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