X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=1aee958e6fdcd42c4addfee161471ccb38780a74;hp=c55aa37a5fcc94b6292ac9b912f93e858bbb5b7f;hb=d16eb37cbc73f67fc29a60645e0b1ac7fe32767e;hpb=999a754d1519a49ca062ee87b22bf601c1ee9f21 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index c55aa37..1aee958 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.67 2001/03/18 18:08: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,33 +688,35 @@ 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)) - return false; + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1.); std::complex** ppcDetValue = new std::complex* [m_nView]; unsigned int iView; for (iView = 0; iView < m_nView; iView++) { ppcDetValue[iView] = new std::complex [m_nDet]; + DetectorValue* detval = pProj->getDetectorArray (iView).detValues(); for (unsigned int iDet = 0; iDet < m_nDet; iDet++) - ppcDetValue[iView][iDet] = std::complex(getDetectorArray (iView).detValues()[iDet], 0); + ppcDetValue[iView][iDet] = std::complex(detval[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, 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; } @@ -722,6 +724,9 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) bool Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad) { +#ifndef HAVE_FFT + return false; +#else unsigned int nx = rIF.nx(); unsigned int ny = rIF.ny(); ImageFileArray v = rIF.getArray(); @@ -737,62 +742,76 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad return false; } -#ifndef HAVE_FFT - return false; -#else Array2d adView (nx, ny); Array2d adDet (nx, ny); double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); std::complex** ppcDetValue = new std::complex* [m_nView]; - unsigned int iView; - double* pdDet = new double [m_nDet]; - fftw_complex* pcIn = new fftw_complex [m_nDet]; - fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE); - for (iView = 0; iView < m_nView; iView++) { - unsigned int iDet; - for (iDet = 0; iDet < m_nDet; iDet++) { - pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet]; + int iNumDetWithZeros = ProcessSignal::addZeropadFactor (m_nDet, iZeropad); + double dZeropadRatio = iNumDetWithZeros / static_cast(m_nDet); + + double* pdDet = new double [iNumDetWithZeros]; + fftw_complex* pcIn = new fftw_complex [iNumDetWithZeros]; + fftw_plan plan = fftw_create_plan (iNumDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE); + + for (unsigned int iView = 0; iView < m_nView; iView++) { + DetectorValue* detval = getDetectorArray(iView).detValues(); + for (unsigned int iDet = 0; iDet < m_nDet; iDet++) { + pcIn[iDet].re = detval[iDet]; pcIn[iDet].im = 0; } + for (unsigned int iDet2 = m_nDet; iDet2 < iNumDetWithZeros; iDet2++) + pcIn[iDet2].re = pcIn[iDet2].im = 0; + fftw_one (plan, pcIn, NULL); - ppcDetValue[iView] = new std::complex [m_nDet]; - for (iDet = 0; iDet < m_nDet; iDet++) - ppcDetValue[iView][iDet] = std::complex (pcIn[iDet].re, pcIn[iDet].im); - Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet); + ppcDetValue[iView] = new std::complex [iNumDetWithZeros]; + for (unsigned int iD = 0; iD < iNumDetWithZeros; iD++) + ppcDetValue[iView][iD] = std::complex (pcIn[iD].re, pcIn[iD].im); +//ppcDetValue[iView][iD] = std::complex (std::abs(std::complex(pcIn[iD].re, pcIn[iD].im)), 0); + Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumDetWithZeros); } fftw_destroy_plan (plan); delete [] pcIn; - bool bError = calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet); + calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumDetWithZeros, dZeropadRatio); + interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumDetWithZeros, + iInterpolationID); - if (! bError) - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); - - for (iView = 0; iView < m_nView; iView++) - delete [] ppcDetValue[iView]; + for (int i = 0; i < m_nView; i++) + delete [] ppcDetValue[i]; delete [] ppcDetValue; - return bError; + return true; #endif } -bool -Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet) +void +Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet, + int iNumDetWithZeros, double dZeropadRatio) { double xMin = -phmLen() / 2; double xMax = xMin + phmLen(); double yMin = -phmLen() / 2; double yMax = yMin + phmLen(); - + double xCent = (xMin + xMax) / 2; + double yCent = (yMin + yMax) / 2; + + xMin = (xMin - xCent) * dZeropadRatio + xCent; + xMax = (xMax - xCent) * dZeropadRatio + xCent; + yMin = (yMin - yCent) * dZeropadRatio + yCent; + yMax = (yMax - yCent) * dZeropadRatio + yCent; + double xInc = (xMax - xMin) / nx; // size of cells double yInc = (yMax - yMin) / ny; - - int iDetCenter = (m_nDet - 1) / 2; // index refering to L=0 projection + + // +1 is correct for frequency data, ndet-1 is correct for projections + int iDetCenter = (iNumDetWithZeros - 1) / 2; // index refering to L=0 projection + if (isEven (iNumDetWithZeros)) + iDetCenter = (iNumDetWithZeros + 1) / 2; // Calculates polar coordinates (view#, det#) for each point on phantom grid double x = xMin + xInc / 2; // Rectang coords of center of pixel @@ -802,26 +821,26 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double double r = ::sqrt (x * x + y * y); double phi = atan2 (y, x); + if (phi < 0) + phi += TWOPI; + if (phi >= PI) { phi -= PI; - } else if (phi < 0) { - phi += PI; - } else r = -r; + } ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; ppdDet[ix][iy] = (r / m_detInc) + iDetCenter; } } - - return true; } void Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, - unsigned int nx, unsigned int ny, std::complex** ppcDetValue, - double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, int iInterpolationID) + unsigned int nx, unsigned int ny, std::complex** ppcDetValue, double** ppdView, + double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID) { + for (unsigned int ix = 0; ix < ny; ix++) { for (unsigned int iy = 0; iy < ny; iy++) { if (iInterpolationID == POLAR_INTERP_NEAREST) { @@ -831,7 +850,7 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, iView = 0; // iDet = m_nDet - iDet; } - if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) { + if (iDet >= 0 && iDet < nDetWithZeros && iView >= 0 && iView < nView) { v[ix][iy] = ppcDetValue[iView][iDet].real(); if (vImag) vImag[ix][iy] = ppcDetValue[iView][iDet].imag(); @@ -853,11 +872,11 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, v2 = ppcDetValue[iFloorView + 1][iFloorDet]; else v2 = ppcDetValue[0][iFloorDet]; - if (iFloorDet < nDet - 1) + if (iFloorDet < nDetWithZeros - 1) v4 = ppcDetValue[iFloorView][iFloorDet+1]; else v4 = v1; - if (iFloorView < nView - 1 && iFloorDet < nDet - 1) + if (iFloorView < nView - 1 && iFloorDet < nDetWithZeros - 1) v3 = ppcDetValue [iFloorView+1][iFloorDet+1]; else if (iFloorView < nView - 1) v3 = v2; @@ -885,94 +904,79 @@ 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 * m_detInc; + m_detStart = -(m_dFanBeamAngle / 2); m_rotInc = TWOPI / static_cast(iNViews); - m_rotStart = HALFPI; + m_rotStart = 0; 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; } +#if 1 + for (int k = iNDets - 2; k >= 0; k--) + detval[k+1] = detval[k]; + detval[0] = 0; +#endif } + 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 +995,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 + if (isEven (nDet)) // 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 +1016,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 +1041,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 +1061,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 +1071,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 +1086,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 +1113,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 +1194,14 @@ ParallelRaysums::getDetPositions (double* pdDetPos) } // locate by bisection, O(log2(n)) +// iLastFloor is used when sequential calls to interpolate with 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 +1220,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])); }