% 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) % recon ustala, czy od razu robiona jest rekonstrukcja - dajemy 0; % natomiast sortow ustala sposób estymacji położeń pomiarów punktowych na % podstawie aproksymacji obrazu z pomiarów kosinusowych - dajemy 2 x = X(:); n = size(X,1); N = n*n; %qmf = MakeONFilter('Coiflet',2); %S = FWT2_PO(X,3,qmf); %s = S(:); ysn = 0; if M ~= 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 if q_step_dct ~= 0, %nieistotna kwantyzacja 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); else ry = cy; quantyz = 0; end if sortow == 1, if q_step_dct ~= 0, cy= ry.*quantyz;%q_step; end ax = round(At_dct2(cy, n, OM)); ndct = realnoiselet(ax)/sqrt(N); %x-ax [sndct,ind] = sort(abs(ndct),'descend'); OMSORT = ind(1:M2); if q_step ~= 0, sn = ndct(ind(1:M2))/q_step; sn = round(sn); else sn = ndct(OMSORT); % ndct(ind(1:M2)); end ysn = sn; elseif sortow == 2, %OPTYMALIZOWANA PROCEDURA ESTYMACJI POŁOŻEŃ POMIARÓW PUNKTOWYCH !!! if q_step_dct ~= 0, cy= ry.*quantyz;%ry*q_step; end ax = round(At_dct2(cy, n, OM)); %PRZYBLIŻENIE OBRAZU Z POMIARÓW KOSINUSOWYCH 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 1 1 1 1 % 1 1 1 1 1 % 1 1 1 1 1 % 1 1 1 1 1 % 1 1 1 1 1]/25; % f4 = [ 1 2 1 % 2 4 2 % 1 2 1 ]/16; f4 = [ 1 1 1 1 1 1 1 1 1 ]/9; % f4 = [ 1 5 1 % 5 10 5 % 1 5 1 ]/34; % 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); pomL = zeros(n); pom(2:n-1,2:n-1) = AX(2:n-1,2:n-1) - Y4; % PROSTY UNSHARP MASKING - WYOSTRZENIE SZCZEGÓŁÓW %pom(3:n-2,3:n-2) = AX(3:n-2,3:n-2) - Y4; pom2 = round(conv2(pom,f0,'valid'));%'valid' %*****DODATKOWE LAPLASJANOWANIE pom(2:n-1,2:n-1) = pom2; %**** % % pomL(100:225,40:130) = pom(100:225,40:130); %umieszczenie tylko ROI filtrowanego obrazu w macierzy transf. noiseletowej % % pom = pomL; %operacja podstawienia pod pom pomx = zeros(n); % % %probka - wybór pomiarów tylko z ROI (120:225,100:130) % % %pomx(120:225,100:130) = 1; %pomx(183:227,117:170) = 1; % % %%pomx(100:220,15:100) = 1; % % pomx(100:225,40:130) = 1; % % px = pomx(:); % % ind = find(px)'; % % % % load RANDOM_STATES % % rand('state', rand_state); % % randn('state', randn_state); % % q = randperm(size(ind,2))'; % % z = q(1:M2); % % zz = ind(z); % % OMSORT = zz; % % %probka %pomx(2:n-1,2:n-1) = pom2; pomx = pom; %***pom %pomx(120:250,100:250) = pom(120:250,100:250); %pomx(100:156,100:156) = pom(100:156,100:156);%pomx(120:250,100:250) = pom(120:250,100:250); %pomx(120:225,100:130) = pom2(120:225,100:130);% % 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(:);%pomx /max(pom(:))*255; ndct = realnoiselet(aax)/sqrt(N); %x-ax PRÓBNY POMIAR NOISELETOWY Z TAK PRZYGOTOWANEJ ESTYMACJI ISTOTNYCH CECH OBRAZOWYCH!!! [sndct,ind] = sort(abs(ndct),'descend'); OMSORT = ind(1:M2); % WSPISANIE DO WEKTORA POŁOŻEŃ NAJWIĘKSZYCH CO DO MODUŁU POMIARÓW Z ESTYMACJI - PRZEYWIDUJEMY, ŻE W TYCH POŁOŻENIACH ZNAJDUJE SIĘ NAWIĘCEJ ENERGII SYGNAŁU %kk = 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); %WEKTOR WARTOŚCI PREDYKCJI POMIARÓW RZECZYWISTYCH W DOBRANYCH ADAPTACYJNIE POŁOŻENIACH NA PODSTAWIE KOSINUSOWEJ ESTYMACJI ORYGINAŁU (teraz nieistotne, bo w kodowaniu pozwalało to zmniejszych bit rate różnicowych wartości) % %dodatek % vn = zeros(N,1); % vn(OMSORT) = ysn; % v = realnoiselet(vn/sqrt(N)); % if q_step ~= 0, %nieistotny wątek kwantyzacji sn = ndct(OMSORT)/q_step; %ind(1:M2) sn = round(sn); else sn = ndct(OMSORT); end else OMSORT = 0; sn = 0; end else %alternatywne losowe ustalenie położeń - to właśnie staramy się poprawić pomx = zeros(n); %probka - wybór pomiarów tylko z ROI (120:225,100:130) %pomx(120:225,100:130) = 1; pomx(183:227,117:170) = 1; %pomx(100:220,15:100) = 1; px = pomx(:); ind = find(px)'; load RANDOM_STATES rand('state', rand_state); randn('state', randn_state); q = randperm(size(ind,2))'; z = q(1:M2); zz = ind(z); OMSORT = zz; %probka ysn = 0; sn = 0; ry=0; OM=0;quantyz=0; end if recon ==1, if q_step_dct ~= 0, cy= ry.*quantyz;%ry*q_step; end ax = round(At_dct2(cy, n, OM)); AX = reshape(ax, n, n); AX(AX>255) = 255; AX(AX<0) = 0; else AX = 0; end;