src/phm2if.cpp

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (C) 1983-2000 Kevin Rosenberg
   4 **
   5 **  $Id: phm2if.cpp,v 1.4 2000/06/07 10:12:05 kevin Exp $
   6 **  $Log: phm2if.cpp,v $
   7 **  Revision 1.4  2000/06/07 10:12:05  kevin
   8 **  Upgraded from MPI to MPI++
   9 **
  10 **  Revision 1.3  2000/06/07 07:43:19  kevin
  11 **  Converted to IF data files and C++
  12 **
  13 **  Revision 1.2  2000/06/07 03:50:27  kevin
  14 **  *** empty log message ***
  15 **
  16 **  Revision 1.1  2000/06/07 02:29:05  kevin
  17 **  Initial C++ versions
  18 **
  19 **  Revision 1.10  2000/05/24 22:50:04  kevin
  20 **  Added support for new SGP library
  21 **
  22 **  Revision 1.9  2000/05/16 04:33:59  kevin
  23 **  Improved option processing
  24 **
  25 **  Revision 1.8  2000/05/11 01:06:30  kevin
  26 **  Changed sprintf to snprintf
  27 **
  28 **  Revision 1.7  2000/05/08 20:02:32  kevin
  29 **  ANSI C changes
  30 **
  31 **  Revision 1.6  2000/05/03 08:49:50  kevin
  32 **  Code cleanup
  33 **
  34 **  Revision 1.5  2000/04/30 19:17:54  kevin
  35 **  *** empty log message ***
  36 **
  37 **  Revision 1.4  2000/04/30 04:06:13  kevin
  38 **  Update Raysum i/o routines
  39 **  Fix MPI bug in ctrec (scatter_raysum) that referenced rs_global
  40 **
  41 **  Revision 1.3  2000/04/28 18:19:01  kevin
  42 **  removed unused files
  43 **
  44 **  Revision 1.2  2000/04/28 13:50:45  kevin
  45 **  Removed Makefile Makefile.in that are automatically generated by autoconf
  46 **
  47 **  Revision 1.1.1.1  2000/04/28 13:02:44  kevin
  48 **  Initial CVS import for first public release
  49 **
  50 **
  51 **
  52 **  This program is free software; you can redistribute it and/or modify
  53 **  it under the terms of the GNU General Public License (version 2) as
  54 **  published by the Free Software Foundation.
  55 **
  56 **  This program is distributed in the hope that it will be useful,
  57 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  58 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  59 **  GNU General Public License for more details.
  60 **
  61 **  You should have received a copy of the GNU General Public License
  62 **  along with this program; if not, write to the Free Software
  63 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  64 ******************************************************************************/
  65
  66 /* FILE
  67  *   phm2sdf.c                  Generate a SDF image from a phantom
  68  *
  69  * HISTORY
  70  *   1984 - Final DOS verion
  71  *   1999 - First UNIX version
  72  */
  73
  74 #include "ct.h"
  75
  76 enum { O_PHANTOM, O_DESC, O_NSAMPLE, O_FILTER, O_TRACE, O_VERBOSE, O_HELP,
  77        O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
  78
  79 static struct option my_options[] =
  80 {
  81   {"phantom", 1, 0, O_PHANTOM},
  82   {"phmfile", 1, 0, O_PHMFILE},
  83   {"desc", 1, 0, O_DESC},
  84   {"nsample", 1, 0, O_NSAMPLE},
  85   {"filter", 1, 0, O_FILTER},
  86   {"filter-domain", 1, 0, O_FILTER_DOMAIN},
  87   {"filter-bw", 1, 0, O_FILTER_BW},
  88   {"filter-param", 1, 0, O_FILTER_PARAM},
  89   {"trace", 1, 0, O_TRACE},
  90   {"verbose", 0, 0, O_VERBOSE},
  91   {"debug", 0, 0, O_DEBUG},
  92   {"help", 0, 0, O_HELP},
  93   {"version", 0, 0, O_VERSION},
  94   {0, 0, 0, 0}
  95 };
  96
  97 void
  98 phm2sdf_usage (const char *program)
  99 {
 100   fprintf(stdout,"phm2sdf_usage: %s outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]\n",kbasename(program));
 101   fprintf(stdout,"Generate phantom image from a predefined --phantom or a --phmfile\n");
 102   fprintf(stdout,"\n");
 103   fprintf(stdout,"     outfile         Name of output file for image\n");
 104   fprintf(stdout,"     nx              Number of pixels X-axis\n");
 105   fprintf(stdout,"     ny              Number of pixels Y-axis\n");
 106   fprintf(stdout,"     --phantom       Phantom to use for projection\n");
 107   fprintf(stdout,"        herman       Herman head phantom\n");
 108   fprintf(stdout,"        rowland      Rowland head phantom\n");
 109   fprintf(stdout,"        browland     Bordered Rowland head phantom\n");
 110   fprintf(stdout,"        unitpulse    Unit pulse phantom\n");
 111   fprintf(stdout,"     --phmfile       Generate Phantom from phantom file\n");
 112   fprintf(stdout,"     --filter        Generate Phantom from a filter function\n");
 113   fprintf(stdout,"       abs_bandlimit Abs * Bandlimiting\n");
 114   fprintf(stdout,"       abs_sinc      Abs * Sinc\n");
 115   fprintf(stdout,"       abs_cos       Abs * Cosine\n");
 116   fprintf(stdout,"       abs_hamming   Abs * Hamming\n");
 117   fprintf(stdout,"       shepp         Shepp-Logan\n");
 118   fprintf(stdout,"       bandlimit     Bandlimiting\n");
 119   fprintf(stdout,"       sinc          Sinc\n");
 120   fprintf(stdout,"       cos           Cosine\n");
 121   fprintf(stdout,"       triangle      Triangle\n");
 122   fprintf(stdout,"       hamming       Hamming\n");
 123   fprintf(stdout,"     --filter-param  Alpha level for Hamming filter\n");
 124   fprintf(stdout,"     --filter-domain Set domain of filter\n");
 125   fprintf(stdout,"         spatial     Spatial domain (default)\n");
 126   fprintf(stdout,"         freq        Frequency domain\n");
 127   fprintf(stdout,"     --filter-bw     Filter bandwidth (default = 1)\n");
 128   fprintf(stdout,"     --desc          Description of raysum\n");
 129   fprintf(stdout,"     --nsample       Number of samples per axis per pixel (default = 1)\n");
 130   fprintf(stdout,"     --trace         Trace level to use\n");
 131   fprintf(stdout,"        none         No tracing (default)\n");
 132   fprintf(stdout,"        text         Trace text level\n");
 133   fprintf(stdout,"        phm          Trace phantom\n");
 134   fprintf(stdout,"        rays         Trace rays\n");
 135   fprintf(stdout,"        plot         Trace plot\n");
 136   fprintf(stdout,"        clipping     Trace clipping\n");
 137   fprintf(stdout,"     --debug         Debug mode\n");
 138   fprintf(stdout,"     --verbose       Verbose mode\n");
 139   fprintf(stdout,"     --version       Print version\n");
 140   fprintf(stdout,"     --help          Print this help message\n");
 141 }
 142
 143 #ifdef HAVE_MPI
 144 void mpi_gather_image (ImageFile& im_global, ImageFile& im_local, const int opt_debug);
 145 #endif
 146
 147 int
 148 phm2sdf_main (const int argc, char *const argv[])
 149 {
 150   ImageFile* im_global = NULL;
 151   PHANTOM *phm = NULL;
 152   int opt_nx = 0, opt_ny = 0;
 153   int opt_nsample = 1;
 154   int opt_phmnum = -1;
 155   FilterType opt_filter = static_cast<FilterType>(-1);
 156   DomainType opt_filter_domain = D_SPATIAL;
 157   char *opt_outfile = NULL;
 158   int opt_debug = 0;
 159   char str[256];
 160   char opt_desc[256], opt_phmfilename[256];
 161   char *endstr, *endptr;
 162   double opt_filter_param = 1;
 163   double opt_filter_bw = 1.;
 164   int opt_trace = TRACE_NONE;
 165   int opt_verbose = 0;
 166   double time_start=0, time_end=0;
 167 #ifdef HAVE_MPI
 168   ImageFile* im_local = NULL;
 169   int mpi_argc = argc;
 170   char **mpi_argv = (char **) argv;
 171   double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
 172
 173   MPI::Init (mpi_argc, mpi_argv);
 174   mpi_ct.comm = MPI::COMM_WORLD.Dup ();
 175   mpi_ct.nproc = mpi_ct.comm.Get_size();
 176   mpi_ct.my_rank = mpi_ct.comm.Get_rank();
 177
 178   if (mpi_ct.nproc > CT_MPI_MAX_PROCESS) {
 179     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
 180               mpi_ct.nproc, CT_MPI_MAX_PROCESS);
 181   }
 182 #endif
 183
 184 #ifdef HAVE_MPI
 185   time_start = MPI::Wtime();
 186 #else
 187   time_start = td_current_sec();
 188 #endif
 189
 190 #ifdef HAVE_MPI
 191   if (mpi_ct.my_rank == 0) {
 192 #endif
 193
 194     strcpy(opt_desc, "");
 195     strcpy(opt_phmfilename, "");
 196     while (1) {
 197       int c = getopt_long(argc, argv, "", my_options, NULL);
 198       char *endptr = NULL;
 199       char *endstr;
 200
 201       if (c == -1)
 202         break;
 203
 204       switch (c) {
 205       case O_PHANTOM:
 206         if ((opt_phmnum = opt_set_phantom(optarg)) < 0) {
 207           phm2sdf_usage(argv[0]);
 208           return (1);
 209         }
 210         break;
 211       case O_PHMFILE:
 212         strncpy(opt_phmfilename, optarg, sizeof(opt_phmfilename));
 213         phm = phm_create_from_file(opt_phmfilename);
 214 #ifdef HAVE_MPI
 215         if (mpi_ct.my_rank == 0)
 216           fprintf(stderr, "Can't use phantom from file in MPI mode\n");
 217         return (1);
 218 #endif
 219         break;
 220       case O_VERBOSE:
 221         opt_verbose = 1;
 222         break;
 223       case O_DEBUG:
 224         opt_debug = 1;
 225         break;
 226       case O_TRACE:
 227         if ((opt_trace = opt_set_trace(optarg)) < 0) {
 228           phm2sdf_usage(argv[0]);
 229           return (1);
 230         }
 231         break;
 232       case O_FILTER:
 233         if ((opt_filter = opt_set_filter(optarg)) < 0) {
 234           phm2sdf_usage (argv[0]);
 235           return (1);
 236         }
 237         break;
 238       case O_FILTER_DOMAIN:
 239         if ((opt_filter_domain = opt_set_filter_domain(optarg)) < 0) {
 240           phm2sdf_usage (argv[0]);
 241           return (1);
 242         }
 243         break;
 244       case O_DESC:
 245         strncpy(opt_desc, optarg, sizeof(opt_desc));
 246         break;
 247       case O_FILTER_BW:
 248         opt_filter_bw = strtod(optarg, &endptr);
 249         endstr = optarg + strlen(optarg);
 250         if (endptr != endstr) {
 251           fprintf(stderr,"Error setting --filter-bw to %s\n", optarg);
 252           phm2sdf_usage(argv[0]);
 253           return (1);
 254         }
 255         break;
 256       case O_FILTER_PARAM:
 257         opt_filter_param = strtod(optarg, &endptr);
 258         endstr = optarg + strlen(optarg);
 259         if (endptr != endstr) {
 260           fprintf(stderr,"Error setting --filter-param to %s\n", optarg);
 261           phm2sdf_usage(argv[0]);
 262           return (1);
 263         }
 264         break;
 265       case O_NSAMPLE:
 266         opt_nsample = strtol(optarg, &endptr, 10);
 267         endstr = optarg + strlen(optarg);
 268         if (endptr != endstr) {
 269           fprintf(stderr,"Error setting --nsample to %s\n", optarg);
 270           phm2sdf_usage(argv[0]);
 271           return (1);
 272         }
 273         break;
 274         case O_VERSION:
 275 #ifdef VERSION
 276           fprintf(stdout, "Version %s\n", VERSION);
 277 #else
 278           fprintf(stderr, "Unknown version number");
 279 #endif
 280       case O_HELP:
 281       case '?':
 282         phm2sdf_usage(argv[0]);
 283         return (0);
 284       default:
 285         phm2sdf_usage(argv[0]);
 286         return (1);
 287       }
 288     }
 289
 290     if (phm == NULL && opt_phmnum == -1 && opt_filter == -1) {
 291       fprintf(stderr, "No phantom defined\n");
 292       phm2sdf_usage(argv[0]);
 293       return (1);
 294     }
 295
 296     if (optind + 3 != argc) {
 297       phm2sdf_usage(argv[0]);
 298       return (1);
 299     }
 300     opt_outfile = argv[optind];
 301     opt_nx = strtol(argv[optind+1], &endptr, 10);
 302     endstr = argv[optind+1] + strlen(argv[optind+1]);
 303     if (endptr != endstr) {
 304       fprintf(stderr,"Error setting nx to %s\n", argv[optind+1]);
 305       phm2sdf_usage(argv[0]);
 306       return (1);
 307     }
 308     opt_ny = strtol(argv[optind+2], &endptr, 10);
 309     endstr = argv[optind+2] + strlen(argv[optind+2]);
 310     if (endptr != endstr) {
 311       fprintf(stderr,"Error setting ny to %s\n", argv[optind+2]);
 312       phm2sdf_usage(argv[0]);
 313       return (1);
 314     }
 315
 316     snprintf(str, sizeof(str), "nx=%d, ny=%d, nsample=%d, ", opt_nx, opt_ny, opt_nsample);
 317     if (opt_phmfilename[0]) {
 318       strncat(str, "phantom=", sizeof(str));
 319       strncat(str, opt_phmfilename, sizeof(str));
 320     }
 321     else if (opt_phmnum != -1) {
 322       strncat(str, "phantom=", sizeof(str));
 323       strncat(str, name_of_phantom(opt_phmnum), sizeof(str));
 324     }
 325     else if (opt_filter != -1) {
 326       strncat(str, "filter=", sizeof(str));
 327       strncat(str, name_of_filter(opt_filter), sizeof(str));
 328       strncat(str, " - ", sizeof(str));
 329       strncat(str, name_of_filter_domain(opt_filter_domain), sizeof(str));
 330     }
 331     if (opt_desc[0]) {
 332       strncat(str, ": ", sizeof(str));
 333       strncat(str, opt_desc, sizeof(str));
 334     }
 335     strncpy(opt_desc, str, sizeof(opt_desc));
 336
 337     if (opt_verbose)
 338       printf("Rasterize Phantom to Image\n\n");
 339 #ifdef HAVE_MPI
 340   }
 341 #endif
 342
 343 #ifdef HAVE_MPI
 344   mpi_t1 = MPI::Wtime();
 345   mpi_ct.comm.Bcast (&opt_verbose, 1, MPI::INT, 0);
 346   mpi_ct.comm.Bcast (&opt_debug, 1, MPI::INT, 0);
 347   mpi_ct.comm.Bcast (&opt_trace, 1, MPI::INT, 0);
 348   mpi_ct.comm.Bcast (&opt_nx, 1, MPI::INT, 0);
 349   mpi_ct.comm.Bcast (&opt_ny, 1, MPI::INT, 0);
 350   mpi_ct.comm.Bcast (&opt_nsample, 1, MPI::INT, 0);
 351   mpi_ct.comm.Bcast (&opt_phmnum, 1, MPI::INT, 0);
 352   mpi_ct.comm.Bcast (&opt_filter, 1, MPI::INT, 0);
 353   mpi_ct.comm.Bcast (&opt_filter_domain, 1, MPI::INT, 0);
 354   mpi_ct.comm.Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
 355   mpi_ct.comm.Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
 356
 357   if (opt_verbose) {
 358     mpi_t2 = MPI::Wtime();
 359     mpi_t = mpi_t2 - mpi_t1;
 360     mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
 361     if (mpi_ct.my_rank == 0)
 362       printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 363   }
 364
 365   mpi_ct_calc_work_units(opt_nx);
 366
 367   if (mpi_ct.my_rank == 0) {
 368     im_global = new ImageFile (opt_outfile, opt_nx, opt_ny);
 369     im_global->adf.fileCreate();
 370   }
 371
 372   im_local = new ImageFile (opt_nx, opt_ny);
 373 #else
 374   im_global = new ImageFile (opt_outfile, opt_nx, opt_ny);
 375   im_global->adf.fileCreate ();
 376 #endif
 377
 378   if (opt_phmnum >= 0)
 379       phm = phm_create (opt_phmnum);
 380
 381 #ifdef HAVE_MPI
 382   else {
 383     if (mpi_ct.my_rank == 0)
 384       fprintf(stderr, "phmnum < 0\n");
 385     exit(1);
 386   }
 387 #endif
 388
 389   ImageFileArray v = im_global->getArray ();
 390   double calctime = 0;
 391
 392 #ifdef HAVE_MPI
 393   if (phm->type == P_UNIT_PULSE) {
 394     if (mpi_ct.my_rank == 0) {
 395       v[opt_nx/2][opt_ny/2] = 1.;
 396       calctime = 0;
 397     }
 398   } else if (opt_filter != -1) {
 399     if (mpi_ct.my_rank == 0) {
 400       image_filter_response (*im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 401       calctime = 0;
 402     }
 403   } else {
 404     if (opt_verbose)
 405       mpi_t1 = MPI::Wtime();
 406     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);
 407     if (opt_verbose) {
 408       mpi_t2 = MPI::Wtime();
 409       mpi_t = mpi_t2 - mpi_t1;
 410       mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
 411       if (mpi_ct.my_rank == 0)
 412         printf("Time to compile phm = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 413     }
 414     mpi_gather_image(*im_global, *im_local, opt_debug);
 415   }
 416 #else
 417   if (phm->type == P_UNIT_PULSE) {
 418     v[opt_nx/2][opt_ny/2] = 1.;
 419     calctime = 0;
 420   } else if (opt_filter != -1) {
 421     image_filter_response (*im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 422     calctime = 0;
 423   } else {
 424 #if HAVE_SGP
 425       if (opt_trace >= TRACE_PHM)
 426         phm_show(phm);
 427 #endif
 428       phm_to_imagefile (phm, *im_global, 0, opt_nx, opt_nsample, opt_trace);
 429   }
 430 #endif
 431
 432 #ifdef HAVE_MPI
 433   time_end = MPI::Wtime();
 434 #else
 435   time_end = td_current_sec();
 436 #endif
 437
 438 #ifdef HAVE_MPI
 439   if (mpi_ct.my_rank == 0)
 440 #endif
 441   {
 442     im_global->adf.arrayDataWrite ();
 443     calctime = time_end - time_start;
 444     im_global->adf.labelAdd (Array2dFileLabel::L_HISTORY, opt_desc, calctime);
 445     im_global->adf.fileClose ();
 446     if (opt_verbose)
 447       fprintf (stdout, "Time to compile image = %.2f sec\n\n", calctime);
 448   }
 449
 450   if (opt_trace >= TRACE_PHM) {
 451     double dmin, dmax;
 452     int nscale;
 453
 454     printf ("Enter display size scale (nominal = 1): ");
 455     scanf ("%d", &nscale);
 456     printf ("Enter minimum and maximum densities (min, max): ");
 457     scanf ("%lf %lf", &dmin, &dmax);
 458     image_display_scale (*im_global, nscale, dmin, dmax);
 459   }
 460
 461   phm_free (phm);
 462
 463   return (0);
 464 }
 465
 466
 467
 468 #ifdef HAVE_MPI
 469 void mpi_gather_image (ImageFile& im_global, ImageFile& im_local, const int opt_debug)
 470 {
 471   ImageFileArray vLocal = im_local.getArray();
 472   ImageFileArray vGlobal = im_global.getArray();
 473   int nyLocal = im_local.ny();
 474   int nyGlobal = im_global.ny();
 475
 476   int end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
 477   for (int iw = 0; iw <= end_work_unit; iw++) {
 478     mpi_ct.comm.Send(vLocal[iw], nyLocal, im_local.getMPIDataType(), 0, 0);
 479   }
 480
 481   if (mpi_ct.my_rank == 0) {
 482     for (int iproc = 0; iproc < mpi_ct.nproc; iproc++) {
 483       end_work_unit = mpi_ct.start_work_unit[iproc] + mpi_ct.local_work_units[iproc] - 1;
 484
 485       if (opt_debug) {
 486         fprintf(stdout, "Recv rs data from process %d (%d-%d)\n", iproc, mpi_ct.start_work_unit[iproc], end_work_unit);
 487         fflush(stdout);
 488       }
 489
 490       for (int iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
 491         MPI::Status status;
 492         mpi_ct.comm.Recv(vGlobal[iw], nyGlobal, MPI::FLOAT, iproc, 0, status);
 493       }
 494     }
 495   }
 496
 497 }
 498 #endif
 499
 500 #ifndef NO_MAIN
 501 int
 502 main (const int argc, char *const argv[])
 503 {
 504   return (phm2sdf_main(argc, argv));
 505 }
 506 #endif