Demo_CS_GBP.rar_DEMO_gbp

function Demo_CS_GBP() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % the DCT basis is selected as the sparse representation dictionary % instead of seting the whole image as a vector, I process the image in the % fashion of column-by-column, so as to reduce the complexity. % Author: Chengfu Huo, roy@mail.ustc.edu.cn, http://home.ustc.edu.cn/~roy % The key code of GBP is provided by P. Huggins % Reference: P. Huggins and S. Zucker, ��Greedy Basis Pursuit,�� 2006. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %------------ read in the image -------------- img=imread('lena.bmp'); % testing image img=double(img); [height,width]=size(img); %------------ form the measurement matrix and base matrix --------------- Phi=randn(floor(height/3),width); % only keep one third of the original data Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(height/3),1]); % normalize each column mat_dct_1d=zeros(256,256); % building the DCT basis (corresponding to each column) for k=0:1:255 dct_1d=cos([0:1:255]'*k*pi/256); if k>0 dct_1d=dct_1d-mean(dct_1d); end; mat_dct_1d(:,k+1)=dct_1d/norm(dct_1d); end %--------- projection --------- img_cs_1d=Phi*img; % treat each column as a independent signal %-------- recover using omp ------------ sparse_rec_1d=[]; Theta_1d=Phi*mat_dct_1d; for i=1:width column_rec=cs_gbp(img_cs_1d(:,i),Theta_1d,height); sparse_rec_1d=[sparse_rec_1d,column_rec']; % sparse representation end img_rec_1d=[mat_dct_1d,-mat_dct_1d]*sparse_rec_1d; % inverse transform %------------ show the results -------------------- figure(1) subplot(2,2,1),imagesc(img),title('original image') subplot(2,2,2),imagesc(Phi),title('measurement mat') subplot(2,2,3),imagesc(mat_dct_1d),title('1d dct mat') psnr = 20*log10(255/sqrt(mean((img(:)-img_rec_1d(:)).^2))) subplot(2,2,4),imagesc(img_rec_1d),title(strcat('1d rec img ',num2str(psnr),'dB')) disp('over') %************************************************************************% function hat_x=cs_gbp(y,T_Mat,m) % y=T_Mat*x, T_Mat is n-by-m % y - measurements % T_Mat - combination of random matrix and sparse representation basis % m - size of the original signal % the sparsity is length(y)/4 cnt=length(y); x=y'; % x is row vector D=[T_Mat,-T_Mat]'; % each row of D is a atom OPT_VERBOSE_ON = 0; % commentary level epsilon = 1.0e-6; % error tolerance [n, d] = size(D); % problem size pinv_epsilon = 1.0e-12; %% epsilon for PINV_ITER exit_flag = 0; % default exit_flag % Compute the initial atom t = 1; % iteration index [alpha_0, k] = max(D*x'); i_list = [k]; alpha_list = [alpha_0]; n_atoms = length(i_list); % number of atoms in the current representation normal = x / sqrt(sum(x.^2)); % the normal to the hyperplane xApprox = alpha_0*D(k,:); % the current approximation to x d_H = dot(D(k,:), normal); % the distance of the hyperplane to the origin xApprox_H = (d_H / dot(xApprox, normal))*xApprox; % xApprox projected onto the hyperplane err_vec = x - xApprox; % v = err_vec - dot(err_vec, normal)*normal; v = v / sqrt(sum(v.^2)); % Psi = [D(k,:)]; Psi_biortho = Psi; % Iterate for times=1:cnt/4 % Exit when error is acceptable error = sqrt(sum(err_vec.^2)); if OPT_VERBOSE_ON, fprintf('t = %i; i_t = %i; nAtoms = %i; error = %f\n', t, k, n_atoms, error); end if error < epsilon, fprintf('... GBP complete.\n'); break; elseif isnan(error), fprintf('GBP -- Error: Exiting.\n'); exit_flag = -2; end % Exit when number of atoms equals dimension if n_atoms==d, fprintf('Basis found.\n'); break; end % Otherwise, next iteration t = t + 1; % Project atoms onto n-v plane indices = setdiff(1:n, i_list); %delete the elements of i_list from 1:n D_n0 = D(indices,:) * normal'; % select the next atom from the residuals D_v0 = D(indices,:) * v'; D_n = D_n0 - repmat(dot(xApprox_H, normal), length(indices), 1); D_v = D_v0 - repmat(dot(xApprox_H, v), length(indices), 1); % Compute angle D_theta = atan2(-D_n, D_v); D_theta = D_theta + (D_theta < 0)*2*pi; % re-center for min() % Order by angle and select the next intersected atom [junk, i_Dsorted] = sort(D_theta); k = indices(i_Dsorted(1)); psi_k = D(k,:); % Compute auxiliary variables and coefficients % Compute orthogonal component of psi_k psi_k_ortho = psi_k; for i = 1:n_atoms, beta(i) = dot(psi_k, Psi_biortho(i,:)); psi_k_ortho = psi_k_ortho - beta(i)*Psi(i,:); end % Compute new biorthogonal vector and adjust existing biorthogonal vectors psi_k_biortho = psi_k_ortho / sum(psi_k_ortho.^2); for i = 1:n_atoms, Psi_biortho(i,:) = Psi_biortho(i,:) - beta(i)*psi_k_biortho; end % Compute new coefficient and adjust existing coefficients alpha_k = dot(x, psi_k_biortho); for i = 1:n_atoms, alpha_list(i) = alpha_list(i) - beta(i)*alpha_k; end % Check if the new representation is convex (with respect to the new coefficient) % If not, there is a numerical error. if alpha_k < 0, fprintf('GBP -- Error: Hyperplane rotation failed. Exiting. \n'); exit_flag = -1; return; end % Update lists Psi = [Psi; psi_k]; Psi_biortho = [Psi_biortho; psi_k_biortho]; alpha_list = [alpha_list, alpha_k]; i_list = [i_list, k]; n_atoms = n_atoms + 1; %Test biorthogonality AAplus = Psi*Psi_biortho'; Ierr = max(max(abs(AAplus - eye(size(AAplus))))); if Ierr > epsilon, % run pinv_iter with epsilon=epsilon and c=1 X = Psi_biortho'; while Ierr > epsilon, X = X * (2*eye(size(AAplus)) - AAplus); AAplus = Psi*X; Ierr = max(max(abs(AAplus - eye(size(AAplus))))); %fprintf('... pinv: I error: %f\n', Ierr); end % pinv_iter complete Psi_biortho = X'; alpha_list = (Psi_biortho*x')'; end % Test if conv(Psi) contains extraneous atoms, % otherwise, continue wrapping while sum(alpha_list<=0)~=0, if OPT_VERBOSE_ON, fprintf('... Deleting atom. \n'); end % Find an extraneous atom [alpha_j, j] = min(alpha_list); % Delete it psi_j_biortho = Psi_biortho(j,:); gamma = Psi_biortho*psi_j_biortho' / sum(psi_j_biortho.^2); Psi = Psi([1:j-1,j+1:n_atoms],:); Psi_biortho = Psi_biortho - gamma*psi_j_biortho; Psi_biortho = Psi_biortho([1:j-1,j+1:n_atoms],:); % Test biorthogonality AAplus = Psi*Psi_biortho'; Ierr = max(max(abs(AAplus - eye(size(AAplus))))); if Ierr > epsilon, % run pinv_iter with epsilon=epsilon and c=1 X = Psi_biortho'; while Ierr > epsilon, X = X * (2*eye(size(AAplus)) - AAplus); AAplus = Psi*X; Ierr = max(max(abs(AAplus - eye(size(AAplus))))); %fprintf('... pinv: I error: %f\n', Ierr); end % pinv_iter complete Psi_biortho = X'; alpha_list = (Psi_biortho*x')'; else alpha_list = alpha_list - alpha_j*gamma'; alpha_list = alpha_list([1:j-1,j+1:n_atoms]); end % Update the representation i_list = i_list([1:j-1,j+1:n_atoms]); n_atoms = length(i_list); end % Compute the new approximation xApprox = alpha_list * Psi; % Next hyperplane normal = (-D_n(i_Dsorted(1)))*v + (D_v(i_Dsorted(1)))*normal; normal = normal / sqrt(sum(normal.^2)); d_H = dot(D(i_list(1),:)