X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=95cbc89a34c211c50e4abacf2f0a8bfafa0e56cb;hp=16b5b901e3bb77e8d7538caa4ab12969c1c14126;hb=663448e3173a19f054952806d8f8eca2fe59ec90;hpb=0f6257b32b46276415f1c6597fc1d992c3787071 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 16b5b90..95cbc89 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.65 2001/03/13 10:35:06 kevin Exp $ +** $Id: projections.cpp,v 1.66 2001/03/13 14:53:44 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 @@ -697,18 +697,18 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - if (! pProj->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(pProj->getDetectorArray (iView).detValues()[iDet], 0); + 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->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]; @@ -724,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(); @@ -739,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); - - if (! bError) - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumDetWithZeros, dZeropadRatio); + interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumDetWithZeros, + 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 (iNumDetWithZeros % 2 == 0) + 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,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) { @@ -833,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(); @@ -855,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;