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

diff --git a/src/ctrec.cpp b/src/ctrec.cpp
index f48fc58..fed86ee 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
+** FILE IDENTIFICATION
 **
-**  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
-**
-**  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.4 2000/06/08 16:43:10 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,7 +91,7 @@ ctrec_usage (const char *program)
 }
 
 
-#ifdef MPI_CT
+#ifdef HAVE_MPI
 static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int debug);
 #endif
 
@@ -147,12 +100,11 @@ static void print_raysum_info(const RAYSUM *rs);
 int 
 ctrec_main (const int argc, char *const 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,35 +115,35 @@ 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);
+  MPI::Init (mpi_argc, mpi_argv);
+  mpi_ct.comm = MPI::COMM_WORLD.Dup();
+  mpi_ct.nproc = mpi_ct.comm.Get_size();
+  mpi_ct.my_rank = mpi_ct.comm.Get_rank();
 
-  if (mpi_ct.nproc > MPI_MAX_PROCESS) {
+  if (mpi_ct.nproc > CT_MPI_MAX_PROCESS) {
     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
-	      mpi_ct.nproc, MPI_MAX_PROCESS);
+	      mpi_ct.nproc, CT_MPI_MAX_PROCESS);
     exit(1);
   }
 #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
+#ifdef HAVE_MPI
   if (mpi_ct.my_rank == 0) {
 #endif
     while (1) {
@@ -279,11 +231,11 @@ 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
+#ifdef HAVE_MPI
   if (mpi_ct.my_rank == 0) {
     rs_global = raysum_open (rs_name);
     raysum_read (rs_global);
@@ -297,26 +249,26 @@ 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();
+  mpi_ct.comm.Bcast (&opt_verbose, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&opt_debug, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&opt_trace, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&opt_filter, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&opt_interp, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
+  mpi_ct.comm.Bcast (&opt_interp_param, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&opt_backproj, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&mpi_ndet, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&mpi_nview, 1, MPI::INT, 0);
+  mpi_ct.comm.Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
+  mpi_ct.comm.Bcast (&mpi_phmlen, 1, MPI::DOUBLE, 0);
+  mpi_ct.comm.Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0);
+  mpi_ct.comm.Bcast (&nx, 1, MPI::INT, 0);
+  mpi_ct.comm.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);
+    mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
     if (mpi_ct.my_rank == 0)
       printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
   }
@@ -332,96 +284,104 @@ ctrec_main (const int argc, char *const argv[])
   rs_local->rot_inc = mpi_rotinc;
 
   if (opt_verbose)
-    mpi_t1 = MPI_Wtime();
+    mpi_t1 = MPI::Wtime();
   mpi_scatter_rs(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);
+    mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
     if (mpi_ct.my_rank == 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);
+    im_global = new ImageFile (im_filename, nx, ny);
+    im_global->adf.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->adf.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);
+  mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
   if (mpi_ct.my_rank == 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();
+
+  int nxLocal = im_local->adf.nx();
+  int nyLocal = im_local->adf.ny();
+  ImageFileArray vLocal = im_local->getArray();
+  ImageFileArray vGlobal = NULL;
+  if (mpi_ct.my_rank == 0)
+    vGlobal = im_global->getArray();
 
-  for (ix = 0; ix < im_local->nx; ix++) {
+  for (int ix = 0; ix < nxLocal; ix++) {
     void *recvbuf = NULL;
     if (mpi_ct.my_rank == 0)
-      recvbuf = im_global->v[ix];
+      recvbuf = vGlobal[ix];
 
-    MPI_Reduce(im_local->v[ix], recvbuf, im_local->ny, MPI_FLOAT, MPI_SUM, 0, mpi_ct.comm);
+    mpi_ct.comm.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);
+    mpi_ct.comm.Reduce (&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
     if (mpi_ct.my_rank == 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();
+      time_end = MPI::Wtime();
 #else
   time_end = td_current_sec();
 #endif
 
 
-#ifdef MPI_CT
+#ifdef HAVE_MPI
   if (mpi_ct.my_rank == 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->adf.arrayDataWrite ();
+      im_global->adf.labelAdd (Array2dFileLabel::L_HISTORY, rs_global->remark, rs_global->calctime);
+      im_global->adf.labelAdd (Array2dFileLabel::L_HISTORY, remark, calctime);
+      im_global->adf.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
+#ifdef HAVE_MPI
 static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
 {
   int iproc;
@@ -433,9 +393,9 @@ static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_d
       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);
+	mpi_ct.comm.Send(&rs_global->view[iw]->ndet, 1, MPI::INT, iproc, 0);
+	mpi_ct.comm.Send(&rs_global->view[iw]->view_angle, 1, MPI::DOUBLE, iproc, 0);
+	mpi_ct.comm.Send(rs_global->view[iw]->detval, rs_global->ndet, MPI::FLOAT, iproc, 0);
       }
     }
   }
@@ -445,15 +405,15 @@ static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_d
 
   end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
   for (iw = 0; iw <= end_work_unit; iw++) {
-    MPI_Status status;
+    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);
+    mpi_ct.comm.Recv(&rs_local->view[iw]->ndet, 1, MPI::INT, 0, 0, status);
+    mpi_ct.comm.Recv(&rs_local->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0, status);
+    mpi_ct.comm.Recv(rs_local->view[iw]->detval, rs_local->ndet, MPI::FLOAT, 0, 0, status);
   }
   rs_local->nview = mpi_ct.local_work_units[mpi_ct.my_rank];
   if (opt_debug) {
-    MPI_Barrier(MPI_COMM_WORLD);
+    mpi_ct.comm.Barrier();
     fprintf(stdout, "Done with mpi_scatter_rs in process %2d\n", mpi_ct.my_rank);
   }
 }