Update function details.

Author: hanyoseob

demo_parallelbeam_ct.m

@@ -37,7 +37,7 @@
 param.dOffsetDctY   = [];
 param.dOffsetDctX   = 20;
 
-param.compute_filtering = 'conv';
+param.compute_filtering = 'fft';
 
 %% Image Object parameters
 % dImgY, dImgX, dImgZ	: Pixel resolution [mm; (float)]

lib/filtering.m

@@ -4,7 +4,7 @@
 % Filtering operator is implemented by both convolution and FFT versions.
 function [pdOut, pdFlt] = filtering(pdIn, param)
 
-pdOut	= zeros(param.nDctX, param.nView);
+pdOut	= zeros(param.nDctX, param.nView, 'like', pdIn);
 
 % Filter is generated based on Ch.3 Equation (3.29)
 pdFlt	= generate_filter(param.dDctX, param.nDctX);
@@ -14,33 +14,30 @@
         % Ch.3 Equation (3.30)
         % CONVOLUTION ver.
         for iview = 0:param.nView-1
-%             pdFltY(:, iview+1) = convolution1d(pdY(:, iview+1), pdFlt, param.nDctX);
+%             pdOut(:, iview+1) = convolution1d(pdIn(:, iview+1), pdFlt, param.nDctX);
             pdOut(:, iview+1) = conv(pdIn(:, iview+1), pdFlt, 'same');
         end
 
     case 'fft'
         % Ch.3 Equation (3.21)
         % FFT ver.
-        pdOut_          = zeros(param.nDctX, 1);
-        
         nFlt            = length(pdFlt);
-        pdFlt_          = pdFlt;
-        pdFlt_(nFlt)    = 0;
-        pdFlt_          = fft(pdFlt_);
+        
+        len            = 2^(ceil(log2(2*param.nDctX)));
+        pdFlt(len)     = 0;
+        pdFlt_          = fft(pdFlt);
 
         for iview = 0:param.nView-1
-            pdOut_(:)  = 0;
 
             pdIn_    	= pdIn(:, iview+1);
-            pdIn_(nFlt)	= 0;
+            pdIn_(len)	= 0;
             pdIn_     	= fft(pdIn_);
 
-            for iflt = 0:nFlt-1
-                pdOut_(iflt+1)	= pdFlt_(iflt+1)*pdIn_(iflt+1);
-            end
+            pdOut_      = pdFlt_.*pdIn_;
 
             pdOut_              = ifft(pdOut_, 'symmetric');
-            pdOut(:, iview+1)	= pdOut_(nFlt - param.nDctX + 1: nFlt);
+%             pdOut(:, iview+1)	= pdOut_(nFlt - param.nDctX + 1: nFlt);
+            pdOut(:, iview+1)	= pdOut_(param.nDctX:2*param.nDctX - 1);
         end
 end