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