diff options
Diffstat (limited to 'tests')
| -rw-r--r-- | tests/python/test_line2d.py | 105 |
1 files changed, 90 insertions, 15 deletions
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py index 3dd60a2..b1bb7d5 100644 --- a/tests/python/test_line2d.py +++ b/tests/python/test_line2d.py @@ -98,6 +98,30 @@ def intersect_ray_vertical_segment(edge1, edge2, x, y_seg): edge2 = [ (a[1], a[0]) for a in edge2 ] return intersect_ray_horizontal_segment(edge1, edge2, x, y_seg) +# weight of the intersection of line with the horizontal segment at y, with horizontal extent x_seg +# using linear interpolation +def intersect_line_horizontal_segment_linear(src, det, y, x_seg, inter_width): + EPS = 1e-5 + x = intersect_line_horizontal(src, det, y) + + assert(x_seg[1] - x_seg[0] + EPS >= inter_width) + if x < x_seg[0] - 0.5*inter_width: + return 0.0 + elif x < x_seg[0] + 0.5*inter_width: + return (x - (x_seg[0] - 0.5*inter_width)) / inter_width + elif x < x_seg[1] - 0.5*inter_width: + return 1.0 + elif x < x_seg[1] + 0.5*inter_width: + return (x_seg[1] + 0.5*inter_width - x) / inter_width + else: + return 0.0 + +def intersect_line_vertical_segment_linear(src, det, x, y_seg, inter_height): + src = ( src[1], src[0] ) + det = ( det[1], det[0] ) + return intersect_line_horizontal_segment_linear(src, det, x, y_seg, inter_height) + + def area_signed(a, b): return a[0] * b[1] - a[1] * b[0] @@ -286,6 +310,9 @@ def gen_random_geometry_parallel_vec(): nloops = 50 seed = 123 +def proj_type_to_fan(t): + return t + '_fanflat' + class Test2DKernel(unittest.TestCase): def single_test(self, type, proj_type): shape = np.random.randint(*range2d, size=2) @@ -314,10 +341,10 @@ class Test2DKernel(unittest.TestCase): projector_id = astra.create_projector(proj_type, pg, vg) elif type == 'fanflat': pg = gen_random_geometry_fanflat() - projector_id = astra.create_projector(proj_type + '_fanflat', pg, vg) + 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 + '_fanflat', pg, vg) + projector_id = astra.create_projector(proj_type_to_fan(proj_type), pg, vg) data = np.zeros((shape[1], shape[0]), dtype=np.float32) @@ -332,6 +359,12 @@ 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) + if proj_type == 'line': a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) @@ -358,11 +391,6 @@ class Test2DKernel(unittest.TestCase): 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 @@ -393,6 +421,44 @@ class Test2DKernel(unittest.TestCase): pylab.figure() pylab.imshow(c-sinogram) pylab.show() + elif proj_type == 'linear': + 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] + 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]) + 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]) + 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]) + 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]) + l += w * length * pixsize[1] * detweight + a[i] = l + i += 1 + a = a.reshape(astra.functions.geom_size(pg)) + 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() + self.assertFalse(x > TOL) elif proj_type == 'distance_driven': a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) i = 0 @@ -406,8 +472,8 @@ class Test2DKernel(unittest.TestCase): 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_max[0]) * pixsize[0], - origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0]) + x_seg = (origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], + origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0]) 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] @@ -415,7 +481,8 @@ class Test2DKernel(unittest.TestCase): i += 1 a = a.reshape(astra.functions.geom_size(pg)) x = np.max(np.abs(sinogram-a)) - if x > 2e-3: + TOL = 2e-3 + if x > TOL: pylab.gray() pylab.imshow(data) pylab.figure() @@ -425,7 +492,7 @@ class Test2DKernel(unittest.TestCase): pylab.figure() pylab.imshow(sinogram-a) pylab.show() - self.assertFalse(x > 2e-3) + 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] @@ -440,9 +507,9 @@ class Test2DKernel(unittest.TestCase): 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: + # TODO: Investigate tolerance for fanflat/strip + TOL = 2e-3 + if False and x > TOL: pylab.gray() pylab.imshow(data) pylab.figure() @@ -452,12 +519,16 @@ class Test2DKernel(unittest.TestCase): pylab.figure() pylab.imshow(sinogram-a) pylab.show() - self.assertFalse(x > 2e-3) + self.assertFalse(x > TOL) def test_par(self): np.random.seed(seed) for _ in range(nloops): self.single_test('parallel', 'line') + def test_par_linear(self): + np.random.seed(seed) + for _ in range(nloops): + self.single_test('parallel', 'linear') def test_par_dd(self): np.random.seed(seed) for _ in range(nloops): @@ -478,6 +549,10 @@ class Test2DKernel(unittest.TestCase): np.random.seed(seed) for _ in range(nloops): self.single_test('parallel_vec', 'line') + def test_parvec_linear(self): + np.random.seed(seed) + for _ in range(nloops): + self.single_test('parallel_vec', 'linear') def test_parvec_dd(self): np.random.seed(seed) for _ in range(nloops): |