X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=src%2Fctrec.cpp;h=39688cfb817e74967519ceff0de6100b60c1385f;hp=19278c98bbeb2443c3de464edf6fdcb5b9bc59e4;hb=2d39e823ba389fc68e5317c422b55be006094252;hpb=a62ba7c54d39fcf1e469c8df169018fa60ae7c3b

diff --git a/src/ctrec.cpp b/src/ctrec.cpp
index 19278c9..39688cf 100644
--- a/src/ctrec.cpp
+++ b/src/ctrec.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctrec.cpp,v 1.7 2000/06/10 22:33:11 kevin Exp $
+**  $Id: ctrec.cpp,v 1.13 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
@@ -26,6 +26,8 @@
 ******************************************************************************/
 
 #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};
 
@@ -43,10 +45,11 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
+
 void 
 ctrec_usage (const char *program)
 {
-    cout << "usage: " << kbasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl;
+  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;
@@ -92,43 +95,36 @@ ctrec_usage (const char *program)
 
 
 #ifdef HAVE_MPI
-static void mpi_scatter_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int debug);
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug);
+static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal);
 #endif
 
 
 int 
 ctrec_main (int argc, char * argv[])
 {
-  ImageFile *im_global = NULL;
-  RAYSUM *rs_global = NULL;
-  char *rs_name, *im_filename = NULL;
-  char remark[MAXREMARK];
-  char filt_name[80];
-  double time_start = 0, time_end = 0;
+  ImageFile *imGlobal = NULL;
+  Projections projGlobal;
+  char *pj_name, *im_filename = NULL;
+  string remark;
   char *endptr;
   int opt_verbose = 0;
   int opt_debug = 0;
   int opt_trace = TRACE_NONE;
   double opt_filter_param = 1;
-  FilterType opt_filter = FILTER_ABS_BANDLIMIT;
-  InterpolationType opt_interp = I_LINEAR;
+  string optFilterName = "abs_bandlimit";
+  string optInterpName = "linear";
+  string optBackprojName = "idiff2";
   int opt_interp_param = 1;
-  BackprojType opt_backproj = O_BPROJ_DIFF2;
   int nx, ny;
 #ifdef HAVE_MPI
-  ImageFile *im_local;
-  RAYSUM *rs_local;
+  ImageFile *imLocal;
   int mpi_nview, mpi_ndet;
   double mpi_detinc, mpi_rotinc, mpi_phmlen;
-  double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
   MPIWorld mpiWorld (argc, argv);
 #endif
 
-#ifdef HAVE_MPI
-  time_start = MPI::Wtime();
-#else
-  time_start = td_current_sec();
-#endif
+  Timer timerProgram;
 
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0) {
@@ -143,22 +139,13 @@ ctrec_main (int argc, char * argv[])
       switch (c)
 	{
 	case O_INTERP:
-	  if ((opt_interp = opt_set_interpolation(optarg)) < 0) {
-	    ctrec_usage(argv[0]);
-	    return (1);
-	  }
+	  optInterpName = optarg;
 	  break;
 	case O_FILTER:
-	  if ((opt_filter = opt_set_filter(optarg)) < 0) {
-	    ctrec_usage(argv[0]);
-	    return (1);
-	  }
+	  optFilterName = optarg;
 	  break;
 	case O_BACKPROJ:
-	  if ((opt_backproj = opt_set_backproj(optarg)) < 0) {
-	      ctrec_usage(argv[0]);
-	      return (1);
-	  }
+	  optBackprojName = optarg;
 	  break;
 	case O_FILTER_PARAM:
 	  opt_filter_param = strtod(optarg, &endptr);
@@ -180,11 +167,11 @@ ctrec_main (int argc, char * argv[])
 	  break;
         case O_VERSION:
 #ifdef VERSION
-	cout <<  "Version " <<  VERSION << endl;
+	  cout <<  "Version " <<  VERSION << endl;
 #else
-        cout << "Unknown version number" << endl;
+          cout << "Unknown version number" << endl;
 #endif
-	  exit(0);
+	  return (0);
 	case O_HELP:
 	case '?':
 	  ctrec_usage(argv[0]);
@@ -200,51 +187,51 @@ ctrec_main (int argc, char * argv[])
       return (1);
     }
 
-    rs_name = argv[optind];
+    pj_name = argv[optind];
   
     im_filename = argv[optind + 1];
   
     nx = strtol(argv[optind + 2], &endptr, 10);
     ny = strtol(argv[optind + 3], &endptr, 10);
   
-    if (opt_filter == FILTER_G_HAMMING || opt_filter == FILTER_ABS_G_HAMMING)
-      snprintf (filt_name, sizeof(filt_name), "%s: alpha = %.2f",
-	       name_of_filter (opt_filter), opt_filter_param); 
+    ostringstream filterDesc;
+    if (opt_filter_param >= 0)
+      filterDesc << optFilterName << ": alpha=" << opt_filter_param; 
     else
-      snprintf (filt_name, sizeof(filt_name), "%s", name_of_filter (opt_filter));
-  
-    snprintf (remark, sizeof(remark), "Reconstruct: %dx%d, %s, %s, %s",
-	     nx, ny, filt_name, name_of_interpolation (opt_interp), name_of_backproj(opt_backproj));
+      filterDesc << optFilterName;
+
+    ostringstream label;
+    label << "Reconstruct: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", " << optBackprojName;
+    remark = label.str();
   
     if (opt_verbose)
-      fprintf (stdout, "%s\n", remark);
+      cout << "Remark: " << remark << endl;
 #ifdef HAVE_MPI
   }
 #endif
 
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0) {
-    rs_global = raysum_open (rs_name);
-    raysum_read (rs_global);
+    projGlobal.read (pj_name);
     if (opt_verbose)
-      raysum_print_info(rs_global);
+      projGlobal.printScanInfo();
 
-    mpi_ndet = rs_global->ndet;
-    mpi_nview = rs_global->nview;
-    mpi_detinc = rs_global->det_inc;
-    mpi_phmlen = rs_global->phmlen;
-    mpi_rotinc = rs_global->rot_inc;
+    mpi_ndet = projGlobal.nDet();
+    mpi_nview = projGlobal.nView();
+    mpi_detinc = projGlobal.detInc();
+    mpi_phmlen = projGlobal.phmLen();
+    mpi_rotinc = projGlobal.rotInc();
   }
 
