X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=df7ccd2b03b380372eea0124e471aa8711066cd0;hp=b54186447141fafa9f21b672dd120820c32bb20d;hb=1e88cf0f7fa4f690ea9f110e8ed3f2b5338d0a10;hpb=53c732778ed19ac5231bb17c7e5bd4d2201d9456

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index b541864..df7ccd2 100644
--- 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.9 2000/07/03 11:02:06 kevin Exp $
+**  $Id: projections.cpp,v 1.22 2000/08/25 15:59:13 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
@@ -81,8 +81,7 @@ Projections::initFromScanner (const Scanner& scanner)
   m_rotInc = scanner.rotInc();
   m_detInc = scanner.detInc();
   m_rotStart = 0;
-  m_detStart =  -scanner.radius() + (scanner.detInc() / 2);
-  m_phmLen = scanner.phmLen();
+  m_detStart =  -(scanner.detLen() / 2) + (scanner.detInc() / 2);
 }
 
 void
@@ -208,7 +207,6 @@ Projections::headerRead (fnetorderstream& fs)
   if (! fs)
       return false;
 
-  off_t testPos;
   fs.readInt16 (_hsize);
   fs.readInt16 (_signature);
   fs.readInt32 (_nView);
@@ -366,7 +364,7 @@ Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int
   darray.setViewAngle (view_angle);
   //  darray.setNDet ( nDet);
   
-  for (int i = 0; i < nDet; i++) {
+  for (unsigned int i = 0; i < nDet; i++) {
       kfloat32 detval;
       fs.readFloat32 (detval);
       detval_ptr[i] = detval;
@@ -412,7 +410,7 @@ Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, co
   fs.writeFloat64 (view_angle);
   fs.writeInt32 (nDet);
 
-  for (int i = 0; i < nDet; i++) {
+  for (unsigned int i = 0; i < nDet; i++) {
     kfloat32 detval = detval_ptr[i];
     fs.writeFloat32 (detval);
   }
@@ -487,17 +485,16 @@ Projections::printScanInfo (void) const
  */
 
 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)
+Projections::reconstruct (ImageFile& im, 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 trace) const
 {
-  int nview = m_nView;
   double detInc = m_detInc;
-  int n_filteredProj = m_nDet;
+  int n_filteredProj = m_nDet * interpFactor;
   double filteredProj [n_filteredProj];   // filtered projections
 
 #ifdef HAVE_BSPLINE_INTERP
   int spline_order = 0, zoom_factor = 0;
   if (interp_type == I_BSPLINE) {
-    zoom_factor = interp_param;
+    zoom_factor = interpFactor;
     spline_order = 3;
     zoom_factor = 3;
     n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1;
@@ -505,9 +502,10 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 #endif
 
   double filterBW = 1. / detInc;
-  SignalFilter filter (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", 0);
-  if (filter.fail()) {
-      sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", filter.failMessage().c_str());
+  ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace);
+
+  if (processSignal.fail()) {
+      sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str());
       return false;
   }
 
@@ -515,28 +513,34 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
 
 #if HAVE_SGP
-  SGP_ID gid;
-  int n_vec_filter = filter.getNFilterPoints();
-  double plot_xaxis [n_vec_filter]; 			// array for plotting 
+  int nVecFilter = processSignal.getNFilterPoints();
+  double plot_xaxis [nVecFilter]; 			// array for plotting 
 
-  if (trace > TRACE_TEXT)  {
+  if (trace > TRACE_TEXT && nVecFilter > 0)  {
     int i;
     double f;
-    double filterInc = filter.getFilterIncrement();
-    for (i = 0, f = filter.getFilterMin(); i < n_vec_filter; i++, f += filterInc)
+    double filterInc = processSignal.getFilterIncrement();
+    for (i = 0, f = processSignal.getFilterMin(); i < nVecFilter; 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 (processSignal.getFilter()) {
+      SGPDriver sgpDriver ("Filter Function");
+      SGP sgp (sgpDriver);
+      EZPlot ezplot (sgp);
+
+      ezplot.ezset ("title Filter Response");
+      ezplot.addCurve (plot_xaxis, processSignal.getFilter(), nVecFilter);
+      ezplot.plot();
+      cio_put_str ("Press any key to continue");
+      cio_kb_getc ();
+    }
   }
   if (trace >= TRACE_TEXT) {
     printf ("nview=%d, ndet=%d, det_start=%.4f, detInc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
   }
 #endif  //HAVE_SGP
 
-  Backprojector bj (*this, im, backprojectName, interpName);
+  Backprojector bj (*this, im, backprojectName, interpName, interpFactor);
   if (bj.fail()) {
     sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", bj.failMessage().c_str());
     return false;
@@ -546,37 +550,10 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     if (trace >= TRACE_TEXT) 
       printf ("Reconstructing view %d (last = %d)\n",  iview, m_nView - 1);
       
-    DetectorArray& darray = getDetectorArray (iview);
-    DetectorValue* detval = darray.detValues();
+    const DetectorArray& darray = getDetectorArray (iview);
+    const DetectorValue* detval = darray.detValues();
 
-    filter.filterSignal (detval, filteredProj);
-
-    //    for (int j = 0; j < m_nDet; j++)
-    //      filteredProj[j] = filter.convolve (detval, detInc, 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 (filteredProj, plot_xaxis, m_nDet);
-    }
-#endif  //HAVE_SGP
+    processSignal.filterSignal (detval, filteredProj);
 
 #ifdef HAVE_BSPLINE_INTERP
     if (interp_type == I_BSPLINE) 
@@ -594,10 +571,34 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 
 #ifdef HAVE_SGP
     if (trace >= TRACE_PLOT) {
+      SGPDriver sgpDriverProj ("Projection");
+      SGP sgpProj (sgpDriverProj);
+      EZPlot ezplotProj (sgpProj);
+
+      ezplotProj.ezset  ("clear");
+      ezplotProj.ezset ("title Filtered Projection");
+      ezplotProj.ezset  ("xticks major 5.");
+      ezplotProj.ezset  ("xlabel ");
+      ezplotProj.ezset  ("ylabel ");
+      ezplotProj.ezset  ("yporigin .5.");
+      ezplotProj.ezset  ("ylength .5.");
+      ezplotProj.ezset  ("box.");
+      ezplotProj.ezset  ("grid.");
+      ezplotProj.addCurve (plot_xaxis, detval, m_nDet);
+      ezplotProj.plot();
+      ezplotProj.ezset  ("clear");
+      ezplotProj.ezset  ("xticks major 5.");
+      ezplotProj.ezset  ("xlabel ");
+      ezplotProj.ezset  ("ylabel ");
+      ezplotProj.ezset  ("ylength .5.");
+      ezplotProj.ezset ("box");
+      ezplotProj.ezset ("grid");
+      ezplotProj.addCurve (plot_xaxis, filteredProj,  n_filteredProj);
+      ezplotProj.plot();
+
+      cout << "Do you want to exit with current pic (y/n)? " << flush;
       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;
       }