r186: *** empty log message ***
[ctsim.git] / tools / pjrec.cpp
index 03ee9e874520770e76e160dc5673a739c0165708..9fcb8903a6de987e3a130f7ae720eeb28883dc45 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pjrec.cpp,v 1.5 2000/07/22 15:45:33 kevin Exp $
+**  $Id: pjrec.cpp,v 1.16 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
@@ -29,7 +29,7 @@
 #include "timer.h"
 
 
-enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+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[] =
 {
@@ -38,6 +38,7 @@ static struct option my_options[] =
   {"filter", 1, 0, O_FILTER},
   {"filter-method", 1, 0, O_FILTER_METHOD},
   {"zeropad", 1, 0, O_ZEROPAD},
+  {"filter-generation", 1, 0, O_FILTER_GENERATION},
   {"filter-param", 1, 0, O_FILTER_PARAM},
   {"backproj", 1, 0, O_BACKPROJ},
   {"trace", 1, 0, O_TRACE},
@@ -48,6 +49,7 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
+static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.16 2000/08/31 08:38:58 kevin Exp $";
 
 void 
 pjrec_usage (const char *program)
@@ -70,12 +72,12 @@ pjrec_usage (const char *program)
   cout << "  --filter       Filter name" << endl;
   cout << "    abs_bandlimit Abs * Bandlimiting (default)" << endl;
   cout << "    abs_sinc      Abs * Sinc" << endl;
-  cout << "    abs_cos       Abs * Cosine" << 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 << "    cos           Cosine" << endl;
+  cout << "    cosine        Cosine" << endl;
   cout << "    triangle      Triangle" << endl;
   cout << "    hamming       Hamming" << endl;
   cout << "  --filter-method  Filter method before backprojections\n";;
@@ -89,6 +91,9 @@ pjrec_usage (const char *program)
 #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;
@@ -98,12 +103,12 @@ pjrec_usage (const char *program)
   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 << "     text      Text level tracing" << endl;
-  cout << "     phm       Trace phantom" << endl;
-  cout << "     rays      Trace allrays" << endl;
-  cout << "     plot      Trace plotting" << endl;
-  cout << "     clipping  Trace clipping" << 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;
@@ -112,7 +117,7 @@ pjrec_usage (const char *program)
 
 
 #ifdef HAVE_MPI
-static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug);
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const bool bDebug);
 static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal);
 #endif
 
@@ -122,21 +127,22 @@ pjrec_main (int argc, char * argv[])
 {
   Projections projGlobal;
   ImageFile* imGlobal = NULL;
-  char* filenameProj = NULL;
-  char* filenameImage = NULL;
-  string remark;
-  char *endptr;
-  int optVerbose = 0;
-  int optDebug = 0;
-  int optZeroPad = 0;
-  int optTrace = TRACE_NONE;
-  double optFilterParam = -1;
-  string optFilterName (SignalFilter::convertFilterIDToName (SignalFilter::FILTER_ABS_BANDLIMIT));
-  string optFilterMethodName (SignalFilter::convertFilterMethodIDToName (SignalFilter::FILTER_METHOD_CONVOLUTION));
-  string optInterpName (Backprojector::convertInterpIDToName (Backprojector::INTERP_LINEAR));
-  string optBackprojName (Backprojector::convertBackprojectIDToName (Backprojector::BPROJ_IDIFF3));
-  int optPreinterpolationFactor = 1;
+  char* pszFilenameProj = NULL;
+  char* pszFilenameImage = NULL;
+  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));
+  int iOptPreinterpolationFactor = 1;
   int nx, ny;
+  char *endptr;
 #ifdef HAVE_MPI
   ImageFile* imLocal;
   int mpi_nview, mpi_ndet;
