From 723a2d3fbe9a7a8c145b5f5ef481dcd4a3799383 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Wed, 24 Jan 2018 17:39:38 +0000
Subject: all Matlab related stuff have been moved to wrappers
---
Wrappers/Matlab/supp/RMSE.m | 7 +++
Wrappers/Matlab/supp/my_red_yellowMAP.mat | Bin 0 -> 1761 bytes
Wrappers/Matlab/supp/sino_add_artifacts.m | 33 +++++++++++
Wrappers/Matlab/supp/studentst.m | 47 +++++++++++++++
Wrappers/Matlab/supp/zing_rings_add.m | 91 ++++++++++++++++++++++++++++++
5 files changed, 178 insertions(+)
create mode 100644 Wrappers/Matlab/supp/RMSE.m
create mode 100644 Wrappers/Matlab/supp/my_red_yellowMAP.mat
create mode 100644 Wrappers/Matlab/supp/sino_add_artifacts.m
create mode 100644 Wrappers/Matlab/supp/studentst.m
create mode 100644 Wrappers/Matlab/supp/zing_rings_add.m
(limited to 'Wrappers/Matlab/supp')
diff --git a/Wrappers/Matlab/supp/RMSE.m b/Wrappers/Matlab/supp/RMSE.m
new file mode 100644
index 0000000..002f776
--- /dev/null
+++ b/Wrappers/Matlab/supp/RMSE.m
@@ -0,0 +1,7 @@
+function err = RMSE(signal1, signal2)
+%RMSE Root Mean Squared Error
+
+err = sum((signal1 - signal2).^2)/length(signal1); % MSE
+err = sqrt(err); % RMSE
+
+end
\ No newline at end of file
diff --git a/Wrappers/Matlab/supp/my_red_yellowMAP.mat b/Wrappers/Matlab/supp/my_red_yellowMAP.mat
new file mode 100644
index 0000000..c2a5b87
Binary files /dev/null and b/Wrappers/Matlab/supp/my_red_yellowMAP.mat differ
diff --git a/Wrappers/Matlab/supp/sino_add_artifacts.m b/Wrappers/Matlab/supp/sino_add_artifacts.m
new file mode 100644
index 0000000..f601914
--- /dev/null
+++ b/Wrappers/Matlab/supp/sino_add_artifacts.m
@@ -0,0 +1,33 @@
+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
new file mode 100644
index 0000000..99fed1e
--- /dev/null
+++ b/Wrappers/Matlab/supp/studentst.m
@@ -0,0 +1,47 @@
+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
new file mode 100644
index 0000000..d197b1f
--- /dev/null
+++ b/Wrappers/Matlab/supp/zing_rings_add.m
@@ -0,0 +1,91 @@
+% 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