X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=e5c79e7364f67a15a0a1ec5ce4b08f01ed0d2df7;hp=43e511f1b469bd09e3d2e42cfe7424e8e73750d0;hb=728e9fcbe0b785e56e7dcd2bd305786c647050af;hpb=209b7fc5d552957cc224265f51f1677ab7071e24 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 43e511f..e5c79e7 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.64 2001/03/13 09:25:18 kevin Exp $ +** $Id: projections.cpp,v 1.72 2001/03/29 21:25:49 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 @@ -50,6 +50,7 @@ const char* const Projections::s_aszInterpTitle[] = const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*); + /* NAME * Projections Constructor for projections matrix storage * @@ -697,20 +698,21 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - if (! pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet)) - return false; - - std::complex** ppcDetValue = new std::complex* [m_nView]; + std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; unsigned int iView; - 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(pProj->getDetectorArray (iView).detValues()[iDet], 0); + for (iView = 0; iView < pProj->m_nView; iView++) { + ppcDetValue[iView] = new std::complex [pProj->m_nDet]; + DetectorValue* detval = pProj->getDetectorArray (iView).detValues(); + for (unsigned int iDet = 0; iDet < pProj->m_nDet; iDet++) + ppcDetValue[iView][iDet] = std::complex(detval[iDet], 0); } - pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, iInterpolationID); + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, pProj->m_nDet, 1., pProj->m_detInc); - for (iView = 0; iView < m_nView; iView++) + pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, + pProj->m_nDet, iInterpolationID); + + for (iView = 0; iView < pProj->m_nView; iView++) delete [] ppcDetValue[iView]; delete [] ppcDetValue; @@ -724,6 +726,10 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) bool Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad) { +#ifndef HAVE_FFTW + rIF.arrayDataClear(); + return false; +#else unsigned int nx = rIF.nx(); unsigned int ny = rIF.ny(); ImageFileArray v = rIF.getArray(); @@ -738,63 +744,95 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad sys_error (ERR_WARNING, "convertFFTPolar supports Parallel only"); return false; } - -#ifndef HAVE_FFT - return false; -#else - Array2d adView (nx, ny); - Array2d adDet (nx, ny); - double** ppdView = adView.getArray(); - double** ppdDet = adDet.getArray(); + int iInterpDet = nx; + int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad); + + double dZeropadRatio = static_cast(iNumInterpDetWithZeros) / static_cast(iInterpDet); + + fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM); + + fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros]; 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]; + double dInterpScale = (m_nDet-1) / static_cast(iInterpDet-1) / SQRT2; + + double dFFTScale = 1. / static_cast(iInterpDet * iInterpDet); + int iMidPoint = iInterpDet / 2; + double dMidPoint = static_cast(iInterpDet) / 2.; + int iZerosAdded = iNumInterpDetWithZeros - iInterpDet; + for (unsigned int iView = 0; iView < m_nView; iView++) { + DetectorValue* detval = getDetectorArray(iView).detValues(); + LinearInterpolator projInterp (detval, m_nDet); + for (unsigned int iDet = 0; iDet < iInterpDet; iDet++) { + double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale; + pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale; pcIn[iDet].im = 0; } + Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet); + if (iZerosAdded > 0) { + for (unsigned int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++) + pcIn[iDet1+iZerosAdded] = pcIn[iDet1]; + for (unsigned int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; 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 [iNumInterpDetWithZeros]; + for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) { + ppcDetValue[iView][iD] = std::complex (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); + } + + Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros); } + delete [] pcIn; fftw_destroy_plan (plan); - delete [] pcIn; - bool bError = calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet); + Array2d adView (nx, ny); + Array2d adDet (nx, ny); + double** ppdView = adView.getArray(); + double** ppdDet = adDet.getArray(); + calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, + m_detInc * dInterpScale); - if (! bError) - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumInterpDetWithZeros, + 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 dDetInc) { - double xMin = -phmLen() / 2; - double xMax = xMin + phmLen(); - double yMin = -phmLen() / 2; - double yMax = yMin + phmLen(); - +// double dLength = viewDiameter(); + double dLength = phmLen(); + double xMin = -dLength / 2; + double xMax = xMin + dLength; + double yMin = -dLength / 2; + double yMax = yMin + dLength; + 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 @@ -804,84 +842,61 @@ 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; + ppdDet[ix][iy] = (r / dDetInc) + 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) { + typedef std::complex complexValue; + + BilinearInterpolator* pBilinear; + if (iInterpolationID == POLAR_INTERP_BILINEAR) + pBilinear = new BilinearInterpolator (ppcDetValue, nView, nDetWithZeros); + + BicubicPolyInterpolator* pBicubic; + if (iInterpolationID == POLAR_INTERP_BICUBIC) + pBicubic = new BicubicPolyInterpolator (ppcDetValue, nView, nDetWithZeros); + for (unsigned int ix = 0; ix < ny; ix++) { for (unsigned int iy = 0; iy < ny; iy++) { + if (iInterpolationID == POLAR_INTERP_NEAREST) { unsigned int iView = nearest (ppdView[ix][iy]); unsigned int iDet = nearest (ppdDet[ix][iy]); if (iView == nView) { iView = 0; - // iDet = m_nDet - iDet; + 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(); - } else { - sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", - ix, iy, ppdView[ix][iy], ppdDet[ix][iy]); + } else v[ix][iy] = 0; - } + } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { - unsigned int iFloorView = static_cast(ppdView[ix][iy]); - double dFracView = ppdView[ix][iy] - iFloorView; - unsigned int iFloorDet = static_cast(ppdDet[ix][iy]); - double dFracDet = ppdDet[ix][iy] - iFloorDet; - - if (iFloorDet >= 0 && iFloorView >= 0) { - std::complex v1 = ppcDetValue[iFloorView][iFloorDet]; - std::complex v2, v3, v4; - if (iFloorView < nView - 1) - v2 = ppcDetValue[iFloorView + 1][iFloorDet]; - else - v2 = ppcDetValue[0][iFloorDet]; - if (iFloorDet < nDet - 1) - v4 = ppcDetValue[iFloorView][iFloorDet+1]; - else - v4 = v1; - if (iFloorView < nView - 1 && iFloorDet < nDet - 1) - v3 = ppcDetValue [iFloorView+1][iFloorDet+1]; - else if (iFloorView < nView - 1) - v3 = v2; - else - v3 = ppcDetValue[0][iFloorDet+1]; - std::complex vInterp = (1 - dFracView) * (1 - dFracDet) * v1 + - dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 + - dFracDet * (1 - dFracView) * v4; - v[ix][iy] = vInterp.real(); - if (vImag) - vImag[ix][iy] = vInterp.imag(); - } else { - sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", - ix, iy, ppdView[ix][iy], ppdDet[ix][iy]); - v[ix][iy] = 0; - if (vImag) - vImag[ix][iy] = 0; - } + std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); + v[ix][iy] = vInterp.real(); + if (vImag) + vImag[ix][iy] = vInterp.imag(); } else if (iInterpolationID == POLAR_INTERP_BICUBIC) { - v[ix][iy] =0; - if (vImag) - vImag[ix][iy] = 0; + std::complex vInterp = pBicubic->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); + v[ix][iy] = vInterp.real(); + if (vImag) + vImag[ix][iy] = vInterp.imag(); } } } @@ -985,7 +1000,7 @@ Projections::interpolateToParallel () const #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::THETA_RANGE_NORMALIZE_TO_TWOPI); @@ -996,14 +1011,14 @@ Projections::interpolateToParallel () const // 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); + 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, &iLastFloor); + LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false); + detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); } } delete pdThetaValuesForT; @@ -1021,13 +1036,13 @@ Projections::interpolateToParallel () const detArray.setViewAngle (dViewAngle); for (int i = 0; i < pProjNew->nDet(); i++) - pdDetValueCopy[i] = detValues[i]; + 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, &iLastFloor); - } + LinearInterpolator interp (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), false); + for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) + detValues[iD] = interp.interpolate (dDetPos, &iLastFloor); } delete pdDetValueCopy; delete pdOriginalDetPositions; @@ -1175,36 +1190,3 @@ ParallelRaysums::getDetPositions (double* pdDetPos) iPos += m_iNumView; } } - -// 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, 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; - 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; - } - - if (iLastFloor) - *iLastFloor = iLower; - return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower])); -} -