r354: Added Projection Polar conversions
authorKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 4 Jan 2001 21:28:41 +0000 (21:28 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 4 Jan 2001 21:28:41 +0000 (21:28 +0000)
17 files changed:
ChangeLog
include/imagefile.h
include/projections.h
include/sgp.h
libctgraphics/ezplot.cpp
libctgraphics/sgp.cpp
libctsim/backprojectors.cpp
libctsim/imagefile.cpp
libctsim/projections.cpp
libctsim/scanner.cpp
msvc/ctsim/ctsim.plg
src/ctsim.h
src/dialogs.cpp
src/dialogs.h
src/dlgprojections.cpp
src/views.cpp
src/views.h

index a64a063..5e0220b 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,14 @@
 3.0.0beta1 - Released 
 
        * ctsim: Fixed initialization of min/max for PlotFiles
+
+       * ctsim: Added reset to full-intensity scale menu item
+
+       * ctsim: Add convert projections to polar plot
+
+       * ezplot: Cleaned up y-tick label placement
+
+       * sgp: Added better support for projection/reconstruction animation
        
 3.0.0alpha3 - Released 1/02/01
 
index eaa69af..aabb903 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: imagefile.h,v 1.28 2001/01/02 16:02:12 kevin Exp $
+**  $Id: imagefile.h,v 1.29 2001/01/04 21:28:41 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
@@ -151,6 +151,8 @@ public:
       : ImageFileBase ()
   {}
 
+  void getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter);
+
   bool convertRealToComplex ();
   bool convertComplexToReal ();
 
index 338f376..4c93304 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.h,v 1.19 2001/01/03 22:00:46 kevin Exp $
+**  $Id: projections.h,v 1.20 2001/01/04 21:28:41 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
 class Scanner;
 class DetectorArray;
 class Array2dFileLabel;
+class fnetorderstream;
 
+#include "array2dfile.h"
+#include "imagefile.h"
+#include <complex>
 
 // Projections
 class Projections
 {
  public:
+
+  static const int POLAR_INTERP_INVALID;
+  static const int POLAR_INTERP_NEAREST;
+  static const int POLAR_INTERP_BILINEAR;
+  static const int POLAR_INTERP_BICUBIC;
+
   Projections (const int nView, const int nDet);
   Projections (const Scanner& scanner);
   Projections ();
   ~Projections ();
 
+  static const int getInterpCount() {return s_iInterpCount;}
+  static const char** getInterpNameArray() {return s_aszInterpName;}
+  static const char** getInterpTitleArray() {return s_aszInterpTitle;}
+  static int convertInterpNameToID (const char* const interpName);
+  static const char* convertInterpIDToName (const int interpID);
+  static const char* convertInterpIDToTitle (const int interpID);
+
   void initFromScanner (const Scanner& scanner);
 
   void printProjectionData (int startView, int endView);
@@ -58,6 +75,9 @@ class Projections
 
   bool convertPolar (ImageFile& rIF, int iInterpolation);
   bool convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad);
+  void calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet);
+  void interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue,
+    double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, int iInterpolate);
 
   bool reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* frequencyFilterName, const char* const interpName, int interp_param, const char* const backprojName, const int trace) const;
 
@@ -124,6 +144,10 @@ class Projections
 
   const static kuint16 m_signature;
 
+  static const char* s_aszInterpName[];
+  static const char* s_aszInterpTitle[];
+  static const int s_iInterpCount;
+
   bool headerWrite (fnetorderstream& fs);
   bool headerRead (fnetorderstream& fs);
   void newProjData ();
index 2d83b2e..611654e 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: sgp.h,v 1.24 2001/01/02 16:02:13 kevin Exp $
+**  $Id: sgp.h,v 1.25 2001/01/04 21:28:41 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
@@ -125,6 +125,7 @@ private:
   double m_dCurrentWorldX;
   double m_dCurrentWorldY;
   double m_dTextAngle;
+  int m_iTextPointSize;
   bool m_bRecalcTransform;
   double m_dPointsPerPixel;  // points (72pt/in) per screen pixel;
   int m_iLinestyle;
@@ -145,6 +146,8 @@ private:
 #if HAVE_WXWINDOWS
   wxPen m_pen;
   wxFont* m_pFont;
+
+  void initFromDC (wxDC* pDC);
 #endif
 
 public:
index 260d5ba..ea29c76 100644 (file)
@@ -6,7 +6,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ezplot.cpp,v 1.27 2001/01/02 09:58:11 kevin Exp $
+**  $Id: ezplot.cpp,v 1.28 2001/01/04 21:28:41 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
@@ -723,18 +723,20 @@ EZPlot::plot (SGP* pSGP)
   if (o_yaxis == NOAXIS || o_ytlabel == FALSE)
     ytl_ofs = 0.0;
   else if (o_yticks == LEFT) {
-    double xExtent, yExtent;
+    double xExtentMin, xExtentMax, yExtent;
     char s[1024];
-    snprintf (s, sizeof(s), y_numfmt, 0);
-    m_pSGP->getTextExtent (s, &xExtent, &yExtent);
-    ytl_ofs = -2.0 * charwidth - xExtent;
+    snprintf (s, sizeof(s), y_numfmt, ymin);
+    m_pSGP->getTextExtent (s, &xExtentMin, &yExtent);
+    snprintf (s, sizeof(s), y_numfmt, ymax);
+    m_pSGP->getTextExtent (s, &xExtentMax, &yExtent);
+    if (xExtentMin > xExtentMax)
+      xExtentMax = xExtentMin;
+    ytl_ofs = -xExtentMax;
   } else if (o_yticks == RIGHT)
     ytl_ofs = 1.5 * charwidth;
   
   xa_max -= 0.7 * x_fldwid * charwidth; // make room for last x tick label
 
-
-
   xt_min = xa_min;
   yt_min = ya_min;
   xt_max = xa_max;
@@ -754,11 +756,11 @@ EZPlot::plot (SGP* pSGP)
   
   if (ytl_ofs != 0.0 && s_xcross == FALSE) {
     if (o_yticks == LEFT) {
-      xa_min += (2 + y_fldwid) * charwidth;
+      xa_min += 2*charwidth - ytl_ofs; // (2 + y_fldwid) * charwidth;
       xt_min = xa_min;
     } else if (o_yticks == RIGHT) {
       xa_min += 0.0;
-      xt_min = xa_min + ytl_ofs + y_fldwid * charwidth;
+      xt_min = xa_min + ytl_ofs; // + y_fldwid * charwidth;
     }
   } else
     xt_min = xa_min;
index ec129d2..263101f 100644 (file)
@@ -7,7 +7,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: sgp.cpp,v 1.25 2001/01/02 16:02:13 kevin Exp $
+**  $Id: sgp.cpp,v 1.26 2001/01/04 21:28:41 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
@@ -91,7 +91,26 @@ SGP::SGP (const SGPDriver& driver)
   ndc_to_mc.setIdentity();
   m_ctm.setIdentity();
 
-#if HAVE_WXWINDOWS
+#ifdef HAVE_WXWINDOWS
+  initFromDC (driver.idWX());
+#endif
+
+  setWindow (0., 0., 1., 1.);
+  setViewport (0., 0., 1., 1.);
+  moveAbs (0., 0.);
+  stylusNDC (0., 0., false);
+  
+  setTextAngle (0.);
+  setTextPointSize (12);
+  setColor (C_BLACK);
+  setLineStyle (LS_SOLID);
+}
+
+
+#ifdef HAVE_WXWINDOWS
+void
+SGP::initFromDC (wxDC* pDC)
+{
   m_pen.SetWidth (1);
 
   if (m_driver.isWX()) {
@@ -112,18 +131,9 @@ SGP::SGP (const SGPDriver& driver)
     m_dPointsPerPixel = iTestPointSize / dTestCharHeight;
        m_driver.idWX()->SetBackground (*wxWHITE_BRUSH);
   }
+}
 #endif
 
-  setWindow (0., 0., 1., 1.);
-  setViewport (0., 0., 1., 1.);
-  moveAbs (0., 0.);
-  stylusNDC (0., 0., false);
-  
-  setTextAngle (0.);
-  setTextPointSize (12);
-  setColor (C_BLACK);
-  setLineStyle (LS_SOLID);
-}
 
 SGP::~SGP()
 {
@@ -519,7 +529,8 @@ SGP::setTextPointSize (double height)
 #endif
 #if HAVE_WXWINDOWS
   if (m_driver.isWX()) {
-      m_pFont->SetPointSize (static_cast<int>(height+0.5));
+    m_iTextPointSize = static_cast<int>(height+0.5);
+      m_pFont->SetPointSize (m_iTextPointSize);
       m_driver.idWX()->SetFont (*m_pFont);
   }
 #endif
@@ -941,7 +952,10 @@ const unsigned char SGP::MARKER_BITMAP[MARK_COUNT][5] =
 void
 SGP::setDC (wxDC* pDC)
 {
-  if (m_driver.isWX())
+  if (m_driver.isWX()) {
     m_driver.setDC(pDC);
+    initFromDC (pDC);
+    setTextPointSize (m_iTextPointSize);
+  }
 }
 #endif
index 1f8a1d2..8343314 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.22 2001/01/02 16:02:13 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.23 2001/01/04 21:28:41 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
@@ -101,7 +101,7 @@ Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char
 {
   m_fail = false;
   m_pBackprojectImplem = NULL;
-
+  
   initBackprojector (proj, im, backprojName, interpName, interpFactor);
 }
 
@@ -142,35 +142,35 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
     m_failMessage = "Invalid interpolation name ";
     m_failMessage += interpName;
   }
-
+  
   if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
     m_fail = true;
     return false;
   }
-
+  
   if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
-      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) 
-      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
   else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
-      if (m_idBackproject == BPROJ_TRIG)
-         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
-      else if (m_idBackproject == BPROJ_TABLE)
-         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
-      else if (m_idBackproject == BPROJ_DIFF)
-         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
-      else if (m_idBackproject == BPROJ_DIFF2)
-         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
-      else if (m_idBackproject == BPROJ_IDIFF2)
-         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
-      else if (m_idBackproject == BPROJ_IDIFF3)
-         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
+    if (m_idBackproject == BPROJ_TRIG)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
+    else if (m_idBackproject == BPROJ_TABLE)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
+    else if (m_idBackproject == BPROJ_DIFF)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
+    else if (m_idBackproject == BPROJ_DIFF2)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
+    else if (m_idBackproject == BPROJ_IDIFF2)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
+    else if (m_idBackproject == BPROJ_IDIFF3)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
   } else {
-      m_fail = true;
-      m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
-      return false;
+    m_fail = true;
+    m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
+    return false;
   }
-
+  
   return true;
 }
 
@@ -179,24 +179,24 @@ int
 Backprojector::convertBackprojectNameToID (const char* const backprojName)
 {
   int backprojID = BPROJ_INVALID;
-
+  
   for (int i = 0; i < s_iBackprojectCount; i++)
-      if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
-         backprojID = i;
-         break;
-      }
-
-  return (backprojID);
+    if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
+      backprojID = i;
+      break;
+    }
+    
+    return (backprojID);
 }
 
 const char*
 Backprojector::convertBackprojectIDToName (int bprojID)
 {
   static const char *bprojName = "";
-
+  
   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
-      return (s_aszBackprojectName[bprojID]);
-
+    return (s_aszBackprojectName[bprojID]);
+  
   return (bprojName);
 }
 
@@ -204,10 +204,10 @@ const char*
 Backprojector::convertBackprojectIDToTitle (const int bprojID)
 {
   static const char *bprojTitle = "";
-
+  
   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
-      return (s_aszBackprojectTitle[bprojID]);
-
+    return (s_aszBackprojectTitle[bprojID]);
+  
   return (bprojTitle);
 }
 
@@ -216,24 +216,24 @@ int
 Backprojector::convertInterpNameToID (const char* const interpName)
 {
   int interpID = INTERP_INVALID;
-
+  
   for (int i = 0; i < s_iInterpCount; i++)
-      if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
-         interpID = i;
-         break;
-      }
-
-  return (interpID);
+    if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
+      interpID = i;
+      break;
+    }
+    
+    return (interpID);
 }
 
 const char*
 Backprojector::convertInterpIDToName (const int interpID)
 {
   static const char *interpName = "";
-
+  
   if (interpID >= 0 && interpID < s_iInterpCount)
-      return (s_aszInterpName[interpID]);
-
+    return (s_aszInterpName[interpID]);
+  
   return (interpName);
 }
 
@@ -241,10 +241,10 @@ const char*
 Backprojector::convertInterpIDToTitle (const int interpID)
 {
   static const char *interpTitle = "";
-
+  
   if (interpID >= 0 && interpID < s_iInterpCount)
-      return (s_aszInterpTitle[interpID]);
-
+    return (s_aszInterpTitle[interpID]);
+  
   return (interpTitle);
 }
 
@@ -257,33 +257,33 @@ Backprojector::convertInterpIDToTitle (const int interpID)
 //   Pure virtual base class for all backprojectors.
 
 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
-    : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
+: proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
 {
   detInc = proj.detInc();
   nDet = proj.nDet();
   iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection 
   rotScale = proj.rotInc();
-
+  
   if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
-       rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
+    rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
-       rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
+    rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
   else
-         sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
-
+    sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
+  
   v = im.getArray();
   nx = im.nx();
   ny = im.ny();
   im.arrayDataClear();
-
+  
   xMin = -proj.phmLen() / 2;      // Retangular coords of phantom
   xMax = xMin + proj.phmLen();
   yMin = -proj.phmLen() / 2;
   yMax = yMin + proj.phmLen();
-
+  
   xInc = (xMax - xMin) / nx;   // size of cells
   yInc = (yMax - yMin) / ny;
-
+  
   m_dFocalLength = proj.focalLength();
 }
 
@@ -300,8 +300,8 @@ Backproject::ScaleImageByRotIncrement ()
 
 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
 {
-    sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
-    errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
+  sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
+  errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
 }
 
 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
@@ -313,7 +313,7 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl
   os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
   os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
   os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
