r117: *** empty log message ***

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 54f8f7d66cce47b12657f9509ac39ad667718ab9..5ce0a8e022a3ef788ff2d5e811972a0ec7d41ee4 100644 (file)
--- a/libctsim/projections.cpp
+++ b/libctsim/projections.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.3 2000/06/20 17:54:51 kevin Exp $
+**  $Id: projections.cpp,v 1.4 2000/06/22 10:17:28 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
@@ -427,3 +427,145 @@ Projections::printScanInfo (void) const
 }
 
 
+
+/* NAME
+ *   Projections::reconstruct      Reconstruct Image from Projections
+ *
+ * SYNOPSIS
+ *   im = proj.reconstruct (im, filt_type, filt_param, interp_type)
+ *   IMAGE *im                 Output image
+ *   int filt_type             Type of convolution filter to use
+ *   double filt_param         Filter specific parameter
+ *                             Currently, used only with Hamming filters
+ *   int interp_type           Type of interpolation method to use
+ *
+ * ALGORITHM
+ *
+ *     Calculate one-dimensional filter in spatial domain
+ *     Allocate & clear (zero) the 2d output image array
+ *      For each projection view
+ *         Convolve raysum array with filter
+ *         Backproject raysums and add (summate) to image array
+ *      end
+ */
+
+bool
+Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const interpName, int interp_param, const char* const backprojectName, const int trace)
+{
+  int nview = m_nView;
+  double det_inc = m_detInc;
+  double detlen = (m_nDet - 1) * det_inc;
+  int n_filtered_proj = m_nDet;
+  double filtered_proj [n_filtered_proj];   // convolved result
+
+#ifdef HAVE_BSPLINE_INTERP
+  int spline_order = 0, zoom_factor = 0;
+  if (interp_type == I_BSPLINE) {
+    zoom_factor = interp_param;
+    spline_order = 3;
+    zoom_factor = 3;
+    n_filtered_proj = (m_nDet - 1) * (zoom_factor + 1) + 1;
+  }
+#endif
+
+  int n_vec_filter = 2 * m_nDet - 1;
+  double filterBW = 1. / det_inc;
+  double filterMin = -detlen;
+  double filterMax = detlen;
+
+  SignalFilter filter (filterName, filterBW, filterMin, filterMax, n_vec_filter, filt_param, "spatial", 0);
+  if (filter.fail())
+    return false;
+
+  if (trace)
+    cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
+
+#if HAVE_SGP
+  SGP_ID gid;
+  double plot_xaxis [n_vec_filter];                    // array for plotting 
+
+  if (trace > TRACE_TEXT)  {
+    int i;
+    double f;
+    double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
+    for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
+      plot_xaxis[i] = f;
+      
+    gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
+    cio_put_str ("Press any key to continue");
+    cio_kb_getc ();
+    sgp2_close (gid);
+  }
+  if (trace >= TRACE_TEXT) {
+    printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
+  }
+#endif  //HAVE_SGP
+
+  Backprojector bj (*this, im, backprojectName, interpName);
+  if (bj.fail())
+    return false;
+
+  for (int iview = 0; iview < m_nView; iview++)  {
+    if (trace >= TRACE_TEXT) 
+      printf ("Reconstructing view %d (last = %d)\n",  iview, m_nView - 1);
+      
+    DetectorArray& darray = getDetectorArray (iview);
+    DetectorValue* detval = darray.detValues();
+
+    for (int j = 0; j < m_nDet; j++)
+      filtered_proj[j] = filter.convolve (detval, det_inc, j, m_nDet);
+
+#ifdef HAVE_SGP
+    if (trace >= TRACE_PLOT)  {
+      ezset  ("clear.");
+      ezset  ("xticks major 5.");
+      ezset  ("xlabel ");
+      ezset  ("ylabel ");
+      ezset  ("xlength .5.");
+      ezset  ("box.");
+      ezset  ("grid.");
+      ezset  ("ufinish yes.");
+      ezplot (detval, plot_xaxis, m_nDet);
+      ezset  ("clear.");
+      ezset  ("xticks major 5.");
+      ezset  ("xlabel ");
+      ezset  ("ylabel ");
+      ezset  ("ustart yes.");
+      ezset  ("xporigin .5.");
+      ezset  ("xlength .5.");
+      ezset ("box");
+      ezset ("grid");
+      gid = ezplot (filtered_proj, plot_xaxis, m_nDet);
+    }
+#endif  //HAVE_SGP
+
+#ifdef HAVE_BSPLINE_INTERP
+    if (interp_type == I_BSPLINE) 
+       bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+    
+#ifdef HAVE_SGP
+    if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
+       bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+      ezplot_1d (filtered_proj, n_filtered_proj);
+    }
+#endif
+#endif
+
+    bj.BackprojectView (filtered_proj, darray.viewAngle());
+
+#ifdef HAVE_SGP
+    if (trace >= TRACE_PLOT) {
+      char str[256];
+      printf ("Do you want to exit with current pic (y/n) -- ");
+      fgets(str, sizeof(str), stdin);
+      sgp2_close (sgp2_get_active_win());
+      if (tolower(str[0]) == 'y') {
+       break;
+      }
+    } 
+#endif  //HAVE_SGP
+  }
+
+  return true;
+}
+