+
+ 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 < pProj->m_nView; iView++)
+ delete [] ppcDetValue[iView];
+ delete [] ppcDetValue;
+
+ if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR)
+ delete pProj;
+
+ return true;
+}
+
+
+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();
+ if (! rIF.isComplex())
+ rIF.convertRealToComplex();
+ ImageFileArray vImag = rIF.getImaginaryArray();
+
+ if (! v || nx == 0 || ny == 0)
+ return false;
+
+ Projections* pProj = this;
+ if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR)
+ pProj = interpolateToParallel();
+
+ int iInterpDet = nx;
+// int iInterpDet = pProj->m_nDet;
+ int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad);
+
+ double dZeropadRatio = static_cast<double>(iNumInterpDetWithZeros) / static_cast<double>(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<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
+ double dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1) / SQRT2;
+
+ double dFFTScale = 1. / static_cast<double>(iInterpDet * iInterpDet);
+ int iMidPoint = iInterpDet / 2;
+ double dMidPoint = static_cast<double>(iInterpDet) / 2.;
+ int iZerosAdded = iNumInterpDetWithZeros - iInterpDet;
+
+ // For each view, interpolate to nx length, shift to center at origin, and FFt transform
+ for (unsigned int iView = 0; iView < m_nView; iView++) {
+ DetectorValue* detval = pProj->getDetectorArray(iView).detValues();
+ LinearInterpolator<DetectorValue> projInterp (detval, pProj->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<double> [iNumInterpDetWithZeros];
+ for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
+ ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale);
+ }
+
+ Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros);