tools/phm2pj.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          phm2pj.cpp
   5 **   Purpose:       Take projections of a phantom object
   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: phm2pj.cpp,v 1.25 2001/02/11 04:56:38 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_NRAY, O_ROTANGLE, O_PHMFILE, O_GEOMETRY, O_FOCAL_LENGTH, O_VIEW_RATIO, O_SCAN_RATIO,
  33        O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
  34
  35 static struct option phm2pj_options[] =
  36 {
  37   {"phantom", 1, 0, O_PHANTOM},
  38   {"phmfile", 1, 0, O_PHMFILE},
  39   {"desc", 1, 0, O_DESC},
  40   {"nray", 1, 0, O_NRAY},
  41   {"rotangle", 1, 0, O_ROTANGLE},
  42   {"geometry", 1, 0, O_GEOMETRY},
  43   {"focal-length", 1, 0, O_FOCAL_LENGTH},
  44   {"view-ratio", 1, 0, O_VIEW_RATIO},
  45   {"scan-ratio", 1, 0, O_SCAN_RATIO},
  46   {"trace", 1, 0, O_TRACE},
  47   {"verbose", 0, 0, O_VERBOSE},
  48   {"help", 0, 0, O_HELP},
  49   {"debug", 0, 0, O_DEBUG},
  50   {"version", 0, 0, O_VERSION},
  51   {0, 0, 0, 0}
  52 };
  53
  54 static const char* g_szIdStr = "$Id: phm2pj.cpp,v 1.25 2001/02/11 04:56:38 kevin Exp $";
  55
  56
  57 void
  58 phm2pj_usage (const char *program)
  59 {
  60   std::cout << "usage: " << fileBasename(program) << " outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]\n";
  61   std::cout << "Calculate (projections) through phantom object, either a predefined --phantom or a --phmfile\n\n";
  62   std::cout << "     outfile          Name of output file for projections\n";
  63   std::cout << "     ndet             Number of detectors\n";
  64   std::cout << "     nview            Number of rotated views\n";
  65   std::cout << "     --phantom        Phantom to use for projection\n";
  66   std::cout << "        herman        Herman head phantom\n";
  67   std::cout << "        herman-b      Herman head phantom (Bordered)\n";
  68   std::cout << "        shepp-logan   Shepp-Logan head phantom\n";
  69   std::cout << "        shepp-logan-b Shepp-Logan head phantom (Bordered)\n";
  70   std::cout << "        unitpulse     Unit pulse phantom\n";
  71   std::cout << "     --phmfile        Get Phantom from phantom file\n";
  72   std::cout << "     --desc           Description of raysum\n";
  73   std::cout << "     --nray           Number of rays per detector (default = 1)\n";
  74   std::cout << "     --rotangle       Degrees to rotate view through (multiple of PI)\n";
  75   std::cout << "                      (default = select appropriate for geometry)\n";
  76   std::cout << "     --geometry       Geometry of scanning\n";
  77   std::cout << "        parallel      Parallel scan beams (default)\n";
  78   std::cout << "        equilinear    Equilinear divergent scan beams\n";
  79   std::cout << "        equiangular   Equiangular divergent scan beams\n";
  80   std::cout << "     --focal-length   Focal length ratio (ratio to radius of phantom)\n";
  81   std::cout << "                      (default = 1)\n";
  82   std::cout << "     --view-ratio     Length to view (view diameter to phantom diameter)\n";
  83   std::cout << "                      (default = 1)\n";
  84   std::cout << "     --scan-ratio     Length to scan (scan diameter to view diameter)\n";
  85   std::cout << "                      (default = 1)\n";
  86   std::cout << "     --trace          Trace level to use\n";
  87   std::cout << "        none          No tracing (default)\n";
  88   std::cout << "        console       Trace text level\n";
  89   std::cout << "        phantom       Trace phantom image\n";
  90   std::cout << "        proj          Trace projections\n";
  91   std::cout << "        plot          Trace plot\n";
  92   std::cout << "        clipping      Trace clipping\n";
  93   std::cout << "     --verbose        Verbose mode\n";
  94   std::cout << "     --debug          Debug mode\n";
  95   std::cout << "     --version        Print version\n";
  96   std::cout << "     --help           Print this help message\n";
  97 }
  98
  99 #ifdef HAVE_MPI
 100 void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug);
 101 #endif
 102
 103 int
 104 phm2pj_main (int argc, char* const argv[])
 105 {
 106   Phantom phm;
 107   std::string optGeometryName = Scanner::convertGeometryIDToName(Scanner::GEOMETRY_PARALLEL);
 108   char *opt_outfile = NULL;
 109   std::string opt_desc;
 110   std::string optPhmName;
 111   std::string optPhmFileName;
 112   int opt_ndet;
 113   int opt_nview;
 114   int opt_nray = 1;
 115   double dOptFocalLength = 1.;
 116   double dOptViewRatio = 1.;
 117   double dOptScanRatio = 1.;
 118   int opt_trace = Trace::TRACE_NONE;
 119   int opt_verbose = 0;
 120   int opt_debug = 0;
 121   double opt_rotangle = -1;
 122   char* endptr = NULL;
 123   char* endstr;
 124
 125 #ifdef HAVE_MPI
 126   MPIWorld mpiWorld (argc, argv);
 127 #endif
 128
 129   Timer timerProgram;
 130
 131 #ifdef HAVE_MPI
 132   if (mpiWorld.getRank() == 0) {
 133 #endif
 134     while (1) {
 135       int c = getopt_long(argc, argv, "", phm2pj_options, NULL);
 136
 137       if (c == -1)
 138         break;
 139
 140       switch (c) {
 141       case O_PHANTOM:
 142         optPhmName = optarg;
 143         break;
 144       case O_PHMFILE:
 145         optPhmFileName = optarg;
 146         break;
 147       case O_VERBOSE:
 148         opt_verbose = 1;
 149         break;
 150       case O_DEBUG:
 151         opt_debug = 1;
 152         break;
 153         break;
 154       case O_TRACE:
 155         if ((opt_trace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
 156           phm2pj_usage(argv[0]);
 157           return (1);
 158         }
 159         break;
 160       case O_DESC:
 161         opt_desc = optarg;
 162         break;
 163       case O_ROTANGLE:
 164         opt_rotangle = strtod(optarg, &endptr);
 165         endstr = optarg + strlen(optarg);
 166         if (endptr != endstr) {
 167           std::cerr << "Error setting --rotangle to " << optarg << std::endl;
 168           phm2pj_usage(argv[0]);
 169           return (1);
 170         }
 171         break;
 172       case O_GEOMETRY:
 173         optGeometryName = optarg;
 174         break;
 175       case O_FOCAL_LENGTH:
 176         dOptFocalLength = strtod(optarg, &endptr);
 177         endstr = optarg + strlen(optarg);
 178         if (endptr != endstr) {
 179           std::cerr << "Error setting --focal-length to " << optarg << std::endl;
 180           phm2pj_usage(argv[0]);
 181           return (1);
 182         }
 183         break;
 184   case O_VIEW_RATIO:
 185         dOptViewRatio = strtod(optarg, &endptr);
 186         endstr = optarg + strlen(optarg);
 187         if (endptr != endstr) {
 188           std::cerr << "Error setting --view-ratio to " << optarg << std::endl;
 189           phm2pj_usage(argv[0]);
 190           return (1);
 191         }
 192         break;
 193   case O_SCAN_RATIO:
 194         dOptScanRatio = strtod(optarg, &endptr);
 195         endstr = optarg + strlen(optarg);
 196         if (endptr != endstr) {
 197           std::cerr << "Error setting --scan-ratio to " << optarg << std::endl;
 198           phm2pj_usage(argv[0]);
 199           return (1);
 200         }
 201         break;
 202       case O_NRAY:
 203         opt_nray = strtol(optarg, &endptr, 10);
 204         endstr = optarg + strlen(optarg);
 205         if (endptr != endstr) {
 206           std::cerr << "Error setting --nray to %s" << optarg << std::endl;
 207           phm2pj_usage(argv[0]);
 208           return (1);
 209         }
 210         break;
 211         case O_VERSION:
 212 #ifdef VERSION
 213           std::cout << "Version: " << VERSION << std::endl << g_szIdStr << std::endl;
 214 #else
 215           std::cout << "Unknown version number\n";
 216 #endif
 217           return (0);
 218       case O_HELP:
 219       case '?':
 220         phm2pj_usage(argv[0]);
 221         return (0);
 222       default:
 223         phm2pj_usage(argv[0]);
 224         return (1);
 225       }
 226     }
 227
 228     if (optPhmName == "" && optPhmFileName == "") {
 229       std::cerr << "No phantom defined\n" << std::endl;
 230       phm2pj_usage(argv[0]);
 231       return (1);
 232     }
 233     if (optind + 3 != argc) {
 234       phm2pj_usage(argv[0]);
 235       return (1);
 236     }
 237
 238     opt_outfile = argv[optind];
 239     opt_ndet = strtol(argv[optind+1], &endptr, 10);
 240     endstr = argv[optind+1] + strlen(argv[optind+1]);
 241     if (endptr != endstr) {
 242       std::cerr << "Error setting --ndet to " << argv[optind+1] << std::endl;
 243       phm2pj_usage(argv[0]);
 244       return (1);
 245     }
 246     opt_nview = strtol(argv[optind+2], &endptr, 10);
 247     endstr = argv[optind+2] + strlen(argv[optind+2]);
 248     if (endptr != endstr) {
 249       std::cerr << "Error setting --nview to " << argv[optind+2] << std::endl;
 250       phm2pj_usage(argv[0]);
 251       return (1);
 252     }
 253
 254     if (opt_rotangle < 0) {
 255       if (optGeometryName.compare ("parallel") == 0)
 256         opt_rotangle = 1;
 257       else
 258         opt_rotangle = 2;
 259     }
 260
 261     std::ostringstream desc;
 262     desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", Geometry=" << optGeometryName << ", ";
 263     if (optPhmFileName.length()) {
 264       desc << "PhantomFile=" << optPhmFileName;
 265     } else if (optPhmName != "") {
 266       desc << "Phantom=" << optPhmName;
 267     }
 268     if (opt_desc.length()) {
 269       desc << ": " << opt_desc;
 270     }
 271     opt_desc = desc.str();
 272
 273     if (optPhmName != "") {
 274       phm.createFromPhantom (optPhmName.c_str());
 275       if (phm.fail()) {
 276         std::cout << phm.failMessage() << std::endl << std::endl;
 277         phm2pj_usage(argv[0]);
 278         return (1);
 279       }
 280     }
 281
 282     if (optPhmFileName != "") {
 283 #ifdef HAVE_MPI
 284       std::cerr << "Can not read phantom from file in MPI mode\n";
 285       return (1);
 286 #endif
 287       phm.createFromFile (optPhmFileName.c_str());
 288     }
 289
 290 #ifdef HAVE_MPI
 291   }
 292 #endif
 293
 294 #ifdef HAVE_MPI
 295   TimerCollectiveMPI timerBcast(mpiWorld.getComm());
 296   mpiWorld.BcastString (optPhmName);
 297   mpiWorld.getComm().Bcast (&opt_rotangle, 1, MPI::DOUBLE, 0);
 298   mpiWorld.getComm().Bcast (&dOptFocalLength, 1, MPI::DOUBLE, 0);
 299   mpiWorld.getComm().Bcast (&dOptViewRatio, 1, MPI::DOUBLE, 0);
 300   mpiWorld.getComm().Bcast (&dOptScanRatio, 1, MPI::DOUBLE, 0);
 301   mpiWorld.getComm().Bcast (&opt_nview, 1, MPI::INT, 0);
 302   mpiWorld.getComm().Bcast (&opt_ndet, 1, MPI::INT, 0);
 303   mpiWorld.getComm().Bcast (&opt_nray, 1, MPI::INT, 0);
 304   mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0);
 305   mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
 306   mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
 307   if (opt_verbose)
 308     timerBcast.timerEndAndReport ("Time to broadcast variables");
 309
 310   if (mpiWorld.getRank() > 0 && optPhmName != "")
 311     phm.createFromPhantom (optPhmName.c_str());
 312 #endif
 313
 314   opt_rotangle *= PI;
 315   Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_nray, opt_rotangle, dOptFocalLength, dOptViewRatio, dOptScanRatio);
 316   if (scanner.fail()) {
 317     std::cout << "Scanner Creation Error: " << scanner.failMessage() << std::endl;
 318     return (1);
 319   }
 320 #ifdef HAVE_MPI
 321   mpiWorld.setTotalWorkUnits (opt_nview);
 322
 323   Projections pjGlobal;
 324   if (mpiWorld.getRank() == 0)
 325     pjGlobal.initFromScanner (scanner);
 326
 327   if (opt_verbose) {
 328     std::ostringstream os;
 329     pjGlobal.printScanInfo(os);
 330     std::cout << os.str();
 331   }
 332
 333   Projections pjLocal (scanner);
 334   pjLocal.setNView (mpiWorld.getMyLocalWorkUnits());
 335
 336   if (opt_debug)
 337     std::cout << "pjLocal->nview = " << pjLocal.nView() << " (process " << mpiWorld.getRank() << ")\n";;
 338
 339   TimerCollectiveMPI timerProject (mpiWorld.getComm());
 340   scanner.collectProjections (pjLocal, phm, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), false, opt_trace);
 341   if (opt_verbose)
 342     timerProject.timerEndAndReport ("Time to collect projections");
 343
 344   TimerCollectiveMPI timerGather (mpiWorld.getComm());
 345   GatherProjectionsMPI (mpiWorld, pjGlobal, pjLocal, opt_debug);
 346   if (opt_verbose)
 347      timerGather.timerEndAndReport ("Time to gather projections");
 348
 349 #else
 350   Projections pjGlobal (scanner);
 351   scanner.collectProjections (pjGlobal, phm, opt_trace);
 352 #endif
 353
 354 #ifdef HAVE_MPI
 355   if (mpiWorld.getRank() == 0)
 356 #endif
 357     {
 358       pjGlobal.setCalcTime (timerProgram.timerEnd());
 359       pjGlobal.setRemark (opt_desc);
 360       pjGlobal.write (opt_outfile);
 361       if (opt_verbose) {
 362         phm.print (std::cout);
 363         std::cout << std::endl;
 364               std::ostringstream os;
 365               pjGlobal.printScanInfo (os);
 366               std::cout << os.str() << std::endl;
 367               std::cout << "  Remark: " << pjGlobal.remark() << std::endl;
 368               std::cout << "Run time: " << pjGlobal.calcTime() << " seconds\n";
 369       }
 370     }
 371
 372   return (0);
 373 }
 374
 375
 376 /* FUNCTION
 377  *    GatherProjectionsMPI
 378  *
 379  * SYNOPSIS
 380  *    Gather's raysums from all processes in pjLocal in pjGlobal in process 0
 381  */
 382
 383 #ifdef HAVE_MPI
 384 void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug)
 385 {
 386   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
 387     DetectorArray& detArray = pjLocal.getDetectorArray(iw);
 388     double viewAngle = detArray.viewAngle();
 389     int nDet = detArray.nDet();
 390     DetectorValue* detval = detArray.detValues();
 391
 392     mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, 0, 0);
 393     mpiWorld.getComm().Send(&nDet, 1, MPI::INT, 0, 0);
 394     mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, 0, 0);
 395   }
 396
 397   if (mpiWorld.getRank() == 0) {
 398     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
 399       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
 400         MPI::Status status;
 401         double viewAngle;
 402         int nDet;
 403         DetectorArray& detArray = pjGlobal.getDetectorArray(iw);
 404         DetectorValue* detval = detArray.detValues();
 405
 406         mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, iProc, 0, status);
 407         mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, iProc, 0, status);
 408         mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, iProc, 0, status);
 409         detArray.setViewAngle (viewAngle);
 410       }
 411     }
 412   }
 413 }
 414 #endif
 415
 416
 417 #ifndef NO_MAIN
 418 int
 419 main (int argc, char* argv[])
 420 {
 421   int retval = 1;
 422
 423   try {
 424     retval = phm2pj_main(argc, argv);
 425 #if HAVE_DMALLOC
 426     //    dmalloc_shutdown();
 427 #endif
 428   } catch (exception e) {
 429     std::cerr << "Exception: " << e.what() << std::endl;
 430   } catch (...) {
 431     std::cerr << "Unknown exception\n";
 432   }
 433
 434   return (retval);
 435 }
 436 #endif
 437