r636: Optimized Rebinning, Added Reconstruct with Rebinning option

[ctsim.git] / libctsim / reconstruct.cpp

diff --git a/libctsim/reconstruct.cpp b/libctsim/reconstruct.cpp
index df456214c7b2a4066e7d9caa986c5eff9591d352..8caed206c2df46997d38eb77b3c3ef674b644db3 100644 (file)
--- a/libctsim/reconstruct.cpp
+++ b/libctsim/reconstruct.cpp
@@ -6,9 +6,9 @@
 **   Date Started: Aug 84
 **
 **  This is part of the CTSim program
-**  Copyright (C) 1983-2000 Kevin Rosenberg
+**  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: reconstruct.cpp,v 1.6 2000/12/29 15:45:06 kevin Exp $
+**  $Id: reconstruct.cpp,v 1.17 2001/03/13 08:24: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
@@ -49,10 +49,17 @@
  */
 
 
-Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int iTrace, SGP* pSGP)
-  : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace), m_bFail(false), m_adPlotXAxis(0)
+Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const char* const filterName, 
+                              double filt_param, const char* const filterMethodName, const int zeropad, 
+                              const char* filterGenerationName, const char* const interpName, 
+                              int interpFactor, const char* const backprojectName, const int iTrace, 
+                              ReconstructionROI* pROI, bool bRebinToParallel, SGP* pSGP)
+  : m_rOriginalProj(rProj), 
+    m_pProj(bRebinToParallel ? m_rOriginalProj.interpolateToParallel() : &m_rOriginalProj),
+    m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), 
+    m_iTrace(iTrace), m_bRebinToParallel(bRebinToParallel), m_bFail(false), m_adPlotXAxis(0)
 {
-  m_nFilteredProjections = m_rProj.nDet() * interpFactor;
+  m_nFilteredProjections = m_pProj->nDet() * interpFactor;
 
 #ifdef HAVE_BSPLINE_INTERP
   int spline_order = 0, zoom_factor = 0;
@@ -64,8 +71,10 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch
   }
 #endif
 
-  double filterBW = 1. / m_rProj.detInc();
-  m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_rProj.detInc(), m_rProj.nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace, m_rProj.geometry(), m_rProj.focalLength(), pSGP);
+  double filterBW = 1. / m_pProj->detInc();
+  m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_pProj->detInc(), 
+    m_pProj->nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace, 
+    m_pProj->geometry(), m_pProj->focalLength(), m_pProj->sourceDetectorLength(), pSGP);
 
   if (m_pProcessSignal->fail()) {
     m_bFail = true;
@@ -75,7 +84,7 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch
     return;
   }
 
-  m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor);
+  m_pBackprojector = new Backprojector (*m_pProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI);
   if (m_pBackprojector->fail()) {
     m_bFail = true;
     m_strFailMessage = "Error creating backprojector: ";
@@ -85,18 +94,21 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch
     return;
   }
 
-#if HAVE_SGP
-  m_adPlotXAxis = new double [m_rProj.nDet()];
-  double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc();
-  double xInc = m_rProj.detInc();
+#ifdef HAVE_SGP
+  m_adPlotXAxis = new double [m_pProj->nDet()];
+  double x = - ((m_pProj->nDet() - 1) / 2) * m_pProj->detInc();
+  double xInc = m_pProj->detInc();
 
-  for (int i = 0; i < m_rProj.nDet(); i++, x += xInc)
+  for (int i = 0; i < m_pProj->nDet(); i++, x += xInc)
     m_adPlotXAxis[i] = x;
 #endif
 }
 
 Reconstructor::~Reconstructor ()
 {
+  if (m_bRebinToParallel)
+    delete m_pProj;
+
   delete m_pBackprojector;
   delete m_pProcessSignal;
   delete m_adPlotXAxis;
@@ -106,7 +118,7 @@ Reconstructor::~Reconstructor ()
 void
 Reconstructor::plotFilter (SGP* pSGP)
 {
-#if HAVE_SGP
+#ifdef HAVE_SGP
   int nVecFilter = m_pProcessSignal->getNFilterPoints();
   double* adPlotXAxis = new double [nVecFilter];
 
@@ -123,7 +135,7 @@ Reconstructor::plotFilter (SGP* pSGP)
       ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter);
       ezplot.plot (pSGP);
     }
