X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=3eab99c4fe9a8a8725bc2108492eb86302e50a64;hp=c55aa37a5fcc94b6292ac9b912f93e858bbb5b7f;hb=4433641931aa27fd6a2b5ecd0102e6c5bbbccc46;hpb=999a754d1519a49ca062ee87b22bf601c1ee9f21 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index c55aa37..3eab99c 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.57 2001/03/11 06:34:37 kevin Exp $ +** $Id: projections.cpp,v 1.63 2001/03/13 08:24:41 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; } @@ -885,94 +887,74 @@ 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_dFanBeamAngle = (iNDets + 1) * m_detInc; + m_detStart = -(m_dFanBeamAngle / 2); m_rotInc = TWOPI / static_cast(iNViews); m_rotStart = HALFPI; 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; + double dCenter = (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] = 1. / cos ((i - dCenter) * 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 + unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p); long lProjNumber = *reinterpret_cast(p); - p = pArgBase+20; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase+20; SwapBytes4IfLittleEndian (p); long lEscale = *reinterpret_cast(p); - p = pArgBase+28; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase+28; SwapBytes4IfLittleEndian (p); long lTime = *reinterpret_cast(p); - p = pArgBase + 4; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase + 4; SwapBytes4IfLittleEndian (p); double dAlpha = *reinterpret_cast(p) + HALFPI; - p = pArgBase+12; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase+12; SwapBytes4IfLittleEndian (p); double dAlign = *reinterpret_cast(p); - p = pArgBase + 16; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase + 16; 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 = 1. / (2294.4871 * ::pow (2.0, -lEscale)); + lDataPos += 32; for (int id = 0; id < iNDets; id++) { - int iV = pData[lDataPos+1] + 256 * pData[lDataPos]; + int iV = pData[lDataPos+1] + (pData[lDataPos] << 8); 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; } Projections* -Projections::interpolateToParallel () +Projections::interpolateToParallel () const { if (m_geometry == Scanner::GEOMETRY_PARALLEL) - return this; + return const_cast(this); int nDet = m_nDet; int nView = m_nView; @@ -991,13 +973,17 @@ Projections::interpolateToParallel () 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); + ParallelRaysums parallel (this, ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI); double* pdThetaValuesForT = new double [pProjNew->nView()]; double* pdRaysumsForT = new double [pProjNew->nView()]; @@ -1008,10 +994,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; @@ -1032,8 +1019,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; @@ -1051,8 +1039,9 @@ Projections::interpolateToParallel () // /////////////////////////////////////////////////////////////////////////////// -ParallelRaysums::ParallelRaysums (Projections* pProjections) -: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()) +ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange) +: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()), + m_iThetaRange (iThetaRange), m_pCoordinates(NULL) { int iGeometry = pProjections->geometry(); double dDetInc = pProjections->detInc(); @@ -1060,9 +1049,10 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) 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] = new ParallelRaysumCoordinate; + m_vecpCoordinates[i] = m_pCoordinates + i; int iCoordinate = 0; for (int iV = 0; iV < m_iNumView; iV++) { @@ -1074,21 +1064,24 @@ 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); } - if (pC->m_dTheta >= PI) { // convert T/Theta to 0-PI interval - pC->m_dTheta -= PI; - pC->m_dT = -pC->m_dT - pProjections->detInc(); + 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; @@ -1098,8 +1091,7 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) ParallelRaysums::~ParallelRaysums() { - for (int i = 0; i < m_iNumCoordinates; i++) - delete m_vecpCoordinates[i]; + delete m_pCoordinates; } ParallelRaysums::CoordinateContainer& @@ -1180,11 +1172,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; @@ -1203,6 +1198,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])); }