X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=60f9467e06cbb2426ed2db355b9c163248ffdbd7;hb=9ff5b5165b2c8871bd4b29ccd5ca794638414615;hp=dc7248e62b9f3e42c86561612e99233c1bd6c245;hpb=de411914da8b157958e9caae917bf1edeafbb713;p=ctsim.git diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index dc7248e..60f9467 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.59 2001/03/11 15:27:30 kevin Exp $ +** $Id: projections.cpp,v 1.62 2001/03/13 04:44:25 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 @@ -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]; @@ -706,15 +705,18 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) for (iView = 0; iView < m_nView; iView++) { ppcDetValue[iView] = new std::complex [m_nDet]; for (unsigned int iDet = 0; iDet < m_nDet; iDet++) - ppcDetValue[iView][iDet] = std::complex(getDetectorArray (iView).detValues()[iDet], 0); + ppcDetValue[iView][iDet] = std::complex(pProj->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, pProj->m_nView, pProj->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; } @@ -890,80 +892,67 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa { init (iNViews, iNDets); m_geometry = Scanner::GEOMETRY_EQUIANGULAR; - m_dFanBeamAngle = iNDets * convertDegreesToRadians (3.06976 / 60); m_dFocalLength = 510; m_dSourceDetectorLength = 890; m_detInc = convertDegreesToRadians (3.06976 / 60); - m_detStart = -m_dFanBeamAngle / 2; + m_detStart = -(m_dFanBeamAngle / 2); m_rotInc = TWOPI / static_cast(iNViews); m_rotStart = HALFPI; + m_dFanBeamAngle = (iNDets + 1) * m_detInc; m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2; - if (iNDets != 1024) - return false; - bool bValid = (iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) || (iNViews == 1500 && lDataLength == 3120000); - if (! bValid) + if (! ((iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) + || (iNViews == 1500 && lDataLength == 3120000))) return false; + int iCenter = (iNDets - 1) / 2; // change from (Nm+1)/2 because of 0 vs. 1 indexing + double* pdCosScale = new double [iNDets]; + for (int i = 0; i < iNDets; i++) + pdCosScale[i] = cos ((i - iCenter) * m_detInc); + long lDataPos = 0; for (int iv = 0; iv < iNViews; iv++) { unsigned char* pArgBase = pData + lDataPos; unsigned char* p = pArgBase+0; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + SwapBytes4IfLittleEndian (p); long lProjNumber = *reinterpret_cast(p); p = pArgBase+20; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + SwapBytes4IfLittleEndian (p); long lEscale = *reinterpret_cast(p); p = pArgBase+28; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + SwapBytes4IfLittleEndian (p); long lTime = *reinterpret_cast(p); p = pArgBase + 4; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + SwapBytes4IfLittleEndian (p); double dAlpha = *reinterpret_cast(p) + HALFPI; p = pArgBase+12; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + SwapBytes4IfLittleEndian (p); double dAlign = *reinterpret_cast(p); p = pArgBase + 16; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + SwapBytes4IfLittleEndian (p); double dMaxValue = *reinterpret_cast(p); - lDataPos += 32; - double dEScale = pow (2.0, -lEscale); - double dBetaInc = convertDegreesToRadians (3.06976 / 60); - int iCenter = (iNDets + 1) / 2; - - DetectorArray& detArray = getDetectorArray( iv ); + DetectorArray& detArray = getDetectorArray (iv); detArray.setViewAngle (dAlpha); DetectorValue* detval = detArray.detValues(); - double dTempScale = 2294.4871 * dEScale; + double dViewScale = 2294.4871 * ::pow (2.0, -lEscale); + lDataPos += 32; for (int id = 0; id < iNDets; id++) { int iV = pData[lDataPos+1] + 256 * pData[lDataPos]; if (iV > 32767) // two's complement signed conversion iV = iV - 65536; - double dCosScale = cos ((id + 1 - iCenter) * dBetaInc); - detval[id] = iV / (dTempScale * dCosScale); + detval[id] = iV / (dViewScale * pdCosScale[id]); lDataPos += 2; } } + delete pdCosScale; return true; } @@ -1000,7 +989,7 @@ Projections::interpolateToParallel () if (nDet % 2 == 0) // even pProjNew->m_detInc = m_dViewDiameter / (nDet - 1); - ParallelRaysums parallel (this); + ParallelRaysums parallel (this, ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI); double* pdThetaValuesForT = new double [pProjNew->nView()]; double* pdRaysumsForT = new double [pProjNew->nView()]; @@ -1011,10 +1000,11 @@ Projections::interpolateToParallel () parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT); double dViewAngle = m_rotStart; + int iLastFloor = -1; 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); + detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle, &iLastFloor); } } delete pdThetaValuesForT; @@ -1035,8 +1025,9 @@ Projections::interpolateToParallel () pdDetValueCopy[i] = detValues[i]; double dDetPos = pProjNew->m_detStart; + int iLastFloor = -1; for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) { - detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos); + detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos, &iLastFloor); } } delete pdDetValueCopy; @@ -1054,8 +1045,9 @@ Projections::interpolateToParallel () // /////////////////////////////////////////////////////////////////////////////// -ParallelRaysums::ParallelRaysums (Projections* pProjections) -: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()) +ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange) +: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()), + m_iThetaRange (iThetaRange) { int iGeometry = pProjections->geometry(); double dDetInc = pProjections->detInc(); @@ -1077,24 +1069,25 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++]; if (iGeometry == Scanner::GEOMETRY_PARALLEL) { - pC->m_dTheta = normalizeAngle (dViewAngle); + pC->m_dTheta = 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 = dViewAngle + dFanAngle; pC->m_dT = dFocalLength * sin(dFanAngle); } else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) { // fan angle is same as dDetPos - pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos); + pC->m_dTheta = 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; + 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; + } } -#endif pC->m_dRaysum = detValues[iD]; dDetPos += dDetInc; } @@ -1185,11 +1178,14 @@ ParallelRaysums::getDetPositions (double* pdDetPos) } // locate by bisection, O(log2(n)) +// iLastFloor is used when sequential calls to interpolate have monotonically increasing dX double -ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX) +ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX, int* iLastFloor) { int iLower = -1; int iUpper = n; + if (iLastFloor && *iLastFloor >= 0 && pdX[*iLastFloor] < dX) + iLower = *iLastFloor; while (iUpper - iLower > 1) { int iMiddle = (iUpper + iLower) >> 1; @@ -1208,6 +1204,8 @@ ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX) return 0; } + if (iLastFloor) + *iLastFloor = iLower; return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower])); }