--- /dev/null
+## Copyright (C) 2016 Loriane Weber
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## SimuPBI_3WiresPhant_func
+
+## Author: Loriane Weber <lweber@gpid16a-1802>
+## Created: 2016-03-30
+
+## Example of use :
+## SimuPBI_3WiresPhant_func('Analytical', 2, 'test', [0 0.1 0.5], 19, 2.7, 300, 180, 1, 0, '', 'gaussian', 35)
+
+## Please note that dir_out, noise_type, noise_amount are optional arguments.
+
+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)
+
+addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_ESRF')
+addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_LW')
+
+################################################
+############### INPUT parameters ###############
+################################################
+
+################################
+######### Parameter related to the computation - and names of the result
+################################
+
+## vers: are the projections calculated analytically or using the Radon transform? use 'Radon' OR 'Analytical'
+%vers='Radon' ;
+
+## Oversampling of the projections : use 2 or 4
+%oversamp=2
+
+## Basename of the result file
+%basename_output='phantom_al_mg_PET_thibaut_3dist';
+
+################################
+######### Parameter related to the physics
+################################
+
+## Distances of propagation (in m)
+%dist=[0 0.28 0.560];
+
+## Energy of the incoming X-ray beam (in keV)
+%energy=19
+
+## Pixel size of the detector (in um)
+%ps=3.5;
+
+## Number of projections
+%nbProj=100;
+
+## Range of the tomography : 180 or 360 degrees
+%range_angle=180;
+
+## Which model do you want to use for the propagation? use 1 or 0
+%model_ctf=0;
+%model_Fresnel=1;
+
+################################
+############# Parameters related to the object
+################################
+
+## height of the object. If height == 1, the simulation is computed in 1D. If height >=2 , the simulation turns into 2D.
+%height=10;
+
+
+################################
+############# OPTIONAL Parameters
+################################
+## Output directory : to be removed
+%dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_PhantTibo';
+%dir_out='';
+
+## Noise-related part
+## Noise addition to the simulated data?
+## Use noise='gaussian' (addition of gaussian noise) or noise='poisson' ( generation of Poisson noise) or '' (no noise is added).
+
+%noise_type='gaussian'
+% default value is '' (no noise)
+
+## Amount of noise
+## If 'gaussian' (additive noise), please specify the Peak-to-peak Signe-to-noise ratio (PPSNR, in dB)
+% noise_amount=40;
+%default value is 35 dB
+## if 'poisson', please specify the amount of noise
+%noise_amount=0.05;
+% default value is 5%.
+## !! note that poisson noise is here not implemented yet
+
+#######################################################
+############### End of INPUT parameters ###############
+#######################################################
+
+#################################################
+############### OUTPUT parameters ###############
+#################################################
+
+## None, file directly save
+
+########################################################
+############### End of OUTPUT parameters ###############
+########################################################
+
+############################# check varargin - optional arguments ###########################
+ if nargin<11
+ error("Some arguments are missing.\n");
+ elseif nargin==11
+ dir_out='';
+ noise_type='';
+ elseif nargin==12
+ noise_type='';
+ elseif nargin==13
+ if strcmp(noise_type, 'gaussian')
+ noise_amount=35;
+ elseif strcmp(noise_type, 'poisson')
+ error("poisson noise is not implemented yet")
+ return; %noise_amount=0.05;
+ else
+ error("Unknown type of noise.\n")
+ endif
+endif
+
+
+############# Unchanged parameters ##############
+if !strcmp(dir_out, '')
+ dir_out=strcat(dir_out, '/');
+end
+
+z1=145; % source to sample distance (in m)
+
+% calcul des angles
+angles=[0:1:nbProj-1];
+delta_angle=range_angle/nbProj;
+angles=angles*delta_angle;
+
+% Wavelength, in meter
+lambda=12.4*10^-10/energy;
+ps=ps*10^-6; % put the pixel size in meter
+
+% Effective prop distance
+D = (z1*dist)./(z1+dist);
+betash = lambda*D; % in meter^2
+
+% pading method
+pad_method='extend';
+
+% Absorption image ('mu', in cm-1)
+mu=edfread('attenuation.edf'); % map of mu in cm-1
+[m n]=size(mu)
+
+absorption=mu*100; % in m-1 , to be consistant with lambda!! /10 to have less attenuation
+beta=(lambda/(4*pi))*absorption; % one slice of beta /!\ units
+
+delta_beta=edfread('delta_beta_map_1200.edf');
+delta=delta_beta.*beta;
+
+## Noise addition
+if(strcmp(noise_type, 'gaussian'))
+ PPSNR=noise_amount;
+ noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
+%elseif(strcmp(noise_type, 'poisson'))
+ %scale_fact=noise_amount;
+ %noise_str=strcat(noise_type, num2str(scale_fact), '_');
+elseif(strcmp(noise_type, ''))
+ noise_str='';
+end
+
+
+if height ==1
+ basename_output=strcat(basename_output, '-1D_');
+else
+ basename_output=strcat(basename_output, '-2D_');
+end
+
+%############# END OF Unchanged parameters ##############
+
+%############# Simulations des donnees ##############
+ps=ps/oversamp;
+
+## Creation of the sinogram
+if(strcmp(vers, 'Radon'))
+ sino_abs=radon(absorption, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
+ sino_delta=radon(delta, angles)*ps*oversamp;
+elseif(strcmp(vers, 'Analytical'))
+ [sino_abs, sino_delta, sino_beta] = PhantTibo_Analytical (angles, ps*oversamp, 1);
+else
+ disp('error in the variable vers')
+ return;
+endif
+
+## Creation of the projections
+abs_proj=cell(1, nbProj);
+B_proj=cell(1, nbProj);
+Phi_proj=cell(1, nbProj);
+
+
+if height ==1
+
+ for j=1:nbProj
+ abs_proj{j}=interp(sino_abs(:,j), oversamp);
+ B_proj{j}=(1/2)*abs_proj{j};
+ Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
+ end
+
+ ## Fourier domain
+ [mp np]=size(Phi_proj{j});
+ [ftot]=FrequencySpace1(2*mp, ps);
+
+ ### CTF model ###
+ if(model_ctf==1)
+ IdCTF=cell(length(D), nbProj);
+ sino_IdCTF=cell(length(D),1);
+
+ for k=1:length(D)
+
+ for j=1:nbProj
+ sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
+ end
+
+ if strcmp(noise_type, 'gaussian')
+ ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ sino_IdCTF{k}=sino_IdCTF{k}+bruit;
+ end
+
+ name_ctf=strcat(dir_out, basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
+ edfwrite(name_ctf, sino_IdCTF{k}, 'float32');
+ end
+ end % end of model CTF
+
+ ### Fresnel model ###
+ if(model_Fresnel==1)
+ IdFresnel=cell(length(D), nbProj);
+ sino_IdFresnel=cell(length(D),1);
+
+ for k=1:length(D)
+ for j=1:nbProj
+ tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
+ IdFresnel{k, j}=tmp;
+ sino_IdFresnel{k}(:,j)=tmp;
+ end
+
+ if strcmp(noise_type, 'gaussian')
+ ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
+ end
+
+ name_fre=strcat(dir_out, basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
+ edfwrite(name_fre, sino_IdFresnel{k}, 'float32');
+ end
+ end % end of Fresnel model
+
+##### case height >=2 (2D) #####
+else
+
+ for j=1:nbProj
+ abs_proj{j}=repmat(interp(sino_abs(:,j), oversamp), [1 oversamp*height]);
+ B_proj{j}=(1/2)*abs_proj{j};
+ Phi_proj{j}=repmat( interp(sino_delta(:,j), oversamp), [1 oversamp*height]);
+ Phi_proj{j}=(-2*pi/lambda)*Phi_proj{j};
+ end
+
+ ## Fourier domain
+ [mp np]=size(Phi_proj{j}); % take into account the oversampling
+ [ftot, gtot]=FrequencySpace(2*mp, 2*np, ps);
+
+ ### CTF model ###
+ if(model_ctf==1)
+ IdCTF=cell(length(D), nbProj);
+
+ for k=1:length(D)
+ nom=strcat(basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_' );
+
+ if(exist(strcat(dir_out, nom))!=7)
+ unix(['mkdir ' dir_out nom]); % create directory
+ else
+ unix(['rm ' dir_out nom '/' basename_output 'CTF*' vers '*edf']); % erase existing proj
+ end % end create directory
+
+ for j=1:nbProj
+ IdCTF{k, j}=CTFPropagation_2D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot, gtot);
+
+ if strcmp(noise_type, 'gaussian')
+ ampsignal=max(max(IdCTF{k, j}))-min(min(IdCTF{k, j})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(IdCTF{k, j})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ IdCTF{k, j}=IdCTF{k, j}+bruit;
+ end
+
+ name_ctf=strcat(dir_out, nom, '/', nom, num2str(j-1, '%0.4d'), '.edf');
+ edfwrite(name_ctf, IdCTF{k,j}, 'float32');
+ end % end proj
+ end % end length
+ end % end ctf
+
+ #### Fresnel model ###
+ if(model_Fresnel==1)
+ IdFresnel=cell(length(D), nbProj);
+
+ for k=1:length(D)
+ nom=strcat(basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_' );
+
+ if(exist(strcat(dir_out, nom))!=7)
+ unix(['mkdir ' dir_out nom]); % create directory
+ else
+ unix(['rm ' dir_out nom '/' basename_output 'Fresnel*' vers '*edf']); % erase existing proj
+ end % end create directory
+
+ for j=1:nbProj
+ IdFresnel{k, j}=FresnelTransform_2D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot, gtot);
+ if strcmp(noise_type, 'gaussian')
+ ampsignal=max(max(IdFresnel{k, j}))-min(min(IdFresnel{k, j})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(IdFresnel{k, j})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ IdFresnel{k, j}=IdFresnel{k, j}+bruit;
+ end
+ name_fre=strcat(dir_out, nom, '/', nom, num2str(j-1, '%0.4d'), '.edf');
+ edfwrite(name_fre, IdFresnel{k,j}, 'float32');
+ end % end proj
+ end % end distances
+
+ end %end Fresnel
+
+end
+
+
+
+endfunction
--- /dev/null
+## Copyright (C) 2016 Loriane Weber
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## SimuPBI_circle_func: fonction that launches the phase-contrast simulation of a wire.
+
+## Author: Loriane Weber <lweber@gpid16a-1802>
+## Created: 2016-02-11
+
+## Example of use :
+## SimuPBI_circle_func('Analytical', 2, 'test', [0 0.1 0.5], 19, 2.7, 300, 180, 1, 1, 0.89, 8.265*10^-7, 100, 500, [200 250], '', 'gaussian', 35)
+
+## Please note that dir_out, noise_type, noise_amount are optional arguments.
+
+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)
+
+addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_ESRF')
+addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_LW')
+
+################################################
+############### INPUT parameters ###############
+################################################
+
+
+################################
+######### Parameter related to the computation - and names of the result
+################################
+## vers: are the projections calculated analytically or using the Radon transform? use the strings 'Radon' OR 'Analytical'
+%vers='Analytical'
+
+## Oversampling of the projections : use 2 or 4
+%oversamp=2
+
+## Basename of the result file. Is a string.
+%basename_output='Test-1D';
+
+################################
+######### Parameter related to the physics
+################################
+
+## Distances of propagation (in m)
+%dist=[0 0.01 0.1 0.20 0.50];
+
+## Energy of the incoming X-ray beam (in keV)
+%energy=19
+
+## Pixel size of the detector (in um)
+%ps=1;
+
+# Number of projections
+%nbProj=360
+
+## Range of the tomography : 180 or 360 degrees
+%range_angle=180;
+
+## Which model do you want to use for the propagation? use 1 or 0
+%model_ctf=1;
+%model_Fresnel=1;
+
+
+################################
+############# Parameters related to the circle
+################################
+
+## material of the circle (here, PET at 19 keV)
+%mu_mat=0.89889 % (in cm-1);
+%delta_mat=8.265*10^-7;
+
+## size and coordinates
+%R_circle=50; % radius of the circle
+%size_image=500; % size of the image (square)
+%circle_center=[200; 250]; % 2D coordinates of the center
+
+################################
+############# OPTIONAL Parameters
+################################
+
+## Output directory
+%dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_Circle'
+% default value if '', i.e. the working directory
+
+## Noise-related part
+## Noise addition to the simulated data?
+## Use noise='gaussian' (addition of gaussian noise) or noise='poisson' ( generation of Poisson noise) or '' (no noise is added).
+
+%noise_type='gaussian'
+% default value is '' (no noise)
+
+## Amount of noise
+## If 'gaussian' (additive noise), please specify the Peak-to-peak Signe-to-noise ratio (PPSNR, in dB)
+% noise_amount=40;
+%default value is 35 dB
+## if 'poisson', please specify the amount of noise
+%noise_amount=0.05;
+% default value is 5%.
+
+#######################################################
+############### End of INPUT parameters ###############
+#######################################################
+
+#################################################
+############### OUTPUT parameters ###############
+#################################################
+
+## None, files are directly save
+
+########################################################
+############### End of OUTPUT parameters ###############
+########################################################
+
+
+############################# check varargin - optional arguments ###########################
+ if nargin<15
+ error("Some arguments are missing.\n");
+ elseif nargin==15
+ dir_out='';
+ noise_type='';
+ elseif nargin==16
+ noise_type='';
+ elseif nargin==17
+ if strcmp(noise_type, 'gaussian')
+ noise_amount=35;
+ else strcmp(noise_type, 'poisson')
+ noise_amount=0.05;
+ endif
+endif
+
+
+############################# Unchanged parameters ###########################
+basename_output=strcat(basename_output, '_');
+
+% calcul des angles
+angles=[0:1:nbProj-1];
+delta_angle=range_angle/nbProj;
+angles=angles*delta_angle;
+
+lambda=12.4*10^-10/energy; % in meter
+ps=ps*10^-6; % put the pixel size in meter
+
+z1=145;
+D = (z1*dist)./(z1+dist);
+betash = lambda*D; % in meter^2
+pad_method='extend';
+
+if !strcmp(dir_out, '')
+ dir_out= strcat(dir_out, '/'); % add underscore at the end
+end
+
+% absorption and phase parameter
+mu_mat=mu_mat*100; % in m-1 , to be consistant with lambda!!
+beta_mat=lambda*mu_mat/(4*pi);
+
+%% Cartes 2D d'attenuation et d/b
+mu=zeros(size_image);
+for(i=1:size_image)
+ for(j=1:size_image)
+ if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
+ mu(i,j)=mu_mat;
+ end
+ end
+end
+namnam=strcat(dir_out, 'attenuationCirclePET.edf')
+edfwrite(namnam, mu, 'float32');
+
+
+delta=zeros(size_image);
+for(i=1:size_image)
+ for(j=1:size_image)
+ if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
+ delta(i,j)=delta_mat;
+ end
+ end
+end
+namnam=strcat(dir_out, 'deltaCirclePET.edf')
+edfwrite(namnam, delta, 'float32');
+
+beta=zeros(size_image);
+for(i=1:size_image)
+ for(j=1:size_image)
+ if( ((i-circle_center(1))^2 + (j-circle_center(2))^2) <= (R_circle^2) )
+ beta(i,j)=beta_mat;
+ end
+ end
+end
+namnam=strcat(dir_out, 'betaCirclePET.edf')
+edfwrite(namnam, beta, 'float32');
+
+
+if(strcmp(noise_type, 'gaussian'))
+ PPSNR=noise_amount;
+ noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
+elseif(strcmp(noise_type, 'poisson'))
+ scale_fact=noise_amount;
+ noise_str=strcat(noise_type, num2str(scale_fact), '_');
+elseif(strcmp(noise_type, ''))
+ noise_str='';
+else
+ error("Unknown type of noise")
+end
+
+############################# END OF Unchanged parameters ###########################
+
+############################# Beginning of simulation ###############################
+ps=ps/oversamp;
+
+ % Create the sinograms
+disp('Create sinograms');
+ if(strcmp(vers, 'Radon'))
+ sino_abs=radon(mu, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
+ sino_delta=radon(delta, angles)*ps*oversamp;
+ sino_beta=radon(beta, angles)*ps*oversamp;
+
+ elseif(strcmp(vers, 'Analytical'))
+ [sino_abs, sino_delta, sino_beta] = CircleMat_Analytical (angles, ps*oversamp, R_circle, circle_center, mu_mat, beta_mat, delta_mat, size_image, 0);
+
+ else
+ disp('error in the variable vers')
+ return;
+ endif
+
+ % creation of projections, mu et B, taking into account the oversampling
+disp('create projections');
+abs_proj=cell(1, nbProj);
+B_proj=cell(1, nbProj);
+Phi_proj=cell(1, nbProj);
+
+ for j=1:nbProj
+ abs_proj{j}=interp(sino_abs(:,j), oversamp);
+ B_proj{j}=(1/2)*abs_proj{j};
+ Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
+ end
+
+% Fourier domain
+[mp np]=size(Phi_proj{j});
+[ftot]=FrequencySpace1(2*mp, ps);
+
+
+###################### CTF model ################
+if(model_ctf==1)
+ disp('CTF model')
+ IdCTF=cell(length(D), nbProj);
+ sino_IdCTF=cell(length(D),1);
+
+ for k=1:length(D)
+ k
+
+ for j=1:nbProj
+ sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
+ end
+
+ if(strcmp(noise_type, 'gaussian'))
+ ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ sino_IdCTF{k}=sino_IdCTF{k}+bruit;
+ elseif(strcmp(noise_type, 'poisson'))
+ noisy=scale_fact*poissrnd(sino_IdCTF{k}/scale_fact);
+ sino_IdCTF{k}=noisy;
+ end
+
+ name_ctf=strcat(basename_output, 'CTF_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
+ edfwrite(strcat(dir_out, name_ctf), sino_IdCTF{k}, 'float32');
+ end
+end % end of model CTF
+
+
+###################### Fresnel model ################
+
+if(model_Fresnel==1)
+ disp('Fresnel model')
+
+ IdFresnel=cell(length(D), nbProj);
+ sino_IdFresnel=cell(length(D),1);
+
+ for k=1:length(D)
+ for j=1:nbProj
+
+ tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
+
+ IdFresnel{k, j}=tmp;
+ sino_IdFresnel{k}(:,j)=tmp;
+ end
+
+ if(strcmp(noise_type, 'gaussian'))
+ ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
+ elseif(strcmp(noise_type, 'poisson'))
+ noisy=scale_fact*poissrnd(sino_IdFresnel{k}/scale_fact);
+ sino_IdFresnel{k}=noisy;
+ end
+
+ name_fre=strcat(basename_output, 'Fresnel_', vers, '_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
+ edfwrite(strcat(dir_out, name_fre), sino_IdFresnel{k}, 'float32');
+ end
+end % end of Fresnel model
+
+############################# End of simulation ###############################
--- /dev/null
+## Copyright (C) 2016 Loriane Weber
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## SimuPBI_unknown_1D: fonction that launches the phase-contrast simulation of a custom object.
+
+## Author: Loriane Weber <lweber@gpid16a-1802>
+## Created: 2016-04-07
+
+## Example of use :
+## 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)
+
+## Please note that dir_out, noise_type, noise_amount are optional arguments.
+
+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)
+addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_ESRF')
+addpath('/users/lweber/Matlab/SimulationsPBI/CreaPhase/utilities_LW')
+
+################################################
+############### INPUT parameters ###############
+################################################
+
+
+################################
+######### Parameter related to the computation - and names of the result
+################################
+## vers does not exist anymore; the projections are calculated using the Radon transform
+
+## Oversampling of the projections : use 2 or 4
+%oversamp=2
+
+## Basename of the result file. Is a string.
+%basename_output='Test-1D';
+
+################################
+######### Parameter related to the physics
+################################
+
+## Distances of propagation (in m)
+%dist=[0 0.01 0.20];
+
+## Energy of the incoming X-ray beam (in keV)
+%energy=19
+
+## Pixel size of the detector (in um)
+%ps=1;
+
+# Number of projections
+%nbProj=360
+
+## Range of the tomography : 180 or 360 degrees
+%range_angle=180;
+
+## Which model do you want to use for the propagation? use 1 or 0
+%model_ctf=1;
+%model_Fresnel=1;
+
+################################
+############# Parameters related to the object
+################################
+
+## attenuation map (in cm^-1)
+%attenuation_file='mu_star.edf'
+
+## Refrcative index decrement map (delta)
+%phase_file='delta_star.edf'
+
+################################
+############# OPTIONAL Parameters
+################################
+
+## Output directory
+% dir_out='/mntdirect/_users/lweber/Matlab/SimulationsPBI/Results_CustomObject'
+% default value if '', i.e. the working directory
+
+## Noise-related part
+## Noise addition to the simulated data?
+## Use noise='gaussian' (addition of gaussian noise) or noise='poisson' ( generation of Poisson noise) or '' (no noise is added).
+
+%noise_type='gaussian'
+% default value is '' (no noise)
+
+## Amount of noise
+## If 'gaussian' (additive noise), please specify the Peak-to-peak Signe-to-noise ratio (PPSNR, in dB)
+% noise_amount=40;
+%default value is 35 dB
+## if 'poisson', please specify the amount of noise
+%noise_amount=0.05;
+% default value is 5%.
+
+%noise=0; % O (no noise) or 1 (noise)
+%PPSNR=24; % noise, in dB
+
+
+#######################################################
+############### End of INPUT parameters ###############
+#######################################################
+
+#################################################
+############### OUTPUT parameters ###############
+#################################################
+
+## None, files are directly save
+
+########################################################
+############### End of OUTPUT parameters ###############
+########################################################
+
+############################# check varargin - optional arguments ###########################
+ if nargin<11
+ error("Some arguments are missing.\n");
+ elseif nargin==11
+ dir_out='';
+ noise_type='';
+ elseif nargin==12
+ noise_type='';
+ elseif nargin==13
+ if strcmp(noise_type, 'gaussian')
+ noise_amount=35;
+ else strcmp(noise_type, 'poisson')
+ noise_amount=0.05;
+ endif
+endif
+
+
+############################# Unchanged parameters ###########################
+basename_output=strcat(basename_output, '_');
+
+z1=145; % source to sample distance, in meter.
+
+% calcul des angles
+angles=[0:1:nbProj-1];
+delta_angle=range_angle/nbProj;
+angles=angles*delta_angle;
+
+lambda=12.4*10^-10/energy; % in meter
+ps=ps*10^-6; % put the pixel size in meter
+
+D = (z1*dist)./(z1+dist);
+betash = lambda*D; % in meter^2
+pad_method='extend';
+
+if !strcmp(dir_out, '')
+ dir_out= strcat(dir_out, '/'); % add underscore at the end
+end
+
+
+% absorption and phase maps
+mu=edfread(attenuation_file); % in cm-1
+mu=mu*100; % convert into m-1
+
+delta=edfread(phase_file);
+
+
+## NOISE
+if(strcmp(noise_type, 'gaussian'))
+ PPSNR=noise_amount;
+ noise_str=strcat(noise_type, num2str(PPSNR), 'dB_');
+elseif(strcmp(noise_type, 'poisson'))
+ scale_fact=noise_amount;
+ noise_str=strcat(noise_type, num2str(scale_fact), '_');
+elseif(strcmp(noise_type, ''))
+ noise_str='';
+else
+ error("Unknown type of noise")
+end
+
+############################# END OF Unchanged parameters ###########################
+
+############################# Beginning of SIMULATIONS ###################################
+ps=ps/oversamp;
+
+%%Create the sinograms
+sino_abs=radon(mu, angles)*ps*oversamp; % to get rif of oversamp in the pixelsize
+sino_delta=radon(delta, angles)*ps*oversamp;
+
+%% creation of projections, mu et B, taking into account the oversampling
+disp('create projections');
+abs_proj=cell(1, nbProj);
+B_proj=cell(1, nbProj);
+Phi_proj=cell(1, nbProj);
+
+ for j=1:nbProj
+ abs_proj{j}=interp(sino_abs(:,j), oversamp);
+ B_proj{j}=(1/2)*abs_proj{j};
+ Phi_proj{j}=(-2*pi/lambda)*interp(sino_delta(:,j), oversamp);
+ end
+
+clear abs_proj
+
+%% Fourier domain
+[mp np]=size(Phi_proj{j});
+[ftot]=FrequencySpace1(2*mp, ps);
+
+###################### CTF model ################
+if(model_ctf==1)
+ disp('CTF model')
+ IdCTF=cell(length(D), nbProj);
+ sino_IdCTF=cell(length(D),1);
+
+ for k=1:length(D)
+ k
+ for j=1:nbProj
+ sino_IdCTF{k}(:,j)=CTFPropagation_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
+ end
+
+ if(strcmp(noise_type, 'gaussian'))
+ ampsignal=max(max(sino_IdCTF{k}))-min(min(sino_IdCTF{k})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(sino_IdCTF{k})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ sino_IdCTF{k}=sino_IdCTF{k}+bruit;
+ elseif(strcmp(noise_type, 'poisson'))
+ noisy=scale_fact*poissrnd(sino_IdCTF{k}/scale_fact);
+ sino_IdCTF{k}=noisy;
+ end
+
+ name_ctf=strcat(basename_output, 'CTF_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
+ edfwrite(strcat(dir_out, name_ctf), sino_IdCTF{k}, 'float32');
+ end
+end % end of model CTF
+
+
+###################### Fresnel model ################
+
+if(model_Fresnel==1)
+ disp('Fresnel model')
+
+ IdFresnel=cell(length(D), nbProj);
+ sino_IdFresnel=cell(length(D),1);
+
+ for k=1:length(D)
+ for j=1:nbProj
+
+ tmp=FresnelTransform_1D(Phi_proj{j}, B_proj{j}, D(k), lambda, oversamp, ftot);
+
+ IdFresnel{k, j}=tmp;
+ sino_IdFresnel{k}(:,j)=tmp;
+ end
+
+ if(strcmp(noise_type, 'gaussian'))
+ ampsignal=max(max(sino_IdFresnel{k}))-min(min(sino_IdFresnel{k})); %add the amount of noise wrt the amplitude at the first dist
+ bruit=randn(size(sino_IdFresnel{k})); % gaussian white noise
+ bruit=bruit/max(max(bruit));
+ bruit=bruit*ampsignal*10^(-PPSNR/20); %/2;
+ sino_IdFresnel{k}=sino_IdFresnel{k}+bruit;
+ elseif(strcmp(noise_type, 'poisson'))
+ noisy=scale_fact*poissrnd(sino_IdFresnel{k}/scale_fact);
+ sino_IdFresnel{k}=noisy;
+ end
+
+ name_fre=strcat(basename_output, 'Fresnel_oversamp', num2str(oversamp), '_', noise_str, num2str(nbProj), 'proj_', num2str(k), '_.edf')
+ edfwrite(name_fre, sino_IdFresnel{k}, 'float32');
+ end
+end % end of Fresnel model
+
+
+endfunction
git://git.creatis.insa-lyon.fr / CreaPhase.git / commitdiff