summaryrefslogtreecommitdiffstats
path: root/include
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2017-10-12 15:17:10 +0200
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2017-10-12 15:17:10 +0200
commit08df6c4167118e0a3b5f128078f9c13f206a171a (patch)
tree0883b5967008e83fb1ae3f0a00bd675bba7116b4 /include
parent5ad3d86b6e2c39de7465186ec8702053a82b6152 (diff)
downloadastra-08df6c4167118e0a3b5f128078f9c13f206a171a.tar.gz
astra-08df6c4167118e0a3b5f128078f9c13f206a171a.tar.bz2
astra-08df6c4167118e0a3b5f128078f9c13f206a171a.tar.xz
astra-08df6c4167118e0a3b5f128078f9c13f206a171a.zip
Fix some warnings
Diffstat (limited to 'include')
-rw-r--r--include/astra/FanFlatBeamLineKernelProjector2D.inl22
-rw-r--r--include/astra/ParallelBeamBlobKernelProjector2D.inl1
-rw-r--r--include/astra/ParallelBeamLineKernelProjector2D.inl179
-rw-r--r--include/astra/ParallelBeamLinearKernelProjector2D.inl151
-rw-r--r--include/astra/ParallelBeamStripKernelProjector2D.inl126
5 files changed, 236 insertions, 243 deletions
diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl
index c1e1e94..927aa09 100644
--- a/include/astra/FanFlatBeamLineKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl
@@ -99,6 +99,9 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in
Ry = proj->fSrcY - Dy;
bool vertical = fabs(Rx) < fabs(Ry);
+ bool isin = false;
+
+ // vertically
if (vertical) {
RxOverRy = Rx/Ry;
lengthPerRow = detSize * pixelLengthX * sqrt(Rx*Rx + Ry*Ry) / abs(Ry);
@@ -106,19 +109,6 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in
S = 0.5f - 0.5f*fabs(RxOverRy);
T = 0.5f + 0.5f*fabs(RxOverRy);
invTminSTimesLengthPerRow = lengthPerRow / (T - S);
- } else {
- RyOverRx = Ry/Rx;
- lengthPerCol = detSize * 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);
- }
-
- bool isin = false;
-
- // vertically
- if (vertical) {
// calculate c for row 0
c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
@@ -163,6 +153,12 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in
// horizontally
else {
+ RyOverRx = Ry/Rx;
+ lengthPerCol = detSize * 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);
// calculate r for col 0
r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.inl b/include/astra/ParallelBeamBlobKernelProjector2D.inl
index 67662ad..ccd2166 100644
--- a/include/astra/ParallelBeamBlobKernelProjector2D.inl
+++ b/include/astra/ParallelBeamBlobKernelProjector2D.inl
@@ -124,7 +124,6 @@ void CParallelBeamBlobKernelProjector2D::projectBlock_internal(int _iProjFrom, i
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
for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl
index e516fe1..7db0a34 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLineKernelProjector2D.inl
@@ -49,94 +49,94 @@ void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int
//----------------------------------------------------------------------------------------
-// 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
-//
+/* 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)
{
@@ -155,7 +155,6 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
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 (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
@@ -163,7 +162,7 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
// 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;
+ int iVolumeIndex, iRayIndex, row, col;
const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl
index d2c529f..61e4973 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -51,80 +51,80 @@ void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, in
//----------------------------------------------------------------------------------------
-// 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
-//
+/* 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)
{
@@ -143,13 +143,12 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
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 (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
// variables
- float32 Dx, Dy, Ex, Ey, x, y, c, r, deltac, deltar, offset;
+ float32 Dx, Dy, Ex, Ey, c, r, deltac, deltar, offset;
float32 RxOverRy, RyOverRx, lengthPerRow, lengthPerCol;
int iVolumeIndex, iRayIndex, row, col, iDetector;
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl
index 4f828f0..fdfcd90 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -47,69 +47,69 @@ 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
-//
+/* 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)
{