src/ctrec.cpp

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