+ p = pArgBase+28; SwapBytes4IfLittleEndian (p);
+ // long lTime = *reinterpret_cast<long*>(p);
+
+ p = pArgBase + 4; SwapBytes4IfLittleEndian (p);
+ double dAlpha = *reinterpret_cast<float*>(p) + HALFPI;
+
+ p = pArgBase+12; SwapBytes4IfLittleEndian (p);
+ // double dAlign = *reinterpret_cast<float*>(p);
+
+ p = pArgBase + 16; SwapBytes4IfLittleEndian (p);
+ // double dMaxValue = *reinterpret_cast<float*>(p);
+
+ DetectorArray& detArray = getDetectorArray (iv);
+ detArray.setViewAngle (dAlpha);
+ DetectorValue* detval = detArray.detValues();
+
+ double dViewScale = 1. / (2294.4871 * ::pow (2.0, -lEscale));
+ lDataPos += 32;
+ for (int id = 0; id < iNDets; id++) {
+ int iV = pData[lDataPos+1] + (pData[lDataPos] << 8);
+ if (iV > 32767) // two's complement signed conversion
+ iV = iV - 65536;
+ detval[id] = iV * dViewScale * pdCosScale[id];
+ lDataPos += 2;
+ }
+#if 1
+ for (int k = iNDets - 2; k >= 0; k--)
+ detval[k+1] = detval[k];
+ detval[0] = 0;
+#endif
+ }
+
+ delete pdCosScale;
+ return true;
+}
+
+Projections*
+Projections::interpolateToParallel () const
+{
+ if (m_geometry == Scanner::GEOMETRY_PARALLEL)
+ return const_cast<Projections*>(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 (isEven (nDet)) // even
+ pProjNew->m_detInc = m_dViewDiameter / (nDet - 1);
+
+ ParallelRaysums parallel (this, ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI);
+
+ 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;
+ int iLastFloor = -1;
+ LinearInterpolator<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false);
+ for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
+ DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
+ detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor);
+ }
+ }
+ 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;
+ int iLastFloor = -1;
+ LinearInterpolator<double> interp (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), false);
+ for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc)
+ detValues[iD] = interp.interpolate (dDetPos, &iLastFloor);
+ }
+ delete pdDetValueCopy;
+ delete pdOriginalDetPositions;
+
+ return pProjNew;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Class ParallelRaysums
+//
+// Used for converting divergent beam raysums into Parallel raysums
+//
+///////////////////////////////////////////////////////////////////////////////
+
+ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange)
+: m_pCoordinates(NULL), m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
+ m_iThetaRange (iThetaRange)
+{
+ int iGeometry = pProjections->geometry();
+ double dDetInc = pProjections->detInc();
+ double dDetStart = pProjections->detStart();
+ double dFocalLength = pProjections->focalLength();
+
+ m_iNumCoordinates = m_iNumView * m_iNumDet;
+ m_pCoordinates = new ParallelRaysumCoordinate [m_iNumCoordinates];
+ m_vecpCoordinates.reserve (m_iNumCoordinates);
+ for (int i = 0; i < m_iNumCoordinates; i++)
+ m_vecpCoordinates[i] = m_pCoordinates + i;
+
+ int iCoordinate = 0;
+ 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 < m_iNumDet; iD++) {
+ ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++];
+
+ if (iGeometry == Scanner::GEOMETRY_PARALLEL) {
+ pC->m_dTheta = dViewAngle;
+ pC->m_dT = dDetPos;
+ } else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+ double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength());
+ pC->m_dTheta = dViewAngle + dFanAngle;
+ pC->m_dT = dFocalLength * sin(dFanAngle);
+
+ } else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+ // fan angle is same as dDetPos
+ pC->m_dTheta = dViewAngle + dDetPos;
+ pC->m_dT = dFocalLength * sin (dDetPos);
+ }
+ if (m_iThetaRange != THETA_RANGE_UNCONSTRAINED) {
+ pC->m_dTheta = normalizeAngle (pC->m_dTheta);
+ if (m_iThetaRange == THETA_RANGE_FOLD_TO_PI && pC->m_dTheta >= PI) {
+ pC->m_dTheta -= PI;
+ pC->m_dT = -pC->m_dT;
+ }
+ }
+ pC->m_dRaysum = detValues[iD];
+ dDetPos += dDetInc;
+ }
+ }
+}
+
+ParallelRaysums::~ParallelRaysums()
+{
+ delete m_pCoordinates;
+}
+
+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);
+ }
+
+ 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_vecpCoordinates[0]->m_dT;
+ *dMinTheta = *dMaxTheta = m_vecpCoordinates[0]->m_dTheta;
+
+ for (int i = 0; i < m_iNumCoordinates; i++) {
+ double dT = m_vecpCoordinates[i]->m_dT;
+ double dTheta = m_vecpCoordinates[i]->m_dTheta;
+
+ if (dT < *dMinT)
+ *dMinT = dT;
+ else if (dT > *dMaxT)
+ *dMaxT = dT;
+
+ if (dTheta < *dMinTheta)
+ *dMinTheta = dTheta;
+ else if (dTheta > *dMaxTheta)
+ *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;
+ }
+}