X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=SimuPBI_unknown_1D_func.m;fp=SimuPBI_unknown_1D_func.m;h=cac813ea3f4db15fb25a72a297eed426cf49d711;hb=e1e85a737673e93b6b4b01bb0b072cd35dfa3c72;hp=0000000000000000000000000000000000000000;hpb=99b9890c11d31f9ae22dd481873e01c708175073;p=CreaPhase.git

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
+## <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