src/phm2if.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          phm2if.cpp
   5 **   Purpose:       Convert an phantom object to an image file
   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: phm2if.cpp,v 1.8 2000/06/13 16:20:31 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_NSAMPLE, O_FILTER, O_TRACE, O_VERBOSE, O_HELP,
  33        O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
  34
  35 static struct option my_options[] =
  36 {
  37   {"phantom", 1, 0, O_PHANTOM},
  38   {"phmfile", 1, 0, O_PHMFILE},
  39   {"desc", 1, 0, O_DESC},
  40   {"nsample", 1, 0, O_NSAMPLE},
  41   {"filter", 1, 0, O_FILTER},
  42   {"filter-domain", 1, 0, O_FILTER_DOMAIN},
  43   {"filter-bw", 1, 0, O_FILTER_BW},
  44   {"filter-param", 1, 0, O_FILTER_PARAM},
  45   {"trace", 1, 0, O_TRACE},
  46   {"verbose", 0, 0, O_VERBOSE},
  47   {"debug", 0, 0, O_DEBUG},
  48   {"help", 0, 0, O_HELP},
  49   {"version", 0, 0, O_VERSION},
  50   {0, 0, 0, 0}
  51 };
  52
  53 void
  54 phm2if_usage (const char *program)
  55 {
  56   cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]" << endl;
  57   cout << "Generate phantom image from a predefined --phantom or a --phmfile" << endl;
  58   cout << endl;
  59   cout << "     outfile         Name of output file for image" << endl;
  60   cout << "     nx              Number of pixels X-axis" << endl;
  61   cout << "     ny              Number of pixels Y-axis" << endl;
  62   cout << "     --phantom       Phantom to use for projection" << endl;
  63   cout << "        herman       Herman head phantom" << endl;
  64   cout << "        rowland      Rowland head phantom" << endl;
  65   cout << "        browland     Bordered Rowland head phantom" << endl;
  66   cout << "        unitpulse    Unit pulse phantom" << endl;
  67   cout << "     --phmfile       Generate Phantom from phantom file" << endl;
  68   cout << "     --filter        Generate Phantom from a filter function" << endl;
  69   cout << "       abs_bandlimit Abs * Bandlimiting" << endl;
  70   cout << "       abs_sinc      Abs * Sinc" << endl;
  71   cout << "       abs_cos       Abs * Cosine" << endl;
  72   cout << "       abs_hamming   Abs * Hamming" << endl;
  73   cout << "       shepp         Shepp-Logan" << endl;
  74   cout << "       bandlimit     Bandlimiting" << endl;
  75   cout << "       sinc          Sinc" << endl;
  76   cout << "       cos           Cosine" << endl;
  77   cout << "       triangle      Triangle" << endl;
  78   cout << "       hamming       Hamming" << endl;
  79   cout << "     --filter-param  Alpha level for Hamming filter" << endl;
  80   cout << "     --filter-domain Set domain of filter" << endl;
  81   cout << "         spatial     Spatial domain (default)" << endl;
  82   cout << "         freq        Frequency domain" << endl;
  83   cout << "     --filter-bw     Filter bandwidth (default = 1)" << endl;
  84   cout << "     --desc          Description of raysum" << endl;
  85   cout << "     --nsample       Number of samples per axis per pixel (default = 1)" << endl;
  86   cout << "     --trace         Trace level to use" << endl;
  87   cout << "        none         No tracing (default)" << endl;
  88   cout << "        text         Trace text level" << endl;
  89   cout << "        phm          Trace phantom" << endl;
  90   cout << "        rays         Trace rays" << endl;
  91   cout << "        plot         Trace plot" << endl;
  92   cout << "        clipping     Trace clipping" << endl;
  93   cout << "     --debug         Debug mode" << endl;
  94   cout << "     --verbose       Verbose mode" << endl;
  95   cout << "     --version       Print version" << endl;
  96   cout << "     --help          Print this help message" << endl;
  97 }
  98
  99 #ifdef HAVE_MPI
 100 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug);
 101 #endif
 102
 103 int
 104 phm2if_main (int argc, char* argv[])
 105 {
 106   ImageFile* imGlobal = NULL;
 107   PHANTOM *phm = NULL;
 108   int opt_nx = 0, opt_ny = 0;
 109   int opt_nsample = 1;
 110   int opt_phmnum = -1;
 111   FilterType opt_filter = static_cast<FilterType>(-1);
 112   DomainType opt_filter_domain = D_SPATIAL;
 113   char *opt_outfile = NULL;
 114   int opt_debug = 0;
 115   char str[256];
 116   char opt_desc[256], opt_phmfilename[256];
 117   char *endstr, *endptr;
 118   double opt_filter_param = 1;
 119   double opt_filter_bw = 1.;
 120   int opt_trace = TRACE_NONE;
 121   int opt_verbose = 0;
 122 #ifdef HAVE_MPI
 123   ImageFile* imLocal = NULL;
 124   MPIWorld mpiWorld (argc, argv);
 125 #endif
 126
 127   Timer timerProgram;
 128
 129 #ifdef HAVE_MPI
 130   if (mpiWorld.getRank() == 0) {
 131 #endif
 132     strcpy(opt_desc, "");
 133     strcpy(opt_phmfilename, "");
 134     while (1) {
 135       int c = getopt_long(argc, argv, "", my_options, NULL);
 136       char *endptr = NULL;
 137       char *endstr;
 138
 139       if (c == -1)
 140         break;
 141
 142       switch (c) {
 143       case O_PHANTOM:
 144         if ((opt_phmnum = opt_set_phantom(optarg)) < 0) {
 145           phm2if_usage(argv[0]);
 146           return (1);
 147         }
 148         break;
 149       case O_PHMFILE:
 150         strncpy(opt_phmfilename, optarg, sizeof(opt_phmfilename));
 151         phm = phm_create_from_file(opt_phmfilename);
 152 #ifdef HAVE_MPI
 153         if (mpiWorld.getRank() == 0)
 154           cerr << "Can't use phantom from file in MPI mode" << endl;
 155         return (1);
 156 #endif
 157         break;
 158       case O_VERBOSE:
 159         opt_verbose = 1;
 160         break;
 161       case O_DEBUG:
 162         opt_debug = 1;
 163         break;
 164       case O_TRACE:
 165         if ((opt_trace = opt_set_trace(optarg)) < 0) {
 166           phm2if_usage(argv[0]);
 167           return (1);
 168         }
 169         break;
 170       case O_FILTER:
 171         if ((opt_filter = opt_set_filter(optarg)) < 0) {
 172           phm2if_usage (argv[0]);
 173           return (1);
 174         }
 175         break;
 176       case O_FILTER_DOMAIN:
 177         if ((opt_filter_domain = opt_set_filter_domain(optarg)) < 0) {
 178           phm2if_usage (argv[0]);
 179           return (1);
 180         }
 181         break;
 182       case O_DESC:
 183         strncpy(opt_desc, optarg, sizeof(opt_desc));
 184         break;
 185       case O_FILTER_BW:
 186         opt_filter_bw = strtod(optarg, &endptr);
 187         endstr = optarg + strlen(optarg);
 188         if (endptr != endstr) {
 189           fprintf(stderr,"Error setting --filter-bw to %s\n", optarg);
 190           phm2if_usage(argv[0]);
 191           return (1);
 192         }
 193         break;
 194       case O_FILTER_PARAM:
 195         opt_filter_param = strtod(optarg, &endptr);
 196         endstr = optarg + strlen(optarg);
 197         if (endptr != endstr) {
 198           fprintf(stderr,"Error setting --filter-param to %s\n", optarg);
 199           phm2if_usage(argv[0]);
 200           return (1);
 201         }
 202         break;
 203       case O_NSAMPLE:
 204         opt_nsample = strtol(optarg, &endptr, 10);
 205         endstr = optarg + strlen(optarg);
 206         if (endptr != endstr) {
 207           fprintf(stderr,"Error setting --nsample to %s\n", optarg);
 208           phm2if_usage(argv[0]);
 209           return (1);
 210         }
 211         break;
 212         case O_VERSION:
 213 #ifdef VERSION
 214           cout <<  "Version " << VERSION << endl;
 215 #else
 216           cerr << "Unknown version number" << endl;
 217 #endif
 218       case O_HELP:
 219       case '?':
 220         phm2if_usage(argv[0]);
 221         return (0);
 222       default:
 223         phm2if_usage(argv[0]);
 224         return (1);
 225       }
 226     }
 227
 228     if (phm == NULL && opt_phmnum == -1 && opt_filter == -1) {
 229       cerr << "No phantom defined" << endl;
 230       phm2if_usage(argv[0]);
 231       return (1);
 232     }
 233
 234     if (optind + 3 != argc) {
 235       phm2if_usage(argv[0]);
 236       return (1);
 237     }
 238     opt_outfile = argv[optind];
 239     opt_nx = strtol(argv[optind+1], &endptr, 10);
 240     endstr = argv[optind+1] + strlen(argv[optind+1]);
 241     if (endptr != endstr) {
 242       fprintf(stderr,"Error setting nx to %s\n", argv[optind+1]);
 243       phm2if_usage(argv[0]);
 244       return (1);
 245     }
 246     opt_ny = strtol(argv[optind+2], &endptr, 10);
 247     endstr = argv[optind+2] + strlen(argv[optind+2]);
 248     if (endptr != endstr) {
 249       fprintf(stderr,"Error setting ny to %s\n", argv[optind+2]);
 250       phm2if_usage(argv[0]);
 251       return (1);
 252     }
 253
 254     snprintf(str, sizeof(str), "nx=%d, ny=%d, nsample=%d, ", opt_nx, opt_ny, opt_nsample);
 255     if (opt_phmfilename[0]) {
 256       strncat(str, "phantom=", sizeof(str));
 257       strncat(str, opt_phmfilename, sizeof(str));
 258     }
 259     else if (opt_phmnum != -1) {
 260       strncat(str, "phantom=", sizeof(str));
 261       strncat(str, name_of_phantom(opt_phmnum), sizeof(str));
 262     }
 263     else if (opt_filter != -1) {
 264       strncat(str, "filter=", sizeof(str));
 265       strncat(str, name_of_filter(opt_filter), sizeof(str));
 266       strncat(str, " - ", sizeof(str));
 267       strncat(str, name_of_filter_domain(opt_filter_domain), sizeof(str));
 268     }
 269     if (opt_desc[0]) {
 270       strncat(str, ": ", sizeof(str));
 271       strncat(str, opt_desc, sizeof(str));
 272     }
 273     strncpy(opt_desc, str, sizeof(opt_desc));
 274
 275     if (opt_verbose)
 276       cout << "Rasterize Phantom to Image" << endl << endl;
 277 #ifdef HAVE_MPI
 278   }
 279 #endif
 280
 281 #ifdef HAVE_MPI
 282   TimerMPI timerBcast (mpiWorld.getComm());
 283   mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0);
 284   mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
 285   mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
 286   mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
 287   mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
 288   mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
 289   mpiWorld.getComm().Bcast (&opt_phmnum, 1, MPI::INT, 0);
 290   mpiWorld.getComm().Bcast (&opt_filter, 1, MPI::INT, 0);
 291   mpiWorld.getComm().Bcast (&opt_filter_domain, 1, MPI::INT, 0);
 292   mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
 293   mpiWorld.getComm().Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
 294
 295   if (opt_verbose)
 296     timerBcast.timerEndAndReport ("Time to broadcast variables");
 297
 298   mpiWorld.setTotalWorkUnits (opt_nx);
 299
 300   if (mpiWorld.getRank() == 0) {
 301     imGlobal = new ImageFile (opt_outfile, opt_nx, opt_ny);
 302     imGlobal->fileCreate();
 303   }
 304   imLocal = new ImageFile (opt_nx, opt_ny);
 305 #else
 306   imGlobal = new ImageFile (opt_outfile, opt_nx, opt_ny);
 307   imGlobal->fileCreate ();
 308 #endif
 309
 310   if (opt_phmnum >= 0)
 311       phm = phm_create (opt_phmnum);
 312
 313 #ifdef HAVE_MPI
 314   else {
 315     if (mpiWorld.getRank() == 0)
 316       cerr << "phmnum < 0" << endl;
 317     exit(1);
 318   }
 319 #endif
 320
 321   ImageFileArray v;
 322 #ifdef HAVE_MPI
 323   if (mpiWorld.getRank() == 0)
 324     v = imGlobal->getArray ();
 325
 326   if (phm->type == P_UNIT_PULSE) {
 327     if (mpiWorld.getRank() == 0) {
 328       v[opt_nx/2][opt_ny/2] = 1.;
 329     }
 330   } else if (opt_filter != -1) {
 331     if (mpiWorld.getRank() == 0) {
 332       image_filter_response (*imGlobal, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 333     }
 334   } else {
 335     TimerMPI timerRasterize (mpiWorld.getComm());
 336     phm_to_imagefile (phm, *imLocal, mpiWorld.getMyStartWorkUnit(),  mpiWorld.getMyLocalWorkUnits(), opt_nsample, opt_trace);
 337     if (opt_verbose)
 338       timerRasterize.timerEndAndReport ("Time to rasterize phantom");
 339
 340     TimerMPI timerGather (mpiWorld.getComm());
 341     mpi_gather_image (mpiWorld, imGlobal, imLocal, opt_debug);
 342     if (opt_verbose)
 343       timerGather.timerEndAndReport ("Time to gather image");
 344   }
 345 #else
 346   v = imGlobal->getArray ();
 347   if (phm->type == P_UNIT_PULSE) {
 348     v[opt_nx/2][opt_ny/2] = 1.;
 349   } else if (opt_filter != -1) {
 350     image_filter_response (*imGlobal, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 351   } else {
 352 #if HAVE_SGP
 353       if (opt_trace >= TRACE_PHM)
 354         phm_show(phm);
 355 #endif
 356       phm_to_imagefile (phm, *imGlobal, 0, opt_nx, opt_nsample, opt_trace);
 357   }
 358 #endif
 359
 360 #ifdef HAVE_MPI
 361   if (mpiWorld.getRank() == 0)
 362 #endif
 363   {
 364     imGlobal->arrayDataWrite ();
 365     double calctime = timerProgram.timerEnd ();
 366     imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, opt_desc, calctime);
 367     imGlobal->fileClose ();
 368     if (opt_verbose)
 369       cout << "Time to rasterized phantom: " << calctime << " seconds" << endl;
 370
 371     if (opt_trace >= TRACE_PHM) {
 372       double dmin, dmax;
 373       int nscale;
 374
 375       printf ("Enter display size scale (nominal = 1): ");
 376       scanf ("%d", &nscale);
 377       printf ("Enter minimum and maximum densities (min, max): ");
 378       scanf ("%lf %lf", &dmin, &dmax);
 379       image_display_scale (*imGlobal, nscale, dmin, dmax);
 380     }
 381   }
 382
 383   phm_free (phm);
 384
 385   return (0);
 386 }
 387
 388
 389
 390 #ifdef HAVE_MPI
 391 void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug)
 392 {
 393   ImageFileArray vLocal = imLocal->getArray();
 394   ImageFileArray vGlobal = NULL;
 395   int nyLocal = imLocal->ny();
 396
 397   if (mpiWorld.getRank() == 0)
 398     vGlobal = imGlobal->getArray();
 399
 400   for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
 401     mpiWorld.getComm().Send(vLocal[iw], nyLocal, imLocal->getMPIDataType(), 0, 0);
 402
 403   if (mpiWorld.getRank() == 0) {
 404     for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
 405       for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
 406         MPI::Status status;
 407         mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, imLocal->getMPIDataType(), iProc, 0, status);
 408       }
 409     }
 410   }
 411
 412 }
 413 #endif
 414
 415 #ifndef NO_MAIN
 416 int
 417 main (int argc, char* argv[])
 418 {
 419   return (phm2if_main(argc, argv));
 420 }
 421 #endif