X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=be8e82941f9a47b30fe1946905fa9fba96bc7f0c;hp=8d3689567550a02ea1a4e670135193410552951a;hb=a05f3cb550877e94aa118cc04b361c0c8fdb3dc3;hpb=99dd1d6ed10db1f669a5fe6af71225a50fc0ddfb

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 8d36895..be8e829 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.1 2000/06/19 02:59:34 kevin Exp $
+**  $Id: projections.cpp,v 1.24 2000/08/31 08:38:58 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
@@ -39,24 +39,20 @@
  */
 
 Projections::Projections (const Scanner& scanner)
+  : m_projData(0)
 {
-  init (scanner.nView(), scanner.nDet());
-
-  m_phmLen = scanner.phmLen();
-  m_rotInc = scanner.rotInc();
-  m_detInc = scanner.detInc();
-  m_rotStart = 0;
-  m_detStart =  -scanner.radius() + (scanner.detInc() / 2);
-  m_phmLen = scanner.phmLen();
+  initFromScanner (scanner);
 }
 
 
 Projections::Projections (const int nView, const int nDet)
+  : m_projData(0)
 {
   init (nView, nDet);
 }
 
 Projections::Projections (void)
+  : m_projData(0)
 {
   init (0, 0);
 }
@@ -72,11 +68,37 @@ Projections::init (const int nView, const int nDet)
 {
   m_nView = nView;
   m_nDet = nDet;
-  m_projData = NULL;
   newProjData ();
-  m_fd = -1;
 }
 
+void
+Projections::initFromScanner (const Scanner& scanner)
+{
+  deleteProjData();
+  init (scanner.nView(), scanner.nDet());
+
+  m_phmLen = scanner.phmLen();
+  m_rotInc = scanner.rotInc();
+  m_detInc = scanner.detInc();
+  m_geometry = scanner.geometry();
+  m_focalLength = scanner.focalLength();
+  m_fieldOfView = scanner.fieldOfView();
+  m_rotStart = 0;
+  m_detStart =  -(scanner.detLen() / 2);
+#if 0
+  if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
+    m_detInc /= 2;
+    cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
+  }
+#endif
+}
+
+void
+Projections::setNView (int nView)  // used by MPI to reduce # of views
+{
+  deleteProjData();
+  init (nView, m_nDet);
+}
 
 // NAME
 // newProjData
@@ -118,120 +140,128 @@ Projections::deleteProjData (void)
 
 
 /* NAME
- *    projections_write_header         Write data header for projections file
+ *    Projections::headerWwrite         Write data header for projections file
  *
  */
 
 bool 
-Projections::headerWrite (void)
+Projections::headerWrite (fnetorderstream& fs)
 {
-  kuint32 _hsize = m_headerSize;
+  kuint16 _hsize = m_headerSize;
+  kuint16 _signature = m_signature;
   kuint32 _nView = m_nView;
   kuint32 _nDet = m_nDet;
   kuint32 _geom = m_geometry;
-  kuint32 _remarksize = m_remark.length();
+  kuint16 _remarksize = m_remark.length();
+  kuint16 _year = m_year;
+  kuint16 _month = m_month;
+  kuint16 _day = m_day;
+  kuint16 _hour = m_hour;
+  kuint16 _minute = m_minute;
+  kuint16 _second = m_second;
+
   kfloat64 _calcTime = m_calcTime;
   kfloat64 _rotStart = m_rotStart;
   kfloat64 _rotInc = m_rotInc;
   kfloat64 _detStart = m_detStart;
   kfloat64 _detInc = m_detInc;
   kfloat64 _phmLen = m_phmLen;
-  
-  if (m_fd < 0) {
-    sys_error (ERR_SEVERE, "m_fd < 0 [projections_write_header]");
-    return false;
-  }
-    
-  if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) {
-    sys_error (ERR_WARNING, "Error seeking to beginning of fiel [projections_write_header]");
-    return false;
-  }
+  kfloat64 _fieldOfView = m_fieldOfView;
+  kfloat64 _focalLength = m_focalLength;
 
