r7061: initial property settings
[ctsim.git] / tools / pjrec.cpp
index 157d6584c1cdd28e8c06d3041817f0ed8a1503c7..ea230489c36429778cb9f133499624885a7d49d8 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pjrec.cpp,v 1.15 2000/08/27 20:32:55 kevin Exp $
+**  $Id$
 **
 **  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
@@ -28,7 +28,6 @@
 #include "ct.h"
 #include "timer.h"
 
-
 enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_FILTER_GENERATION, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
 
 static struct option my_options[] =
@@ -49,70 +48,69 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.15 2000/08/27 20:32:55 kevin Exp $";
+static const char* g_szIdStr = "$Id$";
 
 void 
 pjrec_usage (const char *program)
 {
-  cout << "usage: " << fileBasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl;
-  cout << "Image reconstruction from raysum projections" << endl;
-  cout << endl;
-  cout << "  raysum-file     Input raysum file" << endl;
-  cout << "  image-file      Output image file in SDF2D format" << endl;
-  cout << "  nx-image        Number of columns in output image" << endl;
-  cout << "  ny-image        Number of rows in output image" << endl;
-  cout << "  --interp        Interpolation method during backprojection" << endl;
-  cout << "    nearest     Nearest neighbor interpolation" << endl;
-  cout << "    linear      Linear interpolation" << endl;
+  std::cout << "usage: " << fileBasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << std::endl;
+  std::cout << "Image reconstruction from raysum projections" << std::endl;
+  std::cout << std::endl;
+  std::cout << "  raysum-file     Input raysum file" << std::endl;
+  std::cout << "  image-file      Output image file in SDF2D format" << std::endl;
+  std::cout << "  nx-image        Number of columns in output image" << std::endl;
+  std::cout << "  ny-image        Number of rows in output image" << std::endl;
+  std::cout << "  --interp        Interpolation method during backprojection" << std::endl;
+  std::cout << "    nearest         Nearest neighbor interpolation" << std::endl;
+  std::cout << "    linear          Linear interpolation (default)" << std::endl;
+  std::cout << "    cubic           Cubic interpolation\n";
 #if HAVE_BSPLINE_INTERP
-  cout << "    bspline     B-spline interpolation" << endl;
+  std::cout << "    bspline         B-spline interpolation" << std::endl;
 #endif
-  cout << "  --preinterpolate  Preinterpolation factor (default = 1)\n";
-  cout << "                    Used only with frequency-based filtering\n";
-  cout << "  --filter       Filter name" << endl;
-  cout << "    abs_bandlimit Abs * Bandlimiting (default)" << endl;
-  cout << "    abs_sinc      Abs * Sinc" << endl;
-  cout << "    abs_cosine    Abs * Cosine" << endl;
-  cout << "    abs_hamming   Abs * Hamming" << endl;
-  cout << "    shepp         Shepp-Logan" << endl;
-  cout << "    bandlimit     Bandlimiting" << endl;
-  cout << "    sinc          Sinc" << endl;
-  cout << "    cosine        Cosine" << endl;
-  cout << "    triangle      Triangle" << endl;
-  cout << "    hamming       Hamming" << endl;
-  cout << "  --filter-method  Filter method before backprojections\n";;
-  cout << "    convolution      Spatial filtering (default)\n";
-  cout << "    fourier          Frequency filtering with discete fourier\n";
-  cout << "    fourier_table    Frequency filtering with table lookup fourier\n";
-  cout << "    fft              Fast Fourier Transform\n";
+  std::cout << "  --preinterpolate  Preinterpolation factor (default = 1)\n";
+  std::cout << "                    Used only with frequency-based filtering\n";
+  std::cout << "  --filter       Filter name" << std::endl;
+  std::cout << "    abs_bandlimit  Abs * Bandlimiting (default)" << std::endl;
+  std::cout << "    abs_sinc       Abs * Sinc" << std::endl;
+  std::cout << "    abs_cosine     Abs * Cosine" << std::endl;
+  std::cout << "    abs_hamming    Abs * Hamming" << std::endl;
+  std::cout << "    abs_hanning    Abs * Hanning" << std::endl;
+  std::cout << "    shepp          Shepp-Logan" << std::endl;
+  std::cout << "    bandlimit      Bandlimiting" << std::endl;
+  std::cout << "    sinc           Sinc" << std::endl;
+  std::cout << "    cosine         Cosine" << std::endl;
+  std::cout << "    triangle       Triangle" << std::endl;
+  std::cout << "    hamming        Hamming" << std::endl;
+  std::cout << "    hanning        Hanning" << std::endl;
+  std::cout << "  --filter-method  Filter method before backprojections\n";;
+  std::cout << "    convolution      Spatial filtering (default)\n";
+  std::cout << "    fourier          Frequency filtering with discete fourier\n";
+  std::cout << "    fourier_table    Frequency filtering with table lookup fourier\n";
+  std::cout << "    fft              Fast Fourier Transform\n";
 #if HAVE_FFTW
-  cout << "    fftw             Fast Fourier Transform West library\n";
-  cout << "    rfftw            Fast Fourier Transform West (real-mode) library\n";
+  std::cout << "    fftw             Fast Fourier Transform West library\n";
+  std::cout << "    rfftw            Fast Fourier Transform West (real-mode) library\n";
 #endif
-  cout << "  --zeropad n   Set zeropad level (default = 0)\n";
-  cout << "                set n to number of powers to two to pad\n";
-  cout << "  --filter-generation  Filter Generation mode\n";
-  cout << "    direct       Use direct filter in spatial or frequency domain (default)\n";
-  cout << "    inverse_fourier  Use inverse fourier transform of inverse filter\n";
-  cout << "  --backproj    Backprojection Method" << endl;
-  cout << "    trig        Trigometric functions at every point" << endl;
-  cout << "    table       Trigometric functions with precalculated table" << endl;
-  cout << "    diff        Difference method" << endl;
-  cout << "    diff2       Optimized difference method (default)" << endl;
-  cout << "    idiff2      Optimized difference method with integer math" << endl;
-  cout << "    idiff3      Highly-optimized difference method with integer math" << endl;
-  cout << "  --filter-param Alpha level for Hamming filter" << endl;
-  cout << "  --trace        Set tracing to level" << endl;
-  cout << "     none        No tracing (default)" << endl;
-  cout << "     console     Text level tracing" << endl;
-  cout << "     phantom     Trace phantom" << endl;
-  cout << "     proj        Trace allrays" << endl;
-  cout << "     plot        Trace plotting" << endl;
-  cout << "     clipping    Trace clipping" << endl;
-  cout << "  --verbose      Turn on verbose mode" << endl;
-  cout << "  --debug        Turn on debug mode" << endl;
-  cout << "  --version      Print version" << endl;
-  cout << "  --help         Print this help message" << endl;
+  std::cout << "  --zeropad n   Set zeropad level (default = 0)\n";
+  std::cout << "                set n to number of powers to two to pad\n";
+  std::cout << "  --filter-generation  Filter Generation mode\n";
+  std::cout << "    direct       Use direct filter in spatial or frequency domain (default)\n";
+  std::cout << "    inverse_fourier  Use inverse fourier transform of inverse filter\n";
+  std::cout << "  --backproj    Backprojection Method" << std::endl;
+  std::cout << "    trig        Trigometric functions at every point" << std::endl;
+  std::cout << "    table       Trigometric functions with precalculated table" << std::endl;
+  std::cout << "    diff        Difference method" << std::endl;
+  std::cout << "    diff2       Optimized difference method (default)" << std::endl;
+  std::cout << "    idiff2      Optimized difference method with integer math" << std::endl;
+  std::cout << "    idiff3      Highly-optimized difference method with integer math" << std::endl;
+  std::cout << "  --filter-param Alpha level for Hamming filter" << std::endl;
+  std::cout << "  --trace        Set tracing to level" << std::endl;
+  std::cout << "     none        No tracing (default)" << std::endl;
+  std::cout << "     console     Text level tracing" << std::endl;
+  std::cout << "  --verbose      Turn on verbose mode" << std::endl;
+  std::cout << "  --debug        Turn on debug mode" << std::endl;
+  std::cout << "  --version      Print version" << std::endl;
+  std::cout << "  --help         Print this help message" << std::endl;
 }
 
 
@@ -123,23 +121,23 @@ static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* i
 
 
 int 
-pjrec_main (int argc, char * argv[])
+pjrec_main (int argc, char * const argv[])
 {
   Projections projGlobal;
   ImageFile* imGlobal = NULL;
   char* pszFilenameProj = NULL;
   char* pszFilenameImage = NULL;
-  string sRemark;
+  std::string sRemark;
   bool bOptVerbose = false;
   bool bOptDebug = 1;
   int iOptZeropad = 1;
   int optTrace = Trace::TRACE_NONE;
   double dOptFilterParam = -1;
-  string sOptFilterName (SignalFilter::convertFilterIDToName (SignalFilter::FILTER_ABS_BANDLIMIT));
-  string sOptFilterMethodName (ProcessSignal::convertFilterMethodIDToName (ProcessSignal::FILTER_METHOD_CONVOLUTION));
-  string sOptFilterGenerationName (ProcessSignal::convertFilterGenerationIDToName (ProcessSignal::FILTER_GENERATION_DIRECT));
-  string sOptInterpName (Backprojector::convertInterpIDToName (Backprojector::INTERP_LINEAR));
-  string sOptBackprojectName (Backprojector::convertBackprojectIDToName (Backprojector::BPROJ_IDIFF3));
+  std::string sOptFilterName (SignalFilter::convertFilterIDToName (SignalFilter::FILTER_ABS_BANDLIMIT));
+  std::string sOptFilterMethodName (ProcessSignal::convertFilterMethodIDToName (ProcessSignal::FILTER_METHOD_CONVOLUTION));
+  std::string sOptFilterGenerationName (ProcessSignal::convertFilterGenerationIDToName (ProcessSignal::FILTER_GENERATION_DIRECT));
+  std::string sOptInterpName (Backprojector::convertInterpIDToName (Backprojector::INTERP_LINEAR));
+  std::string sOptBackprojectName (Backprojector::convertBackprojectIDToName (Backprojector::BPROJ_IDIFF));
   int iOptPreinterpolationFactor = 1;
   int nx, ny;
   char *endptr;
@@ -214,9 +212,9 @@ pjrec_main (int argc, char * argv[])
          break;
         case O_VERSION:
 #ifdef VERSION
-         cout <<  "Version " <<  VERSION << endl << g_szIdStr << endl;
+         std::cout <<  "Version " <<  VERSION << std::endl << g_szIdStr << std::endl;
 #else
-          cout << "Unknown version number" << endl;
+          std::cout << "Unknown version number" << std::endl;
 #endif
          return (0);
        case O_HELP:
@@ -241,18 +239,18 @@ pjrec_main (int argc, char * argv[])
     nx = strtol(argv[optind + 2], &endptr, 10);
     ny = strtol(argv[optind + 3], &endptr, 10);
   
-    ostringstream filterDesc;
+    std::ostringstream filterDesc;
     if (dOptFilterParam >= 0)
       filterDesc << sOptFilterName << ": alpha=" << dOptFilterParam; 
     else
       filterDesc << sOptFilterName;
 
-    ostringstream label;
+    std::ostringstream label;
     label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << sOptInterpName << ", preinterpolationFactor=" << iOptPreinterpolationFactor << ", " << sOptBackprojectName;
     sRemark = label.str();
   
     if (bOptVerbose)
-      cout << "SRemark: " << sRemark << endl;
+      std::cout << "SRemark: " << sRemark << std::endl;
 #ifdef HAVE_MPI
   }
 #endif
