## 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('utilities_ESRF')
addpath('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 ###############################
