r355: Polar conversions of projections

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index bdf607369f96d645a7226af85d41a90128828fe2..22171504420064543e7d958bb168e461f9f9f79c 100644 (file)
--- a/libctsim/projections.cpp
+++ b/libctsim/projections.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.37 2001/01/04 21:28:41 kevin Exp $
+**  $Id: projections.cpp,v 1.38 2001/01/06 15:33:15 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
@@ -37,14 +37,14 @@ const char* Projections::s_aszInterpName[] =
 {
   {"nearest"},
   {"bilinear"},
-  {"bicubic"},
+//  {"bicubic"},
 };
 
 const char* Projections::s_aszInterpTitle[] = 
 {
   {"Nearest"},
   {"Bilinear"},
-  {"Bicubic"},
+//  {"Bicubic"},
 };
 
 const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*);
@@ -696,6 +696,58 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID)
   return true;
 }
 
+
+bool 
+Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
+{
+  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;
+  
+  Array2d<double> adView (nx, ny);
+  Array2d<double> adDet (nx, ny);
+  double** ppdView = adView.getArray();
+  double** ppdDet = adDet.getArray();
+
+  std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+  unsigned int iView;
+  double* pdDet = new double [m_nDet];
+  fftw_complex* pcIn = new fftw_complex [m_nDet];
+  fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE);
+
+  for (iView = 0; iView < m_nView; iView++) {
+    ppcDetValue[iView] = new std::complex<double> [m_nDet];
+    unsigned int iDet;
+    for (iDet = 0; iDet < m_nDet; iDet++) {
+      pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet];
+      pcIn[iDet].im = 0;
+    }
+    fftw_one (plan, pcIn, NULL);
+    for (iDet = 0; iDet < m_nDet; iDet++)
+      ppcDetValue[iView][iDet] = std::complex<double> (pcIn[iDet].re, pcIn[iDet].im);
+   Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet);
+  }
+
+  fftw_destroy_plan (plan);  
+  delete [] pcIn;
+  
+  calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
+
+  interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
+
+  for (iView = 0; iView < m_nView; iView++)
+    delete [] ppcDetValue[iView];
+  delete [] ppcDetValue;
+
+  return true;
+}
+
 void
 Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet)
 {
@@ -743,8 +795,10 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
       if (iInterpolationID == POLAR_INTERP_NEAREST) {
         int iView = nearest<int> (ppdView[ix][iy]);
         int iDet = nearest<int> (ppdDet[ix][iy]);
-        if (iView == nView)
-          iView = nView - 1;
+        if (iView == nView) {
+          iView = 0;
+       //   iDet = m_nDet - iDet;
+        }
         if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) {
           v[ix][iy] = ppcDetValue[iView][iDet].real();
           if (vImag)
@@ -766,7 +820,7 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
           if (iFloorView < nView - 1)
             v2 = ppcDetValue[iFloorView + 1][iFloorDet];
           else 
-            v2 = v1;
+            v2 = ppcDetValue[0][iFloorDet];
           if (iFloorDet < nDet - 1) 
             v4 = ppcDetValue[iFloorView][iFloorDet+1];
           else
@@ -774,10 +828,12 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
           if (iFloorView < nView - 1 && iFloorDet < nDet - 1)
             v3 = ppcDetValue [iFloorView+1][iFloorDet+1];
           else if (iFloorView < nView - 1)
-            v4 = v2;
+            v3 = v2;
           else 
-            v4 = v1;
-          std::complex<double> vInterp = v1 + dFracView * (v2 - v1) + dFracDet * (v4 - v1);
+            v3 = ppcDetValue[0][iFloorDet+1];
+          std::complex<double> vInterp = (1 - dFracView) * (1 - dFracDet) * v1 +
+            dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 +
+            dFracDet * (1 - dFracView) * v4;
           v[ix][iy] = vInterp.real();
           if (vImag)
             vImag[ix][iy] = vInterp.imag();
@@ -797,57 +853,4 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
   }
 }
 
-bool 
-Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
-{
-  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;
-  
-  Array2d<double> adView (nx, ny);
-  Array2d<double> adDet (nx, ny);
-  double** ppdView = adView.getArray();
-  double** ppdDet = adDet.getArray();
-
-  std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
-  unsigned int iView;
-  double* pdDet = new double [m_nDet];
-  fftw_complex* pcIn = new fftw_complex [m_nDet];
-  fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE);
-
-  for (iView = 0; iView < m_nView; iView++) {
-    ppcDetValue[iView] = new std::complex<double> [m_nDet];
-    unsigned int iDet;
-    for (iDet = 0; iDet < m_nDet; iDet++) {
-      pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet];
-      pcIn[iDet].im = 0;
-    }
-    fftw_one (plan, pcIn, NULL);
-    for (iDet = 0; iDet < m_nDet; iDet++)
-      ppcDetValue[iView][iDet] = std::complex<double> (pcIn[iDet].re, pcIn[iDet].im);
-    Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet);
-  }
-
-  fftw_destroy_plan (plan);
-  
-  delete [] pcIn;
-  
-  calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
-
-  interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
-
-  for (iView = 0; iView < m_nView; iView++)
-    delete [] ppcDetValue[iView];
-  delete [] ppcDetValue;
-
-  return true;
-}
-
-