X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;ds=sidebyside;f=utilities_LW%2FCTFPropagation_2D.m;fp=utilities_LW%2FCTFPropagation_2D.m;h=f307e63574d48895728212a6c88e72454e069360;hb=d0401c49b6b7511cfdaa0534b78bd3c5b2b0637a;hp=0000000000000000000000000000000000000000;hpb=f8358f5ec65f099f3080043580ef861c3fd3ba2e;p=CreaPhase.git diff --git a/utilities_LW/CTFPropagation_2D.m b/utilities_LW/CTFPropagation_2D.m new file mode 100644 index 0000000..f307e63 --- /dev/null +++ b/utilities_LW/CTFPropagation_2D.m @@ -0,0 +1,82 @@ +## 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 +## . + +## CTFPropagation_2D + +## Author: Loriane Weber +## 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