diff options
| -rw-r--r-- | tests/python/test_line2d.py | 342 |
1 files changed, 191 insertions, 151 deletions
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py index b1bb7d5..e7bca98 100644 --- a/tests/python/test_line2d.py +++ b/tests/python/test_line2d.py @@ -4,10 +4,35 @@ import astra import math import pylab +# Display sinograms with mismatch on test failure +DISPLAY=False + +NONUNITDET=False +OBLIQUE=False +FLEXVOL=False +NONSQUARE=False # non-square pixels not supported yet by most projectors + +# Round interpolation weight to 8 bits to emulate CUDA texture unit precision +CUDA_8BIT_LINEAR=True +CUDA_TOL=2e-2 + +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) + + + # 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 @@ -89,7 +114,6 @@ def intersect_ray_horizontal_segment(edge1, edge2, y, x_seg): (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): @@ -128,7 +152,8 @@ def area_signed(a, b): # 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 + EPS = 1e-5 + return area_signed( (b[0] - a[0], b[1] - a[1]), (c[0] - a[0], c[1] - a[1]) ) > EPS # compute area of rect on left side of line def halfarea_rect_line(src, det, xmin, xmax, ymin, ymax): @@ -172,6 +197,13 @@ def intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax): return abs(s1 - s2) +# width of projection of detector orthogonal to ray direction +# i.e., effective detector width +def effective_detweight(src, det, u): + ray = np.array(det) - np.array(src) + ray = ray / np.linalg.norm(ray, ord=2) + return abs(area_signed(ray, u)) + # LINE GENERATORS # --------------- @@ -264,39 +296,56 @@ 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()) + if not NONUNITDET: + w = 1.0 + else: + w = 0.6 + 0.8 * np.random.random() + pg = astra.create_proj_geom('fanflat', w, 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)) + if not NONUNITDET: + w = 1.0 + else: + w = 0.8 + 0.4 * np.random.random() + pg = astra.create_proj_geom('parallel', w, 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() + if not NONUNITDET: + w = 1.0 + else: + 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() + if OBLIQUE: + angle2 = angle1 + 0.5 * np.random.random() + else: + angle2 = angle1 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() + if not NONUNITDET: + w = 1.0 + else: + 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() + if OBLIQUE: + angle2 = angle1 + 0.5 * np.random.random() + else: + angle2 = angle1 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 ] @@ -307,11 +356,39 @@ def gen_random_geometry_parallel_vec(): -nloops = 50 -seed = 123 - def proj_type_to_fan(t): - return t + '_fanflat' + if t == 'cuda': + return t + else: + return t + '_fanflat' + +def display_mismatch(data, sinogram, a): + pylab.gray() + pylab.imshow(data) + pylab.figure() + pylab.imshow(sinogram) + pylab.figure() + pylab.imshow(a) + pylab.figure() + pylab.imshow(sinogram-a) + pylab.show() + +def display_mismatch_triple(data, sinogram, a, b, c): + 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() class Test2DKernel(unittest.TestCase): def single_test(self, type, proj_type): @@ -319,9 +396,11 @@ class Test2DKernel(unittest.TestCase): # 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() + 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.) @@ -331,7 +410,6 @@ class Test2DKernel(unittest.TestCase): 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() @@ -352,6 +430,8 @@ class Test2DKernel(unittest.TestCase): sinogram_id, sinogram = astra.create_sino(data, projector_id) + self.assertTrue(np.all(np.isfinite(sinogram))) + #print(pg) #print(vg) @@ -359,11 +439,11 @@ 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) + # NB: Flipped y-axis here, since that is how astra interprets 2D volumes + 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] if proj_type == 'line': @@ -371,200 +451,161 @@ class Test2DKernel(unittest.TestCase): 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): + for i, (center, edge1, edge2) in enumerate(gen_lines(pg)): (src, det) = center - #print(src,det) - # NB: Flipped y-axis here, since that is how astra interprets 2D volumes + 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, - 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], + xmin, xmax, ymin, ymax, 1e-3) - a[i] = aw - b[i] = bw - c[i] = cw - i += 1 - a *= detweight - b *= detweight - c *= detweight + a[i] = aw * detweight + b[i] = bw * detweight + c[i] = cw * 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)) + if not np.all(np.isfinite(a)): + raise RuntimeError("Invalid value in reference sinogram") + if not np.all(np.isfinite(b)): + raise RuntimeError("Invalid value in reference sinogram") + if not np.all(np.isfinite(c)): + raise RuntimeError("Invalid value in reference sinogram") + self.assertTrue(np.all(np.isfinite(sinogram))) + # 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 + if DISPLAY and (z < 0 or y < 0): + display_mismatch_triple(data, sinogram, a, b, c) 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 == 'linear': + elif proj_type == 'linear' or proj_type == 'cuda': 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] + 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) - y_seg = (origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1], - origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1]) + 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] w = intersect_line_vertical_segment_linear(center[0], center[1], x, y_seg, pixsize[1]) + # 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 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]) + 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] w = intersect_line_horizontal_segment_linear(center[0], center[1], y, x_seg, pixsize[0]) + # 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 a[i] = l - i += 1 a = a.reshape(astra.functions.geom_size(pg)) + if not np.all(np.isfinite(a)): + raise RuntimeError("Invalid value in reference sinogram") 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() + TOL = 2e-3 if proj_type != 'cuda' else CUDA_TOL + if DISPLAY and x > TOL: + display_mismatch(data, sinogram, a) self.assertFalse(x > TOL) 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): + for i, (center, edge1, edge2) in enumerate(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]) + 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] else: - x_seg = (origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], - origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0]) + 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] a[i] = l - i += 1 a = a.reshape(astra.functions.geom_size(pg)) + if not np.all(np.isfinite(a)): + raise RuntimeError("Invalid value in reference sinogram") x = np.max(np.abs(sinogram-a)) TOL = 2e-3 - if 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() + if DISPLAY and x > TOL: + display_mismatch(data, sinogram, a) 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] - 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): + for i, (center, edge1, edge2) in enumerate(gen_lines(pg)): a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax) - i += 1 a = a.reshape(astra.functions.geom_size(pg)) + if not np.all(np.isfinite(a)): + raise RuntimeError("Invalid value in reference sinogram") x = np.max(np.abs(sinogram-a)) - # TODO: Investigate tolerance for fanflat/strip - TOL = 2e-3 - if False 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() + TOL = 8e-3 + if DISPLAY and x > TOL: + display_mismatch(data, sinogram, a) self.assertFalse(x > TOL) - def test_par(self): + def multi_test(self, type, proj_type): np.random.seed(seed) for _ in range(nloops): - self.single_test('parallel', 'line') + self.single_test(type, proj_type) + + def test_par(self): + self.multi_test('parallel', 'line') def test_par_linear(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel', 'linear') + self.multi_test('parallel', 'linear') + def test_par_cuda(self): + self.multi_test('parallel', 'cuda') def test_par_dd(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel', 'distance_driven') + self.multi_test('parallel', 'distance_driven') def test_par_strip(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel', 'strip') + self.multi_test('parallel', 'strip') def test_fan(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('fanflat', 'line') + self.multi_test('fanflat', 'line') def test_fan_strip(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('fanflat', 'strip') + self.multi_test('fanflat', 'strip') + def test_fan_cuda(self): + self.multi_test('fanflat', 'cuda') def test_parvec(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel_vec', 'line') + self.multi_test('parallel_vec', 'line') def test_parvec_linear(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel_vec', 'linear') + self.multi_test('parallel_vec', 'linear') def test_parvec_dd(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel_vec', 'distance_driven') + self.multi_test('parallel_vec', 'distance_driven') def test_parvec_strip(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('parallel_vec', 'strip') + self.multi_test('parallel_vec', 'strip') + def test_parvec_cuda(self): + self.multi_test('parallel_vec', 'cuda') def test_fanvec(self): - np.random.seed(seed) - for _ in range(nloops): - self.single_test('fanflat_vec', 'line') + self.multi_test('fanflat_vec', 'line') + def test_fanvec_cuda(self): + self.multi_test('fanflat_vec', 'cuda') + @@ -572,4 +613,3 @@ class Test2DKernel(unittest.TestCase): if __name__ == '__main__': unittest.main() -#print(intersect_line_rectangle((0.,-256.),(-27.,0.),11.6368454385 20.173128227 3.18989047649 5.62882841606) |