r645: no message

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 2a80405a43dc1eb0356f1ddaf6d2d9bb6c58ed7e..968ee0ec48b1bdb807362ead34d0e187384da5bf 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.70 2001/03/24 05:28:28 kevin Exp $
+**  $Id: projections.cpp,v 1.71 2001/03/28 16:53:43 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 
 *
@@ -756,24 +757,32 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
   double dInterpScale = (m_nDet-1) / static_cast<double>(iInterpDet-1) / SQRT2;
   
   double dFFTScale = 1. / static_cast<double>(iInterpDet * iInterpDet);
-  dFFTScale /= 2; // Not sure why this scaling is necessary
+  int iMidPoint = iInterpDet / 2;
+  double dMidPoint = static_cast<double>(iInterpDet) / 2.;
+  int iZerosAdded = iNumInterpDetWithZeros - iInterpDet;
   for (unsigned int iView = 0; iView < m_nView; iView++) {
     DetectorValue* detval = getDetectorArray(iView).detValues();
     LinearInterpolator<DetectorValue> 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) * PI * SQRT2;
       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<double> [iNumInterpDetWithZeros];
-    for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++)
+    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);
   }