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, 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 addpath(genpath('octave_packages'))
  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; (in 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 delta_mat=delta_mat*10^-7;
 164 beta_mat=lambda*mu_mat/(4*pi);
 165
 166 %% Cartes 2D d'attenuation et d/b
 167 mu=zeros(size_image);
 168 for(i=1:size_image)
 169         for(j=1:size_image)
 170                 if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
 171                         mu(i,j)=mu_mat;
 172                 end
 173         end
 174 end
 175 namnam=strcat(dir_out, 'attenuationCirclePET.edf')
 176 edfwrite(namnam, mu, 'float32');
 177
 178
 179 delta=zeros(size_image);
 180 for(i=1:size_image)
 181         for(j=1:size_image)
 182                 if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
 183                         delta(i,j)=delta_mat;
 184                 end
 185         end
 186 end
 187 namnam=strcat(dir_out, 'deltaCirclePET.edf')
 188 edfwrite(namnam, delta, 'float32');
 189
 190 beta=zeros(size_image);
 191 for(i=1:size_image)
 192         for(j=1:size_image)
 193                 if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
 194                         beta(i,j)=beta_mat;
 195                 end
 196         end
 197 end
 198 namnam=strcat(dir_out, 'betaCirclePET.edf')
 199 edfwrite(namnam, beta, 'float32');
 200
 201
 202 if(strcmp(noise_type, 'gaussian'))
 203         PPSNR=noise_amount;
 204         noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
 205 elseif(strcmp(noise_type, 'poisson'))
 206         scale_fact=noise_amount;
 207         noise_str=strcat(noise_type, num2str(scale_fact), '_');
 208 elseif(strcmp(noise_type, ''))
 209         noise_str='';
 210 else
 211         error("Unknown type of noise")
 212 end
 213
 214 ############################# END OF Unchanged parameters ###########################
 215
 216 ############################# Beginning of simulation ###############################
 217 ps=ps/oversamp;
 218
 219         % Create the sinograms
 220 disp('Create sinograms');
 221         if(strcmp(vers, 'Radon'))
 222                 sino_abs=radon(mu, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
 223                 sino_delta=radon(delta, angles)*ps*oversamp;
 224                 sino_beta=radon(beta, angles)*ps*oversamp;
 225
 226         elseif(strcmp(vers, 'Analytical'))
 227                 [sino_abs, sino_delta, sino_beta] = CircleMat_Analytical (angles, ps*oversamp, R_circle, circle_center, mu_mat, beta_mat, delta_mat, size_image, 0);
 228
 229         else
 230                 disp('error in the variable vers')
 231                 return;
 232         endif
 233
 234         % creation of projections, mu et B, taking into account the oversampling
 235 disp('create projections');
 236 abs_proj=cell(1, nbProj);
 237 B_proj=cell(1, nbProj);
 238 Phi_proj=cell(1, nbProj);
 239
 240         for j=1:nbProj
 241                 abs_proj{j}=interp(sino_abs(:,j), oversamp);
 242                 B_proj{j}=(1/2)*abs_proj{j};
 243                 Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
 244         end
 245
 246 % Fourier domain
 247 [mp np]=size(Phi_proj{j});
 248 [ftot]=FrequencySpace1(2*mp, ps);
 249
 250
 251 ###################### CTF model ################
 252 if(model_ctf==1)
 253         disp('CTF model')
 254         IdCTF=cell(length(D), nbProj);
 255         sino_IdCTF=cell(length(D),1);
 256
 257         for k=1:length(D)
 258         k
 259
 260                 for j=1:nbProj
 261                         sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 262                 end
 263
 264                 if(strcmp(noise_type, 'gaussian'))
 265                         ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
 266                         bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
 267                         bruit=bruit/max(max(bruit));
 268                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 269                         sino_IdCTF{k}=sino_IdCTF{k}+bruit;
 270                 elseif(strcmp(noise_type, 'poisson'))
 271                         noisy=scale_fact*poissrnd(sino_IdCTF{k}/scale_fact);
 272                         sino_IdCTF{k}=noisy;
 273                 end
 274
 275                 name_ctf=strcat(basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k),  '_.edf')
 276                 edfwrite(strcat(dir_out, name_ctf), sino_IdCTF{k}, 'float32');
 277         end
 278 end % end of model CTF
 279
 280
 281 ###################### Fresnel model ################
 282
 283 if(model_Fresnel==1)
 284         disp('Fresnel model')
 285
 286         IdFresnel=cell(length(D), nbProj);
 287         sino_IdFresnel=cell(length(D),1);
 288
 289         for k=1:length(D)
 290                 for j=1:nbProj
 291
 292                         tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 293
 294                         IdFresnel{k, j}=tmp;
 295                         sino_IdFresnel{k}(:,j)=tmp;
 296                 end
 297
 298                 if(strcmp(noise_type, 'gaussian'))
 299                         ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
 300                         bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
 301                         bruit=bruit/max(max(bruit));
 302                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 303                         sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
 304                 elseif(strcmp(noise_type, 'poisson'))
 305                         noisy=scale_fact*poissrnd(sino_IdFresnel{k}/scale_fact);
 306                         sino_IdFresnel{k}=noisy;
 307                 end
 308
 309                 name_fre=strcat(basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
 310                 edfwrite(strcat(dir_out, name_fre), sino_IdFresnel{k}, 'float32');
 311         end
 312 end % end of Fresnel model
 313
 314 ############################# End of simulation ###############################