Clean up projector unit tests

1 files changed, 191 insertions, 151 deletions

diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index b1bb7d5..e7bca98 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -4,10 +4,35 @@ import astra
 import math
 import pylab
 
+# Display sinograms with mismatch on test failure
+DISPLAY=False
+
+NONUNITDET=False
+OBLIQUE=False
+FLEXVOL=False
+NONSQUARE=False  # non-square pixels not supported yet by most projectors
+
+# Round interpolation weight to 8 bits to emulate CUDA texture unit precision
+CUDA_8BIT_LINEAR=True
+CUDA_TOL=2e-2
+
+nloops = 50
+seed = 123
+
+
+# FAILURES:
+# fan/cuda with flexible volume
+# detweight for fan/cuda
+# fan/strip relatively high numerical errors?
+# parvec/line+linear for oblique
+
+# INCONSISTENCY:
+# effective_detweight vs norm(detu) in line/linear (oblique)
+
+
+
 # return length of intersection of the line through points src = (x,y)
 # and det (x,y), and the rectangle defined by xmin, ymin, xmax, ymax
-#
-# TODO: Generalize from 2D to n-dimensional
 def intersect_line_rectangle(src, det, xmin, xmax, ymin, ymax):
   EPS = 1e-5
 
@@ -89,7 +114,6 @@ def intersect_ray_horizontal_segment(edge1, edge2, y, x_seg):
   (x1, x2) = np.sort(x_seg)
   l = np.max([e1, x1])
   r = np.min([e2, x2])
-  #print(edge1, edge2, y, x_seg, r-l)
   return np.max([r-l, 0.0])
 
 def intersect_ray_vertical_segment(edge1, edge2, x, y_seg):
@@ -128,7 +152,8 @@ def area_signed(a, b):
 
 # is c to the left of ab
 def is_left_of(a, b, c):
-  return area_signed( (b[0] - a[0], b[1] - a[1]), (c[0] - a[0], c[1] - a[1]) ) > 0
+  EPS = 1e-5
+  return area_signed( (b[0] - a[0], b[1] - a[1]), (c[0] - a[0], c[1] - a[1]) ) > EPS
 
 # compute area of rect on left side of line
 def halfarea_rect_line(src, det, xmin, xmax, ymin, ymax):
@@ -172,6 +197,13 @@ def intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax):
   return abs(s1 - s2)
 
 
+# width of projection of detector orthogonal to ray direction
+# i.e., effective detector width
+def effective_detweight(src, det, u):
+  ray = np.array(det) - np.array(src)
+  ray = ray / np.linalg.norm(ray, ord=2)
+  return abs(area_signed(ray, u))
+
 
 # LINE GENERATORS
 # ---------------
@@ -264,39 +296,56 @@ range2d = ( 8, 64 )
 
 
 def gen_random_geometry_fanflat():
-  pg = astra.create_proj_geom('fanflat', 0.6 + 0.8 * np.random.random(), np.random.randint(*range2d), np.linspace(0, 2*np.pi, np.random.randint(*range2d), endpoint=False), 256 * (0.5 + np.random.random()), 256 * np.random.random())
+  if not NONUNITDET:
+    w = 1.0
+  else:
+    w = 0.6 + 0.8 * np.random.random()
+  pg = astra.create_proj_geom('fanflat', w, np.random.randint(*range2d), np.linspace(0, 2*np.pi, np.random.randint(*range2d), endpoint=False), 256 * (0.5 + np.random.random()), 256 * np.random.random())
   return pg
 
 def gen_random_geometry_parallel():
-  pg = astra.create_proj_geom('parallel', 0.8 + 0.4 * np.random.random(), np.random.randint(*range2d), np.linspace(0, 2*np.pi, np.random.randint(*range2d), endpoint=False))
+  if not NONUNITDET:
+    w = 1.0
+  else:
+    w = 0.8 + 0.4 * np.random.random()
+  pg = astra.create_proj_geom('parallel', w, np.random.randint(*range2d), np.linspace(0, 2*np.pi, np.random.randint(*range2d), endpoint=False))
   return pg
 
 def gen_random_geometry_fanflat_vec():
   Vectors = np.zeros([16,6])
   # We assume constant detector width in these tests
