1 ## Copyright (C) 2016 Loriane Weber
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.
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.
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/>.
17 ## SimuPBI_unknown_1D: fonction that launches the phase-contrast simulation of a custom object.
19 ## Author: Loriane Weber <lweber@gpid16a-1802>
20 ## Created: 2016-04-07
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)
25 ## Please note that dir_out, noise_type, noise_amount are optional arguments.
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('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_ESRF')
29 addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_LW')
31 ################################################
32 ############### INPUT parameters ###############
33 ################################################
36 ################################
37 ######### Parameter related to the computation - and names of the result
38 ################################
39 ## vers does not exist anymore; the projections are calculated using the Radon transform
41 ## Oversampling of the projections : use 2 or 4
44 ## Basename of the result file. Is a string.
45 %basename_output='Test-1D';
47 ################################
48 ######### Parameter related to the physics
49 ################################
51 ## Distances of propagation (in m)
54 ## Energy of the incoming X-ray beam (in keV)
57 ## Pixel size of the detector (in um)
60 # Number of projections
63 ## Range of the tomography : 180 or 360 degrees
66 ## Which model do you want to use for the propagation? use 1 or 0
70 ################################
71 ############# Parameters related to the object
72 ################################
74 ## attenuation map (in cm^-1)
75 %attenuation_file='mu_star.edf'
77 ## Refrcative index decrement map (delta)
78 %phase_file='delta_star.edf'
80 ################################
81 ############# OPTIONAL Parameters
82 ################################
85 % dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_CustomObject'
86 % default value if '', i.e. the working directory
89 ## Noise addition to the simulated data?
90 ## Use noise='gaussian' (addition of gaussian noise) or noise='poisson' ( generation of Poisson noise) or '' (no noise is added).
92 %noise_type='gaussian'
93 % default value is '' (no noise)
96 ## If 'gaussian' (additive noise), please specify the Peak-to-peak Signe-to-noise ratio (PPSNR, in dB)
98 %default value is 35 dB
99 ## if 'poisson', please specify the amount of noise
101 % default value is 5%.
103 %noise=0; % O (no noise) or 1 (noise)
104 %PPSNR=24; % noise, in dB
107 #######################################################
108 ############### End of INPUT parameters ###############
109 #######################################################
111 #################################################
112 ############### OUTPUT parameters ###############
113 #################################################
115 ## None, files are directly save
117 ########################################################
118 ############### End of OUTPUT parameters ###############
119 ########################################################
121 ############################# check varargin - optional arguments ###########################
123 error("Some arguments are missing.\n");
130 if strcmp(noise_type, 'gaussian')
132 else strcmp(noise_type, 'poisson')
138 ############################# Unchanged parameters ###########################
139 basename_output=strcat(basename_output, '_');
141 z1=145; % source to sample distance, in meter.
144 angles=[0:1:nbProj-1];
145 delta_angle=range_angle/nbProj;
146 angles=angles*delta_angle;
148 lambda=12.4*10^-10/energy; % in meter
149 ps=ps*10^-6; % put the pixel size in meter
151 D = (z1*dist)./(z1+dist);
152 betash = lambda*D; % in meter^2
155 if !strcmp(dir_out, '')
156 dir_out= strcat(dir_out, '/'); % add underscore at the end
160 % absorption and phase maps
161 mu=edfread(attenuation_file); % in cm-1
162 mu=mu*100; % convert into m-1
164 delta=edfread(phase_file);
168 if(strcmp(noise_type, 'gaussian'))
170 noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
171 elseif(strcmp(noise_type, 'poisson'))
172 scale_fact=noise_amount;
173 noise_str=strcat(noise_type, num2str(scale_fact), '_');
174 elseif(strcmp(noise_type, ''))
177 error("Unknown type of noise")
180 ############################# END OF Unchanged parameters ###########################
182 ############################# Beginning of SIMULATIONS ###################################
185 %%Create the sinograms
186 sino_abs=radon(mu, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
187 sino_delta=radon(delta, angles)*ps*oversamp;
189 %% creation of projections, mu et B, taking into account the oversampling
190 disp('create projections');
191 abs_proj=cell(1, nbProj);
192 B_proj=cell(1, nbProj);
193 Phi_proj=cell(1, nbProj);
196 abs_proj{j}=interp(sino_abs(:,j), oversamp);
197 B_proj{j}=(1/2)*abs_proj{j};
198 Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
204 [mp np]=size(Phi_proj{j});
205 [ftot]=FrequencySpace1(2*mp, ps);
207 ###################### CTF model ################
210 IdCTF=cell(length(D), nbProj);
211 sino_IdCTF=cell(length(D),1);
216 sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
219 if(strcmp(noise_type, 'gaussian'))
220 ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
221 bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
222 bruit=bruit/max(max(bruit));
223 bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
224 sino_IdCTF{k}=sino_IdCTF{k}+bruit;
225 elseif(strcmp(noise_type, 'poisson'))
226 noisy=scale_fact*poissrnd(sino_IdCTF{k}/scale_fact);
230 name_ctf=strcat(basename_output, 'CTF_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
231 edfwrite(strcat(dir_out, name_ctf), sino_IdCTF{k}, 'float32');
233 end % end of model CTF
236 ###################### Fresnel model ################
239 disp('Fresnel model')
241 IdFresnel=cell(length(D), nbProj);
242 sino_IdFresnel=cell(length(D),1);
247 tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
250 sino_IdFresnel{k}(:,j)=tmp;
253 if(strcmp(noise_type, 'gaussian'))
254 ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
255 bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
256 bruit=bruit/max(max(bruit));
257 bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
258 sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
259 elseif(strcmp(noise_type, 'poisson'))
260 noisy=scale_fact*poissrnd(sino_IdFresnel{k}/scale_fact);
261 sino_IdFresnel{k}=noisy;
264 name_fre=strcat(basename_output, 'Fresnel_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
265 edfwrite(name_fre, sino_IdFresnel{k}, 'float32');
267 end % end of Fresnel model