X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=src%2Fctrec.cpp;h=ada1c10ce3c49c3c7a6a943f2c4c2762b0352079;hp=f48fc584de6d60b110b1835060df5ebaf71d44d6;hb=4b16507e92bb80b09575b28bed66810e33d2681f;hpb=f173363bba9997045e5ec825e64d6253ec4da235

diff --git a/src/ctrec.cpp b/src/ctrec.cpp
index f48fc58..ada1c10 100644
--- a/src/ctrec.cpp
+++ b/src/ctrec.cpp
@@ -1,53 +1,15 @@
 /*****************************************************************************
-**  This is part of the CTSim program
-**  Copyright (C) 1983-2000 Kevin Rosenberg
-**
-**  $Id: ctrec.cpp,v 1.1 2000/06/07 02:29:05 kevin Exp $
-**  $Log: ctrec.cpp,v $
-**  Revision 1.1  2000/06/07 02:29:05  kevin
-**  Initial C++ versions
-**
-**  Revision 1.13  2000/06/05 01:33:25  kevin
-**  BSpline changes
-**
-**  Revision 1.12  2000/05/24 22:50:04  kevin
-**  Added support for new SGP library
-**
-**  Revision 1.11  2000/05/16 04:33:59  kevin
-**  Improved option processing
-**
-**  Revision 1.10  2000/05/11 01:06:30  kevin
-**  Changed sprintf to snprintf
-**
-**  Revision 1.9  2000/05/08 20:02:32  kevin
-**  ANSI C changes
-**
-**  Revision 1.8  2000/05/04 18:16:34  kevin
-**  renamed filter definitions
-**
-**  Revision 1.7  2000/05/03 08:49:50  kevin
-**  Code cleanup
-**
-**  Revision 1.6  2000/05/02 15:31:47  kevin
-**  code cleaning
-**
-**  Revision 1.5  2000/04/30 11:41:06  kevin
-**  Cleaned up debugging code
-**
-**  Revision 1.4  2000/04/30 10:13:27  kevin
-**  Fixed MPI bugs
-**
-**  Revision 1.3  2000/04/30 04:06:13  kevin
-**  Update Raysum i/o routines
-**  Fix MPI bug in ctrec (scatter_raysum) that referenced rs_global
+** FILE IDENTIFICATION
 **
-**  Revision 1.2  2000/04/29 23:24:56  kevin
-**  *** empty log message ***
-**
-**  Revision 1.1.1.1  2000/04/28 13:02:44  kevin
-**  Initial CVS import for first public release
+**   Name:          ctrec.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: ctrec.cpp,v 1.6 2000/06/09 11:03:08 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
@@ -63,15 +25,6 @@
 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 ******************************************************************************/
 
-/* FILE
- *   ctrec.c		Reconstruct an image from raysums
- *
- * DATE
- *   Aug 84
- *   Jul 99  -- Converted to ANSI C
- *              Added MPI parallel processing
- */
-
 #include "ct.h"
 
 enum {O_INTERP, O_FILTER, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
@@ -138,21 +91,20 @@ ctrec_usage (const char *program)
 }
 
 
-#ifdef MPI_CT
-static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int debug);
+#ifdef HAVE_MPI
+static void mpi_scatter_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int debug);
 #endif
 
 static void print_raysum_info(const RAYSUM *rs);
 
 int 
-ctrec_main (const int argc, char *const argv[])
+ctrec_main (int argc, char * argv[])
 {
-  IMAGE *im_global = NULL;
+  ImageFile *im_global = NULL;
   RAYSUM *rs_global = NULL;
   char *rs_name, *im_filename = NULL;
   char remark[MAXREMARK];
   char filt_name[80];
-  int nx, ny;
   double time_start = 0, time_end = 0;
   char *endptr;
   int opt_verbose = 0;
@@ -163,36 +115,24 @@ ctrec_main (const int argc, char *const argv[])
   InterpolationType opt_interp = I_LINEAR;
   int opt_interp_param = 1;
   BackprojType opt_backproj = O_BPROJ_DIFF2;
-#ifdef MPI_CT
-  IMAGE *im_local;
+  int nx, ny;
+#ifdef HAVE_MPI
+  ImageFile *im_local;
   RAYSUM *rs_local;
-  int mpi_argc = argc;
-  char **mpi_argv = (char **) argv;
   int mpi_nview, mpi_ndet;
   double mpi_detinc, mpi_rotinc, mpi_phmlen;
   double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
-  int ix;
-
-  MPI_Init(&mpi_argc, &mpi_argv);
-  MPI_Comm_dup (MPI_COMM_WORLD, &mpi_ct.comm);
-  MPI_Comm_size(mpi_ct.comm, &mpi_ct.nproc);
-  MPI_Comm_rank(mpi_ct.comm, &mpi_ct.my_rank);
-
-  if (mpi_ct.nproc > MPI_MAX_PROCESS) {
-    sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
-	      mpi_ct.nproc, MPI_MAX_PROCESS);
-    exit(1);
-  }
+  MPIWorld mpiWorld (argc, argv);
 #endif
 
