--- /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 / blobdiff