SimuPBI_unknown_1D_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_unknown_1D: fonction that launches the phase-contrast simulation of a custom object.
  18
  19 ## Author: Loriane Weber <lweber@gpid16a-1802>
  20 ## Created: 2016-04-07
  21
  22 ## Example of use :
  23 ## SimuPBI_unknown_1D('star', 2, [0 0.1 0.5], 26, 2.7, 300, 180, 1, 1, '/users/lweber/Matlab/SimulationsPBI/mu_star.edf', '/users/lweber/Matlab/SimulationsPBI/delta_star.edf', '', 'gaussian', 35)
  24
  25 ## Please note that dir_out, noise_type, noise_amount are optional arguments.
  26
  27 function [ ret ] = SimuPBI_unknown_1D_func(basename_output, oversamp, dist, energy, ps, nbProj, range_angle, model_ctf, model_Fresnel, attenuation_file, phase_file, dir_out, noise_type, noise_amount)
  28 addpath('utilities_ESRF')
  29 addpath('utilities_LW')
  30 addpath('octave_packages/image-1.0.15')
  31 addpath('octave_packages/signal-1.1.3')
  32
  33 ################################################
  34 ############### INPUT parameters ###############
  35 ################################################
  36
  37
  38 ################################
  39 ######### Parameter related to the computation - and names of the result
  40 ################################
  41 ## vers does not exist anymore; the projections are calculated using the Radon transform
  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.20];
  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 ############# Parameters related to the object
  74 ################################
  75
  76 ## attenuation map (in cm^-1)
  77 %attenuation_file='mu_star.edf'
  78
  79 ## Refrcative index decrement map (delta)
  80 %phase_file='delta_star.edf'
  81
  82 ################################
  83 ############# OPTIONAL Parameters
  84 ################################
  85
  86 ## Output directory
  87 % dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_CustomObject'
  88 % default value if '', i.e. the working directory
  89
  90 ## Noise-related part
  91 ## Noise addition to the simulated data?
  92 ## Use noise='gaussian' (addition of gaussian noise) or noise='poisson' ( generation of Poisson noise) or '' (no noise is added).
  93
  94 %noise_type='gaussian'
  95 % default value is '' (no noise)
  96
  97 ## Amount of noise
  98 ## If 'gaussian' (additive noise), please specify the Peak-to-peak Signe-to-noise ratio (PPSNR, in dB)
  99 % noise_amount=40;
 100 %default value is 35 dB
 101 ## if 'poisson', please specify the amount of noise
 102 %noise_amount=0.05;
 103 % default value is 5%.
 104
 105 %noise=0; % O (no noise) or 1 (noise)
 106 %PPSNR=24; % noise, in dB
 107
 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 ############################# check varargin - optional arguments ###########################
 124  if nargin<11
 125        error("Some arguments are missing.\n");
 126  elseif nargin==11
 127     dir_out='';
 128     noise_type='';
 129  elseif nargin==12
 130         noise_type='';
 131  elseif nargin==13
 132         if strcmp(noise_type, 'gaussian')
 133                 noise_amount=35;
 134         else strcmp(noise_type, 'poisson')
 135         noise_amount=0.05;
 136     endif
 137 endif
 138
 139
 140 ############################# Unchanged parameters ###########################
 141 basename_output=strcat(basename_output, '_');
 142
 143 z1=145; % source to sample distance, in meter.
 144
 145 % calcul des angles
 146 angles=[0:1:nbProj-1];
 147 delta_angle=range_angle/nbProj;
 148 angles=angles*delta_angle;
 149
 150 lambda=12.4*10^-10/energy; % in meter
 151 ps=ps*10^-6; % put the pixel size in meter
 152
 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
 162 % absorption and phase maps
 163 mu=edfread(attenuation_file); % in cm-1
 164 mu=mu*100; % convert into m-1
 165
 166 delta=edfread(phase_file);
 167
 168
 169 ## NOISE
 170 if(strcmp(noise_type, 'gaussian'))
 171         PPSNR=noise_amount;
 172         noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
 173 elseif(strcmp(noise_type, 'poisson'))
 174         scale_fact=noise_amount;
 175         noise_str=strcat(noise_type, num2str(scale_fact), '_');
 176 elseif(strcmp(noise_type, ''))
 177         noise_str='';
 178 else
 179         error("Unknown type of noise")
 180 end
 181
 182 ############################# END OF Unchanged parameters ###########################
 183
 184 ############################# Beginning of SIMULATIONS ###################################
 185 ps=ps/oversamp;
 186
 187 %%Create the sinograms
 188 sino_abs=radon(mu, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
 189 sino_delta=radon(delta, angles)*ps*oversamp;
 190
 191 %% creation of projections, mu et B, taking into account the oversampling
 192 disp('create projections');
 193 abs_proj=cell(1, nbProj);
 194 B_proj=cell(1, nbProj);
 195 Phi_proj=cell(1, nbProj);
 196
 197         for j=1:nbProj
 198                 abs_proj{j}=interp(sino_abs(:,j), oversamp);
 199                 B_proj{j}=(1/2)*abs_proj{j};
 200                 Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
 201         end
 202
 203 clear abs_proj
 204
 205 %% Fourier domain
 206 [mp np]=size(Phi_proj{j});
 207 [ftot]=FrequencySpace1(2*mp, ps);
 208
 209 ###################### CTF model ################
 210 if(model_ctf==1)
 211         disp('CTF model')
 212         IdCTF=cell(length(D), nbProj);
 213         sino_IdCTF=cell(length(D),1);
 214
 215         for k=1:length(D)
 216         k
 217                 for j=1:nbProj
 218                         sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 219                 end
 220
 221                 if(strcmp(noise_type, 'gaussian'))
 222                         ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
 223                         bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
 224                         bruit=bruit/max(max(bruit));
 225                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 226                         sino_IdCTF{k}=sino_IdCTF{k}+bruit;
 227                 elseif(strcmp(noise_type, 'poisson'))
 228                         noisy=scale_fact*poissrnd(sino_IdCTF{k}/scale_fact);
 229                         sino_IdCTF{k}=noisy;
 230                 end
 231
 232                 name_ctf=strcat(basename_output, 'CTF_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k),  '_.edf')
 233                 edfwrite(strcat(dir_out, name_ctf), sino_IdCTF{k}, 'float32');
 234         end
 235 end % end of model CTF
 236
 237
 238 ###################### Fresnel model ################
 239
 240 if(model_Fresnel==1)
 241         disp('Fresnel model')
 242
 243         IdFresnel=cell(length(D), nbProj);
 244         sino_IdFresnel=cell(length(D),1);
 245
 246         for k=1:length(D)
 247                 for j=1:nbProj
 248
 249                         tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
 250
 251                         IdFresnel{k, j}=tmp;
 252                         sino_IdFresnel{k}(:,j)=tmp;
 253                 end
 254
 255                 if(strcmp(noise_type, 'gaussian'))
 256                         ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
 257                         bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
 258                         bruit=bruit/max(max(bruit));
 259                         bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
 260                         sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
 261                 elseif(strcmp(noise_type, 'poisson'))
 262                         noisy=scale_fact*poissrnd(sino_IdFresnel{k}/scale_fact);
 263                         sino_IdFresnel{k}=noisy;
 264                 end
 265
 266                 name_fre=strcat(basename_output, 'Fresnel_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
 267                 edfwrite(name_fre, sino_IdFresnel{k}, 'float32');
 268         end
 269 end % end of Fresnel model
 270
 271
 272 endfunction