r184: *** empty log message ***
[ctsim.git] / tools / phm2if.cpp
index ce522b98ac13c00ffb2f880693f2fe7e4488a247..9ff3b9a1cbe54d0483d81a12287b7a2e0adf65e8 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: phm2if.cpp,v 1.2 2000/07/20 11:17:31 kevin Exp $
+**  $Id: phm2if.cpp,v 1.10 2000/08/25 15:59:13 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
@@ -50,6 +50,8 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
+static const char* g_szIdStr = "$Id: phm2if.cpp,v 1.10 2000/08/25 15:59:13 kevin Exp $";
+
 void 
 phm2if_usage (const char *program)
 {
@@ -61,9 +63,9 @@ phm2if_usage (const char *program)
   cout << "     ny              Number of pixels Y-axis" << endl;
   cout << "     --phantom       Phantom to use for projection" << endl;
   cout << "        herman       Herman head phantom" << endl;
-  cout << "        bherman      Bordered Herman head phantom" << endl;
-  cout << "        rowland      Rowland head phantom" << endl;
-  cout << "        browland     Bordered Rowland head phantom" << endl;
+  cout << "        herman-b     Herman head phantom (Bordered)" << endl;
+  cout << "        shepp-logan  Shepp-Logan head phantom" << endl;
+  cout << "        shepp-logan-b Shepp-Logan head phantom (Bordered)" << endl;
   cout << "        unitpulse    Unit pulse phantom" << endl;
   cout << "     --phmfile       Generate Phantom from phantom file" << endl;
   cout << "     --filter        Generate Phantom from a filter function" << endl;
@@ -98,43 +100,39 @@ phm2if_usage (const char *program)
 }
 
 #ifdef HAVE_MPI
