--- /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 ###############################
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff