Merge branch 'consistent_scaling'

4 files changed, 316 insertions, 329 deletions

diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index d04ffb8..e5d8f2b 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -7,9 +7,9 @@ import pylab
 # Display sinograms with mismatch on test failure
 DISPLAY=False
 
-NONUNITDET=False
-OBLIQUE=False
-FLEXVOL=False
+NONUNITDET=True
+OBLIQUE=True
+FLEXVOL=True
 NONSQUARE=False  # non-square pixels not supported yet by most projectors
 
 # Round interpolation weight to 8 bits to emulate CUDA texture unit precision
@@ -20,15 +20,8 @@ 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)
-
+# KNOWN FAILURES:
+# fan/strip relatively high numerical errors around 45 degrees
 
 
 # return length of intersection of the line through points src = (x,y)
@@ -454,23 +447,15 @@ class Test2DKernel(unittest.TestCase):
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (src, det) = center
 
-          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,
                         xmin, xmax, ymin, ymax,
                         1e-3)
-          a[i] = aw * detweight
-          b[i] = bw * detweight
-          c[i] = cw * detweight
+          a[i] = aw
+          b[i] = bw
+          c[i] = cw
         a = a.reshape(astra.functions.geom_size(pg))
         b = b.reshape(astra.functions.geom_size(pg))
         c = c.reshape(astra.functions.geom_size(pg))
@@ -494,17 +479,9 @@ class Test2DKernel(unittest.TestCase):
         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)
+            length = math.sqrt(1.0 + abs(yd/xd)**2) * pixsize[0]
             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]
@@ -512,9 +489,9 @@ class Test2DKernel(unittest.TestCase):
               # 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
+              l += w * length
           else:
-            length = math.sqrt(1.0 + abs(xd/yd)**2)
+            length = math.sqrt(1.0 + abs(xd/yd)**2) * pixsize[1]
             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]
@@ -522,7 +499,7 @@ class Test2DKernel(unittest.TestCase):
               # 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
+              l += w * length
           a[i] = l
         a = a.reshape(astra.functions.geom_size(pg))
         if not np.all(np.isfinite(a)):
@@ -532,21 +509,26 @@ class Test2DKernel(unittest.TestCase):
         if DISPLAY and x > TOL:
           display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
-      elif proj_type == 'distance_driven':
+      elif proj_type == 'distance_driven' and 'par' in type:
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
-          (xd, yd) = center[1] - center[0]
+          (src, det) = center
+          try:
+            detweight = pg['DetectorWidth']
+          except KeyError:
+            detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
+          (xd, yd) = det - src
           l = 0.0
           if np.abs(xd) > np.abs(yd): # horizontal ray
             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]
+              l += intersect_ray_vertical_segment(edge1, edge2, x, y_seg) * pixsize[0] / detweight
           else:
             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]
+              l += intersect_ray_horizontal_segment(edge1, edge2, y, x_seg) * pixsize[1] / detweight
           a[i] = l
         a = a.reshape(astra.functions.geom_size(pg))
         if not np.all(np.isfinite(a)):
@@ -560,6 +542,7 @@ class Test2DKernel(unittest.TestCase):
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (src, det) = center
+          detweight = effective_detweight(src, det, edge2[1] - edge1[1])
           det_dist = np.linalg.norm(src-det, ord=2)
           l = 0.0
           for j in range(rect_min[0], rect_max[0]):
@@ -570,7 +553,7 @@ class Test2DKernel(unittest.TestCase):
               ymin = origin[1] + (+0.5 * shape[1] - k - 1) * pixsize[1]
               ymax = origin[1] + (+0.5 * shape[1] - k) * pixsize[1]
               ycen = 0.5 * (ymin + ymax)
-              scale = det_dist / np.linalg.norm( src - np.array((xcen,ycen)), ord=2 )
+              scale = det_dist / (np.linalg.norm( src - np.array((xcen,ycen)), ord=2 ) * detweight)
               w = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax)
               l += w * scale
           a[i] = l
@@ -578,14 +561,20 @@ class Test2DKernel(unittest.TestCase):
         if not np.all(np.isfinite(a)):
           raise RuntimeError("Invalid value in reference sinogram")
         x = np.max(np.abs(sinogram-a))
