From 0f7c3261f159882316429d8a13009a8c1f858cbc Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Mon, 25 Mar 2019 14:51:15 +0100
Subject: Clean up projector unit tests
---
tests/python/test_line2d.py | 342 +++++++++++++++++++++++++-------------------
1 file 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)
--
cgit v1.2.3