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.5 2000/06/08 16:43:10 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 (ImageFile* im_global, ImageFile* im_local, const int opt_debug);
 107 #endif
 108
 109 int
 110 phm2sdf_main (const int argc, char *const 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   ImageFile* im_local = NULL;
 131   int mpi_argc = argc;
 132   char **mpi_argv = (char **) argv;
 133   double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
 134
 135   MPI::Init (mpi_argc, mpi_argv);
 136   mpi_ct.comm = MPI::COMM_WORLD.Dup ();
 137   mpi_ct.nproc = mpi_ct.comm.Get_size();
 138   mpi_ct.my_rank = mpi_ct.comm.Get_rank();
 139
 140   if (mpi_ct.nproc > CT_MPI_MAX_PROCESS) {
 141     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
 142               mpi_ct.nproc, CT_MPI_MAX_PROCESS);
 143   }
 144 #endif
 145
 146 #ifdef HAVE_MPI
 147   time_start = MPI::Wtime();
 148 #else
 149   time_start = td_current_sec();
 150 #endif
 151
 152 #ifdef HAVE_MPI
 153   if (mpi_ct.my_rank == 0) {
 154 #endif
 155
 156     strcpy(opt_desc, "");
 157     strcpy(opt_phmfilename, "");
 158     while (1) {
 159       int c = getopt_long(argc, argv, "", my_options, NULL);
 160       char *endptr = NULL;
 161       char *endstr;
 162
 163       if (c == -1)
 164         break;
 165
 166       switch (c) {
 167       case O_PHANTOM:
 168         if ((opt_phmnum = opt_set_phantom(optarg)) < 0) {
 169           phm2sdf_usage(argv[0]);
 170           return (1);
 171         }
 172         break;
 173       case O_PHMFILE:
 174         strncpy(opt_phmfilename, optarg, sizeof(opt_phmfilename));
 175         phm = phm_create_from_file(opt_phmfilename);
 176 #ifdef HAVE_MPI
 177         if (mpi_ct.my_rank == 0)
 178           fprintf(stderr, "Can't use phantom from file in MPI mode\n");
 179         return (1);
 180 #endif
 181         break;
 182       case O_VERBOSE:
 183         opt_verbose = 1;
 184         break;
 185       case O_DEBUG:
 186         opt_debug = 1;
 187         break;
 188       case O_TRACE:
 189         if ((opt_trace = opt_set_trace(optarg)) < 0) {
 190           phm2sdf_usage(argv[0]);
 191           return (1);
 192         }
 193         break;
 194       case O_FILTER:
 195         if ((opt_filter = opt_set_filter(optarg)) < 0) {
 196           phm2sdf_usage (argv[0]);
 197           return (1);
 198         }
 199         break;
 200       case O_FILTER_DOMAIN:
 201         if ((opt_filter_domain = opt_set_filter_domain(optarg)) < 0) {
 202           phm2sdf_usage (argv[0]);
 203           return (1);
 204         }
 205         break;
 206       case O_DESC:
 207         strncpy(opt_desc, optarg, sizeof(opt_desc));
 208         break;
 209       case O_FILTER_BW:
 210         opt_filter_bw = strtod(optarg, &endptr);
 211         endstr = optarg + strlen(optarg);
 212         if (endptr != endstr) {
 213           fprintf(stderr,"Error setting --filter-bw to %s\n", optarg);
 214           phm2sdf_usage(argv[0]);
 215           return (1);
 216         }
 217         break;
 218       case O_FILTER_PARAM:
 219         opt_filter_param = strtod(optarg, &endptr);
 220         endstr = optarg + strlen(optarg);
 221         if (endptr != endstr) {
 222           fprintf(stderr,"Error setting --filter-param to %s\n", optarg);
 223           phm2sdf_usage(argv[0]);
 224           return (1);
 225         }
 226         break;
 227       case O_NSAMPLE:
 228         opt_nsample = strtol(optarg, &endptr, 10);
 229         endstr = optarg + strlen(optarg);
 230         if (endptr != endstr) {
 231           fprintf(stderr,"Error setting --nsample to %s\n", optarg);
 232           phm2sdf_usage(argv[0]);
 233           return (1);
 234         }
 235         break;
 236         case O_VERSION:
 237 #ifdef VERSION
 238           fprintf(stdout, "Version %s\n", VERSION);
 239 #else
 240           fprintf(stderr, "Unknown version number");
 241 #endif
 242       case O_HELP:
 243       case '?':
 244         phm2sdf_usage(argv[0]);
 245         return (0);
 246       default:
 247         phm2sdf_usage(argv[0]);
 248         return (1);
 249       }
 250     }
 251
 252     if (phm == NULL && opt_phmnum == -1 && opt_filter == -1) {
 253       fprintf(stderr, "No phantom defined\n");
 254       phm2sdf_usage(argv[0]);
 255       return (1);
 256     }
 257
 258     if (optind + 3 != argc) {
 259       phm2sdf_usage(argv[0]);
 260       return (1);
 261     }
 262     opt_outfile = argv[optind];
 263     opt_nx = strtol(argv[optind+1], &endptr, 10);
 264     endstr = argv[optind+1] + strlen(argv[optind+1]);
 265     if (endptr != endstr) {
 266       fprintf(stderr,"Error setting nx to %s\n", argv[optind+1]);
 267       phm2sdf_usage(argv[0]);
 268       return (1);
 269     }
 270     opt_ny = strtol(argv[optind+2], &endptr, 10);
 271     endstr = argv[optind+2] + strlen(argv[optind+2]);
 272     if (endptr != endstr) {
 273       fprintf(stderr,"Error setting ny to %s\n", argv[optind+2]);
 274       phm2sdf_usage(argv[0]);
 275       return (1);
 276     }
 277
 278     snprintf(str, sizeof(str), "nx=%d, ny=%d, nsample=%d, ", opt_nx, opt_ny, opt_nsample);
 279     if (opt_phmfilename[0]) {
 280       strncat(str, "phantom=", sizeof(str));
 281       strncat(str, opt_phmfilename, sizeof(str));
 282     }
 283     else if (opt_phmnum != -1) {
 284       strncat(str, "phantom=", sizeof(str));
 285       strncat(str, name_of_phantom(opt_phmnum), sizeof(str));
 286     }
 287     else if (opt_filter != -1) {
 288       strncat(str, "filter=", sizeof(str));
 289       strncat(str, name_of_filter(opt_filter), sizeof(str));
 290       strncat(str, " - ", sizeof(str));
 291       strncat(str, name_of_filter_domain(opt_filter_domain), sizeof(str));
 292     }
 293     if (opt_desc[0]) {
 294       strncat(str, ": ", sizeof(str));
 295       strncat(str, opt_desc, sizeof(str));
 296     }
 297     strncpy(opt_desc, str, sizeof(opt_desc));
 298
 299     if (opt_verbose)
 300       printf("Rasterize Phantom to Image\n\n");
 301 #ifdef HAVE_MPI
 302   }
 303 #endif
 304
 305 #ifdef HAVE_MPI
 306   mpi_t1 = MPI::Wtime();
 307   mpi_ct.comm.Bcast (&opt_verbose, 1, MPI::INT, 0);
 308   mpi_ct.comm.Bcast (&opt_debug, 1, MPI::INT, 0);
 309   mpi_ct.comm.Bcast (&opt_trace, 1, MPI::INT, 0);
 310   mpi_ct.comm.Bcast (&opt_nx, 1, MPI::INT, 0);
 311   mpi_ct.comm.Bcast (&opt_ny, 1, MPI::INT, 0);
 312   mpi_ct.comm.Bcast (&opt_nsample, 1, MPI::INT, 0);
 313   mpi_ct.comm.Bcast (&opt_phmnum, 1, MPI::INT, 0);
 314   mpi_ct.comm.Bcast (&opt_filter, 1, MPI::INT, 0);
 315   mpi_ct.comm.Bcast (&opt_filter_domain, 1, MPI::INT, 0);
 316   mpi_ct.comm.Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
 317   mpi_ct.comm.Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
 318
 319   if (opt_verbose) {
 320     mpi_t2 = MPI::Wtime();
 321     mpi_t = mpi_t2 - mpi_t1;
 322     mpi_ct.comm.Reduce (&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
 323     if (mpi_ct.my_rank == 0)
 324       printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 325   }
 326
 327   mpi_ct_calc_work_units(opt_nx);
 328
 329   if (mpi_ct.my_rank == 0) {
 330     im_global = new ImageFile (opt_outfile, opt_nx, opt_ny);
 331     im_global->adf.fileCreate();
 332   }
 333   im_local = new ImageFile (opt_nx, opt_ny);
 334 #else
 335   im_global = new ImageFile (opt_outfile, opt_nx, opt_ny);
 336   im_global->adf.fileCreate ();
 337 #endif
 338
 339   if (opt_phmnum >= 0)
 340       phm = phm_create (opt_phmnum);
 341
 342 #ifdef HAVE_MPI
 343   else {
 344     if (mpi_ct.my_rank == 0)
 345       fprintf(stderr, "phmnum < 0\n");
 346     exit(1);
 347   }
 348 #endif
 349
 350   double calctime = 0;
 351   ImageFileArray v;
 352
 353 #ifdef HAVE_MPI
 354   if (mpi_ct.my_rank == 0)
 355     v = im_global->getArray ();
 356
 357   if (phm->type == P_UNIT_PULSE) {
 358     if (mpi_ct.my_rank == 0) {
 359       v[opt_nx/2][opt_ny/2] = 1.;
 360       calctime = 0;
 361     }
 362   } else if (opt_filter != -1) {
 363     if (mpi_ct.my_rank == 0) {
 364       image_filter_response (*im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 365       calctime = 0;
 366     }
 367   } else {
 368     mpi_ct.comm.Barrier();
 369
 370     if (opt_verbose)
 371       mpi_t1 = MPI::Wtime();
 372     phm_to_imagefile (phm, *im_local, mpi_ct.start_work_unit[mpi_ct.my_rank],  mpi_ct.local_work_units[mpi_ct.my_rank], opt_nsample, opt_trace);
 373     if (opt_verbose) {
 374       mpi_t2 = MPI::Wtime();
 375       mpi_t = mpi_t2 - mpi_t1;
 376       mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
 377       if (mpi_ct.my_rank == 0)
 378         printf("Time to compile phm = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 379     }
 380     mpi_gather_image(im_global, im_local, opt_debug);
 381   }
 382 #else
 383   v = im_global->getArray ();
 384   if (phm->type == P_UNIT_PULSE) {
 385     v[opt_nx/2][opt_ny/2] = 1.;
 386     calctime = 0;
 387   } else if (opt_filter != -1) {
 388     image_filter_response (*im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 389     calctime = 0;
 390   } else {
 391 #if HAVE_SGP
 392       if (opt_trace >= TRACE_PHM)
 393         phm_show(phm);
 394 #endif
 395       phm_to_imagefile (phm, *im_global, 0, opt_nx, opt_nsample, opt_trace);
 396   }
 397 #endif
 398
 399 #ifdef HAVE_MPI
 400   time_end = MPI::Wtime();
 401 #else
 402   time_end = td_current_sec();
 403 #endif
 404
 405 #ifdef HAVE_MPI
 406   if (mpi_ct.my_rank == 0)
 407 #endif
 408   {
 409     im_global->adf.arrayDataWrite ();
 410     calctime = time_end - time_start;
 411     im_global->adf.labelAdd (Array2dFileLabel::L_HISTORY, opt_desc, calctime);
 412     im_global->adf.fileClose ();
 413     if (opt_verbose)
 414       fprintf (stdout, "Time to compile image = %.2f sec\n\n", calctime);
 415
 416     if (opt_trace >= TRACE_PHM) {
 417       double dmin, dmax;
 418       int nscale;
 419
 420       printf ("Enter display size scale (nominal = 1): ");
 421       scanf ("%d", &nscale);
 422       printf ("Enter minimum and maximum densities (min, max): ");
 423       scanf ("%lf %lf", &dmin, &dmax);
 424       image_display_scale (*im_global, nscale, dmin, dmax);
 425     }
 426   }
 427
 428   phm_free (phm);
 429
 430   return (0);
 431 }
 432
 433
 434
 435 #ifdef HAVE_MPI
 436 void mpi_gather_image (ImageFile* im_global, ImageFile* im_local, const int opt_debug)
 437 {
 438   ImageFileArray vLocal = im_local->getArray();
 439   ImageFileArray vGlobal = NULL;
 440   if (mpi_ct.my_rank == 0)
 441     vGlobal = im_global->getArray();
 442   int nyLocal = im_local->ny();
 443
 444   int end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
 445   for (int iw = 0; iw <= end_work_unit; iw++) {
 446     mpi_ct.comm.Send(vLocal[iw], nyLocal, im_local->getMPIDataType(), 0, 0);
 447   }
 448
 449   if (mpi_ct.my_rank == 0) {
 450     for (int iproc = 0; iproc < mpi_ct.nproc; iproc++) {
 451       end_work_unit = mpi_ct.start_work_unit[iproc] + mpi_ct.local_work_units[iproc] - 1;
 452
 453       if (opt_debug) {
 454         fprintf(stdout, "Recv rs data from process %d (%d-%d)\n", iproc, mpi_ct.start_work_unit[iproc], end_work_unit);
 455         fflush(stdout);
 456       }
 457
 458       for (int iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
 459         MPI::Status status;
 460         mpi_ct.comm.Recv(vGlobal[iw], nyLocal, im_local->getMPIDataType(), iproc, 0, status);
 461       }
 462     }
 463   }
 464
 465 }
 466 #endif
 467
 468 #ifndef NO_MAIN
 469 int
 470 main (const int argc, char *const argv[])
 471 {
 472   return (phm2sdf_main(argc, argv));
 473 }
 474 #endif