-  if (! write_nint32 (&_hsize, m_fd) ||
-      ! write_nint32 (&_nView, m_fd) ||
-      ! write_nint32 (&_nDet, m_fd) ||
-      ! write_nint32 (&_geom, m_fd) ||
-      ! write_nfloat64 (&_calcTime, m_fd) ||
-      ! write_nfloat64 (&_rotStart, m_fd) ||
-      ! write_nfloat64 (&_rotInc, m_fd) ||
-      ! write_nfloat64 (&_detStart, m_fd) ||
-      ! write_nfloat64 (&_detInc, m_fd) ||
-      ! write_nfloat64 (&_phmLen, m_fd) ||
-      ! write_nint32 (&_remarksize, m_fd) ||
-      (::write (m_fd, m_remark.c_str(), _remarksize) != _remarksize)) {
-    sys_error (ERR_SEVERE, "Error writing header information [projections_write_header] %ld", _remarksize);
+  fs.seekp(0);
+  if (! fs)
     return false;
-  }
-  m_headerSize = _hsize = lseek(m_fd, (off_t) 0, SEEK_CUR);
-  if ((lseek(m_fd, 0, SEEK_SET) != 0) || ! write_nint32 (&_hsize, m_fd)) {
-    sys_error (ERR_SEVERE, "Error writing header information [projections_write_header]");
-    return false;
-  }
+
+  fs.writeInt16 (_hsize);
+  fs.writeInt16 (_signature);
+  fs.writeInt32 (_nView);
+  fs.writeInt32 (_nDet);
+  fs.writeInt32 (_geom);
+  fs.writeFloat64 (_calcTime);
+  fs.writeFloat64 (_rotStart);
+  fs.writeFloat64 (_rotInc);
+  fs.writeFloat64 (_detStart);
+  fs.writeFloat64 (_detInc);
+  fs.writeFloat64 (_phmLen);
+  fs.writeFloat64 (_focalLength);
+  fs.writeFloat64 (_fieldOfView);
+  fs.writeInt16 (_year);
+  fs.writeInt16 (_month);
+  fs.writeInt16 (_day);
+  fs.writeInt16 (_hour);
+  fs.writeInt16 (_minute);
+  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 true;
 }
 
-
 /* NAME
  *    projections_read_header         Read data header for projections file
  *
  */
-
 bool
-Projections::headerRead (void)
+Projections::headerRead (fnetorderstream& fs)
 {
-  kuint32 _hsize;
-  kuint32 _nView;
-  kuint32 _nDet;
-  kuint32 _geom;
-  kuint32 _remarksize = 0;
-  kfloat64 _calcTime;
-  kfloat64 _rotStart;
-  kfloat64 _rotInc;
-  kfloat64 _detStart;
-  kfloat64 _detInc;
-  kfloat64 _phmLen;
+  kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
+  kuint32 _nView, _nDet, _geom;
+  kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen, _focalLength, _fieldOfView;
   
-  if (m_fd < 0) {
-    sys_error (ERR_SEVERE, "m_fd < 0 [projections_read_header]");
-    return false;
-  }
-    
-  if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) {
-    sys_error (ERR_WARNING, "Error seeking to beginning of field [projections_read_header]");
-    return false;
+  fs.seekg(0);
+  if (! fs)
+      return false;
+
+  fs.readInt16 (_hsize);
+  fs.readInt16 (_signature);
+  fs.readInt32 (_nView);
+  fs.readInt32 (_nDet);
+  fs.readInt32 (_geom);
+  fs.readFloat64 (_calcTime);
+  fs.readFloat64 (_rotStart);
+  fs.readFloat64 (_rotInc);
+  fs.readFloat64 (_detStart);
+  fs.readFloat64 (_detInc);
+  fs.readFloat64 (_phmLen);
+  fs.readFloat64 (_focalLength);
+  fs.readFloat64 (_fieldOfView);
+  fs.readInt16 (_year);
+  fs.readInt16 (_month);
+  fs.readInt16 (_day);
+  fs.readInt16 (_hour);
+  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;
   }
 
