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.22 2001/02/08 06:25:07 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_VIEW_RATIO, 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   {"view-ratio", 1, 0, O_VIEW_RATIO},
  47   {"verbose", 0, 0, O_VERBOSE},
  48   {"debug", 0, 0, O_DEBUG},
  49   {"help", 0, 0, O_HELP},
  50   {"version", 0, 0, O_VERSION},
  51   {0, 0, 0, 0}
  52 };
  53
  54 static const char* g_szIdStr = "$Id: phm2if.cpp,v 1.22 2001/02/08 06:25:07 kevin Exp $";
  55
  56 void
  57 phm2if_usage (const char *program)
  58 {
  59   std::cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]\n";
  60   std::cout << "Generate phantom image from a predefined --phantom or a --phmfile\n";
  61   std::cout << std::endl;
  62   std::cout << "     outfile         Name of output file for image\n";
  63   std::cout << "     nx              Number of pixels X-axis\n";
  64   std::cout << "     ny              Number of pixels Y-axis\n";
  65   std::cout << "     --view-ratio    View diameter to phantom diameter ratio (default = 1)\n";
  66   std::cout << "     --phantom       Phantom to use for projection\n";
  67   std::cout << "        herman       Herman head phantom\n";
  68   std::cout << "        herman-b     Herman head phantom (Bordered)\n";
  69   std::cout << "        shepp-logan  Shepp-Logan head phantom\n";
  70   std::cout << "        shepp-logan-b Shepp-Logan head phantom (Bordered)\n";
  71   std::cout << "        unitpulse    Unit pulse phantom\n";
  72   std::cout << "     --phmfile       Generate Phantom from phantom file\n";
  73   std::cout << "     --filter        Generate Phantom from a filter function\n";
  74   std::cout << "       abs_bandlimit Abs * Bandlimiting\n";
  75   std::cout << "       abs_sinc      Abs * Sinc\n";
  76   std::cout << "       abs_cos       Abs * Cosine\n";
  77   std::cout << "       abs_hamming   Abs * Hamming\n";
  78   std::cout << "       shepp         Shepp-Logan\n";
  79   std::cout << "       bandlimit     Bandlimiting\n";
  80   std::cout << "       sinc          Sinc\n";
  81   std::cout << "       cos           Cosine\n";
  82   std::cout << "       triangle      Triangle\n";
  83   std::cout << "       hamming       Hamming\n";
  84   std::cout << "     --filter-param  Alpha level for Hamming filter\n";
  85   std::cout << "     --filter-domain Set domain of filter\n";
  86   std::cout << "         spatial     Spatial domain (default)\n";
  87   std::cout << "         freq        Frequency domain\n";
  88   std::cout << "     --filter-bw     Filter bandwidth (default = 1)\n";
  89   std::cout << "     --desc          Description of raysum\n";
  90   std::cout << "     --nsample       Number of samples per axis per pixel (default = 1)\n";
  91   std::cout << "     --trace         Trace level to use\n";
  92   std::cout << "        none         No tracing (default)\n";
  93   std::cout << "        console      Trace text level\n";
  94   std::cout << "        phantom      Trace phantom\n";
  95   std::cout << "        proj         Trace projections\n";
  96   std::cout << "        plot         Trace plot\n";
  97   std::cout << "        clipping     Trace clipping\n";
  98   std::cout << "     --debug         Debug mode\n";
  99   std::cout << "     --verbose       Verbose mode\n";
 100   std::cout << "     --version       Print version\n";
 101   std::cout << "     --help          Print this help message\n";
 102 }
 103
 104 #ifdef HAVE_MPI
 105 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug);
 106 #endif
 107
 108 int
 109 phm2if_main (int argc, char* const argv[])
 110 {
 111   ImageFile* pImGlobal = NULL;
 112   Phantom phm;
 113   std::string optPhmName;
 114   std::string optFilterName;
 115   std::string optDomainName (SignalFilter::convertDomainIDToName (SignalFilter::DOMAIN_SPATIAL));
 116   std::string optOutFilename;
 117   std::string optDesc;
 118   std::string optPhmFilename;
 119   int opt_nx = 0;
 120   int opt_ny = 0;
 121   int opt_nsample = 1;
 122   double optViewRatio = 1.;
 123   double optFilterParam = 1.;
 124   double optFilterBW = 1.;
 125   int optTrace = Trace::TRACE_NONE;
 126   bool optVerbose = false;
 127   bool optDebug = false;
 128   char *endptr = NULL;
 129   char *endstr;
 130   Timer timerProgram;
 131
 132 #ifdef HAVE_MPI
 133   ImageFile* pImLocal = NULL;
 134   MPIWorld mpiWorld (argc, argv);
 135   if (mpiWorld.getRank() == 0) {
 136 #endif
 137     while (1) {
 138       int c = getopt_long(argc, argv, "", my_options, NULL);
 139       if (c == -1)
 140         break;
 141
 142       switch (c) {
 143       case O_PHANTOM:
 144         optPhmName = optarg;
 145         break;
 146       case O_PHMFILE:
 147         optPhmFilename = optarg;
 148         break;
 149       case O_VERBOSE:
 150         optVerbose = true;
 151         break;
 152       case O_DEBUG:
 153         optDebug = true;
 154         break;
 155       case O_TRACE:
 156         if ((optTrace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
 157           phm2if_usage(argv[0]);
 158           return (1);
 159         }
 160         break;
 161       case O_FILTER:
 162         optFilterName = optarg;
 163         break;
 164       case O_FILTER_DOMAIN:
 165         optDomainName = optarg;
 166         break;
 167       case O_DESC:
 168         optDesc =  optarg;
 169         break;
 170       case O_FILTER_BW:
 171         optFilterBW = strtod(optarg, &endptr);
 172         endstr = optarg + strlen(optarg);
 173         if (endptr != endstr) {
 174           sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
 175           phm2if_usage(argv[0]);
 176           return (1);
 177         }
 178         break;
 179       case O_VIEW_RATIO:
 180         optViewRatio = strtod(optarg, &endptr);
 181         endstr = optarg + strlen(optarg);
 182         if (endptr != endstr) {
 183           sys_error(ERR_SEVERE,"Error setting --view-ratio to %s\n", optarg);
 184           phm2if_usage(argv[0]);
 185           return (1);
 186         }
 187         break;
 188       case O_FILTER_PARAM:
 189         optFilterParam = strtod(optarg, &endptr);
 190         endstr = optarg + strlen(optarg);
 191         if (endptr != endstr) {
 192           sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
 193           phm2if_usage(argv[0]);
 194           return (1);
 195         }
 196         break;
 197       case O_NSAMPLE:
 198         opt_nsample = strtol(optarg, &endptr, 10);
 199         endstr = optarg + strlen(optarg);
 200         if (endptr != endstr) {
 201           sys_error(ERR_SEVERE,"Error setting --nsample to %s\n", optarg);
 202           phm2if_usage(argv[0]);
 203           return (1);
 204         }
 205         break;
 206       case O_VERSION:
 207 #ifdef VERSION
 208         std::cout << "Version " << VERSION << std::endl << g_szIdStr << std::endl;
 209 #else
 210         std::cerr << "Unknown version number\n";
 211 #endif
 212       case O_HELP:
 213       case '?':
 214         phm2if_usage(argv[0]);
 215         return (0);
 216       default:
 217         phm2if_usage(argv[0]);
 218         return (1);
 219       }
 220     }
 221
 222     if (optPhmName == "" && optFilterName == "" && optPhmFilename == "") {
 223       std::cerr << "No phantom defined\n" << std::endl;
 224       phm2if_usage(argv[0]);
 225       return (1);
 226     }
 227
 228     if (optind + 3 != argc) {
 229       phm2if_usage(argv[0]);
 230       return (1);
 231     }
 232     optOutFilename = argv[optind];
 233     opt_nx = strtol(argv[optind+1], &endptr, 10);
 234     endstr = argv[optind+1] + strlen(argv[optind+1]);
 235     if (endptr != endstr) {
 236       sys_error(ERR_SEVERE,"Error setting nx to %s\n", argv[optind+1]);
 237       phm2if_usage(argv[0]);
 238       return (1);
 239     }
 240     opt_ny = strtol(argv[optind+2], &endptr, 10);
 241     endstr = argv[optind+2] + strlen(argv[optind+2]);
 242     if (endptr != endstr) {
 243       sys_error(ERR_SEVERE,"Error setting ny to %s\n", argv[optind+2]);
 244       phm2if_usage(argv[0]);
 245       return (1);
 246     }
 247
 248     std::ostringstream oss;
 249     oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", viewRatio=" << optViewRatio << ", nsample=" << opt_nsample << ", ";
 250     if (optPhmFilename != "")
 251       oss << "phantomFile=" << optPhmFilename;
 252     else if (optPhmName != "")
 253       oss << "phantom=" << optPhmName;
 254     else if (optFilterName != "") {
 255       oss << "filter=" << optFilterName << " - " << optDomainName;
 256     }
 257     if (optDesc != "")
 258       oss << ": " << optDesc;
 259     optDesc = oss.str();
 260
 261     if (optPhmName != "") {
 262       phm.createFromPhantom (optPhmName.c_str());
 263       if (phm.fail()) {
 264         std::cout << phm.failMessage() << std::endl << std::endl;
 265         phm2if_usage(argv[0]);
 266         return (1);
 267       }
 268     }
 269
 270     if (optPhmFilename != "") {
 271       phm.createFromFile(optPhmFilename.c_str());
 272 #ifdef HAVE_MPI
 273       if (mpiWorld.getRank() == 0)
 274         std::cerr << "Can't use phantom from file in MPI mode\n";
 275       return (1);
 276 #endif
 277     }
 278
 279     if (optVerbose)
 280       std::cout << "Rasterize Phantom to Image\n" << std::endl;
 281 #ifdef HAVE_MPI
 282   }
 283 #endif
 284
 285 #ifdef HAVE_MPI
 286   TimerCollectiveMPI timerBcast (mpiWorld.getComm());
 287   mpiWorld.BcastString (optPhmName);
 288   mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
 289   mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
 290   mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
 291   mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
 292   mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
 293   mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
 294   mpiWorld.getComm().Bcast *&optViewRatio, 1, MPI::DOUBLE, 0);
 295   mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
 296   mpiWorld.getComm().Bcast (&optFilterBW, 1, MPI::DOUBLE, 0);
 297
 298   mpiWorld.BcastString (optFilterName);
 299   mpiWorld.BcastString (optDomainName);
 300
 301   if (optVerbose)
 302     timerBcast.timerEndAndReport ("Time to broadcast variables");
 303
 304   mpiWorld.setTotalWorkUnits (opt_nx);
 305
 306   if (mpiWorld.getRank() > 0 && optPhmName != "")
 307     phm.createFromPhantom (optPhmName.c_str());
 308
 309   if (mpiWorld.getRank() == 0) {
 310     pImGlobal = new ImageFile (opt_nx, opt_ny);
 311   }
 312   pImLocal = new ImageFile (opt_nx, opt_ny);
 313 #else
 314   pImGlobal = new ImageFile (opt_nx, opt_ny);
 315 #endif
 316
 317   ImageFileArray v = NULL;
 318 #ifdef HAVE_MPI
 319   if (mpiWorld.getRank() == 0)
 320     v = pImGlobal->getArray ();
 321
 322   if (phm.getComposition() == P_UNIT_PULSE) {
 323     if (mpiWorld.getRank() == 0) {
 324       v[opt_nx/2][opt_ny/2] = 1.;
 325     }
 326   } else if (optFilterName != "") {
 327     if (mpiWorld.getRank() == 0) {
 328       pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
 329     }
 330   } else {
 331     TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
 332     phm.convertToImagefile (*pImLocal, optViewRatio, opt_nsample, optTrace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), false);
 333     if (optVerbose)
 334       timerRasterize.timerEndAndReport ("Time to rasterize phantom");
 335
 336     TimerCollectiveMPI timerGather (mpiWorld.getComm());
 337     mpi_gather_image (mpiWorld, pImGlobal, pImLocal, optDebug);
 338     if (optVerbose)
 339       timerGather.timerEndAndReport ("Time to gather image");
 340   }
 341 #else
 342   v = pImGlobal->getArray ();
 343   if (phm.getComposition() == P_UNIT_PULSE) {
 344     v[opt_nx/2][opt_ny/2] = 1.;
 345   } else if (optFilterName != "") {
 346     pImGlobal->filterResponse (optDomainName.c_str(), optFilterBW, optFilterName.c_str(), optFilterParam);
 347   } else {
 348     phm.convertToImagefile (*pImGlobal, optViewRatio, opt_nsample, optTrace);
 349   }
 350 #endif
 351
 352 #ifdef HAVE_MPI
 353   if (mpiWorld.getRank() == 0)
 354 #endif
 355   {
 356     double calctime = timerProgram.timerEnd ();
 357     pImGlobal->labelAdd (Array2dFileLabel::L_HISTORY, optDesc.c_str(), calctime);
 358     pImGlobal->fileWrite (optOutFilename.c_str());
 359     if (optVerbose)
 360       std::cout << "Time to rasterize phantom: " << calctime << " seconds\n";
 361   }
 362
 363   delete pImGlobal;
 364 #ifdef HAVE_MPI
 365   delete pImLocal;
 366 #endif
 367
 368   return (0);
 369 }
 370
 371
 372
 373 #ifdef HAVE_MPI
 374 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* pImGlobal, ImageFile* pImLocal, const int optDebug)
 375 {
 376   ImageFileArray vLocal = pImLocal->getArray();
 377   ImageFileArray vGlobal = NULL;
 378   int nyLocal = pImLocal->ny();
 379
 380   if (mpiWorld.getRank() == 0)
 381     vGlobal = pImGlobal->getArray();
 382
 383   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
 384     mpiWorld.getComm().Send(vLocal[iw], nyLocal, pImLocal->getMPIDataType(), 0, 0);
 385
 386   if (mpiWorld.getRank() == 0) {
 387     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
 388       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
 389         MPI::Status status;
 390         mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, pImLocal->getMPIDataType(), iProc, 0, status);
 391       }
 392     }
 393   }
 394
 395 }
 396 #endif
 397
 398 #ifndef NO_MAIN
 399 int
 400 main (int argc, char* argv[])
 401 {
 402   int retval = 1;
 403
 404   try {
 405     retval = phm2if_main(argc, argv);
 406   } catch (exception e) {
 407     std::cerr << "Exception: " << e.what() << std::endl;
 408   } catch (...) {
 409     std::cerr << "Unknown exception\n";
 410   }
 411
 412   return (retval);
 413 }
 414 #endif