matlab-基于matlab的CS-OMP压缩感知的图像重建算法matlab仿真(含教程)

function Demo_CS_TVeq() %------------ 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=zeros(height,width); Theta_1d=Phi*mat_dct_1d; for i=1:width column_rec=cs_TVeq(img_cs_1d(:,i),Theta_1d,height); sparse_rec_1d(:,i)=column_rec'; % sparse representation end img_rec_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_TVeq(y,T_Mat,m) n=length(y); % length of measurements s=floor(n/4); % sparsity, i.e. number of iterations hat_x_initial=T_Mat'*y; % initialization hat_x=tveq_logbarrier(s, hat_x_initial, T_Mat, [], y); function xp = tveq_logbarrier(s, x0, A, At, b, lbtol, mu, slqtol, slqmaxiter) largescale = isa(A,'function_handle'); if (nargin < 6), lbtol = 1e-3; end if (nargin < 7), mu = 10; end if (nargin < 8), slqtol = 1e-3; end if (nargin < 9), slqmaxiter = 200; end newtontol = lbtol; newtonmaxiter = 3; % there is error if it is too large N = length(x0); n = round(sqrt(N)); % create (sparse) differencing matrices for TV Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); % starting point --- make sure that it is feasible if (largescale) if (norm(A(x0)-b)/norm(b) > slqtol) disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); AAt = @(z) A(At(z)); w = cgsolve(AAt, b, slqtol, slqmaxiter, 0); if (cgres > 1/2) disp('A*At is ill-conditioned: cannot find starting point'); xp = x0; return; end x0 = At(w); end else if (norm(A*x0-b)/norm(b) > slqtol) disp('Starting point infeasible; using x0 = At*inv(AAt)*y.'); opts.POSDEF = true; opts.SYM = true; [w, hcond] = linsolve(A*A', b, opts); if (hcond < 1e-14) disp('A*At is ill-conditioned: cannot find starting point'); xp = x0; return; end x0 = A'*w; end end x = x0; Dhx = Dh*x; Dvx = Dv*x; t = (0.95)*sqrt(Dhx.^2 + Dvx.^2) + (0.1)*max(sqrt(Dhx.^2 + Dvx.^2)); tau = N/sum(sqrt(Dhx.^2+Dvx.^2)); lbiter = ceil((log(N)-log(lbtol)-log(tau))/log(mu)); disp(sprintf('Number of log barrier iterations = %d\n', lbiter)); totaliter = 0; for ii = 1:lbiter [xp, tp, ntiter] = tveq_newton(s, x, t, A, At, b, tau, newtontol, newtonmaxiter, slqtol, slqmaxiter); totaliter = totaliter + ntiter; tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2)); x = xp; t = tp; tau = mu*tau; end function [xp, tp, niter] = tveq_newton(cnt, x0, t0, A, At, b, tau, newtontol, newtonmaxiter, slqtol, slqmaxiter) largescale = isa(A,'function_handle'); alpha = 0.01; beta = 0.5; N = length(x0); n = round(sqrt(N)); K = length(b); % create (sparse) differencing matrices for TV Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); % auxillary matrices for preconditioning Mdv = spdiags([reshape([ones(n-1,n); zeros(1,n)],N,1) ... reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); Mdh = spdiags([reshape([ones(n,n-1) zeros(n,1)],N,1) ... reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); Mmd = reshape([ones(n-1,n-1) zeros(n-1,1); zeros(1,n)],N,1); % initial point x = x0; t = t0; Dhx = Dh*x; Dvx = Dv*x; ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2); f = sum(t) - (1/tau)*(sum(log(-ft))); niter = 0; done = 0; times=0; while ((~done) && (times<cnt)) % the number of iterations is determined by newtonmaxiter % there is error when it is too larg times=times+1; ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx); ntgt = -tau - t./ft; gradf = -(1/tau)*[ntgx; ntgt]; sig22 = 1./ft + (t.^2)./(ft.^2); sig12 = -t./ft.^2; sigb = 1./ft.^2 - (sig12.^2)./sig22; w1p = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt); wp = [w1p; zeros(K,1)]; if (largescale) % diagonal of H11p dg11p = Mdh'*(-1./ft + sigb.*Dhx.^2) + Mdv'*(-1./ft + sigb.*Dvx.^2) + 2*Mmd.*sigb.*Dhx.*Dvx; afac = max(dg11p); hpfun = @(z) Hpeval(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, afac); [dxv,slqflag,slqres,slqiter] = symmlq(hpfun, wp, slqtol, slqmaxiter); if (cgres > 1/2) disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); xp = x; return end else H11p = Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ... Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ... Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ... Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh; afac = max(diag(H11p)); Hp = full([H11p afac*A'; afac*A zeros(K)]); %keyboard opts.SYM = true; [dxv, hcond] = linsolve(Hp, wp, opts); if (hcond < 1e-14) disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); xp = x; tp = t; return end end dx = dxv(1:N); Dhdx = Dh*dx; Dvdx = Dv*dx; dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx)); % minimum step size that stays in the interior aqt = Dhdx.^2 + Dvdx.^2 - dt.^2; bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt); cqt = Dhx.^2 + Dvx.^2 - t.^2; tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ... (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ]; indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0)); smax = min(1, min(tsols(indt))); s = (0.99)*smax; % line search suffdec = 0; backiter = 0; while (~suffdec) xp = x + s*dx; tp = t + s*dt; Dhxp = Dhx + s*Dhdx; Dvxp = Dvx + s*Dvdx; ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2); fp = sum(tp) - (1/tau)*(sum(log(-ftp))); flin = f + alpha*s*(gradf'*[dx; dt]); suffdec = (fp <= flin); s = beta*s; backiter = backiter + 1; if (backiter > 32) disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)'); xp = x; tp = t; return end end % set up for next iteration x = xp; t = tp; Dvx = Dvxp; Dhx = Dhxp; ft = ftp; f = fp; lambda2 = -(gradf'*[dx; dt]); stepsize = s*norm([dx; dt]); niter = niter + 1; done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Implicit application of Hessian function y = Hpeval(z, A, At, Dh, Dv, Dhx, Dvx,