** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.56 2001/03/10 23:56:58 kevin Exp $
+** $Id: projections.cpp,v 1.57 2001/03/11 06:34:37 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
m_rotInc = scanner.rotInc();
m_detInc = scanner.detInc();
+ m_detStart = scanner.detStart();
m_geometry = scanner.geometry();
m_dFocalLength = scanner.focalLength();
m_dSourceDetectorLength = scanner.sourceDetectorLength();
m_dViewDiameter = scanner.viewDiameter();
m_rotStart = 0;
- m_detStart = -(scanner.detLen() / 2);
m_dFanBeamAngle = scanner.fanBeamAngle();
}
return true;
}
+Projections*
+Projections::interpolateToParallel ()
+{
+ if (m_geometry == Scanner::GEOMETRY_PARALLEL)
+ return this;
+
+ int nDet = m_nDet;
+ int nView = m_nView;
+ Projections* pProjNew = new Projections (nView, nDet);
+ pProjNew->m_geometry = Scanner::GEOMETRY_PARALLEL;
+ pProjNew->m_dFocalLength = m_dFocalLength;
+ pProjNew->m_dSourceDetectorLength = m_dSourceDetectorLength;
+ pProjNew->m_dViewDiameter = m_dViewDiameter;
+ pProjNew->m_dFanBeamAngle = m_dFanBeamAngle;
+ pProjNew->m_calcTime = 0;
+ pProjNew->m_remark = m_remark;
+ pProjNew->m_remark += "; Interpolate to Parallel";
+ pProjNew->m_label.setLabelType (Array2dFileLabel::L_HISTORY);
+ pProjNew->m_label.setLabelString (pProjNew->m_remark);
+ pProjNew->m_label.setCalcTime (pProjNew->m_calcTime);
+ pProjNew->m_label.setDateTime (pProjNew->m_year, pProjNew->m_month, pProjNew->m_day, pProjNew->m_hour, pProjNew->m_minute, pProjNew->m_second);
+
+ pProjNew->m_rotStart = 0;
+ pProjNew->m_rotInc = PI / nView;;
+ pProjNew->m_detStart = -m_dViewDiameter / 2;
+ pProjNew->m_detInc = m_dViewDiameter / nDet;
+ if (nDet % 2 == 0) // even
+ pProjNew->m_detInc = m_dViewDiameter / (nDet - 1);
+
+ ParallelRaysums parallel (this);
+
+ double* pdThetaValuesForT = new double [pProjNew->nView()];
+ double* pdRaysumsForT = new double [pProjNew->nView()];
+
+ // interpolate to evenly spaced theta (views)
+ double dDetPos = pProjNew->m_detStart;
+ for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
+ parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT);
+
+ double dViewAngle = m_rotStart;
+ for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
+ DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
+
+ detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle);
+ }
+ }
+ delete pdThetaValuesForT;
+ delete pdRaysumsForT;
+
+ // interpolate to evenly space t (detectors)
+ double* pdOriginalDetPositions = new double [pProjNew->nDet()];
+ parallel.getDetPositions (pdOriginalDetPositions);
+
+ double* pdDetValueCopy = new double [pProjNew->nDet()];
+ double dViewAngle = m_rotStart;
+ for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
+ DetectorArray& detArray = pProjNew->getDetectorArray (iV);
+ DetectorValue* detValues = detArray.detValues();
+ detArray.setViewAngle (dViewAngle);
+
+ for (int i = 0; i < pProjNew->nDet(); i++)
+ pdDetValueCopy[i] = detValues[i];
+
+ double dDetPos = pProjNew->m_detStart;
+ for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
+ detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos);
+ }
+ }
+ delete pdDetValueCopy;
+ delete pdOriginalDetPositions;
+
+ return pProjNew;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Class ParallelRaysums
+//
+// Used for converting divergent beam raysums into Parallel raysums
+//
+///////////////////////////////////////////////////////////////////////////////
ParallelRaysums::ParallelRaysums (Projections* pProjections)
-: m_iNumCoordinates(0)
+: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet())
{
- int nDet = pProjections->nDet();
- int nView = pProjections->nView();
int iGeometry = pProjections->geometry();
double dDetInc = pProjections->detInc();
double dDetStart = pProjections->detStart();
double dFocalLength = pProjections->focalLength();
- m_iNumCoordinates = nDet * nView;
+ m_iNumCoordinates = m_iNumView * m_iNumDet;
m_vecpCoordinates.reserve (m_iNumCoordinates);
for (int i = 0; i < m_iNumCoordinates; i++)
m_vecpCoordinates[i] = new ParallelRaysumCoordinate;
int iCoordinate = 0;
- for (int iV = 0; iV < nView; iV++) {
+ for (int iV = 0; iV < m_iNumView; iV++) {
double dViewAngle = pProjections->getDetectorArray(iV).viewAngle();
+ const DetectorValue* detValues = pProjections->getDetectorArray(iV).detValues();
double dDetPos = dDetStart;
- for (int iD = 0; iD < nDet; iD++) {
+ for (int iD = 0; iD < m_iNumDet; iD++) {
ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++];
if (iGeometry == Scanner::GEOMETRY_PARALLEL) {
pC->m_dTheta = normalizeAngle (dViewAngle);
pC->m_dT = dDetPos;
-
} else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) {
double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength());
pC->m_dTheta = normalizeAngle (dViewAngle + dFanAngle);
pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos);
pC->m_dT = dFocalLength * sin (dDetPos);
}
+ if (pC->m_dTheta >= PI) { // convert T/Theta to 0-PI interval
+ pC->m_dTheta -= PI;
+ pC->m_dT = -pC->m_dT - pProjections->detInc();
+ }
+ pC->m_dRaysum = detValues[iD];
dDetPos += dDetInc;
}
}
ParallelRaysums::getSortedByTheta()
{
if (m_vecpSortedByTheta.size() == 0) {
- m_vecpSortedByTheta.reserve (m_iNumCoordinates);
+ m_vecpSortedByTheta.resize (m_iNumCoordinates);
for (int i = 0; i < m_iNumCoordinates; i++)
m_vecpSortedByTheta[i] = m_vecpCoordinates[i];
std::sort (m_vecpSortedByTheta.begin(), m_vecpSortedByTheta.end(), ParallelRaysumCoordinate::compareByTheta);
ParallelRaysums::getSortedByT()
{
if (m_vecpSortedByT.size() == 0) {
- m_vecpSortedByT.reserve (m_iNumCoordinates);
+ m_vecpSortedByT.resize (m_iNumCoordinates);
for (int i = 0; i < m_iNumCoordinates; i++)
m_vecpSortedByT[i] = m_vecpCoordinates[i];
std::sort (m_vecpSortedByT.begin(), m_vecpSortedByT.end(), ParallelRaysumCoordinate::compareByT);
*dMaxTheta = dTheta;
}
}
+
+void
+ParallelRaysums::getThetaAndRaysumsForT (int iTheta, double* pTheta, double* pRaysum)
+{
+ const CoordinateContainer& coordsT = getSortedByT();
+
+ int iBase = iTheta * m_iNumView;
+ for (int i = 0; i < m_iNumView; i++) {
+ int iPos = iBase + i;
+ pTheta[i] = coordsT[iPos]->m_dTheta;
+ pRaysum[i] = coordsT[iPos]->m_dRaysum;
+ }
+}
+
+void
+ParallelRaysums::getDetPositions (double* pdDetPos)
+{
+ const CoordinateContainer& coordsT = getSortedByT();
+
+ int iPos = 0;
+ for (int i = 0; i < m_iNumDet; i++) {
+ pdDetPos[i] = coordsT[iPos]->m_dT;
+ iPos += m_iNumView;
+ }
+}
+
+// locate by bisection, O(log2(n))
+double
+ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX)
+{
+ int iLower = -1;
+ int iUpper = n;
+
+ while (iUpper - iLower > 1) {
+ int iMiddle = (iUpper + iLower) >> 1;
+ if (dX >= pdX[iMiddle])
+ iLower = iMiddle;
+ else
+ iUpper = iMiddle;
+ }
+ if (dX <= pdX[0])
+ return pdY[0];
+ else if (dX >= pdX[n-1])
+ return pdY[1];
+
+ if (iLower < 0 || iLower >= n) {
+ sys_error (ERR_SEVERE, "Coordinate out of range [locateThetaBase]");
+ return 0;
+ }
+
+ return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower]));
+}
+
projects / ctsim.git / blobdiff