diff options
| -rw-r--r-- | include/astra/FanFlatBeamLineKernelProjector2D.inl | 22 | ||||
| -rw-r--r-- | include/astra/ParallelBeamBlobKernelProjector2D.inl | 1 | ||||
| -rw-r--r-- | include/astra/ParallelBeamLineKernelProjector2D.inl | 179 | ||||
| -rw-r--r-- | include/astra/ParallelBeamLinearKernelProjector2D.inl | 151 | ||||
| -rw-r--r-- | include/astra/ParallelBeamStripKernelProjector2D.inl | 126 | ||||
| -rw-r--r-- | src/CudaFilteredBackProjectionAlgorithm.cpp | 13 |
6 files changed, 236 insertions, 256 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) { diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 2829b7d..f3cca12 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -64,7 +64,6 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) // if already initialized, clear first if (m_bIsInitialized) { -#warning FIXME Necessary? clear(); } @@ -236,18 +235,6 @@ bool CCudaFilteredBackProjectionAlgorithm::check() return true; } -static int calcNextPowerOfTwo(int _iValue) -{ - int iOutput = 1; - - while(iOutput < _iValue) - { - iOutput *= 2; - } - - return iOutput; -} - static bool stringCompareLowerCase(const char * _stringA, const char * _stringB) { int iCmpReturn = 0; |