-
+  
   sys_error (ERR_WARNING, os.str().c_str());
 #endif
 }
@@ -329,7 +329,7 @@ void
 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double theta = view_angle;
-
+  
   double x = xMin + xInc / 2;  // Rectang coords of center of pixel 
   for (int ix = 0; ix < nx; x += xInc, ix++) {
     double y = yMin + yInc / 2;
@@ -337,23 +337,23 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
       double r = sqrt (x * x + y * y);   // distance of cell from center
       double phi = atan2 (y, x);         // angle of cell from center
       double L = r * cos (theta - phi);  // position on detector
-
+      
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
-
-       if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
-           errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
-       else
-         v[ix][iy] += rotScale * filteredProj[iDetPos];
+        int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
+        
+        if (iDetPos < 0 || iDetPos >= nDet)    // check for impossible: index outside of raysum pos 
+          errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
+        else
+          v[ix][iy] += rotScale * filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-         double p = L / detInc;        // position along detector
-         double pFloor = floor (p);
-         int iDetPos = iDetCenter + static_cast<int>(pFloor);
-         double frac = p - pFloor;     // fraction distance from det
-         if (iDetPos < 0 || iDetPos >= nDet - 1)       // check for impossible: index outside of raysum pos 
-           errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
-         else
-           v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+        double p = L / detInc; // position along detector
+        double pFloor = floor (p);
+        int iDetPos = iDetCenter + static_cast<int>(pFloor);
+        double frac = p - pFloor;      // fraction distance from det
+        if (iDetPos < 0 || iDetPos >= nDet - 1)        // check for impossible: index outside of raysum pos 
+          errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
+        else
+          v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }
   }
@@ -367,13 +367,13 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 //   Precalculates trigometric function value for each point in image for backprojection.
 
 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
-  : Backproject (proj, im, interpType, interpFactor)
+: Backproject (proj, im, interpType, interpFactor)
 {
   arrayR.initSetSize (im.nx(), im.ny());
   arrayPhi.initSetSize (im.nx(), im.ny());
   r = arrayR.getArray();
   phi = arrayPhi.getArray();
-
+  
   double x, y;                 // Rectang coords of center of pixel 
   int ix, iy;
   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
@@ -392,29 +392,29 @@ void
 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double theta = view_angle;
-
+  
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
-
+    
     for (int iy = 0; iy < ny; iy++) {
       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
-
+      
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector 
-
-       if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
-         errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
-       else
-         pImCol[iy] += filteredProj[iDetPos];
+        int iDetPos = iDetCenter + nearest<int>(L / detInc);   // calc index in the filtered raysum vector 
+        
+        if (iDetPos < 0 || iDetPos >= nDet)    // check for impossible: index outside of raysum pos 
+          errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
+        else
+          pImCol[iy] += filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-       double dPos = L / detInc;               // position along detector 
-       double dPosFloor = floor (dPos);
-       int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
-       double frac = dPos - dPosFloor; // fraction distance from det 
-       if (iDetPos < 0 || iDetPos >= nDet - 1)
-           errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
-       else
-         pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+        double dPos = L / detInc;              // position along detector 
+        double dPosFloor = floor (dPos);
+        int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+        double frac = dPos - dPosFloor;        // fraction distance from det 
+        if (iDetPos < 0 || iDetPos >= nDet - 1)
+          errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
+        else
+          pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }  // end for y 
   }    // end for x 
@@ -429,14 +429,14 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
 //   Iterates in x & y direction by adding difference in L position
 
 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
-  :  Backproject (proj, im, interpType, interpFactor)
+:  Backproject (proj, im, interpType, interpFactor)
 {
   // calculate center of first pixel v[0][0] 
   double x = xMin + xInc / 2;
   double y = yMin + yInc / 2;
   start_r = sqrt (x * x + y * y);
   start_phi = atan2 (y, x);
-
+  
   im.arrayDataClear();
 }
 
@@ -452,31 +452,31 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
   double det_dx = xInc * cos (theta);
   double det_dy = yInc * sin (theta);
   double lColStart = start_r * cos (theta - start_phi);  // calculate L for first point in image
-       
+  
   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
     double curDetPos = lColStart;
     ImageFileColumn pImCol = v[ix];
-  
+    
     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 #ifdef DEBUG
       printf ("[%2d,%2d]:  %8.5f  ", ix, iy, curDetPos);
 #endif
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc);    // calc index in the filtered raysum vector 
-
-       if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
-           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-       else
-         pImCol[iy] += filteredProj[iDetPos];
+        int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc);   // calc index in the filtered raysum vector 
+        
+        if (iDetPos < 0 || iDetPos >= nDet)    // check for impossible: index outside of raysum pos 
+          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+        else
+          pImCol[iy] += filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-       double detPos = curDetPos / detInc;             // position along detector 
-       double detPosFloor = floor (detPos);
-       int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
-       double frac = detPos - detPosFloor;     // fraction distance from det 
-       if (iDetPos < 0 || iDetPos >= nDet - 1)
-           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-       else
-         pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+        double detPos = curDetPos / detInc;            // position along detector 
+        double detPosFloor = floor (detPos);
+        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
+        double frac = detPos - detPosFloor;    // fraction distance from det 
+        if (iDetPos < 0 || iDetPos >= nDet - 1)
+          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+        else
+          pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }  // end for y 
   }    // end for x 
@@ -493,40 +493,40 @@ void
 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double theta = view_angle;
-
+  
   // Distance betw. detectors for an angle given in units of detectors 
   double det_dx = xInc * cos (theta) / detInc;
   double det_dy = yInc * sin (theta) / detInc;
-
+  
   // calculate detPosition for first point in image (ix=0, iy=0) 
   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
-       
+  
 #ifdef DEBUG
   printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale);
 #endif
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
     double curDetPos = detPosColStart;
     ImageFileColumn pImCol = v[ix];
-
+    
     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 #ifdef DEBUG
       printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
 #endif
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector 
-       
-       if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
-           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-       else
-         *pImCol++ += filteredProj[iDetPos];
+        int iDetPos = iDetCenter + nearest<int> (curDetPos);   // calc index in the filtered raysum vector 
+        
+        if (iDetPos < 0 || iDetPos >= nDet)    // check for impossible: index outside of raysum pos 
+          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+        else
+          *pImCol++ += filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-       double detPosFloor = floor (curDetPos);
-       int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
-       double frac = curDetPos - detPosFloor;  // fraction distance from det 
-       if (iDetPos < 0 || iDetPos >= nDet - 1)
-           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-       else
-         *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
+        double detPosFloor = floor (curDetPos);
+        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
+        double frac = curDetPos - detPosFloor; // fraction distance from det 
+        if (iDetPos < 0 || iDetPos >= nDet - 1)
+          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+        else
+          *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
       }
     }  // end for y
   }    // end for x
@@ -542,43 +542,43 @@ void
 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double theta = view_angle;
-
+  
   static const kint32 scale = 1 << 16;
   static const double dScale = scale;
   static const kint32 halfScale = scale / 2;
-
+  
   const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
   const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
-
+  
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
-       
+  
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
     kint32 curDetPos = detPosColStart;
     ImageFileColumn pImCol = v[ix];
-
+    
     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
-       int iDetPos = iDetCenter + detPosNearest;       // calc index in the filtered raysum vector 
-
-       if (iDetPos < 0 || iDetPos >= nDet)  // check for index outside of raysum pos 
-           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-       else
-         *pImCol++ += filteredProj[iDetPos];
+        int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
+        int iDetPos = iDetCenter + detPosNearest;      // calc index in the filtered raysum vector 
+        
+        if (iDetPos < 0 || iDetPos >= nDet)  // check for index outside of raysum pos 
+          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+        else
+          *pImCol++ += filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-       kint32 detPosFloor = curDetPos / scale;
-       kint32 detPosRemainder = curDetPos % scale;
-       if (detPosRemainder < 0) {
-         detPosFloor--;
-         detPosRemainder += scale;
-       }
-       int iDetPos = iDetCenter + detPosFloor;
-       double frac = detPosRemainder / dScale;
-       if (iDetPos < 0 || iDetPos >= nDet - 1)
-           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-       else
-         *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+        kint32 detPosFloor = curDetPos / scale;
+        kint32 detPosRemainder = curDetPos % scale;
+        if (detPosRemainder < 0) {
+          detPosFloor--;
+          detPosRemainder += scale;
+        }
+        int iDetPos = iDetCenter + detPosFloor;
+        double frac = detPosRemainder / dScale;
+        if (iDetPos < 0 || iDetPos >= nDet - 1)
+          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+        else
+          *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }  // end for y
   }    // end for x
@@ -599,13 +599,13 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
   static const kint32 scaleBitmask = scale - 1;
   static const kint32 halfScale = scale / 2;
   static const double dInvScale = 1. / scale;
-
+  
   const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
   const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
-
+  
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
-       
+  
   // precalculate scaled difference for linear interpolation
   double* deltaFilteredProj = new double [nDet];
   if (interpType == Backprojector::INTERP_LINEAR) {
@@ -613,34 +613,34 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
   }
   deltaFilteredProj[nDet - 1] = 0;  // last detector
-
+  
   int iLastDet = nDet - 1;
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
     kint32 curDetPos = detPosColStart;
     ImageFileColumn pImCol = v[ix];
-
+    
     if (interpType == Backprojector::INTERP_NEAREST) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-       const int iDetPos = (curDetPos + halfScale) >> 16;
-       if (iDetPos >= 0 && iDetPos <= iLastDet)
-         *pImCol++ += filteredProj[iDetPos];
+        const int iDetPos = (curDetPos + halfScale) >> 16;
+        if (iDetPos >= 0 && iDetPos <= iLastDet)
+          *pImCol++ += filteredProj[iDetPos];
       }        // end for iy
     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-       const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
-       if (iDetPos >= 0 && iDetPos <= iLastDet)
-       *pImCol++ += filteredProj[iDetPos];
+        const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
+        if (iDetPos >= 0 && iDetPos <= iLastDet)
+          *pImCol++ += filteredProj[iDetPos];
       }        // end for iy
     } else if (interpType == Backprojector::INTERP_LINEAR) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-       const kint32 iDetPos = curDetPos >> scaleShift;
-       const kint32 detRemainder = curDetPos & scaleBitmask;
-       if (iDetPos >= 0 && iDetPos <= iLastDet)
-         *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+        const kint32 iDetPos = curDetPos >> scaleShift;
+        const kint32 detRemainder = curDetPos & scaleBitmask;
+        if (iDetPos >= 0 && iDetPos <= iLastDet)
+          *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
       }        // end for iy
     } //end linear
   } // end for ix
-
+  
   delete deltaFilteredProj;
 }
 
@@ -649,10 +649,10 @@ void
 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double beta = view_angle;
-
+  
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
-
+    
     for (int iy = 0; iy < ny; iy++) { 
       double dAngleDiff = beta - phi[ix][iy];
       double rcos_t = r[ix][iy] * cos (dAngleDiff);
@@ -660,23 +660,23 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const
       double dFLPlusSin = m_dFocalLength + rsin_t;
       double gamma =  atan (rcos_t / dFLPlusSin);
       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
-
+      
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int iDetPos =iDetCenter + nearest<int>(gamma / detInc); // calc index in the filtered raysum vector 
-
-       if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos  
-         ; //    errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
-       } else
-         pImCol[iy] += filteredProj[iDetPos] / dL2;
+        int iDetPos =iDetCenter + nearest<int>(gamma / detInc);        // calc index in the filtered raysum vector 
+        
+        if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos  
+          ; //   errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+        } else
+          pImCol[iy] += filteredProj[iDetPos] / dL2;
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-       double dPos = gamma / detInc;           // position along detector 
-       double dPosFloor = floor (dPos);
-       int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
-       double frac = dPos - dPosFloor; // fraction distance from det 
-       if (iDetPos < 0 || iDetPos >= nDet - 1) {
-         ; //      errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
-       } else
-         pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
+        double dPos = gamma / detInc;          // position along detector 
+        double dPosFloor = floor (dPos);
+        int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+        double frac = dPos - dPosFloor;        // fraction distance from det 
+        if (iDetPos < 0 || iDetPos >= nDet - 1) {
+          ; //     errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+        } else
+          pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
       }
     }  // end for y 
   }    // end for x 
@@ -686,37 +686,37 @@ void
 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double beta = view_angle;
-
+  
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
-
+    
     for (int iy = 0; iy < ny; iy++) {
       double dAngleDiff = beta - phi[ix][iy];
       double rcos_t = r[ix][iy] * cos (dAngleDiff);
       double rsin_t = r[ix][iy] * sin (dAngleDiff);
-
+      
       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
       double dDetPos =  rcos_t / dU;
       // double to scale for imaginary detector that passes through origin
       // of phantom, see Kak-Slaney Figure 3.22
       dDetPos *= 2; 
-
+      
       if (interpType == Backprojector::INTERP_NEAREST) {
-       int iDetPos = iDetCenter + nearest<int>(dDetPos / detInc);      // calc index in the filtered raysum vector 
-
-       if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
-         ; ///   errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
-       else
-         pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
+        int iDetPos = iDetCenter + nearest<int>(dDetPos / detInc);     // calc index in the filtered raysum vector 
+        
+        if (iDetPos < 0 || iDetPos >= nDet)    // check for impossible: index outside of raysum pos 
+          ; ///          errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+        else
+          pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-       double dPos = dDetPos / detInc;         // position along detector 
-       double dPosFloor = floor (dPos);
-       int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
-       double frac = dPos - dPosFloor; // fraction distance from det 
-       if (iDetPos < 0 || iDetPos >= nDet - 1)
-         ; //      errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
-       else
-         pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU);
+        double dPos = dDetPos / detInc;                // position along detector 
+        double dPosFloor = floor (dPos);
+        int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+        double frac = dPos - dPosFloor;        // fraction distance from det 
+        if (iDetPos < 0 || iDetPos >= nDet - 1)
+          ; //     errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+        else
+          pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU);
       }
     }  // end for y 
   }    // end for x 
