X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=60f9467e06cbb2426ed2db355b9c163248ffdbd7;hp=2d2592e62f8d9382d121d332c29b2cfaa9e42089;hb=9ff5b5165b2c8871bd4b29ccd5ca794638414615;hpb=59f85941934e85d4824c4d76e4d7ed3375bf2e3a diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 2d2592e..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.61 2001/03/11 18:52:03 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 @@ -892,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; } @@ -1002,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()]; @@ -1013,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; @@ -1037,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; @@ -1056,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(); @@ -1079,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; } @@ -1187,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; @@ -1210,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])); }