-  }\r
+  }
   delete adPlotXAxis;
 #endif
 }
@@ -132,40 +144,46 @@ Reconstructor::plotFilter (SGP* pSGP)
 void
 Reconstructor::reconstructAllViews ()
 {
-  reconstructView (0, m_rProj.nView());
-  delete m_pBackprojector; m_pBackprojector = NULL;
+  reconstructView (0, m_pProj->nView());
+  postProcessing();
+}
+
+void
+Reconstructor::postProcessing()
+{
+  m_pBackprojector->PostProcessing();
 }
 
 
 void
-Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView)
+Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth)
 {
   double* adFilteredProj = new double [m_nFilteredProjections];   // filtered projections
 
   if (iViewCount <= 0)
-    iViewCount = m_rProj.nView() - iStartView;
+    iViewCount = m_pProj->nView() - iStartView;
       
   for (int iView = iStartView; iView < (iStartView + iViewCount); iView++)  {
     if (m_iTrace == Trace::TRACE_CONSOLE) 
-               std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
+               std::cout <<"Reconstructing view " << iView << " (last = " << m_pProj->nView() - 1 << ")\n";
       
-    const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
+    const DetectorArray& rDetArray = m_pProj->getDetectorArray (iView);
     const DetectorValue* detval = rDetArray.detValues();
 
     m_pProcessSignal->filterSignal (detval, adFilteredProj);
 
 #ifdef HAVE_BSPLINE_INTERP
     if (interp_type == I_BSPLINE) 
-       bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
+       bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
     
 #ifdef HAVE_SGP
     if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
-       bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
+       bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
       ezplot_1d (adFilteredProj, m_nFilteredProjections);
     }
 #endif
 #endif
-\r
+
        if (bBackprojectView)
       m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
 
@@ -173,31 +191,52 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool
     if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
       EZPlot ezplotProj;
 
+      std::ostringstream osXLength;
+      osXLength << "xlength " << dGraphWidth;
+
       ezplotProj.ezset ("clear");
       ezplotProj.ezset ("title Raw Projection");
       ezplotProj.ezset ("xticks major 5");
+      ezplotProj.ezset ("yticks major 5");
       ezplotProj.ezset ("xlabel ");
       ezplotProj.ezset ("ylabel ");
-      ezplotProj.ezset ("yporigin .5");
-      ezplotProj.ezset ("ylength .5");
+      ezplotProj.ezset ("yporigin 0.55");
+      ezplotProj.ezset ("ylength 0.45");
+      ezplotProj.ezset (osXLength.str().c_str());
       ezplotProj.ezset ("box.");
       ezplotProj.ezset ("grid.");
-      ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
+#if 0  // workaround c++ optimizer bug, now disabled by using /O1 in code
+      double* pdDetval = new double [m_pProj->nDet()];
+      for (unsigned int id = 0; id < m_pProj->nDet(); id++) {
+        pdDetval[id] = detval[id];
+      }
+      ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_pProj->nDet());
+      delete pdDetval;
+#else
+      ezplotProj.addCurve (m_adPlotXAxis, detval, m_pProj->nDet());
+#endif
+      pSGP->setTextPointSize (12);
       ezplotProj.plot (pSGP);
+
       ezplotProj.ezset ("clear");
       ezplotProj.ezset ("title Filtered Projection");
       ezplotProj.ezset ("xticks major 5");
       ezplotProj.ezset ("xlabel ");
       ezplotProj.ezset ("ylabel ");
-      ezplotProj.ezset ("ylength .5");
+      ezplotProj.ezset ("yticks major 5");
+      ezplotProj.ezset ("yporigin 0.10");
+      ezplotProj.ezset ("ylength 0.45");
+      ezplotProj.ezset (osXLength.str().c_str());
       ezplotProj.ezset ("box");
       ezplotProj.ezset ("grid");
       ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj,  m_nFilteredProjections);
+      pSGP->setTextPointSize (12);
       ezplotProj.plot (pSGP);
-    } 
+
+} 
 #endif  //HAVE_SGP
   }
-\r
-  delete adFilteredProj;\r
+
+  delete adFilteredProj;
 }