index df0e51e..946e4db 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: imagefile.cpp,v 1.33 2001/01/02 16:02:13 kevin Exp $
+**  $Id: imagefile.cpp,v 1.34 2001/01/04 21:28:41 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
@@ -84,21 +84,29 @@ F64Image::F64Image (void)
 }
 
 void 
-ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
+ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter)
 {
-  ImageFileArray v = getArray();
-  SignalFilter filter (filterName, domainName, bw, filt_param);
-  
-  int iXCenter, iYCenter;
   if (isEven (m_nx))
     iXCenter = m_nx / 2;
   else
     iXCenter = (m_nx - 1) / 2;
+
   if (isEven (m_ny))
     iYCenter = m_ny / 2;
   else
     iYCenter = (m_ny - 1) / 2;
+}
+
+
+void 
+ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
+{
+  ImageFileArray v = getArray();
+  SignalFilter filter (filterName, domainName, bw, filt_param);
   
+  unsigned int iXCenter, iYCenter;
+  getCenterCoordinates (iXCenter, iYCenter);
+
   for (unsigned int ix = 0; ix < m_nx; ix++)
     for (unsigned int iy = 0; iy < m_ny; iy++) {
       long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter));
index bdcf35e..bdf6073 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.36 2001/01/03 22:00:46 kevin Exp $
+**  $Id: projections.cpp,v 1.37 2001/01/04 21:28:41 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
 
 const kuint16 Projections::m_signature = ('P'*256 + 'J');
 
+const int Projections::POLAR_INTERP_INVALID = -1;
+const int Projections::POLAR_INTERP_NEAREST = 0;
+const int Projections::POLAR_INTERP_BILINEAR = 1;
+const int Projections::POLAR_INTERP_BICUBIC = 2;
+
+const char* Projections::s_aszInterpName[] = 
+{
+  {"nearest"},
+  {"bilinear"},
+  {"bicubic"},
+};
+
+const char* Projections::s_aszInterpTitle[] = 
+{
+  {"Nearest"},
+  {"Bilinear"},
+  {"Bicubic"},
+};
+
+const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*);
+
+
 /* NAME
- *    Projections              Constructor for projections matrix storage 
- *
- * SYNOPSIS
- *    proj = projections_create (filename, nView, nDet)
- *    Projections& proj                Allocated projections structure & matrix
- *    int nView                        Number of rotated view
- *    int nDet                 Number of detectors
- *
- */
+*    Projections               Constructor for projections matrix storage 
+*
+* SYNOPSIS
+*    proj = projections_create (filename, nView, nDet)
+*    Projections& proj         Allocated projections structure & matrix
+*    int nView                 Number of rotated view
+*    int nDet                  Number of detectors
+*
+*/
 
 Projections::Projections (const Scanner& scanner)
-  : m_projData(0)
+: m_projData(0)
 {
   initFromScanner (scanner);
 }
 
 
 Projections::Projections (const int nView, const int nDet)
-  : m_projData(0)
+: m_projData(0)
 {
   init (nView, nDet);
 }
 
 Projections::Projections (void)
-  : m_projData(0)
+: m_projData(0)
 {
   init (0, 0);
 }
@@ -63,6 +85,43 @@ Projections::~Projections (void)
   deleteProjData();
 }
 
+int
+Projections::convertInterpNameToID (const char* const interpName)
+{
+  int interpID = POLAR_INTERP_INVALID;
+  
+  for (int i = 0; i < s_iInterpCount; i++)
+    if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
+      interpID = i;
+      break;
+    }
+    
+    return (interpID);
+}
+
+const char*
+Projections::convertInterpIDToName (const int interpID)
+{
+  static const char *interpName = "";
+  
+  if (interpID >= 0 && interpID < s_iInterpCount)
+    return (s_aszInterpName[interpID]);
+  
+  return (interpName);
+}
+
+const char*
+Projections::convertInterpIDToTitle (const int interpID)
+{
+  static const char *interpTitle = "";
+  
+  if (interpID >= 0 && interpID < s_iInterpCount)
+    return (s_aszInterpTitle[interpID]);
+  
+  return (interpTitle);
+}
+
+
 
 void
 Projections::init (const int nView, const int nDet)
@@ -71,7 +130,7 @@ Projections::init (const int nView, const int nDet)
   m_nView = nView;
   m_nDet = nDet;
   newProjData ();
-
+  
   time_t t = time (NULL);
   tm* lt = localtime (&t);
   m_year = lt->tm_year;
@@ -88,7 +147,7 @@ Projections::initFromScanner (const Scanner& scanner)
   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
   deleteProjData();
   init (scanner.nView(), scanner.nDet());
-
+  
   m_phmLen = scanner.phmLen();
   m_rotInc = scanner.rotInc();
   m_detInc = scanner.detInc();
@@ -97,12 +156,6 @@ Projections::initFromScanner (const Scanner& scanner)
   m_fieldOfView = scanner.fieldOfView();
   m_rotStart = 0;
   m_detStart =  -(scanner.detLen() / 2);
-#if 0
-  if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
-    m_detInc /= 2;
-       std::cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
-  }
-#endif
 }
 
 void
@@ -120,10 +173,10 @@ Projections::newProjData (void)
 {
   if (m_projData)
     sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
-
+  
   if (m_nView > 0 && m_nDet) {
     m_projData = new DetectorArray* [m_nView];
-
+    
     for (int i = 0; i < m_nView; i++)
       m_projData[i] = new DetectorArray (m_nDet);
   }
@@ -131,12 +184,12 @@ Projections::newProjData (void)
 
 
 /* NAME
- *   projections_free                  Free memory allocated to projections
- *
- * SYNOPSIS
- *   projections_free(proj)
- *   Projections& proj                         Projectionss to be deallocated
- */
+*   projections_free                   Free memory allocated to projections
+*
+* SYNOPSIS
+*   projections_free(proj)
+*   Projections& proj                          Projectionss to be deallocated
+*/
 
 void 
 Projections::deleteProjData (void)
@@ -144,7 +197,7 @@ Projections::deleteProjData (void)
   if (m_projData != NULL) {
     for (int i = 0; i < m_nView; i++)
       delete m_projData[i];
-
+    
     delete m_projData;
     m_projData = NULL;
   }
@@ -152,9 +205,9 @@ Projections::deleteProjData (void)
 
 
 /* NAME
- *    Projections::headerWwrite         Write data header for projections file
- *
- */
+*    Projections::headerWwrite         Write data header for projections file
+*
+*/
 
 bool 
 Projections::headerWrite (fnetorderstream& fs)
@@ -171,7 +224,7 @@ Projections::headerWrite (fnetorderstream& fs)
   kuint16 _hour = m_hour;
   kuint16 _minute = m_minute;
   kuint16 _second = m_second;
-
+  
   kfloat64 _calcTime = m_calcTime;
   kfloat64 _rotStart = m_rotStart;
   kfloat64 _rotInc = m_rotInc;
@@ -180,11 +233,11 @@ Projections::headerWrite (fnetorderstream& fs)
   kfloat64 _phmLen = m_phmLen;
   kfloat64 _fieldOfView = m_fieldOfView;
   kfloat64 _focalLength = m_focalLength;
-
+  
   fs.seekp(0);
   if (! fs)
     return false;
-
+  
   fs.writeInt16 (_hsize);
   fs.writeInt16 (_signature);
   fs.writeInt32 (_nView);
@@ -206,21 +259,21 @@ Projections::headerWrite (fnetorderstream& fs)
   fs.writeInt16 (_second);
   fs.writeInt16 (_remarksize);
   fs.write (m_remark.c_str(), _remarksize);
-
+  
   m_headerSize = fs.tellp();
   _hsize = m_headerSize;
   fs.seekp(0);
   fs.writeInt16 (_hsize);
   if (! fs)
-      return false;
+    return false;
   
   return true;
 }
 
 /* NAME
- *    projections_read_header         Read data header for projections file
- *
- */
+*    projections_read_header         Read data header for projections file
+*
+*/
 bool
 Projections::headerRead (fnetorderstream& fs)
 {
@@ -230,8 +283,8 @@ Projections::headerRead (fnetorderstream& fs)
   
   fs.seekg(0);
   if (! fs)
-      return false;
-
+    return false;
+  
   fs.readInt16 (_hsize);
   fs.readInt16 (_signature);
   fs.readInt32 (_nView);
@@ -252,17 +305,17 @@ Projections::headerRead (fnetorderstream& fs)
   fs.readInt16 (_minute);
   fs.readInt16 (_second);
   fs.readInt16 (_remarksize);
-
+  
   if (! fs) {
-      sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
-      return false;
+    sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
+    return false;
   }
-
+  
   if (_signature != m_signature) {
     sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
     return false;
   }
-
+  
   char* pszRemarkStorage = new char [_remarksize+1];
   fs.read (pszRemarkStorage, _remarksize);
   if (! fs) {
@@ -272,11 +325,11 @@ Projections::headerRead (fnetorderstream& fs)
   pszRemarkStorage[_remarksize] = 0;
   m_remark = pszRemarkStorage;
   delete pszRemarkStorage;
-
+  
   off_t _hsizeread = fs.tellg();
   if (!fs || _hsizeread != _hsize) {
-      sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize);
-      return false;
+    sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize);
+    return false;
   }
   
   m_headerSize = _hsize;
@@ -297,12 +350,12 @@ Projections::headerRead (fnetorderstream& fs)
   m_hour = _hour;
   m_minute = _minute;
   m_second = _second;
-
+  
   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
   m_label.setLabelString (m_remark);
   m_label.setCalcTime (m_calcTime);
   m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second);
-
+  
   return true;
 }
 
@@ -322,21 +375,21 @@ Projections::read (const char* filename)
 #else
   frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary | std::ios::nocreate);
 #endif
-
+  
   if (fileRead.fail())
     return false;
-
+  
   if (! headerRead (fileRead))
     return false;
-
+  
   deleteProjData ();
   newProjData();
-
+  
   for (int i = 0; i < m_nView; i++) {
     if (! detarrayRead (fileRead, *m_projData[i], i))
       break;
   }
-
+  
   fileRead.close();
   return true;
 }
@@ -345,7 +398,7 @@ Projections::read (const char* filename)
 bool 
 Projections::copyViewData (const std::string& filename, std::ostream& os, int startView, int endView)
 {
-    return copyViewData (filename.c_str(), os, startView, endView);
+  return copyViewData (filename.c_str(), os, startView, endView);
 }
 
 bool 
@@ -354,56 +407,56 @@ Projections::copyViewData (const char* const filename, std::ostream& os, int sta
   frnetorderstream is (filename, ios::in | ios::binary);
   kuint16 sizeHeader, signature;
   kuint32 _nView, _nDet;
-
+  
   is.readInt16 (sizeHeader);
   is.readInt16 (signature);
   is.readInt32 (_nView);
   is.readInt32 (_nDet);
   int nView = _nView;
   int nDet = _nDet;
-
+  
   if (signature != m_signature) {
     sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
     return false;
   }
-
+  
   if (startView < 0)
-      startView = 0;
+    startView = 0;
   if (startView > nView - 1)
-      startView = nView;
+    startView = nView;
   if (endView < 0 || endView > nView - 1)
-      endView = nView - 1;
-
+    endView = nView - 1;
+  
   if (startView > endView) { // swap if start > end
-      int tempView = endView;
-      endView = startView;
-      startView = tempView;
+    int tempView = endView;
+    endView = startView;
+    startView = tempView;
   }
-
+  
   int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet);
   unsigned char* pViewData = new unsigned char [sizeView];
-
+  
   for (int i = startView; i <= endView; i++) {
-      is.seekg (sizeHeader + i * sizeView);
-      is.read (reinterpret_cast<char*>(pViewData), sizeView);
-      os.write (reinterpret_cast<char*>(pViewData), sizeView);
-      if (is.fail() || os.fail())
-         break;
+    is.seekg (sizeHeader + i * sizeView);
+    is.read (reinterpret_cast<char*>(pViewData), sizeView);
+    os.write (reinterpret_cast<char*>(pViewData), sizeView);
+    if (is.fail() || os.fail())
+      break;
   }
-
+  
   delete pViewData;
   if (is.fail()) 
     sys_error (ERR_FATAL, "Error reading projection file");
   if (os.fail()) 
     sys_error (ERR_FATAL, "Error writing projection file");
-    
+  
   return (! (is.fail() | os.fail()));
 }
 
 bool 
 Projections::copyHeader (const std::string& filename, std::ostream& os)
 {
-    return copyHeader (filename.c_str(), os);
+  return copyHeader (filename.c_str(), os);
 }
 
 bool
@@ -418,20 +471,20 @@ Projections::copyHeader (const char* const filename, std::ostream& os)
     sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
     return false;
   }
-
+  
   unsigned char* pHdrData = new unsigned char [sizeHeader];
   is.read (reinterpret_cast<char*>(pHdrData), sizeHeader);
   if (is.fail()) {
-      sys_error (ERR_FATAL, "Error reading header");
-      return false;
+    sys_error (ERR_FATAL, "Error reading header");
+    return false;
   }
-
+  
   os.write (reinterpret_cast<char*>(pHdrData), sizeHeader);
   if (os.fail()) {
-      sys_error (ERR_FATAL, "Error writing header");
-      return false;
+    sys_error (ERR_FATAL, "Error writing header");
+    return false;
   }
-
+  
   return true;
 }
 
@@ -450,32 +503,32 @@ Projections::write (const char* filename)
     sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
     return false;
   }
-
+  
   if (! headerWrite (fs))
-      return false;
-
+    return false;
+  
   if (m_projData != NULL) {
     for (int i = 0; i < m_nView; i++) {
       if (! detarrayWrite (fs, *m_projData[i], i))
-       break;
+        break;
     }
   }
   if (! fs)
     return false;
-
- fs.close();
-
+  
 fs.close();
+  
   return true;
 }
 
 /* NAME
- *   detarrayRead              Read a Detector Array structure from the disk
- *
- * SYNOPSIS
- *   detarrayRead (proj, darray, view_num)
- *   DETARRAY *darray          Detector array storage location to be filled
- *   int      view_num         View number to read
- */
+*   detarrayRead               Read a Detector Array structure from the disk
+*
+* SYNOPSIS
+*   detarrayRead (proj, darray, view_num)
+*   DETARRAY *darray           Detector array storage location to be filled
+*   int      view_num          View number to read
+*/
 
 bool
 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
