tools/pjrec.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          pjrec.cpp
   5 **   Purpose:       Reconstruct an image from projections
   6 **   Programmer:    Kevin Rosenberg
   7 **   Date Started:  Aug 1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: pjrec.cpp,v 1.2 2000/07/15 08:36:13 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_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
  33
  34 static struct option my_options[] =
  35 {
  36   {"interp", 1, 0, O_INTERP},
  37   {"preinterpolation-factor", 1, 0, O_PREINTERPOLATION_FACTOR},
  38   {"filter", 1, 0, O_FILTER},
  39   {"filter-method", 1, 0, O_FILTER_METHOD},
  40   {"zeropad", 1, 0, O_ZEROPAD},
  41   {"filter-param", 1, 0, O_FILTER_PARAM},
  42   {"backproj", 1, 0, O_BACKPROJ},
  43   {"trace", 1, 0, O_TRACE},
  44   {"debug", 0, 0, O_DEBUG},
  45   {"verbose", 0, 0, O_VERBOSE},
  46   {"help", 0, 0, O_HELP},
  47   {"version", 0, 0, O_VERSION},
  48   {0, 0, 0, 0}
  49 };
  50
  51
  52 void
  53 pjrec_usage (const char *program)
  54 {
  55   cout << "usage: " << fileBasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl;
  56   cout << "Image reconstruction from raysum projections" << endl;
  57   cout << endl;
  58   cout << "  raysum-file     Input raysum file" << endl;
  59   cout << "  image-file      Output image file in SDF2D format" << endl;
  60   cout << "  nx-image        Number of columns in output image" << endl;
  61   cout << "  ny-image        Number of rows in output image" << endl;
  62   cout << "  --interp        Interpolation method during backprojection" << endl;
  63   cout << "    nearest     Nearest neighbor interpolation" << endl;
  64   cout << "    linear      Linear interpolation" << endl;
  65 #if HAVE_BSPLINE_INTERP
  66   cout << "    bspline     B-spline interpolation" << endl;
  67 #endif
  68   cout << "  --preinterpolate  Preinterpolation factor (default = 1)\n";
  69   cout << "                    Used only with frequency-based filtering\n";
  70   cout << "  --filter       Filter name" << endl;
  71   cout << "    abs_bandlimit Abs * Bandlimiting (default)" << endl;
  72   cout << "    abs_sinc      Abs * Sinc" << endl;
  73   cout << "    abs_cos       Abs * Cosine" << endl;
  74   cout << "    abs_hamming   Abs * Hamming" << endl;
  75   cout << "    shepp         Shepp-Logan" << endl;
  76   cout << "    bandlimit     Bandlimiting" << endl;
  77   cout << "    sinc          Sinc" << endl;
  78   cout << "    cos           Cosine" << endl;
  79   cout << "    triangle      Triangle" << endl;
  80   cout << "    hamming       Hamming" << endl;
  81   cout << "  --filter-method  Filter method before backprojections\n";;
  82   cout << "    convolution      Spatial filtering (default)\n";
  83   cout << "    fourier          Frequency filtering with discete fourier\n";
  84   cout << "    fourier_table    Frequency filtering with table lookup fourier\n";
  85   cout << "    fft              Fast Fourier Transform\n";
  86 #if HAVE_FFTW
  87   cout << "    fftw             Fast Fourier Transform West library\n";
  88   cout << "    rfftw            Fast Fourier Transform West (real-mode) library\n";
  89 #endif
  90   cout << "  --zeropad n   Set zeropad level (default = 0)\n";
  91   cout << "                set n to number of powers to two to pad\n";
  92   cout << "  --backproj    Backprojection Method" << endl;
  93   cout << "    trig        Trigometric functions at every point" << endl;
  94   cout << "    table       Trigometric functions with precalculated table" << endl;
  95   cout << "    diff        Difference method" << endl;
  96   cout << "    diff2       Optimized difference method (default)" << endl;
  97   cout << "    idiff2      Optimized difference method with integer math" << endl;
  98   cout << "    idiff3      Highly-optimized difference method with integer math" << endl;
  99   cout << "  --filter-param Alpha level for Hamming filter" << endl;
 100   cout << "  --trace        Set tracing to level" << endl;
 101   cout << "     none      No tracing (default)" << endl;
 102   cout << "     text      Text level tracing" << endl;
 103   cout << "     phm       Trace phantom" << endl;
 104   cout << "     rays      Trace allrays" << endl;
 105   cout << "     plot      Trace plotting" << endl;
 106   cout << "     clipping  Trace clipping" << endl;
 107   cout << "  --verbose      Turn on verbose mode" << endl;
 108   cout << "  --debug        Turn on debug mode" << endl;
 109   cout << "  --version      Print version" << endl;
 110   cout << "  --help         Print this help message" << endl;
 111 }
 112
 113
 114 #ifdef HAVE_MPI
 115 static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug);
 116 static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal);
 117 #endif
 118
 119
 120 int
 121 pjrec_main (int argc, char * argv[])
 122 {
 123   Projections projGlobal;
 124   ImageFile* imGlobal = NULL;
 125   char* filenameProj = NULL;
 126   char* filenameImage = NULL;
 127   string remark;
 128   char *endptr;
 129   int optVerbose = 0;
 130   int optDebug = 0;
 131   int optZeroPad = 0;
 132   int optTrace = TRACE_NONE;
 133   double optFilterParam = -1;
 134   string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR;
 135   string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
 136   string optInterpName = Backprojector::INTERP_LINEAR_STR;
 137   string optBackprojName = Backprojector::BPROJ_IDIFF2_STR;
 138   int optPreinterpolationFactor = 1;
 139   int nx, ny;
 140 #ifdef HAVE_MPI
 141   ImageFile* imLocal;
 142   int mpi_nview, mpi_ndet;
 143   double mpi_detinc, mpi_rotinc, mpi_phmlen;
 144   MPIWorld mpiWorld (argc, argv);
 145 #endif
 146
 147   Timer timerProgram;
 148
 149 #ifdef HAVE_MPI
 150   if (mpiWorld.getRank() == 0) {
 151 #endif
 152     while (1) {
 153       int c = getopt_long(argc, argv, "", my_options, NULL);
 154       char *endptr = NULL;
 155
 156       if (c == -1)
 157         break;
 158
 159       switch (c)
 160         {
 161         case O_INTERP:
 162           optInterpName = optarg;
 163           break;
 164         case O_PREINTERPOLATION_FACTOR:
 165           optPreinterpolationFactor = strtol(optarg, &endptr, 10);
 166           if (endptr != optarg + strlen(optarg)) {
 167             pjrec_usage(argv[0]);
 168             return(1);
 169           }
 170           break;
 171         case O_FILTER:
 172           optFilterName = optarg;
 173           break;
 174         case O_FILTER_METHOD:
 175           optFilterMethodName = optarg;
 176           break;
 177         case O_BACKPROJ:
 178           optBackprojName = optarg;
 179           break;
 180         case O_FILTER_PARAM:
 181           optFilterParam = strtod(optarg, &endptr);
 182           if (endptr != optarg + strlen(optarg)) {
 183             pjrec_usage(argv[0]);
 184             return(1);
 185           }
 186           break;
 187         case O_ZEROPAD:
 188           optZeroPad = strtol(optarg, &endptr, 10);
 189           if (endptr != optarg + strlen(optarg)) {
 190             pjrec_usage(argv[0]);
 191             return(1);
 192           }
 193           break;
 194         case O_VERBOSE:
 195           optVerbose = 1;
 196           break;
 197         case O_DEBUG:
 198           optDebug = 1;
 199           break;
 200         case O_TRACE:
 201           if ((optTrace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
 202             pjrec_usage(argv[0]);
 203             return (1);
 204           }
 205           break;
 206         case O_VERSION:
 207 #ifdef VERSION
 208           cout <<  "Version " <<  VERSION << endl;
 209 #else
 210           cout << "Unknown version number" << endl;
 211 #endif
 212           return (0);
 213         case O_HELP:
 214         case '?':
 215           pjrec_usage(argv[0]);
 216           return (0);
 217         default:
 218           pjrec_usage(argv[0]);
 219           return (1);
 220         }
 221     }
 222
 223     if (optind + 4 != argc) {
 224       pjrec_usage(argv[0]);
 225       return (1);
 226     }
 227
 228     filenameProj = argv[optind];
 229
 230     filenameImage = argv[optind + 1];
 231
 232     nx = strtol(argv[optind + 2], &endptr, 10);
 233     ny = strtol(argv[optind + 3], &endptr, 10);
 234
 235     ostringstream filterDesc;
 236     if (optFilterParam >= 0)
 237       filterDesc << optFilterName << ": alpha=" << optFilterParam;
 238     else
 239       filterDesc << optFilterName;
 240
 241     ostringstream label;
 242     label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolation=" << optPreinterpolationFactor << ", " << optBackprojName;
 243     remark = label.str();
 244
 245     if (optVerbose)
 246       cout << "Remark: " << remark << endl;
 247 #ifdef HAVE_MPI
 248   }
 249 #endif
 250
 251 #ifdef HAVE_MPI
 252   if (mpiWorld.getRank() == 0) {
 253     projGlobal.read (filenameProj);
 254     if (optVerbose)
 255       projGlobal.printScanInfo();
 256
 257     mpi_ndet = projGlobal.nDet();
 258     mpi_nview = projGlobal.nView();
 259     mpi_detinc = projGlobal.detInc();
 260     mpi_phmlen = projGlobal.phmLen();
 261     mpi_rotinc = projGlobal.rotInc();
 262   }
 263
 264   TimerCollectiveMPI timerBcast (mpiWorld.getComm());
 265   mpiWorld.BcastString (optBackprojName);
 266   mpiWorld.BcastString (optFilterName);
 267   mpiWorld.BcastString (optFilterMethodName);
 268   mpiWorld.BcastString (optInterpName);
 269   mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
 270   mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
 271   mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
 272   mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
 273   mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0);
 274   mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 1, MPI::INT, 0);
 275   mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0);
 276   mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
 277   mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
 278   mpiWorld.getComm().Bcast (&mpi_phmlen, 1, MPI::DOUBLE, 0);
 279   mpiWorld.getComm().Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0);
 280   mpiWorld.getComm().Bcast (&nx, 1, MPI::INT, 0);
 281   mpiWorld.getComm().Bcast (&ny, 1, MPI::INT, 0);
 282   if (optVerbose)
 283       timerBcast.timerEndAndReport ("Time to broadcast variables");
 284
 285   mpiWorld.setTotalWorkUnits (mpi_nview);
 286
 287   Projections projLocal (mpiWorld.getMyLocalWorkUnits(), mpi_ndet);
 288   projLocal.setDetInc (mpi_detinc);
 289   projLocal.setPhmLen (mpi_phmlen);
 290   projLocal.setRotInc (mpi_rotinc);
 291
 292   TimerCollectiveMPI timerScatter (mpiWorld.getComm());
 293   ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug);
 294   if (optVerbose)
 295       timerScatter.timerEndAndReport ("Time to scatter projections");
 296
 297   if (mpiWorld.getRank() == 0) {
 298     imGlobal = new ImageFile (nx, ny);
 299   }
 300
 301   imLocal = new ImageFile (nx, ny);
 302 #else
 303   projGlobal.read (filenameProj);
 304   if (optVerbose)
 305     projGlobal.printScanInfo();
 306
 307   imGlobal = new ImageFile (nx, ny);
 308 #endif
 309
 310 #ifdef HAVE_MPI
 311   TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
 312   projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
 313   if (optVerbose)
 314       timerReconstruct.timerEndAndReport ("Time to reconstruct");
 315
 316   TimerCollectiveMPI timerReduce (mpiWorld.getComm());
 317   ReduceImageMPI (mpiWorld, imLocal, imGlobal);
 318   if (optVerbose)
 319       timerReduce.timerEndAndReport ("Time to reduce image");
 320 #else
 321   projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
 322 #endif
 323
 324 #ifdef HAVE_MPI
 325   if (mpiWorld.getRank() == 0)
 326 #endif
 327     {
 328       double calcTime = timerProgram.timerEnd();
 329       imGlobal->labelAdd (projGlobal.getLabel());
 330       imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime);
 331       imGlobal->fileWrite (filenameImage);
 332       if (optVerbose)
 333         cout << "Run time: " << calcTime << " seconds" << endl;
 334     }
 335 #ifdef HAVE_MPI
 336   MPI::Finalize();
 337 #endif
 338
 339   return (0);
 340 }
 341
 342
 343 //////////////////////////////////////////////////////////////////////////////////////
 344 // MPI Support Routines
 345 //
 346 //////////////////////////////////////////////////////////////////////////////////////
 347
 348 #ifdef HAVE_MPI
 349 static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug)
 350 {
 351   if (mpiWorld.getRank() == 0) {
 352     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
 353       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
 354         DetectorArray& detarray = projGlobal.getDetectorArray( iw );
 355         int nDet = detarray.nDet();
 356         DetectorValue* detval = detarray.detValues();
 357
 358         double viewAngle = detarray.viewAngle();
 359         mpiWorld.getComm().Send(&nDet, 1, MPI::INT, iProc, 0);
 360         mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, iProc, 0);
 361         mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, iProc, 0);
 362       }
 363     }
 364   }
 365
 366   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
 367     MPI::Status status;
 368     int nDet;
 369     double viewAngle;
 370     DetectorValue* detval = projLocal.getDetectorArray(iw).detValues();
 371
 372     mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, 0, 0, status);
 373     mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, 0, 0, status);
 374     mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, 0, 0, status);
 375     projLocal.getDetectorArray(iw).setViewAngle( viewAngle );
 376   }
 377 }
 378
 379 static void
 380 ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal)
 381 {
 382   ImageFileArray vLocal = imLocal->getArray();
 383
 384   for (int ix = 0; ix < imLocal->nx(); ix++) {
 385     void *recvbuf = NULL;
 386     if (mpiWorld.getRank() == 0) {
 387       ImageFileArray vGlobal = imGlobal->getArray();
 388       recvbuf = vGlobal[ix];
 389     }
 390     mpiWorld.getComm().Reduce (vLocal[ix], recvbuf, imLocal->ny(), imLocal->getMPIDataType(), MPI::SUM, 0);
 391   }
 392 }
 393
 394 #endif
 395
 396
 397 #ifndef NO_MAIN
 398 int
 399 main (int argc, char* argv[])
 400 {
 401   int retval = 1;
 402
 403   try {
 404     retval = pjrec_main(argc, argv);
 405   } catch (exception e) {
 406     cerr << "Exception: " << e.what() << endl;
 407   } catch (...) {
 408     cerr << "Unknown exception" << endl;
 409   }
 410
 411   return (retval);
 412 }
 413 #endif
 414