tools

tools/admm_stop.m

@@ -0,0 +1,18 @@
+function  [r, s,h]= admm_stop(a,b,c,bold,rho,i_d,h)
+a_length = numel(a);
+
+ABSTOL = 1e-3;
+RELTOL = 1e-3; 
+h.r_norm(i_d) = norm(a(:)-b(:));
+h.s_norm(i_d) = norm(-rho*(b(:)-bold(:)));      
+h.eps_pri(i_d)=sqrt(a_length)*ABSTOL+RELTOL*max(norm(a(:)),norm(b(:)));
+h.eps_dual(i_d)=sqrt(a_length)*ABSTOL+RELTOL*norm(c(:));
+r = h.r_norm(i_d) < h.eps_pri(i_d);
+s = h.s_norm(i_d) < h.eps_dual(i_d);   
+%     if h.r_norm(i_d)>10*h.s_norm(i_d)
+%         rho = rho*2;
+%     elseif h.s_norm(i_d)>10*h.r_norm(i_d)
+%         rho = rho/2;
+%     end 
+
+end

tools/auto_para_apg.m

@@ -0,0 +1,35 @@
+function PARA = auto_para_apg(R,K,psf_s,b,des,precS,use_gpu,tol)
+PARA = [];
+PARA.verbose = des; 
+PARA.precS = precS;
+PARA.gpu = use_gpu;
+PARA.K=K;    
+PARA.psf_s = psf_s; 
+PARA.psf_radius = floor( PARA.psf_s/2 );
+PARA.tol = tol;
+lambda_residual = 1; 
+lambda_l1 = 0.1; 
+%%
+PARA.lambda = [lambda_residual, lambda_l1];
+PARA.max_it = 100;
+PARA.max_it_d=100;
+PARA.max_it_z=100;
+
+PARA.n =1;
+PARA.size_ori_x = size(b);
+PARA.N = size(b,3);
+PARA.size_x = [size(b,1) + 2*PARA.psf_radius(1), size(b,2) + 2*PARA.psf_radius(2),PARA.n];
+PARA.size_z = [PARA.size_x(1), PARA.size_x(2), PARA.K,PARA.n];
+PARA.size_k = [2*PARA.psf_radius(1) + 1, 2*PARA.psf_radius(2) + 1,PARA.K]; 
+PARA.size_k_full = [PARA.size_x(1), PARA.size_x(2), PARA.K]; 
+PARA.kernel_size = [PARA.psf_s,PARA.K];
+%% 
+mul_heur = 5000;%%%%% use 5000 to recover ori.
+gamma_heuristic = mul_heur* 1/(100*max(b(:)));
+PARA.rho_D = gamma_heuristic;
+PARA.rho_Z = gamma_heuristic;
+PARA.max_it_w = 100;
+PARA.rho_W = gamma_heuristic; 
+PARA.R=R; % used to be baseK
+PARA.size_y = [PARA.size_x(1), PARA.size_x(2), PARA.R,PARA.n];
+end

tools/d2dsmall.m

@@ -0,0 +1,25 @@
+function d_small = d2dsmall(d,para)
+if size(para.psf_s,2) == 2
+    d_small = circshift( d, [para.psf_radius,0] ); 
+    %d_small = d_small(1:para.psf_radius*2+1,:);
+
+    if mod(para.psf_s(1),2)==1
+        left_idx = 1:para.psf_radius(1)*2+1;
+    else
+        left_idx = 1:para.psf_radius(1)*2+1;
+    end
+
+    if mod(para.psf_s(2),2)==1
+        right_idx = 1:para.psf_radius(2)*2+1;
+    else
+        right_idx = 1:para.psf_radius(2)*2;
+    end
+
+    d_small = d_small(left_idx,right_idx,:);
+
+else 
+    keyboard
+%     d_small = circshift( d, [para.psf_radius, para.psf_radius,0] ); 
+%     d_small = d_small(1:para.psf_radius*2+1, 1:para.psf_radius*2+1,:);    
+end
+end

tools/dsmall2d.m

@@ -0,0 +1,6 @@
+function d = dsmall2d(d_small,para)
+
+    ndim = length( para.size_z ) - 2;
+    d = padarray( d_small, [para.size_x(1:ndim) - para.psf_s, 0], 0, 'post');
+    d = circshift(d, -[para.psf_radius, 0] ); 
+end

tools/init_SWZ.m

