X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=src%2Fpjrec.cpp;fp=src%2Fpjrec.cpp;h=b5760a79361252a0047751acdc2e2cc5537b2906;hb=2f3d6e2580db607105bb072b13e4aff453ae4495;hp=0000000000000000000000000000000000000000;hpb=08f34bf3ba14d4f436f4d2ef0ee5af1d6eb266ac;p=ctsim.git

diff --git a/src/pjrec.cpp b/src/pjrec.cpp
new file mode 100644
index 0000000..b5760a7
--- /dev/null
+++ b/src/pjrec.cpp
@@ -0,0 +1,376 @@
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+**   Name:          pjrec.cpp
+**   Purpose:       Reconstruct an image from projections
+**   Programmer:    Kevin Rosenberg
+**   Date Started:  Aug 1984
+**
+**  This is part of the CTSim program
+**  Copyright (C) 1983-2000 Kevin Rosenberg
+**
+**  $Id: pjrec.cpp,v 1.1 2000/06/26 21:15:24 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
+**  published by the Free Software Foundation.
+**
+**  This program is distributed in the hope that it will be useful,
+**  but WITHOUT ANY WARRANTY; without even the implied warranty of
+**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+**  GNU General Public License for more details.
+**
+**  You should have received a copy of the GNU General Public License
+**  along with this program; if not, write to the Free Software
+**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+******************************************************************************/
+
+#include "ct.h"
+#include "timer.h"
+
+
+enum {O_INTERP, O_FILTER, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+
+static struct option my_options[] =
+{
+  {"interp", 1, 0, O_INTERP},
+  {"filter", 1, 0, O_FILTER},
+  {"filter-param", 1, 0, O_FILTER_PARAM},
+  {"backproj", 1, 0, O_BACKPROJ},
+  {"trace", 1, 0, O_TRACE},
+  {"debug", 0, 0, O_DEBUG},
+  {"verbose", 0, 0, O_VERBOSE},
+  {"help", 0, 0, O_HELP},
+  {"version", 0, 0, O_VERSION},
+  {0, 0, 0, 0}
+};
+
+
+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;
+#if HAVE_BSPLINE_INTERP
+  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;
+}
+
+
+#ifdef HAVE_MPI
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug);
+static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal);
+#endif
+
+
+int 
+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 optTrace = TRACE_NONE;
+  double optFilterParam = -1;
+  string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR;
+  string optInterpName = Backprojector::INTERP_LINEAR_STR;
+  string optBackprojName = Backprojector::BPROJ_IDIFF2_STR;
+  //  string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
+  int optInterpParam = 1;
+  int nx, ny;
+#ifdef HAVE_MPI
+  ImageFile* imLocal;
+  int mpi_nview, mpi_ndet;
+  double mpi_detinc, mpi_rotinc, mpi_phmlen;
+  MPIWorld mpiWorld (argc, argv);
+#endif
+
+  Timer timerProgram;
+
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0) {
+#endif
+    while (1) {
+      int c = getopt_long(argc, argv, "", my_options, NULL);
+      char *endptr = NULL;
+      
+      if (c == -1)
+	break;
+      
+      switch (c)
+	{
+	case O_INTERP:
+	  optInterpName = optarg;
+	  break;
+	case O_FILTER:
+	  optFilterName = optarg;
+	  break;
+	case O_BACKPROJ:
+	  optBackprojName = optarg;
+	  break;
+	case O_FILTER_PARAM:
+	  optFilterParam = strtod(optarg, &endptr);
+	  if (endptr != optarg + strlen(optarg)) {
+	    pjrec_usage(argv[0]);
+	  }
+	  break;
+	case O_VERBOSE:
+	  optVerbose = 1;
+	  break;
+	case O_DEBUG:
+	  optDebug = 1;
+	  break;
+	case O_TRACE:
+	  if ((optTrace = opt_set_trace(optarg)) < 0) {
+	    pjrec_usage(argv[0]);
+	    return (1);
+	  }
+	  break;
+        case O_VERSION:
+#ifdef VERSION
+	  cout <<  "Version " <<  VERSION << endl;
+#else
+          cout << "Unknown version number" << endl;
+#endif
+	  return (0);
+	case O_HELP:
+	case '?':
+	  pjrec_usage(argv[0]);
+	  return (0);
+	default:
+	  pjrec_usage(argv[0]);
+	  return (1);
+	}
+    }
+  
+    if (optind + 4 != argc) {
+      pjrec_usage(argv[0]);
+      return (1);
+    }
+
+    filenameProj = argv[optind];
+  
+    filenameImage = 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; 
+    else
+      filterDesc << optFilterName;
+
+    ostringstream label;
+    label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", " << optBackprojName;
+    remark = label.str();
+  
+    if (optVerbose)
+      cout << "Remark: " << remark << endl;
+#ifdef HAVE_MPI
+  }
+#endif
+
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0) {
+    projGlobal.read (filenameProj);
+    if (optVerbose)
+      projGlobal.printScanInfo();
+
+    mpi_ndet = projGlobal.nDet();
+    mpi_nview = projGlobal.nView();
+    mpi_detinc = projGlobal.detInc();
+    mpi_phmlen = projGlobal.phmLen();
+    mpi_rotinc = projGlobal.rotInc();
+  }
+
+  TimerCollectiveMPI timerBcast (mpiWorld.getComm());
+  mpiWorld.BcastString (optBackprojName);
+  mpiWorld.BcastString (optFilterName);
+  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 (&mpi_ndet, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
+  mpiWorld.getComm().Bcast (&mpi_phmlen, 1, MPI::DOUBLE, 0);
+  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)
+      timerBcast.timerEndAndReport ("Time to broadcast variables");
+
+  mpiWorld.setTotalWorkUnits (mpi_nview);
+
+  Projections projLocal (mpiWorld.getMyLocalWorkUnits(), mpi_ndet);
+  projLocal.setDetInc (mpi_detinc);
+  projLocal.setPhmLen (mpi_phmlen);
+  projLocal.setRotInc (mpi_rotinc);
+
+  TimerCollectiveMPI timerScatter (mpiWorld.getComm());
+  ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug);
+  if (optVerbose)
+      timerScatter.timerEndAndReport ("Time to scatter projections");
+
+  if (mpiWorld.getRank() == 0) {
+    imGlobal = new ImageFile (nx, ny);
+  }
+
+  imLocal = new ImageFile (nx, ny);
+#else
+  projGlobal.read (filenameProj);
+  if (optVerbose)
+    projGlobal.printScanInfo();
+
+  imGlobal = new ImageFile (nx, ny);
+#endif
+
+#ifdef HAVE_MPI
+  TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
+  projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
+  if (optVerbose)
+      timerReconstruct.timerEndAndReport ("Time to reconstruct");
+
+  TimerCollectiveMPI timerReduce (mpiWorld.getComm());
+  ReduceImageMPI (mpiWorld, imLocal, imGlobal);
+  if (optVerbose)
+      timerReduce.timerEndAndReport ("Time to reduce image");
+#else
+  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
+#endif
+
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0)
+#endif
+    {
+      double calcTime = timerProgram.timerEnd();
+      imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, projGlobal.remark(), projGlobal.calcTime());
+      imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime);
+      imGlobal->fileWrite (filenameImage);
+      if (optVerbose)
+	cout << "Run time: " << calcTime << " seconds" << endl;
+    }
+#ifdef HAVE_MPI
+  MPI::Finalize();
+#endif
+
+  return (0);
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+// MPI Support Routines
+//
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef HAVE_MPI
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug)
+{
+  if (mpiWorld.getRank() == 0) {
+    for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
+      for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
+	DetectorArray& detarray = projGlobal.getDetectorArray( iw );
+	int nDet = detarray.nDet();
+	DetectorValue* detval = detarray.detValues();
+
+	double viewAngle = detarray.viewAngle();
+	mpiWorld.getComm().Send(&nDet, 1, MPI::INT, iProc, 0);
+	mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, iProc, 0);
+	mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, iProc, 0);
+      }
+    }
+  }
+
+  for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
+    MPI::Status status;
+    int nDet;
+    double viewAngle;
+    DetectorValue* detval = projLocal.getDetectorArray(iw).detValues();
+
+    mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, 0, 0, status);
+    mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, 0, 0, status);
+    mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, 0, 0, status);
+    projLocal.getDetectorArray(iw).setViewAngle( viewAngle );
+  }
+}
+
+static void
+ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal)
+{
+  ImageFileArray vLocal = imLocal->getArray();
+
+  for (int ix = 0; ix < imLocal->nx(); ix++) {
+    void *recvbuf = NULL;
+    if (mpiWorld.getRank() == 0) {
+      ImageFileArray vGlobal = imGlobal->getArray();
+      recvbuf = vGlobal[ix];
+    }
+    mpiWorld.getComm().Reduce (vLocal[ix], recvbuf, imLocal->ny(), imLocal->getMPIDataType(), MPI::SUM, 0);
+  }
+}
+
+#endif
+
+
+#ifndef NO_MAIN
+int 
+main (int argc, char* argv[])
+{
+  int retval = 1;
+
+  try {
+    retval = pjrec_main(argc, argv);
+  } catch (exception e) {
+    cerr << "Exception: " << e.what() << endl;
+  } catch (...) {
+    cerr << "Unknown exception" << endl;
+  }
+
+  return (retval);
+}
+#endif
+