clean up ringremove

functions/tomography/ringremove.m

@@ -6,8 +6,7 @@ function sinoCorrected = ringremove(sino,settings)
 % Removes ring artifacts in tomographic reconstructions by reducing the amount of
 % horizontal lines in the given sinogram. Ring removal is either based on the method
 % proposed by Ketcham et al. :cite:`Ketcham_PS_2006`, denoted as 'additive', or on the method proposed
-% by Muench et al. :cite:`Muench_OE_2009`, denoted as 'wavelet'. Large parts of the code for the
-% wavelet-based approach are taken from the original paper
+% by Muench et al. :cite:`Muench_OE_2009`, denoted as 'wavelet'.
 %
 % Parameters
 % ----------
@@ -38,7 +37,7 @@ function sinoCorrected = ringremove(sino,settings)
 %     Highest decomposition level in the wavelet decomposition when using method
 %     'wavelet'.
 % wavelet_Sigma : Default = 1
-%     Standard deviation for the Gaussian filtering of the horizontal part of the
+%     Width of the Gaussian filter kernel for the horizontal part of the
 %     sinogram in Fourier space when using method 'wavelet'.
 %
 % Returns
@@ -123,7 +122,6 @@ settings = completeStruct(settings, defaults);
 
 switch settings.method
     case 'additive'
-        numAngles = size(sino,2);
         % averaging of sinogram along the angle
         sinoMean = mean(sino,2);
         
@@ -140,50 +138,43 @@ switch settings.method
         end
         % identification of horizontal stripes
         correctionProfile = sinoMean - sinoMeanSmoothed;
-        % creation of 2d correction mask
-        correction2D = repmat(correctionProfile,[1 numAngles]);
         % sinogram correction
-        sinoCorrected = sino - settings.additive_Strength*correction2D;
+        sinoCorrected = sino - settings.additive_Strength * correctionProfile;
        
     case 'wavelet'
-        numAngles = size(sino,2);
-        sizeProjs = size(sino,1);
-        
-        % wavelet decomposition
-        for decompIdx = 1:settings.wavelet_DecNum
-            [sino,horzPart{decompIdx},vertPart{decompIdx},diagPart{decompIdx}] = ...
-                dwt2(sino,settings.wavelet_Wname);
+		
+        % perform multilevel wavelet decomposition		
+		H = cell(settings.wavelet_DecNum, 1);
+		V = cell(settings.wavelet_DecNum, 1);
+		D = cell(settings.wavelet_DecNum, 1);
+		S = cell(settings.wavelet_DecNum, 1);
+				
+        for l = 1:settings.wavelet_DecNum
+			S{l} = size(sino);
+            [sino, H{l}, V{l}, D{l}] = dwt2(sino,settings.wavelet_Wname);
         end
         
-        % FFT transform of horizontal frequency bands
-        for decompIdx = 1:settings.wavelet_DecNum
+        % filter horizontal frequency bands
+        for l = 1:settings.wavelet_DecNum
             % FFT
-            horzFourier = fftshift(fft(horzPart{decompIdx}'));
-            [my,mx] = size(horzFourier);
+            horzFourier = fft(H{l}');
             
             % damping of horizontal stripe information
-            dampGauss = 1 - exp(-(-floor(my/2):-floor(my/2) + my - 1).^2/(2*settings.wavelet_Sigma^2));
-            horzFourier = horzFourier.*repmat(dampGauss',1,mx);
+			m = size(horzFourier,1);
+            k = fftfreq(m, 2*pi/m);
+            dampGauss = 1 - exp(-1/2 * k.^2 / settings.wavelet_Sigma^2 );
+            horzFourier = horzFourier .* dampGauss;
             
             % inverse FFT
-            horzPart{decompIdx} = ifft(ifftshift(horzFourier))';
-        end
-        
-        % wavelet reconstruction
-        sinoCorrected = sino;
-        for decompIdx = settings.wavelet_DecNum:-1:1
-            sinoCorrected = sinoCorrected(1:size(vertPart{decompIdx},1),1:size(vertPart{decompIdx},2));
-            sinoCorrected = idwt2(sinoCorrected,horzPart{decompIdx},vertPart{decompIdx},...
-                diagPart{decompIdx},settings.wavelet_Wname);
+            H{l} = ifft(horzFourier)';
         end
         
-        % make sure that corrected sinogram has the same size as the original sinogram
-        if size(sinoCorrected,2) == numAngles + 1
-            sinoCorrected = sinoCorrected(:,1:end - 1);
-        end
-        if size(sinoCorrected,1) == sizeProjs + 1
-            sinoCorrected = sinoCorrected(1:end - 1,:);
-        end
+		% wavelet reconstruction
+		sinoCorrected = sino;
+		for l = settings.wavelet_DecNum:-1:1
+			sinoCorrected = sinoCorrected(1:size(V{l},1),1:size(V{l},2));
+			sinoCorrected = idwt2(sinoCorrected,H{l}, V{l}, D{l},settings.wavelet_Wname, S{l});
+		end	
         
     otherwise
         error('Choose between ringremove based on "additive" or "wavelet" correction for setting "method"');