From 0f4ceb4c7f3f63fddf8fbf44c59fcd8f415e3847 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn 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') 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 Date: Fri, 29 Mar 2019 15:03:57 +0100 Subject: Adjust linear/cuda kernels to line integral scaling --- cuda/2d/astra.cu | 3 ++- cuda/2d/par_fp.cu | 10 ++++------ include/astra/ParallelBeamLinearKernelProjector2D.inl | 6 ++---- tests/python/test_line2d.py | 16 ++++------------ 4 files changed, 12 insertions(+), 23 deletions(-) (limited to 'tests') diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index ec03517..7ff1c95 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -302,7 +302,8 @@ static bool convertAstraGeometry_internal(const CVolumeGeometry2D* pVolGeom, pProjs[i].scale(factor); } // CHECKME: Check factor - fOutputScale *= pVolGeom->getPixelLengthX() * pVolGeom->getPixelLengthY(); + // NB: Only valid for square pixels + fOutputScale *= pVolGeom->getPixelLengthX(); return true; } diff --git a/cuda/2d/par_fp.cu b/cuda/2d/par_fp.cu index 0835301..e03381c 100644 --- a/cuda/2d/par_fp.cu +++ b/cuda/2d/par_fp.cu @@ -115,10 +115,9 @@ __global__ void FPhorizontal_simple(float* D_projData, unsigned int projPitch, u float fSliceStep = cos_theta / sin_theta; float fDistCorr; if (sin_theta > 0.0f) - fDistCorr = -fDetStep; + fDistCorr = outputScale / sin_theta; else - fDistCorr = fDetStep; - fDistCorr *= outputScale; + fDistCorr = -outputScale / sin_theta; float fVal = 0.0f; // project detector on slice @@ -193,10 +192,9 @@ __global__ void FPvertical_simple(float* D_projData, unsigned int projPitch, uns float fSliceStep = sin_theta / cos_theta; float fDistCorr; if (cos_theta < 0.0f) - fDistCorr = -fDetStep; + fDistCorr = -outputScale / cos_theta; else - fDistCorr = fDetStep; - fDistCorr *= outputScale; + fDistCorr = outputScale / cos_theta; float fVal = 0.0f; float fP = (detector - 0.5f*dims.iProjDets + 0.5f - gC_angle_offset[angle]) * fDetStep + (startSlice - 0.5f*dims.iVolHeight + 0.5f) * fSliceStep + 0.5f*dims.iVolWidth - 0.5f + 0.5f; diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl index 10d4892..1acd422 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.inl +++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl @@ -154,16 +154,14 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; - float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY); - const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); if (vertical) { RxOverRy = proj->fRayX/proj->fRayY; - lengthPerRow = detSize * m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); + lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; } else { RyOverRx = proj->fRayY/proj->fRayX; - lengthPerCol = detSize * m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); + lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; } 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 Date: Fri, 29 Mar 2019 16:01:02 +0100 Subject: Adjust distance driven kernels to line integral scaling --- include/astra/ParallelBeamDistanceDrivenProjector2D.inl | 9 ++++++--- tests/python/test_line2d.py | 15 +++++++++++---- 2 files changed, 17 insertions(+), 7 deletions(-) (limited to 'tests') diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl index aedcee9..3a18c6f 100644 --- a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl +++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl @@ -72,7 +72,7 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro const int rowCount = m_pVolumeGeometry->getGridRowCount(); // Performance note: - // This is not a very well optimizated version of the distance driven + // This is not a very well optimized version of the distance driven // projector. The CPU projector model in ASTRA requires ray-driven iteration, // which limits re-use of intermediate computations. @@ -86,6 +86,9 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; const float32 Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f; + const float32 rayWidth = fabs(proj->fDetUX * proj->fRayY - proj->fDetUY * proj->fRayX) / + sqrt(proj->fRayX * proj->fRayX + proj->fRayY * proj->fRayY); + // loop detectors for (int iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { @@ -100,7 +103,7 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro if (vertical) { const float32 RxOverRy = proj->fRayX/proj->fRayY; - const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY() / rayWidth; const float32 deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; const float32 deltad = 0.5f * fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX); @@ -157,7 +160,7 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro } else { const float32 RyOverRx = proj->fRayY/proj->fRayX; - const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY() / rayWidth; const float32 deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; const float32 deltad = 0.5f * fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY); 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 Date: Fri, 29 Mar 2019 19:08:11 +0100 Subject: Adjust strip kernels to line integral scaling --- .../astra/FanFlatBeamStripKernelProjector2D.inl | 12 +++++++---- .../astra/ParallelBeamStripKernelProjector2D.inl | 8 ++++++-- tests/python/test_line2d.py | 24 +++++++++++----------- 3 files changed, 26 insertions(+), 18 deletions(-) (limited to 'tests') diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl index 889d0a3..f5a688c 100644 --- a/include/astra/FanFlatBeamStripKernelProjector2D.inl +++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl @@ -109,8 +109,12 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i // POLICY: RAY PRIOR if (!p.rayPrior(iRayIndex)) continue; - float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + - (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW; + float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW; + dist_srcDetPixSquared = dist_srcDetPixSquared * dist_srcDetPixSquared / (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() * DW * DW); + + + + //float32 InvRayWidthSquared = (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance()) / dist_srcDetPixSquared; float32 sin_theta_left, cos_theta_left; float32 sin_theta_right, cos_theta_right; @@ -257,8 +261,8 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i // POLICY: RAY PRIOR if (!p.rayPrior(iRayIndex)) continue; - float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + - (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW; + float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW; + dist_srcDetPixSquared = dist_srcDetPixSquared * dist_srcDetPixSquared / (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() * DW * DW); // get theta_l = alpha_left + theta and theta_r = alpha_right + theta float32 sin_theta_left, cos_theta_left; diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl index 0d775b3..bed02ab 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.inl +++ b/include/astra/ParallelBeamStripKernelProjector2D.inl @@ -142,6 +142,10 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; + const float32 rayWidth = fabs(proj->fDetUX * proj->fRayY - proj->fDetUY * proj->fRayX) / + sqrt(proj->fRayX * proj->fRayX + proj->fRayY * proj->fRayY); + const float32 relPixelArea = pixelArea / rayWidth; + bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); if (vertical) { RxOverRy = proj->fRayX/proj->fRayY; @@ -219,7 +223,7 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, else if (-S < offsetL) res -= 0.5f + offsetL; else if (-T < offsetL) res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS; - p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res); + p.addWeight(iRayIndex, iVolumeIndex, relPixelArea*res); p.pixelPosterior(iVolumeIndex); } } @@ -272,7 +276,7 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, else if (-S < offsetL) res -= 0.5f + offsetL; else if (-T < offsetL) res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS; - p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res); + p.addWeight(iRayIndex, iVolumeIndex, relPixelArea*res); p.pixelPosterior(iVolumeIndex); } } 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 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') 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 Date: Fri, 29 Mar 2019 21:21:29 +0100 Subject: Adjust adjoint to line integral scaling --- cuda/2d/algo.cu | 12 +++---- cuda/2d/fan_bp.cu | 62 +++++++++++++++++++++++++---------- cuda/2d/par_bp.cu | 20 +++++++++-- cuda/3d/cone_bp.cu | 6 ++-- src/CudaReconstructionAlgorithm2D.cpp | 3 +- tests/python/test_line2d.py | 59 +++++++++++++++++++++++++++++++++ 6 files changed, 132 insertions(+), 30 deletions(-) (limited to 'tests') diff --git a/cuda/2d/algo.cu b/cuda/2d/algo.cu index b4c2864..11422ff 100644 --- a/cuda/2d/algo.cu +++ b/cuda/2d/algo.cu @@ -258,10 +258,10 @@ bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPit if (!ok) return false; - // rescale sinogram to adjust for pixel size - processSino(D_sinoData, fSinogramScale, - //1.0f/(fPixelSize*fPixelSize), - sinoPitch, dims); + // rescale sinogram + if (fSinogramScale != 1.0f) + processSino(D_sinoData, fSinogramScale, + sinoPitch, dims); ok = copyVolumeToDevice(pfReconstruction, iReconstructionPitch, dims, @@ -331,11 +331,11 @@ bool ReconAlgo::callBP(float* D_volumeData, unsigned int volumePitch, if (parProjs) { assert(!fanProjs); return BP(D_volumeData, volumePitch, D_projData, projPitch, - dims, parProjs, outputScale); + dims, parProjs, fOutputScale * outputScale); } else { assert(fanProjs); return FanBP(D_volumeData, volumePitch, D_projData, projPitch, - dims, fanProjs, outputScale); + dims, fanProjs, fOutputScale * outputScale); } } diff --git a/cuda/2d/fan_bp.cu b/cuda/2d/fan_bp.cu index dac3ac2..2072cd4 100644 --- a/cuda/2d/fan_bp.cu +++ b/cuda/2d/fan_bp.cu @@ -48,7 +48,7 @@ const unsigned int g_anglesPerBlock = 16; const unsigned int g_blockSliceSize = 32; const unsigned int g_blockSlices = 16; -const unsigned int g_MaxAngles = 2560; +const unsigned int g_MaxAngles = 2240; __constant__ float gC_SrcX[g_MaxAngles]; __constant__ float gC_SrcY[g_MaxAngles]; @@ -56,6 +56,7 @@ __constant__ float gC_DetSX[g_MaxAngles]; __constant__ float gC_DetSY[g_MaxAngles]; __constant__ float gC_DetUX[g_MaxAngles]; __constant__ float gC_DetUY[g_MaxAngles]; +__constant__ float gC_Scale[g_MaxAngles]; static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int width, unsigned int height, cudaTextureAddressMode mode = cudaAddressModeBorder) @@ -96,8 +97,6 @@ __global__ void devFanBP(float* D_volData, unsigned int volPitch, unsigned int s float fVal = 0.0f; float fA = startAngle + 0.5f; - // TODO: Distance correction? - for (int angle = startAngle; angle < endAngle; ++angle) { const float fSrcX = gC_SrcX[angle]; @@ -106,15 +105,24 @@ __global__ void devFanBP(float* D_volData, unsigned int volPitch, unsigned int s const float fDetSY = gC_DetSY[angle]; const float fDetUX = gC_DetUX[angle]; const float fDetUY = gC_DetUY[angle]; + const float fScale = gC_Scale[angle]; const float fXD = fSrcX - fX; const float fYD = fSrcY - fY; const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX; const float fDen = fDetUX * fYD - fDetUY * fXD; - - const float fT = fNum / fDen; - fVal += tex2D(gT_FanProjTexture, fT, fA); + + // fDen = || u (x-s) || (2x2 determinant) + + // Scale factor is the approximate number of rays traversing this pixel, + // given by the inverse size of a detector pixel scaled by the magnification + // factor of this pixel. + // Magnification factor is || u (d-s) || / || u (x-s) || + + const float fr = __fdividef(1.0f, fDen); + const float fT = fNum * fr; + fVal += tex2D(gT_FanProjTexture, fT, fA) * fScale * fr; fA += 1.0f; } @@ -148,8 +156,6 @@ __global__ void devFanBP_SS(float* D_volData, unsigned int volPitch, unsigned in float fVal = 0.0f; float fA = startAngle + 0.5f; - // TODO: Distance correction? - for (int angle = startAngle; angle < endAngle; ++angle) { const float fSrcX = gC_SrcX[angle]; @@ -158,6 +164,7 @@ __global__ void devFanBP_SS(float* D_volData, unsigned int volPitch, unsigned in const float fDetSY = gC_DetSY[angle]; const float fDetUX = gC_DetUX[angle]; const float fDetUY = gC_DetUY[angle]; + const float fScale = gC_Scale[angle]; // TODO: Optimize these loops... float fX = fXb; @@ -169,9 +176,10 @@ __global__ void devFanBP_SS(float* D_volData, unsigned int volPitch, unsigned in const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX; const float fDen = fDetUX * fYD - fDetUY * fXD; - - const float fT = fNum / fDen; - fVal += tex2D(gT_FanProjTexture, fT, fA); + const float fr = __fdividef(1.0f, fDen); + + const float fT = fNum * fr; + fVal += tex2D(gT_FanProjTexture, fT, fA) * fScale * fr; fY -= fSubStep; } fX += fSubStep; @@ -202,8 +210,6 @@ __global__ void devFanBP_SART(float* D_volData, unsigned int volPitch, const SDi float* volData = (float*)D_volData; - // TODO: Distance correction? - // TODO: Constant memory vs parameters. const float fSrcX = gC_SrcX[0]; const float fSrcY = gC_SrcY[0]; @@ -211,15 +217,17 @@ __global__ void devFanBP_SART(float* D_volData, unsigned int volPitch, const SDi const float fDetSY = gC_DetSY[0]; const float fDetUX = gC_DetUX[0]; const float fDetUY = gC_DetUY[0]; + const float fScale = gC_Scale[0]; const float fXD = fSrcX - fX; const float fYD = fSrcY - fY; const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX; const float fDen = fDetUX * fYD - fDetUY * fXD; - - const float fT = fNum / fDen; - const float fVal = tex2D(gT_FanProjTexture, fT, 0.5f); + const float fr = __fdividef(1.0f, fDen); + + const float fT = fNum * fr; + const float fVal = tex2D(gT_FanProjTexture, fT, 0.5f) * fScale * fr; volData[Y*volPitch+X] += fVal * fOutputScale; } @@ -248,7 +256,7 @@ __global__ void devFanBP_FBPWeighted(float* D_volData, unsigned int volPitch, un float fVal = 0.0f; float fA = startAngle + 0.5f; - // TODO: Distance correction? + // TODO: Update for new projection scaling for (int angle = startAngle; angle < endAngle; ++angle) { @@ -299,6 +307,13 @@ bool FanBP_internal(float* D_volumeData, unsigned int volumePitch, #undef TRANSFER_TO_CONSTANT + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + double detsize = sqrt(angles[i].fDetUX * angles[i].fDetUX + angles[i].fDetUY * angles[i].fDetUY); + double scale = (angles[i].fDetUX * (angles[i].fSrcY - angles[i].fDetSY) - angles[i].fDetUY * (angles[i].fSrcX - angles[i].fDetSX)) / detsize; + tmp[i] = (float)scale; + } + cudaMemcpyToSymbol(gC_Scale, tmp, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); + delete[] tmp; dim3 dimBlock(g_blockSlices, g_blockSliceSize); @@ -346,6 +361,14 @@ bool FanBP_FBPWeighted_internal(float* D_volumeData, unsigned int volumePitch, #undef TRANSFER_TO_CONSTANT + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + double detsize = sqrt(angles[i].fDetUX * angles[i].fDetUX + angles[i].fDetUY * angles[i].fDetUY); + double scale = (angles[i].fDetUX * (angles[i].fSrcY - angles[i].fDetSY) - angles[i].fDetUY * (angles[i].fSrcX - angles[i].fDetSX)) / detsize; + tmp[i] = (float)scale; + } + cudaMemcpyToSymbol(gC_Scale, tmp, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); + + delete[] tmp; dim3 dimBlock(g_blockSlices, g_blockSliceSize); @@ -389,6 +412,11 @@ bool FanBP_SART(float* D_volumeData, unsigned int volumePitch, #undef TRANSFER_TO_CONSTANT + double detsize = sqrt(angles[angle].fDetUX * angles[angle].fDetUX + angles[angle].fDetUY * angles[angle].fDetUY); + double scale = (angles[angle].fDetUX * (angles[angle].fSrcY - angles[angle].fDetSY) - angles[angle].fDetUY * (angles[angle].fSrcX - angles[angle].fDetSX)) / detsize; + float tmp = (float)scale; + cudaMemcpyToSymbol(gC_Scale, &tmp, sizeof(float), 0, cudaMemcpyHostToDevice); + dim3 dimBlock(g_blockSlices, g_blockSliceSize); dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices, (dims.iVolHeight+g_blockSliceSize-1)/g_blockSliceSize); diff --git a/cuda/2d/par_bp.cu b/cuda/2d/par_bp.cu index 09a6554..21484b1 100644 --- a/cuda/2d/par_bp.cu +++ b/cuda/2d/par_bp.cu @@ -53,6 +53,7 @@ const unsigned int g_MaxAngles = 2560; __constant__ float gC_angle_scaled_sin[g_MaxAngles]; __constant__ float gC_angle_scaled_cos[g_MaxAngles]; __constant__ float gC_angle_offset[g_MaxAngles]; +__constant__ float gC_angle_scale[g_MaxAngles]; static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int width, unsigned int height, cudaTextureAddressMode mode = cudaAddressModeBorder) { @@ -70,6 +71,7 @@ static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int wi return true; } +// TODO: Templated version with/without scale? (Or only the global outputscale) __global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int startAngle, const SDimensions dims, float fOutputScale) { const int relX = threadIdx.x; @@ -97,9 +99,10 @@ __global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int star const float scaled_cos_theta = gC_angle_scaled_cos[angle]; const float scaled_sin_theta = gC_angle_scaled_sin[angle]; const float TOffset = gC_angle_offset[angle]; + const float scale = gC_angle_scale[angle]; const float fT = fX * scaled_cos_theta - fY * scaled_sin_theta + TOffset; - fVal += tex2D(gT_projTexture, fT, fA); + fVal += tex2D(gT_projTexture, fT, fA) * scale; fA += 1.0f; } @@ -138,6 +141,7 @@ __global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int s const float cos_theta = gC_angle_scaled_cos[angle]; const float sin_theta = gC_angle_scaled_sin[angle]; const float TOffset = gC_angle_offset[angle]; + const float scale = gC_angle_scale[angle]; float fT = fX * cos_theta - fY * sin_theta + TOffset; @@ -145,7 +149,7 @@ __global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int s float fTy = fT; fT += fSubStep * cos_theta; for (int iSubY = 0; iSubY < dims.iRaysPerPixelDim; ++iSubY) { - fVal += tex2D(gT_projTexture, fTy, fA); + fVal += tex2D(gT_projTexture, fTy, fA) * scale; fTy -= fSubStep * sin_theta; } } @@ -172,6 +176,8 @@ __global__ void devBP_SART(float* D_volData, unsigned int volPitch, float offset const float fT = fX * angle_cos - fY * angle_sin + offset; const float fVal = tex2D(gT_projTexture, fT, 0.5f); + // NB: The 'scale' constant in devBP has been folded into fOutputScale by the caller here + D_volData[Y*volPitch+X] += fVal * fOutputScale; } @@ -186,27 +192,34 @@ bool BP_internal(float* D_volumeData, unsigned int volumePitch, float* angle_scaled_sin = new float[dims.iProjAngles]; float* angle_scaled_cos = new float[dims.iProjAngles]; float* angle_offset = new float[dims.iProjAngles]; + float* angle_scale = new float[dims.iProjAngles]; bindProjDataTexture(D_projData, projPitch, dims.iProjDets, dims.iProjAngles); for (unsigned int i = 0; i < dims.iProjAngles; ++i) { double d = angles[i].fDetUX * angles[i].fRayY - angles[i].fDetUY * angles[i].fRayX; angle_scaled_cos[i] = angles[i].fRayY / d; - angle_scaled_sin[i] = -angles[i].fRayX / d; // TODO: Check signs + angle_scaled_sin[i] = -angles[i].fRayX / d; angle_offset[i] = (angles[i].fDetSY * angles[i].fRayX - angles[i].fDetSX * angles[i].fRayY) / d; + angle_scale[i] = sqrt(angles[i].fRayX * angles[i].fRayX + angles[i].fRayY * angles[i].fRayY) / abs(d); } + //fprintf(stderr, "outputscale in BP_internal: %f, %f\n", fOutputScale, angle_scale[0]); + //fprintf(stderr, "ray in BP_internal: %f,%f (length %f)\n", angles[0].fRayX, angles[0].fRayY, sqrt(angles[0].fRayX * angles[0].fRayX + angles[0].fRayY * angles[0].fRayY)); cudaError_t e1 = cudaMemcpyToSymbol(gC_angle_scaled_sin, angle_scaled_sin, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); cudaError_t e2 = cudaMemcpyToSymbol(gC_angle_scaled_cos, angle_scaled_cos, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); cudaError_t e3 = cudaMemcpyToSymbol(gC_angle_offset, angle_offset, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); + cudaError_t e4 = cudaMemcpyToSymbol(gC_angle_scale, angle_scale, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); assert(e1 == cudaSuccess); assert(e2 == cudaSuccess); assert(e3 == cudaSuccess); + assert(e4 == cudaSuccess); delete[] angle_scaled_sin; delete[] angle_scaled_cos; delete[] angle_offset; + delete[] angle_scale; dim3 dimBlock(g_blockSlices, g_blockSliceSize); dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices, @@ -267,6 +280,7 @@ bool BP_SART(float* D_volumeData, unsigned int volumePitch, float angle_scaled_cos = angles[angle].fRayY / d; float angle_scaled_sin = -angles[angle].fRayX / d; // TODO: Check signs float angle_offset = (angles[angle].fDetSY * angles[angle].fRayX - angles[angle].fDetSX * angles[angle].fRayY) / d; + fOutputScale *= sqrt(angles[angle].fRayX * angles[angle].fRayX + angles[angle].fRayY * angles[angle].fRayY) / abs(d); dim3 dimBlock(g_blockSlices, g_blockSliceSize); dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices, diff --git a/cuda/3d/cone_bp.cu b/cuda/3d/cone_bp.cu index feebda2..2a7ec80 100644 --- a/cuda/3d/cone_bp.cu +++ b/cuda/3d/cone_bp.cu @@ -140,8 +140,10 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng if (FDKWEIGHT) { // The correct factor here is this one: // Z[idx] += (fr*fCd.w)*(fr*fCd.w)*fVal; - // This is the square of the inverse magnification factor - // from fX,fY,fZ to the detector. + // This is the square of the magnification factor + // from fX,fY,fZ to the virtual detector, where the + // virtual detector is the plane through the origin + // parallel to the detector, so spanned by u,v. // Since we are assuming we have a circular cone // beam trajectory, fCd.w is constant, and we instead diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 1e81390..939a026 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -309,8 +309,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) m_bAlgoInit = true; } - float fPixelSize = volgeom.getPixelLengthX(); - float fSinogramScale = 1.0f/(fPixelSize*fPixelSize); + float fSinogramScale = 1.0f; ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(), fSinogramScale, m_pReconstruction->getDataConst(), volgeom.getGridColCount(), 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 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') 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 5ea1bf556419204511195fa5b2bedbd1318b51ff Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 3 Apr 2019 23:17:06 +0200 Subject: Remove C++ projector tests These have been superseded by python versions. --- build/linux/Makefile.in | 2 - tests/test_ParallelBeamLineKernelProjector2D.cpp | 82 ---------- tests/test_ParallelBeamLinearKernelProjector2D.cpp | 170 --------------------- 3 files changed, 254 deletions(-) delete mode 100644 tests/test_ParallelBeamLineKernelProjector2D.cpp delete mode 100644 tests/test_ParallelBeamLinearKernelProjector2D.cpp (limited to 'tests') diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 209206e..f478bb7 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -257,8 +257,6 @@ endif TEST_OBJECTS=\ tests/main.o \ tests/test_AstraObjectManager.o \ - tests/test_ParallelBeamLineKernelProjector2D.o \ - tests/test_ParallelBeamLinearKernelProjector2D.o \ tests/test_Float32Data2D.o \ tests/test_VolumeGeometry2D.o \ tests/test_ParallelProjectionGeometry2D.o \ diff --git a/tests/test_ParallelBeamLineKernelProjector2D.cpp b/tests/test_ParallelBeamLineKernelProjector2D.cpp deleted file mode 100644 index 71130c1..0000000 --- a/tests/test_ParallelBeamLineKernelProjector2D.cpp +++ /dev/null @@ -1,82 +0,0 @@ -/* ------------------------------------------------------------------------ -Copyright: 2010-2018, imec Vision Lab, University of Antwerp - 2014-2018, CWI, Amsterdam - -Contact: astra@astra-toolbox.com -Website: http://www.astra-toolbox.com/ - -This file is part of the ASTRA Toolbox. - - -The ASTRA Toolbox is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -The ASTRA Toolbox is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with the ASTRA Toolbox. If not, see . - ------------------------------------------------------------------------ -*/ - - -#define BOOST_TEST_DYN_LINK -#include -#include -#include - -#include "astra/ParallelBeamLineKernelProjector2D.h" -#include "astra/ParallelProjectionGeometry2D.h" -#include "astra/VolumeGeometry2D.h" - -struct TestParallelBeamLineKernelProjector2D { - TestParallelBeamLineKernelProjector2D() - { - astra::float32 angles[] = { 1.0f }; - BOOST_REQUIRE( projGeom.initialize(1, 9, 1.0f, angles) ); - BOOST_REQUIRE( volGeom.initialize(6, 4) ); - BOOST_REQUIRE( proj.initialize(&projGeom, &volGeom) ); - } - ~TestParallelBeamLineKernelProjector2D() - { - - } - - astra::CParallelBeamLineKernelProjector2D proj; - astra::CParallelProjectionGeometry2D projGeom; - astra::CVolumeGeometry2D volGeom; -}; - -BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_General, TestParallelBeamLineKernelProjector2D ) -{ - -} - -BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_Rectangle, TestParallelBeamLineKernelProjector2D ) -{ - int iMax = proj.getProjectionWeightsCount(0); - BOOST_REQUIRE(iMax > 0); - - astra::SPixelWeight* pPix = new astra::SPixelWeight[iMax]; - BOOST_REQUIRE(pPix); - - int iCount; - proj.computeSingleRayWeights(0, 4, pPix, iMax, iCount); - BOOST_REQUIRE(iCount <= iMax); - - astra::float32 fWeight = 0; - for (int i = 0; i < iCount; ++i) - fWeight += pPix[i].m_fWeight; - - BOOST_CHECK_SMALL(fWeight - 7.13037f, 0.00001f); // 6 / sin(1) - - delete[] pPix; -} - - diff --git a/tests/test_ParallelBeamLinearKernelProjector2D.cpp b/tests/test_ParallelBeamLinearKernelProjector2D.cpp deleted file mode 100644 index 4610aa5..0000000 --- a/tests/test_ParallelBeamLinearKernelProjector2D.cpp +++ /dev/null @@ -1,170 +0,0 @@ -/* ------------------------------------------------------------------------ -Copyright: 2010-2018, imec Vision Lab, University of Antwerp - 2014-2018, CWI, Amsterdam - -Contact: astra@astra-toolbox.com -Website: http://www.astra-toolbox.com/ - -This file is part of the ASTRA Toolbox. - - -The ASTRA Toolbox is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -The ASTRA Toolbox is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with the ASTRA Toolbox. If not, see . - ------------------------------------------------------------------------ -*/ - - -#define BOOST_TEST_DYN_LINK -#include -#include -#include - -#include "astra/ParallelBeamLineKernelProjector2D.h" -#include "astra/ParallelBeamLinearKernelProjector2D.h" -#include "astra/ParallelBeamStripKernelProjector2D.h" -#include "astra/ParallelProjectionGeometry2D.h" -#include "astra/VolumeGeometry2D.h" -#include "astra/ProjectionGeometry2D.h" - -#include - -using astra::float32; - -struct TestParallelBeamLinearKernelProjector2D { - TestParallelBeamLinearKernelProjector2D() - { - astra::float32 angles[] = { 2.6f }; - BOOST_REQUIRE( projGeom.initialize(1, 3, 1.0f, angles) ); - BOOST_REQUIRE( volGeom.initialize(3, 2) ); - BOOST_REQUIRE( proj.initialize(&projGeom, &volGeom) ); - } - ~TestParallelBeamLinearKernelProjector2D() - { - - } - - astra::CParallelBeamLinearKernelProjector2D proj; - astra::CParallelProjectionGeometry2D projGeom; - astra::CVolumeGeometry2D volGeom; -}; - -BOOST_FIXTURE_TEST_CASE( testParallelBeamLinearKernelProjector2D_General, TestParallelBeamLinearKernelProjector2D ) -{ - -} - - -// Compute linear kernel for a single volume pixel/detector pixel combination -float32 compute_linear_kernel(const astra::CProjectionGeometry2D& projgeom, const astra::CVolumeGeometry2D& volgeom, - int iX, int iY, int iDet, float32 fAngle) -{ - // projection of center of volume pixel on detector array - float32 fDetProj = (iX - (volgeom.getGridColCount()-1.0f)/2.0f ) * volgeom.getPixelLengthX() * cos(fAngle) - (iY - (volgeom.getGridRowCount()-1.0f)/2.0f ) * volgeom.getPixelLengthY() * sin(fAngle); - - // start of detector pixel on detector array - float32 fDetStart = projgeom.indexToDetectorOffset(iDet) - 0.5f; - -// printf("(%d,%d,%d): %f in (%f,%f)\n", iX,iY,iDet,fDetProj, fDetStart, fDetStart+1.0f); - - // projection of center of next volume pixel on detector array - float32 fDetStep; - // length of projection ray through volume pixel - float32 fWeight; - - if (fabs(cos(fAngle)) > fabs(sin(fAngle))) { - fDetStep = volgeom.getPixelLengthY() * fabs(cos(fAngle)); - fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthX() * 1.0f / fabs(cos(fAngle)); - } else { - fDetStep = volgeom.getPixelLengthX() * fabs(sin(fAngle)); - fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthY() * 1.0f / fabs(sin(fAngle)); - } - -// printf("step: %f\n weight: %f\n", fDetStep, fWeight); - - // center of detector pixel on detector array - float32 fDetCenter = fDetStart + 0.5f; - - // unweighted contribution of this volume pixel: - // linear interpolation between - // fDetCenter - fDetStep |---> 0 - // fDetCenter |---> 1 - // fDetCenter + fDetStep |---> 0 - float32 fBase; - if (fDetCenter <= fDetProj) { - fBase = (fDetCenter - (fDetProj - fDetStep))/fDetStep; - } else { - fBase = ((fDetProj + fDetStep) - fDetCenter)/fDetStep; - } -// printf("base: %f\n", fBase); - if (fBase < 0) fBase = 0; - return fBase * fWeight; -} - -BOOST_AUTO_TEST_CASE( testParallelBeamLinearKernelProjector2D_Rectangles ) -{ - astra::CParallelBeamLinearKernelProjector2D proj; - astra::CParallelProjectionGeometry2D projGeom; - astra::CVolumeGeometry2D volGeom; - - const unsigned int iRandomTestCount = 100; - - unsigned int iSeed = time(0); - srand(iSeed); - - for (unsigned int iTest = 0; iTest < iRandomTestCount; ++iTest) { - int iDetectorCount = 1 + (rand() % 100); - int iRows = 1 + (rand() % 100); - int iCols = 1 + (rand() % 100); - - - astra::float32 angles[] = { rand() * 2.0f*astra::PI / RAND_MAX }; - projGeom.initialize(1, iDetectorCount, 0.8f, angles); - volGeom.initialize(iCols, iRows); - proj.initialize(&projGeom, &volGeom); - - int iMax = proj.getProjectionWeightsCount(0); - BOOST_REQUIRE(iMax > 0); - - astra::SPixelWeight* pPix = new astra::SPixelWeight[iMax]; - BOOST_REQUIRE(pPix); - - astra::float32 fWeight = 0; - for (int iDet = 0; iDet < projGeom.getDetectorCount(); ++iDet) { - int iCount; - proj.computeSingleRayWeights(0, iDet, pPix, iMax, iCount); - BOOST_REQUIRE(iCount <= iMax); - - astra::float32 fW = 0; - for (int i = 0; i < iCount; ++i) { - float32 fTest = compute_linear_kernel( - projGeom, - volGeom, - pPix[i].m_iIndex % volGeom.getGridColCount(), - pPix[i].m_iIndex / volGeom.getGridColCount(), - iDet, - projGeom.getProjectionAngle(0)); - BOOST_CHECK_SMALL( pPix[i].m_fWeight - fTest, 0.0004f); - fW += pPix[i].m_fWeight; - } - - fWeight += fW; - - } - - delete[] pPix; - } -} - - -- cgit v1.2.3 From f944d9a95d3cd7c70a137b23147368dc15039c7f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn 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') 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 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') 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 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') 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