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