@@ -487,39 +540,39 @@ Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int
   DetectorValue* detval_ptr = darray.detValues();  
   kfloat64 view_angle;
   kuint32 nDet;
-
+  
   fs.seekg (start_data);
-
+  
   fs.readFloat64 (view_angle);
   fs.readInt32 (nDet);
   darray.setViewAngle (view_angle);
   //  darray.setNDet ( nDet);
   
   for (unsigned int i = 0; i < nDet; i++) {
-      kfloat32 detval;
-      fs.readFloat32 (detval);
-      detval_ptr[i] = detval;
+    kfloat32 detval;
+    fs.readFloat32 (detval);
+    detval_ptr[i] = detval;
   }
   if (! fs)
     return false;
-
+  
   return true;
 }
 
 
 /* NAME
- *   detarrayWrite                     Write detector array data to the disk
- *
- * SYNOPSIS
- *   detarrayWrite (darray, view_num)
- *   DETARRAY *darray                  Detector array data to be written
- *   int      view_num                 View number to write
- *
- * DESCRIPTION
- *       This routine writes the detarray data from the disk sequentially to
- *    the file that was opened with open_projections().  Data is written in
- *    binary format.
- */
+*   detarrayWrite                      Write detector array data to the disk
+*
+* SYNOPSIS
+*   detarrayWrite (darray, view_num)
+*   DETARRAY *darray                   Detector array data to be written
+*   int      view_num                  View number to write
+*
+* DESCRIPTION
+*       This routine writes the detarray data from the disk sequentially to
+*    the file that was opened with open_projections().  Data is written in
+*    binary format.
+*/
 
 bool
 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
@@ -537,27 +590,27 @@ Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, co
     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
     return false;
   }
-
+  
   fs.writeFloat64 (view_angle);
   fs.writeInt32 (nDet);
-
+  
   for (unsigned int i = 0; i < nDet; i++) {
     kfloat32 detval = detval_ptr[i];
     fs.writeFloat32 (detval);
   }
-
+  
   if (! fs)
     return (false);
-
+  
   return true;
 }
 
 /* NAME
- *   printProjectionData                       Print projections data
- *
- * SYNOPSIS
- *   printProjectionData ()
- */
+*   printProjectionData                        Print projections data
+*
+* SYNOPSIS
+*   printProjectionData ()
+*/
 
 void
 Projections::printProjectionData ()
@@ -586,7 +639,7 @@ Projections::printProjectionData (int startView, int endView)
       printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
       DetectorValue* detval = m_projData[ir]->detValues();
       for (int id = 0; id < m_projData[ir]->nDet(); id++)
-       printf("%8.4f  ", detval[id]);
+        printf("%8.4f  ", detval[id]);
       printf("\n");
     }
   }
@@ -609,14 +662,191 @@ Projections::printScanInfo (std::ostringstream& os) const
 }
 
 
-bool Projections::convertPolar (ImageFile& rIF, int iInterpolation)
+bool 
+Projections::convertPolar (ImageFile& rIF, int iInterpolationID)
+{
+  unsigned int nx = rIF.nx();
+  unsigned int ny = rIF.ny();
+  ImageFileArray v = rIF.getArray();
+  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();
+
+  calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
+
+  std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+  for (unsigned int iView = 0; iView < m_nView; iView++) {
+    ppcDetValue[iView] = new std::complex<double> [m_nDet];
+    for (unsigned int iDet = 0; iDet < m_nDet; iDet++)
+      ppcDetValue[iView][iDet] = std::complex<double>(getDetectorArray (iView).detValues()[iDet], 0);
+  }
+
+  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)
+{
+  double xMin = -m_phmLen / 2;
+  double xMax = xMin + m_phmLen;
+  double yMin = -m_phmLen / 2;
+  double yMax = yMin + m_phmLen;
+  
+  double xInc = (xMax - xMin) / nx;    // size of cells
+  double yInc = (yMax - yMin) / ny;
+  
+  int iDetCenter = (m_nDet - 1) / 2;   // index refering to L=0 projection 
+
+  if (m_geometry != Scanner::GEOMETRY_PARALLEL) {
+    sys_error (ERR_WARNING, "convertPolar supports Parallel only");
+  }
+  
+  double x = xMin + xInc / 2;  // Rectang coords of center of pixel 
+  for (unsigned int ix = 0; ix < nx; x += xInc, ix++) {
+    double y = yMin + yInc / 2;
+    for (unsigned int iy = 0; iy < ny; y += yInc, iy++) {
+      double r = ::sqrt (x * x + y * y);
+      double phi = atan2 (y, x);
+
+      if (phi >= PI) {
+        phi -= PI;
+      } else if (phi < 0) {
+        phi += PI;
+      } else
+        r = -r;
+      
+      ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc;
+      ppdDet[ix][iy] = (r / m_detInc) + iDetCenter;
+    }
+  }
+}
+
+void
+Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
+     unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue,
+     double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, int iInterpolationID)
 {
-  return false;
+  for (unsigned int ix = 0; ix < ny; ix++) {
+    for (unsigned int iy = 0; iy < ny; iy++) {
+      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 (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) {
+          v[ix][iy] = ppcDetValue[iView][iDet].real();
+          if (vImag)
+            vImag[ix][iy] = ppcDetValue[iView][iDet].imag();
+        } else {
+          sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", 
+            ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
+          v[ix][iy] = 0;
+        }
+      } else if (iInterpolationID == POLAR_INTERP_BILINEAR) {
+        int iFloorView = ppdView[ix][iy];
+        double dFracView = ppdView[ix][iy] - iFloorView;
+        int iFloorDet = ppdDet[ix][iy];
+        double dFracDet = ppdDet[ix][iy] - iFloorDet;
+
+        if (iFloorDet >= 0 && iFloorView >= 0) { 
+          std::complex<double> v1 = ppcDetValue[iFloorView][iFloorDet];
+          std::complex<double> v2, v3, v4;
+          if (iFloorView < nView - 1)
+            v2 = ppcDetValue[iFloorView + 1][iFloorDet];
+          else 
+            v2 = v1;
+          if (iFloorDet < nDet - 1) 
+            v4 = ppcDetValue[iFloorView][iFloorDet+1];
+          else
+            v4 = v1;
+          if (iFloorView < nView - 1 && iFloorDet < nDet - 1)
+            v3 = ppcDetValue [iFloorView+1][iFloorDet+1];
+          else if (iFloorView < nView - 1)
+            v4 = v2;
+          else 
+            v4 = v1;
+          std::complex<double> vInterp = v1 + dFracView * (v2 - v1) + dFracDet * (v4 - v1);
+          v[ix][iy] = vInterp.real();
+          if (vImag)
+            vImag[ix][iy] = vInterp.imag();
+        } else {
+          sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", 
+            ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
+          v[ix][iy] = 0;
+          if (vImag)
+            vImag[ix][iy] = 0;
+        }
+      } else if (iInterpolationID == POLAR_INTERP_BICUBIC) {
+        v[ix][iy] =0;
+          if (vImag)
+            vImag[ix][iy] = 0;
+      }
+    }
+  }
 }
 
-bool Projections::convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad)
+bool 
+Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
 {
-  return false;
+  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;
 }
 
 
index c593b09..e10446f 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.22 2000/12/18 06:32:13 kevin Exp $
+**  $Id: scanner.cpp,v 1.23 2001/01/04 21:28:41 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
@@ -71,21 +71,21 @@ DetectorArray::~DetectorArray (void)
 
 
 /* NAME
- *   Scanner::Scanner          Construct a user specified detector structure
- *
- * SYNOPSIS
- *   Scanner (phm, nDet, nView, nSample)
- *   Phantom& phm              PHANTOM that we are making detector for
- *   int geomety                Geometry of detector
- *   int nDet                  Number of detector along detector array
- *   int nView                 Number of rotated views
- *   int nSample               Number of rays per detector
- */
+*   Scanner::Scanner           Construct a user specified detector structure
+*
+* SYNOPSIS
+*   Scanner (phm, nDet, nView, nSample)
+*   Phantom& phm               PHANTOM that we are making detector for
+*   int geomety                Geometry of detector
+*   int nDet                   Number of detector along detector array
+*   int nView                  Number of rotated views
+*   int nSample                Number of rays per detector
+*/
 
 Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen, const double dFocalLengthRatio, const double dFieldOfViewRatio)
 {
   m_phmLen = phm.maxAxisLength();      // maximal length along an axis
-
+  
   m_fail = false;
   m_idGeometry = convertGeometryNameToID (geometryName);
   if (m_idGeometry == GEOMETRY_INVALID) {
@@ -94,7 +94,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     m_failMessage += geometryName;
     return;
   }
-
+  
   if (nView < 1 || nDet < 1) {
     m_fail = true;
     m_failMessage = "nView & nDet must be greater than 0";
@@ -102,7 +102,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
   }
   if (nSample < 1)
     m_nSample = 1;
-
+  
   m_nDet     = nDet;
   m_nView    = nView;
   m_nSample  = nSample;
@@ -110,7 +110,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
   m_dFieldOfViewRatio = dFieldOfViewRatio;
   m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio;
   m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio;
-
+  
   m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
   m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
   m_rotLen  = rot_anglen;
@@ -145,7 +145,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     double dAngle = atan( dHalfSquare / dFocalPastPhm );
 #endif
     double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
-
+    
     m_detLen = dHalfDetLen * 2;
     m_detInc  = m_detLen / m_nDet;
     if (m_nDet % 2 == 0) // Adjust for Even number of detectors
@@ -181,14 +181,14 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     m_dAngularDetIncrement = m_detInc * 2; // Angular Position 2x gamma angle
     m_dAngularDetLen = m_detLen * 2;
     m_initPos.dAngularDet = -m_dAngularDetLen / 2;
-
+    
     m_initPos.angle = 0;
     m_initPos.xs1 = m_dXCenter;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;;
     m_initPos.xs2 = m_dXCenter;
     m_initPos.ys2 = m_dYCenter + m_dFocalLength;
   }
-
+  
   // Calculate incrementatal rotation matrix 
   GRFMTX_2D temp;
   xlat_mtx2 (m_rotmtxIncrement, -m_dXCenter, -m_dYCenter);
@@ -208,10 +208,10 @@ const char*
 Scanner::convertGeometryIDToName (const int geomID)
 {
   const char *name = "";
-
+  
   if (geomID >= 0 && geomID < s_iGeometryCount)
-      return (s_aszGeometryName[geomID]);
-
+    return (s_aszGeometryName[geomID]);
+  
   return (name);
 }
 
@@ -219,37 +219,37 @@ const char*
 Scanner::convertGeometryIDToTitle (const int geomID)
 {
   const char *title = "";
-
+  
   if (geomID >= 0 && geomID < s_iGeometryCount)
-      return (s_aszGeometryName[geomID]);
-
+    return (s_aszGeometryName[geomID]);
+  
   return (title);
 }
-      
+
 int
 Scanner::convertGeometryNameToID (const char* const geomName) 
 {
   int id = GEOMETRY_INVALID;
-
+  
   for (int i = 0; i < s_iGeometryCount; i++)
-      if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) {
-         id = i;
-         break;
-      }
-
-  return (id);
+    if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) {
+      id = i;
+      break;
+    }
+    
+    return (id);
 }
-  
+
 
 /* NAME
- *   collectProjections                Calculate projections for a Phantom
- *
- * SYNOPSIS
- *   collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace)
- *   Projectrions& proj      Projection storage
- *   Phantom& phm           Phantom for which we collect projections
- *   bool bStoreViewPos      TRUE then storage proj at normal view position
- *   int trace              Trace level
+*   collectProjections         Calculate projections for a Phantom
+*
+* SYNOPSIS
+*   collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace)
+*   Projectrions& proj      Projection storage
+*   Phantom& phm            Phantom for which we collect projections
+*   bool bStoreViewPos      TRUE then storage proj at normal view position
+*   int trace               Trace level
 */
 
 
@@ -264,7 +264,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
 {
   m_trace = trace;
   double start_angle = iStartView * proj.rotInc();
-
+  
   // Calculate initial rotation matrix 
   GRFMTX_2D rotmtx_initial, temp;
   xlat_mtx2 (rotmtx_initial, -m_dXCenter, -m_dYCenter);
@@ -272,15 +272,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
   xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
-
+  
   double xd1=0, yd1=0, xd2=0, yd2=0;
   if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
-      xd1 = m_initPos.xd1;
-      yd1 = m_initPos.yd1;
-      xd2 = m_initPos.xd2;
-      yd2 = m_initPos.yd2;
-      xform_mtx2 (rotmtx_initial, xd1, yd1);      // rotate detector endpoints 
-      xform_mtx2 (rotmtx_initial, xd2, yd2);      // to initial view_angle 
+    xd1 = m_initPos.xd1;
+    yd1 = m_initPos.yd1;
+    xd2 = m_initPos.xd2;
+    yd2 = m_initPos.yd2;
+    xform_mtx2 (rotmtx_initial, xd1, yd1);      // rotate detector endpoints 
+    xform_mtx2 (rotmtx_initial, xd2, yd2);      // to initial view_angle 
   }
   
   double xs1 = m_initPos.xs1;
@@ -289,18 +289,18 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   double ys2 = m_initPos.ys2;
   xform_mtx2 (rotmtx_initial, xs1, ys1);      // rotate source endpoints to
   xform_mtx2 (rotmtx_initial, xs2, ys2);      // initial view angle
