Useful functions for simulations (created by LW, free to use)

[CreaPhase.git] / utilities_LW / PhantTibo_Analytical.m

diff --git a/utilities_LW/PhantTibo_Analytical.m b/utilities_LW/PhantTibo_Analytical.m
new file mode 100644 (file)
index 0000000..b603e62
--- /dev/null
+++ b/utilities_LW/PhantTibo_Analytical.m
@@ -0,0 +1,99 @@
+## Copyright (C) 2015 Loriane Weber
+## 
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## PhantTibo_Analytical
+
+## Author: Loriane Weber <lweber@rnice7-0102>
+## Created: 2015-12-10
+
+%%%%%% INPUT PARAMETERS
+%% ps : is the physical pixel size, in um
+%% RedAttFact : reduced attenuation factor
+%% angles : angles of projections
+%% save : 1 or 0 (1 if you want to save in edf, 0 sinon
+
+function [MuP_tot, Phi_tot, B_tot ] = PhantTibo_Analytical (angles, ps, save)
+
+disp('params');
+
+RedAttFact=1; 
+
+mu_Al=10.7765 * 100/RedAttFact; % *100 to be in m-1 
+mu_Mg= 5.56697 * 100/RedAttFact; 
+mu_PET=0.89889 * 100/RedAttFact; 
+
+delta_Al=15.03*10^-7/RedAttFact;
+delta_Mg=9.917*10^-7/RedAttFact;
+delta_PET=8.265*10^-7/RedAttFact;
+
+beta_Al=5.59*10^-9/RedAttFact;
+beta_Mg=2.89*10^-9/RedAttFact;
+beta_PET=0.46*10^-9/RedAttFact;
+
+R_Al=73/2;
+R_Mg=37/2;
+R_PET=58/2;
+
+cen_Al=[607-600; 585-600]; % 600 to put the center at the center of the image (1200*1200)
+cen_Mg=[679-600; 640-600];
+cen_PET=[730-600; 630-600];
+
+du= 1; % echantillonage de la projection
+
+angles_rad=(angles*pi/180) + pi/2; % in radians % +pi/2 to be coherent with radon transform
+
+dim=1701; % taille de la proj, wrt radon function
+u=[-(dim-1)/2:du:+(dim-1)/2]; 
+disp('proj')
+
+
+for i=1:length(angles)
+
+       u0_Al= cos(angles_rad(i))*cen_Al(1) + sin(angles_rad(i))*cen_Al(2);
+       u0_Mg= cos(angles_rad(i))*cen_Mg(1) + sin(angles_rad(i))*cen_Mg(2);
+       u0_PET= cos(angles_rad(i))*cen_PET(1) + sin(angles_rad(i))*cen_PET(2);
+
+## Proj attenuation
+       MuP_Al=Analytical_1D(mu_Al, u, R_Al, u0_Al);
+       MuP_Mg=Analytical_1D(mu_Mg, u, R_Mg, u0_Mg);
+       MuP_PET=Analytical_1D(mu_PET, u, R_PET, u0_PET);
+
+       MuP_tot(:, i) = ( MuP_Al+ MuP_Mg + MuP_PET ) * ps;
+
+## Proj Phi
+       Phi_Al=Analytical_1D(delta_Al, u, R_Al, u0_Al);
+       Phi_Mg=Analytical_1D(delta_Mg, u, R_Mg, u0_Mg);
+       Phi_PET=Analytical_1D(delta_PET, u, R_PET, u0_PET);
+
+       Phi_tot(:, i) = ( Phi_Al+ Phi_Mg + Phi_PET ) * ps;
+
+## Proj Beta
+       B_Al=Analytical_1D(beta_Al, u, R_Al, u0_Al);
+       B_Mg=Analytical_1D(beta_Mg, u, R_Mg, u0_Mg);
+       B_PET=Analytical_1D(beta_PET, u, R_PET, u0_PET);
+
+       B_tot(:, i) = ( B_Al+ B_Mg + B_PET ) *ps;
+
+end
+
+if save
+       edfwrite('Sino_Analytical_Mu.edf', MuP_tot, 'float32')
+       edfwrite('Sino_Analytical_Phi.edf', Phi_tot, 'float32')
+       edfwrite('Sino_Analytical_B.edf', B_tot, 'float32')
+end
+
+
+endfunction