From 0f4ceb4c7f3f63fddf8fbf44c59fcd8f415e3847 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 29 Mar 2019 14:42:53 +0100
Subject: Adjust line kernels to line integral scaling
---
tests/python/test_line2d.py | 14 +++-----------
1 file changed, 3 insertions(+), 11 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index d04ffb8..755bd27 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -454,23 +454,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))
--
cgit v1.2.3
From 35957b6ef72749cdc520ded67a0eb8cdfd7ea655 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 29 Mar 2019 15:03:57 +0100
Subject: Adjust linear/cuda kernels to line integral scaling
---
tests/python/test_line2d.py | 16 ++++------------
1 file changed, 4 insertions(+), 12 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 755bd27..5647053 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -486,17 +486,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]
@@ -504,9 +496,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]
@@ -514,7 +506,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)):
--
cgit v1.2.3
From 15f53527e2f38e8eeb8bece79376b5e4686f3dd4 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 29 Mar 2019 16:01:02 +0100
Subject: Adjust distance driven kernels to line integral scaling
---
tests/python/test_line2d.py | 15 +++++++++++----
1 file changed, 11 insertions(+), 4 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 5647053..3019277 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -516,21 +516,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)):
@@ -578,6 +583,8 @@ 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)
--
cgit v1.2.3
From 9c869d834ddc681299df9999787eb92bd6010dac Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <wjp@usecode.org>
Date: Fri, 29 Mar 2019 19:08:11 +0100
Subject: Adjust strip kernels to line integral scaling
---
tests/python/test_line2d.py | 24 ++++++++++++------------
1 file changed, 12 insertions(+), 12 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 3019277..ba3f7ea 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -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)
@@ -549,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]):
@@ -559,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
@@ -567,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")
--
cgit v1.2.3
From 87715885f3b4a80693493e37aa8293899a6b987e Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <wjp@usecode.org>
Date: Fri, 29 Mar 2019 21:18:10 +0100
Subject: Dynamically create python test functions
---
tests/python/test_line2d.py | 44 ++++++++++----------------------------------
1 file changed, 10 insertions(+), 34 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index ba3f7ea..4c3d174 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -591,40 +591,16 @@ class Test2DKernel(unittest.TestCase):
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')
-
-
-
-
+__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)
+ setattr(Test2DKernel, 'test_' + k + '_' + v, f(k,v))
if __name__ == '__main__':
unittest.main()
--
cgit v1.2.3
From 3cf63d335ebe392a8c77f0c90395c18150647aeb Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <wjp@usecode.org>
Date: Fri, 29 Mar 2019 21:21:29 +0100
Subject: Adjust adjoint to line integral scaling
---
tests/python/test_line2d.py | 59 +++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 59 insertions(+)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 4c3d174..8b6e200 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -586,10 +586,66 @@ class Test2DKernel(unittest.TestCase):
else:
raise RuntimeError("Unsupported projector")
+ 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' ],
@@ -600,7 +656,10 @@ 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()
--
cgit v1.2.3
From 35052afe6c27119b7d1d58a035d17be043d686f2 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Tue, 2 Apr 2019 15:58:14 +0200
Subject: Add test for reconstruction scaling
---
tests/python/test_rec_scaling.py | 79 ++++++++++++++++++++++++++++++++++++++++
1 file changed, 79 insertions(+)
create mode 100644 tests/python/test_rec_scaling.py
(limited to 'tests/python')
diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py
new file mode 100644
index 0000000..33d09f9
--- /dev/null
+++ b/tests/python/test_rec_scaling.py
@@ -0,0 +1,79 @@
+import numpy as np
+import unittest
+import astra
+import math
+import pylab
+
+DISPLAY=False
+
+def VolumeGeometries():
+ for s in [0.8, 1.0, 1.25]:
+ yield astra.create_vol_geom(128, 128, -64*s, 64*s, -64*s, 64*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)
+
+
+class Test2DRecScale(unittest.TestCase):
+ def single_test(self, geom_type, proj_type, alg, iters):
+ for vg in VolumeGeometries():
+ for pg in ProjectionGeometries(geom_type):
+ vol = np.zeros((128,128))
+ vol[50:70,50:70] = 1
+ proj_id = astra.create_projector(proj_type, pg, vg)
+ sino_id, sinogram = astra.create_sino(vol, proj_id)
+ rec_id = astra.data2d.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)
+
+ astra.algorithm.run(alg_id, iters)
+ rec = astra.data2d.get(rec_id)
+ astra.astra.delete([sino_id, alg_id, alg_id, proj_id])
+ val = np.sum(rec[55:65,55:65]) / 100.
+ TOL = 5e-2
+ if DISPLAY and abs(val-1.0) >= TOL:
+ print(geom_type, proj_type, alg, vg, pg)
+ print(val)
+ pylab.gray()
+ pylab.imshow(rec)
+ pylab.show()
+ self.assertTrue(abs(val-1.0) < TOL)
+
+
+__combinations = {
+ 'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+ 'fanflat': [ 'line_fanflat', 'strip_fanflat', 'cuda' ],
+# 'fanflat': [ 'cuda' ],
+ }
+
+__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
+}
+
+for k, l in __combinations.items():
+ for v in l:
+ A = __algs if v != 'cuda' else __algs_CUDA
+ for a, i in A.items():
+ def f(k, v, a, i):
+ return lambda self: self.single_test(k, v, a, i)
+ setattr(Test2DRecScale, 'test_' + a + '_' + k + '_' + v, f(k,v,a,i))
+
+if __name__ == '__main__':
+ unittest.main()
+
--
cgit v1.2.3
From f944d9a95d3cd7c70a137b23147368dc15039c7f Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 4 Apr 2019 15:58:26 +0200
Subject: Add 3D reconstruction scaling test
---
tests/python/test_rec_scaling.py | 96 ++++++++++++++++++++++++++++++++--------
1 file changed, 78 insertions(+), 18 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py
index 33d09f9..d656edc 100644
--- a/tests/python/test_rec_scaling.py
+++ b/tests/python/test_rec_scaling.py
@@ -6,9 +6,19 @@ import pylab
DISPLAY=False
-def VolumeGeometries():
- for s in [0.8, 1.0, 1.25]:
- yield astra.create_vol_geom(128, 128, -64*s, 64*s, -64*s, 64*s)
+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':
@@ -19,17 +29,40 @@ def ProjectionGeometries(type):
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 == '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)
-class Test2DRecScale(unittest.TestCase):
+class TestRecScale(unittest.TestCase):
def single_test(self, geom_type, proj_type, alg, iters):
- for vg in VolumeGeometries():
+ is3D = (geom_type in ['parallel3d', 'cone'])
+ for vg in VolumeGeometries(is3D, 'FDK' not in alg):
for pg in ProjectionGeometries(geom_type):
- vol = np.zeros((128,128))
- vol[50:70,50:70] = 1
+ 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)
- sino_id, sinogram = astra.create_sino(vol, proj_id)
- rec_id = astra.data2d.create('-vol', vg, 0.0 if 'EM' not in alg else 1.0)
+ 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
@@ -37,24 +70,32 @@ class Test2DRecScale(unittest.TestCase):
cfg['ProjectorId'] = proj_id
alg_id = astra.algorithm.create(cfg)
- astra.algorithm.run(alg_id, iters)
- rec = astra.data2d.get(rec_id)
+ 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])
- val = np.sum(rec[55:65,55:65]) / 100.
+ 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()
- pylab.imshow(rec)
+ if not is3D:
+ pylab.imshow(rec)
+ else:
+ pylab.imshow(rec[:,32,:])
pylab.show()
- self.assertTrue(abs(val-1.0) < TOL)
+ #self.assertTrue(abs(val-1.0) < TOL)
__combinations = {
'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
'fanflat': [ 'line_fanflat', 'strip_fanflat', 'cuda' ],
-# 'fanflat': [ 'cuda' ],
+ 'parallel3d': [ 'cuda3d' ],
+ 'cone': [ 'cuda3d' ],
}
__algs = {
@@ -66,14 +107,33 @@ __algs_CUDA = {
'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:
- A = __algs if v != 'cuda' else __algs_CUDA
+ 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(Test2DRecScale, 'test_' + a + '_' + k + '_' + v, f(k,v,a,i))
-
+ setattr(TestRecScale, 'test_' + a + '_' + k + '_' + v, f(k,v,a,i))
+
if __name__ == '__main__':
unittest.main()
--
cgit v1.2.3
From b595d260193b39981834211682ff41fd1a91e398 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 4 Apr 2019 15:59:02 +0200
Subject: Enable all 2D projector tests
---
tests/python/test_line2d.py | 6 +++---
tests/python/test_rec_scaling.py | 7 +++++--
2 files changed, 8 insertions(+), 5 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 8b6e200..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
diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py
index d656edc..1bd3267 100644
--- a/tests/python/test_rec_scaling.py
+++ b/tests/python/test_rec_scaling.py
@@ -43,6 +43,8 @@ def ProjectionGeometries(type):
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):
@@ -88,7 +90,7 @@ class TestRecScale(unittest.TestCase):
else:
pylab.imshow(rec[:,32,:])
pylab.show()
- #self.assertTrue(abs(val-1.0) < TOL)
+ self.assertTrue(abs(val-1.0) < TOL)
__combinations = {
@@ -99,7 +101,8 @@ __combinations = {
}
__algs = {
- 'SIRT': 50, 'SART': 10*180, 'CGLS': 30, 'FBP': 1
+ 'SIRT': 50, 'SART': 10*180, 'CGLS': 30,
+ 'FBP': 1
}
__algs_CUDA = {
--
cgit v1.2.3
From 0c2482e1dce65ded6215cd5634d86b4c00381e27 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 5 Apr 2019 16:07:48 +0200
Subject: Add unit tests for 3D adjoints
---
tests/python/test_rec_scaling.py | 81 +++++++++++++++++++++++++++++++++++++---
1 file changed, 76 insertions(+), 5 deletions(-)
(limited to 'tests/python')
diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py
index 1bd3267..621fd8a 100644
--- a/tests/python/test_rec_scaling.py
+++ b/tests/python/test_rec_scaling.py
@@ -33,12 +33,46 @@ def ProjectionGeometries(type):
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):
@@ -92,13 +126,44 @@ class TestRecScale(unittest.TestCase):
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' ],
- }
+ '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,
@@ -137,6 +202,12 @@ for k, l in __combinations.items():
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()
--
cgit v1.2.3