** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.55 2001/03/10 23:14:16 kevin Exp $
+** $Id: projections.cpp,v 1.59 2001/03/11 15:27:30 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();
}
}
}
-
bool
Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pData, unsigned long lDataLength)
{
init (iNViews, iNDets);
m_geometry = Scanner::GEOMETRY_EQUIANGULAR;
m_dFanBeamAngle = iNDets * convertDegreesToRadians (3.06976 / 60);
- m_dFocalLength = 51;
- m_dSourceDetectorLength = 89;
+ m_dFocalLength = 510;
+ m_dSourceDetectorLength = 890;
m_detInc = convertDegreesToRadians (3.06976 / 60);
m_detStart = -m_dFanBeamAngle / 2;
m_rotInc = TWOPI / static_cast<double>(iNViews);
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;
+#ifdef CONVERT_PARALLEL_PI
+ pProjNew->m_rotInc = PI / nView;;
+#else
+ pProjNew->m_rotInc = TWOPI / nView;
+#endif
+ 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_ppCoordinates(NULL), 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_ppCoordinates = new ParallelRaysumCoordinate* [m_iNumCoordinates];
+ m_iNumCoordinates = m_iNumView * m_iNumDet;
+ m_vecpCoordinates.reserve (m_iNumCoordinates);
for (int i = 0; i < m_iNumCoordinates; i++)
- m_ppCoordinates[i] = new ParallelRaysumCoordinate;
-
- ParallelRaysumCoordinate** ppCurrentCoordinate = m_ppCoordinates;
+ m_vecpCoordinates[i] = new ParallelRaysumCoordinate;
- for (int iV = 0; iV < nView; iV++) {
+ int iCoordinate = 0;
+ for (int iV = 0; iV < m_iNumView; iV++) {
double dViewAngle = pProjections->getDetectorArray(iV).viewAngle();
- for (int iD = 0; iD < nDet; iD++) {
- ParallelRaysumCoordinate* pC = *ppCurrentCoordinate;
+ const DetectorValue* detValues = pProjections->getDetectorArray(iV).detValues();
+
+ double dDetPos = dDetStart;
+ for (int iD = 0; iD < m_iNumDet; iD++) {
+ ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++];
+
if (iGeometry == Scanner::GEOMETRY_PARALLEL) {
- pC->m_dTheta = dViewAngle;
- pC->m_dT = dDetStart + (iD * dDetInc);
+ pC->m_dTheta = normalizeAngle (dViewAngle);
+ pC->m_dT = dDetPos;
} else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) {
- double dDetPos = dDetStart + (iD * dDetInc);
double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength());
- pC->m_dTheta = dViewAngle + dFanAngle;
+ pC->m_dTheta = normalizeAngle (dViewAngle + dFanAngle);
pC->m_dT = dFocalLength * sin(dFanAngle);
+
} else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
- double dFanAngle = dDetStart + (iD * dDetInc);
- pC->m_dTheta = dViewAngle + dFanAngle;
- pC->m_dT = dFocalLength * sin(dFanAngle);
+ // fan angle is same as dDetPos
+ pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos);
+ pC->m_dT = dFocalLength * sin (dDetPos);
+ }
+#ifdef CONVERT_PARALLEL_PI
+ if (pC->m_dTheta >= PI) { // convert T/Theta to 0-PI interval
+ pC->m_dTheta -= PI;
+ pC->m_dT = -pC->m_dT;
}
- ++ppCurrentCoordinate;
+#endif
+ pC->m_dRaysum = detValues[iD];
+ dDetPos += dDetInc;
}
}
}
ParallelRaysums::~ParallelRaysums()
{
for (int i = 0; i < m_iNumCoordinates; i++)
- delete m_ppCoordinates[i];
+ delete m_vecpCoordinates[i];
+}
+
+ParallelRaysums::CoordinateContainer&
+ParallelRaysums::getSortedByTheta()
+{
+ if (m_vecpSortedByTheta.size() == 0) {
+ 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);
+ }
- delete m_ppCoordinates;
+ return m_vecpSortedByTheta;
}
+ParallelRaysums::CoordinateContainer&
+ParallelRaysums::getSortedByT()
+{
+ if (m_vecpSortedByT.size() == 0) {
+ 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);
+ }
+
+ return m_vecpSortedByT;
+}
+
+
void
ParallelRaysums::getLimits (double* dMinT, double* dMaxT, double* dMinTheta, double* dMaxTheta) const
{
if (m_iNumCoordinates <= 0)
return;
- *dMinT = *dMaxT = m_ppCoordinates[0]->m_dT;
- *dMinTheta = *dMaxTheta = m_ppCoordinates[0]->m_dTheta;
+ *dMinT = *dMaxT = m_vecpCoordinates[0]->m_dT;
+ *dMinTheta = *dMaxTheta = m_vecpCoordinates[0]->m_dTheta;
for (int i = 0; i < m_iNumCoordinates; i++) {
- double dT = m_ppCoordinates[i]->m_dT;
- double dTheta = m_ppCoordinates[i]->m_dTheta;
+ double dT = m_vecpCoordinates[i]->m_dT;
+ double dTheta = m_vecpCoordinates[i]->m_dTheta;
+
if (dT < *dMinT)
*dMinT = dT;
else if (dT > *dMaxT)
*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