Update copyright date; remove old CVS keyword
[ctsim.git] / tools / phm2if.cpp
index f56a8f1a13c102d9d409b2930f047bc41fae8906..ce08c25f0ac236c422d3ca0314dd2fe5c0a8c3a4 100644 (file)
@@ -7,9 +7,7 @@
 **   Date Started:  1984
 **
 **  This is part of the CTSim program
-**  Copyright (C) 1983-2000 Kevin Rosenberg
-**
-**  $Id: phm2if.cpp,v 1.13 2000/11/09 00:12:25 kevin Exp $
+**  Copyright (C) 1983-2009 Kevin Rosenberg
 **
 **  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
 #include "timer.h"
 
 
-enum { O_PHANTOM, O_DESC, O_NSAMPLE, O_FILTER, O_TRACE, O_VERBOSE, O_HELP, 
-       O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
+enum { O_PHANTOM, O_DESC, O_NSAMPLE, O_FILTER, O_VIEW_RATIO, O_TRACE, O_VERBOSE, O_HELP,
+O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
 
-static struct option my_options[] = 
+static struct option my_options[] =
 {
   {"phantom", 1, 0, O_PHANTOM},
   {"phmfile", 1, 0, O_PHMFILE},
@@ -43,6 +41,7 @@ static struct option my_options[] =
   {"filter-bw", 1, 0, O_FILTER_BW},
   {"filter-param", 1, 0, O_FILTER_PARAM},
   {"trace", 1, 0, O_TRACE},
+  {"view-ratio", 1, 0, O_VIEW_RATIO},
   {"verbose", 0, 0, O_VERBOSE},
   {"debug", 0, 0, O_DEBUG},
   {"help", 0, 0, O_HELP},
@@ -50,73 +49,69 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: phm2if.cpp,v 1.13 2000/11/09 00:12:25 kevin Exp $";
+static const char* g_szIdStr = "$Id$";
 
-void 
+void
 phm2if_usage (const char *program)
 {
-  cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]" << endl;
-  cout << "Generate phantom image from a predefined --phantom or a --phmfile" << endl;
-  cout << endl;
-  cout << "     outfile         Name of output file for image" << endl;
-  cout << "     nx              Number of pixels X-axis" << endl;
-  cout << "     ny              Number of pixels Y-axis" << endl;
-  cout << "     --phantom       Phantom to use for projection" << endl;
-  cout << "        herman       Herman 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;
-  cout << "       abs_bandlimit Abs * Bandlimiting" << 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 << "     --filter-param  Alpha level for Hamming filter" << endl;
-  cout << "     --filter-domain Set domain of filter" << endl;
-  cout << "         spatial     Spatial domain (default)" << endl;
-  cout << "         freq        Frequency domain" << endl;
-  cout << "     --filter-bw     Filter bandwidth (default = 1)" << endl;
-  cout << "     --desc          Description of raysum" << endl;
-  cout << "     --nsample       Number of samples per axis per pixel (default = 1)" << endl;
-  cout << "     --trace         Trace level to use" << endl;
-  cout << "        none         No tracing (default)" << endl;
-  cout << "        console      Trace text level" << endl;
-  cout << "        phantom      Trace phantom" << endl;
-  cout << "        proj         Trace projections" << endl;
-  cout << "        plot         Trace plot" << endl;
-  cout << "        clipping     Trace clipping" << endl;
-  cout << "     --debug         Debug mode" << endl;
-  cout << "     --verbose       Verbose mode" << endl;
-  cout << "     --version       Print version" << endl;
-  cout << "     --help          Print this help message" << endl;
+  std::cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]\n";
+  std::cout << "Generate phantom image from a predefined --phantom or a --phmfile\n";
+  std::cout << std::endl;
+  std::cout << "     outfile         Name of output file for image\n";
+  std::cout << "     nx              Number of pixels X-axis\n";
+  std::cout << "     ny              Number of pixels Y-axis\n";
+  std::cout << "     --view-ratio    View diameter to phantom diameter ratio (default = 1)\n";
+  std::cout << "     --phantom       Phantom to use for projection\n";
+  std::cout << "        herman       Herman head phantom\n";
+  std::cout << "        shepp-logan  Shepp-Logan head phantom\n";
+  std::cout << "        unit-pulse    Unit pulse phantom\n";
+  std::cout << "     --phmfile       Generate Phantom from phantom file\n";
+  std::cout << "     --filter        Generate Phantom from a filter function\n";
+  std::cout << "       abs_bandlimit Abs * Bandlimiting\n";
+  std::cout << "       abs_sinc      Abs * Sinc\n";
+  std::cout << "       abs_cos       Abs * Cosine\n";
+  std::cout << "       abs_hamming   Abs * Hamming\n";
+  std::cout << "       shepp         Shepp-Logan\n";
+  std::cout << "       bandlimit     Bandlimiting\n";
+  std::cout << "       sinc          Sinc\n";
+  std::cout << "       cos           Cosine\n";
+  std::cout << "       triangle      Triangle\n";
+  std::cout << "       hamming       Hamming\n";
+  std::cout << "     --filter-param  Alpha level for Hamming filter\n";
+  std::cout << "     --filter-domain Set domain of filter\n";
+  std::cout << "         spatial     Spatial domain (default)\n";
+  std::cout << "         freq        Frequency domain\n";
+  std::cout << "     --filter-bw     Filter bandwidth (default = 1)\n";
+  std::cout << "     --desc          Description of raysum\n";
+  std::cout << "     --nsample       Number of samples per axis per pixel (default = 1)\n";
+  std::cout << "     --trace         Trace level to use\n";
+  std::cout << "        none         No tracing (default)\n";
+  std::cout << "        console      Trace text level\n";
+  std::cout << "     --debug         Debug mode\n";
+  std::cout << "     --verbose       Verbose mode\n";
+  std::cout << "     --version       Print version\n";
+  std::cout << "     --help          Print this help message\n";
 }
 
 #ifdef HAVE_MPI
 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug);
 #endif
 
