Merge pull request #186 from holtzhau/radon_fixes

BUG: Fixed padding in radon, iradon. Tests for small images.

skimage/transform/radon_transform.py

@@ -44,8 +44,8 @@ def radon(image, theta=None):
         theta = np.arange(180)    
     height, width = image.shape
     diagonal = np.sqrt(height ** 2 + width ** 2)
-    heightpad = np.ceil(diagonal - height) + 2
-    widthpad = np.ceil(diagonal - width) + 2
+    heightpad = np.ceil(diagonal - height)
+    widthpad = np.ceil(diagonal - width)
     padded_image = np.zeros((int(height + heightpad),
                              int(width + widthpad)))
     y0, y1 = int(np.ceil(heightpad / 2)), \
@@ -57,14 +57,16 @@ def radon(image, theta=None):
     out = np.zeros((max(padded_image.shape), len(theta)))
 
     h, w = padded_image.shape
-    shift0 = np.array([[1, 0, -w/2.],
-                       [0, 1, -h/2.],
+    dh, dw = h / 2, w / 2
+    shift0 = np.array([[1, 0, -dw],
+                       [0, 1, -dh],
                        [0, 0, 1]])
 
-    shift1 = np.array([[1, 0, w/2.],
-                       [0, 1, h/2.],
+    shift1 = np.array([[1, 0, dw],
+                       [0, 1, dh],
                        [0, 0, 1]])
 
+
     def build_rotation(theta):
         T = -np.deg2rad(theta)
 
@@ -129,7 +131,7 @@ def iradon(radon_image, theta=None, output_size=None,
     th = (np.pi / 180.0) * theta        
     # if output size not specified, estimate from input radon image
     if not output_size:
-        output_size = 2 * np.floor(radon_image.shape[0] / (2 * np.sqrt(2)))
+        output_size = int(np.floor(np.sqrt((radon_image.shape[0]) ** 2 / 2.0)))
     n = radon_image.shape[0]
   
     img = radon_image.copy()
@@ -166,13 +168,14 @@ def iradon(radon_image, theta=None, output_size=None,
     # resize filtered image back to original size
     radon_filtered = radon_filtered[:radon_image.shape[0], :]
     reconstructed = np.zeros((output_size, output_size))
-    mid_index = np.ceil(n/2);
+    mid_index = np.ceil(n / 2.0)
+    
     x = output_size
     y = output_size
     [X, Y] = np.mgrid[0.0:x, 0.0:y]
-    xpr = X - (output_size + 1.0) / 2.0
-    ypr = Y - (output_size + 1.0) / 2.0
-    
+    xpr = X - int(output_size) / 2
+    ypr = Y - int(output_size) / 2 
+
     # reconstruct image by interpolation
     if interpolation == "nearest":   
         for i in range(len(theta)):

skimage/transform/tests/test_radon_transform.py

@@ -10,6 +10,7 @@ def rescale(x):
 
 def test_radon_iradon():
     size = 100
+    debug = False
     image = np.tri(size) + np.tri(size)[::-1]
     for filter_type in ["ramp", "shepp-logan", "cosine", "hamming", "hann"]:
         reconstructed = iradon(radon(image), filter=filter_type)
@@ -18,12 +19,13 @@ def test_radon_iradon():
         reconstructed = rescale(reconstructed)
         delta = np.mean(np.abs(image - reconstructed))
 
-        ## print delta
-        ## import matplotlib.pyplot as plt
-        ## f, (ax1, ax2) = plt.subplots(1, 2)
-        ## ax1.imshow(image, cmap=plt.cm.gray)
-        ## ax2.imshow(reconstructed, cmap=plt.cm.gray)
-        ## plt.show()
+        if debug:
+            print delta
+            import matplotlib.pyplot as plt
+            f, (ax1, ax2) = plt.subplots(1, 2)
+            ax1.imshow(image, cmap=plt.cm.gray)
+            ax2.imshow(reconstructed, cmap=plt.cm.gray)
+            plt.show()
 
         assert delta < 0.05
 
@@ -60,7 +62,34 @@ def test_iradon_angles():
     # Loss of quality when the number of projections is reduced
     assert delta_80 > delta_200
 
+def test_radon_minimal():
+    """
+    Test for small images for various angles
+    """
+    thetas = [np.arange(180)]
+    for theta in thetas:
+        a = np.zeros((3, 3))
+        a[1, 1] = 1
+        p = radon(a, theta)
+        reconstructed = iradon(p, theta)
+        reconstructed /= np.max(reconstructed)
+        assert np.all(abs(a - reconstructed) < 0.3) 
         
+        b = np.zeros((4, 4))
+        b[1:3, 1:3] = 1
+        p = radon(b, theta)
+        reconstructed = iradon(p, theta)
+        reconstructed /= np.max(reconstructed)
+        assert np.all(abs(b - reconstructed) < 0.3) 
+        
+        c = np.zeros((5, 5))
+        c[1:3, 1:3] = 1
+        p = radon(c, theta)
+        reconstructed = iradon(p, theta)
+        reconstructed /= np.max(reconstructed)
+        assert np.all(abs(c - reconstructed) < 0.3) 
+     
+
 if __name__ == "__main__":
     run_module_suite()