libctsim/reconstruct.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:         reconstruct.cpp         Reconstruction class
   5 **   Programmer:   Kevin Rosenberg
   6 **   Date Started: Aug 84
   7 **
   8 **  This is part of the CTSim program
   9 **  Copyright (c) 1983-2001 Kevin Rosenberg
  10 **
  11 **  $Id: reconstruct.cpp,v 1.16 2001/03/11 15:27:30 kevin Exp $
  12 **
  13 **  This program is free software; you can redistribute it and/or modify
  14 **  it under the terms of the GNU General Public License (version 2) as
  15 **  published by the Free Software Foundation.
  16 **
  17 **  This program is distributed in the hope that it will be useful,
  18 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20 **  GNU General Public License for more details.
  21 **
  22 **  You should have received a copy of the GNU General Public License
  23 **  along with this program; if not, write to the Free Software
  24 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  25 ******************************************************************************/
  26
  27 #include "ct.h"
  28
  29
  30 /* NAME
  31  *   Reconstructor::Reconstructor      Reconstruct Image from Projections
  32  *
  33  * SYNOPSIS
  34  *   im = proj.reconstruct (im, filt_type, filt_param, interp_type)
  35  *   IMAGE *im                  Output image
  36  *   int filt_type              Type of convolution filter to use
  37  *   double filt_param          Filter specific parameter
  38  *                              Currently, used only with Hamming filters
  39  *   int interp_type            Type of interpolation method to use
  40  *
  41  * ALGORITHM
  42  *
  43  *      Calculate one-dimensional filter in spatial domain
  44  *      Allocate & clear (zero) the 2d output image array
  45  *      For each projection view
  46  *          Convolve raysum array with filter
  47  *          Backproject raysums and add (summate) to image array
  48  *      end
  49  */
  50
  51
  52 Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const char* const filterName,
  53                               double filt_param, const char* const filterMethodName, const int zeropad,
  54                               const char* filterGenerationName, const char* const interpName,
  55                               int interpFactor, const char* const backprojectName, const int iTrace,
  56                               ReconstructionROI* pROI, SGP* pSGP)
  57   : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace),
  58     m_bFail(false), m_adPlotXAxis(0)
  59 {
  60   m_nFilteredProjections = m_rProj.nDet() * interpFactor;
  61
  62 #ifdef HAVE_BSPLINE_INTERP
  63   int spline_order = 0, zoom_factor = 0;
  64   if (interp_type == I_BSPLINE) {
  65     zoom_factor = interpFactor;
  66     spline_order = 3;
  67     zoom_factor = 3;
  68     m_nFilteredProjections = (m_nDet - 1) * (zoom_factor + 1) + 1;
  69   }
  70 #endif
  71
  72   double filterBW = 1. / m_rProj.detInc();
  73   m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_rProj.detInc(),
  74     m_rProj.nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace,
  75     m_rProj.geometry(), m_rProj.focalLength(), m_rProj.sourceDetectorLength(), pSGP);
  76
  77   if (m_pProcessSignal->fail()) {
  78     m_bFail = true;
  79     m_strFailMessage = "Error creating ProcessSignal: ";
  80     m_strFailMessage += m_pProcessSignal->failMessage();
  81     delete m_pProcessSignal; m_pProcessSignal = NULL;
  82     return;
  83   }
  84
  85   m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI);
  86   if (m_pBackprojector->fail()) {
  87     m_bFail = true;
  88     m_strFailMessage = "Error creating backprojector: ";
  89     m_strFailMessage += m_pBackprojector->failMessage();
  90     delete m_pBackprojector; m_pBackprojector = NULL;
  91     delete m_pProcessSignal; m_pProcessSignal = NULL;
  92     return;
  93   }
  94
  95 #ifdef HAVE_SGP
  96   m_adPlotXAxis = new double [m_rProj.nDet()];
  97   double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc();
  98   double xInc = m_rProj.detInc();
  99
 100   for (int i = 0; i < m_rProj.nDet(); i++, x += xInc)
 101     m_adPlotXAxis[i] = x;
 102 #endif
 103 }
 104
 105 Reconstructor::~Reconstructor ()
 106 {
 107   delete m_pBackprojector;
 108   delete m_pProcessSignal;
 109   delete m_adPlotXAxis;
 110 }
 111
 112
 113 void
 114 Reconstructor::plotFilter (SGP* pSGP)
 115 {
 116 #ifdef HAVE_SGP
 117   int nVecFilter = m_pProcessSignal->getNFilterPoints();
 118   double* adPlotXAxis = new double [nVecFilter];
 119
 120   if (nVecFilter > 0 && pSGP)  {
 121     double f = m_pProcessSignal->getFilterMin();
 122     double filterInc = m_pProcessSignal->getFilterIncrement();
 123     for (int i = 0; i < nVecFilter; i++, f += filterInc)
 124       adPlotXAxis[i] = f;
 125
 126     if (m_pProcessSignal->getFilter()) {
 127       EZPlot ezplot;
 128
 129       ezplot.ezset ("title Filter Response");
 130       ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter);
 131       ezplot.plot (pSGP);
 132     }
 133   }
 134   delete adPlotXAxis;
 135 #endif
 136 }
 137
 138
 139 void
 140 Reconstructor::reconstructAllViews ()
 141 {
 142   reconstructView (0, m_rProj.nView());
 143   postProcessing();
 144 }
 145
 146 void
 147 Reconstructor::postProcessing()
 148 {
 149   m_pBackprojector->PostProcessing();
 150 }
 151
 152
 153 void
 154 Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth)
 155 {
 156   double* adFilteredProj = new double [m_nFilteredProjections];   // filtered projections
 157
 158   if (iViewCount <= 0)
 159     iViewCount = m_rProj.nView() - iStartView;
 160
 161   for (int iView = iStartView; iView < (iStartView + iViewCount); iView++)  {
 162     if (m_iTrace == Trace::TRACE_CONSOLE)
 163                 std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
 164
 165     const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
 166     const DetectorValue* detval = rDetArray.detValues();
 167
 168     m_pProcessSignal->filterSignal (detval, adFilteredProj);
 169
 170 #ifdef HAVE_BSPLINE_INTERP
 171     if (interp_type == I_BSPLINE)
 172         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 173
 174 #ifdef HAVE_SGP
 175     if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
 176         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 177       ezplot_1d (adFilteredProj, m_nFilteredProjections);
 178     }
 179 #endif
 180 #endif
 181
 182         if (bBackprojectView)
 183       m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
 184
 185 #ifdef HAVE_SGP
 186     if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
 187       EZPlot ezplotProj;
 188
 189       std::ostringstream osXLength;
 190       osXLength << "xlength " << dGraphWidth;
 191
 192       ezplotProj.ezset ("clear");
 193       ezplotProj.ezset ("title Raw Projection");
 194       ezplotProj.ezset ("xticks major 5");
 195       ezplotProj.ezset ("yticks major 5");
 196       ezplotProj.ezset ("xlabel ");
 197       ezplotProj.ezset ("ylabel ");
 198       ezplotProj.ezset ("yporigin 0.55");
 199       ezplotProj.ezset ("ylength 0.45");
 200       ezplotProj.ezset (osXLength.str().c_str());
 201       ezplotProj.ezset ("box.");
 202       ezplotProj.ezset ("grid.");
 203 #if 0  // workaround c++ optimizer bug, now disabled by using /O1 in code
 204       double* pdDetval = new double [m_rProj.nDet()];
 205       for (unsigned int id = 0; id < m_rProj.nDet(); id++) {
 206         pdDetval[id] = detval[id];
 207       }
 208       ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_rProj.nDet());
 209       delete pdDetval;
 210 #else
 211       ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
 212 #endif
 213       pSGP->setTextPointSize (12);
 214       ezplotProj.plot (pSGP);
 215
 216       ezplotProj.ezset ("clear");
 217       ezplotProj.ezset ("title Filtered Projection");
 218       ezplotProj.ezset ("xticks major 5");
 219       ezplotProj.ezset ("xlabel ");
 220       ezplotProj.ezset ("ylabel ");
 221       ezplotProj.ezset ("yticks major 5");
 222       ezplotProj.ezset ("yporigin 0.10");
 223       ezplotProj.ezset ("ylength 0.45");
 224       ezplotProj.ezset (osXLength.str().c_str());
 225       ezplotProj.ezset ("box");
 226       ezplotProj.ezset ("grid");
 227       ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj,  m_nFilteredProjections);
 228       pSGP->setTextPointSize (12);
 229       ezplotProj.plot (pSGP);
 230
 231 }
 232 #endif  //HAVE_SGP
 233   }
 234
 235   delete adFilteredProj;
 236 }
 237