diff options
| -rw-r--r-- | tests/python/test_line2d.py | 247 |
1 files changed, 174 insertions, 73 deletions
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py index de68033..7b46458 100644 --- a/tests/python/test_line2d.py +++ b/tests/python/test_line2d.py @@ -58,6 +58,53 @@ def intersect_line_rectangle_interval(src, det, xmin, xmax, ymin, ymax, f): c = intersect_line_rectangle_feather(src, det, xmin, xmax, ymin, ymax, f) return (a,b,c) + +# x-coord of intersection of the line (src, det) with the horizontal line at y +def intersect_line_horizontal(src, det, y): + EPS = 1e-5 + + if np.abs(src[1] - det[1]) < EPS: + return np.nan + + t = (y - src[1]) / (det[1] - src[1]) + + return src[0] + t * (det[0] - src[0]) + + +# length of the intersection of the strip with boundaries edge1, edge2 with the horizontal +# segment at y, with horizontal extent x_seg +def intersect_ray_horizontal_segment(edge1, edge2, y, x_seg): + e1 = intersect_line_horizontal(edge1[0], edge1[1], y) + e2 = intersect_line_horizontal(edge2[0], edge2[1], y) + + if not (np.isfinite(e1) and np.isfinite(e2)): + return np.nan + + (e1, e2) = np.sort([e1, e2]) + (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): + # mirror edge1 and edge2 + edge1 = [ (a[1], a[0]) for a in edge1 ] + edge2 = [ (a[1], a[0]) for a in edge2 ] + return intersect_ray_horizontal_segment(edge1, edge2, x, y_seg) + + + + + + +# LINE GENERATORS +# --------------- +# +# Per ray these yield three lines, at respectively the center and two edges of the detector pixel. +# Each line is given by two points on the line. +# ( ( (p0x, p0y), (q0x, q0y) ), ( (p1x, p1y), (q1x, q1y) ), ( (p2x, p2y), (q2x, q2y) ) ) + def gen_lines_fanflat(proj_geom): angles = proj_geom['ProjectionAngles'] for theta in angles: @@ -76,7 +123,9 @@ def gen_lines_fanflat(proj_geom): detb= detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu for i in range(proj_geom['DetectorCount']): - yield (src, detb + i * detu) + yield ((src, detb + i * detu), + (src, detb + (i - 0.5) * detu), + (src, detb + (i + 0.5) * detu)) def gen_lines_fanflat_vec(proj_geom): v = proj_geom['Vectors'] @@ -87,7 +136,9 @@ def gen_lines_fanflat_vec(proj_geom): detb = detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu for i in range(proj_geom['DetectorCount']): - yield (src, detb + i * detu) + yield ((src, detb + i * detu), + (src, detb + (i - 0.5) * detu), + (src, detb + (i + 0.5) * detu)) def gen_lines_parallel(proj_geom): angles = proj_geom['ProjectionAngles'] @@ -106,7 +157,10 @@ def gen_lines_parallel(proj_geom): 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) + yield ((detb + i * detu - ray, detb + i * detu), + (detb + (i - 0.5) * detu - ray, detb + (i - 0.5) * detu), + (detb + (i + 0.5) * detu - ray, detb + (i + 0.5) * detu)) + def gen_lines_parallel_vec(proj_geom): v = proj_geom['Vectors'] @@ -118,7 +172,9 @@ def gen_lines_parallel_vec(proj_geom): 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) + yield ((detb + i * detu - ray, detb + i * detu), + (detb + (i - 0.5) * detu - ray, detb + (i - 0.5) * detu), + (detb + (i + 0.5) * detu - ray, detb + (i + 0.5) * detu)) def gen_lines(proj_geom): @@ -179,8 +235,8 @@ def gen_random_geometry_parallel_vec(): nloops = 50 seed = 123 -class TestLineKernel(unittest.TestCase): - def single_test(self, type): +class Test2DKernel(unittest.TestCase): + def single_test(self, type, proj_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 ] @@ -201,16 +257,16 @@ class TestLineKernel(unittest.TestCase): if type == 'parallel': pg = gen_random_geometry_parallel() - projector_id = astra.create_projector('line', pg, vg) + projector_id = astra.create_projector(proj_type, pg, vg) elif type == 'parallel_vec': pg = gen_random_geometry_parallel_vec() - projector_id = astra.create_projector('line', pg, vg) + projector_id = astra.create_projector(proj_type, pg, vg) elif type == 'fanflat': pg = gen_random_geometry_fanflat() - projector_id = astra.create_projector('line_fanflat', pg, vg) + projector_id = astra.create_projector(proj_type + '_fanflat', pg, vg) elif type == 'fanflat_vec': pg = gen_random_geometry_fanflat_vec() - projector_id = astra.create_projector('line_fanflat', pg, vg) + projector_id = astra.create_projector(proj_type + '_fanflat', pg, vg) data = np.zeros((shape[1], shape[0]), dtype=np.float32) @@ -225,84 +281,129 @@ class TestLineKernel(unittest.TestCase): 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() + if proj_type == 'line': + + 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 center, edge1, edge2 in gen_lines(pg): + (src, det) = center + #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() + 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): + (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]) + 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_max[0]) * pixsize[0], + origin[0] + (-0.5 * shape[0] + rect_min[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] + a[i] = l + i += 1 + a = a.reshape(astra.functions.geom_size(pg)) + x = np.max(np.abs(sinogram-a)) + if 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) for _ in range(nloops): - self.single_test('parallel') + self.single_test('parallel', 'line') + def test_par_dd(self): + np.random.seed(seed) + for _ in range(nloops): + self.single_test('parallel', 'distance_driven') def test_fan(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('fanflat') + self.single_test('fanflat', 'line') def test_parvec(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('parallel_vec') + self.single_test('parallel_vec', 'line') + def test_parvec_dd(self): + np.random.seed(seed) + for _ in range(nloops): + self.single_test('parallel_vec', 'distance_driven') def test_fanvec(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('fanflat_vec') + self.single_test('fanflat_vec', 'line') + - if __name__ == '__main__': unittest.main() |