X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=utilities_LW%2FCTFPropagation_1D.m;fp=utilities_LW%2FCTFPropagation_1D.m;h=2801aab68bfc12097a1b7f134fddb9a118b8e915;hb=d0401c49b6b7511cfdaa0534b78bd3c5b2b0637a;hp=0000000000000000000000000000000000000000;hpb=f8358f5ec65f099f3080043580ef861c3fd3ba2e;p=CreaPhase.git diff --git a/utilities_LW/CTFPropagation_1D.m b/utilities_LW/CTFPropagation_1D.m new file mode 100644 index 0000000..2801aab --- /dev/null +++ b/utilities_LW/CTFPropagation_1D.m @@ -0,0 +1,65 @@ +## 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 1D + +## Author: Loriane Weber +## Created: 2015-11-20 + +function [ Id ] = CTFPropagation_1D (ph, B, D, lambda, oversamp, ftot) + +[mp, np]=size(ph); + +if(np!=1) + disp('error') + return +end + +%% propagators in 1D (sin and cos) +ssq=2*ssquare1(lambda*D, ftot, 1)'; +csq=2*csquare1(lambda*D, ftot, 1)'; + +tmp_B=zeros(2*mp, 1); +tmp_ph=zeros(2*mp, 1); +tmp_B(1:mp,1)=B; +tmp_ph(1:mp,1)=ph; + +% fourier domain +tmp_B=fft(tmp_B); +tmp_ph=fft(tmp_ph); + +tmp_Id=-csq.*tmp_B + ssq.*tmp_ph; +tmp_Id=1+real(ifft(tmp_Id)); + + +Id=tmp_Id(1:mp, 1); + + if oversamp==2 + ipf = [0.5 1 0.5]; % modelise le binning du detecteur + Idd = conv(Id,ipf,'same')./conv(ones(size(Id)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, :); + Id=Idd; + elseif oversamp==4 + ipf = [0.5 1.0 1.0 1.0 0.5]; % modelise le binning du detecteur + Idd = conv(Id,ipf,'same')./conv(ones(size(Id)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, :); + + + Id=Idd; + end % oversamp + + +endfunction