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