X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=e5c79e7364f67a15a0a1ec5ce4b08f01ed0d2df7;hp=4d4a02d2ed560990419e565490fe8cfd0b3e142c;hb=728e9fcbe0b785e56e7dcd2bd305786c647050af;hpb=3ea498d51ce4597e9649cd21f155b51175ea0bea diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 4d4a02d..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.68 2001/03/21 21:45:31 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,21 +698,21 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - 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 (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 < m_nDet; iDet++) + for (unsigned int iDet = 0; iDet < pProj->m_nDet; iDet++) ppcDetValue[iView][iDet] = std::complex(detval[iDet], 0); } - pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1., m_detInc); + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, pProj->m_nDet, 1., pProj->m_detInc); 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++) + for (iView = 0; iView < pProj->m_nView; iView++) delete [] ppcDetValue[iView]; delete [] ppcDetValue; @@ -755,23 +756,32 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad std::complex** ppcDetValue = new std::complex* [m_nView]; 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 = iInterpDet * dInterpScale; - double dInterpPos = (m_nDet / 2.) + (iDet - iInterpDet/2.) * dInterpScale; - pcIn[iDet].re = projInterp.interpolate (dInterpPos); + double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale; + pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale; pcIn[iDet].im = 0; } - for (unsigned int iDet2 = iInterpDet; iDet2 < iNumInterpDetWithZeros; iDet2++) - pcIn[iDet2].re = pcIn[iDet2].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 [iNumInterpDetWithZeros]; - for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) - ppcDetValue[iView][iD] = std::complex (pcIn[iD].re / iInterpDet / (iInterpDet/2), pcIn[iD].im / iInterpDet / (iInterpDet/2)); + 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); } @@ -834,7 +844,6 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double if (phi < 0) phi += TWOPI; - if (phi >= PI) { phi -= PI; r = -r; @@ -852,7 +861,14 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID) { typedef std::complex complexValue; - BilinearInterpolator bilinear (ppcDetValue, nView, nDetWithZeros); + + 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++) { @@ -872,52 +888,15 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, v[ix][iy] = 0; } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { -#if 1 - std::complex vInterp = bilinear.interpolate (ppdView[ix][iy], ppdDet[ix][iy]); + std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); v[ix][iy] = vInterp.real(); if (vImag) vImag[ix][iy] = vInterp.imag(); -#else - int iFloorView = ::floor (ppdView[ix][iy]); - double dFracView = ppdView[ix][iy] - iFloorView; - int iFloorDet = ::floor (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 < nDetWithZeros - 1) - v4 = ppcDetValue[iFloorView][iFloorDet+1]; - else - v4 = v1; - if (iFloorView < nView - 1 && iFloorDet < nDetWithZeros - 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; - } -#endif } 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(); } } } @@ -1038,7 +1017,7 @@ Projections::interpolateToParallel () const int iLastFloor = -1; for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) { DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues(); - LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView()); + LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false); detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); } } @@ -1057,11 +1036,11 @@ 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; - LinearInterpolator interp (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet()); + 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); }