% Makes CS sensing and reconstruction considering the DWT as the sparse domain
%
% Inputs:
%	M: Number of measurements
%	q_step: quantization step
%	epsilon: optimization parameter
%	ImgName: image file name 
%
% Writen by: Adriana Schulz, Visgraf - IMPA aschulz@impa.br
% Last Update: January 2009

function [AX,ry,OM,OMSORT,sn,quantyz,ysn] = cs_l1dctart(M, M2, q_step_dct, q_step, X, recon, sortow)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%  INITIALIZATION  %%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%addpath ../Measurements
%addpath ../Optimization
%addpath ../Utils
%addpath ../Data

%clear all;

x = X(:);
n = size(X,1);
N = n*n;
%qmf = MakeONFilter('Coiflet',2);
%S = FWT2_PO(X,3,qmf);
%s = S(:);
ysn = 0;
% DCT -- for linear approximation of DCT
lporder = bdct_linapprox_ordering(n, n);
% low pass DCT2, K1 = number of lowpass coefficients
OM = lporder(1:M);
cy = A_dct2(x, n, OM); %wspo³czynniki DCT 
 quantyz = ones(M,1)*q_step_dct;
 przyrost = 3;%9;
 for i=2:M,
    quantyz(i) = quantyz(i) + przyrost*q_step_dct*i/M;
 end
 cy = cy./quantyz; 
 ry = round(cy); 
 
 
 if sortow == 1,
    cy= ry.*quantyz;%q_step;
    ax = round(At_dct2(cy, n, OM));
    ndct = realnoiselet(ax)/sqrt(N); %x-ax
    [sndct,ind] = sort(abs(ndct),'descend');
    OMSORT = ind(1:M2);
    sn = ndct(ind(1:M2))/q_step; 
    sn = round(sn);  
 elseif sortow == 2,
    cy= ry.*quantyz;%ry*q_step;
    ax = round(At_dct2(cy, n, OM));
    AX = reshape(ax, n, n);
    %bb=ncanny(AX,'canny');

    f0 = [
        1  -2  1
        -2  4  -2
        1  -2  1
        ];

     pom1 = round(conv2(AX,f0,'valid'));%'valid'
     
f0 = [
        1  -2  1
        -2  5  -2
        1  -2  1
        ];

     pom2 = round(conv2(AX,f0,'valid'));%'valid'     

     
    %unsharp masking
%     f4=[1 4 7 4 1
%     4  16  26  16  4
%     7  26  41  26  7
%     4  16  26  16  4
%     1  4  7  4  1]/273;

%     f4=[1 3 5 3 1
%     3  8  12  8  3
%     5  12  16  12  5
%     3  8  12  8  3
%     1  3  5  3  1]/144;

    f4 = [ 1  2  1
           2  4  2
           1  2  1 ]/16;
%     f4 = [ 1  1  1
%            1  2  1
%            1  1  1 ]/10;       

%     f4 = [ 1  2  1
%            2  5  2
%            1  2  1 ]/17;

    Y4=round(conv2(AX,f4,'valid'));
    %Y4=conv2(X,f4,'valid');
    alfa = 0.2;
    
   
     %Y5 = (1-alfa)*X(1+px:rozmX-px,1+py:rozmY-py) + alfa*(X(1+px:rozmX-px,1+py:rozmY-py)-Y4);
   % Y5 = AX + alfa*(AX-Y4);
    %srednia = mean2(AX);
    pom = zeros(n);
    pom(2:n-1,2:n-1) = AX(2:n-1,2:n-1) - Y4;
    
    pomx = zeros(n);
    pomx(120:250,100:250) = pom2(120:250,100:250);
%     pom1 = (AX - srednia) + pom; 
%     pom1 = pom1/max(max(abs(pom1)))* (255-srednia);
%     pom1 = pom1 + srednia;
   % pom1 = AX + alfa*pom; 
%     kodow = zeros(256);
%     kodow(2:255,2:255) = pom1;  
    
    aax = pomx(:);%/max(pom(:))*255;    
    ndct = realnoiselet(aax)/sqrt(N); %x-ax    aax
    [sndct,ind] = sort(abs(ndct),'descend');
    OMSORT = ind(1:M2);
    OMSORT = find(aax>0);%*************
    
    pom1 = zeros(n)+ round(mean2(Y4)); pom1(2:n-1,2:n-1) = Y4;
    kodow = ax;%pom1(:);%pom;%ax;
    %kodow(2:255,2:255) = pom2;
    %kodow = round(AX + alfa*pom);
    ndct = realnoiselet(kodow(:))/sqrt(N);
     %ndct = realnoiselet(ax)/sqrt(N);
     ysn = ndct(OMSORT);
     sn = ndct(OMSORT)/q_step; %ind(1:M2)
     sn = round(sn);
 else
    OMSORT = 0; 
    sn = 0;
 end
 
 if recon ==1,
    cy= ry.*quantyz;%ry*q_step;
    ax = round(At_dct2(cy, n, OM));
    AX = reshape(ax, n, n);
    AX(AX>255) = 255;
    AX(AX<0) = 0;
 else
    AX = 0;
 end;
 
 
 
 %disp(sprintf('Entropy (wspó³czynniki DCT: %d %d',M*entropy(ry)/N));
%     h = fspecial('laplacian', 0.2);
%     RXpom = round(abs(imfilter(AXpom,h,'replicate')));
  %  mx = ax;
%     bb=ncanny(AXpom,'canny');
%     mx = bb(:);
%     %mx = RXpom(:);
%     
% %     f0=[0  0  -1  0  0
% %         0  -1  -2  -1  0
% %         -1  -2  16  -2  -1
% %         0  -1  -2  -1  0
% %         0  0  -1  0  0
% %         ]/16;
%     f0=[
%         -1  -2  -1
%         -2  12  -2
%         -1  -2  -1
%         ]
% % 
%      %Y4=uint8(conv2(AXpom,f4,'valid')+128);
%      %bb=bb*255;
% %      K2 = 100;
%      Y5=abs(conv2(bb,f0,'same'));%'valid'
%      mx = Y5(:);
%     ndct = realnoiselet(mx)/sqrt(N); %x-ax
%     [sndct,ind] = sort(abs(ndct),'descend');
% %     posort1 = zeros(1,65536);
% %     posort1(ind(1:K2)) = 1;
%     %posort1 = sort(ind(1:K2));
%     OMDCT = ind(1:K1);%K1);%OMDCT