diff options
Diffstat (limited to 'tests/python')
| -rw-r--r-- | tests/python/test_line2d.py | 180 | ||||
| -rw-r--r-- | tests/python/test_rec_scaling.py | 213 |
2 files changed, 316 insertions, 77 deletions
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py index d04ffb8..e5d8f2b 100644 --- a/tests/python/test_line2d.py +++ b/tests/python/test_line2d.py @@ -7,9 +7,9 @@ import pylab # Display sinograms with mismatch on test failure DISPLAY=False -NONUNITDET=False -OBLIQUE=False -FLEXVOL=False +NONUNITDET=True +OBLIQUE=True +FLEXVOL=True NONSQUARE=False # non-square pixels not supported yet by most projectors # Round interpolation weight to 8 bits to emulate CUDA texture unit precision @@ -20,15 +20,8 @@ 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) - +# KNOWN FAILURES: +# fan/strip relatively high numerical errors around 45 degrees # return length of intersection of the line through points src = (x,y) @@ -454,23 +447,15 @@ class Test2DKernel(unittest.TestCase): for i, (center, edge1, edge2) in enumerate(gen_lines(pg)): (src, det) = center - 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, xmin, xmax, ymin, ymax, 1e-3) - a[i] = aw * detweight - b[i] = bw * detweight - c[i] = cw * detweight + a[i] = aw + b[i] = bw + c[i] = cw a = a.reshape(astra.functions.geom_size(pg)) b = b.reshape(astra.functions.geom_size(pg)) c = c.reshape(astra.functions.geom_size(pg)) @@ -494,17 +479,9 @@ class Test2DKernel(unittest.TestCase): 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) + length = math.sqrt(1.0 + abs(yd/xd)**2) * pixsize[0] 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] @@ -512,9 +489,9 @@ class Test2DKernel(unittest.TestCase): # 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 + l += w * length else: - length = math.sqrt(1.0 + abs(xd/yd)**2) + length = math.sqrt(1.0 + abs(xd/yd)**2) * pixsize[1] 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] @@ -522,7 +499,7 @@ class Test2DKernel(unittest.TestCase): # 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 + l += w * length a[i] = l a = a.reshape(astra.functions.geom_size(pg)) if not np.all(np.isfinite(a)): @@ -532,21 +509,26 @@ class Test2DKernel(unittest.TestCase): if DISPLAY and x > TOL: display_mismatch(data, sinogram, a) self.assertFalse(x > TOL) - elif proj_type == 'distance_driven': + elif proj_type == 'distance_driven' and 'par' in type: a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) for i, (center, edge1, edge2) in enumerate(gen_lines(pg)): - (xd, yd) = center[1] - center[0] + (src, det) = center + try: + detweight = pg['DetectorWidth'] + except KeyError: + detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6]) + (xd, yd) = det - src l = 0.0 if np.abs(xd) > np.abs(yd): # horizontal ray 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] + l += intersect_ray_vertical_segment(edge1, edge2, x, y_seg) * pixsize[0] / detweight else: 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] + l += intersect_ray_horizontal_segment(edge1, edge2, y, x_seg) * pixsize[1] / detweight a[i] = l a = a.reshape(astra.functions.geom_size(pg)) if not np.all(np.isfinite(a)): @@ -560,6 +542,7 @@ class Test2DKernel(unittest.TestCase): a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) for i, (center, edge1, edge2) in enumerate(gen_lines(pg)): (src, det) = center + detweight = effective_detweight(src, det, edge2[1] - edge1[1]) det_dist = np.linalg.norm(src-det, ord=2) l = 0.0 for j in range(rect_min[0], rect_max[0]): @@ -570,7 +553,7 @@ class Test2DKernel(unittest.TestCase): ymin = origin[1] + (+0.5 * shape[1] - k - 1) * pixsize[1] ymax = origin[1] + (+0.5 * shape[1] - k) * pixsize[1] ycen = 0.5 * (ymin + ymax) - scale = det_dist / np.linalg.norm( src - np.array((xcen,ycen)), ord=2 ) + scale = det_dist / (np.linalg.norm( src - np.array((xcen,ycen)), ord=2 ) * detweight) w = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax) l += w * scale a[i] = l @@ -578,14 +561,20 @@ class Test2DKernel(unittest.TestCase): if not np.all(np.isfinite(a)): raise RuntimeError("Invalid value in reference sinogram") x = np.max(np.abs(sinogram-a)) - TOL = 8e-3 + # BUG: Known bug in fan/strip code around 45 degree projections causing larger errors than desirable + TOL = 4e-2 if DISPLAY and x > TOL: display_mismatch(data, sinogram, a) self.assertFalse(x > TOL) elif proj_type == 'strip': a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) for i, (center, edge1, edge2) in enumerate(gen_lines(pg)): - a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax) + (src, det) = center + try: + detweight = pg['DetectorWidth'] + except KeyError: + detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6]) + a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax) / detweight a = a.reshape(astra.functions.geom_size(pg)) if not np.all(np.isfinite(a)): raise RuntimeError("Invalid value in reference sinogram") @@ -594,46 +583,83 @@ class Test2DKernel(unittest.TestCase): if DISPLAY and x > TOL: display_mismatch(data, sinogram, a) self.assertFalse(x > TOL) + else: + raise RuntimeError("Unsupported projector") - def multi_test(self, type, proj_type): - np.random.seed(seed) - for _ in range(nloops): - self.single_test(type, proj_type) - - def test_par(self): - self.multi_test('parallel', 'line') - def test_par_linear(self): - self.multi_test('parallel', 'linear') - def test_par_cuda(self): - self.multi_test('parallel', 'cuda') - def test_par_dd(self): - self.multi_test('parallel', 'distance_driven') - def test_par_strip(self): - self.multi_test('parallel', 'strip') - def test_fan(self): - self.multi_test('fanflat', 'line') - def test_fan_strip(self): - self.multi_test('fanflat', 'strip') - def test_fan_cuda(self): - self.multi_test('fanflat', 'cuda') - def test_parvec(self): - self.multi_test('parallel_vec', 'line') - def test_parvec_linear(self): - self.multi_test('parallel_vec', 'linear') - def test_parvec_dd(self): - self.multi_test('parallel_vec', 'distance_driven') - def test_parvec_strip(self): - self.multi_test('parallel_vec', 'strip') - def test_parvec_cuda(self): - self.multi_test('parallel_vec', 'cuda') - def test_fanvec(self): - self.multi_test('fanflat_vec', 'line') - def test_fanvec_cuda(self): - self.multi_test('fanflat_vec', 'cuda') + def single_test_adjoint(self, type, proj_type): + shape = np.random.randint(*range2d, size=2) + 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.) + 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]) + if type == 'parallel': + pg = gen_random_geometry_parallel() + projector_id = astra.create_projector(proj_type, pg, vg) + elif type == 'parallel_vec': + pg = gen_random_geometry_parallel_vec() + projector_id = astra.create_projector(proj_type, pg, vg) + elif type == 'fanflat': + pg = gen_random_geometry_fanflat() + 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_to_fan(proj_type), pg, vg) + for i in range(5): + X = np.random.random((shape[1], shape[0])) + Y = np.random.random(astra.geom_size(pg)) + + sinogram_id, fX = astra.create_sino(X, projector_id) + bp_id, fTY = astra.create_backprojection(Y, projector_id) + + astra.data2d.delete(sinogram_id) + astra.data2d.delete(bp_id) + + da = np.dot(fX.ravel(), Y.ravel()) + db = np.dot(X.ravel(), fTY.ravel()) + m = np.abs(da - db) + TOL = 1e-3 if 'cuda' not in proj_type else 1e-1 + if m / da >= TOL: + print(vg) + print(pg) + print(m/da, da/db, da, db) + self.assertTrue(m / da < TOL) + astra.projector.delete(projector_id) + def multi_test(self, type, proj_type): + np.random.seed(seed) + for _ in range(nloops): + self.single_test(type, proj_type) + def multi_test_adjoint(self, type, proj_type): + np.random.seed(seed) + for _ in range(nloops): + self.single_test_adjoint(type, proj_type) + +__combinations = { 'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ], + 'parallel_vec': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ], + 'fanflat': [ 'line', 'strip', 'cuda' ], + 'fanflat_vec': [ 'line', 'cuda' ] } + +for k, l in __combinations.items(): + for v in l: + def f(k,v): + return lambda self: self.multi_test(k, v) + def f_adj(k,v): + return lambda self: self.multi_test_adjoint(k, v) + setattr(Test2DKernel, 'test_' + k + '_' + v, f(k,v)) + setattr(Test2DKernel, 'test_' + k + '_' + v + '_adjoint', f_adj(k,v)) if __name__ == '__main__': unittest.main() diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py new file mode 100644 index 0000000..621fd8a --- /dev/null +++ b/tests/python/test_rec_scaling.py @@ -0,0 +1,213 @@ +import numpy as np +import unittest +import astra +import math +import pylab + +DISPLAY=False + +def VolumeGeometries(is3D,noncube): + if not is3D: + for s in [0.8, 1.0, 1.25]: + yield astra.create_vol_geom(128, 128, -64*s, 64*s, -64*s, 64*s) + elif noncube: + for sx in [0.8, 1.0]: + for sy in [0.8, 1.0]: + for sz in [0.8, 1.0]: + yield astra.create_vol_geom(64, 64, 64, -32*sx, 32*sx, -32*sy, 32*sy, -32*sz, 32*sz) + else: + for s in [0.8, 1.0]: + yield astra.create_vol_geom(64, 64, 64, -32*s, 32*s, -32*s, 32*s, -32*s, 32*s) + + +def ProjectionGeometries(type): + if type == 'parallel': + for dU in [0.8, 1.0, 1.25]: + yield astra.create_proj_geom('parallel', dU, 256, np.linspace(0,np.pi,180,False)) + elif type == 'fanflat': + for dU in [0.8, 1.0, 1.25]: + for src in [500, 1000]: + for det in [0, 250, 500]: + yield astra.create_proj_geom('fanflat', dU, 256, np.linspace(0,2*np.pi,180,False), src, det) + elif type == 'parallel3d': + for dU in [0.8, 1.0]: + for dV in [0.8, 1.0]: + yield astra.create_proj_geom('parallel3d', dU, dV, 128, 128, np.linspace(0,np.pi,180,False)) + elif type == 'parallel3d_vec': + for j in range(10): + Vectors = np.zeros([180,12]) + wu = 0.6 + 0.8 * np.random.random() + wv = 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() + angle3 = 0.1*np.pi*np.random.random() + detc = 10 * np.random.random(size=3) + detu = [ math.cos(angle1) * wu, math.sin(angle1) * wu, 0 ] + detv = [ -math.sin(angle1) * math.sin(angle3) * wv, math.cos(angle1) * math.sin(angle3) * wv, math.cos(angle3) * wv ] + ray = [ math.sin(angle2) * l, -math.cos(angle2) * l, 0 ] + Vectors[i, :] = [ ray[0], ray[1], ray[2], detc[0], detc[1], detc[2], detu[0], detu[1], detu[2], detv[0], detv[1], detv[2] ] + pg = astra.create_proj_geom('parallel3d_vec', 128, 128, Vectors) + yield pg + elif type == 'cone': + for dU in [0.8, 1.0]: + for dV in [0.8, 1.0]: + for src in [500, 1000]: + for det in [0, 250]: + yield astra.create_proj_geom('cone', dU, dV, 128, 128, np.linspace(0,2*np.pi,180,False), src, det) + elif type == 'cone_vec': + for j in range(10): + Vectors = np.zeros([180,12]) + wu = 0.6 + 0.8 * np.random.random() + wv = 0.6 + 0.8 * np.random.random() + for i in range(Vectors.shape[0]): + l = 256 * (0.5 * np.random.random()) + angle1 = 2*np.pi*np.random.random() + angle2 = angle1 + 0.5 * np.random.random() + angle3 = 0.1*np.pi*np.random.random() + detc = 10 * np.random.random(size=3) + detu = [ math.cos(angle1) * wu, math.sin(angle1) * wu, 0 ] + detv = [ -math.sin(angle1) * math.sin(angle3) * wv, math.cos(angle1) * math.sin(angle3) * wv, math.cos(angle3) * wv ] + src = [ math.sin(angle2) * l, -math.cos(angle2) * l, 0 ] + Vectors[i, :] = [ src[0], src[1], src[2], detc[0], detc[1], detc[2], detu[0], detu[1], detu[2], detv[0], detv[1], detv[2] ] + pg = astra.create_proj_geom('parallel3d_vec', 128, 128, Vectors) + yield pg + + +class TestRecScale(unittest.TestCase): + def single_test(self, geom_type, proj_type, alg, iters): + if alg == 'FBP' and 'fanflat' in geom_type: + self.skipTest('CPU FBP is parallel-beam only') + is3D = (geom_type in ['parallel3d', 'cone']) + for vg in VolumeGeometries(is3D, 'FDK' not in alg): + for pg in ProjectionGeometries(geom_type): + if not is3D: + vol = np.zeros((128,128),dtype=np.float32) + vol[50:70,50:70] = 1 + else: + vol = np.zeros((64,64,64),dtype=np.float32) + vol[25:35,25:35,25:35] = 1 + proj_id = astra.create_projector(proj_type, pg, vg) + if not is3D: + sino_id, sinogram = astra.create_sino(vol, proj_id) + else: + sino_id, sinogram = astra.create_sino3d_gpu(vol, pg, vg) + if not is3D: + DATA = astra.data2d + else: + DATA = astra.data3d + + rec_id = DATA.create('-vol', vg, 0.0 if 'EM' not in alg else 1.0) + + cfg = astra.astra_dict(alg) + cfg['ReconstructionDataId'] = rec_id + cfg['ProjectionDataId'] = sino_id + cfg['ProjectorId'] = proj_id + alg_id = astra.algorithm.create(cfg) + + for i in range(iters): + astra.algorithm.run(alg_id, 1) + rec = DATA.get(rec_id) + astra.astra.delete([sino_id, alg_id, alg_id, proj_id]) + if not is3D: + val = np.sum(rec[55:65,55:65]) / 100. + else: + val = np.sum(rec[27:32,27:32,27:32]) / 125. + TOL = 5e-2 + if DISPLAY and abs(val-1.0) >= TOL: + print(geom_type, proj_type, alg, vg, pg) + print(val) + pylab.gray() + if not is3D: + pylab.imshow(rec) + else: + pylab.imshow(rec[:,32,:]) + pylab.show() + self.assertTrue(abs(val-1.0) < TOL) + + def single_test_adjoint3D(self, geom_type, proj_type): + for vg in VolumeGeometries(True, True): + for pg in ProjectionGeometries(geom_type): + for i in range(5): + X = np.random.random(astra.geom_size(vg)) + Y = np.random.random(astra.geom_size(pg)) + proj_id, fX = astra.create_sino3d_gpu(X, pg, vg) + bp_id, fTY = astra.create_backprojection3d_gpu(Y, pg, vg) + + astra.data3d.delete([proj_id, bp_id]) + + da = np.dot(fX.ravel(), Y.ravel()) + db = np.dot(X.ravel(), fTY.ravel()) + m = np.abs(da - db) + TOL = 1e-1 + if m / da >= TOL: + print(vg) + print(pg) + print(m/da, da/db, da, db) + self.assertTrue(m / da < TOL) + + + + + +__combinations = { + 'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ], + 'fanflat': [ 'line_fanflat', 'strip_fanflat', 'cuda' ], + 'parallel3d': [ 'cuda3d' ], + 'cone': [ 'cuda3d' ], +} + +__combinations_adjoint = { + 'parallel3d': [ 'cuda3d' ], + 'cone': [ 'cuda3d' ], + 'parallel3d_vec': [ 'cuda3d' ], + 'cone_vec': [ 'cuda3d' ], +} + +__algs = { + 'SIRT': 50, 'SART': 10*180, 'CGLS': 30, + 'FBP': 1 +} + +__algs_CUDA = { + 'SIRT_CUDA': 50, 'SART_CUDA': 10*180, 'CGLS_CUDA': 30, 'EM_CUDA': 50, + 'FBP_CUDA': 1 +} + +__algs_parallel3d = { + 'SIRT3D_CUDA': 200, 'CGLS3D_CUDA': 20, +} + +__algs_cone = { + 'SIRT3D_CUDA': 200, 'CGLS3D_CUDA': 20, + 'FDK_CUDA': 1 +} + + + +for k, l in __combinations.items(): + for v in l: + is3D = (k in ['parallel3d', 'cone']) + if k == 'parallel3d': + A = __algs_parallel3d + elif k == 'cone': + A = __algs_cone + elif v == 'cuda': + A = __algs_CUDA + else: + A = __algs + for a, i in A.items(): + def f(k, v, a, i): + return lambda self: self.single_test(k, v, a, i) + setattr(TestRecScale, 'test_' + a + '_' + k + '_' + v, f(k,v,a,i)) + +for k, l in __combinations_adjoint.items(): + for v in l: + def g(k, v): + return lambda self: self.single_test_adjoint3D(k, v) + setattr(TestRecScale, 'test_adjoint_' + k + '_' + v, g(k,v)) + +if __name__ == '__main__': + unittest.main() + |