-
+  
   int iView;
   double viewAngle;
   for (iView = 0, viewAngle = start_angle;  iView < iNumViews; iView++, viewAngle += proj.rotInc()) {
     int iStoragePosition = iView;
     if (bStoreAtViewPosition)
       iStoragePosition += iStartView;
-
+    
     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
-
+    
 #ifdef HAVE_SGP 
-  if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
+    if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
       m_pSGP = pSGP;
       double dWindowSize = dmax (m_detLen, m_dFocalLength * 2) * SQRT2;
       double dHalfWindowSize = dWindowSize / 2;
@@ -309,53 +309,53 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       m_dYMinWin = m_dYCenter - dHalfWindowSize;
       m_dYMaxWin = m_dYCenter + dHalfWindowSize;
       double dHalfPhmLen = m_phmLen /  2;
-
-    m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
-    m_pSGP->setRasterOp (RO_COPY);
-    m_pSGP->setColor (C_RED);
-    m_pSGP->moveAbs (0., 0.);
-    m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen);
-    m_pSGP->moveAbs (0., 0.);
-    m_pSGP->drawCircle (m_dFocalLength);
-    m_pSGP->setColor (C_BLUE);
-    phm.draw (*m_pSGP);
-    m_dTextHeight = m_pSGP->getCharHeight ();
-
-    traceShowParam ("Phantom:",       "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str());
-    traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry));
-    traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio);
-    traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio);
-    traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet());
-    traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView());
-    traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample);
-    
-    m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
-  }
+      
+      m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
+      m_pSGP->setRasterOp (RO_COPY);
+      m_pSGP->setColor (C_RED);
+      m_pSGP->moveAbs (0., 0.);
+      m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen);
+      m_pSGP->moveAbs (0., 0.);
+      m_pSGP->drawCircle (m_dFocalLength);
+      m_pSGP->setColor (C_BLUE);
+      phm.draw (*m_pSGP);
+      m_dTextHeight = m_pSGP->getCharHeight ();
+      
+      traceShowParam ("Phantom:",       "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str());
+      traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry));
+      traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio);
+      traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio);
+      traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet());
+      traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView());
+      traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample);
+      
+      m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
+    }
 #endif
-
+    
 #ifdef HAVE_SGP
     if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
       m_pSGP->setColor (C_BLACK);
       m_pSGP->setPenWidth (2);
       if (m_idGeometry == GEOMETRY_PARALLEL) {
-       m_pSGP->moveAbs (xs1, ys1);
-       m_pSGP->lineAbs (xs2, ys2);
-       m_pSGP->moveAbs (xd1, yd1);
-       m_pSGP->lineAbs (xd2, yd2);
+        m_pSGP->moveAbs (xs1, ys1);
+        m_pSGP->lineAbs (xs2, ys2);
+        m_pSGP->moveAbs (xd1, yd1);
+        m_pSGP->lineAbs (xd2, yd2);
       } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {        
-       m_pSGP->setPenWidth (4);
-       m_pSGP->moveAbs (xs1, ys1);
-       m_pSGP->lineAbs (xs2, ys2);
-       m_pSGP->setPenWidth (2);
-       m_pSGP->moveAbs (xd1, yd1);
-       m_pSGP->lineAbs (xd2, yd2);
+        m_pSGP->setPenWidth (4);
+        m_pSGP->moveAbs (xs1, ys1);
+        m_pSGP->lineAbs (xs2, ys2);
+        m_pSGP->setPenWidth (2);
+        m_pSGP->moveAbs (xd1, yd1);
+        m_pSGP->lineAbs (xd2, yd2);
       } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-       m_pSGP->setPenWidth (4);
-       m_pSGP->moveAbs (xs1, ys1);
-       m_pSGP->lineAbs (xs2, ys2);
-       m_pSGP->setPenWidth (2);
-       m_pSGP->moveAbs (0., 0.);
-       m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
+        m_pSGP->setPenWidth (4);
+        m_pSGP->moveAbs (xs1, ys1);
+        m_pSGP->lineAbs (xs2, ys2);
+        m_pSGP->setPenWidth (2);
+        m_pSGP->moveAbs (0., 0.);
+        m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
       }
       m_pSGP->setPenWidth (1);
     }
@@ -363,11 +363,11 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast<double>(m_nView) * 100.);
 #endif
     if (m_trace == Trace::TRACE_CONSOLE)
-               std::cout << "Current View: " << iView+iStartView << std::endl;
-
+      std::cout << "Current View: " << iView+iStartView << std::endl;
+    
     projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI);
     detArray.setViewAngle (viewAngle);
-      
+    
 #ifdef HAVE_SGP
     if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
       //       rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta);
@@ -376,124 +376,124 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
     xform_mtx2 (m_rotmtxIncrement, xs1, ys1);
     xform_mtx2 (m_rotmtxIncrement, xs2, ys2);
     if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
-       xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints 
-       xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
+      xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints 
+      xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
     }
   } /* for each iView */
 }
 
 
 /* NAME
- *    rayview                  Calculate raysums for a view at any angle
- *
- * SYNOPSIS
- *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
- *    Phantom& phm             Phantom to scan
- *    DETARRAY *detArray               Storage of values for detector array
- *    Scanner& det             Scanner parameters
- *    double xd1, yd1, xd2, yd2        Beginning & ending detector positions
- *    double xs1, ys1, xs2, ys2        Beginning & ending source positions
- *
- * RAY POSITIONING
- *         For each detector, have there are a variable number of rays traced.
- *     The source of each ray is the center of the source x-ray cell. The
- *     detector positions are equally spaced within the cell
- *
- *         The increments between rays are calculated so that the cells start
- *     at the beginning of a detector cell and they end on the endpoint
- *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
- *         The exception to this is if there is only one ray per detector.
- *     In that case, the detector position is the center of the detector cell.
- */
+*    rayview                   Calculate raysums for a view at any angle
+*
+* SYNOPSIS
+*    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
+*    Phantom& phm              Phantom to scan
+*    DETARRAY *detArray                Storage of values for detector array
+*    Scanner& det              Scanner parameters
+*    double xd1, yd1, xd2, yd2 Beginning & ending detector positions
+*    double xs1, ys1, xs2, ys2 Beginning & ending source positions
+*
+* RAY POSITIONING
+*         For each detector, have there are a variable number of rays traced.
+*     The source of each ray is the center of the source x-ray cell. The
+*     detector positions are equally spaced within the cell
+*
+*         The increments between rays are calculated so that the cells start
+*     at the beginning of a detector cell and they end on the endpoint
+*     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
+*         The exception to this is if there is only one ray per detector.
+*     In that case, the detector position is the center of the detector cell.
+*/
 
 void 
 Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, const double dDetAngle)
 {
-
+  
   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords 
   double sdy = (ys2 - ys1) / detArray.nDet();  // between source
   double xs_maj = xs1 + (sdx / 2);     // put ray source in center of cell 
   double ys_maj = ys1 + (sdy / 2);
-
+  
   double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0;
   double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0;
   if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-      dAngleInc = m_dAngularDetIncrement;
-      dAngleSampleInc = dAngleInc / m_nSample;
-      dAngleSampleOffset = dAngleSampleInc / 2;
-      dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset;
+    dAngleInc = m_dAngularDetIncrement;
+    dAngleSampleInc = dAngleInc / m_nSample;
+    dAngleSampleOffset = dAngleSampleInc / 2;
+    dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset;
   } else {
-      ddx = (xd2 - xd1) / detArray.nDet();  // change in coords 
-      ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
-      ddx2 = ddx / m_nSample;  // Incr. between rays with detector cell
-      ddy2 = ddy / m_nSample;  // Doesn't include detector endpoints 
-      ddx2_ofs = ddx2 / 2;    // offset of 1st ray from start of detector cell
-      ddy2_ofs = ddy2 / 2;
-      
-      xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
-      yd_maj = yd1 + ddy2_ofs;
+    ddx = (xd2 - xd1) / detArray.nDet();  // change in coords 
+    ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
+    ddx2 = ddx / m_nSample;    // Incr. between rays with detector cell
+    ddy2 = ddy / m_nSample;  // Doesn't include detector endpoints 
+    ddx2_ofs = ddx2 / 2;    // offset of 1st ray from start of detector cell
+    ddy2_ofs = ddy2 / 2;
+    
+    xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
+    yd_maj = yd1 + ddy2_ofs;
   }
-
+  
   DetectorValue* detval = detArray.detValues();
-
+  
   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
     for (int d = 0; d < detArray.nDet(); d++)
       if (detArray.nDet() / 2 == d && (d % 2) == 1)
-       detval[d] = 1;
+        detval[d] = 1;
       else
-       detval[d] = 0;
+        detval[d] = 0;
   } else {
-      for (int d = 0; d < detArray.nDet(); d++) {
+    for (int d = 0; d < detArray.nDet(); d++) {
       double xs = xs_maj;
       double ys = ys_maj;
       double xd=0, yd=0, dAngle=0;
       if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-         dAngle = dAngleMajor;
+        dAngle = dAngleMajor;
       } else {
-         xd = xd_maj;
-         yd = yd_maj;
+        xd = xd_maj;
+        yd = yd_maj;
       }
       double sum = 0.0;
       for (unsigned int i = 0; i < m_nSample; i++) {
-       if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-           xd = m_dFocalLength * cos (dAngle);
-           yd = m_dFocalLength * sin (dAngle);
-       }
-
+        if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
+          xd = m_dFocalLength * cos (dAngle);
+          yd = m_dFocalLength * sin (dAngle);
+        }
+        
 #ifdef HAVE_SGP
-       if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
-         m_pSGP->setColor (C_YELLOW);
-         m_pSGP->setRasterOp (RO_AND);
-         m_pSGP->moveAbs (xs, ys);
-         m_pSGP->lineAbs (xd, yd);
-       }
+        if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
+          m_pSGP->setColor (C_YELLOW);
+          m_pSGP->setRasterOp (RO_AND);
+          m_pSGP->moveAbs (xs, ys);
+          m_pSGP->lineAbs (xd, yd);
+        }
 #endif
-
-       sum += projectSingleLine (phm, xd, yd, xs, ys);
-             
+        
+        sum += projectSingleLine (phm, xd, yd, xs, ys);
+        
 #ifdef HAVE_SGP
-       //      if (m_trace >= Trace::TRACE_CLIPPING) {
-       //        traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, "        ");
-       //        traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
-       //      }
+        //     if (m_trace >= Trace::TRACE_CLIPPING) {
+        //       traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, "        ");
+        //       traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
+        //     }
 #endif
-       if (m_idGeometry == GEOMETRY_EQUIANGULAR)
-           dAngle += dAngleSampleInc;
-       else {
-           xd += ddx2;
-           yd += ddy2;
-       }
+        if (m_idGeometry == GEOMETRY_EQUIANGULAR)
+          dAngle += dAngleSampleInc;
+        else {
+          xd += ddx2;
+          yd += ddy2;
+        }
       } // for each sample in detector
-
+      
       detval[d] = sum / m_nSample;
       xs_maj += sdx;
       ys_maj += sdy;
       if (m_idGeometry == GEOMETRY_EQUIANGULAR)
-         dAngleMajor += dAngleInc;
-       else {
-           xd_maj += ddx;
-           yd_maj += ddy;
-       }
+        dAngleMajor += dAngleInc;
+      else {
+        xd_maj += ddx;
+        yd_maj += ddy;
+      }
     } /* for each detector */
   } /* if not unit pulse */
 }
@@ -529,9 +529,9 @@ void
 Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char *fmt, int row, int color, va_list args)
 {  
   char szValue[256];
-
+  
   vsnprintf (szValue, sizeof(szValue), fmt, args);
-
+  
 #ifdef HAVE_SGP
   if (m_pSGP) {
     m_pSGP->setRasterOp (iRasterOp);
@@ -544,24 +544,24 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char
     m_pSGP->drawText (szValue);
   } else 
 #endif
-    {
-      cio_put_str (szLabel);
-      cio_put_str (szValue);
-      cio_put_str ("\n");
-    }
+  {
+    cio_put_str (szLabel);
+    cio_put_str (szValue);
+    cio_put_str ("\n");
+  }
 }
 
 
 
 /* NAME
- *    projectSingleLine                        INTERNAL: Calculates raysum along a line for a Phantom
- *
- * SYNOPSIS
- *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
- *    double rsum              Ray sum of Phantom along given line
- *    Phantom& phm;            Phantom from which to calculate raysum
- *    double *x1, *y1, *x2, y2 Endpoints of ray path (in Phantom coords)
- */
+*    projectSingleLine                 INTERNAL: Calculates raysum along a line for a Phantom
+*
+* SYNOPSIS
+*    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
+*    double rsum               Ray sum of Phantom along given line
+*    Phantom& phm;             Phantom from which to calculate raysum
+*    double *x1, *y1, *x2, y2  Endpoints of ray path (in Phantom coords)
+*/
 
 double 
 Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2)
@@ -570,20 +570,20 @@ Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1
   double rsum = 0.0;
   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
     rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2);
-
+  
   return (rsum);
 }
 
 
 /* NAME
- *   pelem_ray_attenuation             Calculate raysum of an pelem along one line
- *
- * SYNOPSIS
- *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
- *   double rsum               Computed raysum
- *   PhantomElement& pelem             Pelem to scan
- *   double x1, y1, x2, y2     Endpoints of raysum line
- */
+*   pelem_ray_attenuation              Calculate raysum of an pelem along one line
+*
+* SYNOPSIS
+*   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
+*   double rsum                Computed raysum
+*   PhantomElement& pelem              Pelem to scan
+*   double x1, y1, x2, y2      Endpoints of raysum line
+*/
 
 double 
 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2)
@@ -593,7 +593,7 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double
       cio_tone (1000., 0.05);
     return (0.0);
   }
-
+  
 #ifdef HAVE_SGP
   if (m_pSGP && m_trace == Trace::TRACE_CLIPPING) {
     m_pSGP->setRasterOp (RO_XOR);
@@ -605,7 +605,7 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double
     m_pSGP->setRasterOp (RO_SET);
   }
 #endif
-
+  
   double len = lineLength (x1, y1, x2, y2);
   return (len * pelem.atten());
 }
index 27d403b..5ff72f1 100644 (file)
@@ -3,41 +3,44 @@
 <pre>
 <h1>Build Log</h1>
 <h3>