-  mpi_t1 = MPI::Wtime();
+  TimerCollectiveMPI timerBcast (mpiWorld.getComm());
+  mpiWorld.BcastString (optBackprojName);
+  mpiWorld.BcastString (optFilterName);
+  mpiWorld.BcastString (optInterpName);
   mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&opt_filter, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&opt_interp, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
   mpiWorld.getComm().Bcast (&opt_interp_param, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&opt_backproj, 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);
@@ -252,149 +239,139 @@ ctrec_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 (opt_verbose) {
-    mpi_t2 = MPI::Wtime();
-    mpi_t = mpi_t2 - mpi_t1;
-    mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
-    if (mpiWorld.getRank() == 0)
-      printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
-  }
+  if (opt_verbose)
+      timerBcast.timerEndAndReport ("Time to broadcast variables");
 
   mpiWorld.setTotalWorkUnits (mpi_nview);
 
-  rs_local = raysum_create (NULL, mpiWorld.getMyLocalWorkUnits(), mpi_ndet);
-  rs_local->ndet = mpi_ndet;
-  rs_local->nview = mpi_nview;
-  rs_local->det_inc = mpi_detinc;
-  rs_local->phmlen = mpi_phmlen;
-  rs_local->rot_inc = mpi_rotinc;
+  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, opt_debug);
   if (opt_verbose)
-    mpi_t1 = MPI::Wtime();
-  mpi_scatter_rs(mpiWorld, rs_global, rs_local, opt_debug);
-  if (opt_verbose) {
-    mpi_t2 = MPI::Wtime();
-    mpi_t = mpi_t2 - mpi_t1;
-    mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
-    if (mpiWorld.getRank() == 0)
-      printf("Time to scatter rs = %f secs, Max time = %f sec\n", mpi_t, mpi_t_g);
-  }
+      timerScatter.timerEndAndReport ("Time to scatter projections");
 
   if (mpiWorld.getRank() == 0) {
-    im_global = new ImageFile (im_filename, nx, ny);
-    im_global->fileCreate();
+    imGlobal = new ImageFile (im_filename, nx, ny);
+    imGlobal->fileCreate();
   }
 
-  im_local = new ImageFile (nx, ny);
+  imLocal = new ImageFile (nx, ny);
 #else
-  rs_global = raysum_open (rs_name);
-  raysum_read (rs_global);
+  projGlobal.read (pj_name);
   if (opt_verbose)
-    raysum_print_info(rs_global);
+    projGlobal.printScanInfo();
 
-  im_global = new ImageFile (im_filename, nx, ny);
-  im_global->fileCreate();
-#endif
-
-#ifdef HAVE_MPI
-  mpi_t1 = MPI::Wtime();
-  proj_reconst (*im_local, rs_local, opt_filter, opt_filter_param, 
-		 opt_interp, opt_interp_param, opt_backproj, opt_trace);
-
-  mpi_t2 = MPI::Wtime();
-  mpi_t = mpi_t2 - mpi_t1;
-  mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
-  if (mpiWorld.getRank() == 0 && opt_verbose)
-    printf("Time to reconstruct = %f, Max time = %f\n", mpi_t, mpi_t_g);
-#else
-  proj_reconst (*im_global, rs_global, opt_filter, opt_filter_param, 
-		 opt_interp, opt_interp_param, opt_backproj, opt_trace);
+  imGlobal = new ImageFile (im_filename, nx, ny);
+  imGlobal->fileCreate();
 #endif
 
 #ifdef HAVE_MPI
