r4228: Interpolation fixes

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index db8fd0c3d60a63632acbea37bacdd08387d3b8d5..64730029714c60dac9fb70dcc9630ff81448e60d 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.81 2003/03/15 10:27:30 kevin Exp $
+**  $Id: projections.cpp,v 1.82 2003/03/23 18:37:42 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
@@ -26,6 +26,7 @@
 
 #include "ct.h"
 #include <ctime>
+#include "interpolator.h"
 
 const kuint16 Projections::m_signature = ('P'*256 + 'J');
 
@@ -1091,7 +1092,6 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double
                                         int iNumDetWithZeros, double dZeropadRatio, double dDetInc)
 {
   double dLength = viewDiameter();
-//  double dLength = phmLen();
   double xMin = -dLength / 2;
   double xMax = xMin + dLength;
   double yMin = -dLength / 2;
@@ -1120,9 +1120,9 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double
       double r = ::sqrt (x * x + y * y);
       double phi = atan2 (y, x);
 
-      if (phi < 0)
+      if (phi <= -m_rotInc / 2)
         phi += TWOPI;
-      if (phi >= PI) {
+      if (phi >= PI - (m_rotInc / 2)) {
         phi -= PI;
         r = -r;
       }
@@ -1140,30 +1140,29 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
 {
   typedef std::complex<double> complexValue;
 
-  BilinearInterpolator<complexValue>* pBilinear =  NULL;  
+  BilinearPolarInterpolator<complexValue>* pBilinear = NULL;  
+  BicubicPolyInterpolator<complexValue>* pBicubic = NULL;  
   if (iInterpolationID == POLAR_INTERP_BILINEAR)
-    pBilinear = new BilinearInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
-
-  BicubicPolyInterpolator<complexValue>* pBicubic;  
-  if (iInterpolationID == POLAR_INTERP_BICUBIC)
+    pBilinear = new BilinearPolarInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
+  else if (iInterpolationID == POLAR_INTERP_BICUBIC)
     pBicubic = new BicubicPolyInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
-
+  
   for (unsigned int ix = 0; ix < ny; ix++) {
     for (unsigned int iy = 0; iy < ny; iy++) {
-
       if (iInterpolationID == POLAR_INTERP_NEAREST) {
         unsigned int iView = nearest<int> (ppdView[ix][iy]);
         unsigned int iDet = nearest<int> (ppdDet[ix][iy]);
-        if (iView == nView) {
-          iView = 0;
-          iDet = m_nDet - iDet;
-        }
+        if (iView == nView)
+         iView = 0;
         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();
-        } else
+        } else {
           v[ix][iy] = 0;
+          if (vImag)
+            vImag[ix][iy] = 0;
+       }
 
       } else if (iInterpolationID == POLAR_INTERP_BILINEAR) {
         std::complex<double> vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]);