r142: *** empty log message ***
[ctsim.git] / src / pjrec.cpp
index 6f9fb7a3fbc3349a9aa7fa210ee7242ebd2f7951..d6a229caeb35cdc7ffb0ff3ccc24716a9dc1159c 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pjrec.cpp,v 1.3 2000/06/29 12:39:46 kevin Exp $
+**  $Id: pjrec.cpp,v 1.10 2000/07/11 10:32:44 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
 #include "timer.h"
 
 
-enum {O_INTERP, O_FILTER, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+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};
 
 static struct option my_options[] =
 {
   {"interp", 1, 0, O_INTERP},
+  {"preinterpolation-factor", 1, 0, O_PREINTERPOLATION_FACTOR},
   {"filter", 1, 0, O_FILTER},
+  {"filter-method", 1, 0, O_FILTER_METHOD},
+  {"zeropad", 1, 0, O_ZEROPAD},
   {"filter-param", 1, 0, O_FILTER_PARAM},
   {"backproj", 1, 0, O_BACKPROJ},
   {"trace", 1, 0, O_TRACE},
@@ -52,45 +55,59 @@ 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;
+  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;
 #if HAVE_BSPLINE_INTERP
-  cout << "       bspline     B-spline interpolation" << endl;
+  cout << "    bspline     B-spline interpolation" << endl;
 #endif
-  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_hamming   Abs * Hamming" << endl;
-  cout << "       shepp         Shepp-Logan" << endl;
-  cout << "       bandlimit     Bandlimiting" << endl;
-  cout << "       sinc          Sinc" << endl;
-  cout << "       cos           Cosine" << endl;
-  cout << "       triangle      Triangle" << endl;
-  cout << "       hamming       Hamming" << endl;
-  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 << "    --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 << "    --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;
+  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_cos       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 << "    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";
+#if HAVE_FFTW
+  cout << "    fftw             Fast Fourier Transform West library\n";
+  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 << "  --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 << "     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 << "  --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;
 }
 
 
@@ -111,6 +128,7 @@ pjrec_main (int argc, char * argv[])
   char *endptr;
   int optVerbose = 0;
   int optDebug = 0;
+  int optZeroPad = 0;
   int optTrace = TRACE_NONE;
   double optFilterParam = -1;
   string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR;
@@ -118,7 +136,7 @@ pjrec_main (int argc, char * argv[])
   string optInterpName = Backprojector::INTERP_LINEAR_STR;
   string optBackprojName = Backprojector::BPROJ_IDIFF2_STR;
   //  string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
-  int optInterpParam = 1;
+  int optPreinterpolationFactor = 1;
   int nx, ny;
 #ifdef HAVE_MPI
   ImageFile* imLocal;
@@ -144,9 +162,19 @@ pjrec_main (int argc, char * argv[])
        case O_INTERP:
          optInterpName = optarg;
          break;
+       case O_PREINTERPOLATION_FACTOR:
+         optPreinterpolationFactor = strtol(optarg, &endptr, 10);
+         if (endptr != optarg + strlen(optarg)) {
+           pjrec_usage(argv[0]);
+           return(1);
+         }
+         break;
        case O_FILTER:
          optFilterName = optarg;
          break;
+        case O_FILTER_METHOD:
+         optFilterMethodName = optarg;
+          break;
        case O_BACKPROJ:
          optBackprojName = optarg;
          break;
@@ -154,6 +182,14 @@ pjrec_main (int argc, char * argv[])
          optFilterParam = strtod(optarg, &endptr);
          if (endptr != optarg + strlen(optarg)) {
            pjrec_usage(argv[0]);
+           return(1);
+         }
+         break;
+       case O_ZEROPAD:
+         optZeroPad = strtol(optarg, &endptr, 10);
+         if (endptr != optarg + strlen(optarg)) {
+           pjrec_usage(argv[0]);
+           return(1);
          }
          break;
        case O_VERBOSE:
@@ -163,7 +199,7 @@ pjrec_main (int argc, char * argv[])
          optDebug = 1;
          break;
        case O_TRACE:
-         if ((optTrace = opt_set_trace(optarg)) < 0) {
+         if ((optTrace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
            pjrec_usage(argv[0]);
            return (1);
          }
@@ -204,7 +240,7 @@ pjrec_main (int argc, char * argv[])
       filterDesc << optFilterName;
 
     ostringstream label;
-    label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", " << optBackprojName;
+    label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolation=" << optPreinterpolationFactor << ", " << optBackprojName;
     remark = label.str();
   
     if (optVerbose)
@@ -229,12 +265,14 @@ 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.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
-  mpiWorld.getComm().Bcast (&optInterpParam, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 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);
@@ -272,7 +310,7 @@ pjrec_main (int argc, char * argv[])
 
 #ifdef HAVE_MPI
   TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
-  projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
+  projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
   if (optVerbose)
       timerReconstruct.timerEndAndReport ("Time to reconstruct");
 
@@ -281,7 +319,7 @@ pjrec_main (int argc, char * argv[])
   if (optVerbose)
       timerReduce.timerEndAndReport ("Time to reduce image");
 #else
-  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
+  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
 #endif
 
 #ifdef HAVE_MPI