X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=SimuPBI_3WiresPhant_func.m;fp=SimuPBI_3WiresPhant_func.m;h=46b7eaf0bf53002f0cf4678277017481b7e3ca69;hp=0000000000000000000000000000000000000000;hb=e1e85a737673e93b6b4b01bb0b072cd35dfa3c72;hpb=99b9890c11d31f9ae22dd481873e01c708175073
diff --git a/SimuPBI_3WiresPhant_func.m b/SimuPBI_3WiresPhant_func.m
new file mode 100644
index 0000000..46b7eaf
--- /dev/null
+++ b/SimuPBI_3WiresPhant_func.m
@@ -0,0 +1,348 @@
+## 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
+## .
+
+## SimuPBI_3WiresPhant_func
+
+## Author: Loriane Weber
+## 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