-  w = 0.6 + 0.8 * np.random.random()
+  if not NONUNITDET:
+    w = 1.0
+  else:
+    w = 0.6 + 0.8 * np.random.random()
   for i in range(Vectors.shape[0]):
     angle1 = 2*np.pi*np.random.random()
-    angle2 = angle1 + 0.5 * np.random.random()
+    if OBLIQUE:
+      angle2 = angle1 + 0.5 * np.random.random()
+    else:
+      angle2 = angle1
     dist1 = 256 * (0.5 + np.random.random())
     detc = 10 * np.random.random(size=2)
     detu = [ math.cos(angle1) * w, math.sin(angle1) * w ]
     src = [ math.sin(angle2) * dist1, -math.cos(angle2) * dist1 ]
     Vectors[i, :] = [ src[0], src[1], detc[0], detc[1], detu[0], detu[1] ]
   pg = astra.create_proj_geom('fanflat_vec', np.random.randint(*range2d), Vectors)
-
-  # TODO: Randomize more
-  pg = astra.create_proj_geom('fanflat_vec', np.random.randint(*range2d), Vectors)
   return pg
 
 def gen_random_geometry_parallel_vec():
   Vectors = np.zeros([16,6])
   # We assume constant detector width in these tests
-  w = 0.6 + 0.8 * np.random.random()
+  if not NONUNITDET:
+    w = 1.0
+  else:
+    w = 0.6 + 0.8 * np.random.random()
   for i in range(Vectors.shape[0]):
     l = 0.6 + 0.8 * np.random.random()
     angle1 = 2*np.pi*np.random.random()
-    angle2 = angle1 + 0.5 * np.random.random()
+    if OBLIQUE:
+      angle2 = angle1 + 0.5 * np.random.random()
+    else:
+      angle2 = angle1
     detc = 10 * np.random.random(size=2)
     detu = [ math.cos(angle1) * w, math.sin(angle1) * w ]
     ray = [ math.sin(angle2) * l, -math.cos(angle2) * l ]
@@ -307,11 +356,39 @@ def gen_random_geometry_parallel_vec():
 
 
 
-nloops = 50
-seed = 123
-
 def proj_type_to_fan(t):
-  return t + '_fanflat'
+  if t == 'cuda':
+    return t
+  else:
+    return t + '_fanflat'
+
+def display_mismatch(data, sinogram, a):
+  pylab.gray()
+  pylab.imshow(data)
+  pylab.figure()
+  pylab.imshow(sinogram)
+  pylab.figure()
+  pylab.imshow(a)
+  pylab.figure()
+  pylab.imshow(sinogram-a)
+  pylab.show()
+
+def display_mismatch_triple(data, sinogram, a, b, c):
+  pylab.gray()
+  pylab.imshow(data)
+  pylab.figure()
+  pylab.imshow(sinogram)
+  pylab.figure()
+  pylab.imshow(b)
+  pylab.figure()
+  pylab.imshow(a)
+  pylab.figure()
+  pylab.imshow(c)
+  pylab.figure()
+  pylab.imshow(sinogram-a)
+  pylab.figure()
+  pylab.imshow(c-sinogram)
+  pylab.show()
 
 class Test2DKernel(unittest.TestCase):
   def single_test(self, type, proj_type):
@@ -319,9 +396,11 @@ class Test2DKernel(unittest.TestCase):
       # these rectangles are biased, but that shouldn't matter
       rect_min = [ np.random.randint(0, a) for a in shape ]
       rect_max = [ np.random.randint(rect_min[i]+1, shape[i]+1) for i in range(len(shape))]