@@ -159,53 +165,56 @@ pjrec_main (int argc, char * argv[])
       switch (c)
        {
        case O_INTERP:
-         optInterpName = optarg;
+         sOptInterpName = optarg;
          break;
        case O_PREINTERPOLATION_FACTOR:
-         optPreinterpolationFactor = strtol(optarg, &endptr, 10);
+         iOptPreinterpolationFactor = strtol(optarg, &endptr, 10);
          if (endptr != optarg + strlen(optarg)) {
            pjrec_usage(argv[0]);
            return(1);
          }
          break;
        case O_FILTER:
-         optFilterName = optarg;
+         sOptFilterName = optarg;
          break;
         case O_FILTER_METHOD:
-         optFilterMethodName = optarg;
+         sOptFilterMethodName = optarg;
           break;
-       case O_BACKPROJ:
-         optBackprojName = optarg;
+       case O_FILTER_GENERATION:
+         sOptFilterGenerationName = optarg;
          break;
        case O_FILTER_PARAM:
-         optFilterParam = strtod(optarg, &endptr);
+         dOptFilterParam = strtod(optarg, &endptr);
          if (endptr != optarg + strlen(optarg)) {
            pjrec_usage(argv[0]);
            return(1);
          }
          break;
        case O_ZEROPAD:
-         optZeroPad = strtol(optarg, &endptr, 10);
+         iOptZeropad = strtol(optarg, &endptr, 10);
          if (endptr != optarg + strlen(optarg)) {
            pjrec_usage(argv[0]);
            return(1);
          }
          break;
+       case O_BACKPROJ:
+         sOptBackprojectName = optarg;
+         break;
        case O_VERBOSE:
-         optVerbose = 1;
+         bOptVerbose = true;
          break;
        case O_DEBUG:
-         optDebug = 1;
+         bOptDebug = true;
          break;
        case O_TRACE:
-         if ((optTrace = TraceLevel::convertTraceNameToID(optarg)) == TRACE_INVALID) {
+         if ((optTrace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
            pjrec_usage(argv[0]);
            return (1);
          }
          break;
         case O_VERSION:
 #ifdef VERSION
-         cout <<  "Version " <<  VERSION << endl;
+         cout <<  "Version " <<  VERSION << endl << g_szIdStr << endl;
 #else
           cout << "Unknown version number" << endl;
 #endif
@@ -225,34 +234,37 @@ pjrec_main (int argc, char * argv[])
       return (1);
     }
 
-    filenameProj = argv[optind];
+    pszFilenameProj = argv[optind];
   
-    filenameImage = argv[optind + 1];
+    pszFilenameImage = argv[optind + 1];
   
     nx = strtol(argv[optind + 2], &endptr, 10);
     ny = strtol(argv[optind + 3], &endptr, 10);
   
     ostringstream filterDesc;
-    if (optFilterParam >= 0)
-      filterDesc << optFilterName << ": alpha=" << optFilterParam; 
+    if (dOptFilterParam >= 0)
+      filterDesc << sOptFilterName << ": alpha=" << dOptFilterParam; 
     else
-      filterDesc << optFilterName;
+      filterDesc << sOptFilterName;
 
     ostringstream label;
-    label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolationFactor=" << optPreinterpolationFactor << ", " << optBackprojName;
-    remark = label.str();
+    label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << sOptInterpName << ", preinterpolationFactor=" << iOptPreinterpolationFactor << ", " << sOptBackprojectName;
+    sRemark = label.str();
   
-    if (optVerbose)
-      cout << "Remark: " << remark << endl;
+    if (bOptVerbose)
+      cout << "SRemark: " << sRemark << endl;
 #ifdef HAVE_MPI
   }
 #endif
 
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0) {
-    projGlobal.read (filenameProj);
-    if (optVerbose)
-      projGlobal.printScanInfo();
+    projGlobal.read (pszFilenameProj);
+    if (bOptVerbose) {
+      ostringstream os;
+      projGlobal.printScanInfo (os);
+      cout << os.str();
+    }
 
     mpi_ndet = projGlobal.nDet();
     mpi_nview = projGlobal.nView();
@@ -262,16 +274,16 @@ pjrec_main (int argc, char * argv[])
   }
 
   TimerCollectiveMPI timerBcast (mpiWorld.getComm());