+  TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
+  projLocal.reconstruct (*imLocal, optFilterName.c_str(), opt_filter_param, optInterpName.c_str(), opt_interp_param, optBackprojName.c_str(), opt_trace);
   if (opt_verbose)
-    mpi_t1 = MPI::Wtime();
-
-  int nxLocal = im_local->nx();
-  int nyLocal = im_local->ny();
-  ImageFileArray vLocal = im_local->getArray();
-  ImageFileArray vGlobal = NULL;
-  if (mpiWorld.getRank() == 0)
-    vGlobal = im_global->getArray();
-
-  for (int ix = 0; ix < nxLocal; ix++) {
-    void *recvbuf = NULL;
-    if (mpiWorld.getRank() == 0)
-      recvbuf = vGlobal[ix];
-
-    mpiWorld.getComm().Reduce(vLocal[ix], recvbuf, nyLocal, im_local->getMPIDataType(), MPI::SUM, 0);
-  }
+      timerReconstruct.timerEndAndReport ("Time to reconstruct");
 
-  if (opt_verbose) {
-    mpi_t2 = MPI::Wtime();
-    mpi_t = mpi_t2 - mpi_t1;
-    mpiWorld.getComm().Reduce (&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
-    if (mpiWorld.getRank() == 0)
-      printf("Time to reduce image = %f secs, max time = %f\n", mpi_t, mpi_t_g);
-  }
-
-  if (mpiWorld.getRank() == 0)
-      time_end = MPI::Wtime();
+  TimerCollectiveMPI timerReduce (mpiWorld.getComm());
+  ReduceImageMPI (mpiWorld, imLocal, imGlobal);
+  if (opt_verbose)
+      timerReduce.timerEndAndReport ("Time to reduce image");
 #else
-  time_end = td_current_sec();
+  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), opt_filter_param, optInterpName.c_str(), opt_interp_param, optBackprojName.c_str(), opt_trace);
 #endif
 
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0)
 #endif
     {
-      raysum_close (rs_global);
-      double calctime = time_end - time_start;
-      im_global->arrayDataWrite ();
-      im_global->labelAdd (Array2dFileLabel::L_HISTORY, rs_global->remark, rs_global->calctime);
-      im_global->labelAdd (Array2dFileLabel::L_HISTORY, remark, calctime);
-      im_global->fileClose ();
+      double calcTime = timerProgram.timerEnd();
+      imGlobal->arrayDataWrite ();
+      imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, projGlobal.remark(), projGlobal.calcTime());
+      imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime);
+      imGlobal->fileClose ();
       if (opt_verbose)
-	cout << "Time active = " << calctime << " sec" << endl;
+	cout << "Run time: " << calcTime << " seconds" << endl;
     }
-
 #ifdef HAVE_MPI
-	MPI::Finalize();
+  MPI::Finalize();
 #endif
 
-	return (0);
+  return (0);
 }
 
+
+//////////////////////////////////////////////////////////////////////////////////////
+// MPI Support Routines
+//
+//////////////////////////////////////////////////////////////////////////////////////
+
 #ifdef HAVE_MPI
-static void mpi_scatter_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int opt_debug)
 {
   if (mpiWorld.getRank() == 0) {
     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
-	mpiWorld.getComm().Send(&rs_global->view[iw]->ndet, 1, MPI::INT, iProc, 0);
-	mpiWorld.getComm().Send(&rs_global->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0);
-	mpiWorld.getComm().Send(rs_global->view[iw]->detval, rs_global->ndet, MPI::FLOAT, iProc, 0);
+	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 );
+  }
+}
 
-    mpiWorld.getComm().Recv(&rs_local->view[iw]->ndet, 1, MPI::INT, 0, 0, status);
-    mpiWorld.getComm().Recv(&rs_local->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0, status);
-    mpiWorld.getComm().Recv(rs_local->view[iw]->detval, rs_local->ndet, MPI::FLOAT, 0, 0, status);
+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);
   }
-  rs_local->nview = mpiWorld.getMyLocalWorkUnits();
 }
 
 #endif
 
+
 #ifndef NO_MAIN
 int 
 main (int argc, char* argv[])
 {
-  return (ctrec_main(argc, argv));
+  int retval = 1;
+
+  try {
+    retval = ctrec_main(argc, argv);
+  } catch (exception e) {
+    cerr << "Exception: " << e.what() << endl;
+  } catch (...) {
+    cerr << "Unknown exception" << endl;
+  }
+
+  return (retval);
 }
 #endif