diff options
39 files changed, 1712 insertions, 1074 deletions
diff --git a/astra_vc11.vcxproj b/astra_vc11.vcxproj index bc11b23..c8d0ec3 100644 --- a/astra_vc11.vcxproj +++ b/astra_vc11.vcxproj @@ -537,6 +537,7 @@ <ClCompile Include="src\ParallelBeamStripKernelProjector2D.cpp" /> <ClCompile Include="src\ParallelProjectionGeometry2D.cpp" /> <ClCompile Include="src\ParallelProjectionGeometry3D.cpp" /> + <ClCompile Include="src\ParallelVecProjectionGeometry2D.cpp" /> <ClCompile Include="src\ParallelVecProjectionGeometry3D.cpp" /> <ClCompile Include="src\PlatformDepSystemCode.cpp" /> <ClCompile Include="src\ProjectionGeometry2D.cpp" /> @@ -648,6 +649,7 @@ <ClInclude Include="include\astra\ParallelProjectionGeometry2D.h" /> <ClInclude Include="include\astra\ParallelProjectionGeometry3D.h" /> <ClInclude Include="include\astra\ParallelVecProjectionGeometry3D.h" /> + <ClInclude Include="include\astra\ParallelVecProjectionGeometry2D.h" /> <ClInclude Include="include\astra\PlatformDepSystemCode.h" /> <ClInclude Include="include\astra\ProjectionGeometry2D.h" /> <ClInclude Include="include\astra\ProjectionGeometry3D.h" /> diff --git a/astra_vc11.vcxproj.filters b/astra_vc11.vcxproj.filters index a597962..d0e7254 100644 --- a/astra_vc11.vcxproj.filters +++ b/astra_vc11.vcxproj.filters @@ -312,6 +312,10 @@ <ClCompile Include="src\CudaSirtAlgorithm3D.cpp"> <Filter>CUDA\astra source</Filter> </ClCompile> + <ClCompile Include="src\GeometryUtil3D.cpp" /> + <ClCompile Include="src\ParallelVecProjectionGeometry2D.cpp"> + <Filter>Geometries\source</Filter> + </ClCompile> </ItemGroup> <ItemGroup> <ClInclude Include="include\astra\Algorithm.h"> @@ -653,6 +657,11 @@ <ClInclude Include="cuda\3d\util3d.h"> <Filter>CUDA\cuda headers</Filter> </ClInclude> + <ClInclude Include="include\astra\GeometryUtil2D.h" /> + <ClInclude Include="include\astra\GeometryUtil3D.h" /> + <ClInclude Include="include\astra\ParallelVecProjectionGeometry2D.h"> + <Filter>Geometries\headers</Filter> + </ClInclude> </ItemGroup> <ItemGroup> <None Include="include\astra\DataProjectorPolicies.inl"> diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 2d862f2..bc9821a 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -129,6 +129,7 @@ BASE_OBJECTS=\ src/ParallelBeamLineKernelProjector2D.lo \ src/ParallelBeamStripKernelProjector2D.lo \ src/ParallelProjectionGeometry2D.lo \ + src/ParallelVecProjectionGeometry2D.lo \ src/ParallelProjectionGeometry3D.lo \ src/ParallelVecProjectionGeometry3D.lo \ src/PlatformDepSystemCode.lo \ diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index 2f72db0..379c690 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -828,17 +828,12 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, // If there are existing detector offsets, or if we need to translate, // we need to return offsets - if (pProjGeom->getExtraDetectorOffset() || offCenter) + if (offCenter) { float* offset = new float[nth]; - if (pProjGeom->getExtraDetectorOffset()) { - for (int i = 0; i < nth; ++i) - offset[i] = pProjGeom->getExtraDetectorOffset()[i]; - } else { - for (int i = 0; i < nth; ++i) - offset[i] = 0.0f; - } + for (int i = 0; i < nth; ++i) + offset[i] = 0.0f; if (offCenter) { float dx = (pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl index 2f87659..1676cb1 100644 --- a/include/astra/FanFlatBeamLineKernelProjector2D.inl +++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl @@ -51,246 +51,201 @@ 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; - - const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry); - const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry); - - float32 old_theta, theta, alpha; + float32 detX, detY, S, T, I, x, y, c, r, update_c, update_r, offset; + float32 lengthPerRow, lengthPerCol, inv_pixelLengthX, inv_pixelLengthY, invTminSTimesLengthPerRow, invTminSTimesLengthPerCol; + int iVolumeIndex, iRayIndex, row, col, iAngle, iDetector, colCount, rowCount, detCount; const SFanProjection * proj = 0; + // 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 + inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); + inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); + colCount = m_pVolumeGeometry->getGridColCount(); + rowCount = m_pVolumeGeometry->getGridRowCount(); + detCount = pVecProjectionGeometry->getDetectorCount(); + // 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; - } + 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); + + detX = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX; + detY = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY; + + float32 fRayX = proj->fSrcX - detX; + float32 fRayY = proj->fSrcY - detY; + + bool vertical = fabs(fRayX) < fabs(fRayY); + if (vertical) { + lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * sqrt(fRayY*fRayY + fRayX*fRayX) / abs(fRayY); + update_c = -m_pVolumeGeometry->getPixelLengthY() * (fRayX/fRayY) * inv_pixelLengthX; + S = 0.5f - 0.5f*fabs(fRayX/fRayY); + T = 0.5f + 0.5f*fabs(fRayX/fRayY); + invTminSTimesLengthPerRow = lengthPerRow / (T - S); } else { - assert(false); + lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(fRayY*fRayY + fRayX*fRayX) / abs(fRayX); + update_r = -m_pVolumeGeometry->getPixelLengthX() * (fRayY/fRayX) * inv_pixelLengthY; + S = 0.5f - 0.5f*fabs(fRayY/fRayX); + T = 0.5f + 0.5f*fabs(fRayY/fRayX); + 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); - // vertically - if (old_theta <= PIdiv4) { - + if (vertical) { + // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX; + x = detX + (fRayX/fRayY)*(m_pVolumeGeometry->pixelRowToCenterY(0)-detY); + c = (x - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f; // 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 += update_c) { + + col = int(c+0.5f); + offset = c - float32(col); - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue; + if (col <= 0 || col >= colCount-1) continue; // 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 (offset < -S) { + I = (offset + T) * invTminSTimesLengthPerRow; + + iVolumeIndex = row * colCount + col - 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); - } + iVolumeIndex++; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, I); + 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 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); - } + else if (S < offset) { + I = (offset - S) * invTminSTimesLengthPerRow; + + iVolumeIndex = row * colCount + col; + // 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, 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); - } + } + + // centre + else { + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow); + p.pixelPosterior(iVolumeIndex); } } + } } // 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 y for col 0 + y = detY + (fRayY/fRayX)*(m_pVolumeGeometry->pixelColToCenterX(0)-detX); + r = (m_pVolumeGeometry->getWindowMaxY() - y) * inv_pixelLengthY - 0.5f; // 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 += update_r) { + + int row = int(r+0.5f); + offset = r - float32(row); - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue; + if (row <= 0 || row >= rowCount-1) continue; // 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 (offset < -S) { + I = (offset + T) * invTminSTimesLengthPerCol; + + iVolumeIndex = (row-1) * colCount + 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); - } + iVolumeIndex += colCount; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, I); + p.pixelPosterior(iVolumeIndex); } } - // 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 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); - } + else if (S < offset) { + I = (offset - S) * invTminSTimesLengthPerCol; + + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-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); - } + + iVolumeIndex += colCount; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, I); + p.pixelPosterior(iVolumeIndex); + } + } + + // centre + else { + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol); + p.pixelPosterior(iVolumeIndex); } } + } - } // 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 180fe68..8750232 100644 --- a/include/astra/FanFlatProjectionGeometry2D.h +++ b/include/astra/FanFlatProjectionGeometry2D.h @@ -30,6 +30,7 @@ $Id$ #define _INC_ASTRA_FANFLATPROJECTIONGEOMETRY2D #include "ProjectionGeometry2D.h" +#include "FanFlatVecProjectionGeometry2D.h" #include <cmath> @@ -190,6 +191,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 d4ee92e..680cecd 100644 --- a/include/astra/GeometryUtil2D.h +++ b/include/astra/GeometryUtil2D.h @@ -31,6 +31,18 @@ $Id$ 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; +}; + + struct SFanProjection { // the source float fSrcX, fSrcY; diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.h b/include/astra/ParallelBeamBlobKernelProjector2D.h index ecf71ef..a718f56 100644 --- a/include/astra/ParallelBeamBlobKernelProjector2D.h +++ b/include/astra/ParallelBeamBlobKernelProjector2D.h @@ -215,7 +215,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 467e066..1945cdd 100644 --- a/include/astra/ParallelBeamBlobKernelProjector2D.inl +++ b/include/astra/ParallelBeamBlobKernelProjector2D.inl @@ -27,186 +27,151 @@ $Id$ */ - -//---------------------------------------------------------------------------------------- -// 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); - - int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) return; - - // get values - float32 t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector); - float32 theta = m_pProjectionGeometry->getProjectionAngle(_iProjection); - if (theta >= 7*PIdiv4) theta -= 2*PI; - - bool flip = false; + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} - if (theta >= 3*PIdiv4) { - theta -= PI; - t = -t; - flip = true; +//---------------------------------------------------------------------------------------- +// PROJECT BLOCK - vector projection geometry +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); } + // precomputations + 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 detX, detY, x, y, c, r, update_c, update_r, offset, invBlobExtent, weight, d; + 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) { + update_c = -m_pVolumeGeometry->getPixelLengthY() * (proj->fRayX/proj->fRayY) * inv_pixelLengthX; + float32 normR = sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX); + invBlobExtent = m_pVolumeGeometry->getPixelLengthY() / abs(m_fBlobSize * normR / proj->fRayY); + } else { + update_r = -m_pVolumeGeometry->getPixelLengthX() * (proj->fRayY/proj->fRayX) * inv_pixelLengthY; + float32 normR = sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX); + invBlobExtent = m_pVolumeGeometry->getPixelLengthX() / abs(m_fBlobSize * normR / proj->fRayX); + } - if (theta <= PIdiv4) { // -pi/4 <= theta <= pi/4 - - // precalculate sin, cos, 1/cos - float32 sin_theta = sin(theta); - float32 cos_theta = cos(theta); - float32 inv_cos_theta = 1.0f / cos_theta; - - // 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; - - // 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; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { + // loop detectors + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { - // 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)]; - } + iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; - int iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, fWeight); - p.pixelPosterior(iVolumeIndex); + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + + detX = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX; + detY = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY; + + // vertically + if (vertical) { + + // calculate x for row 0 + x = detX + (proj->fRayX/proj->fRayY)*(m_pVolumeGeometry->pixelRowToCenterY(0)-detY); + c = (x - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f; + + // for each row + for (row = 0; row < rowCount; ++row, c += update_c) { + + col_left = int(c - 0.5f - m_fBlobSize); + col_right = int(c + 0.5f + m_fBlobSize); + + if (col_left < 0) col_left = 0; + if (col_right > colCount-1) col_right = colCount-1; + + // for each column + for (col = col_left; col <= col_right; ++col) { + + offset = abs(c - float32(col)) * invBlobExtent; + index = (int)(offset*m_iBlobSampleCount+0.5f); + weight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)]; + + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, weight); + 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)]; - } + // calculate y for col 0 + y = detY + (proj->fRayY/proj->fRayX)*(m_pVolumeGeometry->pixelColToCenterX(0)-detX); + r = (m_pVolumeGeometry->getWindowMaxY() - y) * inv_pixelLengthY - 0.5f; + // for each col + for (col = 0; col < colCount; ++col, r += update_r) { - int iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, fWeight); - p.pixelPosterior(iVolumeIndex); - } - } - } + row_top = int(r - 0.5f - m_fBlobSize); + row_bottom = int(r + 0.5f + m_fBlobSize); - // update P and x - P += updatePerCol; - x += updatePerCol * inv_pixelLengthY; - } + if (row_top < 0) row_top = 0; + if (row_bottom > rowCount-1) row_bottom = rowCount-1; - } - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); + // for each column + for (row = row_top; row <= row_bottom; ++row) { + offset = abs(r - float32(row)) * invBlobExtent; + index = (int)(offset*m_iBlobSampleCount+0.5f); + weight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)]; + + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, weight); + 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 24f43d5..5a10021 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.h +++ b/include/astra/ParallelBeamLineKernelProjector2D.h @@ -30,6 +30,7 @@ $Id$ #define _INC_ASTRA_PARALLELBEAMLINEKERNELPROJECTOR #include "ParallelProjectionGeometry2D.h" +#include "ParallelVecProjectionGeometry2D.h" #include "Float32Data2D.h" #include "Projector2D.h" @@ -180,6 +181,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 199d69e..9c58d60 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.inl +++ b/include/astra/ParallelBeamLineKernelProjector2D.inl @@ -30,7 +30,7 @@ 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,49 +48,54 @@ void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int } + + //---------------------------------------------------------------------------------------- -// PROJECT BLOCK +// PROJECT BLOCK - vector projection geometry 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 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 - 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 detX, detY, S, T, I, x, y, c, r, update_c, update_r, offset; + float32 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) { + lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); + update_c = -m_pVolumeGeometry->getPixelLengthY() * (proj->fRayX/proj->fRayY) * inv_pixelLengthX; + S = 0.5f - 0.5f*fabs(proj->fRayX/proj->fRayY); + T = 0.5f + 0.5f*fabs(proj->fRayX/proj->fRayY); + invTminSTimesLengthPerRow = lengthPerRow / (T - S); + } else { + lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); + update_r = -m_pVolumeGeometry->getPixelLengthX() * (proj->fRayY/proj->fRayX) * inv_pixelLengthY; + S = 0.5f - 0.5f*fabs(proj->fRayY/proj->fRayX); + T = 0.5f + 0.5f*fabs(proj->fRayY/proj->fRayX); + 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); - // loop detectors for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { @@ -99,185 +104,140 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i // POLICY: RAY PRIOR if (!p.rayPrior(iRayIndex)) continue; - // get t - t = m_pProjectionGeometry->indexToDetectorOffset(iDetector); - if (switch_t) t = -t; + detX = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX; + detY = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY; // 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; + if (vertical) { - // get coords - int nextx1 = int((x > 0.0f) ? x : x-1.0f); - float nextx2 = x - nextx1; + // calculate x for row 0 + x = detX + (proj->fRayX/proj->fRayY)*(m_pVolumeGeometry->pixelRowToCenterY(0)-detY); + c = (x - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f; // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { + for (row = 0; row < rowCount; ++row, c += update_c) { - x1 = nextx1; - x2 = nextx2; - - nextx2 += updatePerRow; - while (nextx2 >= 1.0f) { - nextx2 -= 1.0f; - nextx1++; - } - while (nextx2 < 0.0f) { - nextx2 += 1.0f; - nextx1--; - } + col = int(c+0.5f); + offset = c - float32(col); - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue; + if (col <= 0 || col >= colCount-1) continue; // 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) { + I = (offset + T) * invTminSTimesLengthPerRow; + + iVolumeIndex = row * colCount + col - 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); - } + iVolumeIndex++; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, I); + 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) { + I = (offset - S) * invTminSTimesLengthPerRow; + + iVolumeIndex = row * colCount + col; + // 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, I); + 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); - } + } + + // centre + else { + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow); + p.pixelPosterior(iVolumeIndex); } } + } } // 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->getWindowMaxY() - P) * inv_pixelLengthY; + else { - // get coords - int nextx1 = int((x > 0.0f) ? x : x-1.0f); - float nextx2 = x - nextx1; + // calculate y for col 0 + y = detY + (proj->fRayY/proj->fRayX)*(m_pVolumeGeometry->pixelColToCenterX(0)-detX); + r = (m_pVolumeGeometry->getWindowMaxY() - y) * inv_pixelLengthY - 0.5f; // for each col - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { + for (col = 0; col < colCount; ++col, r += update_r) { - x1 = nextx1; - x2 = nextx2; + int row = int(r+0.5f); + offset = r - float32(row); - nextx2 += updatePerCol; - while (nextx2 >= 1.0f) { - nextx2 -= 1.0f; - nextx1++; - } - while (nextx2 < 0.0f) { - nextx2 += 1.0f; - nextx1--; - } - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue; + if (row <= 0 || row >= rowCount-1) continue; // 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); - } + if (offset < -S) { + I = (offset + T) * invTminSTimesLengthPerCol; + + iVolumeIndex = (row-1) * colCount + 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); - } + iVolumeIndex += colCount; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, I); + p.pixelPosterior(iVolumeIndex); } } - // 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->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } + else if (S < offset) { + I = (offset - S) * invTminSTimesLengthPerCol; + + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); + p.pixelPosterior(iVolumeIndex); + } + + iVolumeIndex += colCount; + // 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()) { - 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: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol); + p.pixelPosterior(iVolumeIndex); } } + } - } // end loop col + } // POLICY: RAY POSTERIOR p.rayPosterior(iRayIndex); @@ -286,4 +246,3 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i } // end loop angles } - diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.h b/include/astra/ParallelBeamLinearKernelProjector2D.h index 855093a..8e25c38 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.h +++ b/include/astra/ParallelBeamLinearKernelProjector2D.h @@ -30,6 +30,7 @@ $Id$ #define _INC_ASTRA_PARALLELLINEARKERNELPROJECTOR #include "ParallelProjectionGeometry2D.h" +#include "ParallelVecProjectionGeometry2D.h" #include "Float32Data2D.h" #include "Projector2D.h" @@ -185,7 +186,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 67e0d58..79b82d4 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.inl +++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl @@ -26,12 +26,11 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. $Id$ */ - 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> @@ -48,48 +47,48 @@ void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, in _iDetector, _iDetector + 1, p); } + + //---------------------------------------------------------------------------------------- -// PROJECT BLOCK +// PROJECT BLOCK - vector projection geometry 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, inv_pixelLengthX; - float32 lengthPerCol, updatePerCol, inv_pixelLengthY; - 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 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 - 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 detX, detY, x, y, c, r, update_c, update_r, offset; + float32 lengthPerRow, lengthPerCol; + int iVolumeIndex, iRayIndex, row, col, iDetector; + + const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; + + bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); + if (vertical) { + lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); + update_c = -m_pVolumeGeometry->getPixelLengthY() * (proj->fRayX/proj->fRayY) * inv_pixelLengthX; + } else { + lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX); + update_r = -m_pVolumeGeometry->getPixelLengthX() * (proj->fRayY/proj->fRayX) * 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; - 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(); - // loop detectors for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { @@ -98,79 +97,70 @@ 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; - } + detX = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX; + detY = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY; // vertically - if (theta <= PIdiv4) { - + if (vertical) { + // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = m_pVolumeGeometry->coordXToColF(P) - 0.5f; + x = detX + (proj->fRayX/proj->fRayY)*(m_pVolumeGeometry->pixelRowToCenterY(0)-detY); + c = (x - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f; // 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; - - // 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); - } + for (row = 0; row < rowCount; ++row, c += update_c) { + + col = int(c); + offset = c - float32(col); + + if (col <= 0 || col >= colCount-1) continue; + + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, (1.0f - offset) * 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: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, offset * lengthPerRow); + p.pixelPosterior(iVolumeIndex); } + } } // 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->coordYToRowF(P) - 0.5f; + // calculate y for col 0 + y = detY + (proj->fRayY/proj->fRayX)*(m_pVolumeGeometry->pixelColToCenterX(0)-detX); + r = (m_pVolumeGeometry->getWindowMaxY() - y) * inv_pixelLengthY - 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); - } + // for each col + for (col = 0; col < colCount; ++col, r += update_r) { + + int row = int(r); + offset = r - float32(row); + + if (row <= 0 || row >= rowCount-1) continue; + + iVolumeIndex = row * colCount + col; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, (1.0f - offset) * 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); - } + + iVolumeIndex += colCount; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { + p.addWeight(iRayIndex, iVolumeIndex, offset * lengthPerCol); + p.pixelPosterior(iVolumeIndex); } + } } @@ -179,6 +169,4 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, } // end loop detector } // end loop angles - } - diff --git a/include/astra/ParallelBeamStripKernelProjector2D.h b/include/astra/ParallelBeamStripKernelProjector2D.h index 624bd3c..8fb665f 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.h +++ b/include/astra/ParallelBeamStripKernelProjector2D.h @@ -30,6 +30,7 @@ $Id$ #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 7d73888..85faaa3 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.inl +++ b/include/astra/ParallelBeamStripKernelProjector2D.inl @@ -30,7 +30,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> @@ -52,247 +52,177 @@ void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int 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 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); + const float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); + const float32 pixelArea = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + 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; + #pragma omp parallel for + for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { + + // variables + float32 detLX, detLY, detRX, detRY, S, T, update_c, update_r, offsetL, offsetR, invTminS; + float32 res, fRxOverRy, fRyOverRx; + 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) { + fRxOverRy = proj->fRayX/proj->fRayY; + update_c = -m_pVolumeGeometry->getPixelLengthY() * fRxOverRy * inv_pixelLengthX; + S = 0.5f - 0.5f*fabs(fRxOverRy); + T = 0.5f + 0.5f*fabs(fRxOverRy); + invTminS = 1.0f / (T-S); + } else { + fRyOverRx = proj->fRayY/proj->fRayX; + update_r = -m_pVolumeGeometry->getPixelLengthX() * fRyOverRx * inv_pixelLengthY; + S = 0.5f - 0.5f*fabs(fRyOverRx); + T = 0.5f + 0.5f*fabs(fRyOverRx); + invTminS = 1.0f / (T-S); } - // 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) { + // loop detectors + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { - iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) continue; + iRayIndex = iAngle * detCount + iDetector; - // 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; + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + detLX = proj->fDetSX + iDetector * proj->fDetUX; + detLY = proj->fDetSY + iDetector * proj->fDetUY; + detRX = detLX + proj->fDetUX; + detRY = detLY + proj->fDetUY; + + // vertically + if (vertical) { + + // calculate cL and cR for row 0 + float32 xL = detLX + fRxOverRy*(m_pVolumeGeometry->pixelRowToCenterY(0)-detLY); + float32 cL = (xL - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f; + float32 xR = detRX + fRxOverRy*(m_pVolumeGeometry->pixelRowToCenterY(0)-detRY); + float32 cR = (xR - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f; + + if (cR < cL) { + float32 tmp = cL; + cL = cR; + cR = tmp; + } + // 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 - } else { + for (row = 0; row < rowCount; ++row, cL += update_c, cR += update_c) { - // 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; + col_left = int(cL-0.5f+S); + col_right = int(cR+1.5-S); - updateX = cos_theta * inv_sin_theta; + if (col_left < 0) col_left = 0; + if (col_right > colCount-1) col_right = colCount-1; - // loop detectors - for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { - - iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; + offsetL = cL - float32(col_left); + offsetR = cR - float32(col_left); - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) continue; + // for each column + for (col = col_left; col <= col_right; ++col, offsetL -= 1.0f, offsetR -= 1.0f) { - // get t - t = m_pProjectionGeometry->indexToDetectorOffset(iDetector); - if (switch_t) t = -t; - - // 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; + iVolumeIndex = row * colCount + col; - // 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; + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { - // for each col - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { + // 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; - // get strip extremes in column indices - x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); - x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.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; - // 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; + p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res); + p.pixelPosterior(iVolumeIndex); + } + } + } + } - // for each affected col - for (row = x1L; row <= x1R; ++row) { + // horizontally + else { - if (row < 0 || row >= m_pVolumeGeometry->getGridRowCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + // calculate rL and rR for row 0 + float32 yL = detLY + fRyOverRx*(m_pVolumeGeometry->pixelColToCenterX(0)-detLX); + float32 rL = (m_pVolumeGeometry->getWindowMaxY() - yL) * inv_pixelLengthY - 0.5f; + float32 yR = detRY + fRyOverRx*(m_pVolumeGeometry->pixelColToCenterX(0)-detRX); + float32 rR = (m_pVolumeGeometry->getWindowMaxY() - yR) * inv_pixelLengthY - 0.5f; - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + if (rR < rL) { + float32 tmp = rL; + rL = rR; + rR = tmp; + } - // POLICY: PIXEL PRIOR - if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + // for each column + for (col = 0; col < colCount; ++col, rL += update_r, rR += update_r) { - // 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; } + row_top = int(rL-0.5f+S); + row_bottom = int(rR+1.5-S); - // POLICY: ADD - p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + if (row_top < 0) row_top = 0; + if (row_bottom > rowCount-1) row_bottom = rowCount-1; - // POLICY: PIXEL POSTERIOR - p.pixelPosterior(iVolumeIndex); + offsetL = rL - float32(row_top); + offsetR = rR - float32(row_top); - x2L -= 1.0f; - x2R -= 1.0f; + // for each row + for (row = row_top; row <= row_bottom; ++row, offsetL -= 1.0f, offsetR -= 1.0f) { - } // end row loop + iVolumeIndex = row * colCount + col; - } // end col loop + // POLICY: PIXEL PRIOR + ADD + POSTERIOR + if (p.pixelPrior(iVolumeIndex)) { - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); + // 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; - } // end detector loop + // 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; - } // end theta switch + p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res); + p.pixelPosterior(iVolumeIndex); + } + } + } + } + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } // end loop detector + } // end loop angles - } // end angle loop } - - diff --git a/include/astra/ParallelProjectionGeometry2D.h b/include/astra/ParallelProjectionGeometry2D.h index 2f7d36f..36b4b6f 100644 --- a/include/astra/ParallelProjectionGeometry2D.h +++ b/include/astra/ParallelProjectionGeometry2D.h @@ -30,6 +30,7 @@ $Id$ #define _INC_ASTRA_PARALLELPROJECTIONGEOMETRY2D #include "ProjectionGeometry2D.h" +#include "ParallelVecProjectionGeometry2D.h" namespace astra { @@ -134,7 +135,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. * @@ -151,7 +152,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 b8324e2..d26e6a7 100644 --- a/include/astra/ProjectionGeometry2D.h +++ b/include/astra/ProjectionGeometry2D.h @@ -65,10 +65,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. */ @@ -202,9 +198,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 @@ -273,12 +266,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 { diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m index 3875e6b..c3e737f 100644 --- a/matlab/algorithms/DART/IterativeTomography.m +++ b/matlab/algorithms/DART/IterativeTomography.m @@ -77,7 +77,7 @@ classdef IterativeTomography < matlab.mixin.Copyable function ok = initialize(this) % Initialize this object. Returns 1 if succesful. % >> tomography.initialize(); - + disp('sdfqnlmkqdsfmlkjdfqsjklm'); % create projection geometry with super-resolution if this.superresolution > 1 this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution); diff --git a/matlab/algorithms/DART/examples/example1.m b/matlab/algorithms/DART/examples/example1.m index cb02e0f..791e440 100644 --- a/matlab/algorithms/DART/examples/example1.m +++ b/matlab/algorithms/DART/examples/example1.m @@ -35,36 +35,41 @@ vol_geom = astra_create_vol_geom(det_count, det_count); [sinogram_id, sinogram] = astra_create_sino_cuda(I, proj_geom, vol_geom); astra_mex_data2d('delete', sinogram_id); -% DART -D = DARTalgorithm(sinogram, proj_geom); -D.t0 = 100; -D.t = 10; + % DART + D = DARTalgorithm(sinogram, proj_geom); + D.t0 = 100; + D.t = 10; -D.tomography.method = 'SIRT_CUDA'; -D.tomography.gpu_core = gpu_core; -D.tomography.use_minc = 'yes'; + D.tomography.method = 'SIRT'; + D.tomography.gpu_core = gpu_core; + D.tomography.use_minc = 'yes'; + D.tomography.gpu = 'no'; -D.segmentation.rho = rho; -D.segmentation.tau = tau; + D.segmentation = SegmentationPDM(); + D.segmentation.rho = rho*1.8; + D.segmentation.tau = tau*1.5; + D.segmentation.interval = 5; -D.smoothing.b = 0.1; -D.smoothing.gpu_core = gpu_core; - -D.masking.random = 0.1; -D.masking.gpu_core = gpu_core; + D.smoothing.b = 0.1; + D.smoothing.gpu_core = gpu_core; + D.smoothing.gpu = 'no'; + + D.masking.random = 0.1; + D.masking.gpu_core = gpu_core; + D.masking.gpu = 'no'; + + D.output.directory = outdir; + D.output.pre = [prefix '_']; + D.output.save_images = 'no'; + D.output.save_results = {'stats', 'settings', 'S', 'V'}; + D.output.save_interval = dart_iterations; + D.output.verbose = 'yes'; -D.output.directory = outdir; -D.output.pre = [prefix '_']; -D.output.save_images = 'no'; -D.output.save_results = {'stats', 'settings', 'S', 'V'}; -D.output.save_interval = dart_iterations; -D.output.verbose = 'yes'; + D.statistics.proj_diff = 'no'; -D.statistics.proj_diff = 'no'; + D.initialize(); -D.initialize(); - -D.iterate(dart_iterations); + D.iterate(dart_iterations); % save the reconstruction and the segmentation to file imwritesc(D.S, [outdir '/' prefix '_S.png']); diff --git a/matlab/mex/astra_mex_data2d_c.cpp b/matlab/mex/astra_mex_data2d_c.cpp index 909d229..6863b7a 100644 --- a/matlab/mex/astra_mex_data2d_c.cpp +++ b/matlab/mex/astra_mex_data2d_c.cpp @@ -45,6 +45,7 @@ $Id$ #include "astra/SparseMatrixProjectionGeometry2D.h" #include "astra/FanFlatProjectionGeometry2D.h" #include "astra/FanFlatVecProjectionGeometry2D.h" +#include "astra/ParallelVecProjectionGeometry2D.h" using namespace std; using namespace astra; @@ -160,6 +161,8 @@ void astra_mex_data2d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra pGeometry = new CFanFlatProjectionGeometry2D(); } else if (type == "fanflat_vec") { pGeometry = new CFanFlatVecProjectionGeometry2D(); + } else if (type == "parallel_vec") { + pGeometry = new CParallelVecProjectionGeometry2D(); } else { pGeometry = new CParallelProjectionGeometry2D(); } @@ -449,6 +452,8 @@ void astra_mex_data2d_change_geometry(int nlhs, mxArray* plhs[], int nrhs, const pGeometry = new CFanFlatProjectionGeometry2D(); } else if (type == "fanflat_vec") { pGeometry = new CFanFlatVecProjectionGeometry2D(); + } else if (type == "parallel_vec") { + pGeometry = new CParallelVecProjectionGeometry2D(); } else { pGeometry = new CParallelProjectionGeometry2D(); } diff --git a/matlab/tools/astra_create_proj_geom.m b/matlab/tools/astra_create_proj_geom.m index 862e410..ff7d74d 100644 --- a/matlab/tools/astra_create_proj_geom.m +++ b/matlab/tools/astra_create_proj_geom.m @@ -107,6 +107,19 @@ if strcmp(type,'parallel') 'ProjectionAngles', varargin{3} ... ); +elseif strcmp(type,'parallel_vec') + if numel(varargin) < 2 + error('not enough variables: astra_create_proj_geom(parallel_vec, det_count, V') + end + if size(varargin{2}, 2) ~= 6 + error('V should be a Nx6 matrix, with N the number of projections') + end + proj_geom = struct( ... + 'type', 'parallel_vec', ... + 'DetectorCount', varargin{1}, ... + 'Vectors', varargin{2} ... + ); + elseif strcmp(type,'fanflat') if numel(varargin) < 5 error('not enough variables: astra_create_proj_geom(fanflat, det_width, det_count, angles, source_origin, source_det)'); diff --git a/matlab/tools/astra_geom_2vec.m b/matlab/tools/astra_geom_2vec.m index 0abd07c..e563f47 100644 --- a/matlab/tools/astra_geom_2vec.m +++ b/matlab/tools/astra_geom_2vec.m @@ -1,25 +1,65 @@ -function proj_geom_out = astra_geom_2vec(proj_geom) +function proj_geom_vec = astra_geom_2vec(proj_geom) +%-------------------------------------------------------------------------- +% proj_geom_vec = astra_geom_2vec(proj_geom) +% +% Convert a conventional projection geometry to a corresponding vector-base +% projection geometry +% +% proj_geom: input projection geometry (parallel, fanflat, parallel3d, cone) +% proj_geom_vec: output vector-base projection geometry +%-------------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +%-------------------------------------------------------------------------- +% $Id$ + + % PARALLEL 2D + if strcmp(proj_geom.type,'parallel') + + vectors = zeros(numel(proj_geom.ProjectionAngles), 6); + for i = 1:numel(proj_geom.ProjectionAngles) + + % ray direction + vectors(i,1) = sin(proj_geom.ProjectionAngles(i)); + vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)); + + % center of detector + vectors(i,3) = 0; + vectors(i,4) = 0; + + % vector from detector pixel 0 to 1 + vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth; + vectors(i,6) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth; + + end + + proj_geom_vec = astra_create_proj_geom('parallel_vec', proj_geom.DetectorCount, vectors); % FANFLAT - if strcmp(proj_geom.type,'fanflat') + elseif strcmp(proj_geom.type,'fanflat') vectors = zeros(numel(proj_geom.ProjectionAngles), 6); for i = 1:numel(proj_geom.ProjectionAngles) % source vectors(i,1) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource; - vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource; + vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource; % center of detector vectors(i,3) = -sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector; - vectors(i,4) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector; + vectors(i,4) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector; % vector from detector pixel 0 to 1 vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth; vectors(i,6) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth; end - proj_geom_out = astra_create_proj_geom('fanflat_vec', proj_geom.DetectorCount, vectors); + proj_geom_vec = astra_create_proj_geom('fanflat_vec', proj_geom.DetectorCount, vectors); % CONE elseif strcmp(proj_geom.type,'cone') @@ -29,13 +69,13 @@ function proj_geom_out = astra_geom_2vec(proj_geom) % source vectors(i,1) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource; - vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource; - vectors(i,3) = 0; + vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource; + vectors(i,3) = 0; % center of detector vectors(i,4) = -sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector; - vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector; - vectors(i,6) = 0; + vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector; + vectors(i,6) = 0; % vector from detector pixel (0,0) to (0,1) vectors(i,7) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorSpacingX; @@ -45,13 +85,13 @@ function proj_geom_out = astra_geom_2vec(proj_geom) % vector from detector pixel (0,0) to (1,0) vectors(i,10) = 0; vectors(i,11) = 0; - vectors(i,12) = proj_geom.DetectorSpacingY; + vectors(i,12) = proj_geom.DetectorSpacingY; end - proj_geom_out = astra_create_proj_geom('cone_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors); + proj_geom_vec = astra_create_proj_geom('cone_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors); - % PARALLEL - elseif strcmp(proj_geom.type,'parallel3d') + % PARALLEL 3D + elseif strcmp(proj_geom.type,'parallel3d') vectors = zeros(numel(proj_geom.ProjectionAngles), 12); for i = 1:numel(proj_geom.ProjectionAngles) @@ -59,12 +99,12 @@ function proj_geom_out = astra_geom_2vec(proj_geom) % ray direction vectors(i,1) = sin(proj_geom.ProjectionAngles(i)); vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)); - vectors(i,3) = 0; + vectors(i,3) = 0; % center of detector vectors(i,4) = 0; vectors(i,5) = 0; - vectors(i,6) = 0; + vectors(i,6) = 0; % vector from detector pixel (0,0) to (0,1) vectors(i,7) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorSpacingX; @@ -74,11 +114,13 @@ function proj_geom_out = astra_geom_2vec(proj_geom) % vector from detector pixel (0,0) to (1,0) vectors(i,10) = 0; vectors(i,11) = 0; - vectors(i,12) = proj_geom.DetectorSpacingY; + vectors(i,12) = proj_geom.DetectorSpacingY; end - proj_geom_out = astra_create_proj_geom('parallel3d_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors); + proj_geom_vec = astra_create_proj_geom('parallel3d_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors); else error(['No suitable vector geometry found for type: ' proj_geom.type]) end + +end diff --git a/matlab/tools/astra_geom_postalignment.m b/matlab/tools/astra_geom_postalignment.m index 4115af2..c16f8ed 100644 --- a/matlab/tools/astra_geom_postalignment.m +++ b/matlab/tools/astra_geom_postalignment.m @@ -1,11 +1,36 @@ function proj_geom = astra_geom_postalignment(proj_geom, factor) +%-------------------------------------------------------------------------- +% proj_geom = astra_geom_postalignment(proj_geom, factorU) +% proj_geom = astra_geom_postalignment(proj_geom, [factorU factorV]) +% +% Apply a postalignment to a vector-based projection geometry. Can be used to model the rotation axis offset. +% +% proj_geom: input projection geometry (vector-based only, use astra_geom_2vec to convert conventional projection geometries) +% dim (optional): which dimension +% s: output +%-------------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +%-------------------------------------------------------------------------- +% $Id$ + + if strcmp(proj_geom.type,'fanflat_vec') || strcmp(proj_geom.type,'parallel_vec') + proj_geom.Vectors(:,3:4) = proj_geom.Vectors(:,3:4) + factor(1) * proj_geom.Vectors(:,5:6); - if strcmp(proj_geom.type,'fanflat_vec') - proj_geom.Vectors(:,3:4) = proj_geom.Vectors(:,3:4) + factor * proj_geom.Vectors(:,5:6); - elseif strcmp(proj_geom.type,'cone_vec') || strcmp(proj_geom.type,'parallel3d_vec') - proj_geom.Vectors(:,4:6) = proj_geom.Vectors(:,4:6) + factor * proj_geom.Vectors(:,7:9); + if numel(factor) == 1 + proj_geom.Vectors(:,4:6) = proj_geom.Vectors(:,4:6) + factor * proj_geom.Vectors(:,7:9); + elseif numel(factor) > 1 + proj_geom.Vectors(:,4:6) = proj_geom.Vectors(:,4:6) + factor(1) * proj_geom.Vectors(:,7:9) + factor(2) * proj_geom.Vectors(:,10:12); + end else error('Projection geometry not suited for postalignment correction.') end + +end
\ No newline at end of file diff --git a/matlab/tools/astra_geom_size.m b/matlab/tools/astra_geom_size.m index 7044515..b3c1522 100644 --- a/matlab/tools/astra_geom_size.m +++ b/matlab/tools/astra_geom_size.m @@ -1,4 +1,22 @@ function s = astra_geom_size(geom, dim) +%-------------------------------------------------------------------------- +% s = astra_geom_size(geom, dim) +% +% Get the size of a volume or projection geometry. +% +% geom: volume or projection geometry +% dim (optional): which dimension +% s: output +%-------------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +%-------------------------------------------------------------------------- +% $Id$ if isfield(geom, 'GridSliceCount') % 3D Volume geometry? @@ -6,23 +24,23 @@ function s = astra_geom_size(geom, dim) elseif isfield(geom, 'GridColCount') % 2D Volume geometry? s = [ geom.GridRowCount, geom.GridColCount ]; - elseif strcmp(geom.type,'parallel') || strcmp(geom.type,'fanflat') + elseif strcmp(geom.type,'parallel') || strcmp(geom.type,'fanflat') s = [numel(geom.ProjectionAngles), geom.DetectorCount]; - - elseif strcmp(geom.type,'parallel3d') || strcmp(geom.type,'cone') + + elseif strcmp(geom.type,'parallel3d') || strcmp(geom.type,'cone') s = [geom.DetectorColCount, numel(geom.ProjectionAngles), geom.DetectorRowCount]; - - elseif strcmp(geom.type,'fanflat_vec') + + elseif strcmp(geom.type,'fanflat_vec') || strcmp(geom.type,'parallel_vec') s = [size(geom.Vectors,1), geom.DetectorCount]; - - elseif strcmp(geom.type,'parallel3d_vec') || strcmp(geom.type,'cone_vec') + + elseif strcmp(geom.type,'parallel3d_vec') || strcmp(geom.type,'cone_vec') s = [geom.DetectorColCount, size(geom.Vectors,1), geom.DetectorRowCount]; - + end if nargin == 2 s = s(dim); end - + end diff --git a/matlab/tools/astra_geom_visualize.m b/matlab/tools/astra_geom_visualize.m new file mode 100644 index 0000000..0044844 --- /dev/null +++ b/matlab/tools/astra_geom_visualize.m @@ -0,0 +1,216 @@ +function astra_geom_visualize(proj_geom, vol_geom) + + if strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') ||strcmp(proj_geom.type,'parallel3d') || strcmp(proj_geom.type,'cone') + proj_geom = astra_geom_2vec(proj_geom); + end + + % open window + f = figure('Visible','off'); + hold on + + % display projection 1 + displayProjection(1); + + % label + txt = uicontrol('Style','text', 'Position',[10 10 70 20], 'String','Projection'); + + % slider + anglecount = size(proj_geom.Vectors,1); + sld = uicontrol('Style', 'slider', ... + 'Min', 1, 'Max', anglecount, 'SliderStep', [1 1]/anglecount, 'Value', 1, ... + 'Position', [80 10 200 20], ... + 'Callback', @updateProjection); + + f.Visible = 'on'; + + function updateProjection(source, callbackdata) + displayProjection(floor(source.Value)); + end + + function displayProjection(a) + + colours = get(gca,'ColorOrder'); + + + % set title + title(['projection ' num2str(a)]); + + if strcmp(proj_geom.type,'parallel_vec') + + v = proj_geom.Vectors; + d = proj_geom.DetectorCount; + + if ~isfield(vol_geom, 'option') + minx = -vol_geom.GridRowCount/2; + miny = -vol_geom.GridColCount/2; + minz = -vol_geom.GridSliceCount/2; + maxx = vol_geom.GridRowCount/2; + else + minx = vol_geom.option.WindowMinX; + miny = vol_geom.option.WindowMinY; + maxx = vol_geom.option.WindowMaxX; + maxy = vol_geom.option.WindowMaxY; + end + + % axis + cla + axis([minx maxx miny maxy]*2.25) + axis square + + % volume + plot([minx minx maxx maxx minx], [miny maxy maxy miny miny], 'LineWidth', 1, 'Color', colours(1,:)) + + % ray + s = maxx - minx; + plot([0 v(a,1)]*s*0.33, [0 v(a,2)]*s*0.33, 'LineWidth', 2, 'Color', colours(3,:)) + + % detector + s2 = s*0.75; + plot([-d/2 d/2]*v(a,5) + v(a,3) + s2*v(a,1), [-d/2 d/2]*v(a,6) + v(a,4) + s2*v(a,2), 'LineWidth', 2, 'Color', colours(5,:)) + + elseif strcmp(proj_geom.type,'fanflat_vec') + + v = proj_geom.Vectors; + d = proj_geom.DetectorCount; + + if ~isfield(vol_geom, 'option') + minx = -vol_geom.GridRowCount/2; + miny = -vol_geom.GridColCount/2; + minz = -vol_geom.GridSliceCount/2; + maxx = vol_geom.GridRowCount/2; + else + minx = vol_geom.option.WindowMinX; + miny = vol_geom.option.WindowMinY; + maxx = vol_geom.option.WindowMaxX; + maxy = vol_geom.option.WindowMaxY; + end + + % axis + cla + axis([minx maxx miny maxy]*2.25) + axis square + + % volume + plot([minx minx maxx maxx minx], [miny maxy maxy miny miny], 'LineWidth', 1, 'Color', colours(1,:)) + + % detector + D1 = v(a,3:4) - d/2*v(a,5:6); + D2 = v(a,3:4) + d/2*v(a,5:6); + plot([D1(1) D2(1)], [D1(2) D2(2)], 'LineWidth', 2, 'Color', colours(5,:)) + + % beam + plot([v(a,1) D1(1)], [v(a,2) D1(2)], 'LineWidth', 1, 'Color', colours(3,:)) + plot([v(a,1) D2(1)], [v(a,2) D2(2)], 'LineWidth', 1, 'Color', colours(3,:)) + + elseif strcmp(proj_geom.type,'parallel3d_vec') + + v = proj_geom.Vectors; + d1 = proj_geom.DetectorRowCount; + d2 = proj_geom.DetectorColCount; + + if ~isfield(vol_geom, 'option') + minx = -vol_geom.GridRowCount/2; + miny = -vol_geom.GridColCount/2; + minz = -vol_geom.GridSliceCount/2; + maxx = vol_geom.GridRowCount/2; + maxy = vol_geom.GridColCount/2; + maxz = vol_geom.GridSliceCount/2; + else + minx = vol_geom.option.WindowMinX; + miny = vol_geom.option.WindowMinY; + minz = vol_geom.option.WindowMinZ; + maxx = vol_geom.option.WindowMaxX; + maxy = vol_geom.option.WindowMaxY; + maxz = vol_geom.option.WindowMaxZ; + end + + % axis + windowminx = min(v(:,4)); + windowminy = min(v(:,5)); + windowminz = max(v(:,6)); + windowmaxx = max(v(:,4)); + windowmaxy = max(v(:,5)); + windowmaxz = max(v(:,6)); + cla + axis([minx maxx miny maxy minz maxz]*5.10) + + % volume + plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [minz minz minz minz minz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [maxz maxz maxz maxz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([minx minx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([maxx maxx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([minx minx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([maxx maxx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + + % detector + D1 = v(a,4:6) - d1/2*v(a,7:9) - d2/2*v(a,10:12); + D2 = v(a,4:6) + d1/2*v(a,7:9) - d2/2*v(a,10:12); + D3 = v(a,4:6) + d1/2*v(a,7:9) + d2/2*v(a,10:12); + D4 = v(a,4:6) - d1/2*v(a,7:9) + d2/2*v(a,10:12); + plot3([D1(1) D2(1) D3(1) D4(1) D1(1)], [D1(2) D2(2) D3(2) D4(2) D1(2)], [D1(3) D2(3) D3(3) D4(3) D1(3)], 'LineWidth', 2, 'Color', colours(5,:)) + + % ray + s = maxx - minx; + plot3([0 v(a,1)]*s*0.30, [0 v(a,2)]*s*0.30, [0 v(a,3)]*s*0.30, 'LineWidth', 2, 'Color', colours(3,:)) + + elseif strcmp(proj_geom.type,'cone_vec') + + v = proj_geom.Vectors; + d1 = proj_geom.DetectorRowCount; + d2 = proj_geom.DetectorColCount; + + if ~isfield(vol_geom, 'option') + minx = -vol_geom.GridRowCount/2; + miny = -vol_geom.GridColCount/2; + minz = -vol_geom.GridSliceCount/2; + maxx = vol_geom.GridRowCount/2; + maxy = vol_geom.GridColCount/2; + maxz = vol_geom.GridSliceCount/2; + else + minx = vol_geom.option.WindowMinX; + miny = vol_geom.option.WindowMinY; + minz = vol_geom.option.WindowMinZ; + maxx = vol_geom.option.WindowMaxX; + maxy = vol_geom.option.WindowMaxY; + maxz = vol_geom.option.WindowMaxZ; + end + + % axis + windowminx = min(v(:,4)); + windowminy = min(v(:,5)); + windowminz = max(v(:,6)); + windowmaxx = max(v(:,4)); + windowmaxy = max(v(:,5)); + windowmaxz = max(v(:,6)); + cla + axis([minx maxx miny maxy minz maxz]*5.10) + + % volume + plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [minz minz minz minz minz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [maxz maxz maxz maxz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([minx minx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([maxx maxx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([minx minx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + plot3([maxx maxx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:)) + + % detector + D1 = v(a,4:6) - d1/2*v(a,7:9) - d2/2*v(a,10:12); + D2 = v(a,4:6) + d1/2*v(a,7:9) - d2/2*v(a,10:12); + D3 = v(a,4:6) + d1/2*v(a,7:9) + d2/2*v(a,10:12); + D4 = v(a,4:6) - d1/2*v(a,7:9) + d2/2*v(a,10:12); + plot3([D1(1) D2(1) D3(1) D4(1) D1(1)], [D1(2) D2(2) D3(2) D4(2) D1(2)], [D1(3) D2(3) D3(3) D4(3) D1(3)], 'LineWidth', 2, 'Color', colours(5,:)) + + % beam + plot3([v(a,1) D1(1)], [v(a,2) D1(2)], [v(a,3) D1(3)], 'LineWidth', 1, 'Color', colours(3,:)) + plot3([v(a,1) D2(1)], [v(a,2) D2(2)], [v(a,3) D2(3)], 'LineWidth', 1, 'Color', colours(3,:)) + plot3([v(a,1) D3(1)], [v(a,2) D3(2)], [v(a,3) D3(3)], 'LineWidth', 1, 'Color', colours(3,:)) + plot3([v(a,1) D4(1)], [v(a,2) D4(2)], [v(a,3) D4(3)], 'LineWidth', 1, 'Color', colours(3,:)) + + + else + error('invalid projector type') + + end + end + +end diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 909f58f..540aad4 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -126,6 +126,10 @@ cdef extern from "astra/FanFlatVecProjectionGeometry2D.h" namespace "astra": cdef cppclass CFanFlatVecProjectionGeometry2D: CFanFlatVecProjectionGeometry2D() +cdef extern from "astra/ParallelVecProjectionGeometry2D.h" namespace "astra": + cdef cppclass CParallelVecProjectionGeometry2D: + CParallelVecProjectionGeometry2D() + cdef extern from "astra/ParallelProjectionGeometry2D.h" namespace "astra": cdef cppclass CParallelProjectionGeometry2D: CParallelProjectionGeometry2D() diff --git a/python/astra/creators.py b/python/astra/creators.py index 68bc8a2..4cd7f2d 100644 --- a/python/astra/creators.py +++ b/python/astra/creators.py @@ -140,6 +140,13 @@ This method can be called in a number of ways: :type angles: :class:`numpy.ndarray` :returns: A parallel projection geometry. +``create_proj_geom('parallel_vec', det_count, V)``: + +:param det_count: Number of detector pixels. +:type det_count: :class:`int` +:param V: Vector array. +:type V: :class:`numpy.ndarray` +:returns: A parallel-beam projection geometry. ``create_proj_geom('fanflat', det_width, det_count, angles, source_origin, source_det)``: @@ -227,6 +234,12 @@ This method can be called in a number of ways: raise Exception( 'not enough variables: astra_create_proj_geom(parallel, detector_spacing, det_count, angles)') return {'type': 'parallel', 'DetectorWidth': args[0], 'DetectorCount': args[1], 'ProjectionAngles': args[2]} + elif intype == 'parallel_vec': + if len(args) < 2: + raise Exception('not enough variables: astra_create_proj_geom(parallel_vec, det_count, V)') + if not args[1].shape[1] == 6: + raise Exception('V should be a Nx6 matrix, with N the number of projections') + return {'type':'parallel_vec', 'DetectorCount':args[0], 'Vectors':args[1]} elif intype == 'fanflat': if len(args) < 5: raise Exception('not enough variables: astra_create_proj_geom(fanflat, det_width, det_count, angles, source_origin, source_det)') diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index 4919bf2..42854b3 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -104,6 +104,8 @@ def create(datatype, geometry, data=None, link=False): ppGeometry = <CProjectionGeometry2D * >new CFanFlatProjectionGeometry2D() elif (tpe == 'fanflat_vec'): ppGeometry = <CProjectionGeometry2D * >new CFanFlatVecProjectionGeometry2D() + elif (tpe == 'parallel_vec'): + ppGeometry = <CProjectionGeometry2D * >new CParallelVecProjectionGeometry2D() else: ppGeometry = <CProjectionGeometry2D * >new CParallelProjectionGeometry2D() if not ppGeometry.initialize(cfg[0]): @@ -197,6 +199,8 @@ def change_geometry(i, geom): ppGeometry = <CProjectionGeometry2D * >new CFanFlatProjectionGeometry2D() elif (tpe == 'fanflat_vec'): ppGeometry = <CProjectionGeometry2D * >new CFanFlatVecProjectionGeometry2D() + elif (tpe == 'parallel_vec'): + ppGeometry = <CProjectionGeometry2D * >new CParallelVecProjectionGeometry2D() else: ppGeometry = <CProjectionGeometry2D * >new CParallelProjectionGeometry2D() if not ppGeometry.initialize(cfg[0]): diff --git a/python/astra/functions.py b/python/astra/functions.py index e38b5bc..13ca911 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -172,7 +172,27 @@ def geom_2vec(proj_geom): :param proj_geom: Projection geometry to convert :type proj_geom: :class:`dict` """ - if proj_geom['type'] == 'fanflat': + + if proj_geom['type'] == 'parallel': + angles = proj_geom['ProjectionAngles'] + vectors = np.zeros((len(angles), 6)) + for i in range(len(angles)): + + # source + vectors[i, 0] = np.sin(angles[i]) + vectors[i, 1] = -np.cos(angles[i]) + + # center of detector + vectors[i, 2] = 0 + vectors[i, 3] = 0 + + # vector from detector pixel 0 to 1 + vectors[i, 4] = np.cos(angles[i]) * proj_geom['DetectorWidth'] + vectors[i, 5] = np.sin(angles[i]) * proj_geom['DetectorWidth'] + proj_geom_out = ac.create_proj_geom( + 'parallel_vec', proj_geom['DetectorCount'], vectors) + + elif proj_geom['type'] == 'fanflat': angles = proj_geom['ProjectionAngles'] vectors = np.zeros((len(angles), 6)) for i in range(len(angles)): @@ -251,3 +271,37 @@ def geom_2vec(proj_geom): raise ValueError( 'No suitable vector geometry found for type: ' + proj_geom['type']) return proj_geom_out + + +def geom_postalignment(proj_geom, factor): + """Returns the size of a volume or sinogram, based on the projection or volume geometry. + + :param proj_geom: input projection geometry (vector-based only, use astra.geom_2vec to convert conventional projection geometries) + :type proj_geom: :class:`dict` + :param factor: Optional axis index to return + :type factor: :class:`float` + """ + + if proj_geom['type'] == 'parallel_vec' or proj_geom['type'] == 'fanflat_vec': + for i in range(proj_geom['Vectors'].shape[0]): + proj_geom['Vectors'][i,2] = proj_geom['Vectors'][i,2] + factor * proj_geom['Vectors'][i,4]; + proj_geom['Vectors'][i,3] = proj_geom['Vectors'][i,3] + factor * proj_geom['Vectors'][i,5]; + + elif proj_geom['type'] == 'parallel3d_vec' or proj_geom['type'] == 'cone_vec': + + if len(factor) == 1: + for i in range(proj_geom['Vectors'].shape[0]): + proj_geom['Vectors'][i,3] = proj_geom['Vectors'][i,3] + factor * proj_geom['Vectors'][i,6]; + proj_geom['Vectors'][i,4] = proj_geom['Vectors'][i,4] + factor * proj_geom['Vectors'][i,7]; + proj_geom['Vectors'][i,5] = proj_geom['Vectors'][i,5] + factor * proj_geom['Vectors'][i,8]; + + elif len(factor) > 1: + for i in range(proj_geom['Vectors'].shape[0]): + proj_geom['Vectors'][i,3] = proj_geom['Vectors'][i,3] + factor[0] * proj_geom['Vectors'][i,6] + factor[1] * proj_geom['Vectors'][i, 9]; + proj_geom['Vectors'][i,4] = proj_geom['Vectors'][i,4] + factor[0] * proj_geom['Vectors'][i,7] + factor[1] * proj_geom['Vectors'][i,10]; + proj_geom['Vectors'][i,5] = proj_geom['Vectors'][i,5] + factor[0] * proj_geom['Vectors'][i,8] + factor[1] * proj_geom['Vectors'][i,11]; + else: + raise ValueError('No suitable geometry for postalignment: ' + proj_geom['type']) + + return proj_geom + diff --git a/src/FanFlatProjectionGeometry2D.cpp b/src/FanFlatProjectionGeometry2D.cpp index 32a19bc..91b9578 100644 --- a/src/FanFlatProjectionGeometry2D.cpp +++ b/src/FanFlatProjectionGeometry2D.cpp @@ -215,7 +215,25 @@ Config* CFanFlatProjectionGeometry2D::getConfiguration() const cfg->self.addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); return cfg; } + //---------------------------------------------------------------------------------------- +CFanFlatVecProjectionGeometry2D* CFanFlatProjectionGeometry2D::toVectorGeometry() +{ + SFanProjection* vectors = new SFanProjection[m_iProjectionAngleCount]; + for (int i = 0; i < m_iProjectionAngleCount; ++i) + { + vectors[i].fSrcX = sinf(m_pfProjectionAngles[i]) * m_fOriginSourceDistance; + vectors[i].fSrcY = -cosf(m_pfProjectionAngles[i]) * m_fOriginSourceDistance; + vectors[i].fDetUX = cosf(m_pfProjectionAngles[i]) * m_fDetectorWidth; + vectors[i].fDetUY = sinf(m_pfProjectionAngles[i]) * m_fDetectorWidth; + vectors[i].fDetSX = -sinf(m_pfProjectionAngles[i]) * m_fOriginDetectorDistance - 0.5f * m_iDetectorCount * vectors[i].fDetUX; + vectors[i].fDetSY = cosf(m_pfProjectionAngles[i]) * m_fOriginDetectorDistance - 0.5f * m_iDetectorCount * vectors[i].fDetUY; + } + CFanFlatVecProjectionGeometry2D* vecGeom = new CFanFlatVecProjectionGeometry2D(); + vecGeom->initialize(m_iProjectionAngleCount, m_iDetectorCount, vectors); + delete[] vectors; + return vecGeom; +} } // namespace astra diff --git a/src/ParallelBeamBlobKernelProjector2D.cpp b/src/ParallelBeamBlobKernelProjector2D.cpp index 4559a48..1e6312f 100644 --- a/src/ParallelBeamBlobKernelProjector2D.cpp +++ b/src/ParallelBeamBlobKernelProjector2D.cpp @@ -102,7 +102,7 @@ bool CParallelBeamBlobKernelProjector2D::_check() // check base class ASTRA_CONFIG_CHECK(CProjector2D::_check(), "ParallelBeamBlobKernelProjector2D", "Error in Projector2D initialization"); - ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamBlobKernelProjector2D", "Unsupported projection geometry"); + ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry) || dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamBlobKernelProjector2D", "Unsupported projection geometry"); ASTRA_CONFIG_CHECK(m_iBlobSampleCount > 0, "ParallelBeamBlobKernelProjector2D", "m_iBlobSampleCount should be strictly positive."); ASTRA_CONFIG_CHECK(m_pfBlobValues, "ParallelBeamBlobKernelProjector2D", "Invalid Volume Geometry Object."); @@ -155,17 +155,6 @@ bool CParallelBeamBlobKernelProjector2D::initialize(const Config& _cfg) for (int i = 0; i < m_iBlobSampleCount; i++) { m_pfBlobValues[i] = values[i]; } - - // Required: KernelValues - node2 = node.getSingleNode("KernelValuesNeg"); - ASTRA_CONFIG_CHECK(node2, "BlobProjector", "No Kernel/KernelValuesNeg tag specified."); - vector<float32> values2 = node2.getContentNumericalArray(); - ASTRA_CONFIG_CHECK(values2.size() == (unsigned int)m_iBlobSampleCount, "BlobProjector", "Number of specified values doesn't match SampleCount."); - m_pfBlobValuesNeg = new float32[m_iBlobSampleCount]; - for (int i = 0; i < m_iBlobSampleCount; i++) { - m_pfBlobValuesNeg[i] = values2[i]; - } - } // success @@ -228,44 +217,44 @@ void CParallelBeamBlobKernelProjector2D::computeSingleRayWeights(int _iProjectio // Splat a single point std::vector<SDetector2D> CParallelBeamBlobKernelProjector2D::projectPoint(int _iRow, int _iCol) { - float32 x = m_pVolumeGeometry->pixelColToCenterX(_iCol); - float32 y = m_pVolumeGeometry->pixelRowToCenterY(_iRow); - - std::vector<SDetector2D> res; - // loop projectors and detectors - for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // get projection angle - float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - if (theta >= 7*PIdiv4) theta -= 2*PI; - bool inverse = false; - if (theta >= 3*PIdiv4) { - theta -= PI; - inverse = true; - } - - // calculate distance from the center of the voxel to the ray though the origin - float32 t = x * cos(theta) + y * sin(theta); - if (inverse) t *= -1.0f; - - // calculate the offset on the detectorarray (in indices) - float32 d = m_pProjectionGeometry->detectorOffsetToIndexFloat(t); - int dmin = (int)ceil(d - m_fBlobSize); - int dmax = (int)floor(d + m_fBlobSize); - - // add affected detectors to the list - for (int i = dmin; i <= dmax; ++i) { - if (d >= 0 && d < m_pProjectionGeometry->getDetectorCount()) { - SDetector2D det; - det.m_iAngleIndex = iProjection; - det.m_iDetectorIndex = i; - det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - res.push_back(det); - } - } - } - - // return result vector - return res; - + // float32 x = m_pVolumeGeometry->pixelColToCenterX(_iCol); + // float32 y = m_pVolumeGeometry->pixelRowToCenterY(_iRow); + + // std::vector<SDetector2D> res; + // // loop projectors and detectors + // for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { + + // // get projection angle + // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); + // if (theta >= 7*PIdiv4) theta -= 2*PI; + // bool inverse = false; + // if (theta >= 3*PIdiv4) { + // theta -= PI; + // inverse = true; + // } + + // // calculate distance from the center of the voxel to the ray though the origin + // float32 t = x * cos(theta) + y * sin(theta); + // if (inverse) t *= -1.0f; + + // // calculate the offset on the detectorarray (in indices) + // float32 d = m_pProjectionGeometry->detectorOffsetToIndexFloat(t); + // int dmin = (int)ceil(d - m_fBlobSize); + // int dmax = (int)floor(d + m_fBlobSize); + + // // add affected detectors to the list + // for (int i = dmin; i <= dmax; ++i) { + // if (d >= 0 && d < m_pProjectionGeometry->getDetectorCount()) { + // SDetector2D det; + // det.m_iAngleIndex = iProjection; + // det.m_iDetectorIndex = i; + // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; + // res.push_back(det); + // } + // } + // } + + // // return result vector + // return res; + return std::vector<SDetector2D>(); } diff --git a/src/ParallelBeamLineKernelProjector2D.cpp b/src/ParallelBeamLineKernelProjector2D.cpp index 5a23413..c94c024 100644 --- a/src/ParallelBeamLineKernelProjector2D.cpp +++ b/src/ParallelBeamLineKernelProjector2D.cpp @@ -88,7 +88,7 @@ bool CParallelBeamLineKernelProjector2D::_check() // check base class ASTRA_CONFIG_CHECK(CProjector2D::_check(), "ParallelBeamLineKernelProjector2D", "Error in Projector2D initialization"); - ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamLineKernelProjector2D", "Unsupported projection geometry"); + ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry) || dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamLineKernelProjector2D", "Unsupported projection geometry"); ASTRA_CONFIG_CHECK(m_pVolumeGeometry->getPixelLengthX() == m_pVolumeGeometry->getPixelLengthY(), "ParallelBeamLineKernelProjector2D", "Pixel height must equal pixel width."); @@ -162,61 +162,63 @@ void CParallelBeamLineKernelProjector2D::computeSingleRayWeights(int _iProjectio // Project Point std::vector<SDetector2D> CParallelBeamLineKernelProjector2D::projectPoint(int _iRow, int _iCol) { - float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - std::vector<SDetector2D> res; - // loop projectors and detectors - for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // get projection angle - float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - if (theta >= 7*PIdiv4) theta -= 2*PI; - bool inverse = false; - if (theta >= 3*PIdiv4) { - theta -= PI; - inverse = true; - } - - // calculate distance from the center of the voxel to the ray though the origin - float32 tUL = xUL * cos(theta) + yUL * sin(theta); - float32 tUR = xUR * cos(theta) + yUR * sin(theta); - float32 tLL = xLL * cos(theta) + yLL * sin(theta); - float32 tLR = xLR * cos(theta) + yLR * sin(theta); - if (inverse) { - tUL *= -1.0f; - tUR *= -1.0f; - tLL *= -1.0f; - tLR *= -1.0f; - } - float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); - float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); - - // calculate the offset on the detectorarray (in indices) - int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); - int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); - - // add affected detectors to the list - for (int i = dmin; i <= dmax; ++i) { - if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { - SDetector2D det; - det.m_iAngleIndex = iProjection; - det.m_iDetectorIndex = i; - det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - res.push_back(det); - } - } - } - - // return result vector return res; + // float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; + + // std::vector<SDetector2D> res; + // // loop projectors and detectors + // for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { + + // // get projection angle + // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); + // if (theta >= 7*PIdiv4) theta -= 2*PI; + // bool inverse = false; + // if (theta >= 3*PIdiv4) { + // theta -= PI; + // inverse = true; + // } + + // // calculate distance from the center of the voxel to the ray though the origin + // float32 tUL = xUL * cos(theta) + yUL * sin(theta); + // float32 tUR = xUR * cos(theta) + yUR * sin(theta); + // float32 tLL = xLL * cos(theta) + yLL * sin(theta); + // float32 tLR = xLR * cos(theta) + yLR * sin(theta); + // if (inverse) { + // tUL *= -1.0f; + // tUR *= -1.0f; + // tLL *= -1.0f; + // tLR *= -1.0f; + // } + // float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); + // float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); + + // // calculate the offset on the detectorarray (in indices) + // int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); + // int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); + + // // add affected detectors to the list + // for (int i = dmin; i <= dmax; ++i) { + // if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { + // SDetector2D det; + // det.m_iAngleIndex = iProjection; + // det.m_iDetectorIndex = i; + // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; + // res.push_back(det); + // } + // } + // } + + // // return result vector + // return res; } //---------------------------------------------------------------------------------------- diff --git a/src/ParallelBeamLinearKernelProjector2D.cpp b/src/ParallelBeamLinearKernelProjector2D.cpp index a710664..043db71 100644 --- a/src/ParallelBeamLinearKernelProjector2D.cpp +++ b/src/ParallelBeamLinearKernelProjector2D.cpp @@ -89,7 +89,7 @@ bool CParallelBeamLinearKernelProjector2D::_check() // check base class ASTRA_CONFIG_CHECK(CProjector2D::_check(), "ParallelBeamLinearKernelProjector2D", "Error in Projector2D initialization"); - ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamLinearKernelProjector2D", "Unsupported projection geometry"); + ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry) || dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamLinearKernelProjector2D", "Unsupported projection geometry"); /// TODO: ADD PIXEL H/W LIMITATIONS ASTRA_CONFIG_CHECK(m_pVolumeGeometry->getPixelLengthX() == m_pVolumeGeometry->getPixelLengthY(), "ParallelBeamLinearKernelProjector2D", "Pixel height must equal pixel width."); diff --git a/src/ParallelBeamStripKernelProjector2D.cpp b/src/ParallelBeamStripKernelProjector2D.cpp index 44c6fec..bb693a2 100644 --- a/src/ParallelBeamStripKernelProjector2D.cpp +++ b/src/ParallelBeamStripKernelProjector2D.cpp @@ -88,7 +88,7 @@ bool CParallelBeamStripKernelProjector2D::_check() // check base class ASTRA_CONFIG_CHECK(CProjector2D::_check(), "ParallelBeamStripKernelProjector2D", "Error in Projector2D initialization"); - ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamStripKernelProjector2D", "Unsupported projection geometry"); + ASTRA_CONFIG_CHECK(dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry) || dynamic_cast<CParallelVecProjectionGeometry2D*>(m_pProjectionGeometry), "ParallelBeamStripKernelProjector2D", "Unsupported projection geometry"); ASTRA_CONFIG_CHECK(m_pVolumeGeometry->getPixelLengthX() == m_pVolumeGeometry->getPixelLengthY(), "ParallelBeamStripKernelProjector2D", "Pixel height must equal pixel width."); @@ -164,57 +164,57 @@ void CParallelBeamStripKernelProjector2D::computeSingleRayWeights(int _iProjecti // Splat a single point std::vector<SDetector2D> CParallelBeamStripKernelProjector2D::projectPoint(int _iRow, int _iCol) { - float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; + // float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; + // float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; std::vector<SDetector2D> res; - // loop projectors and detectors - for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // get projection angle - float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - if (theta >= 7*PIdiv4) theta -= 2*PI; - bool inverse = false; - if (theta >= 3*PIdiv4) { - theta -= PI; - inverse = true; - } - - // calculate distance from the center of the voxel to the ray though the origin - float32 tUL = xUL * cos(theta) + yUL * sin(theta); - float32 tUR = xUR * cos(theta) + yUR * sin(theta); - float32 tLL = xLL * cos(theta) + yLL * sin(theta); - float32 tLR = xLR * cos(theta) + yLR * sin(theta); - if (inverse) { - tUL *= -1.0f; - tUR *= -1.0f; - tLL *= -1.0f; - tLR *= -1.0f; - } - float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); - float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); - - // calculate the offset on the detectorarray (in indices) - int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); - int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); - - // add affected detectors to the list - for (int i = dmin; i <= dmax; ++i) { - if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { - SDetector2D det; - det.m_iAngleIndex = iProjection; - det.m_iDetectorIndex = i; - det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - res.push_back(det); - } - } - } + // // loop projectors and detectors + // for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { + + // // get projection angle + // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); + // if (theta >= 7*PIdiv4) theta -= 2*PI; + // bool inverse = false; + // if (theta >= 3*PIdiv4) { + // theta -= PI; + // inverse = true; + // } + + // // calculate distance from the center of the voxel to the ray though the origin + // float32 tUL = xUL * cos(theta) + yUL * sin(theta); + // float32 tUR = xUR * cos(theta) + yUR * sin(theta); + // float32 tLL = xLL * cos(theta) + yLL * sin(theta); + // float32 tLR = xLR * cos(theta) + yLR * sin(theta); + // if (inverse) { + // tUL *= -1.0f; + // tUR *= -1.0f; + // tLL *= -1.0f; + // tLR *= -1.0f; + // } + // float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); + // float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); + + // // calculate the offset on the detectorarray (in indices) + // int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); + // int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); + + // // add affected detectors to the list + // for (int i = dmin; i <= dmax; ++i) { + // if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { + // SDetector2D det; + // det.m_iAngleIndex = iProjection; + // det.m_iDetectorIndex = i; + // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; + // res.push_back(det); + // } + // } + // } // return result vector return res; diff --git a/src/ParallelProjectionGeometry2D.cpp b/src/ParallelProjectionGeometry2D.cpp index 699e141..8632b17 100644 --- a/src/ParallelProjectionGeometry2D.cpp +++ b/src/ParallelProjectionGeometry2D.cpp @@ -67,8 +67,7 @@ CParallelProjectionGeometry2D::CParallelProjectionGeometry2D(const CParallelProj initialize(_projGeom.m_iProjectionAngleCount, _projGeom.m_iDetectorCount, _projGeom.m_fDetectorWidth, - _projGeom.m_pfProjectionAngles, - _projGeom.m_pfExtraDetectorOffset); + _projGeom.m_pfProjectionAngles); } //---------------------------------------------------------------------------------------- @@ -182,8 +181,8 @@ Config* CParallelProjectionGeometry2D::getConfiguration() const cfg->self.addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); return cfg; } -//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- CVector3D CParallelProjectionGeometry2D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex /* = 0 */) { CVector3D vOutput; @@ -197,4 +196,24 @@ CVector3D CParallelProjectionGeometry2D::getProjectionDirection(int _iProjection return vOutput; } +//---------------------------------------------------------------------------------------- +CParallelVecProjectionGeometry2D* CParallelProjectionGeometry2D::toVectorGeometry() +{ + SParProjection* vectors = new SParProjection[m_iProjectionAngleCount]; + for (int i = 0; i < m_iProjectionAngleCount; ++i) + { + vectors[i].fRayX = sinf(m_pfProjectionAngles[i]); + vectors[i].fRayY = -cosf(m_pfProjectionAngles[i]); + vectors[i].fDetUX = cosf(m_pfProjectionAngles[i]) * m_fDetectorWidth; + vectors[i].fDetUY = sinf(m_pfProjectionAngles[i]) * m_fDetectorWidth; + vectors[i].fDetSX = -0.5f * m_iDetectorCount * vectors[i].fDetUX; + vectors[i].fDetSY = -0.5f * m_iDetectorCount * vectors[i].fDetUY; + } + + CParallelVecProjectionGeometry2D* vecGeom = new CParallelVecProjectionGeometry2D(); + vecGeom->initialize(m_iProjectionAngleCount, m_iDetectorCount, vectors); + delete[] vectors; + return vecGeom; +} + } // end namespace astra diff --git a/src/ParallelVecProjectionGeometry2D.cpp b/src/ParallelVecProjectionGeometry2D.cpp new file mode 100644 index 0000000..c1ed1c4 --- /dev/null +++ b/src/ParallelVecProjectionGeometry2D.cpp @@ -0,0 +1,239 @@ +/* +----------------------------------------------------------------------- +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$ +*/ + +#include "astra/ParallelVecProjectionGeometry2D.h" + +#include <cstring> +#include <sstream> +#include <boost/lexical_cast.hpp> + +using namespace std; + +namespace astra +{ + +//---------------------------------------------------------------------------------------- +// Default constructor. Sets all variables to zero. +CParallelVecProjectionGeometry2D::CParallelVecProjectionGeometry2D() +{ + _clear(); + m_pProjectionAngles = 0; +} + +//---------------------------------------------------------------------------------------- +// Constructor. +CParallelVecProjectionGeometry2D::CParallelVecProjectionGeometry2D(int _iProjectionAngleCount, + int _iDetectorCount, + const SParProjection* _pProjectionAngles) +{ + this->initialize(_iProjectionAngleCount, + _iDetectorCount, + _pProjectionAngles); +} + +//---------------------------------------------------------------------------------------- +// Copy Constructor +CParallelVecProjectionGeometry2D::CParallelVecProjectionGeometry2D(const CParallelVecProjectionGeometry2D& _projGeom) +{ + _clear(); + this->initialize(_projGeom.m_iProjectionAngleCount, + _projGeom.m_iDetectorCount, + _projGeom.m_pProjectionAngles); +} + +//---------------------------------------------------------------------------------------- +// Destructor. +CParallelVecProjectionGeometry2D::~CParallelVecProjectionGeometry2D() +{ + // TODO + delete[] m_pProjectionAngles; +} + + +//---------------------------------------------------------------------------------------- +// Initialization. +bool CParallelVecProjectionGeometry2D::initialize(int _iProjectionAngleCount, + int _iDetectorCount, + const SParProjection* _pProjectionAngles) +{ + m_iProjectionAngleCount = _iProjectionAngleCount; + m_iDetectorCount = _iDetectorCount; + m_pProjectionAngles = new SParProjection[m_iProjectionAngleCount]; + for (int i = 0; i < m_iProjectionAngleCount; ++i) + m_pProjectionAngles[i] = _pProjectionAngles[i]; + + // TODO: check? + + // success + m_bInitialized = _check(); + return m_bInitialized; +} + +//---------------------------------------------------------------------------------------- +// Initialization with a Config object +bool CParallelVecProjectionGeometry2D::initialize(const Config& _cfg) +{ + ASTRA_ASSERT(_cfg.self); + ConfigStackCheck<CProjectionGeometry2D> CC("ParallelVecProjectionGeometry2D", this, _cfg); + + XMLNode* node; + + // TODO: Fix up class hierarchy... this class doesn't fit very well. + // initialization of parent class + //CProjectionGeometry2D::initialize(_cfg); + + // Required: DetectorCount + node = _cfg.self->getSingleNode("DetectorCount"); + ASTRA_CONFIG_CHECK(node, "ParallelVecProjectionGeometry2D", "No DetectorRowCount tag specified."); + m_iDetectorCount = boost::lexical_cast<int>(node->getContent()); + ASTRA_DELETE(node); + CC.markNodeParsed("DetectorCount"); + + // Required: Vectors + node = _cfg.self->getSingleNode("Vectors"); + ASTRA_CONFIG_CHECK(node, "ParallelVecProjectionGeometry2D", "No Vectors tag specified."); + vector<float32> data = node->getContentNumericalArray(); + CC.markNodeParsed("Vectors"); + ASTRA_DELETE(node); + ASTRA_CONFIG_CHECK(data.size() % 6 == 0, "ParallelVecProjectionGeometry2D", "Vectors doesn't consist of 6-tuples."); + m_iProjectionAngleCount = data.size() / 6; + m_pProjectionAngles = new SParProjection[m_iProjectionAngleCount]; + + for (int i = 0; i < m_iProjectionAngleCount; ++i) { + SParProjection& p = m_pProjectionAngles[i]; + p.fRayX = data[6*i + 0]; + p.fRayY = data[6*i + 1]; + p.fDetUX = data[6*i + 4]; + p.fDetUY = data[6*i + 5]; + + // The backend code currently expects the corner of the detector, while + // the matlab interface supplies the center + p.fDetSX = data[6*i + 2] - 0.5f * m_iDetectorCount * p.fDetUX; + p.fDetSY = data[6*i + 3] - 0.5f * m_iDetectorCount * p.fDetUY; + } + + + + // success + m_bInitialized = _check(); + return m_bInitialized; +} + +//---------------------------------------------------------------------------------------- +// Clone +CProjectionGeometry2D* CParallelVecProjectionGeometry2D::clone() +{ + return new CParallelVecProjectionGeometry2D(*this); +} + +//---------------------------------------------------------------------------------------- +// is equal +bool CParallelVecProjectionGeometry2D::isEqual(CProjectionGeometry2D* _pGeom2) const +{ + if (_pGeom2 == NULL) return false; + + // try to cast argument to CParallelVecProjectionGeometry2D + CParallelVecProjectionGeometry2D* pGeom2 = dynamic_cast<CParallelVecProjectionGeometry2D*>(_pGeom2); + if (pGeom2 == NULL) return false; + + // both objects must be initialized + if (!m_bInitialized || !pGeom2->m_bInitialized) return false; + + // check all values + if (m_iProjectionAngleCount != pGeom2->m_iProjectionAngleCount) return false; + if (m_iDetectorCount != pGeom2->m_iDetectorCount) return false; + + for (int i = 0; i < m_iProjectionAngleCount; ++i) { + if (memcmp(&m_pProjectionAngles[i], &pGeom2->m_pProjectionAngles[i], sizeof(m_pProjectionAngles[i])) != 0) return false; + } + + return true; +} + +//---------------------------------------------------------------------------------------- +// Is of type +bool CParallelVecProjectionGeometry2D::isOfType(const std::string& _sType) +{ + return (_sType == "parallel_vec"); +} + +//---------------------------------------------------------------------------------------- + +CVector3D CParallelVecProjectionGeometry2D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex /* = 0 */) +{ + CVector3D vOutput(0.0f, 0.0f, 0.0f); + + // not implemented + ASTRA_ASSERT(false); + + return vOutput; +} + +//---------------------------------------------------------------------------------------- + +void CParallelVecProjectionGeometry2D::getRayParams(int _iRow, int _iColumn, float32& _fT, float32& _fTheta) const +{ + // not implemented + ASTRA_ASSERT(false); +} + +//---------------------------------------------------------------------------------------- + +bool CParallelVecProjectionGeometry2D::_check() +{ + // TODO + return true; +} + + +//---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CParallelVecProjectionGeometry2D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry2D"); + cfg->self->addAttribute("type", "parallel_vec"); + cfg->self->addChildNode("DetectorCount", getDetectorCount()); + std::string vectors = ""; + for (int i = 0; i < m_iProjectionAngleCount; ++i) { + SParProjection& p = m_pProjectionAngles[i]; + vectors += boost::lexical_cast<string>(p.fRayX) + ","; + vectors += boost::lexical_cast<string>(p.fRayY) + ","; + vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f * m_iDetectorCount * p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f * m_iDetectorCount * p.fDetUY) + ","; + vectors += boost::lexical_cast<string>(p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetUY); + if (i < m_iProjectionAngleCount-1) vectors += ';'; + } + cfg->self->addChildNode("Vectors", vectors); + return cfg; +} +//---------------------------------------------------------------------------------------- + + +} // namespace astra diff --git a/src/ProjectionGeometry2D.cpp b/src/ProjectionGeometry2D.cpp index b89605b..da93b9c 100644 --- a/src/ProjectionGeometry2D.cpp +++ b/src/ProjectionGeometry2D.cpp @@ -72,7 +72,6 @@ void CProjectionGeometry2D::_clear() m_iDetectorCount = 0; m_fDetectorWidth = 0.0f; m_pfProjectionAngles = NULL; - m_pfExtraDetectorOffset = NULL; m_bInitialized = false; } @@ -85,10 +84,8 @@ void CProjectionGeometry2D::clear() m_fDetectorWidth = 0.0f; if (m_bInitialized){ delete[] m_pfProjectionAngles; - delete[] m_pfExtraDetectorOffset; } m_pfProjectionAngles = NULL; - m_pfExtraDetectorOffset = NULL; m_bInitialized = false; } @@ -147,19 +144,6 @@ bool CProjectionGeometry2D::initialize(const Config& _cfg) } CC.markNodeParsed("ProjectionAngles"); - vector<float32> offset = _cfg.self.getOptionNumericalArray("ExtraDetectorOffset"); - m_pfExtraDetectorOffset = new float32[m_iProjectionAngleCount]; - if (offset.size() == (size_t)m_iProjectionAngleCount) { - for (int i = 0; i < m_iProjectionAngleCount; i++) { - m_pfExtraDetectorOffset[i] = offset[i]; - } - } else { - for (int i = 0; i < m_iProjectionAngleCount; i++) { - m_pfExtraDetectorOffset[i] = 0.0f; - } - } - CC.markOptionParsed("ExtraDetectorOffset"); - // some checks ASTRA_CONFIG_CHECK(m_iDetectorCount > 0, "ProjectionGeometry2D", "DetectorCount should be positive."); ASTRA_CONFIG_CHECK(m_fDetectorWidth > 0.0f, "ProjectionGeometry2D", "DetectorWidth should be positive."); @@ -186,10 +170,8 @@ bool CProjectionGeometry2D::_initialize(int _iProjectionAngleCount, m_iDetectorCount = _iDetectorCount; m_fDetectorWidth = _fDetectorWidth; m_pfProjectionAngles = new float32[m_iProjectionAngleCount]; - m_pfExtraDetectorOffset = new float32[m_iProjectionAngleCount]; for (int i = 0; i < m_iProjectionAngleCount; i++) { m_pfProjectionAngles[i] = _pfProjectionAngles[i]; - m_pfExtraDetectorOffset[i] = _pfExtraDetectorOffsets ? _pfExtraDetectorOffsets[i]:0; } // Interface class, so don't set m_bInitialized to true diff --git a/src/Projector2D.cpp b/src/Projector2D.cpp index cf233a0..c8743b5 100644 --- a/src/Projector2D.cpp +++ b/src/Projector2D.cpp @@ -28,6 +28,7 @@ $Id$ #include "astra/Projector2D.h" +#include "astra/ParallelVecProjectionGeometry2D.h" #include "astra/FanFlatProjectionGeometry2D.h" #include "astra/FanFlatVecProjectionGeometry2D.h" #include "astra/SparseMatrixProjectionGeometry2D.h" @@ -131,6 +132,10 @@ bool CProjector2D::initialize(const Config& _cfg) CFanFlatVecProjectionGeometry2D* pFanFlatVecProjectionGeometry = new CFanFlatVecProjectionGeometry2D(); pFanFlatVecProjectionGeometry->initialize(Config(node)); m_pProjectionGeometry = pFanFlatVecProjectionGeometry; + } else if (type == "parallel_vec") { + CParallelVecProjectionGeometry2D* pParallelVecProjectionGeometry = new CParallelVecProjectionGeometry2D(); + pParallelVecProjectionGeometry->initialize(Config(node)); + m_pProjectionGeometry = pParallelVecProjectionGeometry; } else { m_pProjectionGeometry = new CParallelProjectionGeometry2D(); m_pProjectionGeometry->initialize(Config(node)); |