summaryrefslogtreecommitdiffstats
path: root/include/astra
diff options
context:
space:
mode:
Diffstat (limited to 'include/astra')
-rw-r--r--include/astra/CudaProjector3D.h2
-rw-r--r--include/astra/FanFlatBeamLineKernelProjector2D.inl6
-rw-r--r--include/astra/FanFlatBeamStripKernelProjector2D.inl12
-rw-r--r--include/astra/Features.h16
-rw-r--r--include/astra/ParallelBeamDistanceDrivenProjector2D.inl9
-rw-r--r--include/astra/ParallelBeamLineKernelProjector2D.inl31
-rw-r--r--include/astra/ParallelBeamLinearKernelProjector2D.inl19
-rw-r--r--include/astra/ParallelBeamStripKernelProjector2D.inl33
-rw-r--r--include/astra/cuda/2d/algo.h9
-rw-r--r--include/astra/cuda/2d/cgls.h2
-rw-r--r--include/astra/cuda/2d/fbp.h6
-rw-r--r--include/astra/cuda/3d/fdk.h2
-rw-r--r--include/astra/cuda/3d/mem3d.h2
-rw-r--r--include/astra/cuda/3d/util3d.h47
14 files changed, 132 insertions, 64 deletions
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;
};
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/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/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
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/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl
index a9f1aa0..903ebb6 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLineKernelProjector2D.inl
@@ -166,24 +166,7 @@ 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);
- 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);
- 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;
@@ -204,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;
@@ -248,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 10d4892..53451e5 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -154,18 +154,7 @@ 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);
- 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);
- deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
- }
Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
@@ -186,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;
@@ -209,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 0d775b3..2031560 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -142,20 +142,11 @@ 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;
- 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;
@@ -176,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;
@@ -219,7 +216,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);
}
}
@@ -229,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;
@@ -272,7 +275,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/include/astra/cuda/2d/algo.h b/include/astra/cuda/2d/algo.h
index 3fccdef..b670b8b 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);
@@ -133,7 +136,7 @@ protected:
SDimensions dims;
SParProjection* parProjs;
SFanProjection* fanProjs;
- float fOutputScale;
+ float fProjectorScale;
bool freeGPUMemory;
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);
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;
};
}
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);
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);
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-11 05:22:29 +0000