---------------------Configuration: ctsim - Win32 Release--------------------
+--------------------Configuration: ctsim - Win32 Debug--------------------
 </h3>
 <h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2E.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A0.tmp" with contents
 [
-/nologo /G6 /MT /W3 /GR /GX /O2 /I "." /I "..\..\include" /I "..\..\getopt" /I "..\..\..\lpng108" /I "..\..\..\zlib" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2\include" /D "NDEBUG" /D "__WXWIN__" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.0.0beta1\" /FR"Release/" /Fp"Release/ctsim.pch" /YX /Fo"Release/" /Fd"Release/" /FD /c 
-"C:\ctsim\src\views.cpp"
+/nologo /G6 /MTd /W3 /Gm /GR /GX /ZI /Od /I "\wx2\include" /I "." /I "..\..\include" /I "..\..\getopt" /I "..\..\..\lpng108" /I "..\..\..\zlib" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /D VERSION=\"2.5.0\" /D "_DEBUG" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.0.0beta1\" /FR"Debug/" /Fp"Debug/ctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c 
+"C:\ctsim\src\dialogs.cpp"
 ]
-Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2E.tmp" 
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2F.tmp" with contents
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A0.tmp" 
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A1.tmp" with contents
 [
-kernel32.lib user32.lib wsock32.lib comctl32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib libpng.lib zlib.lib /nologo /subsystem:windows /incremental:no /pdb:"Release/ctsim.pdb" /machine:I386 /out:"Release/ctsim.exe" /libpath:"..\..\..\lpng108\msvc\win32\libpng\lib" /libpath:"..\..\..\lpng108\msvc\win32\zlib\lib" 
-.\Release\ctsim.obj
-.\Release\dialogs.obj
-.\Release\dlgprojections.obj
-.\Release\dlgreconstruct.obj
-.\Release\docs.obj
-.\Release\views.obj
-\ctsim\msvc\libctsim\Release\libctsim.lib
-"\fftw-2.1.3\Win32\FFTW2st\Release\FFTW2st.lib"
-"\fftw-2.1.3\Win32\RFFTW2st\Release\RFFTW2st.lib"
-\wx2\lib\wx.lib
+comctl32.lib winmm.lib rpcrt4.lib ws2_32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ../libctsim/Debug/libctsim.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ..\..\..\lpng108\msvc\win32\libpng\lib_dbg\libpng.lib ..\..\..\lpng108\msvc\win32\zlib\lib_dbg\zlib.lib libcmtd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib ../../../wx2/lib/wxd.lib /nologo /subsystem:windows /incremental:yes /pdb:"Debug/ctsim.pdb" /debug /machine:I386 /nodefaultlib:"libcd.lib" /nodefaultlib:"libcid.lib" /nodefaultlib:"msvcrtd.lib" /out:"Debug/ctsim.exe" /pdbtype:sept /libpath:"..\..\..\lpng108\msvc\win32\libpng\lib" /libpath:"..\..\..\lpng108\msvc\win32\zlib\lib" 
+.\Debug\ctsim.obj
+.\Debug\dialogs.obj
+.\Debug\dlgprojections.obj
+.\Debug\dlgreconstruct.obj
+.\Debug\docs.obj
+.\Debug\views.obj
+.\Debug\wx.res
+\ctsim\msvc\libctsim\Debug\libctsim.lib
+"\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib"
+"\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib"
+\wx2\lib\wxd.lib
 ]
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2F.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A1.tmp"
 <h3>Output Window</h3>
 Compiling...
-views.cpp
+dialogs.cpp
 Linking...
-LINK : warning LNK4089: all references to "ADVAPI32.dll" discarded by /OPT:REF
-LINK : warning LNK4089: all references to "WSOCK32.dll" discarded by /OPT:REF
+Creating command line "bscmake.exe /nologo /o"Debug/ctsim.bsc"  .\Debug\ctsim.sbr .\Debug\dialogs.sbr .\Debug\dlgprojections.sbr .\Debug\dlgreconstruct.sbr .\Debug\docs.sbr .\Debug\views.sbr"
+Creating browse info file...
+BSCMAKE: warning BK4503 : minor error in .SBR file '.\Debug\dialogs.sbr' ignored
+<h3>Output Window</h3>
 
 
 
 <h3>Results</h3>
-ctsim.exe - 0 error(s), 2 warning(s)
+ctsim.exe - 0 error(s), 1 warning(s)
 </pre>
 </body>
 </html>
