src/phm2rs.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          phm2rs.cpp
   5 **   Purpose:       Take projections of a phantom object
   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: phm2rs.cpp,v 1.3 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
  29 #include "ct.h"
  30
  31 enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE,
  32        O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
  33
  34 static struct option phm2rs_options[] =
  35 {
  36   {"phantom", 1, 0, O_PHANTOM},
  37   {"phmfile", 1, 0, O_PHMFILE},
  38   {"desc", 1, 0, O_DESC},
  39   {"nray", 1, 0, O_NRAY},
  40   {"rotangle", 1, 0, O_ROTANGLE},
  41   {"trace", 1, 0, O_TRACE},
  42   {"verbose", 0, 0, O_VERBOSE},
  43   {"help", 0, 0, O_HELP},
  44   {"debug", 0, 0, O_DEBUG},
  45   {"version", 0, 0, O_VERSION},
  46   {0, 0, 0, 0}
  47 };
  48
  49
  50 void
  51 phm2rs_usage (const char *program)
  52 {
  53 #ifdef HAVE_MPI
  54 if (mpi_ct.my_rank == 0)
  55   {
  56 #endif
  57   fprintf(stdout,"usage: %s outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]\n", kbasename(program));
  58   fprintf(stdout,"Calculate raysums (projections) through phantom object, either\n");
  59   fprintf(stdout,"a predefined --phantom or a --phmfile\n");
  60   fprintf(stdout,"\n");
  61   fprintf(stdout,"     outfile      Name of output file for raysums\n");
  62   fprintf(stdout,"     ndet         Number of detectors\n");
  63   fprintf(stdout,"     nview        Number of rotated views\n");
  64   fprintf(stdout,"     --phantom    Phantom to use for projection\n");
  65   fprintf(stdout,"        herman    Herman head phantom\n");
  66   fprintf(stdout,"        rowland   Rowland head phantom\n");
  67   fprintf(stdout,"        browland  Bordered Rowland head phantom\n");
  68   fprintf(stdout,"        unitpulse Unit pulse phantom\n");
  69   fprintf(stdout,"     --phmfile    Get Phantom from phantom file\n");
  70   fprintf(stdout,"     --desc       Description of raysum\n");
  71   fprintf(stdout,"     --nray       Number of rays per detector (default = 1)\n");
  72   fprintf(stdout,"     --rotangle   Degrees to rotate view through, multiple of PI (default = 1)\n");
  73   fprintf(stdout,"     --trace      Trace level to use\n");
  74   fprintf(stdout,"        none      No tracing (default)\n");
  75   fprintf(stdout,"        text      Trace text level\n");
  76   fprintf(stdout,"        phm       Trace phantom image\n");
  77   fprintf(stdout,"        rays      Trace rays\n");
  78   fprintf(stdout,"        plot      Trace plot\n");
  79   fprintf(stdout,"        clipping  Trace clipping\n");
  80   fprintf(stdout,"     --verbose    Verbose mode\n");
  81   fprintf(stdout,"     --debug      Debug mode\n");
  82   fprintf(stdout,"     --version    Print version\n");
  83   fprintf(stdout,"     --help       Print this help message\n");
  84 #ifdef HAVE_MPI
  85   }
  86   mpi_ct.comm.Abort(1);
  87 #endif
  88 }
  89
  90 #ifdef HAVE_MPI
  91 void mpi_gather_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug);
  92 #endif
  93
  94 int
  95 phm2rs_main (const int argc, char *const argv[])
  96 {
  97   DETECTOR *det;
  98   PHANTOM *phm = NULL;
  99   RAYSUM *rs_global = NULL;
 100   char str[256], *opt_outfile = NULL, opt_desc[MAXREMARK+1];
 101   char opt_phmfilename[256];
 102   char *endptr, *endstr;
 103   int opt_ndet, opt_nview;
 104   int opt_nray = 1;
 105   int opt_trace = 0, opt_phmnum = -1;
 106   int opt_verbose = 0;
 107   int opt_debug = 0;
 108   double opt_rotangle = 1;
 109   double time_start = 0, time_end = 0;
 110
 111 #ifndef HAVE_MPI
 112   time_start = td_current_sec();
 113 #else
 114   RAYSUM *rs_local;
 115   int mpi_argc = argc;
 116   char **mpi_argv = (char **) argv;
 117   double mpi_t1 = 0, mpi_t2 = 0, mpi_t, mpi_t_g;
 118
 119   MPI::Init (mpi_argc, mpi_argv);
 120   mpi_ct.comm = MPI::COMM_WORLD.Dup ();
 121   mpi_ct.nproc = mpi_ct.comm.Get_size();
 122   mpi_ct.my_rank = mpi_ct.comm.Get_rank();
 123
 124   if (mpi_ct.nproc > CT_MPI_MAX_PROCESS) {
 125     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
 126               mpi_ct.nproc, CT_MPI_MAX_PROCESS);
 127   }
 128
 129   time_start = MPI::Wtime();
 130 #endif
 131
 132 #ifdef HAVE_MPI
 133   if (mpi_ct.my_rank == 0) {
 134 #endif
 135     strcpy(opt_desc, "");
 136     strcpy(opt_phmfilename, "");
 137     while (1) {
 138       int c = getopt_long(argc, argv, "", phm2rs_options, NULL);
 139       char *endptr = NULL;
 140       char *endstr;
 141
 142       if (c == -1)
 143         break;
 144
 145       switch (c) {
 146       case O_PHANTOM:
 147         if ((opt_phmnum = opt_set_phantom (optarg)) < 0) {
 148           phm2rs_usage(argv[0]);
 149           return (1);
 150         }
 151         phm = phm_create (opt_phmnum);
 152         break;
 153       case O_PHMFILE:
 154 #ifdef HAVE_MPI
 155         if (mpi_ct.my_rank == 0)
 156           fprintf(stderr, "Can not read phantom from file in MPI mode\n");
 157         return (1);
 158 #endif
 159         strncpy(opt_phmfilename, optarg, sizeof(opt_phmfilename));
 160         phm = phm_create_from_file(opt_phmfilename);
 161         break;
 162       case O_VERBOSE:
 163         opt_verbose = 1;
 164         break;
 165       case O_DEBUG:
 166         opt_debug = 1;
 167         break;
 168         break;
 169       case O_TRACE:
 170         if ((opt_trace = opt_set_trace(optarg)) < 0) {
 171           phm2rs_usage(argv[0]);
 172           return (1);
 173         }
 174         break;
 175       case O_DESC:
 176         strncpy(opt_desc, optarg, sizeof(opt_desc));
 177         break;
 178       case O_ROTANGLE:
 179         opt_rotangle = strtod(optarg, &endptr);
 180         endstr = optarg + strlen(optarg);
 181         if (endptr != endstr) {
 182           fprintf(stderr,"Error setting --rotangle to %s\n", optarg);
 183           phm2rs_usage(argv[0]);
 184           return (1);
 185         }
 186         break;
 187       case O_NRAY:
 188         opt_nray = strtol(optarg, &endptr, 10);
 189         endstr = optarg + strlen(optarg);
 190         if (endptr != endstr) {
 191           fprintf(stderr,"Error setting --nray to %s\n", optarg);
 192           phm2rs_usage(argv[0]);
 193           return (1);
 194         }
 195         break;
 196         case O_VERSION:
 197 #ifdef VERSION
 198           fprintf(stdout, "Version %s\n", VERSION);
 199 #else
 200           fprintf(stderr, "Unknown version number");
 201 #endif
 202           exit(0);
 203       case O_HELP:
 204       case '?':
 205         phm2rs_usage(argv[0]);
 206         return (0);
 207       default:
 208         phm2rs_usage(argv[0]);
 209         return (1);
 210       }
 211     }
 212
 213     if (phm == NULL) {
 214       fprintf(stderr, "No phantom defined\n");
 215       phm2rs_usage(argv[0]);
 216       return (1);
 217     }
 218     if (optind + 3 != argc) {
 219       phm2rs_usage(argv[0]);
 220       return (1);
 221     }
 222
 223     opt_outfile = argv[optind];
 224     opt_ndet = strtol(argv[optind+1], &endptr, 10);
 225     endstr = argv[optind+1] + strlen(argv[optind+1]);
 226     if (endptr != endstr) {
 227       fprintf(stderr,"Error setting --ndet to %s\n", argv[optind+1]);
 228       phm2rs_usage(argv[0]);
 229       return (1);
 230     }
 231     opt_nview = strtol(argv[optind+2], &endptr, 10);
 232     endstr = argv[optind+2] + strlen(argv[optind+2]);
 233     if (endptr != endstr) {
 234       fprintf(stderr,"Error setting --nview to %s\n", argv[optind+2]);
 235       phm2rs_usage(argv[0]);
 236       return (1);
 237     }
 238
 239     snprintf(str, sizeof(str),
 240              "Raysum_Collect: NDet=%d, Nview=%d, NRay=%d, RotAngle=%.2f, ",
 241              opt_ndet, opt_nview, opt_nray, opt_rotangle);
 242     if (opt_phmfilename[0]) {
 243       strncat(str, "Phantom=", sizeof(str));
 244       strncat(str, opt_phmfilename, sizeof(str));
 245     } else if (opt_phmnum != -1) {
 246       strncat(str, "Phantom=", sizeof(str));
 247       strncat(str, name_of_phantom(opt_phmnum), sizeof(str));
 248     }
 249     if (opt_desc[0]) {
 250       strncat(str, ": ", sizeof(str));
 251       strncat(str, opt_desc, sizeof(str));
 252     }
 253     strncpy(opt_desc, str, sizeof(opt_desc));
 254 #ifdef HAVE_MPI
 255   }
 256 #endif
 257
 258 #ifdef HAVE_MPI
 259   mpi_ct.comm.Barrier ();
 260   mpi_ct.comm.Bcast (&opt_rotangle, 1, MPI::DOUBLE, 0);
 261   mpi_ct.comm.Bcast (&opt_nview, 1, MPI::INT, 0);
 262   mpi_ct.comm.Bcast (&opt_ndet, 1, MPI::INT, 0);
 263   mpi_ct.comm.Bcast (&opt_nray, 1, MPI::INT, 0);
 264   mpi_ct.comm.Bcast (&opt_phmnum, 1, MPI::INT, 0);
 265   mpi_ct.comm.Bcast (&opt_verbose, 1, MPI::INT, 0);
 266   mpi_ct.comm.Bcast (&opt_debug, 1, MPI::INT, 0);
 267   mpi_ct.comm.Bcast (&opt_trace, 1, MPI::INT, 0);
 268
 269   if (mpi_ct.my_rank > 0 && opt_phmnum >= 0)
 270     phm = phm_create (opt_phmnum);
 271 #endif
 272
 273   opt_rotangle *= PI;
 274   det = detector_create (phm, DETECTOR_PARALLEL, opt_ndet, opt_nview, opt_nray, opt_rotangle);
 275
 276 #ifdef HAVE_MPI
 277   mpi_ct_calc_work_units(opt_nview);
 278   if (mpi_ct.my_rank == 0) {
 279     rs_global = raysum_create_from_det (opt_outfile, det);
 280     raysum_alloc_views(rs_global);
 281   }
 282
 283   rs_local = raysum_create_from_det (NULL, det);
 284   rs_local->nview = mpi_ct.local_work_units[mpi_ct.my_rank];
 285   if (opt_debug)
 286     printf("rs_local->nview = %d (process %d)\n", rs_local->nview, mpi_ct.my_rank);
 287
 288   if (opt_verbose)
 289     mpi_t1 = MPI::Wtime();
 290   raysum_collect (rs_local, det, phm, mpi_ct.start_work_unit[mpi_ct.my_rank], opt_trace, FALSE);
 291   if (opt_verbose) {
 292     mpi_t2 = MPI::Wtime();
 293     mpi_t = mpi_t2 - mpi_t1;
 294     mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
 295     if (mpi_ct.my_rank == 0)
 296       printf("Time to collect rs = %f secs, Max = %f secs\n", mpi_t, mpi_t_g);
 297   }
 298
 299   if (opt_verbose)
 300     mpi_t1 = MPI::Wtime();
 301   mpi_gather_rs(rs_global, rs_local, opt_debug);
 302   if (opt_verbose) {
 303     mpi_t2 = MPI::Wtime();
 304     mpi_t = mpi_t2 - mpi_t1;
 305     mpi_ct.comm.Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0);
 306     if (mpi_ct.my_rank == 0)
 307       printf("Time to gather rs = %f secs, Max = %f secs\n", mpi_t, mpi_t_g);
 308   }
 309 #else
 310   rs_global = raysum_create_from_det (opt_outfile, det);
 311   raysum_collect (rs_global, det, phm, 0, opt_trace, FALSE);
 312 #endif
 313
 314 #ifndef HAVE_MPI
 315   time_end = td_current_sec();
 316 #else
 317   time_end = MPI::Wtime();
 318   if (mpi_ct.my_rank == 0)
 319 #endif
 320     {
 321       rs_global->calctime = time_end - time_start;
 322       strncpy (rs_global->remark, opt_desc, sizeof(rs_global->remark));
 323       if (opt_verbose) {
 324         fprintf(stdout, "Remark: %s\n", rs_global->remark);
 325         fprintf(stdout, "Time active = %.2f secs\n", rs_global->calctime);
 326       }
 327     }
 328
 329 #ifdef HAVE_MPI
 330   if (mpi_ct.my_rank == 0)
 331 #endif
 332     {
 333       raysum_write (rs_global);
 334       raysum_close (rs_global);
 335     }
 336
 337   phm_free (phm);
 338   detector_free (det);
 339
 340   return (0);
 341 }
 342
 343
 344 /* FUNCTION
 345  *    mpi_gather_rs
 346  *
 347  * SYNOPSIS
 348  *    Gather's raysums from all processes in rs_local in rs_global in process 0
 349  */
 350
 351 #ifdef HAVE_MPI
 352 void mpi_gather_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
 353 {
 354   int iproc;
 355   int end_work_unit;
 356   int iw;
 357
 358   end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
 359   for (iw = 0; iw <= end_work_unit; iw++) {
 360     mpi_ct.comm.Send(&rs_local->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0);
 361     mpi_ct.comm.Send(&rs_local->view[iw]->ndet, 1, MPI::INT, 0, 0);
 362     mpi_ct.comm.Send(rs_local->view[iw]->detval, rs_local->ndet, MPI::FLOAT, 0, 0);
 363   }
 364
 365   if (mpi_ct.my_rank == 0) {
 366     for (iproc = 0; iproc < mpi_ct.nproc; iproc++) {
 367       end_work_unit = mpi_ct.start_work_unit[iproc]
 368         + mpi_ct.local_work_units[iproc]
 369         - 1;
 370
 371       if (opt_debug)
 372         fprintf(stdout, "Recv rs data from process %d\n", iproc);
 373
 374       for (iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
 375         MPI::Status status;
 376
 377         mpi_ct.comm.Recv(&rs_global->view[iw]->view_angle, 1, MPI::DOUBLE, iproc, 0, status);
 378         mpi_ct.comm.Recv(&rs_global->view[iw]->ndet, 1, MPI::INT, iproc, 0, status);
 379         mpi_ct.comm.Recv(rs_global->view[iw]->detval, rs_global->ndet, MPI::FLOAT, iproc, 0, status);
 380       }
 381     }
 382   }
 383
 384 }
 385 #endif
 386
 387
 388 #ifndef NO_MAIN
 389 int
 390 main (const int argc, char *const argv[])
 391 {
 392   return (phm2rs_main(argc, argv));
 393 }
 394 #endif
 395