summaryrefslogtreecommitdiffstats
path: root/include/astra
diff options
context:
space:
mode:
Diffstat (limited to 'include/astra')
-rw-r--r--include/astra/CudaFilteredBackProjectionAlgorithm.h28
-rw-r--r--include/astra/FanFlatBeamLineKernelProjector2D.inl302
-rw-r--r--include/astra/FanFlatProjectionGeometry2D.h5
-rw-r--r--include/astra/GeometryUtil2D.h69
-rw-r--r--include/astra/ParallelBeamBlobKernelProjector2D.h7
-rw-r--r--include/astra/ParallelBeamBlobKernelProjector2D.inl312
-rw-r--r--include/astra/ParallelBeamLineKernelProjector2D.h2
-rw-r--r--include/astra/ParallelBeamLineKernelProjector2D.inl403
-rw-r--r--include/astra/ParallelBeamLinearKernelProjector2D.h2
-rw-r--r--include/astra/ParallelBeamLinearKernelProjector2D.inl255
-rw-r--r--include/astra/ParallelBeamStripKernelProjector2D.h1
-rw-r--r--include/astra/ParallelBeamStripKernelProjector2D.inl466
-rw-r--r--include/astra/ParallelProjectionGeometry2D.h16
-rw-r--r--include/astra/ParallelVecProjectionGeometry2D.h163
-rw-r--r--include/astra/ProjectionGeometry2D.h19
15 files changed, 1140 insertions, 910 deletions
diff --git a/include/astra/CudaFilteredBackProjectionAlgorithm.h b/include/astra/CudaFilteredBackProjectionAlgorithm.h
index 057843e..55191ef 100644
--- a/include/astra/CudaFilteredBackProjectionAlgorithm.h
+++ b/include/astra/CudaFilteredBackProjectionAlgorithm.h
@@ -25,30 +25,26 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-----------------------------------------------------------------------
*/
-#ifndef CUDAFILTEREDBACKPROJECTIONALGORITHM2_H
-#define CUDAFILTEREDBACKPROJECTIONALGORITHM2_H
+#ifndef CUDAFILTEREDBACKPROJECTIONALGORITHM_H
+#define CUDAFILTEREDBACKPROJECTIONALGORITHM_H
#ifdef ASTRA_CUDA
#include <astra/Float32ProjectionData2D.h>
#include <astra/Float32VolumeData2D.h>
-#include <astra/ReconstructionAlgorithm2D.h>
+#include <astra/CudaReconstructionAlgorithm2D.h>
#include "../../cuda/2d/astra.h"
namespace astra
{
-class _AstraExport CCudaFilteredBackProjectionAlgorithm : public CReconstructionAlgorithm2D
+class _AstraExport CCudaFilteredBackProjectionAlgorithm : public CCudaReconstructionAlgorithm2D
{
public:
static std::string type;
private:
- CFloat32ProjectionData2D * m_pSinogram;
- CFloat32VolumeData2D * m_pReconstruction;
- int m_iGPUIndex;
- int m_iPixelSuperSampling;
E_FBPFILTER m_eFilter;
float * m_pfFilter;
int m_iFilterWidth; // number of elements per projection direction in filter
@@ -65,15 +61,6 @@ public:
virtual bool initialize(const Config& _cfg);
bool initialize(CFloat32ProjectionData2D * _pSinogram, CFloat32VolumeData2D * _pReconstruction, E_FBPFILTER _eFilter, const float * _pfFilter = NULL, int _iFilterWidth = 0, int _iGPUIndex = -1, float _fFilterParameter = -1.0f);
- virtual void run(int _iNrIterations = 0);
-
- static int calcIdealRealFilterWidth(int _iDetectorCount);
- static int calcIdealFourierFilterWidth(int _iDetectorCount);
-
- //debug
- static void testGenFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount);
- static int getGPUCount();
-
/** Get a description of the class.
*
* @return description string
@@ -83,12 +70,7 @@ public:
protected:
bool check();
- AstraFBP* m_pFBP;
-
- bool m_bAstraFBPInit;
-
- void initializeFromProjector();
- virtual bool requiresProjector() const { return false; }
+ virtual void initCUDAAlgorithm();
};
// inline functions
diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl
index d967844..51dc878 100644
--- a/include/astra/FanFlatBeamLineKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl
@@ -25,6 +25,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-----------------------------------------------------------------------
*/
+#define policy_weight(p,rayindex,volindex,weight) if (p.pixelPrior(volindex)) { p.addWeight(rayindex, volindex, weight); p.pixelPosterior(volindex); }
template <typename Policy>
void CFanFlatBeamLineKernelProjector2D::project(Policy& p)
@@ -48,246 +49,165 @@ void CFanFlatBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _
}
//----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+// PROJECT BLOCK - vector projection geometry
template <typename Policy>
void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
{
- // variables
- float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2;
- float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY;
- int iVolumeIndex, iRayIndex, row, col, iAngle, iDetector, x1;
- bool switch_t;
+ // get vector geometry
+ const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry;
+ if (dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry)) {
+ pVecProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry();
+ } else {
+ pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry);
+ }
+
+ // precomputations
+ const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX();
+ const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY();
+ const float32 inv_pixelLengthX = 1.0f / pixelLengthX;
+ const float32 inv_pixelLengthY = 1.0f / pixelLengthY;
+ const int colCount = m_pVolumeGeometry->getGridColCount();
+ const int rowCount = m_pVolumeGeometry->getGridRowCount();
+ const int detCount = pVecProjectionGeometry->getDetectorCount();
+ const float32 Ex = m_pVolumeGeometry->getWindowMinY() + pixelLengthX*0.5f;
+ const float32 Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
- const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry);
- const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry);
+ // loop angles
+ #pragma omp parallel for
+ for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
- float32 old_theta, theta, alpha;
- const SFanProjection * proj = 0;
+ // variables
+ float32 Dx, Dy, Rx, Ry, S, T, weight, c, r, deltac, deltar, offset, RxOverRy, RyOverRx;
+ float32 lengthPerRow, lengthPerCol, invTminSTimesLengthPerRow, invTminSTimesLengthPerCol;
+ int iVolumeIndex, iRayIndex, row, col, iDetector;
- // loop angles
- for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
- // get theta
- if (pProjectionGeometry) {
- old_theta = pProjectionGeometry->getProjectionAngle(iAngle);
- }
- else if (pVecProjectionGeometry) {
- proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
- old_theta = atan2(-proj->fSrcX, proj->fSrcY);
- if (old_theta < 0) old_theta += 2*PI;
- } else {
- assert(false);
- }
-
- switch_t = false;
- if (old_theta >= 7*PIdiv4) old_theta -= 2*PI;
- if (old_theta >= 3*PIdiv4) {
- old_theta -= PI;
- switch_t = true;
- }
+ const SFanProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
// loop detectors
for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
- iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
+ iRayIndex = iAngle * detCount + iDetector;
// POLICY: RAY PRIOR
if (!p.rayPrior(iRayIndex)) continue;
-
- // get values
- if (pProjectionGeometry) {
- t = -pProjectionGeometry->indexToDetectorOffset(iDetector);
- alpha = atan(t / pProjectionGeometry->getSourceDetectorDistance());
- t = sin(alpha) * pProjectionGeometry->getOriginSourceDistance();
- }
- else if (pVecProjectionGeometry) {
- float32 detX = proj->fDetSX + proj->fDetUX*(0.5f + iDetector);
- float32 detY = proj->fDetSY + proj->fDetUY*(0.5f + iDetector);
- alpha = angleBetweenVectors(-proj->fSrcX, -proj->fSrcY, detX - proj->fSrcX, detY - proj->fSrcY);
- t = sin(alpha) * sqrt(proj->fSrcX*proj->fSrcX + proj->fSrcY*proj->fSrcY);
+
+ Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+ Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
+
+ Rx = proj->fSrcX - Dx;
+ Ry = proj->fSrcY - Dy;
+
+ bool vertical = fabs(Rx) < fabs(Ry);
+ if (vertical) {
+ RxOverRy = Rx/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);
+ invTminSTimesLengthPerRow = lengthPerRow / (T - S);
} else {
- assert(false);
+ RyOverRx = Ry/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);
+ invTminSTimesLengthPerCol = lengthPerCol / (T - S);
}
- if (switch_t) t = -t;
- theta = old_theta + alpha;
-
- // precalculate sin, cos, 1/cos
- sin_theta = sin(theta);
- cos_theta = cos(theta);
- inv_sin_theta = 1.0f / sin_theta;
- inv_cos_theta = 1.0f / cos_theta;
-
- // precalculate kernel limits
- lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
- updatePerRow = sin_theta * inv_cos_theta;
- inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
-
- // precalculate kernel limits
- lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
- updatePerCol = cos_theta * inv_sin_theta;
- inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
-
- // precalculate S and T
- S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
- T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
+ bool isin = false;
// vertically
- if (old_theta <= PIdiv4) {
-
- // calculate x for row 0
- P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
- x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX;
+ if (vertical) {
+
+ // calculate c for row 0
+ c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
// for each row
- for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-
- // get coords
- x1 = int((x > 0.0f) ? x : x-1.0f);
- x2 = x - x1;
- x += updatePerRow;
+ for (row = 0; row < rowCount; ++row, c += deltac) {
- if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
+ col = int(c+0.5f);
+ if (col <= 0 || col >= colCount-1) { if (!isin) continue; else break; }
+ offset = c - float32(col);
// left
- if (x2 < 0.5f-S) {
- I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
- if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridColCount()*/) {//x1 is always less than or equal to gridColCount because of the "continue" in the beginning of the for-loop
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
-
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- }
+ if (offset < -S) {
+ weight = (offset + T) * invTminSTimesLengthPerRow;
+
+ iVolumeIndex = row * colCount + col - 1;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow-weight)
- // center
- else if (x2 <= 0.5f+S) {
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ iVolumeIndex++;
+ policy_weight(p, iRayIndex, iVolumeIndex, weight)
}
// right
- else if (x2 <= 1.0f) {
- I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridColCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ else if (S < offset) {
+ weight = (offset - S) * invTminSTimesLengthPerRow;
+
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow-weight)
+
+ iVolumeIndex++;
+ policy_weight(p, iRayIndex, iVolumeIndex, weight)
+ }
+
+ // centre
+ else {
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow)
}
+ isin = true;
}
}
// horizontally
- //else if (PIdiv4 <= old_theta && old_theta <= 3*PIdiv4) {
else {
- // calculate point P
- P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
- x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY;
+ // calculate r for col 0
+ r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
// for each col
- for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
- // get coords
- x1 = int((x > 0.0f) ? x : x-1.0f);
- x2 = x - x1;
- x += updatePerCol;
+ for (col = 0; col < colCount; ++col, r += deltar) {
- if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
+ int row = int(r+0.5f);
+ if (row <= 0 || row >= rowCount-1) { if (!isin) continue; else break; }
+ offset = r - float32(row);
// up
- if (x2 < 0.5f-T) {
- I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
- if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridRowCount()*/) {//x1 is always less than or equal to gridRowCount because of the "continue" in the beginning of the for-loop
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
-
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- }
+ if (offset < -S) {
+ weight = (offset + T) * invTminSTimesLengthPerCol;
+
+ iVolumeIndex = (row-1) * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight)
- // center
- else if (x2 <= 0.5f+T) {
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ iVolumeIndex += colCount;
+ policy_weight(p, iRayIndex, iVolumeIndex, weight)
}
// down
- else if (x2 <= 1.0f) {
- I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridRowCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ else if (S < offset) {
+ weight = (offset - S) * invTminSTimesLengthPerCol;
+
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight)
+
+ iVolumeIndex += colCount;
+ policy_weight(p, iRayIndex, iVolumeIndex, weight)
+ }
+
+ // centre
+ else {
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol)
}
+ isin = true;
}
- } // end loop col
+ }
// POLICY: RAY POSTERIOR
p.rayPosterior(iRayIndex);
} // end loop detector
} // end loop angles
+
}
diff --git a/include/astra/FanFlatProjectionGeometry2D.h b/include/astra/FanFlatProjectionGeometry2D.h
index e9a0535..c45eb78 100644
--- a/include/astra/FanFlatProjectionGeometry2D.h
+++ b/include/astra/FanFlatProjectionGeometry2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#define _INC_ASTRA_FANFLATPROJECTIONGEOMETRY2D
#include "ProjectionGeometry2D.h"
+#include "FanFlatVecProjectionGeometry2D.h"
#include <cmath>
@@ -189,6 +190,10 @@ public:
* @return a unit vector describing the direction
*/
virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex);
+
+ /** Create a vector geom
+ */
+ CFanFlatVecProjectionGeometry2D* toVectorGeometry();
};
diff --git a/include/astra/GeometryUtil2D.h b/include/astra/GeometryUtil2D.h
index 6434d3c..4d79353 100644
--- a/include/astra/GeometryUtil2D.h
+++ b/include/astra/GeometryUtil2D.h
@@ -30,17 +30,74 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
namespace astra {
+struct SParProjection {
+ // the ray direction
+ float fRayX, fRayY;
+
+ // the start of the (linear) detector
+ float fDetSX, fDetSY;
+
+ // the length of a single detector pixel
+ float fDetUX, fDetUY;
+
+
+ void translate(double dx, double dy) {
+ fDetSX += dx;
+ fDetSY += dy;
+ }
+ void scale(double factor) {
+ fRayX *= factor;
+ fRayY *= factor;
+ fDetSX *= factor;
+ fDetSY *= factor;
+ fDetUX *= factor;
+ fDetUY *= factor;
+ }
+};
+
+
struct SFanProjection {
- // the source
- float fSrcX, fSrcY;
+ // the source
+ float fSrcX, fSrcY;
- // the start of the (linear) detector
- float fDetSX, fDetSY;
+ // the start of the (linear) detector
+ float fDetSX, fDetSY;
- // the length of a single detector pixel
- float fDetUX, fDetUY;
+ // the length of a single detector pixel
+ float fDetUX, fDetUY;
+
+ void translate(double dx, double dy) {
+ fSrcX += dx;
+ fSrcY += dy;
+ fDetSX += dx;
+ fDetSY += dy;
+ }
+ void scale(double factor) {
+ fSrcX *= factor;
+ fSrcY *= factor;
+ fDetSX *= factor;
+ fDetSY *= factor;
+ fDetUX *= factor;
+ fDetUY *= factor;
+ }
};
+
+
+SParProjection* genParProjections(unsigned int iProjAngles,
+ unsigned int iProjDets,
+ double fDetSize,
+ const float *pfAngles,
+ const float *pfExtraOffsets);
+
+SFanProjection* genFanProjections(unsigned int iProjAngles,
+ unsigned int iProjDets,
+ double fOriginSource, double fOriginDetector,
+ double fDetSize,
+ const float *pfAngles);
+
+bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset);
+
}
#endif
diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.h b/include/astra/ParallelBeamBlobKernelProjector2D.h
index 529cc10..12bee5f 100644
--- a/include/astra/ParallelBeamBlobKernelProjector2D.h
+++ b/include/astra/ParallelBeamBlobKernelProjector2D.h
@@ -214,7 +214,12 @@ protected:
float32 m_fBlobSampleRate; //< At which interval are the inserted blob values evaluated?
int m_iBlobSampleCount; //< Number of evaluated blob samples
float32* m_pfBlobValues; //< Evaluated blob values
- float32* m_pfBlobValuesNeg; //< Evaluated blob values
+
+ /** Internal policy-based projection of a range of angles and range.
+ * (_i*From is inclusive, _i*To exclusive) */
+ template <typename Policy>
+ void projectBlock_internal(int _iProjFrom, int _iProjTo,
+ int _iDetFrom, int _iDetTo, Policy& _policy);
};
diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.inl b/include/astra/ParallelBeamBlobKernelProjector2D.inl
index c5040b7..e11c4d4 100644
--- a/include/astra/ParallelBeamBlobKernelProjector2D.inl
+++ b/include/astra/ParallelBeamBlobKernelProjector2D.inl
@@ -26,186 +26,206 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
*/
-
-//----------------------------------------------------------------------------------------
-// PROJECT ALL
template <typename Policy>
void CParallelBeamBlobKernelProjector2D::project(Policy& p)
{
- for (int iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
- for (int iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
- projectSingleRay(iAngle, iDetector, p);
- }
- }
+ projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+ 0, m_pProjectionGeometry->getDetectorCount(), p);
}
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
template <typename Policy>
void CParallelBeamBlobKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
{
- for (int iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
- projectSingleRay(_iProjection, iDetector, p);
- }
+ projectBlock_internal(_iProjection, _iProjection + 1,
+ 0, m_pProjectionGeometry->getDetectorCount(), p);
}
-
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
template <typename Policy>
void CParallelBeamBlobKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
{
- ASTRA_ASSERT(m_bIsInitialized);
+ projectBlock_internal(_iProjection, _iProjection + 1,
+ _iDetector, _iDetector + 1, p);
+}
- int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
+//----------------------------------------------------------------------------------------
+// PROJECT BLOCK - vector projection geometry
+//
+// Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+//
+// For each angle/detector pair:
+//
+// Let D=(Dx,Dy) denote the centre of the detector (point) in volume coordinates, and
+// let R=(Rx,Ry) denote the direction of the ray (vector).
+//
+// For mainly vertical rays (|Rx|<=|Ry|),
+// let E=(Ex,Ey) denote the centre of the most upper left pixel:
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// and let F=(Fx,Fy) denote a vector to the next pixel
+// F = (PixelLengthX, 0)
+//
+// The intersection of the ray (D+aR) with the centre line of the upper row of pixels (E+bF) is
+// { Dx + a*Rx = Ex + b*Fx
+// { Dy + a*Ry = Ey + b*Fy
+// Solving for (a,b) results in:
+// a = (Ey + b*Fy - Dy)/Ry
+// = (Ey - Dy)/Ry
+// b = (Dx + a*Rx - Ex)/Fx
+// = (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+//
+// Define c as the x-value of the intersection of the ray with the upper row in pixel coordinates.
+// c = b
+//
+// The intersection of the ray (D+aR) with the centre line of the second row of pixels (E'+bF) with
+// E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+// expressed in x-value pixel coordinates is
+// c' = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx.
+// And thus:
+// deltac = c' - c = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx - (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+// = [(Ey' - Dy)*Rx/Ry - (Ey - Dy)*Rx/Ry]/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = -PixelLengthY*(Rx/Ry)/Fx.
+//
+// Given c on a certain row, its pixel directly on its left (col), and the distance (offset) to it, can be found:
+// col = floor(c)
+// offset = c - col
+//
+// The index of this pixel is
+// volumeIndex = row * colCount + col
+//
+//
+// Mainly horizontal rays (|Rx|<=|Ry|) are handled in a similar fashion:
+//
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// F = (0, -PixelLengthX)
+//
+// a = (Ex + b*Fx - Dx)/Rx = (Ex - Dx)/Rx
+// b = (Dy + a*Ry - Ey)/Fy = (Dy + (Ex - Dx)*Ry/Rx - Ey)/Fy
+// r = b
+// deltar = PixelLengthX*(Ry/Rx)/Fy.
+// row = floor(r+1/2)
+// offset = r - row
+//
+template <typename Policy>
+void CParallelBeamBlobKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
+{
+ // get vector geometry
+ const CParallelVecProjectionGeometry2D* pVecProjectionGeometry;
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)) {
+ pVecProjectionGeometry = dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry();
+ } else {
+ pVecProjectionGeometry = dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry);
+ }
- // POLICY: RAY PRIOR
- if (!p.rayPrior(iRayIndex)) return;
+ // precomputations
+ const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX();
+ const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY();
+ const float32 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
+ const float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
+ const int colCount = m_pVolumeGeometry->getGridColCount();
+ const int rowCount = m_pVolumeGeometry->getGridRowCount();
+ const int detCount = pVecProjectionGeometry->getDetectorCount();
+
+ // loop angles
+ #pragma omp parallel for
+ for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+ // variables
+ float32 Dx, Dy, Ex, Ey, c, r, deltac, deltar, offset, invBlobExtent, RxOverRy, RyOverRx;
+ int iVolumeIndex, iRayIndex, row, col, iDetector;
+ int col_left, col_right, row_top, row_bottom, index;
+
+ const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+ bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
+ if (vertical) {
+ RxOverRy = proj->fRayX/proj->fRayY;
+ deltac = -m_pVolumeGeometry->getPixelLengthY() * (proj->fRayX/proj->fRayY) * inv_pixelLengthX;
+ invBlobExtent = m_pVolumeGeometry->getPixelLengthY() / abs(m_fBlobSize * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / proj->fRayY);
+ } else {
+ RyOverRx = proj->fRayY/proj->fRayX;
+ deltar = -m_pVolumeGeometry->getPixelLengthX() * (proj->fRayY/proj->fRayX) * inv_pixelLengthY;
+ invBlobExtent = m_pVolumeGeometry->getPixelLengthX() / abs(m_fBlobSize * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / proj->fRayX);
+ }
- // get values
- float32 t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector);
- float32 theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
- if (theta >= 7*PIdiv4) theta -= 2*PI;
+ Ex = m_pVolumeGeometry->getWindowMinY() + pixelLengthX*0.5f;
+ Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
- bool flip = false;
+ // loop detectors
+ for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
+
+ iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
- if (theta >= 3*PIdiv4) {
- theta -= PI;
- t = -t;
- flip = true;
- }
+ // POLICY: RAY PRIOR
+ if (!p.rayPrior(iRayIndex)) continue;
+
+ Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+ Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
+ // vertically
+ if (vertical) {
- if (theta <= PIdiv4) { // -pi/4 <= theta <= pi/4
+ // calculate c for row 0
+ c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
- // precalculate sin, cos, 1/cos
- float32 sin_theta = sin(theta);
- float32 cos_theta = cos(theta);
- float32 inv_cos_theta = 1.0f / cos_theta;
+ // loop rows
+ for (row = 0; row < rowCount; ++row, c += deltac) {
- // precalculate other stuff
- float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
- float32 updatePerRow = sin_theta * lengthPerRow;
- float32 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
- float32 pixelLengthX_over_blobSize = m_pVolumeGeometry->getPixelLengthX() / m_fBlobSize;
-
- // some variables
- int row, col, xmin, xmax;
- float32 P, x, d;
+ col_left = int(c - 0.5f - m_fBlobSize);
+ col_right = int(c + 0.5f + m_fBlobSize);
- // calculate P and x for row 0
- P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
- x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f;
+ if (col_left < 0) col_left = 0;
+ if (col_right > colCount-1) col_right = colCount-1;
- // for each row
- for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-
- // calculate extent
- xmin = (int)ceil((P - m_fBlobSize - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f);
- xmax = (int)floor((P + m_fBlobSize - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f);
-
- // add pixels
- for (col = xmin; col <= xmax; col++) {
- if (col >= 0 && col < m_pVolumeGeometry->getGridColCount()) {
- //d = abs(x - col) * pixelLengthX_over_blobSize;
- //index = (int)(d*m_iBlobSampleCount+0.5f);
- //float32 fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)] * lengthPerRow;
-
- float32 fWeight;
- int index;
- if ((x >= col) ^ flip) {
- d = abs(x - col) * pixelLengthX_over_blobSize * cos_theta;
- index = (int)(d*m_iBlobSampleCount+0.5f);
- fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)];
- } else {
- d = abs(x - col) * pixelLengthX_over_blobSize * cos_theta;
- index = (int)(d*m_iBlobSampleCount+0.5f);
- fWeight = m_pfBlobValuesNeg[min(index,m_iBlobSampleCount-1)];
- }
+ // loop columns
+ for (col = col_left; col <= col_right; ++col) {
- int iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, fWeight);
- p.pixelPosterior(iVolumeIndex);
+ iVolumeIndex = row * colCount + col;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ offset = abs(c - float32(col)) * invBlobExtent;
+ index = (int)(offset*m_iBlobSampleCount+0.5f);
+ p.addWeight(iRayIndex, iVolumeIndex, m_pfBlobValues[min(index,m_iBlobSampleCount-1)]);
+ p.pixelPosterior(iVolumeIndex);
+ }
}
}
}
- // update P and x
- P += updatePerRow;
- x += updatePerRow * inv_pixelLengthX;
- }
+ // horizontally
+ else {
- } else { // pi/4 < theta < 3pi/4
-
- // precalculate sin cos
- float32 sin_90_theta = sin(PIdiv2-theta);
- float32 cos_90_theta = cos(PIdiv2-theta);
- float32 inv_cos_90_theta = 1.0f / cos_90_theta;
-
- // precalculate other stuff
- float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_cos_90_theta;
- float32 updatePerCol = sin_90_theta * lengthPerCol;
- float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
- float32 pixelLengthY_over_blobSize = m_pVolumeGeometry->getPixelLengthY() / m_fBlobSize;
-
- // some variables
- int row, col, xmin, xmax;
- float32 P,x, d;
-
- // calculate P and x for col 0
- P = (sin_90_theta * m_pVolumeGeometry->pixelColToCenterX(0) - t) * inv_cos_90_theta;
- x = (P - m_pVolumeGeometry->getWindowMinY()) * inv_pixelLengthY - 0.5f;
-
- // for each col
- for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
- // calculate extent
- xmin = (int)ceil((P - m_fBlobSize - m_pVolumeGeometry->getWindowMinY()) * inv_pixelLengthY - 0.5f);
- xmax = (int)floor((P + m_fBlobSize - m_pVolumeGeometry->getWindowMinY()) * inv_pixelLengthY - 0.5f);
-
- // add pixels
- for (row = xmin; row <= xmax; row++) {
- if (row >= 0 && row < m_pVolumeGeometry->getGridRowCount()) {
- //d = abs(x - row) * pixelLengthY_over_blobSize;
- //int index = (int)(d*m_iBlobSampleCount+0.5f);
- //float32 fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)] * lengthPerCol;
-
- float32 fWeight;
- int index;
- if ((x <= row) ^ flip) {
- d = abs(x - row) * pixelLengthY_over_blobSize * cos_90_theta;
- index = (int)(d*m_iBlobSampleCount+0.5f);
- fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)];
- } else {
- d = abs(x - row) * pixelLengthY_over_blobSize * cos_90_theta;
- index = (int)(d*m_iBlobSampleCount+0.5f);
- fWeight = m_pfBlobValuesNeg[min(index,m_iBlobSampleCount-1)];
- }
-
-
- int iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, fWeight);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- }
+ // calculate r for col 0
+ r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
- // update P and x
- P += updatePerCol;
- x += updatePerCol * inv_pixelLengthY;
- }
+ // loop columns
+ for (col = 0; col < colCount; ++col, r += deltar) {
- }
+ row_top = int(r - 0.5f - m_fBlobSize);
+ row_bottom = int(r + 0.5f + m_fBlobSize);
- // POLICY: RAY POSTERIOR
- p.rayPosterior(iRayIndex);
+ if (row_top < 0) row_top = 0;
+ if (row_bottom > rowCount-1) row_bottom = rowCount-1;
+ // loop rows
+ for (row = row_top; row <= row_bottom; ++row) {
+ iVolumeIndex = row * colCount + col;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ offset = abs(r - float32(row)) * invBlobExtent;
+ index = (int)(offset*m_iBlobSampleCount+0.5f);
+ p.addWeight(iRayIndex, iVolumeIndex, m_pfBlobValues[min(index,m_iBlobSampleCount-1)]);
+ p.pixelPosterior(iVolumeIndex);
+ }
+ }
+ }
+ }
+
+ // POLICY: RAY POSTERIOR
+ p.rayPosterior(iRayIndex);
+
+ } // end loop detector
+ } // end loop angles
}
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.h b/include/astra/ParallelBeamLineKernelProjector2D.h
index 0f25d83..e0b7b46 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.h
+++ b/include/astra/ParallelBeamLineKernelProjector2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#define _INC_ASTRA_PARALLELBEAMLINEKERNELPROJECTOR
#include "ParallelProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
#include "Float32Data2D.h"
#include "Projector2D.h"
@@ -179,6 +180,7 @@ protected:
template <typename Policy>
void projectBlock_internal(int _iProjFrom, int _iProjTo,
int _iDetFrom, int _iDetTo, Policy& _policy);
+
};
inline std::string CParallelBeamLineKernelProjector2D::getType()
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl
index dcd11bd..83c16d7 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLineKernelProjector2D.inl
@@ -24,12 +24,13 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-----------------------------------------------------------------------
*/
+#define policy_weight(p,rayindex,volindex,weight) if (p.pixelPrior(volindex)) { p.addWeight(rayindex, volindex, weight); p.pixelPosterior(volindex); }
template <typename Policy>
void CParallelBeamLineKernelProjector2D::project(Policy& p)
{
projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
- 0, m_pProjectionGeometry->getDetectorCount(), p);
+ 0, m_pProjectionGeometry->getDetectorCount(), p);
}
template <typename Policy>
@@ -48,235 +49,265 @@ void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int
//----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+// PROJECT BLOCK - vector projection geometry
+//
+// Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+//
+// For each angle/detector pair:
+//
+// Let D=(Dx,Dy) denote the centre of the detector (point) in volume coordinates, and
+// let R=(Rx,Ry) denote the direction of the ray (vector).
+//
+// For mainly vertical rays (|Rx|<=|Ry|),
+// let E=(Ex,Ey) denote the centre of the most upper left pixel:
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// and let F=(Fx,Fy) denote a vector to the next pixel
+// F = (PixelLengthX, 0)
+//
+// The intersection of the ray (D+aR) with the centre line of the upper row of pixels (E+bF) is
+// { Dx + a*Rx = Ex + b*Fx
+// { Dy + a*Ry = Ey + b*Fy
+// Solving for (a,b) results in:
+// a = (Ey + b*Fy - Dy)/Ry
+// = (Ey - Dy)/Ry
+// b = (Dx + a*Rx - Ex)/Fx
+// = (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+//
+// Define c as the x-value of the intersection of the ray with the upper row in pixel coordinates.
+// c = b
+//
+// The intersection of the ray (D+aR) with the centre line of the second row of pixels (E'+bF) with
+// E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+// expressed in x-value pixel coordinates is
+// c' = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx.
+// And thus:
+// deltac = c' - c = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx - (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+// = [(Ey' - Dy)*Rx/Ry - (Ey - Dy)*Rx/Ry]/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = -PixelLengthY*(Rx/Ry)/Fx.
+//
+// Given c on a certain row, its closest pixel (col), and the distance (offset) to it, can be found:
+// col = floor(c+1/2)
+// offset = c - col
+//
+// The index of this pixel is
+// volumeIndex = row * colCount + col
+//
+// The projection kernel is defined by
+//
+// _____ LengthPerRow
+// /| | |\
+// / | | | \
+// __/ | | | \__ 0
+// -T -S 0 S T
+//
+// with S = 1/2 - 1/2*|Rx/Ry|, T = 1/2 + 1/2*|Rx/Ry|, and LengthPerRow = pixelLengthX * sqrt(Rx^2+Ry^2) / |Ry|
+//
+// And thus
+// { (offset+T)/(T-S) * LengthPerRow if -T <= offset < S
+// W_(rayIndex,volIndex) = { LengthPerRow if -S <= offset <= S
+// { (offset-S)/(T-S) * LengthPerRow if S < offset <= T
+//
+// If -T <= offset < S, the weight for the pixel directly to the left is
+// W_(rayIndex,volIndex-1) = LengthPerRow - (offset+T)/(T-S) * LengthPerRow,
+// and if S < offset <= T, the weight for the pixel directly to the right is
+// W_(rayIndex,volIndex+1) = LengthPerRow - (offset-S)/(T-S) * LengthPerRow.
+//
+//
+// Mainly horizontal rays (|Rx|<=|Ry|) are handled in a similar fashion:
+//
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// F = (0, -PixelLengthX)
+//
+// a = (Ex + b*Fx - Dx)/Rx = (Ex - Dx)/Rx
+// b = (Dy + a*Ry - Ey)/Fy = (Dy + (Ex - Dx)*Ry/Rx - Ey)/Fy
+// r = b
+// deltar = PixelLengthX*(Ry/Rx)/Fy.
+// row = floor(r+1/2)
+// offset = r - row
+// S = 1/2 - 1/2*|Ry/Rx|
+// T = 1/2 + 1/2*|Ry/Rx|
+// LengthPerCol = pixelLengthY * sqrt(Rx^2+Ry^2) / |Rx|
+//
+// { (offset+T)/(T-S) * LengthPerCol if -T <= offset < S
+// W_(rayIndex,volIndex) = { LengthPerCol if -S <= offset <= S
+// { (offset-S)/(T-S) * LengthPerCol if S < offset <= T
+//
+// W_(rayIndex,volIndex-colcount) = LengthPerCol - (offset+T)/(T-S) * LengthPerCol
+// W_(rayIndex,volIndex+colcount) = LengthPerCol - (offset-S)/(T-S) * LengthPerCol
+//
template <typename Policy>
void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
{
- // variables
- float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2;
- float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY;
- int iVolumeIndex, iRayIndex, row, col, iAngle, iDetector, x1;
- bool switch_t;
+ // get vector geometry
+ const CParallelVecProjectionGeometry2D* pVecProjectionGeometry;
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)) {
+ pVecProjectionGeometry = dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry();
+ } else {
+ pVecProjectionGeometry = dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry);
+ }
+
+ // precomputations
+ const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX();
+ const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY();
+ const float32 inv_pixelLengthX = 1.0f / pixelLengthX;
+ const float32 inv_pixelLengthY = 1.0f / pixelLengthY;
+ const int colCount = m_pVolumeGeometry->getGridColCount();
+ const int rowCount = m_pVolumeGeometry->getGridRowCount();
+ const int detCount = pVecProjectionGeometry->getDetectorCount();
// loop angles
- for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
- // get theta
- theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
- switch_t = false;
- if (theta >= 7*PIdiv4) theta -= 2*PI;
- if (theta >= 3*PIdiv4) {
- theta -= PI;
- switch_t = true;
+ #pragma omp parallel for
+ for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+ // variables
+ float32 Dx, Dy, Ex, Ey, S, T, weight, c, r, deltac, deltar, offset;
+ float32 RxOverRy, RyOverRx, lengthPerRow, lengthPerCol, invTminSTimesLengthPerRow, invTminSTimesLengthPerCol;
+ int iVolumeIndex, iRayIndex, row, col, iDetector;
+
+ 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);
}
- // precalculate sin, cos, 1/cos
- sin_theta = sin(theta);
- cos_theta = cos(theta);
- inv_sin_theta = 1.0f / sin_theta;
- inv_cos_theta = 1.0f / cos_theta;
-
- // precalculate kernel limits
- lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
- updatePerRow = sin_theta * inv_cos_theta;
- inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
-
- // precalculate kernel limits
- lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
- updatePerCol = cos_theta * inv_sin_theta;
- inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
-
- // precalculate S and T
- S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
- T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
+ Ex = m_pVolumeGeometry->getWindowMinY() + pixelLengthX*0.5f;
+ Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
// loop detectors
- for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
-
+ for (int iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
+
iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
// POLICY: RAY PRIOR
if (!p.rayPrior(iRayIndex)) continue;
-
- // get t
- t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
- if (switch_t) t = -t;
-
- // vertically
- if (theta <= PIdiv4) {
-
- // calculate x for row 0
- P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
- x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX;
- // get coords
- int nextx1 = int((x > 0.0f) ? x : x-1.0f);
- float nextx2 = x - nextx1;
+ Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+ Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
- // for each row
- for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
+ bool isin = false;
+
+ // vertically
+ if (vertical) {
- x1 = nextx1;
- x2 = nextx2;
+ // calculate c for row 0
+ c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
- nextx2 += updatePerRow;
- while (nextx2 >= 1.0f) {
- nextx2 -= 1.0f;
- nextx1++;
- }
- while (nextx2 < 0.0f) {
- nextx2 += 1.0f;
- nextx1--;
- }
+ // loop rows
+ for (row = 0; row < rowCount; ++row, c += deltac) {
- if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
+ col = int(c+0.5f);
+ if (col <= 0 || col >= colCount-1) { if (!isin) continue; else break; }
+ offset = c - float32(col);
// left
- if (x2 < 0.5f-S) {
- I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
- if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
- p.pixelPosterior(iVolumeIndex);
- }
+ if (offset < -S) {
+ weight = (offset + T) * invTminSTimesLengthPerRow;
+
+ iVolumeIndex = row * colCount + col - 1;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-weight);
+ p.pixelPosterior(iVolumeIndex);
}
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
+ iVolumeIndex++;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ p.addWeight(iRayIndex, iVolumeIndex, weight);
+ p.pixelPosterior(iVolumeIndex);
}
}
- // center
- else if (x2 <= 0.5f+S) {
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- }
-
// right
- else {
- I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
+ else if (S < offset) {
+ weight = (offset - S) * invTminSTimesLengthPerRow;
+
+ iVolumeIndex = row * colCount + col;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-weight);
+ p.pixelPosterior(iVolumeIndex);
}
- if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
- p.pixelPosterior(iVolumeIndex);
- }
+
+ iVolumeIndex++;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ p.addWeight(iRayIndex, iVolumeIndex, weight);
+ p.pixelPosterior(iVolumeIndex);
}
}
+
+ // centre
+ else {
+ iVolumeIndex = row * colCount + col;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+ p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow);
+ p.pixelPosterior(iVolumeIndex);
+ }
+ }
+ isin = true;
}
}
// horizontally
- else if (PIdiv4 <= theta && theta <= 3*PIdiv4) {
+ else {
- // calculate point P
- P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
- x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY;
+ // calculate r for col 0
+ r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
- // get coords
- int nextx1 = int((x > 0.0f) ? x : x-1.0f);
- float nextx2 = x - nextx1;
+ // loop columns
+ for (col = 0; col < colCount; ++col, r += deltar) {
- // for each col
- for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
+ row = int(r+0.5f);
+ if (row <= 0 || row >= rowCount-1) { if (!isin) continue; else break; }
+ offset = r - float32(row);
- x1 = nextx1;
- x2 = nextx2;
+ // up
+ if (offset < -S) {
+ weight = (offset + T) * invTminSTimesLengthPerCol;
- nextx2 += updatePerCol;
- while (nextx2 >= 1.0f) {
- nextx2 -= 1.0f;
- nextx1++;
- }
- while (nextx2 < 0.0f) {
- nextx2 += 1.0f;
- nextx1--;
+ iVolumeIndex = (row-1) * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight)
+
+ iVolumeIndex += colCount;
+ policy_weight(p, iRayIndex, iVolumeIndex, weight)
}
- if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
+ // down
+ else if (S < offset) {
+ weight = (offset - S) * invTminSTimesLengthPerCol;
- // up
- if (x2 < 0.5f-T) {
- I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
- if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight)
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ iVolumeIndex += colCount;
+ policy_weight(p, iRayIndex, iVolumeIndex, weight)
}
- // center
- else if (x2 <= 0.5f+T) {
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- }
-
- // down
- else {
- I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ // centre
+ else {
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol)
}
+ isin = true;
}
- } // end loop col
+ }
// POLICY: RAY POSTERIOR
p.rayPosterior(iRayIndex);
@@ -284,5 +315,7 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
} // end loop detector
} // end loop angles
-}
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+ delete pVecProjectionGeometry;
+}
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.h b/include/astra/ParallelBeamLinearKernelProjector2D.h
index 7b40628..3e81fa3 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.h
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#define _INC_ASTRA_PARALLELLINEARKERNELPROJECTOR
#include "ParallelProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
#include "Float32Data2D.h"
#include "Projector2D.h"
@@ -184,7 +185,6 @@ protected:
void projectBlock_internal(int _iProjFrom, int _iProjTo,
int _iDetFrom, int _iDetTo, Policy& _policy);
-
};
//----------------------------------------------------------------------------------------
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl
index 5dd4781..2619a12 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -25,12 +25,13 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-----------------------------------------------------------------------
*/
+#define policy_weight(p,rayindex,volindex,weight) if (p.pixelPrior(volindex)) { p.addWeight(rayindex, volindex, weight); p.pixelPosterior(volindex); }
template <typename Policy>
void CParallelBeamLinearKernelProjector2D::project(Policy& p)
{
projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
- 0, m_pProjectionGeometry->getDetectorCount(), p);
+ 0, m_pProjectionGeometry->getDetectorCount(), p);
}
template <typename Policy>
@@ -47,45 +48,127 @@ void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, in
_iDetector, _iDetector + 1, p);
}
+
+
//----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+// PROJECT BLOCK - vector projection geometry
+//
+// Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+//
+// For each angle/detector pair:
+//
+// Let D=(Dx,Dy) denote the centre of the detector (point) in volume coordinates, and
+// let R=(Rx,Ry) denote the direction of the ray (vector).
+//
+// For mainly vertical rays (|Rx|<=|Ry|),
+// let E=(Ex,Ey) denote the centre of the most upper left pixel:
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// and let F=(Fx,Fy) denote a vector to the next pixel
+// F = (PixelLengthX, 0)
+//
+// The intersection of the ray (D+aR) with the centre line of the upper row of pixels (E+bF) is
+// { Dx + a*Rx = Ex + b*Fx
+// { Dy + a*Ry = Ey + b*Fy
+// Solving for (a,b) results in:
+// a = (Ey + b*Fy - Dy)/Ry
+// = (Ey - Dy)/Ry
+// b = (Dx + a*Rx - Ex)/Fx
+// = (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+//
+// Define c as the x-value of the intersection of the ray with the upper row in pixel coordinates.
+// c = b
+//
+// The intersection of the ray (D+aR) with the centre line of the second row of pixels (E'+bF) with
+// E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+// expressed in x-value pixel coordinates is
+// c' = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx.
+// And thus:
+// deltac = c' - c = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx - (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+// = [(Ey' - Dy)*Rx/Ry - (Ey - Dy)*Rx/Ry]/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = -PixelLengthY*(Rx/Ry)/Fx.
+//
+// Given c on a certain row, its pixel directly on its left (col), and the distance (offset) to it, can be found:
+// col = floor(c)
+// offset = c - col
+//
+// The index of this pixel is
+// volumeIndex = row * colCount + col
+//
+// The projection kernel is defined by
+//
+// LengthPerRow
+// /|\
+// / | \
+// __/ | \__ 0
+// p0 p1 p2
+//
+// And thus
+// W_(rayIndex,volIndex) = (1 - offset) * lengthPerRow
+// W_(rayIndex,volIndex+1) = offset * lengthPerRow
+//
+//
+// Mainly horizontal rays (|Rx|<=|Ry|) are handled in a similar fashion:
+//
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// F = (0, -PixelLengthX)
+//
+// a = (Ex + b*Fx - Dx)/Rx = (Ex - Dx)/Rx
+// b = (Dy + a*Ry - Ey)/Fy = (Dy + (Ex - Dx)*Ry/Rx - Ey)/Fy
+// r = b
+// deltar = PixelLengthX*(Ry/Rx)/Fy.
+// row = floor(r+1/2)
+// offset = r - row
+// LengthPerCol = pixelLengthY * sqrt(Rx^2+Ry^2) / |Rx|
+//
+// W_(rayIndex,volIndex) = (1 - offset) * lengthPerCol
+// W_(rayIndex,volIndex+colcount) = offset * lengthPerCol
+//
template <typename Policy>
void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
{
- // variables
- float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, t;
- float32 lengthPerRow, updatePerRow;
- float32 lengthPerCol, updatePerCol;
- bool switch_t;
- int iAngle, iDetector, iVolumeIndex, iRayIndex;
- int row, col, x1;
- float32 P,x,x2;
+ // get vector geometry
+ const CParallelVecProjectionGeometry2D* pVecProjectionGeometry;
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)) {
+ pVecProjectionGeometry = dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry();
+ } else {
+ pVecProjectionGeometry = dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry);
+ }
+
+ // precomputations
+ const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX();
+ const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY();
+ const float32 inv_pixelLengthX = 1.0f / pixelLengthX;
+ const float32 inv_pixelLengthY = 1.0f / pixelLengthY;
+ const int colCount = m_pVolumeGeometry->getGridColCount();
+ const int rowCount = m_pVolumeGeometry->getGridRowCount();
+ const int detCount = pVecProjectionGeometry->getDetectorCount();
// loop angles
- for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
- // get theta
- theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
- switch_t = false;
- if (theta >= 7*PIdiv4) theta -= 2*PI;
- if (theta >= 3*PIdiv4) {
- theta -= PI;
- switch_t = true;
+ #pragma omp parallel for
+ for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+ // variables
+ float32 Dx, Dy, Ex, Ey, x, y, c, r, deltac, deltar, offset;
+ float32 RxOverRy, RyOverRx, lengthPerRow, lengthPerCol;
+ int iVolumeIndex, iRayIndex, row, col, iDetector;
+
+ const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+ 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;
}
- // precalculate sin, cos, 1/cos
- sin_theta = sin(theta);
- cos_theta = cos(theta);
- inv_cos_theta = 1.0f / cos_theta;
- inv_sin_theta = 1.0f / sin_theta;
-
- // precalculate kernel limits
- lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
- updatePerRow = sin_theta * inv_cos_theta;
-
- // precalculate kernel limits
- lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
- updatePerCol = cos_theta * inv_sin_theta;
+ Ex = m_pVolumeGeometry->getWindowMinY() + pixelLengthX*0.5f;
+ Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
// loop detectors
for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
@@ -95,79 +178,54 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
// POLICY: RAY PRIOR
if (!p.rayPrior(iRayIndex)) continue;
- // get t
- t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
- if (switch_t) {
- t = -t;
- }
+ Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+ Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
+
+ bool isin = false;
// vertically
- if (theta <= PIdiv4) {
-
- // calculate x for row 0
- P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
- x = m_pVolumeGeometry->coordXToColF(P) - 0.5f;
+ if (vertical) {
+
+ // calculate c for row 0
+ c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
+
+ // loop rows
+ for (row = 0; row < rowCount; ++row, c += deltac) {
- // for each row
- for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
+ col = int(c);
+ if (col <= 0 || col >= colCount-1) { if (!isin) continue; else break; }
+ offset = c - float32(col);
+
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - offset) * lengthPerRow)
- // get coords
- x1 = int((x > 0.0f) ? x : x-1.0f);
- x2 = x - x1;
- x += updatePerRow;
-
- // add weights
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerRow);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, (x2) * lengthPerRow);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ iVolumeIndex++;
+ policy_weight(p, iRayIndex, iVolumeIndex, offset * lengthPerRow)
+
+ isin = true;
}
}
// horizontally
- else if (PIdiv4 <= theta && theta <= 3*PIdiv4) {
-
- // calculate point P
- P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
- x = m_pVolumeGeometry->coordYToRowF(P) - 0.5f;
-
- // for each row
- for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
- // get coords
- x1 = int((x > 0.0f) ? x : x-1.0f);
- x2 = x - x1;
- x += updatePerCol;
-
- // add weights
- if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerCol);
- p.pixelPosterior(iVolumeIndex);
- }
- }
- if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridRowCount()) {
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col);
- // POLICY: PIXEL PRIOR + ADD + POSTERIOR
- if (p.pixelPrior(iVolumeIndex)) {
- p.addWeight(iRayIndex, iVolumeIndex, x2 * lengthPerCol);
- p.pixelPosterior(iVolumeIndex);
- }
- }
+ else {
+
+ // calculate r for col 0
+ r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
+
+ // loop columns
+ for (col = 0; col < colCount; ++col, r += deltar) {
+
+ row = int(r);
+ if (row <= 0 || row >= rowCount-1) { if (!isin) continue; else break; }
+ offset = r - float32(row);
+
+ iVolumeIndex = row * colCount + col;
+ policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - offset) * lengthPerCol)
+
+ iVolumeIndex += colCount;
+ policy_weight(p, iRayIndex, iVolumeIndex, offset * lengthPerCol)
+
+ isin = true;
}
}
@@ -177,5 +235,6 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
} // end loop detector
} // end loop angles
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+ delete pVecProjectionGeometry;
}
-
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.h b/include/astra/ParallelBeamStripKernelProjector2D.h
index 373c7fa..908df1f 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.h
+++ b/include/astra/ParallelBeamStripKernelProjector2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#define _INC_ASTRA_PARALLELBEAMSTROKEKERNELPROJECTOR
#include "ParallelProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
#include "Float32Data2D.h"
#include "Projector2D.h"
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl
index e8e3739..3a21ed6 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -29,7 +29,7 @@ template <typename Policy>
void CParallelBeamStripKernelProjector2D::project(Policy& p)
{
projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
- 0, m_pProjectionGeometry->getDetectorCount(), p);
+ 0, m_pProjectionGeometry->getDetectorCount(), p);
}
template <typename Policy>
@@ -48,250 +48,244 @@ void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int
//----------------------------------------------------------------------------------------
// PROJECT BLOCK
+//
+// Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+//
+// For each angle/detector pair:
+//
+// Let DL=(DLx,DLy) denote the left of the detector (point) in volume coordinates, and
+// Let DR=(DRx,DRy) denote the right of the detector (point) in volume coordinates, and
+// let R=(Rx,Ry) denote the direction of the ray (vector).
+//
+// For mainly vertical rays (|Rx|<=|Ry|),
+// let E=(Ex,Ey) denote the centre of the most upper left pixel:
+// E = (WindowMinX + PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// and let F=(Fx,Fy) denote a vector to the next pixel
+// F = (PixelLengthX, 0)
+//
+// The intersection of the left edge of the strip (DL+aR) with the centre line of the upper row of pixels (E+bF) is
+// { DLx + a*Rx = Ex + b*Fx
+// { DLy + a*Ry = Ey + b*Fy
+// Solving for (a,b) results in:
+// a = (Ey + b*Fy - DLy)/Ry
+// = (Ey - DLy)/Ry
+// b = (DLx + a*Rx - Ex)/Fx
+// = (DLx + (Ey - DLy)*Rx/Ry - Ex)/Fx
+//
+// Define cL as the x-value of the intersection of the left edge of the strip with the upper row in pixel coordinates.
+// cL = b
+//
+// cR, the x-value of the intersection of the right edge of the strip with the upper row in pixel coordinates can be found similarly.
+//
+// The intersection of the ray (DL+aR) with the left line of the second row of pixels (E'+bF) with
+// E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+// expressed in x-value pixel coordinates is
+// cL' = (DLx + (Ey' - DLy)*Rx/Ry - Ex)/Fx.
+// And thus:
+// deltac = cL' - cL = (DLx + (Ey' - DLy)*Rx/Ry - Ex)/Fx - (DLx + (Ey - DLy)*Rx/Ry - Ex)/Fx
+// = [(Ey' - DLy)*Rx/Ry - (Ey - DLy)*Rx/Ry]/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = [Ey' - Ey]*(Rx/Ry)/Fx
+// = -PixelLengthY*(Rx/Ry)/Fx.
+//
+// The projection weight for a certain pixel is defined by the area between two points of
+//
+// _____ LengthPerRow
+// /| | |\
+// / | | | \
+// __/ | | | \__ 0
+// -T -S 0 S T
+// with S = 1/2 - 1/2*|Rx/Ry|, T = 1/2 + 1/2*|Rx/Ry|, and LengthPerRow = pixelLengthX * sqrt(Rx^2+Ry^2) / |Ry|
+//
+// For a certain row, all columns that are 'hit' by this kernel lie in the interval
+// (col_left, col_right) = (floor(cL-1/2+S), floor(cR+3/2-S))
+//
+// The offsets for both is
+// (offsetL, offsetR) = (cL - floor(col_left), cR - floor(col_left))
+//
+// The projection weight is found by the difference between the integrated values of the kernel
+// offset <= -T Kernel = 0
+// -T < offset <= -S Kernel = PixelArea/2*(T+offset)^2/(T-S)
+// -S < offset <= S Kernel = PixelArea/2 + offset
+// S < offset <= T Kernel = PixelArea - PixelArea/2*(T-offset)^2/(T-S)
+// T <= offset: Kernel = PixelArea
+//
template <typename Policy>
void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
{
- ASTRA_ASSERT(m_bIsInitialized);
-
- // Some variables
- float32 theta, t;
- int row, col;
- int iAngle;
- int iDetector;
- float32 res;
- float32 PL, PLimitL, PLimitR;
- float32 xL, xR, XLimitL, XLimitR;
- int x1L,x1R;
- float32 x2L, x2R, updateX;
- int iVolumeIndex, iRayIndex;
-
- float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta;
- float32 fabs_sin_theta, fabs_cos_theta, fabs_inv_sin_theta, fabs_inv_cos_theta;
- float32 PW, PH, DW, inv_PW, inv_PH;
- float32 S, T, U, V, inv_4T;
+ // get vector geometry
+ const CParallelVecProjectionGeometry2D* pVecProjectionGeometry;
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)) {
+ pVecProjectionGeometry = dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry();
+ } else {
+ pVecProjectionGeometry = dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry);
+ }
+
+ // precomputations
+ const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX();
+ const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY();
+ const float32 pixelArea = pixelLengthX * pixelLengthY;
+ const float32 inv_pixelLengthX = 1.0f / pixelLengthX;
+ const float32 inv_pixelLengthY = 1.0f / pixelLengthY;
+ const int colCount = m_pVolumeGeometry->getGridColCount();
+ const int rowCount = m_pVolumeGeometry->getGridRowCount();
+ const int detCount = pVecProjectionGeometry->getDetectorCount();
// loop angles
- for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
- // get values
- theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
- bool switch_t = false;
- if (theta >= 7*PIdiv4) theta -= 2*PI;
- if (theta >= 3*PIdiv4) {
- theta -= PI;
- switch_t = true;
- }
-
- // Precalculate sin, cos, 1/cos
- sin_theta = sin(theta);
- cos_theta = cos(theta);
- inv_cos_theta = 1.0f / cos_theta;
- inv_sin_theta = 1.0f / sin_theta;
-
- fabs_sin_theta = (sin_theta < 0.0f) ? -sin_theta : sin_theta;
- fabs_cos_theta = (cos_theta < 0.0f) ? -cos_theta : cos_theta;
- fabs_inv_cos_theta = (inv_cos_theta < 0.0f) ? -inv_cos_theta : inv_cos_theta;
- fabs_inv_sin_theta = (inv_sin_theta < 0.0f) ? -inv_sin_theta : inv_sin_theta;
-
- // Other precalculations
- PW = m_pVolumeGeometry->getPixelLengthX();
- PH = m_pVolumeGeometry->getPixelLengthY();
- DW = m_pProjectionGeometry->getDetectorWidth();
- inv_PW = 1.0f / PW;
- inv_PH = 1.0f / PH;
-
- // [-45?,45?] and [135?,225?]
- if (theta < PIdiv4) {
-
- // Precalculate kernel limits
- S = -0.5f * fabs_sin_theta * fabs_inv_cos_theta;
- T = -S;
- U = 1.0f + S;
- V = 1.0f - S;
- inv_4T = 0.25f / T;
-
- updateX = sin_theta * inv_cos_theta;
-
- // loop detectors
- for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
-
- iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
- // POLICY: RAY PRIOR
- if (!p.rayPrior(iRayIndex)) continue;
-
- // get t
- t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
- if (switch_t) t = -t;
-
- // calculate left strip extremes (volume coordinates)
- PL = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0) - DW*0.5f) * inv_cos_theta;
- PLimitL = PL - 0.5f * fabs_sin_theta * fabs_inv_cos_theta * PH;
- PLimitR = PLimitL + DW * inv_cos_theta + PH * fabs_sin_theta * fabs_inv_cos_theta;
-
- // calculate strip extremes (pixel coordinates)
- XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
- XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
- xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
- xR = xL + (DW * inv_cos_theta) * inv_PW;
-
- // for each row
- for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-
- // get strip extremes in column indices
- x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f);
- x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f);
-
- // get coords w.r.t leftmost column hit by strip
- x2L = xL - x1L;
- x2R = xR - x1L;
-
- // update strip extremes for the next row
- XLimitL += updateX;
- XLimitR += updateX;
- xL += updateX;
- xR += updateX;
-
- // for each affected col
- for (col = x1L; col <= x1R; ++col) {
-
- if (col < 0 || col >= m_pVolumeGeometry->getGridColCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
- // POLICY: PIXEL PRIOR
- if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- // right
- if (x2R >= V) res = 1.0f;
- else if (x2R > U) res = x2R - (x2R-U)*(x2R-U)*inv_4T;
- else if (x2R >= T) res = x2R;
- else if (x2R > S) res = (x2R-S)*(x2R-S) * inv_4T;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- // left
- if (x2L <= S) {} // - 0.0f
- else if (x2L < T) res -= (x2L-S)*(x2L-S) * inv_4T;
- else if (x2L <= U) res -= x2L;
- else if (x2L < V) res -= x2L - (x2L-U)*(x2L-U)*inv_4T;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- // POLICY: ADD
- p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res);
-
- // POLICY: PIXEL POSTERIOR
- p.pixelPosterior(iVolumeIndex);
-
- x2L -= 1.0f;
- x2R -= 1.0f;
-
- } // end col loop
-
- } // end row loop
-
- // POLICY: RAY POSTERIOR
- p.rayPosterior(iRayIndex);
-
- } // end detector loop
-
- // [45?,135?] and [225?,315?]
- // horizontaly
+ #pragma omp parallel for
+ for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+ // variables
+ float32 DLx, DLy, DRx, DRy, Ex, Ey, S, T, deltac, deltar, offsetL, offsetR, invTminS;
+ float32 res, RxOverRy, RyOverRx, cL, cR, rL, rR;
+ int iVolumeIndex, iRayIndex, iDetector;
+ int row, row_top, row_bottom, col, col_left, col_right;
+
+ const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+ 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);
+ }
- // Precalculate kernel limits
- S = -0.5f * fabs_cos_theta * fabs_inv_sin_theta;
- T = -S;
- U = 1.0f + S;
- V = 1.0f - S;
- inv_4T = 0.25f / T;
-
- updateX = cos_theta * inv_sin_theta;
-
- // loop detectors
- for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
-
- iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
- // POLICY: RAY PRIOR
- if (!p.rayPrior(iRayIndex)) continue;
+ Ex = m_pVolumeGeometry->getWindowMinY() + pixelLengthX*0.5f;
+ Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
- // get t
- t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
- if (switch_t) t = -t;
+ // loop detectors
+ for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
- // calculate left strip extremes (volume coordinates)
- PL = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0) + DW*0.5f) * inv_sin_theta;
- PLimitL = PL + 0.5f * fabs_cos_theta * fabs_inv_sin_theta * PW;
- PLimitR = PLimitL - DW * inv_sin_theta - PH * fabs_cos_theta * fabs_inv_sin_theta;
-
- // calculate strip extremes (pixel coordinates)
- XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH;
- XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH;
- xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH;
- xR = xL + (DW * fabs_inv_sin_theta) * inv_PH;
-
- // for each col
- for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
- // get strip extremes in column indices
- x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f);
- x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f);
-
- // get coords w.r.t leftmost column hit by strip
- x2L = xL - x1L;
- x2R = xR - x1L;
-
- // update strip extremes for the next row
- XLimitL += updateX;
- XLimitR += updateX;
- xL += updateX;
- xR += updateX;
+ iRayIndex = iAngle * detCount + iDetector;
- // for each affected col
- for (row = x1L; row <= x1R; ++row) {
-
- if (row < 0 || row >= m_pVolumeGeometry->getGridRowCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
-
- // POLICY: PIXEL PRIOR
- if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- // right
- if (x2R >= V) res = 1.0f;
- else if (x2R > U) res = x2R - (x2R-U)*(x2R-U)*inv_4T;
- else if (x2R >= T) res = x2R;
- else if (x2R > S) res = (x2R-S)*(x2R-S) * inv_4T;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- // left
- if (x2L <= S) {} // - 0.0f
- else if (x2L < T) res -= (x2L-S)*(x2L-S) * inv_4T;
- else if (x2L <= U) res -= x2L;
- else if (x2L < V) res -= x2L - (x2L-U)*(x2L-U)*inv_4T;
- else { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
- // POLICY: ADD
- p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res);
-
- // POLICY: PIXEL POSTERIOR
- p.pixelPosterior(iVolumeIndex);
-
- x2L -= 1.0f;
- x2R -= 1.0f;
-
- } // end row loop
-
- } // end col loop
-
- // POLICY: RAY POSTERIOR
- p.rayPosterior(iRayIndex);
-
- } // end detector loop
-
-
- } // end theta switch
+ // POLICY: RAY PRIOR
+ if (!p.rayPrior(iRayIndex)) continue;
+
+ DLx = proj->fDetSX + iDetector * proj->fDetUX;
+ DLy = proj->fDetSY + iDetector * proj->fDetUY;
+ DRx = DLx + proj->fDetUX;
+ DRy = DLy + proj->fDetUY;
+
+ // vertically
+ if (vertical) {
+
+ // calculate cL and cR for row 0
+ cL = (DLx + (Ey - DLy)*RxOverRy - Ex) * inv_pixelLengthX;
+ cR = (DRx + (Ey - DRy)*RxOverRy - Ex) * inv_pixelLengthX;
+
+ if (cR < cL) {
+ float32 tmp = cL;
+ cL = cR;
+ cR = tmp;
+ }
+
+ // loop rows
+ for (row = 0; row < rowCount; ++row, cL += deltac, cR += deltac) {
+
+ col_left = int(cL-0.5f+S);
+ col_right = int(cR+1.5-S);
+
+ if (col_left < 0) col_left = 0;
+ if (col_right > colCount-1) col_right = colCount-1;
+
+ float32 tmp = float32(col_left);
+ offsetL = cL - tmp;
+ offsetR = cR - tmp;
+
+ // loop columns
+ for (col = col_left; col <= col_right; ++col, offsetL -= 1.0f, offsetR -= 1.0f) {
+
+ iVolumeIndex = row * colCount + col;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+
+ // right ray edge
+ if (T <= offsetR) res = 1.0f;
+ else if (S < offsetR) res = 1.0f - 0.5f*(T-offsetR)*(T-offsetR)*invTminS;
+ else if (-S < offsetR) res = 0.5f + offsetR;
+ else if (-T < offsetR) res = 0.5f*(offsetR+T)*(offsetR+T)*invTminS;
+ else res = 0.0f;
+
+ // left ray edge
+ if (T <= offsetL) res -= 1.0f;
+ else if (S < offsetL) res -= 1.0f - 0.5f*(T-offsetL)*(T-offsetL)*invTminS;
+ 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.pixelPosterior(iVolumeIndex);
+ }
+ }
+ }
+ }
+
+ // horizontally
+ else {
+
+ // calculate rL and rR for row 0
+ rL = -(DLy + (Ex - DLx)*RyOverRx - Ey) * inv_pixelLengthY;
+ rR = -(DRy + (Ex - DRx)*RyOverRx - Ey) * inv_pixelLengthY;
+
+ if (rR < rL) {
+ float32 tmp = rL;
+ rL = rR;
+ rR = tmp;
+ }
+
+ // loop columns
+ for (col = 0; col < colCount; ++col, rL += deltar, rR += deltar) {
+
+ row_top = int(rL-0.5f+S);
+ row_bottom = int(rR+1.5-S);
+
+ if (row_top < 0) row_top = 0;
+ if (row_bottom > rowCount-1) row_bottom = rowCount-1;
+
+ float32 tmp = float32(row_top);
+ offsetL = rL - tmp;
+ offsetR = rR - tmp;
+
+ // loop rows
+ for (row = row_top; row <= row_bottom; ++row, offsetL -= 1.0f, offsetR -= 1.0f) {
+
+ iVolumeIndex = row * colCount + col;
+ // POLICY: PIXEL PRIOR + ADD + POSTERIOR
+ if (p.pixelPrior(iVolumeIndex)) {
+
+ // right ray edge
+ if (T <= offsetR) res = 1.0f;
+ else if (S < offsetR) res = 1.0f - 0.5f*(T-offsetR)*(T-offsetR)*invTminS;
+ else if (-S < offsetR) res = 0.5f + offsetR;
+ else if (-T < offsetR) res = 0.5f*(offsetR+T)*(offsetR+T)*invTminS;
+ else res = 0.0f;
+
+ // left ray edge
+ if (T <= offsetL) res -= 1.0f;
+ else if (S < offsetL) res -= 1.0f - 0.5f*(T-offsetL)*(T-offsetL)*invTminS;
+ 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.pixelPosterior(iVolumeIndex);
+ }
+ }
+ }
+ }
+
+ // POLICY: RAY POSTERIOR
+ p.rayPosterior(iRayIndex);
+
+ } // end loop detector
+ } // end loop angles
- } // end angle loop
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+ delete pVecProjectionGeometry;
}
-
-
diff --git a/include/astra/ParallelProjectionGeometry2D.h b/include/astra/ParallelProjectionGeometry2D.h
index 963bb06..9625d35 100644
--- a/include/astra/ParallelProjectionGeometry2D.h
+++ b/include/astra/ParallelProjectionGeometry2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#define _INC_ASTRA_PARALLELPROJECTIONGEOMETRY2D
#include "ProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
namespace astra
{
@@ -82,8 +83,7 @@ public:
CParallelProjectionGeometry2D(int _iProjectionAngleCount,
int _iDetectorCount,
float32 _fDetectorWidth,
- const float32* _pfProjectionAngles,
- const float32* _pfExtraDetectorOffsets = 0);
+ const float32* _pfProjectionAngles);
/** Copy constructor.
*/
@@ -115,8 +115,7 @@ public:
bool initialize(int _iProjectionAngleCount,
int _iDetectorCount,
float32 _fDetectorWidth,
- const float32* _pfProjectionAngles,
- const float32* _pfExtraDetectorOffsets = 0);
+ const float32* _pfProjectionAngles);
/** Create a hard copy.
*/
@@ -133,7 +132,7 @@ public:
* @param _sType geometry type to compare to.
* @return true if _sType == "parallel".
*/
- virtual bool isOfType(const std::string& _sType);
+ virtual bool isOfType(const std::string& _sType);
/** Get all settings in a Config object.
*
@@ -150,7 +149,12 @@ public:
*
* @return a unit vector describing the direction
*/
- virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex = 0);
+ virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex = 0);
+
+ /** Create a vector geom
+ */
+ CParallelVecProjectionGeometry2D* toVectorGeometry();
+
};
} // namespace astra
diff --git a/include/astra/ParallelVecProjectionGeometry2D.h b/include/astra/ParallelVecProjectionGeometry2D.h
new file mode 100644
index 0000000..96f8a54
--- /dev/null
+++ b/include/astra/ParallelVecProjectionGeometry2D.h
@@ -0,0 +1,163 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+ 2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+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 <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#ifndef _INC_ASTRA_PARALLELVECPROJECTIONGEOMETRY2D
+#define _INC_ASTRA_PARALLELVECPROJECTIONGEOMETRY2D
+
+#include "ProjectionGeometry2D.h"
+#include "GeometryUtil2D.h"
+
+namespace astra
+{
+
+/**
+ * This class defines a 2D parallel beam projection geometry.
+ *
+ * \par XML Configuration
+ * \astra_xml_item{DetectorCount, int, Number of detectors for each projection.}
+ *
+ * \par MATLAB example
+ * \astra_code{
+ * proj_geom = astra_struct('parallel_vec');\n
+ * proj_geom.DetectorCount = 512;\n
+ * proj_geom.Vectors = V;\n
+ * }
+ *
+ * \par Vectors
+ * Vectors is a matrix containing the actual geometry. Each row corresponds
+ * to a single projection, and consists of:
+ * ( rayX, rayY, dX, dY, uX, uY)
+ * ray: the ray direction
+ * d : the centre of the detector line
+ * u : the vector from detector pixel (0) to (1)
+ */
+class _AstraExport CParallelVecProjectionGeometry2D : public CProjectionGeometry2D
+{
+protected:
+
+ SParProjection *m_pProjectionAngles;
+
+public:
+
+ /** Default constructor. Sets all variables to zero. Note that this constructor leaves the object in an unusable state and must
+ * be followed by a call to init().
+ */
+ CParallelVecProjectionGeometry2D();
+
+ /** Constructor.
+ *
+ * @param _iProjectionAngleCount Number of projection angles.
+ * @param _iDetectorCount Number of detectors, i.e., the number of detector measurements for each projection angle.
+ * @param _pfProjectionAngles Pointer to an array of projection angles. The angles will be copied from this array.
+ */
+ CParallelVecProjectionGeometry2D(int _iProjectionAngleCount,
+ int _iDetectorCount,
+ const SParProjection* _pfProjectionAngles);
+
+ /** Copy constructor.
+ */
+ CParallelVecProjectionGeometry2D(const CParallelVecProjectionGeometry2D& _projGeom);
+
+ /** Assignment operator.
+ */
+ CParallelVecProjectionGeometry2D& operator=(const CParallelVecProjectionGeometry2D& _other);
+
+ /** Destructor.
+ */
+ virtual ~CParallelVecProjectionGeometry2D();
+
+ /** Initialize the geometry with a config object.
+ *
+ * @param _cfg Configuration Object
+ * @return initialization successful?
+ */
+ virtual bool initialize(const Config& _cfg);
+
+ /** Initialization. This function MUST be called after using the default constructor and MAY be called to
+ * reset a previously initialized object.
+ *
+ * @param _iProjectionAngleCount Number of projection angles.
+ * @param _iDetectorCount Number of detectors, i.e., the number of detector measurements for each projection angle.
+ * @param _pfProjectionAngles Pointer to an array of projection angles. The angles will be copied from this array.
+ */
+ bool initialize(int _iProjectionAngleCount,
+ int _iDetectorCount,
+ const SParProjection* _pfProjectionAngles);
+
+ virtual bool _check();
+
+ /** Create a hard copy.
+ */
+ virtual CProjectionGeometry2D* clone();
+
+ /** Returns true if the type of geometry defined in this class is the one specified in _sType.
+ *
+ * @param _sType geometry type to compare to.
+ * @return true if _sType == "fanflat_vec".
+ */
+ virtual bool isOfType(const std::string& _sType);
+
+ /** Return true if this geometry instance is the same as the one specified.
+ *
+ * @return true if this geometry instance is the same as the one specified.
+ */
+ virtual bool isEqual(CProjectionGeometry2D*) const;
+
+ /** Get all settings in a Config object.
+ *
+ * @return Configuration Object.
+ */
+ virtual Config* getConfiguration() const;
+
+
+ /** Get the value for t and theta, based upon the row and column index.
+ *
+ * @param _iRow row index
+ * @param _iColumn column index
+ * @param _fT output: value of t
+ * @param _fTheta output: value of theta, always lies within the [0,pi[ interval.
+ */
+ virtual void getRayParams(int _iRow, int _iColumn, float32& _fT, float32& _fTheta) const;
+
+ /**
+ * Returns a vector describing the direction of a ray belonging to a certain detector
+ *
+ * @param _iProjectionIndex index of projection
+ * @param _iProjectionIndex index of detector
+ *
+ * @return a unit vector describing the direction
+ */
+ virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex);
+
+ const SParProjection* getProjectionVectors() const { return m_pProjectionAngles; }
+
+};
+
+} // namespace astra
+
+#endif /* _INC_ASTRA_PARALLELVECPROJECTIONGEOMETRY2D */
diff --git a/include/astra/ProjectionGeometry2D.h b/include/astra/ProjectionGeometry2D.h
index 393db77..504f588 100644
--- a/include/astra/ProjectionGeometry2D.h
+++ b/include/astra/ProjectionGeometry2D.h
@@ -64,10 +64,6 @@ protected:
*/
float32 m_fDetectorWidth;
- /** An array of m_iProjectionAngleCount elements containing an extra detector offset for each projection.
- */
- float32* m_pfExtraDetectorOffset;
-
/** Dynamically allocated array of projection angles. All angles are represented in radians and lie in
* the [0,2pi[ interval.
*/
@@ -93,8 +89,7 @@ protected:
CProjectionGeometry2D(int _iProjectionAngleCount,
int _iDetectorCount,
float32 _fDetectorWidth,
- const float32* _pfProjectionAngles,
- const float32* _pfExtraDetectorOffsets = 0);
+ const float32* _pfProjectionAngles);
/** Copy constructor.
*/
@@ -120,8 +115,7 @@ protected:
bool _initialize(int _iProjectionAngleCount,
int _iDetectorCount,
float32 _fDetectorWidth,
- const float32* _pfProjectionAngles,
- const float32* _pfExtraDetectorOffsets = 0);
+ const float32* _pfProjectionAngles);
public:
@@ -201,9 +195,6 @@ public:
*/
float32 getProjectionAngleDegrees(int _iProjectionIndex) const;
- float32 getExtraDetectorOffset(int iAngle) const;
- const float32* getExtraDetectorOffset() const { return m_pfExtraDetectorOffset; }
-
/** Get the index coordinate of a point on a detector array.
*
* @param _fOffset distance between the center of the detector array and a certain point
@@ -272,12 +263,6 @@ public:
//----------------------------------------------------------------------------------------
-inline float32 CProjectionGeometry2D::getExtraDetectorOffset(int _iAngle) const
-{
- return m_pfExtraDetectorOffset ? m_pfExtraDetectorOffset[_iAngle] : 0.0f;
-}
-
-
// Get the initialization state.
inline bool CProjectionGeometry2D::isInitialized() const
{
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 11:26:45 +0000