X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=8674647146a8c4f64f160c4323ef7be974b304d9;hb=a2bb680d159717dbd81e3ad179cb28adbf3a3e9a;hp=612ab6bf9161663dd42d6894e3bb1cfc371fc5de;hpb=ab78de6e33fff94c37847621458db3ad430a8bb3;p=ctsim.git diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 612ab6b..8674647 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.80 2002/06/27 03:19:23 kevin Exp $ +** $Id$ ** ** 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 @@ -25,7 +25,8 @@ ******************************************************************************/ #include "ct.h" -#include +#include +#include "interpolator.h" const kuint16 Projections::m_signature = ('P'*256 + 'J'); @@ -36,15 +37,15 @@ const int Projections::POLAR_INTERP_BICUBIC = 2; const char* const Projections::s_aszInterpName[] = { - {"nearest"}, - {"bilinear"}, + "nearest", + "bilinear", // {"bicubic"}, }; const char* const Projections::s_aszInterpTitle[] = { - {"Nearest"}, - {"Bilinear"}, + "Nearest", + "Bilinear", // {"Bicubic"}, }; @@ -1016,51 +1017,53 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) pProj = interpolateToParallel(); - int iInterpDet = nx; -// int iInterpDet = pProj->m_nDet; + int iInterpDet = static_cast(static_cast(sqrt(nx*nx+ny*ny))); int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad); - + double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.05); 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 = static_cast (fftw_malloc (sizeof(fftw_complex) * iNumInterpDetWithZeros)); + fftw_plan plan = fftw_plan_dft_1d (iNumInterpDetWithZeros, pcIn, pcIn, FFTW_FORWARD, FFTW_ESTIMATE); - fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros]; std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; - double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1) / SQRT2; + //double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1); + double dInterpScale = pProj->m_nDet / static_cast(iInterpDet); double dFFTScale = 1. / static_cast(iInterpDet * iInterpDet); int iMidPoint = iInterpDet / 2; double dMidPoint = static_cast(iInterpDet) / 2.; int iZerosAdded = iNumInterpDetWithZeros - iInterpDet; - // For each view, interpolate to nx length, shift to center at origin, and FFt transform + // For each view, interpolate, shift to center at origin, and FFT for (int iView = 0; iView < m_nView; iView++) { DetectorValue* detval = pProj->getDetectorArray(iView).detValues(); LinearInterpolator projInterp (detval, pProj->m_nDet); for (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; + pcIn[iDet][0] = projInterp.interpolate (dInterpPos) * dProjScale; + pcIn[iDet][1] = 0; } Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet); if (iZerosAdded > 0) { - for (int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++) - pcIn[iDet1+iZerosAdded] = pcIn[iDet1]; + for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--) { + pcIn[iDet1+iZerosAdded][0] = pcIn[iDet1][0]; + pcIn[iDet1+iZerosAdded][1] = pcIn[iDet1][1]; + } for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) - pcIn[iDet2].re = pcIn[iDet2].im = 0; + pcIn[iDet2][0] = pcIn[iDet2][1] = 0; } - fftw_one (plan, pcIn, NULL); + fftw_execute (plan); ppcDetValue[iView] = new std::complex [iNumInterpDetWithZeros]; for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) { - ppcDetValue[iView][iD] = std::complex (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); + ppcDetValue[iView][iD] = std::complex (pcIn[iD][0] * dFFTScale, pcIn[iD][1] * dFFTScale); } Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros); } - delete [] pcIn; + fftw_free(pcIn) ; fftw_destroy_plan (plan); @@ -1091,7 +1094,6 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double int iNumDetWithZeros, double dZeropadRatio, double dDetInc) { double dLength = viewDiameter(); -// double dLength = phmLen(); double xMin = -dLength / 2; double xMax = xMin + dLength; double yMin = -dLength / 2; @@ -1107,10 +1109,10 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double double xInc = (xMax - xMin) / nx; // size of cells double yInc = (yMax - yMin) / ny; + double dDetCenter = (iNumDetWithZeros - 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; + // if (isEven (iNumDetWithZeros)) + // dDetCenter = (iNumDetWithZeros + 0) / 2; // Calculates polar coordinates (view#, det#) for each point on phantom grid double x = xMin + xInc / 2; // Rectang coords of center of pixel @@ -1120,15 +1122,15 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double double r = ::sqrt (x * x + y * y); double phi = atan2 (y, x); - if (phi < 0) + if (phi <= -m_rotInc / 2) phi += TWOPI; - if (phi >= PI) { + if (phi >= PI - (m_rotInc / 2)) { phi -= PI; r = -r; } ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; - ppdDet[ix][iy] = (r / dDetInc) + iDetCenter; + ppdDet[ix][iy] = (r / dDetInc) + dDetCenter; } } } @@ -1140,30 +1142,29 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, { typedef std::complex complexValue; - BilinearInterpolator* pBilinear = NULL; + BilinearPolarInterpolator* pBilinear = NULL; + BicubicPolyInterpolator* pBicubic = NULL; if (iInterpolationID == POLAR_INTERP_BILINEAR) - pBilinear = new BilinearInterpolator (ppcDetValue, nView, nDetWithZeros); - - BicubicPolyInterpolator* pBicubic; - if (iInterpolationID == POLAR_INTERP_BICUBIC) + pBilinear = new BilinearPolarInterpolator (ppcDetValue, nView, nDetWithZeros); + else 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; - } + if (iView == nView) + iView = 0; 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 + } else { v[ix][iy] = 0; + if (vImag) + vImag[ix][iy] = 0; + } } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]);