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