@@ -260,8 +258,11 @@ pjrec_main (int argc, char * argv[])
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0) {
     projGlobal.read (pszFilenameProj);
-    if (bOptVerbose)
-      projGlobal.printScanInfo();
+    if (bOptVerbose) {
+      ostringstream os;
+      projGlobal.printScanInfo (os);
+      std::cout << os.str();
+    }
 
     mpi_ndet = projGlobal.nDet();
     mpi_nview = projGlobal.nView();
@@ -310,15 +311,25 @@ pjrec_main (int argc, char * argv[])
   imLocal = new ImageFile (nx, ny);
 #else
   projGlobal.read (pszFilenameProj);
-  if (bOptVerbose)
-    projGlobal.printScanInfo();
+  if (bOptVerbose) {
+    std::ostringstream os;
+    projGlobal.printScanInfo(os);
+    std::cout << os.str();
+  }
 
   imGlobal = new ImageFile (nx, ny);
 #endif
 
 #ifdef HAVE_MPI
   TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
-  projLocal.reconstruct (*imLocal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace);
+
+  Reconstructor reconstruct (projLocal, *imLocal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace);
+  if (reconstruct.fail()) {
+    std::cout << reconstruct.failMessage();
+    return (1);
+  }
+  reconstruct.reconstructAllViews();
+  
   if (bOptVerbose)
       timerReconstruct.timerEndAndReport ("Time to reconstruct");
 
@@ -327,7 +338,12 @@ pjrec_main (int argc, char * argv[])
   if (bOptVerbose)
       timerReduce.timerEndAndReport ("Time to reduce image");
 #else
-  projGlobal.reconstruct (*imGlobal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace);
+  Reconstructor reconstruct (projGlobal, *imGlobal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace);
+  if (reconstruct.fail()) {
+    std::cout << reconstruct.failMessage();
+    return (1);
+  }
+  reconstruct.reconstructAllViews();
 #endif
 
 #ifdef HAVE_MPI
@@ -339,7 +355,7 @@ pjrec_main (int argc, char * argv[])
       imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, sRemark.c_str(), dCalcTime);
       imGlobal->fileWrite (pszFilenameImage);
       if (bOptVerbose)
-       cout << "Run time: " << dCalcTime << " seconds" << endl;
+       std::cout << "Run time: " << dCalcTime << " seconds" << std::endl;
     }
 #ifdef HAVE_MPI
   MPI::Finalize();
@@ -412,9 +428,9 @@ main (int argc, char* argv[])
   try {
     retval = pjrec_main(argc, argv);
   } catch (exception e) {
-    cerr << "Exception: " << e.what() << endl;
+         std::cerr << "Exception: " << e.what() << std::endl;
   } catch (...) {
-    cerr << "Unknown exception" << endl;
+         std::cerr << "Unknown exception" << std::endl;
   }
 
   return (retval);