-      if True:
-          #pixsize = 0.5 + np.random.random(size=2)
-          pixsize = np.array([0.5, 0.5]) + np.random.random()
+      if FLEXVOL:
+          if not NONSQUARE:
+            pixsize = np.array([0.5, 0.5]) + np.random.random()
+          else:
+            pixsize = 0.5 + np.random.random(size=2)
           origin = 10 * np.random.random(size=2)
       else:
           pixsize = (1.,1.)
@@ -331,7 +410,6 @@ class Test2DKernel(unittest.TestCase):
                                  origin[0] + 0.5 * shape[0] * pixsize[0],
                                  origin[1] - 0.5 * shape[1] * pixsize[1],
                                  origin[1] + 0.5 * shape[1] * pixsize[1])
-      #print(vg)
 
       if type == 'parallel':
         pg = gen_random_geometry_parallel()
@@ -352,6 +430,8 @@ class Test2DKernel(unittest.TestCase):
 
       sinogram_id, sinogram = astra.create_sino(data, projector_id)
 
+      self.assertTrue(np.all(np.isfinite(sinogram)))
+
       #print(pg)
       #print(vg)
 
@@ -359,11 +439,11 @@ class Test2DKernel(unittest.TestCase):
 
       astra.projector.delete(projector_id)
 
-      # Weight for pixel / voxel size
-      try:
-        detweight = pg['DetectorWidth']
-      except KeyError:
-        detweight = np.sqrt(pg['Vectors'][0,4]**2 + pg['Vectors'][0,5]**2)
+      # NB: Flipped y-axis here, since that is how astra interprets 2D volumes
+      xmin = origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0]
+      xmax = origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0]
+      ymin = origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1]
+      ymax = origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1]
 
       if proj_type == 'line':
 
@@ -371,200 +451,161 @@ class Test2DKernel(unittest.TestCase):
         b = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         c = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
 
-        i = 0
-        #print( origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0], origin[1] + (-0.5 * shape[1] + rect_min[1]) * pixsize[1], origin[1] + (-0.5 * shape[1] + rect_max[1]) * pixsize[1])
-        for center, edge1, edge2 in gen_lines(pg):
+        for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (src, det) = center
-          #print(src,det)
 
-          # NB: Flipped y-axis here, since that is how astra interprets 2D volumes
+          try:
+            detweight = pg['DetectorWidth']
+          except KeyError:
+            if 'fan' not in type:
+              detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
+            else:
+              detweight = np.linalg.norm(pg['Vectors'][i//pg['DetectorCount'],4:6], ord=2)
+
           # We compute line intersections with slightly bigger (cw) and
           # smaller (aw) rectangles, and see if the kernel falls
           # between these two values.
           (aw,bw,cw) = intersect_line_rectangle_interval(src, det,
-                        origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0],
-                        origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0],
-                        origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1],
-                        origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1],
+                        xmin, xmax, ymin, ymax,
                         1e-3)
-          a[i] = aw
-          b[i] = bw
-          c[i] = cw
-          i += 1
-        a *= detweight
-        b *= detweight
-        c *= detweight
+          a[i] = aw * detweight
+          b[i] = bw * detweight
+          c[i] = cw * detweight
         a = a.reshape(astra.functions.geom_size(pg))
         b = b.reshape(astra.functions.geom_size(pg))
         c = c.reshape(astra.functions.geom_size(pg))
 
+        if not np.all(np.isfinite(a)):
+          raise RuntimeError("Invalid value in reference sinogram")
+        if not np.all(np.isfinite(b)):
+          raise RuntimeError("Invalid value in reference sinogram")
+        if not np.all(np.isfinite(c)):
+          raise RuntimeError("Invalid value in reference sinogram")
+        self.assertTrue(np.all(np.isfinite(sinogram)))
+
         # Check if sinogram lies between a and c
         y = np.min(sinogram-a)
         z = np.min(c-sinogram)
-        x = np.max(np.abs(sinogram-b)) # ideally this is small, but can be large
-                                       # due to discontinuities in line kernel
+        if DISPLAY and (z < 0 or y < 0):
+          display_mismatch_triple(data, sinogram, a, b, c)
         self.assertFalse(z < 0 or y < 0)