-int 
-phm2if_main (int argc, char* argv[])
+int
+phm2if_main (int argc, char* const argv[])
 {
   ImageFile* pImGlobal = NULL;
   Phantom phm;
-  string optPhmName;
-  string optFilterName;
-  string optDomainName (SignalFilter::convertDomainIDToName (SignalFilter::DOMAIN_SPATIAL));
-  string optOutFilename;
-  string optDesc;
-  string optPhmFilename;
+  std::string optPhmName;
+  std::string optFilterName;
+  std::string optDomainName (SignalFilter::convertDomainIDToName (SignalFilter::DOMAIN_SPATIAL));
+  std::string optOutFilename;
+  std::string optDesc;
+  std::string optPhmFilename;
   int opt_nx = 0;
   int opt_ny = 0;
   int opt_nsample = 1;
+  double optViewRatio = 1.;
   double optFilterParam = 1.;
   double optFilterBW = 1.;
   int optTrace = Trace::TRACE_NONE;
@@ -134,81 +129,90 @@ phm2if_main (int argc, char* argv[])
     while (1) {
       int c = getopt_long(argc, argv, "", my_options, NULL);
       if (c == -1)
-       break;
-      
+        break;
+
       switch (c) {
       case O_PHANTOM:
-       optPhmName = optarg;
-       break;
+        optPhmName = optarg;
+        break;
       case O_PHMFILE:
-       optPhmFilename = optarg;
-       break;
+        optPhmFilename = optarg;
+        break;
       case O_VERBOSE:
-       optVerbose = true;
-       break;
+        optVerbose = true;
+        break;
       case O_DEBUG:
-       optDebug = true;
-       break;
+        optDebug = true;
+        break;
       case O_TRACE:
-       if ((optTrace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
-         phm2if_usage(argv[0]);
-         return (1);
-       }
-       break;
+        if ((optTrace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
+          phm2if_usage(argv[0]);
+          return (1);
+        }
+        break;
       case O_FILTER:
-       optFilterName = optarg;
-       break;
+        optFilterName = optarg;
+        break;
       case O_FILTER_DOMAIN:
-       optDomainName = optarg;
-       break;
+        optDomainName = optarg;
+        break;
       case O_DESC:
-       optDesc =  optarg;
-       break;
+        optDesc =  optarg;
+        break;
       case O_FILTER_BW:
-       optFilterBW = strtod(optarg, &endptr);
-       endstr = optarg + strlen(optarg);
-       if (endptr != endstr) {
-         sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
-         phm2if_usage(argv[0]);
-         return (1);
-       }
-       break;
+        optFilterBW = strtod(optarg, &endptr);
+        endstr = optarg + strlen(optarg);
+        if (endptr != endstr) {
+          sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
+          phm2if_usage(argv[0]);
+          return (1);
+        }
+        break;
+      case O_VIEW_RATIO:
+        optViewRatio = strtod(optarg, &endptr);
+        endstr = optarg + strlen(optarg);
+        if (endptr != endstr) {
+          sys_error(ERR_SEVERE,"Error setting --view-ratio to %s\n", optarg);
+          phm2if_usage(argv[0]);
+          return (1);
+        }
+        break;
       case O_FILTER_PARAM:
-       optFilterParam = strtod(optarg, &endptr);
-       endstr = optarg + strlen(optarg);
-       if (endptr != endstr) {
-         sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
-         phm2if_usage(argv[0]);
-         return (1);
-       }
-       break;
+        optFilterParam = strtod(optarg, &endptr);
+        endstr = optarg + strlen(optarg);
+        if (endptr != endstr) {
+          sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
+          phm2if_usage(argv[0]);
+          return (1);
+        }
+        break;
       case O_NSAMPLE:
-       opt_nsample = strtol(optarg, &endptr, 10);
-       endstr = optarg + strlen(optarg);
-       if (endptr != endstr) {
-         sys_error(ERR_SEVERE,"Error setting --nsample to %s\n", optarg);
-         phm2if_usage(argv[0]);
-         return (1);
-       }
-       break;
-        case O_VERSION:
+        opt_nsample = strtol(optarg, &endptr, 10);
+        endstr = optarg + strlen(optarg);
+        if (endptr != endstr) {
+          sys_error(ERR_SEVERE,"Error setting --nsample to %s\n", optarg);
+          phm2if_usage(argv[0]);
+          return (1);
+        }
+        break;
+      case O_VERSION:
 #ifdef VERSION
-         cout << "Version " << VERSION << endl << g_szIdStr << endl;
+        std::cout << "Version " << VERSION << std::endl << g_szIdStr << std::endl;
 #else
-          cerr << "Unknown version number" << endl;
+        std::cerr << "Unknown version number\n";
 #endif
       case O_HELP:
       case '?':
-       phm2if_usage(argv[0]);
-       return (0);
+        phm2if_usage(argv[0]);
+        return (0);
       default:
-       phm2if_usage(argv[0]);
-       return (1);
+        phm2if_usage(argv[0]);
+        return (1);
       }
     }
-    
+
     if (optPhmName == "" && optFilterName == "" && optPhmFilename == "") {
-      cerr << "No phantom defined" << endl << endl;
+      std::cerr << "No phantom defined\n" << std::endl;
       phm2if_usage(argv[0]);
       return (1);
     }
@@ -232,9 +236,9 @@ phm2if_main (int argc, char* argv[])
       phm2if_usage(argv[0]);
       return (1);
     }
-    
-    ostringstream oss;
-    oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", ";
+
+    std::ostringstream oss;
+    oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", viewRatio=" << optViewRatio << ", nsample=" << opt_nsample << ", ";
     if (optPhmFilename != "")
       oss << "phantomFile=" << optPhmFilename;
     else if (optPhmName != "")
@@ -245,27 +249,27 @@ phm2if_main (int argc, char* argv[])
     if (optDesc != "")
       oss << ": " << optDesc;
     optDesc = oss.str();
-    
+
     if (optPhmName != "") {
       phm.createFromPhantom (optPhmName.c_str());
       if (phm.fail()) {
-       cout << phm.failMessage() << endl << endl;
-       phm2if_usage(argv[0]);
-       return (1);
+        std::cout << phm.failMessage() << std::endl << std::endl;
+        phm2if_usage(argv[0]);
+        return (1);
       }
     }
 
     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;
+      if (mpiWorld.getRank() == 0)
+        std::cerr << "Can't use phantom from file in MPI mode\n";
       return (1);
 #endif
     }
 
     if (optVerbose)
-      cout << "Rasterize Phantom to Image" << endl << endl;
+      std::cout << "Rasterize Phantom to Image\n" << std::endl;
 #ifdef HAVE_MPI
   }
 #endif
