r127: *** empty log message ***

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 3a071b478bceb2f25761068260791512d3f32ca7..6e21b3568285ded1049dccf2d2ca9f7052e26a8f 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.7 2000/06/29 12:39:46 kevin Exp $
+**  $Id: projections.cpp,v 1.8 2000/07/02 18:21:39 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
@@ -490,10 +490,10 @@ bool
 Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, 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
+  double detInc = m_detInc;
+  double detlen = (m_nDet - 1) * detInc;
+  int n_filteredProj = m_nDet;
+  double filteredProj [n_filteredProj];   // filtered projections
 
 #ifdef HAVE_BSPLINE_INTERP
   int spline_order = 0, zoom_factor = 0;
@@ -501,16 +501,13 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     zoom_factor = interp_param;
     spline_order = 3;
     zoom_factor = 3;
-    n_filtered_proj = (m_nDet - 1) * (zoom_factor + 1) + 1;
+    n_filteredProj = (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;
+  double filterBW = 1. / detInc;
 
-  SignalFilter filter (filterName, filterBW, filterMin, filterMax, n_vec_filter, filt_param, "spatial", 0);
+  SignalFilter filter (filterName, filterMethodName, filterBW, detlen, m_nDet, filt_param, "spatial", 0);
   if (filter.fail()) {
       sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", filter.failMessage().c_str());
       return false;
@@ -521,13 +518,14 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 
 #if HAVE_SGP
   SGP_ID gid;
+  int n_vec_filter = filter.getNFilterPoints();
   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)
+    double filterInc = (detlen * 2) / (n_vec_filter - 1);
+    for (i = 0, f = -detlen; i < n_vec_filter; i++, f += filterInc)
       plot_xaxis[i] = f;
       
     gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
@@ -536,7 +534,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     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);
+    printf ("nview=%d, ndet=%d, det_start=%.4f, detInc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
   }
 #endif  //HAVE_SGP
 
@@ -553,8 +551,10 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     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);
+    filter.filterSignal (detval, filteredProj, detInc, m_nDet);
+
+    //    for (int j = 0; j < m_nDet; j++)
+    //      filteredProj[j] = filter.convolve (detval, detInc, j, m_nDet);
 
 #ifdef HAVE_SGP
     if (trace >= TRACE_PLOT)  {
@@ -576,23 +576,23 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
       ezset  ("xlength .5.");
       ezset ("box");
       ezset ("grid");
-      gid = ezplot (filtered_proj, plot_xaxis, m_nDet);
+      gid = ezplot (filteredProj, 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);
+       bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
     
 #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);
+       bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
+      ezplot_1d (filteredProj, n_filteredProj);
     }
 #endif
 #endif
 
-    bj.BackprojectView (filtered_proj, darray.viewAngle());
+    bj.BackprojectView (filteredProj, darray.viewAngle());
 
 #ifdef HAVE_SGP
     if (trace >= TRACE_PLOT) {