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 --- include/astra/FanFlatBeamLineKernelProjector2D.inl | 6 ++---- include/astra/ParallelBeamLineKernelProjector2D.inl | 6 ++---- 2 files changed, 4 insertions(+), 8 deletions(-) (limited to 'include') diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl index 58dec61..8f2e673 100644 --- a/include/astra/FanFlatBeamLineKernelProjector2D.inl +++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl @@ -82,8 +82,6 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in const SFanProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; - float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY); - // loop detectors for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { @@ -104,7 +102,7 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in // vertically if (vertical) { RxOverRy = Rx/Ry; - lengthPerRow = detSize * pixelLengthX * sqrt(Rx*Rx + Ry*Ry) / abs(Ry); + lengthPerRow = pixelLengthX * sqrt(Rx*Rx + Ry*Ry) / abs(Ry); deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; S = 0.5f - 0.5f*fabs(RxOverRy); T = 0.5f + 0.5f*fabs(RxOverRy); @@ -154,7 +152,7 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in // horizontally else { RyOverRx = Ry/Rx; - lengthPerCol = detSize * pixelLengthY * sqrt(Rx*Rx + Ry*Ry) / abs(Rx); + lengthPerCol = pixelLengthY * sqrt(Rx*Rx + Ry*Ry) / abs(Rx); deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; S = 0.5f - 0.5f*fabs(RyOverRx); T = 0.5f + 0.5f*fabs(RyOverRx); diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl index a9f1aa0..1fc67a8 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.inl +++ b/include/astra/ParallelBeamLineKernelProjector2D.inl @@ -166,19 +166,17 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; - float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY); - bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); if (vertical) { RxOverRy = proj->fRayX/proj->fRayY; - lengthPerRow = detSize * pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); + lengthPerRow = pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; S = 0.5f - 0.5f*fabs(RxOverRy); T = 0.5f + 0.5f*fabs(RxOverRy); invTminSTimesLengthPerRow = lengthPerRow / (T - S); } else { RyOverRx = proj->fRayY/proj->fRayX; - lengthPerCol = detSize * pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); + lengthPerCol = pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; S = 0.5f - 0.5f*fabs(RyOverRx); T = 0.5f + 0.5f*fabs(RyOverRx); -- 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 --- include/astra/ParallelBeamLinearKernelProjector2D.inl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) (limited to 'include') 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; } -- 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 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) (limited to 'include') 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); -- 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 --- include/astra/FanFlatBeamStripKernelProjector2D.inl | 12 ++++++++---- include/astra/ParallelBeamStripKernelProjector2D.inl | 8 ++++++-- 2 files changed, 14 insertions(+), 6 deletions(-) (limited to 'include') 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); } } -- cgit v1.2.3 From 48f4e7b165404a0375db300b9fe59da92edf1cce Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Sat, 30 Mar 2019 20:49:42 +0100 Subject: Adjust FBP to line integral scaling --- include/astra/cuda/2d/algo.h | 7 +++++-- include/astra/cuda/2d/cgls.h | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) (limited to 'include') diff --git a/include/astra/cuda/2d/algo.h b/include/astra/cuda/2d/algo.h index 3fccdef..b6ec231 100644 --- a/include/astra/cuda/2d/algo.h +++ b/include/astra/cuda/2d/algo.h @@ -56,6 +56,10 @@ public: bool setSuperSampling(int raysPerDet, int raysPerPixelDim); + // Scale the final reconstruction. + // May be called at any time after setGeometry and before iterate(). Multiple calls stack. + bool setReconstructionScale(float fScale); + virtual bool enableVolumeMask(); virtual bool enableSinogramMask(); @@ -88,8 +92,7 @@ public: // sinogram, reconstruction, volume mask and sinogram mask in system RAM, // respectively. The corresponding pitch variables give the pitches // of these buffers, measured in floats. - // The sinogram is multiplied by fSinogramScale after uploading it. - virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale, + virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, const float* pfReconstruction, unsigned int iReconstructionPitch, const float* pfVolMask, unsigned int iVolMaskPitch, const float* pfSinoMask, unsigned int iSinoMaskPitch); diff --git a/include/astra/cuda/2d/cgls.h b/include/astra/cuda/2d/cgls.h index 375a425..a854a74 100644 --- a/include/astra/cuda/2d/cgls.h +++ b/include/astra/cuda/2d/cgls.h @@ -47,7 +47,7 @@ public: virtual bool setBuffers(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch); - virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale, + virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, const float* pfReconstruction, unsigned int iReconstructionPitch, const float* pfVolMask, unsigned int iVolMaskPitch, const float* pfSinoMask, unsigned int iSinoMaskPitch); -- cgit v1.2.3 From 11114e33fb504b0b74f3d239e77fe248a027cc23 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Sat, 30 Mar 2019 21:12:01 +0100 Subject: Clean up outputscale naming confusion in cuda::algo --- include/astra/cuda/2d/algo.h | 2 +- include/astra/cuda/2d/fbp.h | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) (limited to 'include') diff --git a/include/astra/cuda/2d/algo.h b/include/astra/cuda/2d/algo.h index b6ec231..b670b8b 100644 --- a/include/astra/cuda/2d/algo.h +++ b/include/astra/cuda/2d/algo.h @@ -136,7 +136,7 @@ protected: SDimensions dims; SParProjection* parProjs; SFanProjection* fanProjs; - float fOutputScale; + float fProjectorScale; bool freeGPUMemory; diff --git a/include/astra/cuda/2d/fbp.h b/include/astra/cuda/2d/fbp.h index 1adf3b1..3aa4741 100644 --- a/include/astra/cuda/2d/fbp.h +++ b/include/astra/cuda/2d/fbp.h @@ -79,6 +79,11 @@ public: bool setShortScan(bool ss) { m_bShortScan = ss; return true; } + // Scale the final reconstruction. + // May be called at any time before iterate(). + bool setReconstructionScale(float fScale); + + virtual bool init(); virtual bool iterate(unsigned int iterations); @@ -90,6 +95,7 @@ protected: void* D_filter; // cufftComplex* bool m_bShortScan; + float fReconstructionScale; }; } -- cgit v1.2.3 From 027dc56ffe325652314b4f9e5d804b3b20fd0db3 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Sat, 30 Mar 2019 22:43:26 +0100 Subject: Work around some warnings --- .../astra/ParallelBeamLineKernelProjector2D.inl | 29 +++++++++++----------- .../astra/ParallelBeamLinearKernelProjector2D.inl | 17 ++++++------- .../astra/ParallelBeamStripKernelProjector2D.inl | 25 +++++++++---------- 3 files changed, 34 insertions(+), 37 deletions(-) (limited to 'include') diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl index 1fc67a8..903ebb6 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.inl +++ b/include/astra/ParallelBeamLineKernelProjector2D.inl @@ -167,21 +167,6 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); - if (vertical) { - RxOverRy = proj->fRayX/proj->fRayY; - lengthPerRow = pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); - deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; - S = 0.5f - 0.5f*fabs(RxOverRy); - T = 0.5f + 0.5f*fabs(RxOverRy); - invTminSTimesLengthPerRow = lengthPerRow / (T - S); - } else { - RyOverRx = proj->fRayY/proj->fRayX; - lengthPerCol = pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); - deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; - S = 0.5f - 0.5f*fabs(RyOverRx); - T = 0.5f + 0.5f*fabs(RyOverRx); - invTminSTimesLengthPerCol = lengthPerCol / (T - S); - } Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f; @@ -202,6 +187,13 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i // vertically if (vertical) { + RxOverRy = proj->fRayX/proj->fRayY; + lengthPerRow = pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); + deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; + S = 0.5f - 0.5f*fabs(RxOverRy); + T = 0.5f + 0.5f*fabs(RxOverRy); + invTminSTimesLengthPerRow = lengthPerRow / (T - S); + // calculate c for row 0 c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX; @@ -246,6 +238,13 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i // horizontally else { + RyOverRx = proj->fRayY/proj->fRayX; + lengthPerCol = pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); + deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; + S = 0.5f - 0.5f*fabs(RyOverRx); + T = 0.5f + 0.5f*fabs(RyOverRx); + invTminSTimesLengthPerCol = lengthPerCol / (T - S); + // calculate r for col 0 r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY; diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl index 1acd422..53451e5 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.inl +++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl @@ -155,15 +155,6 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); - if (vertical) { - RxOverRy = proj->fRayX/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 = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); - deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; - } Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f; @@ -184,6 +175,10 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, // vertically if (vertical) { + RxOverRy = proj->fRayX/proj->fRayY; + lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); + deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; + // calculate c for row 0 c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX; @@ -207,6 +202,10 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, // horizontally else { + RyOverRx = proj->fRayY/proj->fRayX; + lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); + deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; + // calculate r for col 0 r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY; diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl index bed02ab..2031560 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.inl +++ b/include/astra/ParallelBeamStripKernelProjector2D.inl @@ -147,19 +147,6 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, const float32 relPixelArea = pixelArea / rayWidth; bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); - if (vertical) { - RxOverRy = proj->fRayX/proj->fRayY; - deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX; - S = 0.5f - 0.5f*fabs(RxOverRy); - T = 0.5f + 0.5f*fabs(RxOverRy); - invTminS = 1.0f / (T-S); - } else { - RyOverRx = proj->fRayY/proj->fRayX; - deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY; - S = 0.5f - 0.5f*fabs(RyOverRx); - T = 0.5f + 0.5f*fabs(RyOverRx); - invTminS = 1.0f / (T-S); - } Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f; @@ -180,6 +167,12 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, // vertically if (vertical) { + RxOverRy = proj->fRayX/proj->fRayY; + deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX; + S = 0.5f - 0.5f*fabs(RxOverRy); + T = 0.5f + 0.5f*fabs(RxOverRy); + invTminS = 1.0f / (T-S); + // calculate cL and cR for row 0 cL = (DLx + (Ey - DLy)*RxOverRy - Ex) * inv_pixelLengthX; cR = (DRx + (Ey - DRy)*RxOverRy - Ex) * inv_pixelLengthX; @@ -233,6 +226,12 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, // horizontally else { + RyOverRx = proj->fRayY/proj->fRayX; + deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY; + S = 0.5f - 0.5f*fabs(RyOverRx); + T = 0.5f + 0.5f*fabs(RyOverRx); + invTminS = 1.0f / (T-S); + // calculate rL and rR for row 0 rL = -(DLy + (Ex - DLx)*RyOverRx - Ey) * inv_pixelLengthY; rR = -(DRy + (Ex - DRx)*RyOverRx - Ey) * inv_pixelLengthY; -- cgit v1.2.3 From 4b576ee6ad461b162dd64538f54df47b65b37972 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 3 Apr 2019 22:31:55 +0200 Subject: Remove obsolete DensityWeighting option --- include/astra/CudaProjector3D.h | 2 -- 1 file changed, 2 deletions(-) (limited to 'include') diff --git a/include/astra/CudaProjector3D.h b/include/astra/CudaProjector3D.h index 60df7bb..9b4ff1f 100644 --- a/include/astra/CudaProjector3D.h +++ b/include/astra/CudaProjector3D.h @@ -117,7 +117,6 @@ public: int getVoxelSuperSampling() const { return m_iVoxelSuperSampling; } int getDetectorSuperSampling() const { return m_iDetectorSuperSampling; } int getGPUIndex() const { return m_iGPUIndex; } - bool getDensityWeighting() const { return m_bDensityWeighting; } protected: @@ -125,7 +124,6 @@ protected: int m_iVoxelSuperSampling; int m_iDetectorSuperSampling; int m_iGPUIndex; - bool m_bDensityWeighting; }; -- cgit v1.2.3 From 744605ed826196aa6e0e3e3b5e945e50d830ed3a Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 3 Apr 2019 22:32:11 +0200 Subject: Add feature flags for changed scaling behaviour --- include/astra/Features.h | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) (limited to 'include') diff --git a/include/astra/Features.h b/include/astra/Features.h index d88ae71..c7ef98c 100644 --- a/include/astra/Features.h +++ b/include/astra/Features.h @@ -38,10 +38,22 @@ _AstraExport bool hasFeature(const std::string &feature); FEATURES: -cuda: is cuda support compiled in? +cuda + is cuda support compiled in? NB: To check if there is also actually a usable GPU, use cudaAvailable() -mex_link: is there support for the matlab command astra_mex_data3d('link')? +mex_link + is there support for the matlab command astra_mex_data3d('link')? + +projectors_scaled_as_line_integrals + This is set since all 2D and 3D, CPU and GPU projectors scale their outputs + to approximate line integrals. (Previously, some 2D projectors were scaled + as area integrals.) + +fan_cone_BP_density_weighting_by_default + This is set since fan beam and cone beam BP operations perform ray density + weighting by default to more closely approximate the true mathematical adjoint. + The DensityWeighting cuda3d projector option is removed. For future backward-incompatible changes, extra features will be added here -- cgit v1.2.3 From 9bf1863c585506fda6b2fd23de08164bb4fb530d Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 3 Apr 2019 23:10:26 +0200 Subject: Add missing header --- include/astra/cuda/3d/mem3d.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'include') diff --git a/include/astra/cuda/3d/mem3d.h b/include/astra/cuda/3d/mem3d.h index 8c3956e..fff1490 100644 --- a/include/astra/cuda/3d/mem3d.h +++ b/include/astra/cuda/3d/mem3d.h @@ -110,7 +110,7 @@ bool copyIntoArray(MemHandle3D handle, MemHandle3D subdata, const SSubDimensions bool FP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iDetectorSuperSampling, astra::Cuda3DProjectionKernel projKernel); -bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling, bool bFDKWeighting); +bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling); bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, bool bShortScan, const float *pfFilter = 0); -- cgit v1.2.3 From 9950fc9bf91073c0168c8847a8f6a8814690f97c Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 4 Apr 2019 12:09:07 +0200 Subject: Small clean up of factors --- include/astra/cuda/3d/fdk.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'include') diff --git a/include/astra/cuda/3d/fdk.h b/include/astra/cuda/3d/fdk.h index 3b1a9e1..6817154 100644 --- a/include/astra/cuda/3d/fdk.h +++ b/include/astra/cuda/3d/fdk.h @@ -35,7 +35,7 @@ namespace astraCUDA3d { bool FDK_PreWeight(cudaPitchedPtr D_projData, float fSrcOrigin, float fDetOrigin, float fZShift, - float fDetUSize, float fDetVSize, float fVoxSize, + float fDetUSize, float fDetVSize, bool bShortScan, const SDimensions3D& dims, const float* angles); -- cgit v1.2.3 From 2e334ee443a8fd3ce0a6218dd877117a67163533 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 5 Apr 2019 16:07:21 +0200 Subject: Adjust par3d adjoint scaling, and clean up --- include/astra/cuda/3d/util3d.h | 47 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) (limited to 'include') diff --git a/include/astra/cuda/3d/util3d.h b/include/astra/cuda/3d/util3d.h index 0146d2d..200abfc 100644 --- a/include/astra/cuda/3d/util3d.h +++ b/include/astra/cuda/3d/util3d.h @@ -64,6 +64,53 @@ float dotProduct3D(cudaPitchedPtr data, unsigned int x, unsigned int y, unsigned int calcNextPowerOfTwo(int _iValue); +struct Vec3 { + double x; + double y; + double z; + Vec3(double x_, double y_, double z_) : x(x_), y(y_), z(z_) { } + Vec3 operator+(const Vec3 &b) const { + return Vec3(x + b.x, y + b.y, z + b.z); + } + Vec3 operator-(const Vec3 &b) const { + return Vec3(x - b.x, y - b.y, z - b.z); + } + Vec3 operator-() const { + return Vec3(-x, -y, -z); + } + double norm() const { + return sqrt(x*x + y*y + z*z); + } +}; + +static double det3x(const Vec3 &b, const Vec3 &c) { + return (b.y * c.z - b.z * c.y); +} +static double det3y(const Vec3 &b, const Vec3 &c) { + return -(b.x * c.z - b.z * c.x); +} + +static double det3z(const Vec3 &b, const Vec3 &c) { + return (b.x * c.y - b.y * c.x); +} + +static double det3(const Vec3 &a, const Vec3 &b, const Vec3 &c) { + return a.x * det3x(b,c) + a.y * det3y(b,c) + a.z * det3z(b,c); +} + +static Vec3 cross3(const Vec3 &a, const Vec3 &b) { + return Vec3(det3x(a,b), det3y(a,b), det3z(a,b)); +} + +static Vec3 scaled_cross3(const Vec3 &a, const Vec3 &b, const Vec3 &sc) { + Vec3 ret = cross3(a, b); + ret.x *= sc.y * sc.z; + ret.y *= sc.x * sc.z; + ret.z *= sc.x * sc.y; + return ret; +} + + } #endif -- cgit v1.2.3