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-2009 Kevin Rosenberg
  11 **
  12 **  This program is free software; you can redistribute it and/or modify
  13 **  it under the terms of the GNU General Public License (version 2) as
  14 **  published by the Free Software Foundation.
  15 **
  16 **  This program is distributed in the hope that it will be useful,
  17 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 **  GNU General Public License for more details.
  20 **
  21 **  You should have received a copy of the GNU General Public License
  22 **  along with this program; if not, write to the Free Software
  23 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  24 ******************************************************************************/
  25
  26 #include "ct.h"
  27 #include "timer.h"
  28
  29
  30 enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE, O_IMAGEFILE, O_GEOMETRY, O_FOCAL_LENGTH, O_CENTER_DETECTOR_LENGTH,
  31 O_VIEW_RATIO, O_SCAN_RATIO, O_OFFSETVIEW, O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
  32
  33 static struct option phm2pj_options[] =
  34 {
  35   {"phantom", 1, 0, O_PHANTOM},
  36   {"phmfile", 1, 0, O_PHMFILE},
  37   {"imagefile", 1, 0, O_IMAGEFILE},
  38   {"desc", 1, 0, O_DESC},
  39   {"nray", 1, 0, O_NRAY},
  40   {"rotangle", 1, 0, O_ROTANGLE},
  41   {"geometry", 1, 0, O_GEOMETRY},
  42   {"focal-length", 1, 0, O_FOCAL_LENGTH},
  43   {"center-detector-length", 1, 0, O_CENTER_DETECTOR_LENGTH},
  44   {"offsetview", 1, 0, O_OFFSETVIEW},
  45   {"view-ratio", 1, 0, O_VIEW_RATIO},
  46   {"scan-ratio", 1, 0, O_SCAN_RATIO},
  47   {"trace", 1, 0, O_TRACE},
  48   {"verbose", 0, 0, O_VERBOSE},
  49   {"help", 0, 0, O_HELP},
  50   {"debug", 0, 0, O_DEBUG},
  51   {"version", 0, 0, O_VERSION},
  52   {0, 0, 0, 0}
  53 };
  54
  55 static const char* g_szIdStr = "$Id$";
  56
  57
  58 void
  59 phm2pj_usage (const char *program)
  60 {
  61   std::cout << "usage: " << fileBasename(program) << " outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]\n";
  62   std::cout << "Calculate (projections) through phantom object, either a predefined --phantom or a --phmfile\n\n";
  63   std::cout << "     outfile          Name of output file for projections\n";
  64   std::cout << "     ndet             Number of detectors\n";
  65   std::cout << "     nview            Number of rotated views\n";
  66   std::cout << "     SELECT PHANTOM   Use either --phantom, --phmfile, or --imagefile\n";
  67   std::cout << "     --phantom        Phantom to use for projection\n";
  68   std::cout << "        herman        Herman head phantom\n";
  69   std::cout << "        shepp-logan   Shepp-Logan head phantom\n";
  70   std::cout << "        unit-pulse     Unit pulse phantom\n";
  71   std::cout << "     --phmfile        Get Phantom from phantom file\n";
  72   std::cout << "     --imagefile      Get Phantom from imagefile\n";
  73   std::cout << "     --desc           Add description for projection file\n";
  74   std::cout << "     --nray           Number of rays per detector (default = 1)\n";
  75   std::cout << "     --rotangle       Angle to rotate view through (fraction of a circle)\n";
  76   std::cout << "                      (default = select appropriate for geometry)\n";
  77   std::cout << "     --geometry       Geometry of scanning\n";
  78   std::cout << "        parallel      Parallel scan beams (default)\n";
  79   std::cout << "        equilinear    Equilinear divergent scan beams\n";
  80   std::cout << "        equiangular   Equiangular divergent scan beams\n";
  81   std::cout << "     --focal-length   Focal length ratio (ratio to radius of view area)\n";
  82   std::cout << "                      (default = 2)\n";
  83   std::cout << "     --center-detector-length  Distance from center of phantom to detector array\n";
  84   std::cout << "                      (ratio to radius of view area) (default = 2)\n";
  85   std::cout << "     --view-ratio     Length to view (view diameter to phantom diameter)\n";
  86   std::cout << "                      (default = 1)\n";
  87   std::cout << "     --scan-ratio     Length to scan (scan diameter to view diameter)\n";
  88   std::cout << "                      (default = 1)\n";
  89   std::cout << "     --offsetview     Initial gantry offset in 'views' (default = 0)\n";
  90   std::cout << "     --trace          Trace level to use\n";
  91   std::cout << "        none          No tracing (default)\n";
  92   std::cout << "        console       Trace text level\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   std::string optImageFileName;
 113   int opt_ndet;
 114   int opt_nview;
 115   int opt_offsetview = 0;
 116   int opt_nray = 1;
 117   double dOptFocalLength = 2.;
 118   double dOptCenterDetectorLength = 2.;
 119   double dOptViewRatio = 1.;
 120   double dOptScanRatio = 1.;
 121   int opt_trace = Trace::TRACE_NONE;
 122   int opt_verbose = 0;
 123   int opt_debug = 0;
 124   double opt_rotangle = -1;
 125   char* endptr = NULL;
 126   char* endstr;
 127   UNUSED(opt_debug);
 128
 129 #ifdef HAVE_MPI
 130   MPIWorld mpiWorld (argc, argv);
 131 #endif
 132
 133   Timer timerProgram;
 134
 135 #ifdef HAVE_MPI
 136   if (mpiWorld.getRank() == 0) {
 137 #endif
 138     while (1) {
 139       int c = getopt_long(argc, argv, "", phm2pj_options, NULL);
 140
 141       if (c == -1)
 142         break;
 143
 144       switch (c) {
 145       case O_PHANTOM:
 146         optPhmName = optarg;
 147         break;
 148       case O_PHMFILE:
 149         optPhmFileName = optarg;
 150         break;
 151       case O_IMAGEFILE:
 152         optImageFileName = optarg;
 153         break;
 154       case O_VERBOSE:
 155         opt_verbose = 1;
 156         break;
 157       case O_DEBUG:
 158         opt_debug = 1;
 159         break;
 160         break;
 161       case O_TRACE:
 162         if ((opt_trace = Trace::convertTraceNameToID(optarg)) == Trace::TRACE_INVALID) {
 163           phm2pj_usage(argv[0]);
 164           return (1);
 165         }
 166         break;
 167       case O_DESC:
 168         opt_desc = optarg;
 169         break;
 170       case O_ROTANGLE:
 171         opt_rotangle = strtod(optarg, &endptr);
 172         endstr = optarg + strlen(optarg);
 173         if (endptr != endstr) {
 174           std::cerr << "Error setting --rotangle to " << optarg << std::endl;
 175           phm2pj_usage(argv[0]);
 176           return (1);
 177         }
 178         break;
 179       case O_GEOMETRY:
 180         optGeometryName = optarg;
 181         break;
 182       case O_FOCAL_LENGTH:
 183         dOptFocalLength = strtod(optarg, &endptr);
 184         endstr = optarg + strlen(optarg);
 185         if (endptr != endstr) {
 186           std::cerr << "Error setting --focal-length to " << optarg << std::endl;
 187           phm2pj_usage(argv[0]);
 188           return (1);
 189         }
 190         break;
 191       case O_CENTER_DETECTOR_LENGTH:
 192         dOptCenterDetectorLength = strtod(optarg, &endptr);
 193         endstr = optarg + strlen(optarg);
 194         if (endptr != endstr) {
 195           std::cerr << "Error setting --center-detector-length to " << optarg << std::endl;
 196           phm2pj_usage(argv[0]);
 197           return (1);
 198         }
 199         break;
 200       case O_VIEW_RATIO:
 201         dOptViewRatio = strtod(optarg, &endptr);
 202         endstr = optarg + strlen(optarg);
 203         if (endptr != endstr) {
 204           std::cerr << "Error setting --view-ratio to " << optarg << std::endl;
 205           phm2pj_usage(argv[0]);
 206           return (1);
 207         }
 208         break;
 209       case O_SCAN_RATIO:
 210         dOptScanRatio = strtod(optarg, &endptr);
 211         endstr = optarg + strlen(optarg);
 212         if (endptr != endstr) {
 213           std::cerr << "Error setting --scan-ratio to " << optarg << std::endl;
 214           phm2pj_usage(argv[0]);
 215           return (1);
 216         }
 217         break;
 218       case O_NRAY:
 219         opt_nray = strtol(optarg, &endptr, 10);
 220         endstr = optarg + strlen(optarg);
 221         if (endptr != endstr) {
 222           std::cerr << "Error setting --nray to %s" << optarg << std::endl;
 223           phm2pj_usage(argv[0]);
 224           return (1);
 225         }
 226         break;
 227           case O_OFFSETVIEW:
 228                 opt_offsetview = strtol(optarg, &endptr, 10);
 229                 endstr = optarg + strlen(optarg);
 230                 if (endptr != endstr) {
 231                   std::cerr << "Error setting --offsetview to %s" << optarg << std::endl;
 232                   phm2pj_usage(argv[0]);
 233                   return (1);
 234                 }
 235                 break;
 236
 237       case O_VERSION:
 238 #ifdef VERSION
 239         std::cout << "Version: " << VERSION << std::endl << g_szIdStr << std::endl;
 240 #else
 241         std::cout << "Unknown version number\n";
 242 #endif
 243         return (0);
 244       case O_HELP:
 245       case '?':
 246         phm2pj_usage(argv[0]);
 247         return (0);
 248       default:
 249         phm2pj_usage(argv[0]);
 250         return (1);
 251       }
 252     }
 253
 254     if (optPhmName.empty() && optPhmFileName.empty() && optImageFileName.empty()) {
 255       std::cerr << "No phantom defined\n" << std::endl;
 256       phm2pj_usage(argv[0]);
 257       return (1);
 258     }
 259     if ((! optPhmName.empty() && !optPhmFileName.empty()) ||
 260         (! optPhmName.empty() && !optImageFileName.empty()) ||
 261         (! optPhmFileName.empty() && !optImageFileName.empty())) {
 262       std::cerr << "Multiple phantoms define. Use only --phantom, --phmfile, or --imagefile\n" << std::endl;
 263       phm2pj_usage(argv[0]);
 264       return (1);
 265     }
 266
 267     if (optind + 3 != argc) {
 268       phm2pj_usage(argv[0]);
 269       return (1);
 270     }
 271
 272     opt_outfile = argv[optind];
 273     opt_ndet = strtol(argv[optind+1], &endptr, 10);
 274     endstr = argv[optind+1] + strlen(argv[optind+1]);
 275     if (endptr != endstr) {
 276       std::cerr << "Error setting --ndet to " << argv[optind+1] << std::endl;
 277       phm2pj_usage(argv[0]);
 278       return (1);
 279     }
 280     opt_nview = strtol(argv[optind+2], &endptr, 10);
 281     endstr = argv[optind+2] + strlen(argv[optind+2]);
 282     if (endptr != endstr) {
 283       std::cerr << "Error setting --nview to " << argv[optind+2] << std::endl;
 284       phm2pj_usage(argv[0]);
 285       return (1);
 286     }
 287
 288     if (opt_rotangle < 0) {
 289       if (optGeometryName.compare ("parallel") == 0)
 290         opt_rotangle = 0.5;
 291       else
 292         opt_rotangle = 1.0;
 293     }
 294
 295     std::ostringstream desc;
 296     desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", OffsetView =" << opt_offsetview << ", Geometry=" << optGeometryName << ", ";
 297
 298     if (!optPhmFileName.empty()) {
 299       desc << "PhantomFile=" << optPhmFileName;
 300     } else if (!optPhmName.empty()) {
 301       desc << "Phantom=" << optPhmName;
 302     } else if (!optImageFileName.empty()) {
 303       desc << "Imagefile=" << optImageFileName;
 304     }
 305     if (opt_desc.length()) {
 306       desc << ": " << opt_desc;
 307     }
 308     opt_desc = desc.str();
 309
 310     if (!optPhmName.empty()) {
 311       phm.createFromPhantom (optPhmName.c_str());
 312       if (phm.fail()) {
 313         std::cout << phm.failMessage() << std::endl << std::endl;
 314         phm2pj_usage(argv[0]);
 315         return (1);
 316       }
 317     } else if (!optPhmFileName.empty()) {
 318 #ifdef HAVE_MPI
 319       std::cerr << "Can not read phantom from file in MPI mode\n";
 320       return (1);
 321 #endif
 322       phm.createFromPhmFile (optPhmFileName.c_str());
 323     } else if (!optImageFileName.empty()) {
 324 #ifdef HAVE_MPI
 325       std::cerr << "Can not read image from file in MPI mode\n";
 326       return (1);
 327 #endif
 328       phm.createFromImageFile (optImageFileName.c_str());
 329     }
 330
 331 #ifdef HAVE_MPI
 332   }
 333 #endif
 334
 335 #ifdef HAVE_MPI
 336   TimerCollectiveMPI timerBcast(mpiWorld.getComm());
 337   mpiWorld.BcastString (optPhmName);
 338   mpiWorld.getComm().Bcast (&opt_rotangle, 1, MPI::DOUBLE, 0);
 339   mpiWorld.getComm().Bcast (&dOptFocalLength, 1, MPI::DOUBLE, 0);
 340   mpiWorld.getComm().Bcast (&dOptCenterDetectorLength, 1, MPI::DOUBLE, 0);
 341   mpiWorld.getComm().Bcast (&dOptViewRatio, 1, MPI::DOUBLE, 0);
 342   mpiWorld.getComm().Bcast (&dOptScanRatio, 1, MPI::DOUBLE, 0);
 343   mpiWorld.getComm().Bcast (&opt_nview, 1, MPI::INT, 0);
 344   mpiWorld.getComm().Bcast (&opt_ndet, 1, MPI::INT, 0);
 345   mpiWorld.getComm().Bcast (&opt_nray, 1, MPI::INT, 0);
 346   mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0);
 347   mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
 348   mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
 349   if (opt_verbose)
 350     timerBcast.timerEndAndReport ("Time to broadcast variables");
 351
 352   if (mpiWorld.getRank() > 0 && optPhmName != "")
 353     phm.createFromPhantom (optPhmName.c_str());
 354 #endif
 355
 356   opt_rotangle *= TWOPI;
 357   Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_offsetview, opt_nray,
 358                 opt_rotangle, dOptFocalLength, dOptCenterDetectorLength, dOptViewRatio, dOptScanRatio);
 359   if (scanner.fail()) {
 360     std::cout << "Scanner Creation Error: " << scanner.failMessage() << std::endl;
 361     return (1);
 362   }
 363 #ifdef HAVE_MPI
 364   mpiWorld.setTotalWorkUnits (opt_nview);
 365
 366   Projections pjGlobal;
 367   if (mpiWorld.getRank() == 0)
 368     pjGlobal.initFromScanner (scanner);
 369
 370   if (opt_verbose) {
 371     std::ostringstream os;
 372     pjGlobal.printScanInfo(os);
 373     std::cout << os.str();
 374   }
 375
 376   Projections pjLocal (scanner);
 377   pjLocal.setNView (mpiWorld.getMyLocalWorkUnits());
 378
 379   if (opt_debug)
 380     std::cout << "pjLocal->nview = " << pjLocal.nView() << " (process " << mpiWorld.getRank() << ")\n";;
 381
 382   TimerCollectiveMPI timerProject (mpiWorld.getComm());
 383   scanner.collectProjections (pjLocal, phm, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), false, opt_trace);
 384   if (opt_verbose)
 385     timerProject.timerEndAndReport ("Time to collect projections");
 386
 387   TimerCollectiveMPI timerGather (mpiWorld.getComm());
 388   GatherProjectionsMPI (mpiWorld, pjGlobal, pjLocal, opt_debug);
 389   if (opt_verbose)
 390     timerGather.timerEndAndReport ("Time to gather projections");
 391
 392 #else
 393   Projections pjGlobal (scanner);
 394   scanner.collectProjections (pjGlobal, phm, 0, opt_nview, opt_offsetview, true, opt_trace);
 395 #endif
 396
 397 #ifdef HAVE_MPI
 398   if (mpiWorld.getRank() == 0)
 399 #endif
 400   {
 401     pjGlobal.setCalcTime (timerProgram.timerEnd());
 402     pjGlobal.setRemark (opt_desc);
 403     pjGlobal.write (opt_outfile);
 404     if (opt_verbose) {
 405       phm.print (std::cout);
 406       std::cout << std::endl;
 407       std::ostringstream os;
 408       pjGlobal.printScanInfo (os);
 409       std::cout << os.str() << std::endl;
 410       std::cout << "  Remark: " << pjGlobal.remark() << std::endl;
 411       std::cout << "Run time: " << pjGlobal.calcTime() << " seconds\n";
 412     }
 413   }
 414
 415   return (0);
 416 }
 417
 418
 419 /* FUNCTION
 420 *    GatherProjectionsMPI
 421 *
 422 * SYNOPSIS
 423 *    Gather's raysums from all processes in pjLocal in pjGlobal in process 0
 424 */
 425
 426 #ifdef HAVE_MPI
 427 void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug)
 428 {
 429   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
 430     DetectorArray& detArray = pjLocal.getDetectorArray(iw);
 431     double viewAngle = detArray.viewAngle();
 432     int nDet = detArray.nDet();
 433     DetectorValue* detval = detArray.detValues();
 434
 435     mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, 0, 0);
 436     mpiWorld.getComm().Send(&nDet, 1, MPI::INT, 0, 0);
 437     mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, 0, 0);
 438   }
 439
 440   if (mpiWorld.getRank() == 0) {
 441     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
 442       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
 443         MPI::Status status;
 444         double viewAngle;
 445         int nDet;
 446         DetectorArray& detArray = pjGlobal.getDetectorArray(iw);
 447         DetectorValue* detval = detArray.detValues();
 448
 449         mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, iProc, 0, status);
 450         mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, iProc, 0, status);
 451         mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, iProc, 0, status);
 452         detArray.setViewAngle (viewAngle);
 453       }
 454     }
 455   }
 456 }
 457 #endif
 458
 459
 460 #ifndef NO_MAIN
 461 int
 462 main (int argc, char* argv[])
 463 {
 464   int retval = 1;
 465
 466   try {
 467     retval = phm2pj_main(argc, argv);
 468 #if HAVE_DMALLOC
 469     //    dmalloc_shutdown();
 470 #endif
 471   } catch (exception e) {
 472     std::cerr << "Exception: " << e.what() << std::endl;
 473   } catch (...) {
 474     std::cerr << "Unknown exception\n";
 475   }
 476
 477   return (retval);
 478 }
 479 #endif
 480