From 52ed0437bfbf5bb8a6fdafd65628c4460184fd2d Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 3 Apr 2018 12:34:10 +0100 Subject: cleaning stage, leaving only FGP/ROF and related compilers, demos --- Wrappers/Matlab/supp/sino_add_artifacts.m | 33 ----------- Wrappers/Matlab/supp/studentst.m | 47 ---------------- Wrappers/Matlab/supp/zing_rings_add.m | 91 ------------------------------- 3 files changed, 171 deletions(-) delete mode 100644 Wrappers/Matlab/supp/sino_add_artifacts.m delete mode 100644 Wrappers/Matlab/supp/studentst.m delete mode 100644 Wrappers/Matlab/supp/zing_rings_add.m (limited to 'Wrappers/Matlab/supp') diff --git a/Wrappers/Matlab/supp/sino_add_artifacts.m b/Wrappers/Matlab/supp/sino_add_artifacts.m deleted file mode 100644 index f601914..0000000 --- a/Wrappers/Matlab/supp/sino_add_artifacts.m +++ /dev/null @@ -1,33 +0,0 @@ -function sino_artifacts = sino_add_artifacts(sino,artifact_type) -% function to add various distortions to the sinogram space, current -% version includes: random rings and zingers (streaks) -% Input: -% 1. sinogram -% 2. artifact type: 'rings' or 'zingers' (streaks) - - -[Detectors, anglesNumb, SlicesZ] = size(sino); -fprintf('%s %i %s %i %s %i %s \n', 'Sinogram has a dimension of', Detectors, 'detectors;', anglesNumb, 'projections;', SlicesZ, 'vertical slices.'); - -sino_artifacts = sino; - -if (strcmp(artifact_type,'rings')) - fprintf('%s \n', 'Adding rings...'); - NumRings = round(Detectors/20); % Number of rings relatively to the size of Detectors - IntenOff = linspace(0.05,0.5,NumRings); % the intensity of rings in the selected range - - for k = 1:SlicesZ - % generate random indices to propagate rings - RandInd = randperm(Detectors,Detectors); - for jj = 1:NumRings - ind_c = RandInd(jj); - sino_artifacts(ind_c,1:end,k) = sino_artifacts(ind_c,1:end,k) + IntenOff(jj).*sino_artifacts(ind_c,1:end,k); % generate a constant offset - end - - end -elseif (strcmp(artifact_type,'zingers')) - fprintf('%s \n', 'Adding zingers...'); -else - fprintf('%s \n', 'Nothing selected, the same sinogram returned...'); -end -end \ No newline at end of file diff --git a/Wrappers/Matlab/supp/studentst.m b/Wrappers/Matlab/supp/studentst.m deleted file mode 100644 index 99fed1e..0000000 --- a/Wrappers/Matlab/supp/studentst.m +++ /dev/null @@ -1,47 +0,0 @@ -function [f,g,h,s,k] = studentst(r,k,s) -% Students T penalty with 'auto-tuning' -% -% use: -% [f,g,h,{k,{s}}] = studentst(r) - automatically fits s and k -% [f,g,h,{k,{s}}] = studentst(r,k) - automatically fits s -% [f,g,h,{k,{s}}] = studentst(r,k,s) - use given s and k -% -% input: -% r - residual as column vector -% s - scale (optional) -% k - degrees of freedom (optional) -% -% output: -% f - misfit (scalar) -% g - gradient (column vector) -% h - positive approximation of the Hessian (column vector, Hessian is a diagonal matrix) -% s,k - scale and degrees of freedom -% -% Tristan van Leeuwen, 2012. -% tleeuwen@eos.ubc.ca - -% fit both s and k -if nargin == 1 - opts = optimset('maxFunEvals',1e2); - tmp = fminsearch(@(x)st(r,x(1),x(2)),[1;2],opts); - s = tmp(1); - k = tmp(2); -end - - -if nargin == 2 - opts = optimset('maxFunEvals',1e2); - tmp = fminsearch(@(x)st(r,x,k),[1],opts); - s = tmp(1); -end - -% evaulate penalty -[f,g,h] = st(r,s,k); - - -function [f,g,h] = st(r,s,k) -n = length(r); -c = -n*(gammaln((k+1)/2) - gammaln(k/2) - .5*log(pi*s*k)); -f = c + .5*(k+1)*sum(log(1 + conj(r).*r/(s*k))); -g = (k+1)*r./(s*k + conj(r).*r); -h = (k+1)./(s*k + conj(r).*r); diff --git a/Wrappers/Matlab/supp/zing_rings_add.m b/Wrappers/Matlab/supp/zing_rings_add.m deleted file mode 100644 index d197b1f..0000000 --- a/Wrappers/Matlab/supp/zing_rings_add.m +++ /dev/null @@ -1,91 +0,0 @@ -% uncomment this part of script to generate data with different noise characterisitcs - -fprintf('%s\n', 'Generating Projection Data...'); - -% Creating RHS (b) - the sinogram (using a strip projection model) -% vol_geom = astra_create_vol_geom(N, N); -% proj_geom = astra_create_proj_geom('parallel', 1.0, P, theta_rad); -% proj_id_temp = astra_create_projector('strip', proj_geom, vol_geom); -% [sinogram_id, sinogramIdeal] = astra_create_sino(phantom, proj_id_temp); -% astra_mex_data2d('delete',sinogram_id); -% astra_mex_algorithm('delete',proj_id_temp); - -%% -% inverse crime data generation -[sino_id, sinogramIdeal] = astra_create_sino3d_cuda(phantom, proj_geom, vol_geom); -astra_mex_data3d('delete', sino_id); - -% [id,x] = astra_create_backprojection3d_cuda(sinogramIdeal, proj_geom, vol_geom); -% astra_mex_data3d('delete', id); -%% -% -% % adding Gaussian noise -% eta = 0.04; % Relative noise level -% E = randn(size(sinogram)); -% sinogram = sinogram + eta*norm(sinogram,'fro')*E/norm(E,'fro'); % adding noise to the sinogram -% sinogram(sinogram<0) = 0; -% clear E; - -%% -% adding zingers -val_offset = 0; -sino_zing = sinogramIdeal'; -vec1 = [60, 80, 80, 70, 70, 90, 90, 40, 130, 145, 155, 125]; -vec2 = [350, 450, 190, 500, 250, 530, 330, 230, 550, 250, 450, 195]; -for jj = 1:length(vec1) - for i1 = -2:2 - for j1 = -2:2 - sino_zing(vec1(jj)+i1, vec2(jj)+j1) = val_offset; - end - end -end - -% adding stripes into the signogram -sino_zing_rings = sino_zing; -coeff = linspace2(0.01,0.15,180); -vmax = max(sinogramIdeal(:)); -sino_zing_rings(1:180,120) = sino_zing_rings(1:180,120) + vmax*0.13; -sino_zing_rings(80:180,209) = sino_zing_rings(80:180,209) + vmax*0.14; -sino_zing_rings(50:110,210) = sino_zing_rings(50:110,210) + vmax*0.12; -sino_zing_rings(1:180,211) = sino_zing_rings(1:180,211) + vmax*0.14; -sino_zing_rings(1:180,300) = sino_zing_rings(1:180,300) + vmax*coeff(:); -sino_zing_rings(1:180,301) = sino_zing_rings(1:180,301) + vmax*0.14; -sino_zing_rings(10:100,302) = sino_zing_rings(10:100,302) + vmax*0.15; -sino_zing_rings(90:180,350) = sino_zing_rings(90:180,350) + vmax*0.11; -sino_zing_rings(60:140,410) = sino_zing_rings(60:140,410) + vmax*0.12; -sino_zing_rings(1:180,411) = sino_zing_rings(1:180,411) + vmax*0.14; -sino_zing_rings(1:180,412) = sino_zing_rings(1:180,412) + vmax*coeff(:); -sino_zing_rings(1:180,413) = sino_zing_rings(1:180,413) + vmax*coeff(:); -sino_zing_rings(1:180,500) = sino_zing_rings(1:180,500) - vmax*0.12; -sino_zing_rings(1:180,501) = sino_zing_rings(1:180,501) - vmax*0.12; -sino_zing_rings(1:180,550) = sino_zing_rings(1:180,550) + vmax*0.11; -sino_zing_rings(1:180,551) = sino_zing_rings(1:180,551) + vmax*0.11; -sino_zing_rings(1:180,552) = sino_zing_rings(1:180,552) + vmax*0.11; - -sino_zing_rings(sino_zing_rings < 0) = 0; -%% - -% adding Poisson noise -dose = 50000; -multifactor = 0.002; - -dataExp = dose.*exp(-sino_zing_rings*multifactor); % noiseless raw data -dataPnoise = astra_add_noise_to_sino(dataExp, dose); % pre-log noisy raw data (weights) -sino_zing_rings = log(dose./max(dataPnoise,1))/multifactor; %log corrected data -> sinogram -Dweights = dataPnoise'; % statistical weights -sino_zing_rings = sino_zing_rings'; -clear dataPnoise dataExp - -% w = dose./exp(sinogram*multifactor); % getting back raw data from log-cor - -% figure(1); -% set(gcf, 'Position', get(0,'Screensize')); -% subplot(1,2,1); imshow(phantom,[0 0.6]); title('Ideal Phantom'); colorbar; -% subplot(1,2,2); imshow(sinogram,[0 180]); title('Noisy Sinogram'); colorbar; -% colormap(cmapnew); - -% figure; -% set(gcf, 'Position', get(0,'Screensize')); -% subplot(1,2,1); imshow(sinogramIdeal,[0 180]); title('Ideal Sinogram'); colorbar; -% imshow(sino_zing_rings,[0 180]); title('Noisy Sinogram with zingers and stripes'); colorbar; -% colormap(cmapnew); \ No newline at end of file -- cgit v1.2.3