@@ -279,6 +283,7 @@ phm2if_main (int argc, char* argv[])
   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 (&optViewRatio, 1, MPI::DOUBLE, 0);
   mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
   mpiWorld.getComm().Bcast (&optFilterBW, 1, MPI::DOUBLE, 0);
 
@@ -291,7 +296,7 @@ phm2if_main (int argc, char* argv[])
   mpiWorld.setTotalWorkUnits (opt_nx);
 
   if (mpiWorld.getRank() > 0 && optPhmName != "")
-      phm.createFromPhantom (optPhmName.c_str());
+    phm.createFromPhantom (optPhmName.c_str());
 
   if (mpiWorld.getRank() == 0) {
     pImGlobal = new ImageFile (opt_nx, opt_ny);
@@ -301,7 +306,7 @@ phm2if_main (int argc, char* argv[])
   pImGlobal = new ImageFile (opt_nx, opt_ny);
 #endif
 
-  ImageFileArray v;
+  ImageFileArray v = NULL;
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0)
     v = pImGlobal->getArray ();
@@ -316,7 +321,7 @@ phm2if_main (int argc, char* argv[])
     }
   } else {
     TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
-    phm.convertToImagefile (*pImLocal, opt_nsample, optTrace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), false);
+    phm.convertToImagefile (*pImLocal, optViewRatio, opt_nsample, optTrace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), false);
     if (optVerbose)
       timerRasterize.timerEndAndReport ("Time to rasterize phantom");
 
