src/ctrec.c

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (C) 1983-2000 Kevin Rosenberg
   4 **
   5 **  $Id: ctrec.c,v 1.11 2000/05/16 04:33:59 kevin Exp $
   6 **  $Log: ctrec.c,v $
   7 **  Revision 1.11  2000/05/16 04:33:59  kevin
   8 **  Improved option processing
   9 **
  10 **  Revision 1.10  2000/05/11 01:06:30  kevin
  11 **  Changed sprintf to snprintf
  12 **
  13 **  Revision 1.9  2000/05/08 20:02:32  kevin
  14 **  ANSI C changes
  15 **
  16 **  Revision 1.8  2000/05/04 18:16:34  kevin
  17 **  renamed filter definitions
  18 **
  19 **  Revision 1.7  2000/05/03 08:49:50  kevin
  20 **  Code cleanup
  21 **
  22 **  Revision 1.6  2000/05/02 15:31:47  kevin
  23 **  code cleaning
  24 **
  25 **  Revision 1.5  2000/04/30 11:41:06  kevin
  26 **  Cleaned up debugging code
  27 **
  28 **  Revision 1.4  2000/04/30 10:13:27  kevin
  29 **  Fixed MPI bugs
  30 **
  31 **  Revision 1.3  2000/04/30 04:06:13  kevin
  32 **  Update Raysum i/o routines
  33 **  Fix MPI bug in ctrec (scatter_raysum) that referenced rs_global
  34 **
  35 **  Revision 1.2  2000/04/29 23:24:56  kevin
  36 **  *** empty log message ***
  37 **
  38 **  Revision 1.1.1.1  2000/04/28 13:02:44  kevin
  39 **  Initial CVS import for first public release
  40 **
  41 **
  42 **
  43 **  This program is free software; you can redistribute it and/or modify
  44 **  it under the terms of the GNU General Public License (version 2) as
  45 **  published by the Free Software Foundation.
  46 **
  47 **  This program is distributed in the hope that it will be useful,
  48 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  49 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  50 **  GNU General Public License for more details.
  51 **
  52 **  You should have received a copy of the GNU General Public License
  53 **  along with this program; if not, write to the Free Software
  54 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  55 ******************************************************************************/
  56
  57 /* FILE
  58  *   ctrec.c            Reconstruct an image from raysums
  59  *
  60  * DATE
  61  *   Aug 84
  62  *   Jul 99  -- Converted to ANSI C
  63  *              Added MPI parallel processing
  64  */
  65
  66 #include "ct.h"
  67
  68 enum {O_INTERP, O_FILTER, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
  69
  70 static struct option my_options[] =
  71 {
  72   {"interp", 1, 0, O_INTERP},
  73   {"filter", 1, 0, O_FILTER},
  74   {"filter-param", 1, 0, O_FILTER_PARAM},
  75   {"backproj", 1, 0, O_BACKPROJ},
  76   {"trace", 1, 0, O_TRACE},
  77   {"debug", 0, 0, O_DEBUG},
  78   {"verbose", 0, 0, O_VERBOSE},
  79   {"help", 0, 0, O_HELP},
  80   {"version", 0, 0, O_VERSION},
  81   {0, 0, 0, 0}
  82 };
  83
  84 void
  85 ctrec_usage (const char *program)
  86 {
  87   fprintf(stdout,"usage: %s raysum-file image-file nx-image ny-image [OPTIONS]\n", kbasename(program));
  88   fprintf(stdout,"Image reconstruction from raysum projections\n");
  89   fprintf(stdout,"\n");
  90   fprintf(stdout,"   raysum-file     Input raysum file\n");
  91   fprintf(stdout,"   image-file      Output image file in SDF2D format\n");
  92   fprintf(stdout,"   nx-image        Number of columns in output image\n");
  93   fprintf(stdout,"   ny-image        Number of rows in output image\n");
  94   fprintf(stdout,"   --interp        Interpolation method during backprojection\n");
  95   fprintf(stdout,"       nearest     Nearest neighbor interpolation\n");
  96   fprintf(stdout,"       linear      Linear interpolation\n");
  97   fprintf(stdout,"       bspline     B-spline interpolation\n");
  98   fprintf(stdout,"    --filter       Filter name\n");
  99   fprintf(stdout,"       abs_bandlimit Abs * Bandlimiting (default)\n");
 100   fprintf(stdout,"       abs_sinc      Abs * Sinc\n");
 101   fprintf(stdout,"       abs_cos       Abs * Cosine\n");
 102   fprintf(stdout,"       abs_hamming   Abs * Hamming\n");
 103   fprintf(stdout,"       shepp         Shepp-Logan\n");
 104   fprintf(stdout,"       bandlimit     Bandlimiting\n");
 105   fprintf(stdout,"       sinc          Sinc\n");
 106   fprintf(stdout,"       cos           Cosine\n");
 107   fprintf(stdout,"       triangle      Triangle\n");
 108   fprintf(stdout,"       hamming       Hamming\n");
 109   fprintf(stdout,"    --backproj     Backprojection Method\n");
 110   fprintf(stdout,"       trig        Trigometric functions at every point\n");
 111   fprintf(stdout,"       table       Trigometric functions with precalculated table\n");
 112   fprintf(stdout,"       diff        Difference method\n");
 113   fprintf(stdout,"       diff2       Optimized difference method (default)\n");
 114   fprintf(stdout,"       idiff2      Optimized difference method with integer math\n");
 115   fprintf(stdout,"    --filter-param Alpha level for Hamming filter\n");
 116   fprintf(stdout,"    --trace        Set tracing to level\n");
 117   fprintf(stdout,"         none      No tracing (default)\n");
 118   fprintf(stdout,"         status    Status tracing\n");
 119   fprintf(stdout,"         phm       Trace phantom\n");
 120   fprintf(stdout,"         rays      Trace allrays\n");
 121   fprintf(stdout,"         plot      Trace plotting\n");
 122   fprintf(stdout,"         clipping  Trace clipping\n");
 123   fprintf(stdout,"    --verbose      Turn on verbose mode\n");
 124   fprintf(stdout,"    --debug        Turn on debug mode\n");
 125   fprintf(stdout,"    --version      Print version\n");
 126   fprintf(stdout,"    --help         Print this help message\n");
 127 }
 128
 129
 130 #ifdef MPI_CT
 131 static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int debug);
 132 #endif
 133
 134 static void print_raysum_info(const RAYSUM *rs);
 135
 136 int
 137 ctrec_main (const int argc, char *const argv[])
 138 {
 139   IMAGE *im_global = NULL;
 140   RAYSUM *rs_global = NULL;
 141   char *rs_name, *im_filename = NULL;
 142   char remark[MAXREMARK];
 143   char filt_name[80];
 144   int nx, ny;
 145   double time_start = 0, time_end = 0;
 146   char *endptr;
 147   int opt_verbose = 0;
 148   int opt_debug = 0;
 149   int opt_trace = TRACE_NONE;
 150   double opt_filter_param = 1;
 151   int opt_filter = FILTER_ABS_BANDLIMIT;
 152   int opt_interp = I_LINEAR;
 153   int opt_interp_param = 1;
 154   BackprojType opt_backproj = O_BPROJ_DIFF2;
 155 #ifdef MPI_CT
 156   IMAGE *im_local;
 157   RAYSUM *rs_local;
 158   int mpi_argc = argc;
 159   char **mpi_argv = (char **) argv;
 160   int mpi_nview, mpi_ndet;
 161   double mpi_detinc, mpi_rotinc, mpi_phmlen;
 162   double mpi_t1, mpi_t2, mpi_t, mpi_t_g;
 163   int ix;
 164
 165   MPI_Init(&mpi_argc, &mpi_argv);
 166   MPI_Comm_dup (MPI_COMM_WORLD, &mpi_ct.comm);
 167   MPI_Comm_size(mpi_ct.comm, &mpi_ct.nproc);
 168   MPI_Comm_rank(mpi_ct.comm, &mpi_ct.my_rank);
 169
 170   if (mpi_ct.nproc > MPI_MAX_PROCESS) {
 171     sys_error(ERR_FATAL, "Number of mpi processes (%d) exceeds max processes (%d)",
 172               mpi_ct.nproc, MPI_MAX_PROCESS);
 173     exit(1);
 174   }
 175 #endif
 176
 177 #ifdef MPI_CT
 178   time_start = MPI_Wtime();
 179 #else
 180   time_start = td_current_sec();
 181 #endif
 182
 183 #ifdef MPI_CT
 184   if (mpi_ct.my_rank == 0) {
 185 #endif
 186     while (1) {
 187       int c = getopt_long(argc, argv, "", my_options, NULL);
 188       char *endptr = NULL;
 189
 190       if (c == -1)
 191         break;
 192
 193       switch (c)
 194         {
 195         case O_INTERP:
 196           if ((opt_interp = opt_set_interpolation(optarg)) < 0) {
 197             ctrec_usage(argv[0]);
 198             return (1);
 199           }
 200           break;
 201         case O_FILTER:
 202           if ((opt_filter = opt_set_filter(optarg)) < 0) {
 203             ctrec_usage(argv[0]);
 204             return (1);
 205           }
 206           break;
 207         case O_BACKPROJ:
 208           if ((opt_backproj = opt_set_backproj(optarg)) < 0) {
 209               ctrec_usage(argv[0]);
 210               return (1);
 211           }
 212           break;
 213         case O_FILTER_PARAM:
 214           opt_filter_param = strtod(optarg, &endptr);
 215           if (endptr != optarg + strlen(optarg)) {
 216             ctrec_usage(argv[0]);
 217           }
 218           break;
 219         case O_VERBOSE:
 220           opt_verbose = 1;
 221           break;
 222         case O_DEBUG:
 223           opt_debug = 1;
 224           break;
 225         case O_TRACE:
 226           if ((opt_trace = opt_set_trace(optarg)) < 0) {
 227             ctrec_usage(argv[0]);
 228             return (1);
 229           }
 230           break;
 231         case O_VERSION:
 232 #ifdef VERSION
 233           fprintf(stdout, "Version %s\n", VERSION);
 234 #else
 235           fprintf(stderr, "Unknown version number");
 236 #endif
 237           exit(0);
 238         case O_HELP:
 239         case '?':
 240           ctrec_usage(argv[0]);
 241           return (0);
 242         default:
 243           ctrec_usage(argv[0]);
 244           return (1);
 245         }
 246     }
 247
 248     if (optind + 4 != argc) {
 249       ctrec_usage(argv[0]);
 250       return (1);
 251     }
 252
 253     rs_name = argv[optind];
 254
 255     im_filename = argv[optind + 1];
 256
 257     nx = strtol(argv[optind + 2], &endptr, 10);
 258     ny = strtol(argv[optind + 3], &endptr, 10);
 259
 260     if (opt_filter == FILTER_G_HAMMING || opt_filter == FILTER_ABS_G_HAMMING)
 261       snprintf (filt_name, sizeof(filt_name), "%s: alpha = %.2f",
 262                name_of_filter (opt_filter), opt_filter_param);
 263     else
 264       snprintf (filt_name, sizeof(filt_name), "%s", name_of_filter (opt_filter));
 265
 266     snprintf (remark, sizeof(remark), "Reconstruct: %dx%d, %s, %s, %s",
 267              nx, ny, filt_name, name_of_interpolation (opt_interp), name_of_backproj(opt_backproj));
 268
 269     if (opt_verbose)
 270       fprintf (stdout, "%s\n", remark);
 271 #ifdef MPI_CT
 272   }
 273 #endif
 274
 275 #ifdef MPI_CT
 276   if (mpi_ct.my_rank == 0) {
 277     rs_global = raysum_open (rs_name);
 278     raysum_read (rs_global);
 279     if (opt_verbose)
 280       print_raysum_info(rs_global);
 281
 282     mpi_ndet = rs_global->ndet;
 283     mpi_nview = rs_global->nview;
 284     mpi_detinc = rs_global->det_inc;
 285     mpi_phmlen = rs_global->phmlen;
 286     mpi_rotinc = rs_global->rot_inc;
 287   }
 288
 289   mpi_t1 = MPI_Wtime();
 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   MPI_Bcast(&opt_filter, 1, MPI_INT, 0, mpi_ct.comm);
 294   MPI_Bcast(&opt_interp, 1, MPI_INT, 0, mpi_ct.comm);
 295   MPI_Bcast(&opt_filter_param, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 296   MPI_Bcast(&opt_interp_param, 1, MPI_INT, 0, mpi_ct.comm);
 297   MPI_Bcast(&opt_backproj, 1, MPI_INT, 0, mpi_ct.comm);
 298   MPI_Bcast(&mpi_ndet, 1, MPI_INT, 0, mpi_ct.comm);
 299   MPI_Bcast(&mpi_nview, 1, MPI_INT, 0, mpi_ct.comm);
 300   MPI_Bcast(&mpi_detinc, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 301   MPI_Bcast(&mpi_phmlen, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 302   MPI_Bcast(&mpi_rotinc, 1, MPI_DOUBLE, 0, mpi_ct.comm);
 303   MPI_Bcast(&nx, 1, MPI_INT, 0, mpi_ct.comm);
 304   MPI_Bcast(&ny, 1, MPI_INT, 0, mpi_ct.comm);
 305   if (opt_verbose) {
 306     mpi_t2 = MPI_Wtime();
 307     mpi_t = mpi_t2 - mpi_t1;
 308     MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 309     if (mpi_ct.my_rank == 0)
 310       printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g);
 311   }
 312
 313   mpi_ct_calc_work_units(mpi_nview);
 314
 315   rs_local = raysum_create (NULL, mpi_ct.local_work_units[mpi_ct.my_rank], mpi_ndet);
 316
 317   rs_local->ndet = mpi_ndet;
 318   rs_local->nview = mpi_nview;
 319   rs_local->det_inc = mpi_detinc;
 320   rs_local->phmlen = mpi_phmlen;
 321   rs_local->rot_inc = mpi_rotinc;
 322
 323   if (opt_verbose)
 324     mpi_t1 = MPI_Wtime();
 325   mpi_scatter_rs(rs_global, rs_local, opt_debug);
 326   if (opt_verbose) {
 327     mpi_t2 = MPI_Wtime();
 328     mpi_t = mpi_t2 - mpi_t1;
 329     MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 330     if (mpi_ct.my_rank == 0)
 331       printf("Time to scatter rs = %f secs, Max time = %f sec\n", mpi_t, mpi_t_g);
 332   }
 333
 334   if (mpi_ct.my_rank == 0) {
 335     im_global = image_create (im_filename, nx, ny);
 336     sdf_add_label (LT_HISTORY, rs_global->remark, rs_global->calctime, im_global->dfp_2d->dfp);
 337     sdf_add_empty_label (im_global->dfp_2d->dfp);
 338   }
 339   im_local = image_create (NULL, nx, ny);
 340
 341 #else
 342   rs_global = raysum_open (rs_name);
 343   raysum_read (rs_global);
 344   if (opt_verbose)
 345     print_raysum_info(rs_global);
 346
 347   im_global = image_create (im_filename, nx, ny);
 348   sdf_add_label (LT_HISTORY, rs_global->remark, rs_global->calctime, im_global->dfp_2d->dfp);
 349   sdf_add_empty_label (im_global->dfp_2d->dfp);
 350 #endif
 351
 352 #ifdef MPI_CT
 353   mpi_t1 = MPI_Wtime();
 354   image_reconst (im_local, rs_local, opt_filter, opt_filter_param,
 355                  opt_interp, opt_interp_param, opt_backproj, opt_trace);
 356   mpi_t2 = MPI_Wtime();
 357   mpi_t = mpi_t2 - mpi_t1;
 358   MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 359   if (mpi_ct.my_rank == 0 && opt_verbose)
 360     printf("Time to reconstruct = %f, Max time = %f\n", mpi_t, mpi_t_g);
 361 #else
 362   image_reconst (im_global, rs_global, opt_filter, opt_filter_param,
 363                  opt_interp, opt_interp_param, opt_backproj, opt_trace);
 364 #endif
 365
 366 #ifdef MPI_CT
 367   if (opt_verbose)
 368     mpi_t1 = MPI_Wtime();
 369
 370   for (ix = 0; ix < im_local->nx; ix++) {
 371     void *recvbuf = NULL;
 372     if (mpi_ct.my_rank == 0)
 373       recvbuf = im_global->v[ix];
 374
 375     MPI_Reduce(im_local->v[ix], recvbuf, im_local->ny, MPI_FLOAT, MPI_SUM, 0, mpi_ct.comm);
 376   }
 377
 378   if (opt_verbose) {
 379     mpi_t2 = MPI_Wtime();
 380     mpi_t = mpi_t2 - mpi_t1;
 381     MPI_Reduce(&mpi_t, &mpi_t_g, 1, MPI_DOUBLE, MPI_MAX, 0, mpi_ct.comm);
 382     if (mpi_ct.my_rank == 0)
 383       printf("Time to reduce image = %f secs, max time = %f\n", mpi_t, mpi_t_g);
 384   }
 385
 386   if (mpi_ct.my_rank == 0)
 387       time_end = MPI_Wtime();
 388 #else
 389   time_end = td_current_sec();
 390 #endif
 391
 392
 393 #ifdef MPI_CT
 394   if (mpi_ct.my_rank == 0)
 395 #endif
 396     {
 397       raysum_close (rs_global);
 398       strncpy (im_global->remark, remark, sizeof(im_global->remark));
 399       im_global->calctime = time_end - time_start;
 400       image_save (im_global);
 401       if (opt_verbose)
 402           fprintf (stdout, "Time active = %.2f\n", im_global->calctime);
 403     }
 404
 405 #ifdef MPI_CT
 406         MPI_Finalize();
 407 #endif
 408
 409         return (0);
 410 }
 411
 412
 413 #ifdef MPI_CT
 414 static void mpi_scatter_rs (RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug)
 415 {
 416   int iproc;
 417   int end_work_unit;
 418   int iw = 0;
 419
 420   if (mpi_ct.my_rank == 0) {
 421     for (iproc = 0; iproc < mpi_ct.nproc; iproc++) {
 422       end_work_unit = mpi_ct.start_work_unit[iproc] + mpi_ct.local_work_units[iproc]  - 1;
 423
 424       for (iw = mpi_ct.start_work_unit[iproc]; iw <= end_work_unit; iw++) {
 425         MPI_Send(&rs_global->view[iw]->ndet, 1, MPI_INT, iproc, 0, mpi_ct.comm);
 426         MPI_Send(&rs_global->view[iw]->view_angle, 1, MPI_DOUBLE, iproc, 0, mpi_ct.comm);
 427         MPI_Send(rs_global->view[iw]->detval, rs_global->ndet, MPI_FLOAT, iproc, 0, mpi_ct.comm);
 428       }
 429     }
 430   }
 431
 432   if (opt_debug)
 433     fprintf(stdout, "Receiving rs data in process %d\n", mpi_ct.my_rank);
 434
 435   end_work_unit = mpi_ct.local_work_units[mpi_ct.my_rank] - 1;
 436   for (iw = 0; iw <= end_work_unit; iw++) {
 437     MPI_Status status;
 438
 439     MPI_Recv(&rs_local->view[iw]->ndet, 1, MPI_INT, 0, 0, mpi_ct.comm, &status);
 440     MPI_Recv(&rs_local->view[iw]->view_angle, 1, MPI_DOUBLE, 0, 0, mpi_ct.comm, &status);
 441     MPI_Recv(rs_local->view[iw]->detval, rs_local->ndet, MPI_FLOAT, 0, 0, mpi_ct.comm, &status);
 442   }
 443   rs_local->nview = mpi_ct.local_work_units[mpi_ct.my_rank];
 444   if (opt_debug) {
 445     MPI_Barrier(MPI_COMM_WORLD);
 446     fprintf(stdout, "Done with mpi_scatter_rs in process %2d\n", mpi_ct.my_rank);
 447   }
 448 }
 449
 450 #endif
 451
 452 static void print_raysum_info(const RAYSUM *rs)
 453 {
 454   printf ("Number of detectors: %d\n", rs->ndet);
 455   printf ("    Number of views: %d\n", rs->nview);
 456   printf ("             Remark: %s\n", rs->remark);
 457   printf ("             phmlen: %f\n", rs->phmlen);
 458   printf ("          det_start: %f\n", rs->det_start);
 459   printf ("            det_inc: %f\n", rs->det_inc);
 460   printf ("          rot_start: %f\n", rs->rot_start);
 461   printf ("            rot_inc: %f\n", rs->rot_inc);
 462 }
 463
 464 #ifndef NO_MAIN
 465 int
 466 main (const int argc, char *const argv[])
 467 {
 468   return (ctrec_main(argc, argv));
 469 }
 470 #endif
 471