-void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug);
+void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug);
 #endif
 
 int 
 phm2if_main (int argc, char* argv[])
 {
-  ImageFile* imGlobal = NULL;
+  ImageFile* pImGlobal = NULL;
   Phantom phm;
-  int opt_nx = 0, opt_ny = 0;
-  int opt_nsample = 1;
-  string optPhmName = Phantom::PHM_HERMAN_STR;
+  string optPhmName (Phantom::convertPhantomIDToName(Phantom::PHM_HERMAN));
   string optFilterName;
-  string optDomainName = SignalFilter::DOMAIN_SPATIAL_STR;
-  char *opt_outfile = NULL;
-  int opt_debug = 0;
-  string opt_desc;
-  string opt_phmFileName;
-  char *endstr, *endptr;
-  double opt_filter_param = 1;
-  double opt_filter_bw = 1.;
-  int opt_trace = TRACE_NONE;
-  bool opt_verbose = false;
-#ifdef HAVE_MPI
-  ImageFile* imLocal = NULL;
-  MPIWorld mpiWorld (argc, argv);
-#endif
-
+  string optDomainName (SignalFilter::convertDomainIDToName (SignalFilter::DOMAIN_SPATIAL));
+  string optOutFilename;
+  string optDesc;
+  string optPhmFilename;
+  int opt_nx = 0;
+  int opt_ny = 0;
+  int opt_nsample = 1;
+  double optFilterParam = 1.;
+  double optFilterBW = 1.;
+  int optTrace = TRACE_NONE;
+  bool optVerbose = false;
+  bool optDebug = false;
+  char *endptr = NULL;
+  char *endstr;
   Timer timerProgram;
-  
+
 #ifdef HAVE_MPI
+  ImageFile* pImLocal = NULL;
+  MPIWorld mpiWorld (argc, argv);
   if (mpiWorld.getRank() == 0) {
 #endif
     while (1) {
       int c = getopt_long(argc, argv, "", my_options, NULL);
-      char *endptr = NULL;
-      char *endstr;
-      
       if (c == -1)
        break;
       
@@ -143,16 +141,16 @@ phm2if_main (int argc, char* argv[])
        optPhmName = optarg;
        break;
       case O_PHMFILE:
-       opt_phmFileName = optarg;
+       optPhmFilename = optarg;
        break;
       case O_VERBOSE:
-       opt_verbose = true;
+       optVerbose = true;
        break;
       case O_DEBUG:
-       opt_debug = 1;
+       optDebug = true;
        break;
       case O_TRACE:
-       if ((opt_trace = TraceLevel::convertTraceNameToID(optarg)) == TRACE_INVALID) {
+       if ((optTrace = TraceLevel::convertTraceNameToID(optarg)) == TRACE_INVALID) {
          phm2if_usage(argv[0]);
          return (1);
        }
@@ -164,10 +162,10 @@ phm2if_main (int argc, char* argv[])
        optDomainName = optarg;
        break;
       case O_DESC:
-       opt_desc =  optarg;
+       optDesc =  optarg;
        break;
       case O_FILTER_BW:
-       opt_filter_bw = strtod(optarg, &endptr);
+       optFilterBW = strtod(optarg, &endptr);
        endstr = optarg + strlen(optarg);
        if (endptr != endstr) {
          sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
@@ -176,7 +174,7 @@ phm2if_main (int argc, char* argv[])
        }
        break;
       case O_FILTER_PARAM:
-       opt_filter_param = strtod(optarg, &endptr);
+       optFilterParam = strtod(optarg, &endptr);
        endstr = optarg + strlen(optarg);
        if (endptr != endstr) {
          sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
@@ -195,7 +193,7 @@ phm2if_main (int argc, char* argv[])
        break;
         case O_VERSION:
 #ifdef VERSION
-         cout <<  "Version " << VERSION << endl;
+         cout << "Version " << VERSION << endl << g_szIdStr << endl;
 #else
           cerr << "Unknown version number" << endl;
 #endif
@@ -209,7 +207,7 @@ phm2if_main (int argc, char* argv[])
       }
     }
     
-    if (optPhmName == "" && optFilterName == "" && opt_phmFileName == "") {
+    if (optPhmName == "" && optFilterName == "" && optPhmFilename == "") {
       cerr << "No phantom defined" << endl << endl;
       phm2if_usage(argv[0]);
       return (1);
@@ -219,7 +217,7 @@ phm2if_main (int argc, char* argv[])
       phm2if_usage(argv[0]);
       return (1);
     }
-    opt_outfile = argv[optind];
+    optOutFilename = argv[optind];
     opt_nx = strtol(argv[optind+1], &endptr, 10);
     endstr = argv[optind+1] + strlen(argv[optind+1]);
     if (endptr != endstr) {
@@ -237,16 +235,16 @@ phm2if_main (int argc, char* argv[])
     
     ostringstream oss;
     oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", ";
-    if (opt_phmFileName != "")
-      oss << "phantom=" << opt_phmFileName;
+    if (optPhmFilename != "")
+      oss << "phantomFile=" << optPhmFilename;
     else if (optPhmName != "")
       oss << "phantom=" << optPhmName;
     else if (optFilterName != "") {
       oss << "filter=" << optFilterName << " - " << optDomainName;
     }
-    if (opt_desc != "")
-      oss << ": " << opt_desc;
-    opt_desc = oss.str();
+    if (optDesc != "")
+      oss << ": " << optDesc;
+    optDesc = oss.str();
     
     if (optPhmName != "") {
       phm.createFromPhantom (optPhmName.c_str());
@@ -257,8 +255,8 @@ phm2if_main (int argc, char* argv[])
       }
     }
 
-    if (opt_phmFileName != "") {
-      phm.createFromFile(opt_phmFileName.c_str());
+    if (optPhmFilename != "") {
+      phm.createFromFile(optPhmFilename.c_str());
 #ifdef HAVE_MPI
       if (mpiWorld.getRank() == 0) 
        cerr << "Can't use phantom from file in MPI mode" << endl;
@@ -266,7 +264,7 @@ phm2if_main (int argc, char* argv[])
 #endif
     }
 
-    if (opt_verbose)
+    if (optVerbose)
       cout << "Rasterize Phantom to Image" << endl << endl;
 #ifdef HAVE_MPI
   }
@@ -275,19 +273,19 @@ phm2if_main (int argc, char* argv[])
 #ifdef HAVE_MPI
   TimerCollectiveMPI timerBcast (mpiWorld.getComm());
   mpiWorld.BcastString (optPhmName);
-  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 (&optVerbose, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
-  mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
-  mpiWorld.getComm().Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
+  mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
+  mpiWorld.getComm().Bcast (&optFilterBW, 1, MPI::DOUBLE, 0);
 
   mpiWorld.BcastString (optFilterName);
   mpiWorld.BcastString (optDomainName);
 
-  if (opt_verbose)
+  if (optVerbose)
     timerBcast.timerEndAndReport ("Time to broadcast variables");
 
   mpiWorld.setTotalWorkUnits (opt_nx);
@@ -296,17 +294,17 @@ phm2if_main (int argc, char* argv[])
       phm.createFromPhantom (optPhmName.c_str());
 
   if (mpiWorld.getRank() == 0) {
-    imGlobal = new ImageFile (opt_nx, opt_ny);
+    pImGlobal = new ImageFile (opt_nx, opt_ny);
   }
-  imLocal = new ImageFile (opt_nx, opt_ny);
+  pImLocal = new ImageFile (opt_nx, opt_ny);
 #else
-  imGlobal = new ImageFile (opt_nx, opt_ny);
+  pImGlobal = new ImageFile (opt_nx, opt_ny);
 #endif
 
   ImageFileArray v;
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0)
-    v = imGlobal->getArray ();
+    v = pImGlobal->getArray ();
 
   if (phm.getComposition() == P_UNIT_PULSE) {
     if (mpiWorld.getRank() == 0) {
@@ -314,31 +312,31 @@ phm2if_main (int argc, char* argv[])
     }
   } else if (optFilterName != "") {
     if (mpiWorld.getRank() == 0) {
-      imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param);
+      pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
     }
   } else {
     TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
-    phm.convertToImagefile (*imLocal, opt_nsample, opt_trace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits());
-    if (opt_verbose)
+    phm.convertToImagefile (*pImLocal, opt_nsample, optTrace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits());
+    if (optVerbose)
       timerRasterize.timerEndAndReport ("Time to rasterize phantom");
 
     TimerCollectiveMPI timerGather (mpiWorld.getComm());
-    mpi_gather_image (mpiWorld, imGlobal, imLocal, opt_debug);
-    if (opt_verbose)
+    mpi_gather_image (mpiWorld, pImGlobal, pImLocal, optDebug);
+    if (optVerbose)
       timerGather.timerEndAndReport ("Time to gather image");
   }
 #else
-  v = imGlobal->getArray ();
+  v = pImGlobal->getArray ();
   if (phm.getComposition() == P_UNIT_PULSE) {
     v[opt_nx/2][opt_ny/2] = 1.;
   } else if (optFilterName != "") {
-    imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param);
+    pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
   } else {
 #if HAVE_SGP
-      if (opt_trace >= TRACE_PHM)
+      if (optTrace >= TRACE_PHM)
        phm.show();
 #endif
-      phm.convertToImagefile (*imGlobal, opt_nsample, opt_trace);
+      phm.convertToImagefile (*pImGlobal, opt_nsample, optTrace);
   }
 #endif
 
@@ -347,12 +345,12 @@ phm2if_main (int argc, char* argv[])
 #endif
   {
     double calctime = timerProgram.timerEnd ();
-    imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, opt_desc.c_str(), calctime);
-    imGlobal->fileWrite (opt_outfile);
-    if (opt_verbose)
+    pImGlobal->labelAdd (Array2dFileLabel::L_HISTORY, optDesc.c_str(), calctime);
+    pImGlobal->fileWrite (optOutFilename.c_str());
+    if (optVerbose)
       cout << "Time to rasterized phantom: " << calctime << " seconds" << endl;
 
-    if (opt_trace >= TRACE_PHM) {
+    if (optTrace >= TRACE_PHM) {
       double dmin, dmax;
       int nscale;
       
@@ -360,33 +358,38 @@ phm2if_main (int argc, char* argv[])
       scanf ("%d", &nscale);
       printf ("Enter minimum and maximum densities (min, max): ");
       scanf ("%lf %lf", &dmin, &dmax);
-      imGlobal->displayScaling (nscale, dmin, dmax);
+      pImGlobal->displayScaling (nscale, dmin, dmax);
     }
   }
 
+  delete pImGlobal;
+#ifdef HAVE_MPI
+  delete pImLocal;
+#endif
+
   return (0);
 }
 
 
 
 #ifdef HAVE_MPI
-void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug)
+void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug)
 {
-  ImageFileArray vLocal = imLocal->getArray();
+  ImageFileArray vLocal = pImLocal->getArray();
   ImageFileArray vGlobal = NULL;
-  int nyLocal = imLocal->ny();
+  int nyLocal = pImLocal->ny();
 
   if (mpiWorld.getRank() == 0)
-    vGlobal = imGlobal->getArray();
+    vGlobal = pImGlobal->getArray();
   
   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
-    mpiWorld.getComm().Send(vLocal[iw], nyLocal, imLocal->getMPIDataType(), 0, 0);
+    mpiWorld.getComm().Send(vLocal[iw], nyLocal, pImLocal->getMPIDataType(), 0, 0);
 
   if (mpiWorld.getRank() == 0) {
     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
        MPI::Status status;
-       mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, imLocal->getMPIDataType(), iProc, 0, status);
+       mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, pImLocal->getMPIDataType(), iProc, 0, status);
       }
     }
   }