X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=5ce0a8e022a3ef788ff2d5e811972a0ec7d41ee4;hp=8010912f32046b54eac5617250c1999470ac1cb4;hb=2d39e823ba389fc68e5317c422b55be006094252;hpb=5cf6874680f80d238bf34535d711dc223813f951

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 8010912..5ce0a8e 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.2 2000/06/19 17:58:20 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
@@ -137,7 +137,7 @@ Projections::deleteProjData (void)
  */
 
 bool 
-Projections::headerWrite (void)
+Projections::headerWrite (fnetorderstream& fs)
 {
   kuint32 _hsize = m_headerSize;
   kuint32 _nView = m_nView;
@@ -151,48 +151,39 @@ Projections::headerWrite (void)
   kfloat64 _detInc = m_detInc;
   kfloat64 _phmLen = m_phmLen;
   
-  if (m_fd < 0) {
-    sys_error (ERR_SEVERE, "m_fd < 0 [projections_write_header]");
+  fs.seekp(0);
+  if (! fs)
     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;
-  }
 
-  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);
-    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.writeInt32 (_hsize);
+  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.writeInt32 (_remarksize);
+  fs.write (m_remark.c_str(), _remarksize);
+
+  m_headerSize = fs.tellp();
+  _hsize = m_headerSize;
+  fs.seekp(0);
+  fs.writeInt32 (_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;
@@ -206,46 +197,43 @@ Projections::headerRead (void)
   kfloat64 _detInc;
   kfloat64 _phmLen;
   
-  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;
 
   off_t testPos;
-  if (! read_nint32 (&_hsize, m_fd)) {
-    sys_error (ERR_SEVERE, "Error reading header size , _hsize=%d [projections_read_header]", _hsize);
+  fs.readInt32 (_hsize);
+  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.readInt32 (_remarksize);
+  if (_remarksize > 100000) {
+    sys_error (ERR_SEVERE, "Insane _remarksize %d [projections::headerRead]", _remarksize);
     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;
-      }
+
+  if (! fs) {
+      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;
   }
@@ -264,31 +252,26 @@ Projections::headerRead (void)
   return true;
 }
 
-
 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);
+
+  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;
 }
 
@@ -296,30 +279,25 @@ Projections::read (const char* filename)
 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);
 
-    headerWrite ();
+  if (! fs) {
+    sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
+    return false;
+  }
+  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;
 }
@@ -334,7 +312,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 */;
@@ -344,32 +322,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; 
-  } 
+  fs.seekg (start_data);
 
-  if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
-    sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayRead]"); 
-    return false; 
-  }
-
-  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++) {
       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;
 }
@@ -390,7 +356,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 */;
@@ -400,36 +366,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;
-      }
+    kfloat32 detval = detval_ptr[i];
+    fs.writeFloat32 (detval);
   }
 
+  if (! fs)
+    return (false);
+
   return true;
 }
 
-
 /* NAME
  *   prt_projections			Print projections data
  *
@@ -471,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;
+}
+