-#ifdef MPI_CT
-  time_start = MPI_Wtime();
+#ifdef HAVE_MPI
+  time_start = MPI::Wtime();
 #else
   time_start = td_current_sec();
 #endif
 
-#ifdef MPI_CT
-  if (mpi_ct.my_rank == 0) {
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0) {
 #endif
     while (1) {
       int c = getopt_long(argc, argv, "", my_options, NULL);
@@ -279,12 +219,12 @@ ctrec_main (const int argc, char *const argv[])
   
     if (opt_verbose)
       fprintf (stdout, "%s\n", remark);
-#ifdef MPI_CT
+#ifdef HAVE_MPI
   }
 #endif
 
-#ifdef MPI_CT
-  if (mpi_ct.my_rank == 0) {
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0) {
     rs_global = raysum_open (rs_name);
     raysum_read (rs_global);
     if (opt_verbose)
@@ -297,33 +237,33 @@ ctrec_main (const int argc, char *const argv[])
     mpi_rotinc = rs_global->rot_inc;
   }
 
-  mpi_t1 = MPI_Wtime();
-  MPI_Bcast(&opt_verbose, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_debug, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_trace, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_filter, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_interp, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_filter_param, 1, MPI_DOUBLE, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_interp_param, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&opt_backproj, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&mpi_ndet, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&mpi_nview, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&mpi_detinc, 1, MPI_DOUBLE, 0, mpi_ct.comm);
-  MPI_Bcast(&mpi_phmlen, 1, MPI_DOUBLE, 0, mpi_ct.comm);
-  MPI_Bcast(&mpi_rotinc, 1, MPI_DOUBLE, 0, mpi_ct.comm);
-  MPI_Bcast(&nx, 1, MPI_INT, 0, mpi_ct.comm);
-  MPI_Bcast(&ny, 1, MPI_INT, 0, mpi_ct.comm);
+  mpi_t1 = MPI::Wtime();
+  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);
+  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 (opt_verbose) {
-    mpi_t2 = MPI_Wtime();
+    mpi_t2 = MPI::Wtime();
     mpi_t = mpi_t2 - mpi_t1;
-    MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
-    if (mpi_ct.my_rank == 0)
+    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);
   }
 
-  mpi_ct_calc_work_units(mpi_nview);
+  mpiWorld.setTotalWorkUnits (mpi_nview);
 
-  rs_local = raysum_create (NULL, mpi_ct.local_work_units[mpi_ct.my_rank], mpi_ndet);
+  rs_local = raysum_create (NULL, mpiWorld.getMyLocalWorkUnits(), mpi_ndet);
 
   rs_local->ndet = mpi_ndet;
   rs_local->nview = mpi_nview;
@@ -332,130 +272,124 @@ ctrec_main (const int argc, char *const argv[])
   rs_local->rot_inc = mpi_rotinc;
 
   if (opt_verbose)
-    mpi_t1 = MPI_Wtime();
-  mpi_scatter_rs(rs_global, rs_local, opt_debug);
+    mpi_t1 = MPI::Wtime();
+  mpi_scatter_rs(mpiWorld, rs_global, rs_local, opt_debug);
   if (opt_verbose) {
-    mpi_t2 = MPI_Wtime();
+    mpi_t2 = MPI::Wtime();
     mpi_t = mpi_t2 - mpi_t1;
-    MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
-    if (mpi_ct.my_rank == 0)
+    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);
   }
 
