diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/ConeVecProjectionGeometry3D.cpp | 8 | ||||
| -rw-r--r-- | src/FanFlatProjectionGeometry2D.cpp | 18 | ||||
| -rw-r--r-- | src/ParallelBeamBlobKernelProjector2D.cpp | 92 | ||||
| -rw-r--r-- | src/ParallelBeamLineKernelProjector2D.cpp | 108 | ||||
| -rw-r--r-- | src/ParallelBeamLinearKernelProjector2D.cpp | 2 | ||||
| -rw-r--r-- | src/ParallelBeamStripKernelProjector2D.cpp | 100 | ||||
| -rw-r--r-- | src/ParallelProjectionGeometry2D.cpp | 25 | ||||
| -rw-r--r-- | src/ParallelVecProjectionGeometry2D.cpp | 239 | ||||
| -rw-r--r-- | src/ParallelVecProjectionGeometry3D.cpp | 8 | ||||
| -rw-r--r-- | src/ProjectionGeometry2D.cpp | 18 | ||||
| -rw-r--r-- | src/Projector2D.cpp | 5 |
11 files changed, 439 insertions, 184 deletions
diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp index 9dc273d..a4bd22d 100644 --- a/src/ConeVecProjectionGeometry3D.cpp +++ b/src/ConeVecProjectionGeometry3D.cpp @@ -208,7 +208,7 @@ Config* CConeVecProjectionGeometry3D::getConfiguration() const Config* cfg = new Config(); cfg->initialize("ProjectionGeometry3D"); - cfg->self->addAttribute("type", "cone_vec"); + cfg->self->addAttribute("type", "cone"); cfg->self->addChildNode("DetectorRowCount", m_iDetectorRowCount); cfg->self->addChildNode("DetectorColCount", m_iDetectorColCount); @@ -218,9 +218,9 @@ Config* CConeVecProjectionGeometry3D::getConfiguration() const vectors += boost::lexical_cast<string>(p.fSrcX) + ","; vectors += boost::lexical_cast<string>(p.fSrcY) + ","; vectors += boost::lexical_cast<string>(p.fSrcZ) + ","; - vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f*m_iDetectorRowCount*p.fDetVX + 0.5f*m_iDetectorColCount*p.fDetUX) + ","; - vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f*m_iDetectorRowCount*p.fDetVY + 0.5f*m_iDetectorColCount*p.fDetUY) + ","; - vectors += boost::lexical_cast<string>(p.fDetSZ + 0.5f*m_iDetectorRowCount*p.fDetVZ + 0.5f*m_iDetectorColCount*p.fDetUZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f*m_iDetectorRowCount*p.fDetUX + 0.5f*m_iDetectorColCount*p.fDetVX) + ","; + vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f*m_iDetectorRowCount*p.fDetUY + 0.5f*m_iDetectorColCount*p.fDetVY) + ","; + vectors += boost::lexical_cast<string>(p.fDetSZ + 0.5f*m_iDetectorRowCount*p.fDetUZ + 0.5f*m_iDetectorColCount*p.fDetVZ) + ","; vectors += boost::lexical_cast<string>(p.fDetUX) + ","; vectors += boost::lexical_cast<string>(p.fDetUY) + ","; vectors += boost::lexical_cast<string>(p.fDetUZ) + ","; diff --git a/src/FanFlatProjectionGeometry2D.cpp b/src/FanFlatProjectionGeometry2D.cpp index d757f18..9851d40 100644 --- a/src/FanFlatProjectionGeometry2D.cpp +++ b/src/FanFlatProjectionGeometry2D.cpp @@ -217,7 +217,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 3253f88..d030576 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."); @@ -156,16 +156,6 @@ bool CParallelBeamBlobKernelProjector2D::initialize(const Config& _cfg) 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 +218,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 cac8f30..42ed0a0 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/ParallelVecProjectionGeometry3D.cpp b/src/ParallelVecProjectionGeometry3D.cpp index dc325e9..cfac485 100644 --- a/src/ParallelVecProjectionGeometry3D.cpp +++ b/src/ParallelVecProjectionGeometry3D.cpp @@ -208,7 +208,7 @@ Config* CParallelVecProjectionGeometry3D::getConfiguration() const Config* cfg = new Config(); cfg->initialize("ProjectionGeometry3D"); - cfg->self->addAttribute("type", "parallel3d_vec"); + cfg->self->addAttribute("type", "parallel3d"); cfg->self->addChildNode("DetectorRowCount", m_iDetectorRowCount); cfg->self->addChildNode("DetectorColCount", m_iDetectorColCount); @@ -218,9 +218,9 @@ Config* CParallelVecProjectionGeometry3D::getConfiguration() const vectors += boost::lexical_cast<string>(p.fRayX) + ","; vectors += boost::lexical_cast<string>(p.fRayY) + ","; vectors += boost::lexical_cast<string>(p.fRayZ) + ","; - vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f*m_iDetectorRowCount*p.fDetVX + 0.5f*m_iDetectorColCount*p.fDetUX) + ","; - vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f*m_iDetectorRowCount*p.fDetVY + 0.5f*m_iDetectorColCount*p.fDetUY) + ","; - vectors += boost::lexical_cast<string>(p.fDetSZ + 0.5f*m_iDetectorRowCount*p.fDetVZ + 0.5f*m_iDetectorColCount*p.fDetUZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f*m_iDetectorRowCount*p.fDetUX + 0.5f*m_iDetectorColCount*p.fDetVX) + ","; + vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f*m_iDetectorRowCount*p.fDetUY + 0.5f*m_iDetectorColCount*p.fDetVY) + ","; + vectors += boost::lexical_cast<string>(p.fDetSZ + 0.5f*m_iDetectorRowCount*p.fDetUZ + 0.5f*m_iDetectorColCount*p.fDetVZ) + ","; vectors += boost::lexical_cast<string>(p.fDetUX) + ","; vectors += boost::lexical_cast<string>(p.fDetUY) + ","; vectors += boost::lexical_cast<string>(p.fDetUZ) + ","; diff --git a/src/ProjectionGeometry2D.cpp b/src/ProjectionGeometry2D.cpp index 89b5fe0..a51d86b 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; } @@ -150,19 +147,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."); @@ -189,10 +173,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 32a2956..b0ace92 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)); |