SimuPBI_circle_func.m

   1 ## Copyright (C) 2016 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 ## SimuPBI_circle_func: fonction that launches the phase-contrast simulation of a wire.
  18
  19 ## Author: Loriane Weber <lweber@gpid16a-1802>
  20 ## Created: 2016-02-11
  21
  22 ## Example of use :
  23 ## SimuPBI_circle_func('Analytical', 2, 'test', [0 0.1 0.5], 19, 2.7, 300, 180, 1, 1, 0.89, 8.265*10^-7, 100, 500, [200 250], '', 'gaussian', 35)
  24
  25 ## Please note that dir_out, noise_type, noise_amount are optional arguments.
  26
  27 function [] = SimuPBI_circle_func(vers, oversamp, basename_output, dist, energy, ps, nbProj, range_angle, model_ctf, model_Fresnel, mu_mat, delta_mat, R_circle, size_image, circle_center, dir_out, noise_type, noise_amount)
  28
  29 addpath('utilities_ESRF')
  30 addpath('utilities_LW')
  31
  32 ################################################
  33 ############### INPUT parameters ###############
  34 ################################################
  35
  36
  37 ################################
  38 ######### Parameter related to the computation - and names of the result
  39 ################################
  40 ## vers: are the projections calculated analytically or using the Radon transform? use the strings 'Radon' OR 'Analytical'
  41 %vers='Analytical'
  42
  43 ## Oversampling of the projections : use 2 or 4
  44 %oversamp=2
  45
  46 ## Basename of the result file. Is a string.
  47 %basename_output='Test-1D';
  48
  49 ################################
  50 ######### Parameter related to the physics
  51 ################################
  52
  53 ## Distances of propagation (in m)
  54 %dist=[0 0.01 0.1 0.20 0.50];
  55
  56 ## Energy of the incoming X-ray beam (in keV)
  57 %energy=19
  58
  59 ## Pixel size of the detector (in um)
  60 %ps=1;
  61
  62 # Number of projections
  63 %nbProj=360
  64
  65 ## Range of the tomography : 180 or 360 degrees
  66 %range_angle=180;
  67
  68 ## Which model do you want to use for the propagation? use 1 or 0
  69 %model_ctf=1;
  70 %model_Fresnel=1;
  71
  72
  73 ################################
  74 ############# Parameters related to the circle
  75 ################################
  76
  77 ## material of the circle (here, PET at 19 keV)
  78 %mu_mat=0.89889 % (in cm-1);
  79 %delta_mat=8.265*10^-7;
  80
  81 ## size and coordinates
  82 %R_circle=50; % radius of the circle
  83 %size_image=500; % size of the image (square)
  84 %circle_center=[200; 250]; % 2D coordinates of the center
  85
  86 ################################
  87 ############# OPTIONAL Parameters
  88 ################################
  89
  90 ## Output directory
  91 %dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_Circle'
  92 % default value if '', i.e. the working directory
  93
  94 ## Noise-related part
  95 ## Noise addition to the simulated data?
  96 ## Use noise='gaussian' (addition of gaussian noise) or noise='poisson' ( generation of Poisson noise) or '' (no noise is added).
  97
  98 %noise_type='gaussian'
  99 % default value is '' (no noise)
 100
 101 ## Amount of noise
 102 ## If 'gaussian' (additive noise), please specify the Peak-to-peak Signe-to-noise ratio (PPSNR, in dB)
 103 % noise_amount=40;
 104 %default value is 35 dB
 105 ## if 'poisson', please specify the amount of noise
 106 %noise_amount=0.05;
 107 % default value is 5%.
 108
 109 #######################################################
 110 ############### End of INPUT parameters ###############
 111 #######################################################
 112
 113 #################################################
 114 ############### OUTPUT parameters ###############
 115 #################################################
 116
 117 ## None, files are directly save
 118
 119 ########################################################
 120 ############### End of OUTPUT parameters ###############
 121 ########################################################
 122
 123
 124 ############################# check varargin - optional arguments ###########################
 125  if nargin<15
 126        error("Some arguments are missing.\n");
 127  elseif nargin==15
 128     dir_out='';
 129     noise_type='';
 130  elseif nargin==16
 131         noise_type='';
 132  elseif nargin==17
 133         if strcmp(noise_type, 'gaussian')
 134                 noise_amount=35;
 135         else strcmp(noise_type, 'poisson')
 136         noise_amount=0.05;
 137     endif
 138 endif
 139
 140
 141 ############################# Unchanged parameters ###########################
 142 basename_output=strcat(basename_output, '_');
 143
 144 % calcul des angles
 145 angles=[0:1:nbProj-1];
 146 delta_angle=range_angle/nbProj;
 147 angles=angles*delta_angle;
 148
 149 lambda=12.4*10^-10/energy; % in meter
 150 ps=ps*10^-6; % put the pixel size in meter
 151
 152 z1=145;
 153 D = (z1*dist)./(z1+dist);
 154 betash = lambda*D; % in meter^2
 155 pad_method='extend';
 156
 157 if !strcmp(dir_out, '')
 158         dir_out= strcat(dir_out, '/'); % add underscore at the end
 159 end
 160
 161 % absorption and phase parameter
 162 mu_mat=mu_mat*100; % in m-1 , to be consistant with lambda!!
 163 beta_mat=lambda*mu_mat/(4*pi);
 164
 165 %% Cartes 2D d'attenuation et d/b
 166 mu=zeros(size_image);
 167 for(i=1:size_image)
 168         for(j=1:size_image)
 169                 if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
 170                         mu(i,j)=mu_mat;
 171                 end
 172         end
 173 end
 174 namnam=strcat(dir_out, 'attenuationCirclePET.edf')
 175 edfwrite(namnam, mu, 'float32');
 176
 177
 178 delta=zeros(size_image);
 179 for(i=1:size_image)
 180         for(j=1:size_image)
 181                 if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
 182                         delta(i,j)=delta_mat;
 183                 end
 184         end
 185 end
 186 namnam=strcat(dir_out, 'deltaCirclePET.edf')
 187 edfwrite(namnam, delta, 'float32');
 188
 189 beta=zeros(size_image);
 190 for(i=1:size_image)
 191         for(j=1:size_image)
 192                 if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
 193                         beta(i,j)=beta_mat;
 194                 end
 195         end
 196 end
 197 namnam=strcat(dir_out, 'betaCirclePET.edf')
 198 edfwrite(namnam, beta, 'float32');
 199
 200
 201 if(strcmp(noise_type, 'gaussian'))
 202         PPSNR=noise_amount;
 203         noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
 204 elseif(strcmp(noise_type, 'poisson'))
 205         scale_fact=noise_amount;
 206         noise_str=strcat(noise_type, num2str(scale_fact), '_');
 207 elseif(strcmp(noise_type, ''))
 208         noise_str='';
 209 else
 210         error("Unknown type of noise")
 211 end
 212
 213 ############################# END OF Unchanged parameters ###########################
 214
 215 ############################# Beginning of simulation ###############################
 216 ps=ps/oversamp;
 217
 218         % Create the sinograms
 219 disp('Create sinograms');
 220         if(strcmp(vers, 'Radon'))
 221                 sino_abs=radon(mu, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
 222                 sino_delta=radon(delta, angles)*ps*oversamp;
 223                 sino_beta=radon(beta, angles)*ps*oversamp;
 224
 225         elseif(strcmp(vers, 'Analytical'))
 226                 [sino_abs, sino_delta, sino_beta] = CircleMat_Analytical (angles, ps*oversamp, R_circle, circle_center, mu_mat, beta_mat, delta_mat, size_image, 0);
 227
 228         else
 229                 disp('error in the variable vers')
 230                 return;
 231         endif
 232
 233         % creation of projections, mu et B, taking into account the oversampling
 234 disp('create projections');
 235 abs_proj=cell(1, nbProj);
 236 B_proj=cell(1, nbProj);
 237 Phi_proj=cell(1, nbProj);
 238
 239         for j=1:nbProj
 240                 abs_proj{j}=interp(sino_abs(:,j), oversamp);
 241                 B_proj{j}=(1/2)*abs_proj{j};
 242                 Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
 243         end
 244
 245 % Fourier domain
 246 [mp np]=size(Phi_proj{j});
 247 [ftot]=FrequencySpace1(2*mp, ps);
 248
 249
 250 ###################### CTF model ################
 251 if(model_ctf==1)
 252         disp('CTF model')
 253         IdCTF=cell(length(D), nbProj);
 254         sino_IdCTF=cell(length(D),1);
 255
 256         for k=1:length(D)
 257         k
 258
 259                 for j=1:nbProj
 260                         sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 261                 end
 262
 263                 if(strcmp(noise_type, 'gaussian'))
 264                         ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
 265                         bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
 266                         bruit=bruit/max(max(bruit));
 267                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 268                         sino_IdCTF{k}=sino_IdCTF{k}+bruit;
 269                 elseif(strcmp(noise_type, 'poisson'))
 270                         noisy=scale_fact*poissrnd(sino_IdCTF{k}/scale_fact);
 271                         sino_IdCTF{k}=noisy;
 272                 end
 273
 274                 name_ctf=strcat(basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k),  '_.edf')
 275                 edfwrite(strcat(dir_out, name_ctf), sino_IdCTF{k}, 'float32');
 276         end
 277 end % end of model CTF
 278
 279
 280 ###################### Fresnel model ################
 281
 282 if(model_Fresnel==1)
 283         disp('Fresnel model')
 284
 285         IdFresnel=cell(length(D), nbProj);
 286         sino_IdFresnel=cell(length(D),1);
 287
 288         for k=1:length(D)
 289                 for j=1:nbProj
 290
 291                         tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 292
 293                         IdFresnel{k, j}=tmp;
 294                         sino_IdFresnel{k}(:,j)=tmp;
 295                 end
 296
 297                 if(strcmp(noise_type, 'gaussian'))
 298                         ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
 299                         bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
 300                         bruit=bruit/max(max(bruit));
 301                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 302                         sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
 303                 elseif(strcmp(noise_type, 'poisson'))
 304                         noisy=scale_fact*poissrnd(sino_IdFresnel{k}/scale_fact);
 305                         sino_IdFresnel{k}=noisy;
 306                 end
 307
 308                 name_fre=strcat(basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
 309                 edfwrite(strcat(dir_out, name_fre), sino_IdFresnel{k}, 'float32');
 310         end
 311 end % end of Fresnel model
 312
 313 ############################# End of simulation ###############################