X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=27d5b69ffa7516d5435ce0406068b61fd94bbeea;hb=08c1ec110dc7936c2bbd1c619bd2cf3618c6b4cc;hp=b3d6d35b2df647b5595d60c01bebcf30771af957;hpb=b361677a2f7d5b443641faec70b057f2a84dc77e;p=ctsim.git diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index b3d6d35..27d5b69 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** 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.60 2001/03/11 17:55:29 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 @@ -150,12 +150,12 @@ Projections::initFromScanner (const Scanner& scanner) 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(); } @@ -688,17 +688,16 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) if (! v || nx == 0 || ny == 0) return false; - if (m_geometry != Scanner::GEOMETRY_PARALLEL) { - sys_error (ERR_WARNING, "convertPolar supports Parallel only"); - return false; - } + Projections* pProj = this; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + pProj = interpolateToParallel(); Array2d adView (nx, ny); Array2d adDet (nx, ny); double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - if (! calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet)) + if (! pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet)) return false; std::complex** ppcDetValue = new std::complex* [m_nView]; @@ -709,12 +708,15 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) ppcDetValue[iView][iDet] = std::complex(getDetectorArray (iView).detValues()[iDet], 0); } - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); for (iView = 0; iView < m_nView; iView++) delete [] ppcDetValue[iView]; delete [] ppcDetValue; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + delete pProj; + return true; } @@ -885,15 +887,14 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, } } - 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(iNViews); @@ -968,42 +969,136 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa 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); } - ++ppCurrentCoordinate; +#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; + } +#endif + pC->m_dRaysum = detValues[iD]; + dDetPos += dDetInc; } } } @@ -1011,23 +1106,49 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) ParallelRaysums::~ParallelRaysums() { for (int i = 0; i < m_iNumCoordinates; i++) - delete m_ppCoordinates[i]; + delete m_vecpCoordinates[i]; +} - delete m_ppCoordinates; +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_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) @@ -1039,3 +1160,56 @@ ParallelRaysums::getLimits (double* dMinT, double* dMaxT, double* dMinTheta, dou *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])); +} +