-  off_t testPos;
-  if (! read_nint32 (&_hsize, m_fd)) {
-    sys_error (ERR_SEVERE, "Error reading header size , _hsize=%d [projections_read_header]", _hsize);
+  if (_signature != m_signature) {
+    sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
     return false;
   }
-  if ( ! read_nint32 (&_nView, m_fd) ||
-      ! read_nint32 (&_nDet, m_fd) ||
-      ! read_nint32 (&_geom, m_fd) ||
-      ! read_nfloat64 (&_calcTime, m_fd) ||
-      ! read_nfloat64 (&_rotStart, m_fd) ||
-      ! read_nfloat64 (&_rotInc, m_fd) ||
-      ! read_nfloat64 (&_detStart, m_fd) ||
-      ! read_nfloat64 (&_detInc, m_fd) ||
-      ! read_nfloat64 (&_phmLen, m_fd) ||
-       ! read_nint32 (&_remarksize, m_fd))
-      {
-	  sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
-	  return false;
-      }
 
   char remarkStorage[_remarksize+1];
-  if (::read (m_fd, &remarkStorage, _remarksize) != _remarksize) {
+  fs.read (remarkStorage, _remarksize);
+  if (! fs) {
     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
     return false;
   }
   remarkStorage[_remarksize] = 0;
   m_remark = remarkStorage;
 
-  off_t _hsizeread = lseek(m_fd, 0, SEEK_CUR);
-  if (_hsizeread != _hsize) {
+  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;
   }
@@ -246,66 +276,93 @@ Projections::headerRead (void)
   m_detStart = _detStart;
   m_detInc = _detInc;
   m_phmLen = _phmLen;
+  m_focalLength = _focalLength;
+  m_fieldOfView = _fieldOfView;
+  m_year = _year;
+  m_month = _month;
+  m_day = _day;
+  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;
 }
 
+bool
+Projections::read (const string& filename)
+{
+  return read (filename.c_str());
+}
+
 
 bool
 Projections::read (const char* filename)
 {
-  if ((m_fd = open(filename, OPEN_RDONLY | O_BINARY)) < 0) {
-    sys_error (ERR_SEVERE, "Error opening projections file %s for input [open_projections]", filename);
+  frnetorderstream fileRead (filename, ios::in | ios::binary);
+  m_filename = filename;
+
+  if (! fileRead)
     return false;
-  }
 
-  if (! headerRead()) {
-    sys_error (ERR_SEVERE, "Error reading projections file %s [open_projections]", filename);
+  if (! headerRead (fileRead))
     return false;
-  }
 
   deleteProjData ();
-  newProjData ();
+  newProjData();
 
   for (int i = 0; i < m_nView; i++) {
-    if (! detarrayRead (*m_projData[i], i))
+    if (! detarrayRead (fileRead, *m_projData[i], i))
       break;
   }
 
-  close (m_fd);
-  m_fd = -1;
-
+  fileRead.close();
   return true;
 }
 
 
+bool
+Projections::write (const string& filename)
+{
+  return write (filename.c_str());
+}
+
 bool
 Projections::write (const char* filename)
 {
-    if ((m_fd = open (filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY 
-#ifndef _WIN32
-			, S_IWRITE | S_IREAD | S_IRGRP | S_IROTH
-#endif
-			)) < 0) {
-      sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
-      return false;
-    }
-    if (! headerWrite ()) {
-      sys_error(ERR_SEVERE, "Error writing to file %s [create_projections]", filename);
-      return false;
-    }
+  frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
+  m_filename = filename;
+  if (! fs) {
+    sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
+    return false;
+  }
 
-    headerWrite ();
+  time_t t = time(NULL);
+  tm* lt = localtime(&t);
+  m_year = lt->tm_year;
+  m_month = lt->tm_mon;
+  m_day = lt->tm_mday;
+  m_hour = lt->tm_hour;
+  m_minute = lt->tm_min;
+  m_second = lt->tm_sec;
+
+  if (! headerWrite (fs))
+      return false;
 
   if (m_projData != NULL) {
     for (int i = 0; i < m_nView; i++) {
-      if (! detarrayWrite (*m_projData[i], i))
+      if (! detarrayWrite (fs, *m_projData[i], i))
 	break;
     }
   }
+  if (! fs)
+    return false;
 
-  close (m_fd);
-  m_fd = -1;
+ fs.close();
 
   return true;
 }
@@ -320,7 +377,7 @@ Projections::write (const char* filename)
  */
 
 bool
-Projections::detarrayRead (DetectorArray& darray, const int iview)
+Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
 {
   const int detval_bytes = darray.nDet() * sizeof(kfloat32);
   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
@@ -330,32 +387,20 @@ Projections::detarrayRead (DetectorArray& darray, const int iview)
   kfloat64 view_angle;
   kuint32 nDet;
 
-  if (m_fd < 0) { 
-    sys_error (ERR_SEVERE, "Tried to read detector array on closed file [detarrayRead]"); 
-    return false; 
-  } 
-
-  if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
-    sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayRead]"); 
-    return false; 
-  }
+  fs.seekg (start_data);
 