-        TOL = 8e-3
+        # BUG: Known bug in fan/strip code around 45 degree projections causing larger errors than desirable
+        TOL = 4e-2
         if DISPLAY and x > TOL:
           display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
       elif proj_type == 'strip':
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
-          a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax)
+          (src, det) = center
+          try:
+            detweight = pg['DetectorWidth']
+          except KeyError:
+            detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
+          a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax) / detweight
         a = a.reshape(astra.functions.geom_size(pg))
         if not np.all(np.isfinite(a)):
           raise RuntimeError("Invalid value in reference sinogram")
@@ -594,46 +583,83 @@ class Test2DKernel(unittest.TestCase):
         if DISPLAY and x > TOL:
           display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
+      else:
+        raise RuntimeError("Unsupported projector")
 
-  def multi_test(self, type, proj_type):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test(type, proj_type)
-
-  def test_par(self):
-    self.multi_test('parallel', 'line')
-  def test_par_linear(self):
-    self.multi_test('parallel', 'linear')
-  def test_par_cuda(self):
-    self.multi_test('parallel', 'cuda')
-  def test_par_dd(self):
-    self.multi_test('parallel', 'distance_driven')
-  def test_par_strip(self):
-    self.multi_test('parallel', 'strip')
-  def test_fan(self):
-    self.multi_test('fanflat', 'line')
-  def test_fan_strip(self):
-    self.multi_test('fanflat', 'strip')
-  def test_fan_cuda(self):
-    self.multi_test('fanflat', 'cuda')
-  def test_parvec(self):
-    self.multi_test('parallel_vec', 'line')
-  def test_parvec_linear(self):
-    self.multi_test('parallel_vec', 'linear')
-  def test_parvec_dd(self):
-    self.multi_test('parallel_vec', 'distance_driven')
-  def test_parvec_strip(self):
-    self.multi_test('parallel_vec', 'strip')
-  def test_parvec_cuda(self):
-    self.multi_test('parallel_vec', 'cuda')
-  def test_fanvec(self):
-    self.multi_test('fanflat_vec', 'line')
-  def test_fanvec_cuda(self):
-    self.multi_test('fanflat_vec', 'cuda')
+  def single_test_adjoint(self, type, proj_type):
+      shape = np.random.randint(*range2d, size=2)
+      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.)
+          origin = (0.,0.)
+      vg = astra.create_vol_geom(shape[1], shape[0],
+                                 origin[0] - 0.5 * shape[0] * pixsize[0],
+                                 origin[0] + 0.5 * shape[0] * pixsize[0],
+                                 origin[1] - 0.5 * shape[1] * pixsize[1],
+                                 origin[1] + 0.5 * shape[1] * pixsize[1])
 
+      if type == 'parallel':
+        pg = gen_random_geometry_parallel()
+        projector_id = astra.create_projector(proj_type, pg, vg)
+      elif type == 'parallel_vec':
+        pg = gen_random_geometry_parallel_vec()
+        projector_id = astra.create_projector(proj_type, pg, vg)
+      elif type == 'fanflat':
+        pg = gen_random_geometry_fanflat()
+        projector_id = astra.create_projector(proj_type_to_fan(proj_type), pg, vg)
+      elif type == 'fanflat_vec':
+        pg = gen_random_geometry_fanflat_vec()
+        projector_id = astra.create_projector(proj_type_to_fan(proj_type), pg, vg)
 
