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(-) 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