@@ -0,0 +1,18 @@
+function  [W,s_small,z] = init_SWZ(para)
+
+    [W,s_small] = initWR(para); %% for test only
+     z = randn(para.size_z);
+
+    if (para.precS ==1)
+        s_small = single(s_small);
+        W = single(W);
+        z = single(z);
+    end
+    if (para.gpu ==1)
+        s_small = gpuArray(s_small);
+        W = gpuArray(W);
+        z = gpuArray(z);
+    end  
+
+
+end

tools/my_psnr.m

@@ -0,0 +1,22 @@
+function [ps,rmse] = my_psnr(b,Dz,par)
+ndim = length(par.size_z)-2;
+size_Dz = size(Dz);
+it = sqrt(prod(size_Dz(1:ndim)));
+p1 =20*log10(it);
+
+ndim = length(par.size_z)-2;
+pp = prod(par.size_ori_x(1:ndim));
+
+b = reshape(b,pp,[]) ; 
+Dz = reshape(Dz,pp,[]) ; 
+
+for i=1:size(b,2)
+        b_tmp = b(:,i);
+        Dz_tmp = Dz(:,i);
+        tmp =  norm((b_tmp(:) - Dz_tmp(:)));
+        p2 = 20*log10(tmp);
+        p(i) = p1 - p2;
+end
+ps=mean(p);
+rmse = sqrt( 1/length(b(:)) * (tmp^2)); 
+end

tools/plot_cmp.m

@@ -0,0 +1,12 @@
+function plot_cmp(x,reconstruction)
+figure()
+if size(x,3)~=1
+    x=x(:,:,1);
+    reconstruction = reconstruction(:,:,1);
+ end   
+
+subplot(1,2,1), imshow(x);  axis image, colormap gray, title('Orig');
+subplot(1,2,2), imshow(reconstruction); axis image, colormap gray; title(sprintf('Rec'));
+%psnr_see = psnr(x,reconstruction);
+
+end

tools/rec.m

@@ -0,0 +1,8 @@
+function [rec]=rec(z,d,par)
+rec = zeros([par.size_x(1),par.size_x(2),par.n]);
+for i_n = 1:size(z,4)
+    for i_k = 1:size(z,3)
+        rec(:,:,i_n) = rec(:,:,i_n)+conv2(z(:,:,i_k,i_n),d(:,:,i_k),'full');
+    end
+end
+end

tools/rec_fre.m

@@ -0,0 +1,5 @@
+function [ Dz_small,Dz] = rec_fre(d_hat,z_hat,par)
+Dz = real(ifft2( reshape(sum(d_hat .* z_hat, 3), par.size_x) ));
+Dz_small = x2xsmall(Dz,par);
+end
+

tools/solve_conv_term_D.m

@@ -0,0 +1,15 @@
+function dd_hat = solve_conv_term_D(sub, Ahi,Bhi, par,hk)
+    sy = par.size_z(1); sx = par.size_z(2); 
+    left = Bhi+permute( reshape(sub, sx * sy, 1, hk), [3,2,1] );
+    clear Bhi ss_hat yy_hat
+    if par.gpu==1
+        x2 = pagefun(@mtimes,Ahi,left);
+    else
+%         for i = 1:size(Ahi,3)
+%             x2(:,:,i) = Ahi(:,:,i)*left(:,:,i);
+%         end
+        x2 = mtimesx(Ahi,left);
+    end
+    clear Ahi
+    dd_hat = reshape(permute(x2,[3,1,2]),sy,sx,[]);
+end

tools/sparse_level.m

@@ -0,0 +1,5 @@
+function  sp = sparse_level(z,par)     
+z_length = prod(par.size_z);
+temp1 = z(:);
+sp = length( temp1(abs(temp1)<1e-5))/z_length *100;         
+end

tools/x2xsmall.m

@@ -0,0 +1,7 @@
+function x = x2xsmall(x,par)
+ndim = length(par.size_z)-2;
+if ndim == 2
+x = x(1 + par.psf_radius(1):end - par.psf_radius(1),1 + par.psf_radius(2):end - par.psf_radius(2),:);
+end
+
+end

tools/z0u1.m

@@ -0,0 +1,10 @@
+function [nb,veam,vstd] = z0u1(b)
+% zero mean and unit variance
+size_b=size(b);
+p=prod(size_b(3:end));
+vmsi = permute(reshape(double(b), [], p), [2 1]);
+veam = mean(vmsi,2);
+vstd = std(vmsi, 0, 2);
+vmsi = (vmsi-repmat(veam, [1 size(vmsi,2)])) ./ repmat(vstd, [1 size(vmsi,2)]);
+nb = reshape(vmsi', size_b);    
+end