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