From: Loriane Weber Date: Wed, 20 Apr 2016 08:53:21 +0000 (+0200) Subject: First three function for simulations X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=e1e85a737673e93b6b4b01bb0b072cd35dfa3c72;hp=99b9890c11d31f9ae22dd481873e01c708175073;p=CreaPhase.git First three function for simulations ( circle, 3 wires phantom, and custom object in 1D) --- 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 diff --git a/SimuPBI_circle_func.m b/SimuPBI_circle_func.m new file mode 100644 index 0000000..2ddcbfb --- /dev/null +++ b/SimuPBI_circle_func.m @@ -0,0 +1,313 @@ +## 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_circle_func: fonction that launches the phase-contrast simulation of a wire. + +## Author: Loriane Weber +## 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 ############################### diff --git a/SimuPBI_unknown_1D_func.m b/SimuPBI_unknown_1D_func.m new file mode 100644 index 0000000..cac813e --- /dev/null +++ b/SimuPBI_unknown_1D_func.m @@ -0,0 +1,270 @@ +## 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_unknown_1D: fonction that launches the phase-contrast simulation of a custom object. + +## Author: Loriane Weber +## 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