+      for i in range(5):
+        X = np.random.random((shape[1], shape[0]))
+        Y = np.random.random(astra.geom_size(pg))
+
+        sinogram_id, fX = astra.create_sino(X, projector_id)
+        bp_id, fTY = astra.create_backprojection(Y, projector_id)
+
+        astra.data2d.delete(sinogram_id)
+        astra.data2d.delete(bp_id)
+
+        da = np.dot(fX.ravel(), Y.ravel())
+        db = np.dot(X.ravel(), fTY.ravel())
+        m = np.abs(da - db)
+        TOL = 1e-3 if 'cuda' not in proj_type else 1e-1
+        if m / da >= TOL:
+          print(vg)
+          print(pg)
+          print(m/da, da/db, da, db)
+        self.assertTrue(m / da < TOL)
+      astra.projector.delete(projector_id)
 
 
+  def multi_test(self, type, proj_type):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test(type, proj_type)
+  def multi_test_adjoint(self, type, proj_type):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test_adjoint(type, proj_type)
+
+__combinations = { 'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+                   'parallel_vec': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+                   'fanflat': [ 'line', 'strip', 'cuda' ],
+                   'fanflat_vec': [ 'line', 'cuda' ] }
+
+for k, l in __combinations.items():
+  for v in l:
+    def f(k,v):
+      return lambda self: self.multi_test(k, v)
+    def f_adj(k,v):
+      return lambda self: self.multi_test_adjoint(k, v)
+    setattr(Test2DKernel, 'test_' + k + '_' + v, f(k,v))
+    setattr(Test2DKernel, 'test_' + k + '_' + v + '_adjoint', f_adj(k,v))
 
 if __name__ == '__main__':
   unittest.main()
diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py
new file mode 100644
index 0000000..621fd8a
--- /dev/null
+++ b/tests/python/test_rec_scaling.py
@@ -0,0 +1,213 @@
+import numpy as np
+import unittest
+import astra
+import math
+import pylab
+
+DISPLAY=False
+
+def VolumeGeometries(is3D,noncube):
+  if not is3D:
+    for s in [0.8, 1.0, 1.25]:
+      yield astra.create_vol_geom(128, 128, -64*s, 64*s, -64*s, 64*s)
+  elif noncube:
+    for sx in [0.8, 1.0]:
+      for sy in [0.8, 1.0]:
+        for sz in [0.8, 1.0]:
+          yield astra.create_vol_geom(64, 64, 64, -32*sx, 32*sx, -32*sy, 32*sy, -32*sz, 32*sz)
+  else:
+    for s in [0.8, 1.0]:
+      yield astra.create_vol_geom(64, 64, 64, -32*s, 32*s, -32*s, 32*s, -32*s, 32*s)
+
+
+def ProjectionGeometries(type):
+  if type == 'parallel':
+    for dU in [0.8, 1.0, 1.25]:
+      yield astra.create_proj_geom('parallel', dU, 256, np.linspace(0,np.pi,180,False))
+  elif type == 'fanflat':
+    for dU in [0.8, 1.0, 1.25]:
+      for src in [500, 1000]:
+        for det in [0, 250, 500]:
+          yield astra.create_proj_geom('fanflat', dU, 256, np.linspace(0,2*np.pi,180,False), src, det)
+  elif type == 'parallel3d':
+    for dU in [0.8, 1.0]:
+      for dV in [0.8, 1.0]:
+        yield astra.create_proj_geom('parallel3d', dU, dV, 128, 128, np.linspace(0,np.pi,180,False))
+  elif type == 'parallel3d_vec':
+    for j in range(10):
+       Vectors = np.zeros([180,12])
+       wu = 0.6 + 0.8 * np.random.random()
+       wv = 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()
+         angle3 = 0.1*np.pi*np.random.random()
+         detc = 10 * np.random.random(size=3)
+         detu = [ math.cos(angle1) * wu, math.sin(angle1) * wu, 0 ]
+         detv = [ -math.sin(angle1) * math.sin(angle3) * wv, math.cos(angle1) * math.sin(angle3) * wv, math.cos(angle3) * wv ]
+         ray = [ math.sin(angle2) * l, -math.cos(angle2) * l, 0 ]
+         Vectors[i, :] = [ ray[0], ray[1], ray[2], detc[0], detc[1], detc[2], detu[0], detu[1], detu[2], detv[0], detv[1], detv[2] ]
+       pg = astra.create_proj_geom('parallel3d_vec', 128, 128, Vectors)
+       yield pg
+  elif type == 'cone':
+    for dU in [0.8, 1.0]:
+      for dV in [0.8, 1.0]:
+        for src in [500, 1000]:
+          for det in [0, 250]:
+            yield astra.create_proj_geom('cone', dU, dV, 128, 128, np.linspace(0,2*np.pi,180,False), src, det)
+  elif type == 'cone_vec':
+    for j in range(10):
+       Vectors = np.zeros([180,12])
+       wu = 0.6 + 0.8 * np.random.random()
+       wv = 0.6 + 0.8 * np.random.random()
+       for i in range(Vectors.shape[0]):
+         l = 256 * (0.5 * np.random.random())
+         angle1 = 2*np.pi*np.random.random()
+         angle2 = angle1 + 0.5 * np.random.random()
+         angle3 = 0.1*np.pi*np.random.random()
+         detc = 10 * np.random.random(size=3)
+         detu = [ math.cos(angle1) * wu, math.sin(angle1) * wu, 0 ]
+         detv = [ -math.sin(angle1) * math.sin(angle3) * wv, math.cos(angle1) * math.sin(angle3) * wv, math.cos(angle3) * wv ]
+         src = [ math.sin(angle2) * l, -math.cos(angle2) * l, 0 ]
+         Vectors[i, :] = [ src[0], src[1], src[2], detc[0], detc[1], detc[2], detu[0], detu[1], detu[2], detv[0], detv[1], detv[2] ]
+       pg = astra.create_proj_geom('parallel3d_vec', 128, 128, Vectors)
+       yield pg
+
+
+class TestRecScale(unittest.TestCase):
+  def single_test(self, geom_type, proj_type, alg, iters):
+    if alg == 'FBP' and 'fanflat' in geom_type:
+      self.skipTest('CPU FBP is parallel-beam only')
+    is3D = (geom_type in ['parallel3d', 'cone'])
+    for vg in VolumeGeometries(is3D, 'FDK' not in alg):
+      for pg in ProjectionGeometries(geom_type):
+        if not is3D:
+          vol = np.zeros((128,128),dtype=np.float32)
+          vol[50:70,50:70] = 1
+        else:
+          vol = np.zeros((64,64,64),dtype=np.float32)
+          vol[25:35,25:35,25:35] = 1
+        proj_id = astra.create_projector(proj_type, pg, vg)
+        if not is3D:
+          sino_id, sinogram = astra.create_sino(vol, proj_id)
+        else:
+          sino_id, sinogram = astra.create_sino3d_gpu(vol, pg, vg)
+        if not is3D:
+          DATA = astra.data2d
+        else:
+          DATA = astra.data3d
+
+        rec_id = DATA.create('-vol', vg, 0.0 if 'EM' not in alg else 1.0)
+
+        cfg = astra.astra_dict(alg)
+        cfg['ReconstructionDataId'] = rec_id
+        cfg['ProjectionDataId'] = sino_id
+        cfg['ProjectorId'] = proj_id
+        alg_id = astra.algorithm.create(cfg)
+
+        for i in range(iters):
+            astra.algorithm.run(alg_id, 1)
+        rec = DATA.get(rec_id)
+        astra.astra.delete([sino_id, alg_id, alg_id, proj_id])
+        if not is3D:
+          val = np.sum(rec[55:65,55:65]) / 100.
+        else:
+          val = np.sum(rec[27:32,27:32,27:32]) / 125.
+        TOL = 5e-2
+        if DISPLAY and abs(val-1.0) >= TOL:
+          print(geom_type, proj_type, alg, vg, pg)
+          print(val)
+          pylab.gray()
+          if not is3D:
+            pylab.imshow(rec)
+          else:
+            pylab.imshow(rec[:,32,:])
+          pylab.show()
+        self.assertTrue(abs(val-1.0) < TOL)
+
+  def single_test_adjoint3D(self, geom_type, proj_type):
+    for vg in VolumeGeometries(True, True):
+      for pg in ProjectionGeometries(geom_type):
+        for i in range(5):
+          X = np.random.random(astra.geom_size(vg))
+          Y = np.random.random(astra.geom_size(pg))
+          proj_id, fX = astra.create_sino3d_gpu(X, pg, vg)
+          bp_id, fTY = astra.create_backprojection3d_gpu(Y, pg, vg)
+
+          astra.data3d.delete([proj_id, bp_id])
+
+          da = np.dot(fX.ravel(), Y.ravel())
+          db = np.dot(X.ravel(), fTY.ravel())
+          m = np.abs(da - db)
+          TOL = 1e-1
+          if m / da >= TOL:
+            print(vg)
+            print(pg)
+            print(m/da, da/db, da, db)
+          self.assertTrue(m / da < TOL)
+
+
+
+
+
+__combinations = {
+  'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+  'fanflat': [ 'line_fanflat', 'strip_fanflat', 'cuda' ],
+  'parallel3d': [ 'cuda3d' ],
+  'cone': [ 'cuda3d' ],
+}
+
+__combinations_adjoint = {
+  'parallel3d': [ 'cuda3d' ],
+  'cone': [ 'cuda3d' ],
+  'parallel3d_vec': [ 'cuda3d' ],
+  'cone_vec': [ 'cuda3d' ],
+}
+
+__algs = {
+   'SIRT': 50, 'SART': 10*180, 'CGLS': 30,
+   'FBP': 1
+}
+
+__algs_CUDA = {
+  'SIRT_CUDA': 50, 'SART_CUDA': 10*180, 'CGLS_CUDA': 30, 'EM_CUDA': 50,
+  'FBP_CUDA': 1
+}
+
+__algs_parallel3d = {
+  'SIRT3D_CUDA': 200, 'CGLS3D_CUDA': 20,
+}
+
+__algs_cone = {
+  'SIRT3D_CUDA': 200, 'CGLS3D_CUDA': 20,
+  'FDK_CUDA': 1
+}
+
+
+
+for k, l in __combinations.items():
+  for v in l:
+    is3D = (k in ['parallel3d', 'cone'])
+    if k == 'parallel3d':
+      A = __algs_parallel3d
+    elif k == 'cone':
+      A = __algs_cone
+    elif v == 'cuda':
+      A = __algs_CUDA
+    else:
+      A = __algs
+    for a, i in A.items():
+      def f(k, v, a, i):
+        return lambda self: self.single_test(k, v, a, i)
+      setattr(TestRecScale, 'test_' + a + '_' + k + '_' + v, f(k,v,a,i))
+
+for k, l in __combinations_adjoint.items():
+  for v in l:
+    def g(k, v):
+      return lambda self: self.single_test_adjoint3D(k, v)
+    setattr(TestRecScale, 'test_adjoint_' + k + '_' + v, g(k,v))
+
+if __name__ == '__main__':
+  unittest.main()
+
diff --git a/tests/test_ParallelBeamLineKernelProjector2D.cpp b/tests/test_ParallelBeamLineKernelProjector2D.cpp
deleted file mode 100644
index 71130c1..0000000
--- a/tests/test_ParallelBeamLineKernelProjector2D.cpp
+++ /dev/null
@@ -1,82 +0,0 @@
-/*
------------------------------------------------------------------------
-Copyright: 2010-2018, imec Vision Lab, University of Antwerp
-           2014-2018, CWI, Amsterdam
-
-Contact: astra@astra-toolbox.com
-Website: http://www.astra-toolbox.com/
-
-This file is part of the ASTRA Toolbox.
-
-
-The ASTRA Toolbox is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-The ASTRA Toolbox is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-
------------------------------------------------------------------------
-*/
-
-
-#define BOOST_TEST_DYN_LINK
-#include <boost/test/unit_test.hpp>
-#include <boost/test/auto_unit_test.hpp>
-#include <boost/test/floating_point_comparison.hpp>
-
-#include "astra/ParallelBeamLineKernelProjector2D.h"
-#include "astra/ParallelProjectionGeometry2D.h"
-#include "astra/VolumeGeometry2D.h"
-
-struct TestParallelBeamLineKernelProjector2D {
-        TestParallelBeamLineKernelProjector2D()
-	{
-		astra::float32 angles[] = { 1.0f };
-		BOOST_REQUIRE( projGeom.initialize(1, 9, 1.0f, angles) );
-		BOOST_REQUIRE( volGeom.initialize(6, 4) );
-		BOOST_REQUIRE( proj.initialize(&projGeom, &volGeom) );
-	}
-        ~TestParallelBeamLineKernelProjector2D()
-	{
-
-	}
-
-	astra::CParallelBeamLineKernelProjector2D proj;
-	astra::CParallelProjectionGeometry2D projGeom;
-	astra::CVolumeGeometry2D volGeom;
-};
-
-BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_General, TestParallelBeamLineKernelProjector2D )
-{
-
-}
-
-BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_Rectangle, TestParallelBeamLineKernelProjector2D )
-{
-	int iMax = proj.getProjectionWeightsCount(0);
-	BOOST_REQUIRE(iMax > 0);
-
-	astra::SPixelWeight* pPix = new astra::SPixelWeight[iMax];
-	BOOST_REQUIRE(pPix);
-
-	int iCount;
-	proj.computeSingleRayWeights(0, 4, pPix, iMax, iCount); 
-	BOOST_REQUIRE(iCount <= iMax);
-
-	astra::float32 fWeight = 0;
-	for (int i = 0; i < iCount; ++i)
-		fWeight += pPix[i].m_fWeight;
-
-	BOOST_CHECK_SMALL(fWeight - 7.13037f, 0.00001f); // 6 / sin(1)
-
-	delete[] pPix;
-}
-
-
diff --git a/tests/test_ParallelBeamLinearKernelProjector2D.cpp b/tests/test_ParallelBeamLinearKernelProjector2D.cpp
deleted file mode 100644
index 4610aa5..0000000
--- a/tests/test_ParallelBeamLinearKernelProjector2D.cpp
+++ /dev/null
@@ -1,170 +0,0 @@
-/*
------------------------------------------------------------------------
-Copyright: 2010-2018, imec Vision Lab, University of Antwerp
-           2014-2018, CWI, Amsterdam
-
-Contact: astra@astra-toolbox.com
-Website: http://www.astra-toolbox.com/
-
-This file is part of the ASTRA Toolbox.
-
-
-The ASTRA Toolbox is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-The ASTRA Toolbox is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-
------------------------------------------------------------------------
-*/
-
-
-#define BOOST_TEST_DYN_LINK
-#include <boost/test/unit_test.hpp>
-#include <boost/test/auto_unit_test.hpp>
-#include <boost/test/floating_point_comparison.hpp>
-
-#include "astra/ParallelBeamLineKernelProjector2D.h"
-#include "astra/ParallelBeamLinearKernelProjector2D.h"
-#include "astra/ParallelBeamStripKernelProjector2D.h"
-#include "astra/ParallelProjectionGeometry2D.h"
-#include "astra/VolumeGeometry2D.h"
-#include "astra/ProjectionGeometry2D.h"
-
-#include <ctime>
-
-using astra::float32;
-
-struct TestParallelBeamLinearKernelProjector2D {
-        TestParallelBeamLinearKernelProjector2D()
-	{
-		astra::float32 angles[] = { 2.6f };
-		BOOST_REQUIRE( projGeom.initialize(1, 3, 1.0f, angles) );
-		BOOST_REQUIRE( volGeom.initialize(3, 2) );
-		BOOST_REQUIRE( proj.initialize(&projGeom, &volGeom) );
-	}
-        ~TestParallelBeamLinearKernelProjector2D()
-	{
-
-	}
-
-	astra::CParallelBeamLinearKernelProjector2D proj;
-	astra::CParallelProjectionGeometry2D projGeom;
-	astra::CVolumeGeometry2D volGeom;
-};
-
-BOOST_FIXTURE_TEST_CASE( testParallelBeamLinearKernelProjector2D_General, TestParallelBeamLinearKernelProjector2D )
-{
-
-}
-
-
-// Compute linear kernel for a single volume pixel/detector pixel combination
-float32 compute_linear_kernel(const astra::CProjectionGeometry2D& projgeom, const astra::CVolumeGeometry2D& volgeom,
-                  int iX, int iY, int iDet, float32 fAngle)
-{
-	// projection of center of volume pixel on detector array
-	float32 fDetProj = (iX - (volgeom.getGridColCount()-1.0f)/2.0f ) * volgeom.getPixelLengthX() * cos(fAngle) - (iY - (volgeom.getGridRowCount()-1.0f)/2.0f ) * volgeom.getPixelLengthY() * sin(fAngle);
-
-	// start of detector pixel on detector array
-	float32 fDetStart = projgeom.indexToDetectorOffset(iDet) - 0.5f;
-
-//	printf("(%d,%d,%d): %f in (%f,%f)\n", iX,iY,iDet,fDetProj, fDetStart, fDetStart+1.0f);
-
-	// projection of center of next volume pixel on detector array
-	float32 fDetStep;
-	// length of projection ray through volume pixel
-	float32 fWeight;
-
-	if (fabs(cos(fAngle)) > fabs(sin(fAngle))) {
-		fDetStep = volgeom.getPixelLengthY() * fabs(cos(fAngle));
-		fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthX() * 1.0f / fabs(cos(fAngle));
-	} else {
-		fDetStep = volgeom.getPixelLengthX() * fabs(sin(fAngle));
-		fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthY() * 1.0f / fabs(sin(fAngle));
-	}
-
-//	printf("step: %f\n   weight: %f\n", fDetStep, fWeight);
-
-	// center of detector pixel on detector array
-	float32 fDetCenter = fDetStart + 0.5f;
-
-	// unweighted contribution of this volume pixel:
-	// linear interpolation between
-	//  fDetCenter - fDetStep    |---> 0
-	//  fDetCenter               |---> 1
-	//  fDetCenter + fDetStep    |---> 0
-	float32 fBase;
-	if (fDetCenter <= fDetProj) {
-		fBase = (fDetCenter - (fDetProj - fDetStep))/fDetStep;
-	} else {
-		fBase = ((fDetProj + fDetStep) - fDetCenter)/fDetStep;
-	}
-//	printf("base: %f\n", fBase);
-	if (fBase < 0) fBase = 0;
-	return fBase * fWeight;
-}
-
-BOOST_AUTO_TEST_CASE( testParallelBeamLinearKernelProjector2D_Rectangles )
-{
-	astra::CParallelBeamLinearKernelProjector2D proj;
-	astra::CParallelProjectionGeometry2D projGeom;
-	astra::CVolumeGeometry2D volGeom;
-
-	const unsigned int iRandomTestCount = 100;
-
-	unsigned int iSeed = time(0);
-	srand(iSeed);
-
-	for (unsigned int iTest = 0; iTest < iRandomTestCount; ++iTest) {
-		int iDetectorCount = 1 + (rand() % 100);
-		int iRows = 1 + (rand() % 100);
-		int iCols = 1 + (rand() % 100);
-		
-		
-		astra::float32 angles[] = { rand() * 2.0f*astra::PI / RAND_MAX };
-		projGeom.initialize(1, iDetectorCount, 0.8f, angles);
-		volGeom.initialize(iCols, iRows);
-		proj.initialize(&projGeom, &volGeom);
-
-		int iMax = proj.getProjectionWeightsCount(0);
-		BOOST_REQUIRE(iMax > 0);
-
-		astra::SPixelWeight* pPix = new astra::SPixelWeight[iMax];
-		BOOST_REQUIRE(pPix);
-
-		astra::float32 fWeight = 0;
-		for (int iDet = 0; iDet < projGeom.getDetectorCount(); ++iDet) {
-			int iCount;
-			proj.computeSingleRayWeights(0, iDet, pPix, iMax, iCount); 
-			BOOST_REQUIRE(iCount <= iMax);
-
-			astra::float32 fW = 0;
-			for (int i = 0; i < iCount; ++i) {
-				float32 fTest = compute_linear_kernel(
-				            projGeom,
-				            volGeom,
-				            pPix[i].m_iIndex % volGeom.getGridColCount(),
-				            pPix[i].m_iIndex / volGeom.getGridColCount(),
-				            iDet,
-				            projGeom.getProjectionAngle(0));
-				BOOST_CHECK_SMALL( pPix[i].m_fWeight - fTest, 0.0004f);
-				fW += pPix[i].m_fWeight;
-			}
-
-			fWeight += fW;
-
-		}
-
-		delete[] pPix;
-	}
-}
-
-