summaryrefslogtreecommitdiffstats
path: root/tests
diff options
context:
space:
mode:
Diffstat (limited to 'tests')
-rw-r--r--tests/python/test_line2d.py89
1 files changed, 89 insertions, 0 deletions
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 7b46458..3dd60a2 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -70,6 +70,11 @@ def intersect_line_horizontal(src, det, y):
return src[0] + t * (det[0] - src[0])
+# y-coord of intersection of the line (src, det) with the vertical line at x
+def intersect_line_vertical(src, det, x):
+ src = ( src[1], src[0] )
+ det = ( det[1], det[0] )
+ return intersect_line_horizontal(src, det, x)
# length of the intersection of the strip with boundaries edge1, edge2 with the horizontal
# segment at y, with horizontal extent x_seg
@@ -94,7 +99,53 @@ def intersect_ray_vertical_segment(edge1, edge2, x, y_seg):
return intersect_ray_horizontal_segment(edge1, edge2, x, y_seg)
+def area_signed(a, b):
+ return a[0] * b[1] - a[1] * b[0]
+# 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
+
+# compute area of rect on left side of line
+def halfarea_rect_line(src, det, xmin, xmax, ymin, ymax):
+ pts = ( (xmin,ymin), (xmin,ymax), (xmax,ymin), (xmax,ymax) )
+ pts_left = list(filter( lambda p: is_left_of(src, det, p), pts ))
+ npts_left = len(pts_left)
+ if npts_left == 0:
+ return 0.0
+ elif npts_left == 1:
+ # triangle
+ p = pts_left[0]
+ xd = intersect_line_horizontal(src, det, p[1]) - p[0]
+ yd = intersect_line_vertical(src, det, p[0]) - p[1]
+ ret = 0.5 * abs(xd) * abs(yd)
+ return ret
+ elif npts_left == 2:
+ p = pts_left[0]
+ q = pts_left[1]
+ if p[0] == q[0]:
+ # vertical intersection
+ x1 = intersect_line_horizontal(src, det, p[1]) - p[0]
+ x2 = intersect_line_horizontal(src, det, q[1]) - q[0]
+ ret = 0.5 * (ymax - ymin) * (abs(x1) + abs(x2))
+ return ret
+ else:
+ assert(p[1] == q[1])
+ # horizontal intersection
+ y1 = intersect_line_vertical(src, det, p[0]) - p[1]
+ y2 = intersect_line_vertical(src, det, q[0]) - q[1]
+ ret = 0.5 * (xmax - xmin) * (abs(y1) + abs(y2))
+ return ret
+ else:
+ # mirror and invert
+ ret = ((xmax - xmin) * (ymax - ymin)) - halfarea_rect_line(det, src, xmin, xmax, ymin, ymax)
+ return ret
+
+# area of intersection of the strip with boundaries edge1, edge2 with rectangle
+def intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax):
+ s1 = halfarea_rect_line(edge1[0], edge1[1], xmin, xmax, ymin, ymax)
+ s2 = halfarea_rect_line(edge2[0], edge2[1], xmin, xmax, ymin, ymax)
+ return abs(s1 - s2)
@@ -375,7 +426,33 @@ class Test2DKernel(unittest.TestCase):
pylab.imshow(sinogram-a)
pylab.show()
self.assertFalse(x > 2e-3)
+ 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):
+ a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax)
+ i += 1
+
+ a = a.reshape(astra.functions.geom_size(pg))
+ x = np.max(np.abs(sinogram-a))
+ if x > 2e-3:
+ print(x)
+ if False and x > 2e-3:
+ pylab.gray()
+ pylab.imshow(data)
+ pylab.figure()
+ pylab.imshow(sinogram)
+ pylab.figure()
+ pylab.imshow(a)
+ pylab.figure()
+ pylab.imshow(sinogram-a)
+ pylab.show()
+ self.assertFalse(x > 2e-3)
def test_par(self):
np.random.seed(seed)
@@ -385,10 +462,18 @@ class Test2DKernel(unittest.TestCase):
np.random.seed(seed)
for _ in range(nloops):
self.single_test('parallel', 'distance_driven')
+ def test_par_strip(self):
+ np.random.seed(seed)
+ for _ in range(nloops):
+ self.single_test('parallel', 'strip')
def test_fan(self):
np.random.seed(seed)
for _ in range(nloops):
self.single_test('fanflat', 'line')
+ def test_fan_strip(self):
+ np.random.seed(seed)
+ for _ in range(nloops):
+ self.single_test('fanflat', 'strip')
def test_parvec(self):
np.random.seed(seed)
for _ in range(nloops):
@@ -397,6 +482,10 @@ class Test2DKernel(unittest.TestCase):
np.random.seed(seed)
for _ in range(nloops):
self.single_test('parallel_vec', 'distance_driven')
+ def test_parvec_strip(self):
+ np.random.seed(seed)
+ for _ in range(nloops):
+ self.single_test('parallel_vec', 'strip')
def test_fanvec(self):
np.random.seed(seed)
for _ in range(nloops):