-        if z < 0 or y < 0:
-          print(y,z,x)
-          pylab.gray()
-          pylab.imshow(data)
-          pylab.figure()
-          pylab.imshow(sinogram)
-          pylab.figure()
-          pylab.imshow(b)
-          pylab.figure()
-          pylab.imshow(a)
-          pylab.figure()
-          pylab.imshow(c)
-          pylab.figure()
-          pylab.imshow(sinogram-a)
-          pylab.figure()
-          pylab.imshow(c-sinogram)
-          pylab.show()
-      elif proj_type == 'linear':
+      elif proj_type == 'linear' or proj_type == 'cuda':
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
-        i = 0
-        for center, edge1, edge2 in gen_lines(pg):
-          (xd, yd) = center[1] - center[0]
+        for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
+          (src, det) = center
+          (xd, yd) = det - src
+          try:
+            detweight = pg['DetectorWidth']
+          except KeyError:
+            if 'fan' not in type:
+              detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
+            else:
+              detweight = np.linalg.norm(pg['Vectors'][i//pg['DetectorCount'],4:6], ord=2)
+
           l = 0.0
           if np.abs(xd) > np.abs(yd): # horizontal ray
             length = math.sqrt(1.0 + abs(yd/xd)**2)
-            y_seg = (origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1],
-                     origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1])
+            y_seg = (ymin, ymax)
             for j in range(rect_min[0], rect_max[0]):
               x = origin[0] + (-0.5 * shape[0] + j + 0.5) * pixsize[0]
               w = intersect_line_vertical_segment_linear(center[0], center[1], x, y_seg, pixsize[1])
+              # limited interpolation precision with cuda
+              if CUDA_8BIT_LINEAR and proj_type == 'cuda':
+                w = np.round(w * 256.0) / 256.0
               l += w * length * pixsize[0] * detweight
           else:
             length = math.sqrt(1.0 + abs(xd/yd)**2)
-            x_seg = (origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0],
-                     origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0])
+            x_seg = (xmin, xmax)
             for j in range(rect_min[1], rect_max[1]):
               y = origin[1] + (+0.5 * shape[1] - j - 0.5) * pixsize[1]
               w = intersect_line_horizontal_segment_linear(center[0], center[1], y, x_seg, pixsize[0])
+              # limited interpolation precision with cuda
+              if CUDA_8BIT_LINEAR and proj_type == 'cuda':
+                w = np.round(w * 256.0) / 256.0
               l += w * length * pixsize[1] * detweight
           a[i] = l
-          i += 1
         a = a.reshape(astra.functions.geom_size(pg))
+        if not np.all(np.isfinite(a)):
+          raise RuntimeError("Invalid value in reference sinogram")
         x = np.max(np.abs(sinogram-a))
-        TOL = 2e-3
-        if True and x > TOL:
-          pylab.gray()
-          pylab.imshow(data)
-          pylab.figure()
-          pylab.imshow(sinogram)
-          pylab.figure()
-          pylab.imshow(a)
-          pylab.figure()
-          pylab.imshow(sinogram-a)
-          pylab.show()
+        TOL = 2e-3 if proj_type != 'cuda' else CUDA_TOL
+        if DISPLAY and x > TOL:
+          display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
       elif proj_type == 'distance_driven':
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
-        i = 0
-        for center, edge1, edge2 in gen_lines(pg):
+        for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (xd, yd) = center[1] - center[0]
           l = 0.0
           if np.abs(xd) > np.abs(yd): # horizontal ray
-            y_seg = (origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1],
-                     origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1])
+            y_seg = (ymin, ymax)
             for j in range(rect_min[0], rect_max[0]):
               x = origin[0] + (-0.5 * shape[0] + j + 0.5) * pixsize[0]
               l += intersect_ray_vertical_segment(edge1, edge2, x, y_seg) * pixsize[0]
           else:
