src/phm2sdf.c

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (C) 1983-2000 Kevin Rosenberg
   4 **
   5 **  $Id: phm2sdf.c,v 1.2 2000/04/28 13:50:45 kevin Exp $
   6 **  $Log: phm2sdf.c,v $
   7 **  Revision 1.2  2000/04/28 13:50:45  kevin
   8 **  Removed Makefile Makefile.in that are automatically generated by autoconf
   9 **
  10 **  Revision 1.1.1.1  2000/04/28 13:02:44  kevin
  11 **  Initial CVS import for first public release
  12 **
  13 **
  14 **
  15 **  This program is free software; you can redistribute it and/or modify
  16 **  it under the terms of the GNU General Public License (version 2) as
  17 **  published by the Free Software Foundation.
  18 **
  19 **  This program is distributed in the hope that it will be useful,
  20 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 **  GNU General Public License for more details.
  23 **
  24 **  You should have received a copy of the GNU General Public License
  25 **  along with this program; if not, write to the Free Software
  26 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  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 #define O_PHANTOM       1
  39 #define O_DESC          2
  40 #define O_NSAMPLE       3
  41 #define O_FILTER        4
  42 #define O_TRACE         5
  43 #define O_VERBOSE       6
  44 #define O_HELP          7
  45 #define O_PICFILE       8
  46 #define O_FILTER_DOMAIN 9
  47 #define O_FILTER_BW     10
  48 #define O_FILTER_PARAM  11
  49 #define O_DEBUG         12
  50 #define O_VERSION       13
  51
  52 static struct option my_options[] =
  53 {
  54   {"phantom", 1, 0, O_PHANTOM},
  55   {"picfile", 1, 0, O_PICFILE},
  56   {"desc", 1, 0, O_DESC},
  57   {"nsample", 1, 0, O_NSAMPLE},
  58   {"filter", 1, 0, O_FILTER},
  59   {"filter-domain", 1, 0, O_FILTER_DOMAIN},
  60   {"filter-bw", 1, 0, O_FILTER_BW},
  61   {"filter-param", 1, 0, O_FILTER_PARAM},
  62   {"trace", 1, 0, O_TRACE},
  63   {"verbose", 0, 0, O_VERBOSE},
  64   {"debug", 0, 0, O_DEBUG},
  65   {"help", 0, 0, O_HELP},
  66   {"version", 0, 0, O_VERSION},
  67   {0, 0, 0, 0}
  68 };
  69
  70 void
  71 usage (const char *program)
  72 {
  73   fprintf(stdout,"usage: %s outfile nx ny [--phantom phantom-name] [--picfile filename] [--filter filter-name] [OPTIONS]\n",kbasename(program));
  74   fprintf(stdout,"Generate phantom image from a predefined --phantom or a --picfile\n");
  75   fprintf(stdout,"\n");
  76   fprintf(stdout,"     outfile         Name of output file for image\n");
  77   fprintf(stdout,"     nx              Number of pixels X-axis\n");
  78   fprintf(stdout,"     ny              Number of pixels Y-axis\n");
  79   fprintf(stdout,"     --phantom       Phantom to use for projection\n");
  80   fprintf(stdout,"        herman       Herman head phantom\n");
  81   fprintf(stdout,"        rowland      Rowland head phantom\n");
  82   fprintf(stdout,"        browland     Bordered Rowland head phantom\n");
  83   fprintf(stdout,"        unitpulse    Unit pulse phantom\n");
  84   fprintf(stdout,"     --picfile       Generate Phantom from picture file\n");
  85   fprintf(stdout,"     --filter        Generate Phantom from a filter function\n");
  86   fprintf(stdout,"       abs_bandlimit Abs * Bandlimiting\n");
  87   fprintf(stdout,"       abs_sinc      Abs * Sinc\n");
  88   fprintf(stdout,"       abs_cos       Abs * Cosine\n");
  89   fprintf(stdout,"       abs_hamming   Abs * Hamming\n");
  90   fprintf(stdout,"       shepp         Shepp-Logan\n");
  91   fprintf(stdout,"       bandlimit     Bandlimiting\n");
  92   fprintf(stdout,"       sinc          Sinc\n");
  93   fprintf(stdout,"       cos           Cosine\n");
  94   fprintf(stdout,"       triangle      Triangle\n");
  95   fprintf(stdout,"       hamming       Hamming\n");
  96   fprintf(stdout,"     --filter-param  Alpha level for Hamming filter\n");
  97   fprintf(stdout,"     --filter-domain Set domain of filter\n");
  98   fprintf(stdout,"         spatial     Spatial domain (default)\n");
  99   fprintf(stdout,"         freq        Frequency domain\n");
 100   fprintf(stdout,"     --filter-bw     Filter bandwidth (default = 1)\n");
 101   fprintf(stdout,"     --desc          Description of raysum\n");
 102   fprintf(stdout,"     --nsample       Number of samples per axis per pixel (default = 1)\n");
 103   fprintf(stdout,"     --trace         Trace level to use\n");
 104   fprintf(stdout,"        none         No tracing (default)\n");
 105   fprintf(stdout,"        text         Trace text level\n");
 106   fprintf(stdout,"        pic          Trace picture\n");
 107   fprintf(stdout,"        rays         Trace rays\n");
 108   fprintf(stdout,"        plot         Trace plot\n");
 109   fprintf(stdout,"        clipping     Trace clipping\n");
 110   fprintf(stdout,"     --debug         Debug mode\n");
 111   fprintf(stdout,"     --verbose       Verbose mode\n");
 112   fprintf(stdout,"     --version       Print version\n");
 113   fprintf(stdout,"     --help          Print this help message\n");
 114   exit(1);
 115 }
 116
 117 #ifdef MPI_CT
 118 void mpi_gather_image (IMAGE *im_global, IMAGE *im_local, const int opt_debug);
 119 #endif
 120
 121 int
 122 main (const int argc, char *const argv[])
 123 {
 124   IMAGE *im_global = NULL;
 125   PICTURE *pic = NULL;
 126   int opt_nx = 0, opt_ny = 0;
 127   int opt_nsample = 1;
 128   int opt_picnum = -1;
 129   int opt_filter = -1;
 130   int opt_filter_domain = D_SPATIAL;
 131   char *opt_outfile;
 132   int opt_debug = 0;
 133   char str[256];
 134   char opt_desc[256], opt_picfilename[256];
 135   char *endstr, *endptr;
 136   double opt_filter_param;
 137   double opt_filter_bw = 1.;
 138   int opt_trace = TRACE_NONE;
 139   int opt_verbose = 0;
 140   double time_start, time_end;
 141 #ifdef MPI_CT
 142   IMAGE *im_local = NULL;
 143   int mpi_argc = argc;
 144   char **mpi_argv = (char **) argv;
 145   double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
 146
 147   MPI_Init(&mpi_argc, &mpi_argv);
 148   MPI_Comm_dup (MPI_COMM_WORLD, &mpi_ct.comm);
 149   MPI_Comm_size(mpi_ct.comm, &mpi_ct.nproc);
 150   MPI_Comm_rank(mpi_ct.comm, &mpi_ct.my_rank);
 151
 152   if (mpi_ct.nproc > MPI_MAX_PROCESS) {
 153     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
 154               mpi_ct.nproc, MPI_MAX_PROCESS);
 155     exit(1);
 156   }
 157 #endif
 158
 159 #ifdef MPI_CT
 160   time_start = MPI_Wtime();
 161 #else
 162   time_start = td_current_sec();
 163 #endif
 164
 165 #ifdef MPI_CT
 166   if (mpi_ct.my_rank == 0) {
 167 #endif
 168
 169     strcpy(opt_desc, "");
 170     strcpy(opt_picfilename, "");
 171     while (1) {
 172       int c = getopt_long(argc, argv, "", my_options, NULL);
 173       char *endptr = NULL;
 174       char *endstr;
 175
 176       if (c == -1)
 177         break;
 178
 179       switch (c) {
 180       case O_PHANTOM:
 181         opt_picnum = opt_set_picture(optarg, argv[0]);
 182         break;
 183       case O_PICFILE:
 184         strncpy(opt_picfilename, optarg, sizeof(opt_picfilename));
 185         pic = create_pic_from_file(opt_picfilename);
 186 #ifdef MPI_CT
 187         if (mpi_ct.my_rank == 0)
 188           fprintf(stderr, "Can't use picture from file in MPI mode\n");
 189         exit(1);
 190 #endif
 191         break;
 192       case O_VERBOSE:
 193         opt_verbose = 1;
 194         break;
 195       case O_DEBUG:
 196         opt_debug = 1;
 197         break;
 198       case O_TRACE:
 199         opt_trace = opt_set_trace(optarg, argv[0]);
 200         break;
 201       case O_FILTER:
 202         opt_filter = opt_set_filter(optarg, argv[0]);
 203         break;
 204       case O_FILTER_DOMAIN:
 205         opt_filter_domain = opt_set_filter_domain(optarg, argv[0]);
 206         break;
 207       case O_DESC:
 208         strncpy(opt_desc, optarg, sizeof(opt_desc));
 209         break;
 210       case O_FILTER_BW:
 211         opt_filter_bw = strtod(optarg, &endptr);
 212         endstr = optarg + strlen(optarg);
 213         if (endptr != endstr) {
 214           fprintf(stderr,"Error setting --filter-bw to %s\n", optarg);
 215           usage(argv[0]);
 216           exit(1);
 217         }
 218         break;
 219       case O_FILTER_PARAM:
 220         opt_filter_param = strtod(optarg, &endptr);
 221         endstr = optarg + strlen(optarg);
 222         if (endptr != endstr) {
 223           fprintf(stderr,"Error setting --filter-param to %s\n", optarg);
 224           usage(argv[0]);
 225           exit(1);
 226         }
 227         break;
 228       case O_NSAMPLE:
 229         opt_nsample = strtol(optarg, &endptr, 10);
 230         endstr = optarg + strlen(optarg);
 231         if (endptr != endstr) {
 232           fprintf(stderr,"Error setting --nsample to %s\n", optarg);
 233           usage(argv[0]);
 234           exit(1);
 235         }
 236         break;
 237         case O_VERSION:
 238 #ifdef VERSION
 239           fprintf(stdout, "Version %s\n", VERSION);
 240 #else
 241           fprintf(stderr, "Unknown version number");
 242 #endif
 243       case O_HELP:
 244       case '?':
 245         usage(argv[0]);
 246         exit(0);
 247       default:
 248         usage(argv[0]);
 249         exit(1);
 250       }
 251     }
 252
 253     if (pic == NULL && opt_picnum == -1 && opt_filter == -1) {
 254       fprintf(stderr, "No phantom defined\n");
 255       usage(argv[0]);
 256       exit(1);
 257     }
 258
 259     if (opt_trace >= TRACE_PIC)
 260       show_pic(pic);
 261
 262     if (optind + 3 != argc) {
 263       usage(argv[0]);
 264       exit(1);
 265     }
 266     opt_outfile = argv[optind];
 267     opt_nx = strtol(argv[optind+1], &endptr, 10);
 268     endstr = argv[optind+1] + strlen(argv[optind+1]);
 269     if (endptr != endstr) {
 270       fprintf(stderr,"Error setting nx to %s\n", argv[optind+1]);
 271       usage(argv[0]);
 272       exit(1);
 273     }
 274     opt_ny = strtol(argv[optind+2], &endptr, 10);
 275     endstr = argv[optind+2] + strlen(argv[optind+2]);
 276     if (endptr != endstr) {
 277       fprintf(stderr,"Error setting ny to %s\n", argv[optind+2]);
 278       usage(argv[0]);
 279       exit(1);
 280     }
 281
 282     snprintf(str, sizeof(str), "nx=%d, ny=%d, nsample=%d, ", opt_nx, opt_ny, opt_nsample);
 283     if (opt_picfilename[0]) {
 284       strncat(str, "phantom=", sizeof(str));
 285       strncat(str, opt_picfilename, sizeof(str));
 286     }
 287     else if (opt_picnum != -1) {
 288       strncat(str, "phantom=", sizeof(str));
 289       strncat(str, name_of_picture(opt_picnum), sizeof(str));
 290     }
 291     else if (opt_filter != -1) {
 292       strncat(str, "filter=", sizeof(str));
 293       strncat(str, name_of_filter(opt_filter), sizeof(str));
 294       strncat(str, " - ", sizeof(str));
 295       strncat(str, name_of_filter_domain(opt_filter_domain), sizeof(str));
 296     }
 297     if (opt_desc[0]) {
 298       strncat(str, ": ", sizeof(str));
 299       strncat(str, opt_desc, sizeof(str));
 300     }
 301     strncpy(opt_desc, str, sizeof(opt_desc));
 302
 303     if (opt_verbose)
 304       printf("Compile Phantom to Image\n\n");
 305 #ifdef MPI_CT
 306   }
 307 #endif
 308
 309 #ifdef MPI_CT
 310   mpi_t1 = MPI_Wtime();
 311   MPI_Bcast(&opt_verbose, 1, MPI_INT, 0, mpi_ct.comm);
 312   MPI_Bcast(&opt_debug, 1, MPI_INT, 0, mpi_ct.comm);
 313   MPI_Bcast(&opt_trace, 1, MPI_INT, 0, mpi_ct.comm);
 314   MPI_Bcast(&opt_nx, 1, MPI_INT, 0, mpi_ct.comm);
 315   MPI_Bcast(&opt_ny, 1, MPI_INT, 0, mpi_ct.comm);
 316   MPI_Bcast(&opt_nsample, 1, MPI_INT, 0, mpi_ct.comm);
 317   MPI_Bcast(&opt_picnum, 1, MPI_INT, 0, mpi_ct.comm);
 318   MPI_Bcast(&opt_filter, 1, MPI_INT, 0, mpi_ct.comm);
 319   MPI_Bcast(&opt_filter_domain, 1, MPI_INT, 0, mpi_ct.comm);
 320   MPI_Bcast(&opt_filter_param, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 321   MPI_Bcast(&opt_filter_bw, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 322
 323   if (opt_verbose) {
 324     mpi_t2 = MPI_Wtime();
 325     mpi_t = mpi_t2 - mpi_t1;
 326     MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 327     if (mpi_ct.my_rank == 0)
 328       printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 329   }
 330
 331   mpi_ct_calc_work_units(opt_nx);
 332
 333   if (mpi_ct.my_rank == 0)
 334     im_global = image_create (opt_outfile, opt_nx, opt_ny);
 335
 336   im_local = image_create(NULL, opt_nx, opt_ny);
 337 #else
 338   im_global = image_create (opt_outfile, opt_nx, opt_ny);
 339 #endif
 340
 341   if (opt_picnum > 0)
 342     pic = create_pic(opt_picnum);
 343 #ifdef MPI_CT
 344   else {
 345     if (mpi_ct.my_rank == 0)
 346       fprintf(stderr, "picnum < 0\n");
 347     exit(1);
 348   }
 349 #endif
 350
 351 #ifdef MPI_CT
 352   if (pic->type == P_UNIT_PULSE) {
 353     if (mpi_ct.my_rank == 0) {
 354       im_global->v[opt_nx/2][opt_ny/2] = 1.;
 355       im_global->calctime = 0;
 356     }
 357   } else if (opt_filter != -1) {
 358     if (mpi_ct.my_rank == 0) {
 359       image_filter_init (im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 360       im_global->calctime = 0;
 361     }
 362   } else {
 363     if (opt_verbose)
 364       mpi_t1 = MPI_Wtime();
 365     pic_to_image (pic, im_local, mpi_ct.start_work_unit[mpi_ct.my_rank],
 366                   mpi_ct.local_work_units[mpi_ct.my_rank], opt_nsample, opt_trace);
 367     if (opt_verbose) {
 368       mpi_t2 = MPI_Wtime();
 369       mpi_t = mpi_t2 - mpi_t1;
 370       MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 371       if (mpi_ct.my_rank == 0)
 372         printf("Time to compile pic = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 373     }
 374     mpi_gather_image(im_global, im_local, opt_debug);
 375   }
 376 #else
 377   if (pic->type == P_UNIT_PULSE) {
 378     im_global->v[opt_nx/2][opt_ny/2] = 1.;
 379     im_global->calctime = 0;
 380   } else if (opt_filter != -1) {
 381     image_filter_init (im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
 382     im_global->calctime = 0;
 383   } else {
 384     pic_to_image (pic, im_global, 0, opt_nx, opt_nsample, opt_trace);
 385   }
 386 #endif
 387
 388 #ifdef MPI_CT
 389   time_end = MPI_Wtime();
 390   if (mpi_ct.my_rank == 0) {
 391     im_global->calctime = time_end - time_start;
 392     strncpy (im_global->remark, opt_desc, sizeof(im_global->remark));
 393   }
 394
 395   if (mpi_ct.my_rank == 0) {
 396     image_save (im_global);
 397     if (opt_verbose)
 398       fprintf (stdout, "Time to compile image = %.2f sec\n\n", im_global->calctime);
 399   }
 400 #else
 401   time_end = td_current_sec();
 402   im_global->calctime = time_end - time_start;
 403   strncpy (im_global->remark, opt_desc, sizeof(im_global->remark));
 404   image_save (im_global);
 405
 406   if (opt_trace >= TRACE_PIC)
 407     {
 408       crt_set_cpos (1, 1);
 409     }
 410
 411   if (opt_verbose)
 412     fprintf (stdout, "Time to compile image = %.2f sec\n\n", im_global->calctime);
 413 #endif
 414
 415   if (opt_trace >= TRACE_PIC)
 416     {
 417       double dmin, dmax;
 418       int nx, ny;
 419
 420       printf ("Enter size of cell (nx, ny): ");
 421       scanf ("%d %d", &nx, &ny);
 422       printf ("Enter minimum and maximum densities (min, max): ");
 423       scanf ("%lf %lf", &dmin, &dmax);
 424       //      image_paint (CRTDEV, im_global, 0, 0, nx, ny, dmin, dmax, 0);
 425       WAITKEY();
 426
 427       crt_set_palette (0, crt_calc_color (0., 0., 2.));
 428       crt_put_stra ("Finished", 8 + C_WHITE);
 429       crt_set_cpos (1, 2);
 430     }
 431
 432   return (0);
 433 }
 434
 435
 436
 437 #ifdef MPI_CT
 438 void mpi_gather_image (IMAGE *im_global, IMAGE *im_local, const int opt_debug)
 439 {
 440   int iproc;
 441   int end_work_unit;
 442   int iw;
 443
 444   end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
 445   for (iw = 0; iw <= end_work_unit; iw++) {
 446     MPI_Send(im_local->v[iw], im_local->ny, MPI_FLOAT, 0, 0, mpi_ct.comm);
 447   }
 448
 449   if (mpi_ct.my_rank == 0) {
 450     for (iproc = 0; iproc < mpi_ct.nproc; iproc++) {
 451       end_work_unit = mpi_ct.start_work_unit[iproc]
 452         + mpi_ct.local_work_units[iproc]
 453         - 1;
 454
 455       if (opt_debug) {
 456         fprintf(stdout, "Recv rs data from process %d (%d-%d)\n", iproc,
 457                 mpi_ct.start_work_unit[iproc], end_work_unit);
 458         fflush(stdout);
 459       }
 460
 461
 462       for (iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
 463         MPI_Status status;
 464
 465         MPI_Recv(im_global->v[iw], im_global->ny, MPI_FLOAT, iproc, 0, mpi_ct.comm, &status);
 466       }
 467     }
 468   }
 469
 470 }
 471 #endif