-  if (! read_nfloat64 (&view_angle, m_fd) ||
-      ! read_nint32 (&nDet, m_fd)) {
-      sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]");
-      return false;
-  }
+  fs.readFloat64 (view_angle);
+  fs.readInt32 (nDet);
   darray.setViewAngle (view_angle);
   //  darray.setNDet ( nDet);
   
-  for (int i = 0; i < nDet; i++) {
+  for (unsigned int i = 0; i < nDet; i++) {
       kfloat32 detval;
-      if (! read_nfloat32 (&detval, m_fd)) {
-	  sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]");
-	  return false;
-      }
+      fs.readFloat32 (detval);
       detval_ptr[i] = detval;
   }
+  if (! fs)
+    return false;
 
   return true;
 }
@@ -376,7 +421,7 @@ Projections::detarrayRead (DetectorArray& darray, const int iview)
  */
 
 bool
-Projections::detarrayWrite (const DetectorArray& darray, const int iview)
+Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
 {
   const int detval_bytes = darray.nDet() * sizeof(float);
   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
@@ -386,36 +431,26 @@ Projections::detarrayWrite (const DetectorArray& darray, const int iview)
   kfloat64 view_angle = darray.viewAngle();
   kuint32 nDet = darray.nDet();
   
-  if (m_fd < 0) {
-    sys_error (ERR_SEVERE, "Tried to write detector array with no file open [detarrayWrite]");
+  fs.seekp (start_data);
+  if (! fs) {
+    sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
     return false;
   }
 
-  if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
-      sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
-      return false;
-  }
+  fs.writeFloat64 (view_angle);
+  fs.writeInt32 (nDet);
 
-  if ( ! write_nfloat64 (&view_angle, m_fd)) {
-      sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
-      return false;
-  }
-  if ( ! write_nint32 (&nDet, m_fd)) {
-      sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
-      return false;
-  }
-  for (int i = 0; i < nDet; i++) {
-      kfloat32 detval = detval_ptr[i];
-      if ( ! write_nfloat32 (&detval, m_fd)) {
-	  sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
-	  return false;
-      }
+  for (unsigned int i = 0; i < nDet; i++) {
+    kfloat32 detval = detval_ptr[i];
+    fs.writeFloat32 (detval);
   }
 
+  if (! fs)
+    return (false);
+
   return true;
 }
 
-
 /* NAME
  *   prt_projections			Print projections data
  *
@@ -427,14 +462,16 @@ Projections::detarrayWrite (const DetectorArray& darray, const int iview)
 void
 Projections::printProjectionData (void)
 {
-  printf("Projections Print\n\n");
+  printf("Projections Data\n\n");
   printf("Description: %s\n", m_remark.c_str());
-  printf("nView = %d  nDet = %d\n", m_nView, m_nDet);
-  printf("rotStart = %8.4f   rotInc = %8.4f\n", m_rotStart, m_rotInc);
-  printf("detStart = %8.4f   detInc = %8.4f\n", m_detStart, m_detInc);
-
+  printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
+  printf("nView       = %8d           nDet = %8d\n", m_nView, m_nDet);
+  printf("focalLength = %8.4f  fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
+  printf("rotStart    = %8.4f       rotInc = %8.4f\n", m_rotStart, m_rotInc);
+  printf("detStart    = %8.4f       detInc = %8.4f\n", m_detStart, m_detInc);
   if (m_projData != NULL) {
       for (int ir = 0; ir < m_nView; ir++) {
+	  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]);
@@ -444,16 +481,166 @@ Projections::printProjectionData (void)
 }
 
 void 
-Projections::printScanInfo (void) const
+Projections::printScanInfo (ostringstream& os) const
 {
-  printf ("Number of detectors: %d\n", m_nDet);
-  printf ("    Number of views: %d\n", m_nView);
-  printf ("             Remark: %s\n", m_remark.c_str());
-  printf ("             phmLen: %f\n", m_phmLen);
-  printf ("           detStart: %f\n", m_detStart);
-  printf ("             detInc: %f\n", m_detInc);
-  printf ("           rotStart: %f\n", m_rotStart);
-  printf ("             rotInc: %f\n", m_rotInc);
+  os << "Number of detectors: " << m_nDet << "\n";
+  os << "    Number of views: " << m_nView<< "\n";
+  os << "             Remark: " << m_remark.c_str()<< "\n";
+  os << "           Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
+  os << "       Focal Length: " << m_focalLength<< "\n";
+  os << "      Field Of View: " << m_fieldOfView<< "\n";
+  os << "             phmLen: " << m_phmLen<< "\n";
+  os << "           detStart: " << m_detStart<< "\n";
+  os << "             detInc: " << m_detInc<< "\n";
+  os << "           rotStart: " << m_rotStart<< "\n";
+  os << "             rotInc: " << m_rotInc<< "\n";
 }
 
 
+
+/* 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 filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace) const
+{
+  double detInc = m_detInc;
+  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 = interpFactor;
+    spline_order = 3;
+    zoom_factor = 3;
+    n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1;
+  }
+#endif
+
+  double filterBW = 1. / detInc;
+  ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace, m_geometry, m_focalLength);
+
+  if (processSignal.fail()) {
+      sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str());
+      return false;
+  }
+
+  if (trace)
+    cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
+
+#if HAVE_SGP
+  int nVecFilter = processSignal.getNFilterPoints();
+  double plot_xaxis [nVecFilter]; 			// array for plotting 
+
+  if (trace > Trace::TRACE_CONSOLE && nVecFilter > 0)  {
+    int i;
+    double f;
+    double filterInc = processSignal.getFilterIncrement();
+    for (i = 0, f = processSignal.getFilterMin(); i < nVecFilter; i++, f += filterInc)
+      plot_xaxis[i] = f;
+
+    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::TRACE_CONSOLE) {
+    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, interpFactor);
+  if (bj.fail()) {
+    sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", bj.failMessage().c_str());
+    return false;
+  }
+
+  for (int iview = 0; iview < m_nView; iview++)  {
+    if (trace >= Trace::TRACE_CONSOLE) 
+      cout <<"Reconstructing view " << iview << "(last = " << m_nView - 1 << ")\n";
+      
+    const DetectorArray& darray = getDetectorArray (iview);
+    const DetectorValue* detval = darray.detValues();
+
+    processSignal.filterSignal (detval, filteredProj);
+
+#ifdef HAVE_BSPLINE_INTERP
+    if (interp_type == I_BSPLINE) 
+	bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
+    
+#ifdef HAVE_SGP
+    if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE) {
+	bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
+      ezplot_1d (filteredProj, n_filteredProj);
+    }
+#endif
+#endif
+
+    bj.BackprojectView (filteredProj, darray.viewAngle());
+
+#ifdef HAVE_SGP
+    if (trace >= 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];
+      fgets(str, sizeof(str), stdin);
+      if (tolower(str[0]) == 'y') {
+	break;
+      }
+    } 
+#endif  //HAVE_SGP
+  }
+
+  return true;
+}
+