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