utilities_LW/CTFPropagation_2D.m

   1 ## Copyright (C) 2015 Loriane Weber
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with Octave; see the file COPYING.  If not, see
  15 ## <http://www.gnu.org/licenses/>.
  16
  17 ## CTFPropagation_2D
  18
  19 ## Author: Loriane Weber <lweber@gpid16a-1802>
  20 ## Created: 2015-12-11
  21
  22 function [ IdCTF ] = CTFPropagation_2D (ph, B, D, lambda, oversamp, ftot, gtot)
  23
  24
  25 [mp, np]=size(ph);
  26 [mp_a, np_a] = size(B);
  27
  28 % check dimensions
  29         if (mp_a != mp || np_a != np)
  30                 disp('error in the dimensions !! Phi and B should have the same dimension')
  31         return;
  32         end
  33
  34 % propagators in 2D (sin and cos)
  35 ssq=2*ssquare2(lambda*D, ftot', gtot', 1);
  36 csq=2*csquare2(lambda*D, ftot', gtot', 1);
  37 %keyboard
  38 % initialisation
  39 tmp_B=zeros(2*mp, 2*np);
  40 tmp_Phi=zeros(2*mp, 2*np);
  41
  42 %tmp_B(1:mp,1:np)=B;
  43 %tmp_ph(1:mp,1:np)=ph;
  44
  45 tmp_B=im_pad(B, 2*np, 2*mp, "zeros");
  46 tmp_Phi=im_pad(ph, 2*np, 2*mp, "zeros");
  47
  48 tmp_B=fft2(tmp_B);
  49 tmp_Phi=fft2(tmp_Phi);
  50
  51 tmp_Id=-csq.*tmp_B + ssq.*tmp_Phi;
  52 tmp_Id=1+real(ifft2(tmp_Id));
  53
  54 IdCTF=cut(tmp_Id, mp, np);
  55
  56         if oversamp==2
  57
  58                 ipf = [ 0.25 0.5 0.25; % modelise le binning du detecteur
  59                                 0.50 1.0 0.50;
  60                                 0.25 0.5 0.25];
  61                 Idd = conv2(IdCTF,ipf,'same')./conv2(ones(size(IdCTF)),ipf,'same');
  62                 Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end);
  63                 IdCTF=Idd;
  64
  65         elseif oversamp==4
  66                 ipf = [0.25 0.50 0.50 0.50 0.25;
  67                            0.50 1.00 1.00 1.00 0.50;
  68                            0.50 1.00 1.00 1.00 0.50;
  69                            0.50 1.00 1.00 1.00 0.50;
  70                            0.25 0.50 0.50 0.50 0.25];
  71
  72                 Idd = conv2(IdCTF,ipf,'same')./conv2(ones(size(IdCTF)),ipf,'same');
  73                 Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end);
  74                 IdCTF=Idd;
  75
  76         end % end oversampling
  77
  78
  79
  80
  81
  82 endfunction