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

## Author: Loriane Weber <lweber@gpid16a-1802>
## 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('utilities_ESRF')
addpath('utilities_LW')
addpath(genpath('octave_packages'))

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