@@ -332,34 +337,19 @@ phm2if_main (int argc, char* argv[])
   } else if (optFilterName != "") {
     pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
   } else {
-#if HAVE_SGP
-      if (optTrace >= Trace::TRACE_PHANTOM)
-       phm.show();
-#endif
-      phm.convertToImagefile (*pImGlobal, opt_nsample, optTrace);
+    phm.convertToImagefile (*pImGlobal, optViewRatio, opt_nsample, optTrace);
   }
 #endif
 
 #ifdef HAVE_MPI
-  if (mpiWorld.getRank() == 0) 
+  if (mpiWorld.getRank() == 0)
 #endif
   {
     double calctime = timerProgram.timerEnd ();
     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 (optTrace >= Trace::TRACE_PHANTOM) {
-      double dmin, dmax;
-      int nscale;
-      
-      printf ("Enter display size scale (nominal = 1): ");
-      scanf ("%d", &nscale);
-      printf ("Enter minimum and maximum densities (min, max): ");
-      scanf ("%lf %lf", &dmin, &dmax);
-      pImGlobal->displayScaling (nscale, dmin, dmax);
-    }
+      std::cout << "Time to rasterize phantom: " << calctime << " seconds\n";
   }
 
   delete pImGlobal;
@@ -381,15 +371,15 @@ void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImL
 
   if (mpiWorld.getRank() == 0)
     vGlobal = pImGlobal->getArray();
-  
+
   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
     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, pImLocal->getMPIDataType(), iProc, 0, status);
+        MPI::Status status;
+        mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, pImLocal->getMPIDataType(), iProc, 0, status);
       }
     }
   }
@@ -398,7 +388,7 @@ void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImL
 #endif
 
 #ifndef NO_MAIN
-int 
+int
 main (int argc, char* argv[])
 {
   int retval = 1;
@@ -406,9 +396,9 @@ main (int argc, char* argv[])
   try {
     retval = phm2if_main(argc, argv);
   } catch (exception e) {
-    cerr << "Exception: " << e.what() << endl;
+    std::cerr << "Exception: " << e.what() << std::endl;
   } catch (...) {
-    cerr << "Unknown exception" << endl;
+    std::cerr << "Unknown exception\n";
   }
 
   return (retval);