-  mpiWorld.BcastString (optBackprojName);
-  mpiWorld.BcastString (optFilterName);
-  mpiWorld.BcastString (optFilterMethodName);
-  mpiWorld.BcastString (optInterpName);
-  mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
+  mpiWorld.BcastString (sOptBackprojectName);
+  mpiWorld.BcastString (sOptFilterName);
+  mpiWorld.BcastString (sOptFilterMethodName);
+  mpiWorld.BcastString (sOptInterpName);
+  mpiWorld.getComm().Bcast (&bOptVerbose, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&bOptDebug, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
-  mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&dOptFilterParam, 1, MPI::DOUBLE, 0);
+  mpiWorld.getComm().Bcast (&iOptZeropad, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&iOptPreinterpolationFactor, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
@@ -279,7 +291,7 @@ pjrec_main (int argc, char * argv[])
   mpiWorld.getComm().Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0);
   mpiWorld.getComm().Bcast (&nx, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&ny, 1, MPI::INT, 0);
-  if (optVerbose)
+  if (bOptVerbose)
       timerBcast.timerEndAndReport ("Time to broadcast variables");
 
   mpiWorld.setTotalWorkUnits (mpi_nview);
@@ -290,8 +302,8 @@ pjrec_main (int argc, char * argv[])
   projLocal.setRotInc (mpi_rotinc);
 
   TimerCollectiveMPI timerScatter (mpiWorld.getComm());
-  ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug);
-  if (optVerbose)
+  ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, bOptDebug);
+  if (bOptVerbose)
       timerScatter.timerEndAndReport ("Time to scatter projections");
 
   if (mpiWorld.getRank() == 0) {
@@ -300,37 +312,40 @@ pjrec_main (int argc, char * argv[])
 
   imLocal = new ImageFile (nx, ny);
 #else
-  projGlobal.read (filenameProj);
-  if (optVerbose)
-    projGlobal.printScanInfo();
+  projGlobal.read (pszFilenameProj);
+  if (bOptVerbose) {
+    ostringstream os;
+    projGlobal.printScanInfo(os);
+    cout << os.str();
+  }
 
   imGlobal = new ImageFile (nx, ny);
 #endif
 
 #ifdef HAVE_MPI
   TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
-  projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
-  if (optVerbose)
+  projLocal.reconstruct (*imLocal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace);
+  if (bOptVerbose)
       timerReconstruct.timerEndAndReport ("Time to reconstruct");
 
   TimerCollectiveMPI timerReduce (mpiWorld.getComm());
   ReduceImageMPI (mpiWorld, imLocal, imGlobal);
-  if (optVerbose)
+  if (bOptVerbose)
       timerReduce.timerEndAndReport ("Time to reduce image");
 #else
-  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
+  projGlobal.reconstruct (*imGlobal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace);
 #endif
 
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0)
 #endif
     {
-      double calcTime = timerProgram.timerEnd();
+      double dCalcTime = timerProgram.timerEnd();
       imGlobal->labelAdd (projGlobal.getLabel());
-      imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime);
-      imGlobal->fileWrite (filenameImage);
-      if (optVerbose)
-       cout << "Run time: " << calcTime << " seconds" << endl;
+      imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, sRemark.c_str(), dCalcTime);
+      imGlobal->fileWrite (pszFilenameImage);
+      if (bOptVerbose)
+       cout << "Run time: " << dCalcTime << " seconds" << endl;
     }
 #ifdef HAVE_MPI
   MPI::Finalize();
@@ -346,7 +361,7 @@ pjrec_main (int argc, char * argv[])
 //////////////////////////////////////////////////////////////////////////////////////
 
 #ifdef HAVE_MPI
-static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug)
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const bool bOptDebug)
 {
   if (mpiWorld.getRank() == 0) {
     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
@@ -381,7 +396,7 @@ ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal)
 {
   ImageFileArray vLocal = imLocal->getArray();
 
-  for (int ix = 0; ix < imLocal->nx(); ix++) {
+  for (unsigned int ix = 0; ix < imLocal->nx(); ix++) {
     void *recvbuf = NULL;
     if (mpiWorld.getRank() == 0) {
       ImageFileArray vGlobal = imGlobal->getArray();