src/phm2rs.c

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (C) 1983-2000 Kevin Rosenberg
   4 **
   5 **  $Id: phm2rs.c,v 1.5 2000/04/30 10:13:27 kevin Exp $
   6 **  $Log: phm2rs.c,v $
   7 **  Revision 1.5  2000/04/30 10:13:27  kevin
   8 **  Fixed MPI bugs
   9 **
  10 **  Revision 1.4  2000/04/30 04:06:13  kevin
  11 **  Update Raysum i/o routines
  12 **  Fix MPI bug in ctrec (scatter_raysum) that referenced rs_global
  13 **
  14 **  Revision 1.3  2000/04/29 23:24:56  kevin
  15 **  *** empty log message ***
  16 **
  17 **  Revision 1.2  2000/04/28 13:50:45  kevin
  18 **  Removed Makefile Makefile.in that are automatically generated by autoconf
  19 **
  20 **  Revision 1.1.1.1  2000/04/28 13:02:44  kevin
  21 **  Initial CVS import for first public release
  22 **
  23 **
  24 **
  25 **  This program is free software; you can redistribute it and/or modify
  26 **  it under the terms of the GNU General Public License (version 2) as
  27 **  published by the Free Software Foundation.
  28 **
  29 **  This program is distributed in the hope that it will be useful,
  30 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  31 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  32 **  GNU General Public License for more details.
  33 **
  34 **  You should have received a copy of the GNU General Public License
  35 **  along with this program; if not, write to the Free Software
  36 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  37 ******************************************************************************/
  38 /* FILE
  39  *   phm2rs.c           Generate raysum projections from phantom object
  40  *
  41  * HISTORY
  42  *   1984 -- Final DOS version
  43  *   1999 -- First UNIX version
  44  */
  45
  46 #include "ct.h"
  47
  48 #define O_PHANTOM    1
  49 #define O_DESC       2
  50 #define O_NRAY       3
  51 #define O_ROTANGLE   4
  52 #define O_TRACE      5
  53 #define O_VERBOSE    6
  54 #define O_HELP       7
  55 #define O_PICFILE    8
  56 #define O_DEBUG     10
  57 #define O_VERSION   11
  58
  59 static struct option my_options[] =
  60 {
  61   {"phantom", 1, 0, O_PHANTOM},
  62   {"picfile", 1, 0, O_PICFILE},
  63   {"desc", 1, 0, O_DESC},
  64   {"nray", 1, 0, O_NRAY},
  65   {"rotangle", 1, 0, O_ROTANGLE},
  66   {"trace", 1, 0, O_TRACE},
  67   {"verbose", 0, 0, O_VERBOSE},
  68   {"help", 0, 0, O_HELP},
  69   {"debug", 0, 0, O_DEBUG},
  70   {"version", 0, 0, O_VERSION},
  71   {0, 0, 0, 0}
  72 };
  73
  74 void
  75 usage (const char *program)
  76 {
  77 #ifdef MPI_CT
  78 if (mpi_ct.my_rank == 0)
  79   {
  80 #endif
  81   fprintf(stdout,"usage: %s outfile ndet nview [--phantom phantom-name] [--picfile filename] [OPTIONS]\n", kbasename(program));
  82   fprintf(stdout,"Calculate raysums (projections) through phantom object, either\n");
  83   fprintf(stdout,"a predefined --phantom or a --picfile\n");
  84   fprintf(stdout,"\n");
  85   fprintf(stdout,"     outfile      Name of output file for raysums\n");
  86   fprintf(stdout,"     ndet         Number of detectors\n");
  87   fprintf(stdout,"     nview        Number of rotated views\n");
  88   fprintf(stdout,"     --phantom    Phantom to use for projection\n");
  89   fprintf(stdout,"        herman    Herman head phantom\n");
  90   fprintf(stdout,"        rowland   Rowland head phantom\n");
  91   fprintf(stdout,"        browland  Bordered Rowland head phantom\n");
  92   fprintf(stdout,"        unitpulse Unit pulse phantom\n");
  93   fprintf(stdout,"     --picfile    Get Phantom from picture file\n");
  94   fprintf(stdout,"     --desc       Description of raysum\n");
  95   fprintf(stdout,"     --nray       Number of rays per detector (default = 1)\n");
  96   fprintf(stdout,"     --rotangle   Degrees to rotate view through, multiple of PI (default = 1)\n");
  97   fprintf(stdout,"     --trace      Trace level to use\n");
  98   fprintf(stdout,"        none      No tracing (default)\n");
  99   fprintf(stdout,"        text      Trace text level\n");
 100   fprintf(stdout,"        pic       Trace picture\n");
 101   fprintf(stdout,"        rays      Trace rays\n");
 102   fprintf(stdout,"        plot      Trace plot\n");
 103   fprintf(stdout,"        clipping  Trace clipping\n");
 104   fprintf(stdout,"     --verbose    Verbose mode\n");
 105   fprintf(stdout,"     --debug      Debug mode\n");
 106   fprintf(stdout,"     --version    Print version\n");
 107   fprintf(stdout,"     --help       Print this help message\n");
 108 #ifdef MPI_CT
 109   }
 110   MPI_Abort(mpi_ct.comm, 1);
 111 #endif
 112   exit(1);
 113 }
 114
 115
 116 void mpi_gather_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug);
 117
 118 int
 119 main (const int argc, char *const argv[])
 120 {
 121   DETECTOR *det;
 122   PICTURE *pic = NULL;
 123   RAYSUM *rs_global;
 124   char str[256], *opt_outfile, opt_desc[MAXREMARK+1];
 125   char opt_picfilename[256];
 126   char *endptr, *endstr;
 127   int opt_ndet, opt_nview;
 128   int opt_nray = 1;
 129   int opt_trace = 0, opt_picnum = -1;
 130   int opt_verbose = 0;
 131   int opt_debug = 0;
 132   double opt_rotangle = 1;
 133   double time_start, time_end;
 134 #ifdef MPI_CT
 135   RAYSUM *rs_local;
 136   int mpi_argc = argc;
 137   char **mpi_argv = (char **) argv;
 138   double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
 139
 140   MPI_Init(&mpi_argc, &mpi_argv);
 141   MPI_Comm_dup (MPI_COMM_WORLD, &mpi_ct.comm);
 142   MPI_Comm_size(mpi_ct.comm, &mpi_ct.nproc);
 143   MPI_Comm_rank(mpi_ct.comm, &mpi_ct.my_rank);
 144
 145   if (mpi_ct.nproc > MPI_MAX_PROCESS) {
 146     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
 147               mpi_ct.nproc, MPI_MAX_PROCESS);
 148     exit(1);
 149   }
 150   time_start = MPI_Wtime();
 151 #else
 152   time_start = td_current_sec();
 153 #endif
 154
 155 #ifdef MPI_CT
 156   if (mpi_ct.my_rank == 0) {
 157 #endif
 158     strcpy(opt_desc, "");
 159     strcpy(opt_picfilename, "");
 160     while (1) {
 161       int c = getopt_long(argc, argv, "", my_options, NULL);
 162       char *endptr = NULL;
 163       char *endstr;
 164
 165       if (c == -1)
 166         break;
 167
 168       switch (c) {
 169       case O_PHANTOM:
 170         opt_picnum = opt_set_picture(optarg, argv[0]);
 171         pic = create_pic(opt_picnum);
 172         break;
 173       case O_PICFILE:
 174 #ifdef MPI_CT
 175         if (mpi_ct.my_rank == 0)
 176           fprintf(stderr, "Can not read picture from file in MPI mode\n");
 177         exit(1);
 178 #endif
 179         strncpy(opt_picfilename, optarg, sizeof(opt_picfilename));
 180         pic = create_pic_from_file(opt_picfilename);
 181         break;
 182       case O_VERBOSE:
 183         opt_verbose = 1;
 184         break;
 185       case O_DEBUG:
 186         opt_debug = 1;
 187         break;
 188         break;
 189       case O_TRACE:
 190         opt_trace = opt_set_trace(optarg, argv[0]);
 191         break;
 192       case O_DESC:
 193         strncpy(opt_desc, optarg, sizeof(opt_desc));
 194         break;
 195       case O_ROTANGLE:
 196         opt_rotangle = strtod(optarg, &endptr);
 197         endstr = optarg + strlen(optarg);
 198         if (endptr != endstr) {
 199           fprintf(stderr,"Error setting --rotangle to %s\n", optarg);
 200           usage(argv[0]);
 201           exit(1);
 202         }
 203         break;
 204       case O_NRAY:
 205         opt_nray = strtol(optarg, &endptr, 10);
 206         endstr = optarg + strlen(optarg);
 207         if (endptr != endstr) {
 208           fprintf(stderr,"Error setting --nray to %s\n", optarg);
 209           usage(argv[0]);
 210           exit(1);
 211         }
 212         break;
 213         case O_VERSION:
 214 #ifdef VERSION
 215           fprintf(stdout, "Version %s\n", VERSION);
 216 #else
 217           fprintf(stderr, "Unknown version number");
 218 #endif
 219           exit(0);
 220       case O_HELP:
 221       case '?':
 222         usage(argv[0]);
 223         exit(0);
 224       default:
 225         usage(argv[0]);
 226         exit(1);
 227       }
 228     }
 229
 230     if (pic == NULL) {
 231       fprintf(stderr, "No phantom defined\n");
 232       usage(argv[0]);
 233       exit(1);
 234     }
 235     if (optind + 3 != argc) {
 236       usage(argv[0]);
 237       exit(1);
 238     }
 239     opt_outfile = argv[optind];
 240     opt_ndet = strtol(argv[optind+1], &endptr, 10);
 241     endstr = argv[optind+1] + strlen(argv[optind+1]);
 242     if (endptr != endstr) {
 243       fprintf(stderr,"Error setting --ndet to %s\n", argv[optind+1]);
 244       usage(argv[0]);
 245       exit(1);
 246     }
 247     opt_nview = strtol(argv[optind+2], &endptr, 10);
 248     endstr = argv[optind+2] + strlen(argv[optind+2]);
 249     if (endptr != endstr) {
 250       fprintf(stderr,"Error setting --nview to %s\n", argv[optind+2]);
 251       usage(argv[0]);
 252       exit(1);
 253     }
 254
 255     snprintf(str, sizeof(str),
 256              "Raysum_Collect: NDet=%d, Nview=%d, NRay=%d, RotAngle=%.2f, ",
 257              opt_ndet, opt_nview, opt_nray, opt_rotangle);
 258     if (opt_picfilename[0]) {
 259       strncat(str, "Phantom=", sizeof(str));
 260       strncat(str, opt_picfilename, sizeof(str));
 261     } else if (opt_picnum != -1) {
 262       strncat(str, "Phantom=", sizeof(str));
 263       strncat(str, name_of_picture(opt_picnum), sizeof(str));
 264     }
 265     if (opt_desc[0]) {
 266       strncat(str, ": ", sizeof(str));
 267       strncat(str, opt_desc, sizeof(str));
 268     }
 269     strncpy(opt_desc, str, sizeof(opt_desc));
 270 #ifdef MPI_CT
 271   }
 272 #endif
 273
 274 #ifdef MPI_CT
 275   MPI_Barrier(mpi_ct.comm);
 276   MPI_Bcast(&opt_rotangle, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 277   MPI_Bcast(&opt_nview, 1, MPI_INT, 0, mpi_ct.comm);
 278   MPI_Bcast(&opt_ndet, 1, MPI_INT, 0, mpi_ct.comm);
 279   MPI_Bcast(&opt_nray, 1, MPI_INT, 0, mpi_ct.comm);
 280   MPI_Bcast(&opt_picnum, 1, MPI_INT, 0, mpi_ct.comm);
 281   MPI_Bcast(&opt_verbose, 1, MPI_INT, 0, mpi_ct.comm);
 282   MPI_Bcast(&opt_debug, 1, MPI_INT, 0, mpi_ct.comm);
 283   MPI_Bcast(&opt_trace, 1, MPI_INT, 0, mpi_ct.comm);
 284
 285   if (mpi_ct.my_rank > 0 && opt_picnum > 0)
 286     pic = create_pic(opt_picnum);
 287 #endif
 288
 289   opt_rotangle *= PI;
 290   det = detect_create (pic, opt_ndet, opt_nview, opt_nray, opt_rotangle);
 291
 292 #ifdef MPI_CT
 293   mpi_ct_calc_work_units(opt_nview);
 294   if (mpi_ct.my_rank == 0) {
 295     rs_global = raysum_create_from_det (opt_outfile, det);
 296     raysum_alloc_views(rs_global);
 297   }
 298
 299   rs_local = raysum_create_from_det (NULL, det);
 300   rs_local->nview = mpi_ct.local_work_units[mpi_ct.my_rank];
 301   if (opt_debug)
 302     printf("rs_local->nview = %d (process %d)\n", rs_local->nview, mpi_ct.my_rank);
 303
 304   if (opt_verbose)
 305     mpi_t1 = MPI_Wtime();
 306   raysum_collect (rs_local, det, pic, mpi_ct.start_work_unit[mpi_ct.my_rank], opt_trace, FALSE);
 307   if (opt_verbose) {
 308     mpi_t2 = MPI_Wtime();
 309     mpi_t = mpi_t2 - mpi_t1;
 310     MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 311     if (mpi_ct.my_rank == 0)
 312       printf("Time to collect rs = %f secs, Max = %f secs\n", mpi_t, mpi_t_g);
 313   }
 314
 315   if (opt_verbose)
 316     mpi_t1 = MPI_Wtime();
 317   mpi_gather_rs(rs_global, rs_local, opt_debug);
 318   if (opt_verbose) {
 319     mpi_t2 = MPI_Wtime();
 320     mpi_t = mpi_t2 - mpi_t1;
 321     MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 322     if (mpi_ct.my_rank == 0)
 323       printf("Time to gather rs = %f secs, Max = %f secs\n", mpi_t, mpi_t_g);
 324   }
 325 #else
 326   rs_global = raysum_create_from_det (opt_outfile, det);
 327   raysum_collect (rs_global, det, pic, 0, opt_trace, FALSE);
 328 #endif
 329
 330 #ifdef MPI_CT
 331   time_end = MPI_Wtime();
 332   if (mpi_ct.my_rank == 0) {
 333     rs_global->calctime = time_end - time_start;
 334     strncpy (rs_global->remark, opt_desc, sizeof(rs_global->remark));
 335     if (opt_verbose) {
 336       fprintf(stdout, "Remark: %s\n", rs_global->remark);
 337       fprintf(stdout, "Time active = %.2f secs\n", rs_global->calctime);
 338     }
 339   }
 340 #else
 341   time_end = td_current_sec();
 342   rs_global->calctime = time_end - time_start;
 343   strncpy (rs_global->remark, opt_desc, sizeof(rs_global->remark));
 344   if (opt_verbose) {
 345     fprintf(stdout, "Remark: %s\n", rs_global->remark);
 346     fprintf(stdout, "Time active = %.2f secs\n", rs_global->calctime);
 347   }
 348 #endif
 349
 350 #ifdef MPI_CT
 351   if (mpi_ct.my_rank == 0) {
 352     raysum_write (rs_global);
 353     raysum_close (rs_global);
 354   }
 355 #else
 356   raysum_write (rs_global);
 357   raysum_close (rs_global);
 358 #endif
 359
 360   detect_free (det);
 361
 362   if (opt_trace >= TRACE_PIC)
 363     {
 364       crt_set_cpos (1, 1);
 365       printf("Finished\n");
 366     }
 367
 368   return (0);
 369 }
 370
 371
 372 #ifdef MPI_CT
 373 void mpi_gather_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
 374 {
 375   int iproc;
 376   int end_work_unit;
 377   int iw;
 378
 379   end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
 380   for (iw = 0; iw <= end_work_unit; iw++) {
 381     MPI_Send(&rs_local->view[iw]->view_angle, 1, MPI_DOUBLE, 0, 0, mpi_ct.comm);
 382     MPI_Send(&rs_local->view[iw]->ndet, 1, MPI_INT, 0, 0, mpi_ct.comm);
 383     MPI_Send(rs_local->view[iw]->detval, rs_local->ndet, MPI_FLOAT, 0, 0, mpi_ct.comm);
 384   }
 385
 386   if (mpi_ct.my_rank == 0) {
 387     for (iproc = 0; iproc < mpi_ct.nproc; iproc++) {
 388       end_work_unit = mpi_ct.start_work_unit[iproc]
 389         + mpi_ct.local_work_units[iproc]
 390         - 1;
 391
 392       if (opt_debug)
 393         fprintf(stdout, "Recv rs data from process %d\n", iproc);
 394
 395       for (iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
 396         MPI_Status status;
 397
 398         MPI_Recv(&rs_global->view[iw]->view_angle, 1, MPI_DOUBLE, iproc, 0, mpi_ct.comm, &status);
 399         MPI_Recv(&rs_global->view[iw]->ndet, 1, MPI_INT, iproc, 0, mpi_ct.comm, &status);
 400         MPI_Recv(rs_global->view[iw]->detval, rs_global->ndet, MPI_FLOAT, iproc, 0, mpi_ct.comm, &status);
 401       }
 402     }
 403   }
 404
 405 }
 406 #endif