tools/phm2if.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          phm2if.cpp
   5 **   Purpose:       Convert an phantom object to an image file
   6 **   Programmer:    Kevin Rosenberg
   7 **   Date Started:  1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: phm2if.cpp,v 1.11 2000/08/27 20:32:55 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29 #include "timer.h"
  30
  31
  32 enum { O_PHANTOM, O_DESC, O_NSAMPLE, O_FILTER, O_TRACE, O_VERBOSE, O_HELP,
  33        O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
  34
  35 static struct option my_options[] =
  36 {
  37   {"phantom", 1, 0, O_PHANTOM},
  38   {"phmfile", 1, 0, O_PHMFILE},
  39   {"desc", 1, 0, O_DESC},
  40   {"nsample", 1, 0, O_NSAMPLE},
  41   {"filter", 1, 0, O_FILTER},
  42   {"filter-domain", 1, 0, O_FILTER_DOMAIN},
  43   {"filter-bw", 1, 0, O_FILTER_BW},
  44   {"filter-param", 1, 0, O_FILTER_PARAM},
  45   {"trace", 1, 0, O_TRACE},
  46   {"verbose", 0, 0, O_VERBOSE},
  47   {"debug", 0, 0, O_DEBUG},
  48   {"help", 0, 0, O_HELP},
  49   {"version", 0, 0, O_VERSION},
  50   {0, 0, 0, 0}
  51 };
  52
  53 static const char* g_szIdStr = "$Id: phm2if.cpp,v 1.11 2000/08/27 20:32:55 kevin Exp $";
  54
  55 void
  56 phm2if_usage (const char *program)
  57 {
  58   cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]" << endl;
  59   cout << "Generate phantom image from a predefined --phantom or a --phmfile" << endl;
  60   cout << endl;
  61   cout << "     outfile         Name of output file for image" << endl;
  62   cout << "     nx              Number of pixels X-axis" << endl;
  63   cout << "     ny              Number of pixels Y-axis" << endl;
  64   cout << "     --phantom       Phantom to use for projection" << endl;
  65   cout << "        herman       Herman head phantom" << endl;
  66   cout << "        herman-b     Herman head phantom (Bordered)" << endl;
  67   cout << "        shepp-logan  Shepp-Logan head phantom" << endl;
  68   cout << "        shepp-logan-b Shepp-Logan head phantom (Bordered)" << endl;
  69   cout << "        unitpulse    Unit pulse phantom" << endl;
  70   cout << "     --phmfile       Generate Phantom from phantom file" << endl;
  71   cout << "     --filter        Generate Phantom from a filter function" << endl;
  72   cout << "       abs_bandlimit Abs * Bandlimiting" << endl;
  73   cout << "       abs_sinc      Abs * Sinc" << endl;
  74   cout << "       abs_cos       Abs * Cosine" << endl;
  75   cout << "       abs_hamming   Abs * Hamming" << endl;
  76   cout << "       shepp         Shepp-Logan" << endl;
  77   cout << "       bandlimit     Bandlimiting" << endl;
  78   cout << "       sinc          Sinc" << endl;
  79   cout << "       cos           Cosine" << endl;
  80   cout << "       triangle      Triangle" << endl;
  81   cout << "       hamming       Hamming" << endl;
  82   cout << "     --filter-param  Alpha level for Hamming filter" << endl;
  83   cout << "     --filter-domain Set domain of filter" << endl;
  84   cout << "         spatial     Spatial domain (default)" << endl;
  85   cout << "         freq        Frequency domain" << endl;
  86   cout << "     --filter-bw     Filter bandwidth (default = 1)" << endl;
  87   cout << "     --desc          Description of raysum" << endl;
  88   cout << "     --nsample       Number of samples per axis per pixel (default = 1)" << endl;
  89   cout << "     --trace         Trace level to use" << endl;
  90   cout << "        none         No tracing (default)" << endl;
  91   cout << "        console      Trace text level" << endl;
  92   cout << "        phantom      Trace phantom" << endl;
  93   cout << "        proj         Trace projections" << endl;
  94   cout << "        plot         Trace plot" << endl;
  95   cout << "        clipping     Trace clipping" << endl;
  96   cout << "     --debug         Debug mode" << endl;
  97   cout << "     --verbose       Verbose mode" << endl;
  98   cout << "     --version       Print version" << endl;
  99   cout << "     --help          Print this help message" << endl;
 100 }
 101
 102 #ifdef HAVE_MPI
 103 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug);
 104 #endif
 105
 106 int
 107 phm2if_main (int argc, char* argv[])
 108 {
 109   ImageFile* pImGlobal = NULL;
 110   Phantom phm;
 111   string optPhmName (Phantom::convertPhantomIDToName(Phantom::PHM_HERMAN));
 112   string optFilterName;
 113   string optDomainName (SignalFilter::convertDomainIDToName (SignalFilter::DOMAIN_SPATIAL));
 114   string optOutFilename;
 115   string optDesc;
 116   string optPhmFilename;
 117   int opt_nx = 0;
 118   int opt_ny = 0;
 119   int opt_nsample = 1;
 120   double optFilterParam = 1.;
 121   double optFilterBW = 1.;
 122   int optTrace = Trace::TRACE_NONE;
 123   bool optVerbose = false;
 124   bool optDebug = false;
 125   char *endptr = NULL;
 126   char *endstr;
 127   Timer timerProgram;
 128
 129 #ifdef HAVE_MPI
 130   ImageFile* pImLocal = NULL;
 131   MPIWorld mpiWorld (argc, argv);
 132   if (mpiWorld.getRank() == 0) {
 133 #endif
 134     while (1) {
 135       int c = getopt_long(argc, argv, "", my_options, NULL);
 136       if (c == -1)
 137         break;
 138
 139       switch (c) {
 140       case O_PHANTOM:
 141         optPhmName = optarg;
 142         break;
 143       case O_PHMFILE:
 144         optPhmFilename = optarg;
 145         break;
 146       case O_VERBOSE:
 147         optVerbose = true;
 148         break;
 149       case O_DEBUG:
 150         optDebug = true;
 151         break;
 152       case O_TRACE:
 153         if ((optTrace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
 154           phm2if_usage(argv[0]);
 155           return (1);
 156         }
 157         break;
 158       case O_FILTER:
 159         optFilterName = optarg;
 160         break;
 161       case O_FILTER_DOMAIN:
 162         optDomainName = optarg;
 163         break;
 164       case O_DESC:
 165         optDesc =  optarg;
 166         break;
 167       case O_FILTER_BW:
 168         optFilterBW = strtod(optarg, &endptr);
 169         endstr = optarg + strlen(optarg);
 170         if (endptr != endstr) {
 171           sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
 172           phm2if_usage(argv[0]);
 173           return (1);
 174         }
 175         break;
 176       case O_FILTER_PARAM:
 177         optFilterParam = strtod(optarg, &endptr);
 178         endstr = optarg + strlen(optarg);
 179         if (endptr != endstr) {
 180           sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
 181           phm2if_usage(argv[0]);
 182           return (1);
 183         }
 184         break;
 185       case O_NSAMPLE:
 186         opt_nsample = strtol(optarg, &endptr, 10);
 187         endstr = optarg + strlen(optarg);
 188         if (endptr != endstr) {
 189           sys_error(ERR_SEVERE,"Error setting --nsample to %s\n", optarg);
 190           phm2if_usage(argv[0]);
 191           return (1);
 192         }
 193         break;
 194         case O_VERSION:
 195 #ifdef VERSION
 196           cout << "Version " << VERSION << endl << g_szIdStr << endl;
 197 #else
 198           cerr << "Unknown version number" << endl;
 199 #endif
 200       case O_HELP:
 201       case '?':
 202         phm2if_usage(argv[0]);
 203         return (0);
 204       default:
 205         phm2if_usage(argv[0]);
 206         return (1);
 207       }
 208     }
 209
 210     if (optPhmName == "" && optFilterName == "" && optPhmFilename == "") {
 211       cerr << "No phantom defined" << endl << endl;
 212       phm2if_usage(argv[0]);
 213       return (1);
 214     }
 215
 216     if (optind + 3 != argc) {
 217       phm2if_usage(argv[0]);
 218       return (1);
 219     }
 220     optOutFilename = argv[optind];
 221     opt_nx = strtol(argv[optind+1], &endptr, 10);
 222     endstr = argv[optind+1] + strlen(argv[optind+1]);
 223     if (endptr != endstr) {
 224       sys_error(ERR_SEVERE,"Error setting nx to %s\n", argv[optind+1]);
 225       phm2if_usage(argv[0]);
 226       return (1);
 227     }
 228     opt_ny = strtol(argv[optind+2], &endptr, 10);
 229     endstr = argv[optind+2] + strlen(argv[optind+2]);
 230     if (endptr != endstr) {
 231       sys_error(ERR_SEVERE,"Error setting ny to %s\n", argv[optind+2]);
 232       phm2if_usage(argv[0]);
 233       return (1);
 234     }
 235
 236     ostringstream oss;
 237     oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", ";
 238     if (optPhmFilename != "")
 239       oss << "phantomFile=" << optPhmFilename;
 240     else if (optPhmName != "")
 241       oss << "phantom=" << optPhmName;
 242     else if (optFilterName != "") {
 243       oss << "filter=" << optFilterName << " - " << optDomainName;
 244     }
 245     if (optDesc != "")
 246       oss << ": " << optDesc;
 247     optDesc = oss.str();
 248
 249     if (optPhmName != "") {
 250       phm.createFromPhantom (optPhmName.c_str());
 251       if (phm.fail()) {
 252         cout << phm.failMessage() << endl << endl;
 253         phm2if_usage(argv[0]);
 254         return (1);
 255       }
 256     }
 257
 258     if (optPhmFilename != "") {
 259       phm.createFromFile(optPhmFilename.c_str());
 260 #ifdef HAVE_MPI
 261       if (mpiWorld.getRank() == 0)
 262         cerr << "Can't use phantom from file in MPI mode" << endl;
 263       return (1);
 264 #endif
 265     }
 266
 267     if (optVerbose)
 268       cout << "Rasterize Phantom to Image" << endl << endl;
 269 #ifdef HAVE_MPI
 270   }
 271 #endif
 272
 273 #ifdef HAVE_MPI
 274   TimerCollectiveMPI timerBcast (mpiWorld.getComm());
 275   mpiWorld.BcastString (optPhmName);
 276   mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
 277   mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
 278   mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
 279   mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
 280   mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
 281   mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
 282   mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
 283   mpiWorld.getComm().Bcast (&optFilterBW, 1, MPI::DOUBLE, 0);
 284
 285   mpiWorld.BcastString (optFilterName);
 286   mpiWorld.BcastString (optDomainName);
 287
 288   if (optVerbose)
 289     timerBcast.timerEndAndReport ("Time to broadcast variables");
 290
 291   mpiWorld.setTotalWorkUnits (opt_nx);
 292
 293   if (mpiWorld.getRank() > 0 && optPhmName != "")
 294       phm.createFromPhantom (optPhmName.c_str());
 295
 296   if (mpiWorld.getRank() == 0) {
 297     pImGlobal = new ImageFile (opt_nx, opt_ny);
 298   }
 299   pImLocal = new ImageFile (opt_nx, opt_ny);
 300 #else
 301   pImGlobal = new ImageFile (opt_nx, opt_ny);
 302 #endif
 303
 304   ImageFileArray v;
 305 #ifdef HAVE_MPI
 306   if (mpiWorld.getRank() == 0)
 307     v = pImGlobal->getArray ();
 308
 309   if (phm.getComposition() == P_UNIT_PULSE) {
 310     if (mpiWorld.getRank() == 0) {
 311       v[opt_nx/2][opt_ny/2] = 1.;
 312     }
 313   } else if (optFilterName != "") {
 314     if (mpiWorld.getRank() == 0) {
 315       pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
 316     }
 317   } else {
 318     TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
 319     phm.convertToImagefile (*pImLocal, opt_nsample, optTrace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits());
 320     if (optVerbose)
 321       timerRasterize.timerEndAndReport ("Time to rasterize phantom");
 322
 323     TimerCollectiveMPI timerGather (mpiWorld.getComm());
 324     mpi_gather_image (mpiWorld, pImGlobal, pImLocal, optDebug);
 325     if (optVerbose)
 326       timerGather.timerEndAndReport ("Time to gather image");
 327   }
 328 #else
 329   v = pImGlobal->getArray ();
 330   if (phm.getComposition() == P_UNIT_PULSE) {
 331     v[opt_nx/2][opt_ny/2] = 1.;
 332   } else if (optFilterName != "") {
 333     pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
 334   } else {
 335 #if HAVE_SGP
 336       if (optTrace >= Trace::TRACE_PHANTOM)
 337         phm.show();
 338 #endif
 339       phm.convertToImagefile (*pImGlobal, opt_nsample, optTrace);
 340   }
 341 #endif
 342
 343 #ifdef HAVE_MPI
 344   if (mpiWorld.getRank() == 0)
 345 #endif
 346   {
 347     double calctime = timerProgram.timerEnd ();
 348     pImGlobal->labelAdd (Array2dFileLabel::L_HISTORY, optDesc.c_str(), calctime);
 349     pImGlobal->fileWrite (optOutFilename.c_str());
 350     if (optVerbose)
 351       cout << "Time to rasterized phantom: " << calctime << " seconds" << endl;
 352
 353     if (optTrace >= Trace::TRACE_PHANTOM) {
 354       double dmin, dmax;
 355       int nscale;
 356
 357       printf ("Enter display size scale (nominal = 1): ");
 358       scanf ("%d", &nscale);
 359       printf ("Enter minimum and maximum densities (min, max): ");
 360       scanf ("%lf %lf", &dmin, &dmax);
 361       pImGlobal->displayScaling (nscale, dmin, dmax);
 362     }
 363   }
 364
 365   delete pImGlobal;
 366 #ifdef HAVE_MPI
 367   delete pImLocal;
 368 #endif
 369
 370   return (0);
 371 }
 372
 373
 374
 375 #ifdef HAVE_MPI
 376 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug)
 377 {
 378   ImageFileArray vLocal = pImLocal->getArray();
 379   ImageFileArray vGlobal = NULL;
 380   int nyLocal = pImLocal->ny();
 381
 382   if (mpiWorld.getRank() == 0)
 383     vGlobal = pImGlobal->getArray();
 384
 385   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
 386     mpiWorld.getComm().Send(vLocal[iw], nyLocal, pImLocal->getMPIDataType(), 0, 0);
 387
 388   if (mpiWorld.getRank() == 0) {
 389     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
 390       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
 391         MPI::Status status;
 392         mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, pImLocal->getMPIDataType(), iProc, 0, status);
 393       }
 394     }
 395   }
 396
 397 }
 398 #endif
 399
 400 #ifndef NO_MAIN
 401 int
 402 main (int argc, char* argv[])
 403 {
 404   int retval = 1;
 405
 406   try {
 407     retval = phm2if_main(argc, argv);
 408   } catch (exception e) {
 409     cerr << "Exception: " << e.what() << endl;
 410   } catch (...) {
 411     cerr << "Unknown exception" << endl;
 412   }
 413
 414   return (retval);
 415 }
 416 #endif