libctsim/reconstr.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **      Name:         rec_t.c           Image Reconstruction Procedures
   5 **      Programmer:   Kevin Rosenberg
   6 **      Date Started: Aug 84
   7 **
   8 ** GLOBAL FUNCTIONS
   9 **      proj_reconst()          Reconstruct Image from Projections
  10 **
  11 **  This is part of the CTSim program
  12 **  Copyright (C) 1983-2000 Kevin Rosenberg
  13 **
  14 **  $Id: reconstr.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
  15 **
  16 **  This program is free software; you can redistribute it and/or modify
  17 **  it under the terms of the GNU General Public License (version 2) as
  18 **  published by the Free Software Foundation.
  19 **
  20 **  This program is distributed in the hope that it will be useful,
  21 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  22 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  23 **  GNU General Public License for more details.
  24 **
  25 **  You should have received a copy of the GNU General Public License
  26 **  along with this program; if not, write to the Free Software
  27 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  28 ******************************************************************************/
  29
  30 #include "ct.h"
  31
  32
  33 /* NAME
  34  *   proj_reconst                       Reconstruct Image from Projections
  35  *
  36  * SYNOPSIS
  37  *   im = proj_reconst (im, proj, filt_type, filt_param, interp_type)
  38  *   IMAGE *im                          Output image
  39  *   RAYSUM *proj                               Raysum data if in memory
  40  *   int filt_type                      Type of convolution filter to use
  41  *   double filt_param                  Filter specific parameter
  42  *                                      Currently, used only with Hamming filters
  43  *   int interp_type                    Type of interpolation method to use
  44  *
  45  * ALGORITHM
  46  *
  47  *      Calculate one-dimensional filter in spatial domain
  48  *      Allocate & clear (zero) the 2d output image array
  49  *      For each projection view
  50  *          Convolve raysum array with filter
  51  *          Backproject raysums and add (summate) to image array
  52  *      end
  53  */
  54
  55 ImageFile&
  56 proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, const int trace)
  57 {
  58   int ndet = proj.nDet();
  59   int nview = proj.nView();
  60   double det_inc = proj.detInc();
  61   double detlen = (ndet - 1) * det_inc;
  62   double *vec_proj = new double [ndet];   // projection data
  63   int n_filtered_proj = ndet;
  64   double* filtered_proj = new double [n_filtered_proj];   // convolved result
  65
  66 #ifdef HAVE_BSPLINE_INTERP
  67   int spline_order = 0, zoom_factor = 0;
  68   if (interp_type == I_BSPLINE) {
  69     zoom_factor = interp_param;
  70     spline_order = 3;
  71     zoom_factor = 3;
  72     n_filtered_proj = (ndet - 1) * (zoom_factor + 1) + 1;
  73   }
  74 #endif
  75
  76   int n_vec_filter = 2 * ndet - 1;
  77   double filterBW = 1. / det_inc;
  78   double filterMin = -detlen;
  79   double filterMax = detlen;
  80   double* vec_filter = filter_generate (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0);
  81
  82 #if HAVE_SGP
  83   double *plot_xaxis = NULL;                    /* array for plotting */
  84   SGP_ID gid;
  85   plot_xaxis = new double [n_vec_filter];
  86
  87   if (trace > TRACE_TEXT)  {
  88     int i;
  89     double f;
  90     double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
  91     for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
  92       plot_xaxis[i] = f;
  93
  94     gid = ezplot (plot_xaxis, vec_filter, n_vec_filter);
  95     cio_put_str ("Press any key to continue");
  96     cio_kb_getc ();
  97     sgp2_close (gid);
  98   }
  99   if (trace >= TRACE_TEXT) {
 100     printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", proj.nView(), proj.nDet(), proj.detStart(), proj.detInc());
 101   }
 102 #endif  //HAVE_SGP
 103
 104   Backproject* bj = selectBackprojector (backproj_type, proj, im, interp_type);
 105
 106   for (int iview = 0; iview < nview; iview++)  {
 107     if (trace >= TRACE_TEXT)
 108       printf ("Reconstructing view %d (last = %d)\n",  iview, nview - 1);
 109
 110     DetectorArray& darray = proj.getDetectorArray (iview);
 111     DetectorValue* detval = darray.detValues();
 112
 113     for (int j = 0; j < ndet; j++)
 114       vec_proj[j] = detval[j];
 115
 116     for (int j = 0; j < ndet; j++)
 117       filtered_proj[j] = convolve (vec_proj, vec_filter, det_inc, j, ndet, FUNC_BOTH);
 118
 119 #ifdef HAVE_SGP
 120     if (trace >= TRACE_PLOT)  {
 121       ezset  ("clear.");
 122       ezset  ("xticks major 5.");
 123       ezset  ("xlabel ");
 124       ezset  ("ylabel ");
 125       ezset  ("xlength .5.");
 126       ezset  ("box.");
 127       ezset  ("grid.");
 128       ezset  ("ufinish yes.");
 129       ezplot (vec_proj, plot_xaxis, proj.nDet());
 130       ezset  ("clear.");
 131       ezset  ("xticks major 5.");
 132       ezset  ("xlabel ");
 133       ezset  ("ylabel ");
 134       ezset  ("ustart yes.");
 135       ezset  ("xporigin .5.");
 136       ezset  ("xlength .5.");
 137       ezset ("box");
 138       ezset ("grid");
 139       gid = ezplot (filtered_proj, plot_xaxis, proj.nDet());
 140     }
 141 #endif  //HAVE_SGP
 142
 143 #ifdef HAVE_BSPLINE_INTERP
 144     if (interp_type == I_BSPLINE)
 145         bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
 146
 147 #ifdef HAVE_SGP
 148     if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
 149         bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
 150       ezplot_1d (filtered_proj, n_filtered_proj);
 151     }
 152 #endif
 153 #endif
 154
 155     bj->BackprojectView (filtered_proj, darray.viewAngle());
 156
 157 #ifdef HAVE_SGP
 158     if (trace >= TRACE_PLOT) {
 159       char str[256];
 160       printf ("Do you want to exit with current pic (y/n) -- ");
 161       fgets(str, sizeof(str), stdin);
 162       sgp2_close (sgp2_get_active_win());
 163       if (tolower(str[0]) == 'y') {
 164         break;
 165       }
 166     }
 167 #endif  //HAVE_SGP
 168   }
 169
 170   delete bj;
 171   delete vec_proj;
 172   delete filtered_proj;
 173   delete vec_filter;
 174
 175 #ifdef HAVE_SGP
 176   delete plot_xaxis;
 177 #endif
 178
 179   return (im);
 180 }
 181