r648: no message

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index e5c79e7364f67a15a0a1ec5ce4b08f01ed0d2df7..8440bd49c9eac4c66b47846c9042764485b851cb 100644 (file)
--- 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.72 2001/03/29 21:25:49 kevin Exp $
+**  $Id: projections.cpp,v 1.73 2001/03/30 19:17:32 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
@@ -740,12 +740,12 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
   if (! v || nx == 0 || ny == 0)
     return false;
   
-  if (m_geometry != Scanner::GEOMETRY_PARALLEL) {
-    sys_error (ERR_WARNING, "convertFFTPolar supports Parallel only");
-    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);
@@ -753,21 +753,24 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
   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>* [m_nView];
-  double dInterpScale = (m_nDet-1) / static_cast<double>(iInterpDet-1) / SQRT2;
+  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 = getDetectorArray(iView).detValues();
-    LinearInterpolator<DetectorValue> projInterp (detval, m_nDet);
+    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++)
@@ -793,11 +796,14 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
   Array2d<double> adDet (nx, ny);
   double** ppdView = adView.getArray();
   double** ppdDet = adDet.getArray();
-  calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, 
-    m_detInc * dInterpScale);
+  pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, 
+    pProj->m_detInc * dInterpScale);
+
+  pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, 
+    iNumInterpDetWithZeros, iInterpolationID);
 
-  interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumInterpDetWithZeros,
-    iInterpolationID);
+  if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR)
+    delete pProj;
 
   for (int i = 0; i < m_nView; i++)
     delete [] ppcDetValue[i];