Merge branch 'parallel_vec'

3 files changed, 313 insertions, 3 deletions

diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
new file mode 100644
index 0000000..de68033
--- /dev/null
+++ b/tests/python/test_line2d.py
@@ -0,0 +1,310 @@
+import numpy as np
+import unittest
+import astra
+import math
+import pylab
+
+# 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
+
+  if np.abs(src[0] - det[0]) < EPS:
+    if src[0] >= xmin and src[0] < xmax:
+      return ymax - ymin
+    else:
+      return 0.0
+  if np.abs(src[1] - det[1]) < EPS:
+    if src[1] >= ymin and src[1] < ymax:
+      return xmax - xmin
+    else:
+      return 0.0
+
+  n = np.sqrt((det[0] - src[0]) ** 2 + (det[1] - src[1]) ** 2)
+
+  check = [ (-(xmin - src[0]), -(det[0] - src[0]) / n ),
+            (xmax - src[0], (det[0] - src[0]) / n ),
+            (-(ymin - src[1]), -(det[1] - src[1]) / n ),
+            (ymax - src[1], (det[1] - src[1]) / n ) ]
+
+  pre = [ -np.Inf ]
+  post = [ np.Inf ]
+
+  for p, q in check:
+    r = p / (1.0 * q)
+    if q > 0:
+      post.append(r)    # exiting half-plane
+    else:
+      pre.append(r)     # entering half-plane
+
+  end_r = np.min(post)
+  start_r = np.max(pre)
+
+  if end_r > start_r:
+    return end_r - start_r
+  else:
+    return 0.0
+
+def intersect_line_rectangle_feather(src, det, xmin, xmax, ymin, ymax, feather):
+  return intersect_line_rectangle(src, det,
+                                  xmin-feather, xmax+feather,
+                                  ymin-feather, ymax+feather)
+
+def intersect_line_rectangle_interval(src, det, xmin, xmax, ymin, ymax, f):
+  a = intersect_line_rectangle_feather(src, det, xmin, xmax, ymin, ymax, -f)
+  b = intersect_line_rectangle(src, det, xmin, xmax, ymin, ymax)
+  c = intersect_line_rectangle_feather(src, det, xmin, xmax, ymin, ymax, f)
+  return (a,b,c)
+
+def gen_lines_fanflat(proj_geom):
+  angles = proj_geom['ProjectionAngles']
+  for theta in angles:
+    #theta = -theta
+    src = ( math.sin(theta) * proj_geom['DistanceOriginSource'],
+           -math.cos(theta) * proj_geom['DistanceOriginSource'] )
+    detc= (-math.sin(theta) * proj_geom['DistanceOriginDetector'],
+            math.cos(theta) * proj_geom['DistanceOriginDetector'] )
+    detu= ( math.cos(theta) * proj_geom['DetectorWidth'],
+            math.sin(theta) * proj_geom['DetectorWidth'] )
+
+    src = np.array(src, dtype=np.float64)
+    detc= np.array(detc, dtype=np.float64)
+    detu= np.array(detu, dtype=np.float64)
+
+    detb= detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu
+
+    for i in range(proj_geom['DetectorCount']):
+      yield (src, detb + i * detu)
+
+def gen_lines_fanflat_vec(proj_geom):
+  v = proj_geom['Vectors']
+  for i in range(v.shape[0]):
+    src = v[i,0:2]
+    detc = v[i,2:4]
+    detu = v[i,4:6]
+
+    detb = detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu
+    for i in range(proj_geom['DetectorCount']):
+      yield (src, detb + i * detu)
+
+def gen_lines_parallel(proj_geom):
+  angles = proj_geom['ProjectionAngles']
+  for theta in angles:
+    ray = ( math.sin(theta),
+           -math.cos(theta) )
+    detc= (0, 0 )
+    detu= ( math.cos(theta) * proj_geom['DetectorWidth'],
+            math.sin(theta) * proj_geom['DetectorWidth'] )
+
+    ray = np.array(ray, dtype=np.float64)
+    detc= np.array(detc, dtype=np.float64)
+    detu= np.array(detu, dtype=np.float64)
+
+
+    detb= detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu
+
+    for i in range(proj_geom['DetectorCount']):
+      yield (detb + i * detu - ray, detb + i * detu)
+
+def gen_lines_parallel_vec(proj_geom):
+  v = proj_geom['Vectors']
+  for i in range(v.shape[0]):
+    ray = v[i,0:2]
+    detc = v[i,2:4]
+    detu = v[i,4:6]
+
+    detb = detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu
+
+    for i in range(proj_geom['DetectorCount']):
+      yield (detb + i * detu - ray, detb + i * detu)
+
+
+def gen_lines(proj_geom):
+  g = { 'fanflat': gen_lines_fanflat,
+        'fanflat_vec': gen_lines_fanflat_vec,
+        'parallel': gen_lines_parallel,
+        'parallel_vec': gen_lines_parallel_vec }
+  for l in g[proj_geom['type']](proj_geom):
+    yield l
+
+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())
+  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))
+  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()
+  for i in range(Vectors.shape[0]):
+    angle1 = 2*np.pi*np.random.random()
+    angle2 = angle1 + 0.5 * np.random.random()
+    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()
+  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()
+    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 ]
+    Vectors[i, :] = [ ray[0], ray[1], detc[0], detc[1], detu[0], detu[1] ]
+  pg = astra.create_proj_geom('parallel_vec', np.random.randint(*range2d), Vectors)
+  return pg
+
+
+
+
+nloops = 50
+seed = 123
+
+class TestLineKernel(unittest.TestCase):
+  def single_test(self, type):
+      shape = np.random.randint(*range2d, size=2)
+      # 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()
+          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])
+      #print(vg)
+
+      if type == 'parallel':
+        pg = gen_random_geometry_parallel()
+        projector_id = astra.create_projector('line', pg, vg)
+      elif type == 'parallel_vec':
+        pg = gen_random_geometry_parallel_vec()
+        projector_id = astra.create_projector('line', pg, vg)
+      elif type == 'fanflat':
+        pg = gen_random_geometry_fanflat()
+        projector_id = astra.create_projector('line_fanflat', pg, vg)
+      elif type == 'fanflat_vec':
+        pg = gen_random_geometry_fanflat_vec()
+        projector_id = astra.create_projector('line_fanflat', pg, vg)
+
+
+      data = np.zeros((shape[1], shape[0]), dtype=np.float32)
+      data[rect_min[1]:rect_max[1],rect_min[0]:rect_max[0]] = 1
+
+      sinogram_id, sinogram = astra.create_sino(data, projector_id)
+
+      #print(pg)
+      #print(vg)
+
+      astra.data2d.delete(sinogram_id)
+
+      astra.projector.delete(projector_id)
+
+      a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
+      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 src, det in gen_lines(pg):
+        #print(src,det)
+
+        # NB: Flipped y-axis here, since that is how astra interprets 2D volumes
+        # 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],
+                      1e-3)
+        a[i] = aw
+        b[i] = bw
+        c[i] = cw
+        i += 1
+      # Add weight for pixel / voxel size
+      try:
+        detweight = pg['DetectorWidth']
+      except KeyError:
+        detweight = np.sqrt(pg['Vectors'][0,4]*pg['Vectors'][0,4] + pg['Vectors'][0,5]*pg['Vectors'][0,5] )
+      a *= detweight
+      b *= detweight
+      c *= 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))
+
+      # 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
+      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()
+
+  def test_par(self):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test('parallel')
+  def test_fan(self):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test('fanflat')
+  def test_parvec(self):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test('parallel_vec')
+  def test_fanvec(self):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test('fanflat_vec')
+
+
+ 
+
+if __name__ == '__main__':
+  unittest.main()
+
+#print(intersect_line_rectangle((0.,-256.),(-27.,0.),11.6368454385 20.173128227 3.18989047649 5.62882841606)
diff --git a/tests/test_ParallelBeamLineKernelProjector2D.cpp b/tests/test_ParallelBeamLineKernelProjector2D.cpp
index 58d511e..575799c 100644
--- a/tests/test_ParallelBeamLineKernelProjector2D.cpp
+++ b/tests/test_ParallelBeamLineKernelProjector2D.cpp
@@ -74,7 +74,7 @@ BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_Rectangle, TestPa
 	for (int i = 0; i < iCount; ++i)
 		fWeight += pPix[i].m_fWeight;
 
-	BOOST_CHECK_SMALL(fWeight - 7.13037f, 0.00001f);
+	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
index a34c82c..c62f670 100644
--- a/tests/test_ParallelBeamLinearKernelProjector2D.cpp
+++ b/tests/test_ParallelBeamLinearKernelProjector2D.cpp
@@ -85,10 +85,10 @@ float32 compute_linear_kernel(const astra::CProjectionGeometry2D& projgeom, cons
 
 	if (fabs(cos(fAngle)) > fabs(sin(fAngle))) {
 		fDetStep = volgeom.getPixelLengthY() * fabs(cos(fAngle));
-		fWeight = volgeom.getPixelLengthX() * 1.0f / fabs(cos(fAngle));
+		fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthX() * 1.0f / fabs(cos(fAngle));
 	} else {
 		fDetStep = volgeom.getPixelLengthX() * fabs(sin(fAngle));
-		fWeight = volgeom.getPixelLengthY() * 1.0f / fabs(sin(fAngle));
+		fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthY() * 1.0f / fabs(sin(fAngle));
 	}
 
 //	printf("step: %f\n   weight: %f\n", fDetStep, fWeight);