utilities_LW/PhantTibo_Analytical.m

   1 ## Copyright (C) 2015 Loriane Weber
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with Octave; see the file COPYING.  If not, see
  15 ## <http://www.gnu.org/licenses/>.
  16
  17 ## PhantTibo_Analytical
  18
  19 ## Author: Loriane Weber <lweber@rnice7-0102>
  20 ## Created: 2015-12-10
  21
  22 %%%%%% INPUT PARAMETERS
  23 %% ps : is the physical pixel size, in um
  24 %% RedAttFact : reduced attenuation factor
  25 %% angles : angles of projections
  26 %% save : 1 or 0 (1 if you want to save in edf, 0 sinon
  27
  28 function [MuP_tot, Phi_tot, B_tot ] = PhantTibo_Analytical (angles, ps, save)
  29
  30 disp('params');
  31
  32 RedAttFact=1;
  33
  34 mu_Al=10.7765 * 100/RedAttFact; % *100 to be in m-1
  35 mu_Mg= 5.56697 * 100/RedAttFact;
  36 mu_PET=0.89889 * 100/RedAttFact;
  37
  38 delta_Al=15.03*10^-7/RedAttFact;
  39 delta_Mg=9.917*10^-7/RedAttFact;
  40 delta_PET=8.265*10^-7/RedAttFact;
  41
  42 beta_Al=5.59*10^-9/RedAttFact;
  43 beta_Mg=2.89*10^-9/RedAttFact;
  44 beta_PET=0.46*10^-9/RedAttFact;
  45
  46 R_Al=73/2;
  47 R_Mg=37/2;
  48 R_PET=58/2;
  49
  50 cen_Al=[607-600; 585-600]; % 600 to put the center at the center of the image (1200*1200)
  51 cen_Mg=[679-600; 640-600];
  52 cen_PET=[730-600; 630-600];
  53
  54 du= 1; % echantillonage de la projection
  55
  56 angles_rad=(angles*pi/180) + pi/2; % in radians % +pi/2 to be coherent with radon transform
  57
  58 dim=1701; % taille de la proj, wrt radon function
  59 u=[-(dim-1)/2:du:+(dim-1)/2];
  60 disp('proj')
  61
  62
  63 for i=1:length(angles)
  64
  65         u0_Al= cos(angles_rad(i))*cen_Al(1) + sin(angles_rad(i))*cen_Al(2);
  66         u0_Mg= cos(angles_rad(i))*cen_Mg(1) + sin(angles_rad(i))*cen_Mg(2);
  67         u0_PET= cos(angles_rad(i))*cen_PET(1) + sin(angles_rad(i))*cen_PET(2);
  68
  69 ## Proj attenuation
  70         MuP_Al=Analytical_1D(mu_Al, u, R_Al, u0_Al);
  71         MuP_Mg=Analytical_1D(mu_Mg, u, R_Mg, u0_Mg);
  72         MuP_PET=Analytical_1D(mu_PET, u, R_PET, u0_PET);
  73
  74         MuP_tot(:, i) = ( MuP_Al+ MuP_Mg + MuP_PET ) * ps;
  75
  76 ## Proj Phi
  77         Phi_Al=Analytical_1D(delta_Al, u, R_Al, u0_Al);
  78         Phi_Mg=Analytical_1D(delta_Mg, u, R_Mg, u0_Mg);
  79         Phi_PET=Analytical_1D(delta_PET, u, R_PET, u0_PET);
  80
  81         Phi_tot(:, i) = ( Phi_Al+ Phi_Mg + Phi_PET ) * ps;
  82
  83 ## Proj Beta
  84         B_Al=Analytical_1D(beta_Al, u, R_Al, u0_Al);
  85         B_Mg=Analytical_1D(beta_Mg, u, R_Mg, u0_Mg);
  86         B_PET=Analytical_1D(beta_PET, u, R_PET, u0_PET);
  87
  88         B_tot(:, i) = ( B_Al+ B_Mg + B_PET ) *ps;
  89
  90 end
  91
  92 if save
  93         edfwrite('Sino_Analytical_Mu.edf', MuP_tot, 'float32')
  94         edfwrite('Sino_Analytical_Phi.edf', Phi_tot, 'float32')
  95         edfwrite('Sino_Analytical_B.edf', B_tot, 'float32')
  96 end
  97
  98
  99 endfunction