-            x_seg = (origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0],
-                     origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0])
+            x_seg = (xmin, xmax)
             for j in range(rect_min[1], rect_max[1]):
               y = origin[1] + (+0.5 * shape[1] - j - 0.5) * pixsize[1]
               l += intersect_ray_horizontal_segment(edge1, edge2, y, x_seg) * pixsize[1]
           a[i] = l
-          i += 1
         a = a.reshape(astra.functions.geom_size(pg))
+        if not np.all(np.isfinite(a)):
+          raise RuntimeError("Invalid value in reference sinogram")
         x = np.max(np.abs(sinogram-a))
         TOL = 2e-3
-        if x > TOL:
-          pylab.gray()
-          pylab.imshow(data)
-          pylab.figure()
-          pylab.imshow(sinogram)
-          pylab.figure()
-          pylab.imshow(a)
-          pylab.figure()
-          pylab.imshow(sinogram-a)
-          pylab.show()
+        if DISPLAY and x > TOL:
+          display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
       elif proj_type == 'strip':
-        xmin = origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0]
-        xmax = origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0]
-        ymin = origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1]
-        ymax = origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1]
-
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
-        i = 0
-        for center, edge1, edge2 in gen_lines(pg):
+        for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax)
-          i += 1
 
         a = a.reshape(astra.functions.geom_size(pg))
+        if not np.all(np.isfinite(a)):
+          raise RuntimeError("Invalid value in reference sinogram")
         x = np.max(np.abs(sinogram-a))
-        # TODO: Investigate tolerance for fanflat/strip
-        TOL = 2e-3
-        if False and x > TOL:
-          pylab.gray()
-          pylab.imshow(data)
-          pylab.figure()
-          pylab.imshow(sinogram)
-          pylab.figure()
-          pylab.imshow(a)
-          pylab.figure()
-          pylab.imshow(sinogram-a)
-          pylab.show()
+        TOL = 8e-3
+        if DISPLAY and x > TOL:
+          display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
 
-  def test_par(self):
+  def multi_test(self, type, proj_type):
     np.random.seed(seed)
     for _ in range(nloops):
-      self.single_test('parallel', 'line')
+      self.single_test(type, proj_type)
+
+  def test_par(self):
+    self.multi_test('parallel', 'line')
   def test_par_linear(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel', 'linear')
+    self.multi_test('parallel', 'linear')
+  def test_par_cuda(self):
+    self.multi_test('parallel', 'cuda')
   def test_par_dd(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel', 'distance_driven')
+    self.multi_test('parallel', 'distance_driven')
   def test_par_strip(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel', 'strip')
+    self.multi_test('parallel', 'strip')
   def test_fan(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('fanflat', 'line')
+    self.multi_test('fanflat', 'line')
   def test_fan_strip(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('fanflat', 'strip')
+    self.multi_test('fanflat', 'strip')
+  def test_fan_cuda(self):
+    self.multi_test('fanflat', 'cuda')
   def test_parvec(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel_vec', 'line')
+    self.multi_test('parallel_vec', 'line')
   def test_parvec_linear(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel_vec', 'linear')
+    self.multi_test('parallel_vec', 'linear')
   def test_parvec_dd(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel_vec', 'distance_driven')
+    self.multi_test('parallel_vec', 'distance_driven')
   def test_parvec_strip(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('parallel_vec', 'strip')
+    self.multi_test('parallel_vec', 'strip')
+  def test_parvec_cuda(self):
+    self.multi_test('parallel_vec', 'cuda')
   def test_fanvec(self):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test('fanflat_vec', 'line')
+    self.multi_test('fanflat_vec', 'line')
+  def test_fanvec_cuda(self):
+    self.multi_test('fanflat_vec', 'cuda')
+
 
 
 
@@ -572,4 +613,3 @@ class Test2DKernel(unittest.TestCase):
 if __name__ == '__main__':
   unittest.main()
 
-#print(intersect_line_rectangle((0.,-256.),(-27.,0.),11.6368454385 20.173128227 3.18989047649 5.62882841606)