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