X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=dc7248e62b9f3e42c86561612e99233c1bd6c245;hb=de411914da8b157958e9caae917bf1edeafbb713;hp=7c092f5990af43a9e267092a1f4a92042083ecc8;hpb=7cb74b55a21d0ab57a584ce3c910c280beb7a329;p=ctsim.git diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 7c092f5..dc7248e 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.52 2001/03/02 05:10:22 kevin Exp $ +** $Id: projections.cpp,v 1.59 2001/03/11 15:27:30 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(); } @@ -687,6 +687,11 @@ 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; + } Array2d adView (nx, ny); Array2d adDet (nx, ny); @@ -727,6 +732,11 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad if (! v || nx == 0 || ny == 0) return false; + if (m_geometry != Scanner::GEOMETRY_PARALLEL) { + sys_error (ERR_WARNING, "convertFFTPolar supports Parallel only"); + return false; + } + #ifndef HAVE_FFT return false; #else @@ -784,11 +794,6 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double int iDetCenter = (m_nDet - 1) / 2; // index refering to L=0 projection - if (m_geometry != Scanner::GEOMETRY_PARALLEL) { - sys_error (ERR_WARNING, "convertPolar supports Parallel only"); - return false; - } - // Calculates polar coordinates (view#, det#) for each point on phantom grid double x = xMin + xInc / 2; // Rectang coords of center of pixel for (unsigned int ix = 0; ix < nx; x += xInc, ix++) { @@ -820,8 +825,8 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, for (unsigned int ix = 0; ix < ny; ix++) { for (unsigned int iy = 0; iy < ny; iy++) { if (iInterpolationID == POLAR_INTERP_NEAREST) { - int iView = nearest (ppdView[ix][iy]); - int iDet = nearest (ppdDet[ix][iy]); + unsigned int iView = nearest (ppdView[ix][iy]); + unsigned int iDet = nearest (ppdDet[ix][iy]); if (iView == nView) { iView = 0; // iDet = m_nDet - iDet; @@ -836,9 +841,9 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, v[ix][iy] = 0; } } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { - int iFloorView = static_cast(ppdView[ix][iy]); + unsigned int iFloorView = static_cast(ppdView[ix][iy]); double dFracView = ppdView[ix][iy] - iFloorView; - int iFloorDet = static_cast(ppdDet[ix][iy]); + unsigned int iFloorDet = static_cast(ppdDet[ix][iy]); double dFracDet = ppdDet[ix][iy] - iFloorDet; if (iFloorDet >= 0 && iFloorView >= 0) { @@ -880,19 +885,18 @@ 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); - m_rotStart = 0; + m_rotStart = HALFPI; m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2; if (iNDets != 1024) @@ -903,41 +907,42 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa long lDataPos = 0; for (int iv = 0; iv < iNViews; iv++) { - long* pLong = reinterpret_cast(pData + lDataPos+0); + unsigned char* pArgBase = pData + lDataPos; + unsigned char* p = pArgBase+0; #ifndef WORDS_BIGENDIAN - SwapBytes4 (pLong); + SwapBytes4 (p); #endif - long lProjNumber = *pLong; + long lProjNumber = *reinterpret_cast(p); - pLong = reinterpret_cast(pData + lDataPos+20); + p = pArgBase+20; #ifndef WORDS_BIGENDIAN - SwapBytes4 (pLong); + SwapBytes4 (p); #endif - long lEscale = *pLong; + long lEscale = *reinterpret_cast(p); - pLong = reinterpret_cast(pData + lDataPos+28); + p = pArgBase+28; #ifndef WORDS_BIGENDIAN - SwapBytes4 (pLong); + SwapBytes4 (p); #endif - long lTime = *pLong; + long lTime = *reinterpret_cast(p); - float* pFloat = reinterpret_cast(pData + lDataPos+4); + p = pArgBase + 4; #ifndef WORDS_BIGENDIAN - SwapBytes4 (pFloat); + SwapBytes4 (p); #endif - double dAlpha = *pFloat; + double dAlpha = *reinterpret_cast(p) + HALFPI; - pFloat = reinterpret_cast(pData + lDataPos+12); + p = pArgBase+12; #ifndef WORDS_BIGENDIAN - SwapBytes4 (pFloat); + SwapBytes4 (p); #endif - double dAlign = *pFloat; + double dAlign = *reinterpret_cast(p); - pFloat = reinterpret_cast(pData + lDataPos+16); + p = pArgBase + 16; #ifndef WORDS_BIGENDIAN - SwapBytes4 (pFloat); + SwapBytes4 (p); #endif - double dMaxValue = *pFloat; + double dMaxValue = *reinterpret_cast(p); lDataPos += 32; double dEScale = pow (2.0, -lEscale); @@ -951,7 +956,7 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa double dTempScale = 2294.4871 * dEScale; for (int id = 0; id < iNDets; id++) { int iV = pData[lDataPos+1] + 256 * pData[lDataPos]; - if (iV > 32767) + if (iV > 32767) // two's complement signed conversion iV = iV - 65536; double dCosScale = cos ((id + 1 - iCenter) * dBetaInc); detval[id] = iV / (dTempScale * dCosScale); @@ -962,3 +967,247 @@ 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_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()) +{ + int iGeometry = pProjections->geometry(); + double dDetInc = pProjections->detInc(); + double dDetStart = pProjections->detStart(); + double dFocalLength = pProjections->focalLength(); + + m_iNumCoordinates = m_iNumView * m_iNumDet; + m_vecpCoordinates.reserve (m_iNumCoordinates); + for (int i = 0; i < m_iNumCoordinates; i++) + m_vecpCoordinates[i] = new ParallelRaysumCoordinate; + + 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 = normalizeAngle (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_dT = dFocalLength * sin(dFanAngle); + + } else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) { + // fan angle is same as dDetPos + pC->m_dTheta = normalizeAngle (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; + } +#endif + pC->m_dRaysum = detValues[iD]; + dDetPos += dDetInc; + } + } +} + +ParallelRaysums::~ParallelRaysums() +{ + for (int i = 0; i < m_iNumCoordinates; i++) + delete m_vecpCoordinates[i]; +} + +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; + } +} + +// 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])); +} +