-  if (mpi_ct.my_rank == 0) {
-    im_global = image_create (im_filename, nx, ny);
-    sdf_add_label (LT_HISTORY, rs_global->remark, rs_global->calctime, im_global->dfp_2d->dfp);
-    sdf_add_empty_label (im_global->dfp_2d->dfp);
+  if (mpiWorld.getRank() == 0) {
+    im_global = new ImageFile (im_filename, nx, ny);
+    im_global->fileCreate();
   }
-  im_local = image_create (NULL, nx, ny);
 
+  im_local = new ImageFile (nx, ny);
 #else
   rs_global = raysum_open (rs_name);
   raysum_read (rs_global);
   if (opt_verbose)
     print_raysum_info(rs_global);
 
-  im_global = image_create (im_filename, nx, ny);
-  sdf_add_label (LT_HISTORY, rs_global->remark, rs_global->calctime, im_global->dfp_2d->dfp);
-  sdf_add_empty_label (im_global->dfp_2d->dfp);
+  im_global = new ImageFile (im_filename, nx, ny);
+  im_global->fileCreate();
 #endif
 
-#ifdef MPI_CT
-  mpi_t1 = MPI_Wtime();
-  image_reconst (im_local, rs_local, opt_filter, opt_filter_param, 
+#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_t2 = MPI::Wtime();
   mpi_t = mpi_t2 - mpi_t1;
-  MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
-  if (mpi_ct.my_rank == 0 && opt_verbose)
+  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
-  image_reconst (im_global, rs_global, opt_filter, opt_filter_param, 
+  proj_reconst (*im_global, rs_global, opt_filter, opt_filter_param, 
 		 opt_interp, opt_interp_param, opt_backproj, opt_trace);
 #endif
 
-#ifdef MPI_CT
+#ifdef HAVE_MPI
   if (opt_verbose)
-    mpi_t1 = MPI_Wtime();
+    mpi_t1 = MPI::Wtime();
 
-  for (ix = 0; ix < im_local->nx; ix++) {
+  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 (mpi_ct.my_rank == 0)
-      recvbuf = im_global->v[ix];
+    if (mpiWorld.getRank() == 0)
+      recvbuf = vGlobal[ix];
 
-    MPI_Reduce(im_local->v[ix], recvbuf, im_local->ny, MPI_FLOAT, MPI_SUM, 0, mpi_ct.comm);
+    mpiWorld.getComm().Reduce(vLocal[ix], recvbuf, nyLocal, im_local->getMPIDataType(), MPI::SUM, 0);
   }
 
   if (opt_verbose) {
-    mpi_t2 = MPI_Wtime();
+    mpi_t2 = MPI::Wtime();
     mpi_t = mpi_t2 - mpi_t1;
-    MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
-    if (mpi_ct.my_rank == 0)
+    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 (mpi_ct.my_rank == 0)
-      time_end = MPI_Wtime();
+  if (mpiWorld.getRank() == 0)
+      time_end = MPI::Wtime();
 #else
   time_end = td_current_sec();
 #endif
 
 
-#ifdef MPI_CT
-  if (mpi_ct.my_rank == 0)
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0)
 #endif
     {
       raysum_close (rs_global);
-      strncpy (im_global->remark, remark, sizeof(im_global->remark));
-      im_global->calctime = time_end - time_start;
-      image_save (im_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 ();
       if (opt_verbose)
-	  fprintf (stdout, "Time active = %.2f\n", im_global->calctime);
+	cout << "Time active = " << calctime << " sec" << endl;
     }
 
-#ifdef MPI_CT
-	MPI_Finalize();
+#ifdef HAVE_MPI
+	MPI::Finalize();
 #endif
 
 	return (0);
 }
 
 
-#ifdef MPI_CT
-static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
+#ifdef HAVE_MPI
+static void mpi_scatter_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
 {
-  int iproc;
-  int end_work_unit;
-  int iw = 0;
-
-  if (mpi_ct.my_rank == 0) {
-    for (iproc = 0; iproc < mpi_ct.nproc; iproc++) {
-      end_work_unit = mpi_ct.start_work_unit[iproc] + mpi_ct.local_work_units[iproc]  - 1;
-
-      for (iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
-	MPI_Send(&rs_global->view[iw]->ndet, 1, MPI_INT, iproc, 0, mpi_ct.comm);
-	MPI_Send(&rs_global->view[iw]->view_angle, 1, MPI_DOUBLE, iproc, 0, mpi_ct.comm);
-	MPI_Send(rs_global->view[iw]->detval, rs_global->ndet, MPI_FLOAT, iproc, 0, mpi_ct.comm);
+  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);
       }
     }
   }
 
-  if (opt_debug)
-    fprintf(stdout, "Receiving rs data in process %d\n", mpi_ct.my_rank);
-
-  end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
-  for (iw = 0; iw <= end_work_unit; iw++) {
-    MPI_Status status;
+  for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
+    MPI::Status status;
 
-    MPI_Recv(&rs_local->view[iw]->ndet, 1, MPI_INT, 0, 0, mpi_ct.comm, &status);
-    MPI_Recv(&rs_local->view[iw]->view_angle, 1, MPI_DOUBLE, 0, 0, mpi_ct.comm, &status);
-    MPI_Recv(rs_local->view[iw]->detval, rs_local->ndet, MPI_FLOAT, 0, 0, mpi_ct.comm, &status);
-  }
-  rs_local->nview = mpi_ct.local_work_units[mpi_ct.my_rank];
-  if (opt_debug) {
-    MPI_Barrier(MPI_COMM_WORLD);
-    fprintf(stdout, "Done with mpi_scatter_rs in process %2d\n", mpi_ct.my_rank);
+    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);
   }
+  rs_local->nview = mpiWorld.getMyLocalWorkUnits();
 }
 
 #endif
@@ -474,7 +408,7 @@ static void print_raysum_info(const RAYSUM *rs)
 
 #ifndef NO_MAIN
 int 
-main (const int argc, char *const argv[])
+main (int argc, char* argv[])
 {
   return (ctrec_main(argc, argv));
 }