index 31d02a4..e1055ba 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctsim.h,v 1.20 2001/01/03 22:00:46 kevin Exp $
+**  $Id: ctsim.h,v 1.21 2001/01/04 21:28:41 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
@@ -163,7 +163,6 @@ enum {
   MAINMENU_FILE_EXIT,
 
   IFMENU_FILE_PROPERTIES,
-  PHMMENU_FILE_PROPERTIES,
 
   PJMENU_FILE_PROPERTIES,
   PJMENU_RECONSTRUCT_FBP,
@@ -179,6 +178,7 @@ enum {
 
   IFMENU_VIEW_SCALE_AUTO,
   IFMENU_VIEW_SCALE_MINMAX,
+  IFMENU_VIEW_SCALE_FULL,
 
        IFMENU_COMPARE_IMAGES,
        IFMENU_COMPARE_ROW,
@@ -203,11 +203,13 @@ enum {
   IFMENU_FILTER_SHUFFLENATURALTOFOURIERORDER,
   IFMENU_FILTER_SHUFFLEFOURIERTONATURALORDER,
 
+  PHMMENU_FILE_PROPERTIES,
   PHMMENU_PROCESS_RASTERIZE,
   PHMMENU_PROCESS_PROJECTIONS,
 
        PLOTMENU_VIEW_SCALE_MINMAX,
        PLOTMENU_VIEW_SCALE_AUTO,
+       PLOTMENU_VIEW_SCALE_FULL,
   MAINMENU_WINDOW_BASE,
 };
 
index f0ece54..b5b32a5 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.cpp,v 1.24 2001/01/02 16:02:13 kevin Exp $
+**  $Id: dialogs.cpp,v 1.25 2001/01/04 21:28:41 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
@@ -84,7 +84,7 @@ StringValueAndTitleListBox::StringValueAndTitleListBox (wxDialog* pParent, int n
 };
 
 const char*
-StringValueAndTitleListBox::getSelectionStringValue (void) const
+StringValueAndTitleListBox::getSelectionStringValue () const
 {
   return m_ppszValues[GetSelection()];
 }
@@ -125,7 +125,7 @@ DialogGetPhantom::DialogGetPhantom (wxFrame* pParent, int iDefaultPhantom)
 }
 
 const char*
-DialogGetPhantom::getPhantom(void)
+DialogGetPhantom::getPhantom()
 {
   return m_pListBoxPhantom->getSelectionStringValue();
 }
@@ -183,7 +183,7 @@ DialogGetComparisonImage::DialogGetComparisonImage (wxFrame* pParent, const char
 }
 
 ImageFileDocument*
-DialogGetComparisonImage::getImageFileDocument(void)
+DialogGetComparisonImage::getImageFileDocument()
 {
   return m_rVecIF[ m_pListBoxImageChoices->GetSelection() ];
 }
@@ -241,12 +241,12 @@ DialogGetMinMax::DialogGetMinMax (wxFrame* pParent, const char* const pszTitle,
   pTopSizer->SetSizeHints (this);
 }
 
-DialogGetMinMax::~DialogGetMinMax (void)
+DialogGetMinMax::~DialogGetMinMax ()
 {
 }
 
 double
-DialogGetMinMax::getMinimum (void)
+DialogGetMinMax::getMinimum ()
 {
   wxString strCtrl = m_pTextCtrlMin->GetValue();
   double dValue;
@@ -257,7 +257,7 @@ DialogGetMinMax::getMinimum (void)
 }
 
 double
-DialogGetMinMax::getMaximum (void)
+DialogGetMinMax::getMaximum ()
 {
   wxString strCtrl = m_pTextCtrlMax->GetValue();
   double dValue;
@@ -401,13 +401,13 @@ DialogGetRasterParameters::DialogGetRasterParameters (wxFrame* pParent, int iDef
   pTopSizer->SetSizeHints (this);
 }
 
-DialogGetRasterParameters::~DialogGetRasterParameters (void)
+DialogGetRasterParameters::~DialogGetRasterParameters ()
 {
 }
 
 
 unsigned int
-DialogGetRasterParameters::getXSize (void)
+DialogGetRasterParameters::getXSize ()
 {
   wxString strCtrl = m_pTextCtrlXSize->GetValue();
   unsigned long lValue;
@@ -418,7 +418,7 @@ DialogGetRasterParameters::getXSize (void)
 }
 
 unsigned int
-DialogGetRasterParameters::getYSize (void)
+DialogGetRasterParameters::getYSize ()
 {
   wxString strCtrl = m_pTextCtrlYSize->GetValue();
   unsigned long lValue;
@@ -430,7 +430,7 @@ DialogGetRasterParameters::getYSize (void)
 
 
 unsigned int
-DialogGetRasterParameters::getNSamples (void)
+DialogGetRasterParameters::getNSamples ()
 {
   wxString strCtrl = m_pTextCtrlNSamples->GetValue();
   unsigned long lValue;
@@ -529,13 +529,13 @@ DialogGetProjectionParameters::DialogGetProjectionParameters (wxFrame* pParent,
   pTopSizer->SetSizeHints (this);
 }
 
-DialogGetProjectionParameters::~DialogGetProjectionParameters (void)
+DialogGetProjectionParameters::~DialogGetProjectionParameters ()
 {
 }
 
 
 unsigned int
-DialogGetProjectionParameters::getNDet (void)
+DialogGetProjectionParameters::getNDet ()
 {
   wxString strCtrl = m_pTextCtrlNDet->GetValue();
   unsigned long lValue;
@@ -546,7 +546,7 @@ DialogGetProjectionParameters::getNDet (void)
 }
 
 unsigned int
-DialogGetProjectionParameters::getNView (void)
+DialogGetProjectionParameters::getNView ()
 {
   wxString strCtrl = m_pTextCtrlNView->GetValue();
   unsigned long lValue;
@@ -558,7 +558,7 @@ DialogGetProjectionParameters::getNView (void)
 
 
 unsigned int
-DialogGetProjectionParameters::getNSamples (void)
+DialogGetProjectionParameters::getNSamples ()
 {
   wxString strCtrl = m_pTextCtrlNSamples->GetValue();
   unsigned long lValue;
@@ -569,7 +569,7 @@ DialogGetProjectionParameters::getNSamples (void)
 }
 
 double
-DialogGetProjectionParameters::getRotAngle (void)
+DialogGetProjectionParameters::getRotAngle ()
 {
   wxString strCtrl = m_pTextCtrlRotAngle->GetValue();
   double dValue;
@@ -580,7 +580,7 @@ DialogGetProjectionParameters::getRotAngle (void)
 }
 
 double
-DialogGetProjectionParameters::getFocalLengthRatio (void)
+DialogGetProjectionParameters::getFocalLengthRatio ()
 {
   wxString strCtrl = m_pTextCtrlFocalLength->GetValue();
   double dValue;
@@ -591,7 +591,7 @@ DialogGetProjectionParameters::getFocalLengthRatio (void)
 }
 
 double
-DialogGetProjectionParameters::getFieldOfViewRatio (void)
+DialogGetProjectionParameters::getFieldOfViewRatio ()
 {
   wxString strCtrl = m_pTextCtrlFieldOfView->GetValue();
   double dValue;
@@ -602,13 +602,13 @@ DialogGetProjectionParameters::getFieldOfViewRatio (void)
 }
 
 const char*
-DialogGetProjectionParameters::getGeometry (void)
+DialogGetProjectionParameters::getGeometry ()
 {
   return m_pListBoxGeometry->getSelectionStringValue();
 }
 
 int
-DialogGetProjectionParameters::getTrace (void)
+DialogGetProjectionParameters::getTrace ()
 {
   return Trace::convertTraceNameToID(m_pListBoxTrace->getSelectionStringValue());
 }
@@ -713,13 +713,13 @@ DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* p
   pTopSizer->SetSizeHints (this);
 }
 
-DialogGetReconstructionParameters::~DialogGetReconstructionParameters (void)
+DialogGetReconstructionParameters::~DialogGetReconstructionParameters ()
 {
 }
 
 
 unsigned int
-DialogGetReconstructionParameters::getXSize (void)
+DialogGetReconstructionParameters::getXSize ()
 {
   wxString strCtrl = m_pTextCtrlXSize->GetValue();
   unsigned long lValue;
@@ -730,7 +730,7 @@ DialogGetReconstructionParameters::getXSize (void)
 }
 
 unsigned int
-DialogGetReconstructionParameters::getYSize (void)
+DialogGetReconstructionParameters::getYSize ()
 {
   wxString strCtrl = m_pTextCtrlYSize->GetValue();
   unsigned long lValue;
@@ -741,7 +741,7 @@ DialogGetReconstructionParameters::getYSize (void)
 }
 
 unsigned int
-DialogGetReconstructionParameters::getZeropad (void)
+DialogGetReconstructionParameters::getZeropad ()
 {
   wxString strCtrl = m_pTextCtrlZeropad->GetValue();
   unsigned long lValue;
@@ -753,7 +753,7 @@ DialogGetReconstructionParameters::getZeropad (void)
 
 
 unsigned int
-DialogGetReconstructionParameters::getInterpParam (void)
+DialogGetReconstructionParameters::getInterpParam ()
 {
   wxString strCtrl = m_pTextCtrlInterpParam->GetValue();
   unsigned long lValue;
@@ -764,7 +764,7 @@ DialogGetReconstructionParameters::getInterpParam (void)
 }
 
 double
-DialogGetReconstructionParameters::getFilterParam (void)
+DialogGetReconstructionParameters::getFilterParam ()
 {
   wxString strCtrl = m_pTextCtrlFilterParam->GetValue();
   double dValue;
@@ -775,25 +775,25 @@ DialogGetReconstructionParameters::getFilterParam (void)
 }
 
 const char*
-DialogGetReconstructionParameters::getFilterName (void)
+DialogGetReconstructionParameters::getFilterName ()
 {
   return m_pListBoxFilter->getSelectionStringValue();
 }
 
 const char*
-DialogGetReconstructionParameters::getFilterMethodName (void)
+DialogGetReconstructionParameters::getFilterMethodName ()
 {
   return m_pListBoxFilterMethod->getSelectionStringValue();
 }
 
 const char*
-DialogGetReconstructionParameters::getInterpName (void)
+DialogGetReconstructionParameters::getInterpName ()
 {
   return m_pListBoxInterp->getSelectionStringValue();
 }
 
 int
-DialogGetReconstructionParameters::getTrace (void)
+DialogGetReconstructionParameters::getTrace ()
 {
   int iTrace = 0;
   if (strcmp("full", m_pListBoxTrace->getSelectionStringValue()) == 0)
@@ -802,13 +802,13 @@ DialogGetReconstructionParameters::getTrace (void)
 }
 
 const char*
-DialogGetReconstructionParameters::getBackprojectName (void)
+DialogGetReconstructionParameters::getBackprojectName ()
 {
   return m_pListBoxBackproject->getSelectionStringValue();
 }
 
 const char*
-DialogGetReconstructionParameters::getFilterGenerationName (void)
+DialogGetReconstructionParameters::getFilterGenerationName ()
 {
   return m_pListBoxFilterGeneration->getSelectionStringValue();
 }
@@ -821,6 +821,7 @@ DialogGetReconstructionParameters::getFilterGenerationName (void)
 /////////////////////////////////////////////////////////////////////
 
 
+
 DialogGetFilterParameters::DialogGetFilterParameters (wxFrame* pParent, int iDefaultXSize, int iDefaultYSize, int iDefaultFilterID, double dDefaultFilterParam,  double dDefaultBandwidth, int iDefaultDomainID, double dDefaultInputScale, double dDefaultOutputScale)
 : wxDialog (pParent, -1, "Set Filter Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)
 {
@@ -892,13 +893,13 @@ DialogGetFilterParameters::DialogGetFilterParameters (wxFrame* pParent, int iDef
   pTopSizer->SetSizeHints (this);
 }
 
-DialogGetFilterParameters::~DialogGetFilterParameters (void)
+DialogGetFilterParameters::~DialogGetFilterParameters ()
 {
 }
 
 
 unsigned int
-DialogGetFilterParameters::getXSize (void)
+DialogGetFilterParameters::getXSize ()
 {
   wxString strCtrl = m_pTextCtrlXSize->GetValue();
   unsigned long lValue;
@@ -909,7 +910,7 @@ DialogGetFilterParameters::getXSize (void)
 }
 
 unsigned int
-DialogGetFilterParameters::getYSize (void)
+DialogGetFilterParameters::getYSize ()
 {
   wxString strCtrl = m_pTextCtrlYSize->GetValue();
   unsigned long lValue;
@@ -920,7 +921,7 @@ DialogGetFilterParameters::getYSize (void)
 }
 
 double
-DialogGetFilterParameters::getBandwidth (void)
+DialogGetFilterParameters::getBandwidth ()
 {
   wxString strCtrl = m_pTextCtrlBandwidth->GetValue();
   double dValue;
@@ -931,7 +932,7 @@ DialogGetFilterParameters::getBandwidth (void)
 }
 
 double
-DialogGetFilterParameters::getFilterParam (void)
+DialogGetFilterParameters::getFilterParam ()
 {
   wxString strCtrl = m_pTextCtrlFilterParam->GetValue();
   double dValue;
@@ -942,7 +943,7 @@ DialogGetFilterParameters::getFilterParam (void)
 }
 
 double
-DialogGetFilterParameters::getInputScale (void)
+DialogGetFilterParameters::getInputScale ()
 {
   wxString strCtrl = m_pTextCtrlInputScale->GetValue();
   double dValue;
@@ -953,7 +954,7 @@ DialogGetFilterParameters::getInputScale (void)
 }
 
 double
-DialogGetFilterParameters::getOutputScale (void)
+DialogGetFilterParameters::getOutputScale ()
 {
   wxString strCtrl = m_pTextCtrlOutputScale->GetValue();
   double dValue;
@@ -964,13 +965,13 @@ DialogGetFilterParameters::getOutputScale (void)
 }
 
 const char*
-DialogGetFilterParameters::getFilterName (void)
+DialogGetFilterParameters::getFilterName ()
 {
   return m_pListBoxFilter->getSelectionStringValue();
 }
 
 const char*
-DialogGetFilterParameters::getDomainName (void)
+DialogGetFilterParameters::getDomainName ()
 {
   return m_pListBoxDomain->getSelectionStringValue();
 }
@@ -1011,7 +1012,7 @@ DialogExportParameters::DialogExportParameters (wxFrame* pParent, int iDefaultFo
 }
 
 const char*
-DialogExportParameters::getFormatName(void)
+DialogExportParameters::getFormatName()
 {
   return m_pListBoxFormat->getSelectionStringValue();
 }
@@ -1063,12 +1064,12 @@ DialogGetXYSize::DialogGetXYSize (wxFrame* pParent, const char* const pszTitle,
   pTopSizer->SetSizeHints (this);
 }
 
-DialogGetXYSize::~DialogGetXYSize (void)
+DialogGetXYSize::~DialogGetXYSize ()
 {
 }
 
 unsigned int
-DialogGetXYSize::getXSize (void)
+DialogGetXYSize::getXSize ()
 {
   wxString strCtrl = m_pTextCtrlXSize->GetValue();
   long lValue;
@@ -1079,7 +1080,7 @@ DialogGetXYSize::getXSize (void)
 }
 
 unsigned int
-DialogGetXYSize::getYSize (void)
+DialogGetXYSize::getYSize ()
 {
   wxString strCtrl = m_pTextCtrlYSize->GetValue();
   long lValue;
@@ -1089,3 +1090,115 @@ DialogGetXYSize::getYSize (void)
     return (m_iDefaultYSize);
 }
 
+
+
+/////////////////////////////////////////////////////////////////////
+// CLASS IDENTIFICATION
+//
+// DialogGetConvertPolarParameters
+/////////////////////////////////////////////////////////////////////
+
+DialogGetConvertPolarParameters::DialogGetConvertPolarParameters (wxFrame* pParent, const char* const pszTitle, 
+       int iDefaultXSize, int iDefaultYSize, int iDefaultInterpolationID, int iDefaultZeropad)
+: wxDialog (pParent, -1, pszTitle, wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)
+{
+  m_iDefaultXSize = iDefaultXSize;
+  m_iDefaultYSize = iDefaultYSize;
+  m_iDefaultZeropad = iDefaultZeropad;
+
+  wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL);
+  
+  pTopSizer->Add (new wxStaticText (this, -1, pszTitle), 0, wxALIGN_CENTER | wxTOP | wxLEFT | wxRIGHT, 5);
+  
+  pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);
+  
+  std::ostringstream os;
+  os << iDefaultXSize;
+  m_pTextCtrlXSize = new wxTextCtrl (this, -1, os.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+  std::ostringstream osYSize;
+  osYSize << iDefaultYSize;
+  m_pTextCtrlYSize = new wxTextCtrl (this, -1, osYSize.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+  
+  wxFlexGridSizer* pGridSizer = new wxFlexGridSizer (2);
+  
+  m_pListBoxInterpolation = new StringValueAndTitleListBox (this, Projections::getInterpCount(), Projections::getInterpTitleArray(), Projections::getInterpNameArray());
+  m_pListBoxInterpolation->SetSelection (iDefaultInterpolationID);
+  pGridSizer->Add (new wxStaticText (this, -1, "Interpolation"), 0, wxALIGN_CENTER_VERTICAL | wxALIGN_RIGHT | wxALL, 5);
+  pGridSizer->Add (m_pListBoxInterpolation, 0, wxALL | wxALIGN_LEFT | wxEXPAND);
+  
+  pGridSizer->Add (new wxStaticText (this, -1, "X Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (m_pTextCtrlXSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (new wxStaticText (this, -1, "Y Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+  pGridSizer->Add (m_pTextCtrlYSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  if (iDefaultZeropad >= 0) {
+    std::ostringstream osZeropad;
+    osZeropad << iDefaultZeropad;
+    m_pTextCtrlZeropad = new wxTextCtrl (this, -1, osZeropad.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+    pGridSizer->Add (new wxStaticText (this, -1, "Zeropad"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+    pGridSizer->Add (m_pTextCtrlZeropad, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL);
+  }
+  
+  pTopSizer->Add (pGridSizer, 1, wxALL, 3);
+  
+  pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);
+  
+  wxBoxSizer* pButtonSizer = new wxBoxSizer (wxHORIZONTAL);
+  wxButton* pButtonOk = new wxButton (this, wxID_OK, "Okay");
+  wxButton* pButtonCancel = new wxButton (this, wxID_CANCEL, "Cancel");
+  pButtonSizer->Add (pButtonOk, 0, wxEXPAND | wxALL, 10);
+  pButtonSizer->Add (pButtonCancel, 0, wxEXPAND | wxALL, 10);
+  
+  pTopSizer->Add (pButtonSizer, 0, wxALIGN_CENTER);
+  
+  SetAutoLayout (true);
+  SetSizer (pTopSizer);
+  pTopSizer->Layout();
+  pTopSizer->Fit (this);
+  pTopSizer->SetSizeHints (this);
+}
+
+
+DialogGetConvertPolarParameters::~DialogGetConvertPolarParameters ()
+{
+}
+
+
+unsigned int
+DialogGetConvertPolarParameters::getXSize ()
+{
+  wxString strCtrl = m_pTextCtrlXSize->GetValue();
+  unsigned long lValue;
+  if (strCtrl.ToULong (&lValue))
+    return lValue;
+  else
+    return (m_iDefaultXSize);
+}
+
+unsigned int
+DialogGetConvertPolarParameters::getYSize ()
+{
+  wxString strCtrl = m_pTextCtrlYSize->GetValue();
+  unsigned long lValue;
+  if (strCtrl.ToULong (&lValue))
+    return lValue;
+  else
+    return (m_iDefaultYSize);
+}
+
+unsigned int
+DialogGetConvertPolarParameters::getZeropad ()
+{
+  wxString strCtrl = m_pTextCtrlZeropad->GetValue();
+  unsigned long lValue;
+  if (strCtrl.ToULong (&lValue))
+    return lValue;
+  else
+    return (m_iDefaultZeropad);
+}
+
+const char*
+DialogGetConvertPolarParameters::getInterpolationName ()
+{
+  return m_pListBoxInterpolation->getSelectionStringValue();
+}
+
index 0553b1a..d16981c 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.h,v 1.19 2001/01/03 22:01:50 kevin Exp $
+**  $Id: dialogs.h,v 1.20 2001/01/04 21:28:41 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
@@ -36,6 +36,7 @@
 #include "phantom.h"
 #include "procsignal.h"
 #include "filter.h"
+#include "projections.h"
 
 
 // CLASS StringValueAndTitleListBox
@@ -291,5 +292,31 @@ class DialogGetXYSize : public wxDialog
 };
 
 
+class DialogGetConvertPolarParameters : public wxDialog
+{
+ public:
+   DialogGetConvertPolarParameters (wxFrame* pParent, const char* const pszTitle, int iDefaultXSize = 0, 
+     int iDefaultYSize = 0, int iDefaultInterpolationID = Projections::POLAR_INTERP_BILINEAR, 
+     int iDefaultZeropad = 1);
+   virtual ~DialogGetConvertPolarParameters ();
+
+    unsigned int getXSize();
+    unsigned int getYSize();
+    const char* getInterpolationName();
+    unsigned int getZeropad();
+
+ private:
+    wxTextCtrl* m_pTextCtrlXSize;
+    wxTextCtrl* m_pTextCtrlYSize;
+    wxTextCtrl* m_pTextCtrlZeropad;
+
+    StringValueAndTitleListBox* m_pListBoxInterpolation;
+
+    int m_iDefaultXSize;
+    int m_iDefaultYSize;
+    int m_iDefaultZeropad;
+};
+
+
 #endif
 
index fb2f8cd..1d8c5bf 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dlgprojections.cpp,v 1.15 2001/01/02 16:02:13 kevin Exp $
+**  $Id: dlgprojections.cpp,v 1.16 2001/01/04 21:28:41 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
@@ -68,7 +68,8 @@ IMPLEMENT_CLASS(ProjectionsDialog, wxDialog)
 
 
 ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, const Phantom& rPhantom, const int iTrace, wxWindow *parent)
-: wxDialog(parent, -1, "Collect Projections", wxDefaultPosition), m_rScanner(rScanner), m_rProjections(rProj), m_rPhantom(rPhantom), m_pSGPDriver(NULL), m_pSGP(NULL), m_iTrace(iTrace), m_pDC(NULL), m_btnAbort(0), m_btnPause(0), m_btnStep(0)
+: wxDialog(parent, -1, "Collect Projections", wxDefaultPosition), m_rScanner(rScanner), m_rProjections(rProj), m_rPhantom(rPhantom), 
+  m_pSGPDriver(NULL), m_pSGP(NULL), m_iTrace(iTrace), m_pDC(NULL), m_btnAbort(0), m_btnPause(0), m_btnStep(0)
 {
     m_state = Continue;
     m_iLastView = -1;
@@ -135,7 +136,7 @@ ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, con
        
     wxYield();     // Update the display
        
-    m_pSGPDriver->idWX()->SetFont(*wxSWISS_FONT);
+    m_pSGP->setTextPointSize(10);
 #ifdef __WXMAC__
     MacUpdateImmediately();
 #endif
@@ -228,7 +229,6 @@ ProjectionsDialog::OnPause (wxCommandEvent& event)
        } else if (m_state == Continue) {
                m_memoryDC.SelectObject (m_bitmap);       // in memoryDC
                m_pSGP->setDC (&m_memoryDC);
-               m_memoryDC.SetFont (*wxSWISS_FONT);
                showView (m_iLastView);
                m_state = Paused;
                m_btnPause->SetLabel (wxString("Resume"));
@@ -248,7 +248,6 @@ ProjectionsDialog::OnStep (wxCommandEvent& event)
        } else if (m_state == Continue) {
                m_memoryDC.SelectObject (m_bitmap);       // in memoryDC
                m_pSGP->setDC (&m_memoryDC);
-               m_memoryDC.SetFont (*wxSWISS_FONT);
                showView (m_iLastView);
                // m_rScanner.collectProjections (m_rProjections, m_rPhantom, m_iLastView, 1, true, m_iTrace, m_pSGP);
                m_state = Paused;
@@ -259,7 +258,6 @@ ProjectionsDialog::OnStep (wxCommandEvent& event)
        } else if (m_state == Paused) {
                m_memoryDC.SelectObject (m_bitmap);       // in memoryDC
                m_pSGP->setDC (&m_memoryDC);
-               m_memoryDC.SetFont (*wxSWISS_FONT);
                projectView (m_iLastView + 1);
                m_pSGP->setDC (m_pDC);
                m_memoryDC.SelectObject(wxNullBitmap);
index f7d8a1f..aba89cc 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.50 2001/01/03 22:00:46 kevin Exp $
+**  $Id: views.cpp,v 1.51 2001/01/04 21:28:41 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
@@ -170,6 +170,7 @@ EVT_MENU(IFMENU_FILE_EXPORT, ImageFileView::OnExport)
 EVT_MENU(IFMENU_FILE_PROPERTIES, ImageFileView::OnProperties)
 EVT_MENU(IFMENU_VIEW_SCALE_MINMAX, ImageFileView::OnScaleMinMax)
 EVT_MENU(IFMENU_VIEW_SCALE_AUTO, ImageFileView::OnScaleAuto)
+EVT_MENU(IFMENU_VIEW_SCALE_FULL, ImageFileView::OnScaleFull)
 EVT_MENU(IFMENU_COMPARE_IMAGES, ImageFileView::OnCompare)
 EVT_MENU(IFMENU_COMPARE_ROW, ImageFileView::OnCompareRow)
 EVT_MENU(IFMENU_COMPARE_COL, ImageFileView::OnCompareCol)
@@ -286,6 +287,16 @@ ImageFileView::OnScaleMinMax (wxCommandEvent& event)
   }
 }
 
+void 
+ImageFileView::OnScaleFull (wxCommandEvent& event)
+{
+  if (m_bMinSpecified || m_bMaxSpecified) {
+    m_bMinSpecified = false;
+    m_bMaxSpecified = false;
+    OnUpdate (this, NULL);
+  }
+}
+
 void
 ImageFileView::OnCompare (wxCommandEvent& event)
 {
@@ -304,7 +315,7 @@ ImageFileView::OnCompare (wxCommandEvent& event)
       std::ostringstream os;
       double min, max, mean, mode, median, stddev;
       rIF.statistics (min, max, mean, mode, median, stddev);
-      os << rIF.getFilename() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";
+      os << GetFrame()->GetTitle().c_str() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";
       rCompareIF.statistics (min, max, mean, mode, median, stddev);
       os << pCompareDoc->GetFirstView()->GetFrame()->GetTitle().c_str() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n";
       os << "\n";
@@ -326,7 +337,7 @@ ImageFileView::OnCompare (wxCommandEvent& event)
           return;
         }
         
-        wxString s = GetDocument()->GetFirstView()->GetFrame()->GetTitle() + ": ";
+        wxString s = GetFrame()->GetTitle() + ": ";
         differenceImage.labelsCopy (rIF, s.c_str());
         s = pCompareDoc->GetFirstView()->GetFrame()->GetTitle() + ": ";
         differenceImage.labelsCopy (rCompareIF, s.c_str());
@@ -561,7 +572,7 @@ ImageFileView::OnDivide (wxCommandEvent& event)
 }
 
 
-#ifdef HAVE_FFTW
+#ifdef HAVE_FFT
 void
 ImageFileView::OnFFT (wxCommandEvent& event)
 {
@@ -604,6 +615,7 @@ ImageFileView::OnFourier (wxCommandEvent& event)
     GetDocument()->Modify(TRUE);
   GetDocument()->UpdateAllViews(this);
 }
+
 void
 ImageFileView::OnInverseFourier (wxCommandEvent& event)
 {
@@ -717,6 +729,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   wxMenu *view_menu = new wxMenu;
   view_menu->Append(IFMENU_VIEW_SCALE_MINMAX, "Display Scale &Set...");
   view_menu->Append(IFMENU_VIEW_SCALE_AUTO, "Display Scale &Auto...");
+  view_menu->Append(IFMENU_VIEW_SCALE_FULL, "Display &Full Scale");
   
   wxMenu* filter_menu = new wxMenu;
   filter_menu->Append (IFMENU_FILTER_INVERTVALUES, "&Invert Values");
@@ -725,7 +738,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   filter_menu->Append (IFMENU_FILTER_LOG, "&Log");
   filter_menu->Append (IFMENU_FILTER_EXP, "&Exp");
   filter_menu->AppendSeparator();
-#ifdef HAVE_FFTW
+#ifdef HAVE_FFT
   filter_menu->Append (IFMENU_FILTER_FFT, "&FFT");
   filter_menu->Append (IFMENU_FILTER_IFFT, "&IFFT");
   filter_menu->Append (IFMENU_FILTER_FOURIER, "F&ourier");
@@ -1283,7 +1296,7 @@ PhantomView::PhantomView(void)
   m_iDefaultNDet = 367;
   m_iDefaultNView = 320;
   m_iDefaultNSample = 2;
-  m_dDefaultRotation = 2;
+  m_dDefaultRotation = 1;
   m_dDefaultFocalLength = 2;
   m_dDefaultFieldOfView = 1;
   m_iDefaultGeometry = Scanner::GEOMETRY_PARALLEL;
@@ -1617,6 +1630,11 @@ ProjectionFileView::ProjectionFileView(void)
   m_iDefaultInterpolation = Backprojector::INTERP_LINEAR;
   m_iDefaultInterpParam = 1;
   m_iDefaultTrace = Trace::TRACE_NONE;
+
+  m_iDefaultPolarNX = 256;
+  m_iDefaultPolarNY = 256;
+  m_iDefaultPolarInterpolation = Projections::POLAR_INTERP_BILINEAR;
+  m_iDefaultPolarZeropad = 1;
 }
 
 ProjectionFileView::~ProjectionFileView(void)
@@ -1639,23 +1657,67 @@ void
 ProjectionFileView::OnConvertPolar (wxCommandEvent& event)
 {
   Projections& rProj = GetDocument()->getProjections();
+  DialogGetConvertPolarParameters dialogPolar (m_frame, "Convert Polar", m_iDefaultPolarNX, m_iDefaultPolarNY,
+    m_iDefaultPolarInterpolation, -1);
+  if (dialogPolar.ShowModal() == wxID_OK) {
+    wxString strInterpolation (dialogPolar.getInterpolationName());
+    m_iDefaultPolarNX = dialogPolar.getXSize();
+    m_iDefaultPolarNY = dialogPolar.getYSize();
+    ImageFileDocument* pPolarDoc = dynamic_cast<ImageFileDocument*>(theApp->getDocManager()->CreateDocument("untitled.if", wxDOC_SILENT));
+    ImageFile& rIF = pPolarDoc->getImageFile();
+    if (! pPolarDoc) {
+      sys_error (ERR_SEVERE, "Unable to create image file");
+      return;
+    }
+    rIF.setArraySize (m_iDefaultPolarNX, m_iDefaultPolarNY);
+    m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.c_str());
+    rProj.convertPolar (rIF, m_iDefaultPolarInterpolation);
+    rIF.labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime());
+    std::ostringstream os;
+    os << "Convert projection file " << GetFrame()->GetTitle().c_str() << " to polar image: xSize=" 
+      << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation=" 
+      << strInterpolation.c_str();
+    *theApp->getLog() << os.str().c_str() << "\n";
+    rIF.labelAdd (os.str().c_str());
+    if (theApp->getSetModifyNewDocs())
+      pPolarDoc->Modify(true);
+    pPolarDoc->UpdateAllViews();
+    pPolarDoc->GetFirstView()->OnUpdate (this, NULL);
+  }
 }
 
 void
 ProjectionFileView::OnConvertFFTPolar (wxCommandEvent& event)
 {
   Projections& rProj = GetDocument()->getProjections();
-  wxMessageBox ("Polar conversion not yet implemented", "Unimplemented function");
-#if 0
-  rProj.convertPolar ();
-  if (theApp->getSetModifyNewDocs())
-    GetDocument()->Modify(true);
-  GetDocument()->UpdateAllViews();
-#ifndef HAVE_FFT
-  wxMessageBox ("FFT support has not been compiled into this version of CTSim", "Error");
-#endif
-#endif
-}
+  DialogGetConvertPolarParameters dialogPolar (m_frame, "Convert to FFT Polar", m_iDefaultPolarNX, m_iDefaultPolarNY,
+    m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad);
+  if (dialogPolar.ShowModal() == wxID_OK) {
+    wxString strInterpolation (dialogPolar.getInterpolationName());
+    m_iDefaultPolarNX = dialogPolar.getXSize();
+    m_iDefaultPolarNY = dialogPolar.getYSize();
+    m_iDefaultPolarZeropad = dialogPolar.getZeropad();
+    ImageFileDocument* pPolarDoc = dynamic_cast<ImageFileDocument*>(theApp->getDocManager()->CreateDocument("untitled.if", wxDOC_SILENT));
+    ImageFile& rIF = pPolarDoc->getImageFile();
+    if (! pPolarDoc) {
+      sys_error (ERR_SEVERE, "Unable to create image file");
+      return;
+    }
+    rIF.setArraySize (m_iDefaultPolarNX, m_iDefaultPolarNY);
+    m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.c_str());
+    rProj.convertFFTPolar (rIF, m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad);
+    rIF.labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime());
+    std::ostringstream os;
+    os << "Convert projection file " << GetFrame()->GetTitle().c_str() << " to FFT polar image: xSize=" 
+      << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation=" 
+      << strInterpolation.c_str() << ", zeropad=" << m_iDefaultPolarZeropad;
+    *theApp->getLog() << os.str().c_str() << "\n";
+    rIF.labelAdd (os.str().c_str());
+    if (theApp->getSetModifyNewDocs())
+      pPolarDoc->Modify(true);
+    pPolarDoc->UpdateAllViews();
+    pPolarDoc->GetFirstView()->OnUpdate (this, NULL);
+  }}
 
 void
 ProjectionFileView::OnReconstructFourier (wxCommandEvent& event)
@@ -1947,6 +2009,7 @@ BEGIN_EVENT_TABLE(PlotFileView, wxView)
 EVT_MENU(PJMENU_FILE_PROPERTIES, PlotFileView::OnProperties)
 EVT_MENU(PLOTMENU_VIEW_SCALE_MINMAX, PlotFileView::OnScaleMinMax)
 EVT_MENU(PLOTMENU_VIEW_SCALE_AUTO, PlotFileView::OnScaleAuto)
+EVT_MENU(PLOTMENU_VIEW_SCALE_FULL, PlotFileView::OnScaleFull)
 END_EVENT_TABLE()
 
 PlotFileView::PlotFileView(void) 
@@ -2026,6 +2089,16 @@ PlotFileView::OnScaleMinMax (wxCommandEvent& event)
   }
 }
 
+void 
+PlotFileView::OnScaleFull (wxCommandEvent& event)
+{
+  if (m_bMinSpecified || m_bMaxSpecified) {
+    m_bMinSpecified = false;
+    m_bMaxSpecified = false;
+    OnUpdate (this, NULL);
+  }
+}
+
 
 PlotFileCanvas* 
 PlotFileView::CreateCanvas (wxView *view, wxFrame *parent)
@@ -2067,6 +2140,7 @@ PlotFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   wxMenu *view_menu = new wxMenu;
   view_menu->Append(PLOTMENU_VIEW_SCALE_MINMAX, "Display Scale &Set...");
   view_menu->Append(PLOTMENU_VIEW_SCALE_AUTO, "Display Scale &Auto...");
+  view_menu->Append(PLOTMENU_VIEW_SCALE_FULL, "Display &Full Scale");
   
   wxMenu *help_menu = new wxMenu;
   help_menu->Append(MAINMENU_HELP_CONTENTS, "&Contents");
index 4463689..7ca2e62 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.h,v 1.24 2001/01/03 22:00:46 kevin Exp $
+**  $Id: views.h,v 1.25 2001/01/04 21:28:41 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
@@ -87,7 +87,7 @@ public:
   void OnShuffleNaturalToFourierOrder (wxCommandEvent& event);
   void OnShuffleFourierToNaturalOrder (wxCommandEvent& event);
   
-#ifdef HAVE_FFTW
+#ifdef HAVE_FFT
   void OnFFT (wxCommandEvent& event);
   void OnIFFT (wxCommandEvent& event);
 #endif
@@ -97,6 +97,7 @@ public:
   
   void OnScaleAuto (wxCommandEvent& event);
   void OnScaleMinMax (wxCommandEvent& event);
+  void OnScaleFull (wxCommandEvent& event);
   void OnPlotRow (wxCommandEvent& event);
   void OnPlotCol (wxCommandEvent& event);
   void OnPlotHistogram (wxCommandEvent& event);
@@ -159,6 +160,11 @@ private:
   int m_iDefaultBackprojector;
   int m_iDefaultTrace;
   
+  int m_iDefaultPolarNX;
+  int m_iDefaultPolarNY;
+  int m_iDefaultPolarInterpolation;
+  int m_iDefaultPolarZeropad;
+
 public:
   ProjectionFileView(void);
   virtual ~ProjectionFileView(void);
@@ -266,9 +272,11 @@ public:
   void OnDraw(wxDC* dc);
   void OnUpdate(wxView *sender, wxObject *hint = NULL);
   bool OnClose (bool deleteWindow = true);
+
   void OnProperties (wxCommandEvent& event);
-  void OnScaleAuto (wxCommandEvent& event);
   void OnScaleMinMax (wxCommandEvent& event);
+  void OnScaleAuto (wxCommandEvent& event);
+  void OnScaleFull (wxCommandEvent& event);
   
   wxFrame* getFrame ()
   { return m_frame; }