r4203: *** empty log message ***

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 612ab6bf9161663dd42d6894e3bb1cfc371fc5de..db8fd0c3d60a63632acbea37bacdd08387d3b8d5 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.80 2002/06/27 03:19:23 kevin Exp $
+**  $Id: projections.cpp,v 1.81 2003/03/15 10:27:30 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
@@ -25,7 +25,7 @@
 ******************************************************************************/
 
 #include "ct.h"
-#include <ctime>\r
+#include <ctime>
 
 const kuint16 Projections::m_signature = ('P'*256 + 'J');
 
@@ -1016,36 +1016,36 @@ 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<int>(static_cast<double>(sqrt(nx*nx+ny*ny)));
   int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad);
-
+  double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.5);
   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 dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1);
+  double dInterpScale = pProj->m_nDet / static_cast<double>(iInterpDet);
   
   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 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<DetectorValue> 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].re = projInterp.interpolate (dInterpPos) * dProjScale;
       pcIn[iDet].im = 0;
     }
 
     Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet);
     if (iZerosAdded > 0) {
-      for (int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++)
+       for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--)
         pcIn[iDet1+iZerosAdded] = pcIn[iDet1];
       for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) 
         pcIn[iDet2].re = pcIn[iDet2].im = 0;
@@ -1107,10 +1107,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 
@@ -1128,7 +1128,7 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double
       }
       
       ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc;
-      ppdDet[ix][iy] = (r / dDetInc) + iDetCenter;
+      ppdDet[ix][iy] = (r / dDetInc) + dDetCenter;
     }
   }
 }