SimuPBI_3WiresPhant_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_3WiresPhant_func
  18
  19 ## Author: Loriane Weber <lweber@gpid16a-1802>
  20 ## Created: 2016-03-30
  21
  22 ## Example of use :
  23 ## SimuPBI_3WiresPhant_func('Analytical', 2, 'test', [0 0.1 0.5], 19, 2.7, 300, 180, 1, 0, '', 'gaussian', 35)
  24
  25 ## Please note that dir_out, noise_type, noise_amount are optional arguments.
  26
  27 function [] = SimuPBI_3WiresPhant_func (vers, oversamp, basename_output, dist, energy, ps, nbProj, range_angle, model_ctf, model_Fresnel, height, 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 ######### Parameter related to the computation - and names of the result
  40 ################################
  41
  42 ## vers: are the projections calculated analytically or using the Radon transform? use 'Radon' OR 'Analytical'
  43 %vers='Radon' ;
  44
  45 ## Oversampling of the projections : use 2 or 4
  46 %oversamp=2
  47
  48 ## Basename of the result file
  49 %basename_output='phantom_al_mg_PET_thibaut_3dist';
  50
  51 ################################
  52 ######### Parameter related to the physics
  53 ################################
  54
  55 ## Distances of propagation (in m)
  56 %dist=[0 0.28 0.560];
  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=3.5;
  63
  64 ## Number of projections
  65 %nbProj=100;
  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 refers to the Contrast Transfer Function propagation model. use 1 if you want to simulate propagation with the CTF model, or 0 otherwise.
  72 ## should be equal to 0 or 1
  73 % model_ctf = 1
  74 % model_ctf = 0
  75
  76 ## model_Fresnel refers to the Fresnel propagation model. Use 1 if you want to simulate propagation with the fresnel model, or 0 otherwise.
  77 ## should be equal to 0 or 1
  78 % model_Fresnel = 1
  79 % model_Fresnel = 0
  80
  81 ################################
  82 ############# Parameters related to the object
  83 ################################
  84
  85 ## height of the object. If height == 1, the simulation is computed in 1D. If height >=2 , the simulation turns into 2D.
  86 %height=10;
  87
  88
  89 ################################
  90 ############# OPTIONAL Parameters
  91 ################################
  92 ## Output directory : to be removed
  93 %dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_PhantTibo';
  94 %dir_out='';
  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 ## !! note that poisson noise is here not implemented yet
 111
 112 #######################################################
 113 ############### End of INPUT parameters ###############
 114 #######################################################
 115
 116 #################################################
 117 ############### OUTPUT parameters ###############
 118 #################################################
 119
 120 ## None, file directly save
 121
 122 ########################################################
 123 ############### End of OUTPUT parameters ###############
 124 ########################################################
 125
 126 ############################# check varargin - optional arguments ###########################
 127  if nargin<11
 128        error("Some arguments are missing.\n");
 129  elseif nargin==11
 130     dir_out='';
 131     noise_type='';
 132  elseif nargin==12
 133         noise_type='';
 134  elseif nargin==13
 135         if strcmp(noise_type, 'gaussian')
 136                 noise_amount=35;
 137         elseif strcmp(noise_type, 'poisson')
 138         error("poisson noise is not implemented yet")
 139         return; %noise_amount=0.05;
 140         else
 141                 error("Unknown type of noise.\n")
 142     endif
 143 endif
 144
 145
 146 ############# Unchanged parameters ##############
 147 if !strcmp(dir_out, '')
 148         dir_out=strcat(dir_out, '/');
 149 end
 150
 151 z1=145; % source to sample distance (in m)
 152
 153 % calcul des angles
 154 angles=[0:1:nbProj-1];
 155 delta_angle=range_angle/nbProj;
 156 angles=angles*delta_angle;
 157
 158 % Wavelength, in meter
 159 lambda=12.4*10^-10/energy;
 160 ps=ps*10^-6; % put the pixel size in meter
 161
 162 % Effective prop distance
 163 D = (z1*dist)./(z1+dist);
 164 betash = lambda*D; % in meter^2
 165
 166 % pading method
 167 pad_method='extend';
 168
 169 % Absorption image ('mu', in cm-1)
 170 mu=edfread('images/attenuation.edf'); % map of mu in cm-1
 171 [m n]=size(mu)
 172
 173 absorption=mu*100; % in m-1 , to be consistant with lambda!! /10  to have less attenuation
 174 beta=(lambda/(4*pi))*absorption; % one slice of beta /!\ units
 175
 176 delta_beta=edfread('images/delta_beta_map_1200.edf');
 177 delta=delta_beta.*beta;
 178
 179 ## Noise addition
 180 if(strcmp(noise_type, 'gaussian'))
 181         PPSNR=noise_amount;
 182         noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
 183 %elseif(strcmp(noise_type, 'poisson'))
 184         %scale_fact=noise_amount;
 185         %noise_str=strcat(noise_type, num2str(scale_fact), '_');
 186 elseif(strcmp(noise_type, ''))
 187         noise_str='';
 188 end
 189
 190
 191 if height ==1
 192         basename_output=strcat(basename_output, '-1D_');
 193 else
 194         basename_output=strcat(basename_output, '-2D_');
 195 end
 196
 197 %############# END OF Unchanged parameters ##############
 198
 199 %############# Simulations des donnees ##############
 200 ps=ps/oversamp;
 201
 202 ## Creation of the sinogram
 203 if(strcmp(vers, 'Radon'))
 204         sino_abs=radon(absorption, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
 205         sino_delta=radon(delta, angles)*ps*oversamp;
 206 elseif(strcmp(vers, 'Analytical'))
 207         [sino_abs, sino_delta, sino_beta] = PhantTibo_Analytical (angles, ps*oversamp, 1);
 208 else
 209         disp('error in the variable vers')
 210         return;
 211 endif
 212
 213 ## Creation of the projections
 214 abs_proj=cell(1, nbProj);
 215 B_proj=cell(1, nbProj);
 216 Phi_proj=cell(1, nbProj);
 217
 218
 219 if height ==1
 220
 221         for j=1:nbProj
 222                 abs_proj{j}=interp(sino_abs(:,j), oversamp);
 223                 B_proj{j}=(1/2)*abs_proj{j};
 224                 Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
 225         end
 226
 227         ## Fourier domain
 228         [mp np]=size(Phi_proj{j});
 229         [ftot]=FrequencySpace1(2*mp, ps);
 230
 231         ### CTF model ###
 232         if(model_ctf==1)
 233                 IdCTF=cell(length(D), nbProj);
 234                 sino_IdCTF=cell(length(D),1);
 235
 236                 for k=1:length(D)
 237
 238                         for j=1:nbProj
 239                                 sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 240                         end
 241
 242                         if strcmp(noise_type, 'gaussian')
 243                                 ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
 244                                 bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
 245                                 bruit=bruit/max(max(bruit));
 246                                 bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 247                                 sino_IdCTF{k}=sino_IdCTF{k}+bruit;
 248                         end
 249
 250                         name_ctf=strcat(dir_out, basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k),  '_.edf')
 251                         edfwrite(name_ctf, sino_IdCTF{k}, 'float32');
 252                 end
 253         end % end of model CTF
 254
 255         ### Fresnel model ###
 256         if(model_Fresnel==1)
 257                 IdFresnel=cell(length(D), nbProj);
 258                 sino_IdFresnel=cell(length(D),1);
 259
 260                 for k=1:length(D)
 261                         for j=1:nbProj
 262                                 tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 263                                 IdFresnel{k, j}=tmp;
 264                                 sino_IdFresnel{k}(:,j)=tmp;
 265                         end
 266
 267                         if strcmp(noise_type, 'gaussian')
 268                                 ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
 269                                 bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
 270                                 bruit=bruit/max(max(bruit));
 271                                 bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 272                                 sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
 273                         end
 274
 275                         name_fre=strcat(dir_out, basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
 276                         edfwrite(name_fre, sino_IdFresnel{k}, 'float32');
 277                 end
 278         end % end of Fresnel model
 279
 280 ##### case height >=2 (2D) #####
 281 else
 282
 283         for j=1:nbProj
 284                 abs_proj{j}=repmat(interp(sino_abs(:,j), oversamp), [1 oversamp*height]);
 285                 B_proj{j}=(1/2)*abs_proj{j};
 286                 Phi_proj{j}=repmat( interp(sino_delta(:,j), oversamp), [1 oversamp*height]);
 287                 Phi_proj{j}=(-2*pi/lambda)*Phi_proj{j};
 288         end
 289
 290         ## Fourier domain
 291         [mp np]=size(Phi_proj{j}); % take into account the oversampling
 292         [ftot, gtot]=FrequencySpace(2*mp, 2*np, ps);
 293
 294         ### CTF model ###
 295         if(model_ctf==1)
 296                 IdCTF=cell(length(D), nbProj);
 297
 298                 for k=1:length(D)
 299                         nom=strcat(basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_' );
 300
 301                                 if(exist(strcat(dir_out, nom))!=7)
 302                                         unix(['mkdir ' dir_out nom]); % create directory
 303                                 else
 304                                         unix(['rm ' dir_out nom '/' basename_output 'CTF*' vers '*edf']); % erase existing proj
 305                                 end % end create directory
 306
 307                         for j=1:nbProj
 308                                 IdCTF{k, j}=CTFPropagation_2D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot, gtot);
 309
 310                                 if strcmp(noise_type, 'gaussian')
 311                                         ampsignal=max(max(IdCTF{k, j}))-min(min(IdCTF{k, j})); %add the amount of noise wrt the amplitude at the first dist
 312                                         bruit=randn(size(IdCTF{k, j})); % gaussian white noise
 313                                         bruit=bruit/max(max(bruit));
 314                                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 315                                         IdCTF{k, j}=IdCTF{k, j}+bruit;
 316                                 end
 317
 318                                 name_ctf=strcat(dir_out, nom, '/', nom, num2str(j-1, '%0.4d'), '.edf');
 319                                 edfwrite(name_ctf, IdCTF{k,j}, 'float32');
 320                         end % end proj
 321                 end % end length
 322         end % end ctf
 323
 324         #### Fresnel model ###
 325         if(model_Fresnel==1)
 326                 IdFresnel=cell(length(D), nbProj);
 327
 328                 for k=1:length(D)
 329                         nom=strcat(basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_' );
 330
 331                         if(exist(strcat(dir_out, nom))!=7)
 332                                 unix(['mkdir ' dir_out nom]); % create directory
 333                         else
 334                                 unix(['rm ' dir_out nom '/' basename_output 'Fresnel*' vers '*edf']); % erase existing proj
 335                         end % end create directory
 336
 337                         for j=1:nbProj
 338                                 IdFresnel{k, j}=FresnelTransform_2D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot, gtot);
 339                                 if strcmp(noise_type, 'gaussian')
 340                                         ampsignal=max(max(IdFresnel{k, j}))-min(min(IdFresnel{k, j})); %add the amount of noise wrt the amplitude at the first dist
 341                                         bruit=randn(size(IdFresnel{k, j})); % gaussian white noise
 342                                         bruit=bruit/max(max(bruit));
 343                                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 344                                         IdFresnel{k, j}=IdFresnel{k, j}+bruit;
 345                                 end
 346                                 name_fre=strcat(dir_out, nom, '/', nom, num2str(j-1, '%0.4d'), '.edf');
 347                                 edfwrite(name_fre, IdFresnel{k,j}, 'float32');
 348                         end     % end proj
 349                 end % end distances
 350
 351         end %end Fresnel
 352
 353 end
 354
 355
 356
 357 endfunction