## Copyright (C) 2015 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/>.

## CTFPropagation_2D

## Author: Loriane Weber <lweber@gpid16a-1802>
## Created: 2015-12-11

function [ IdCTF ] = CTFPropagation_2D (ph, B, D, lambda, oversamp, ftot, gtot)


[mp, np]=size(ph);
[mp_a, np_a] = size(B);

% check dimensions 
	if (mp_a != mp || np_a != np)
		disp('error in the dimensions !! Phi and B should have the same dimension')
	return;
	end

% propagators in 2D (sin and cos)
ssq=2*ssquare2(lambda*D, ftot', gtot', 1);
csq=2*csquare2(lambda*D, ftot', gtot', 1);
%keyboard
% initialisation
tmp_B=zeros(2*mp, 2*np);
tmp_Phi=zeros(2*mp, 2*np);

%tmp_B(1:mp,1:np)=B;
%tmp_ph(1:mp,1:np)=ph;

tmp_B=im_pad(B, 2*np, 2*mp, "zeros");
tmp_Phi=im_pad(ph, 2*np, 2*mp, "zeros");
		
tmp_B=fft2(tmp_B);
tmp_Phi=fft2(tmp_Phi);

tmp_Id=-csq.*tmp_B + ssq.*tmp_Phi;
tmp_Id=1+real(ifft2(tmp_Id)); 

IdCTF=cut(tmp_Id, mp, np);

	if oversamp==2
	
		ipf = [ 0.25 0.5 0.25; % modelise le binning du detecteur
				0.50 1.0 0.50;
				0.25 0.5 0.25];
		Idd = conv2(IdCTF,ipf,'same')./conv2(ones(size(IdCTF)),ipf,'same');
		Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end);
		IdCTF=Idd; 
		
	elseif oversamp==4
		ipf = [0.25 0.50 0.50 0.50 0.25;
			   0.50 1.00 1.00 1.00 0.50;
			   0.50 1.00 1.00 1.00 0.50;
			   0.50 1.00 1.00 1.00 0.50;
			   0.25 0.50 0.50 0.50 0.25];
			   
		Idd = conv2(IdCTF,ipf,'same')./conv2(ones(size(IdCTF)),ipf,'same');
		Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end);
		IdCTF=Idd;  
			
	end % end oversampling
			




endfunction
