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.15 2001/03/01 07:30:49 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, SGP* pSGP)
  56   : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace),
  57     m_bFail(false), m_adPlotXAxis(0)
  58 {
  59   m_nFilteredProjections = m_rProj.nDet() * interpFactor;
  60
  61 #ifdef HAVE_BSPLINE_INTERP
  62   int spline_order = 0, zoom_factor = 0;
  63   if (interp_type == I_BSPLINE) {
  64     zoom_factor = interpFactor;
  65     spline_order = 3;
  66     zoom_factor = 3;
  67     m_nFilteredProjections = (m_nDet - 1) * (zoom_factor + 1) + 1;
  68   }
  69 #endif
  70
  71   double filterBW = 1. / m_rProj.detInc();
  72   m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_rProj.detInc(),
  73     m_rProj.nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace,
  74     m_rProj.geometry(), m_rProj.focalLength(), m_rProj.sourceDetectorLength(), pSGP);
  75
  76   if (m_pProcessSignal->fail()) {
  77     m_bFail = true;
  78     m_strFailMessage = "Error creating ProcessSignal: ";
  79     m_strFailMessage += m_pProcessSignal->failMessage();
  80     delete m_pProcessSignal; m_pProcessSignal = NULL;
  81     return;
  82   }
  83
  84   m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor);
  85   if (m_pBackprojector->fail()) {
  86     m_bFail = true;
  87     m_strFailMessage = "Error creating backprojector: ";
  88     m_strFailMessage += m_pBackprojector->failMessage();
  89     delete m_pBackprojector; m_pBackprojector = NULL;
  90     delete m_pProcessSignal; m_pProcessSignal = NULL;
  91     return;
  92   }
  93
  94 #ifdef HAVE_SGP
  95   m_adPlotXAxis = new double [m_rProj.nDet()];
  96   double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc();
  97   double xInc = m_rProj.detInc();
  98
  99   for (int i = 0; i < m_rProj.nDet(); i++, x += xInc)
 100     m_adPlotXAxis[i] = x;
 101 #endif
 102 }
 103
 104 Reconstructor::~Reconstructor ()
 105 {
 106   delete m_pBackprojector;
 107   delete m_pProcessSignal;
 108   delete m_adPlotXAxis;
 109 }
 110
 111
 112 void
 113 Reconstructor::plotFilter (SGP* pSGP)
 114 {
 115 #ifdef HAVE_SGP
 116   int nVecFilter = m_pProcessSignal->getNFilterPoints();
 117   double* adPlotXAxis = new double [nVecFilter];
 118
 119   if (nVecFilter > 0 && pSGP)  {
 120     double f = m_pProcessSignal->getFilterMin();
 121     double filterInc = m_pProcessSignal->getFilterIncrement();
 122     for (int i = 0; i < nVecFilter; i++, f += filterInc)
 123       adPlotXAxis[i] = f;
 124
 125     if (m_pProcessSignal->getFilter()) {
 126       EZPlot ezplot;
 127
 128       ezplot.ezset ("title Filter Response");
 129       ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter);
 130       ezplot.plot (pSGP);
 131     }
 132   }
 133   delete adPlotXAxis;
 134 #endif
 135 }
 136
 137
 138 void
 139 Reconstructor::reconstructAllViews ()
 140 {
 141   reconstructView (0, m_rProj.nView());
 142   postProcessing();
 143 }
 144
 145 void
 146 Reconstructor::postProcessing()
 147 {
 148   m_pBackprojector->PostProcessing();
 149 }
 150
 151
 152 void
 153 Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth)
 154 {
 155   double* adFilteredProj = new double [m_nFilteredProjections];   // filtered projections
 156
 157   if (iViewCount <= 0)
 158     iViewCount = m_rProj.nView() - iStartView;
 159
 160   for (int iView = iStartView; iView < (iStartView + iViewCount); iView++)  {
 161     if (m_iTrace == Trace::TRACE_CONSOLE)
 162                 std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
 163
 164     const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
 165     const DetectorValue* detval = rDetArray.detValues();
 166
 167     m_pProcessSignal->filterSignal (detval, adFilteredProj);
 168
 169 #ifdef HAVE_BSPLINE_INTERP
 170     if (interp_type == I_BSPLINE)
 171         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 172
 173 #ifdef HAVE_SGP
 174     if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
 175         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 176       ezplot_1d (adFilteredProj, m_nFilteredProjections);
 177     }
 178 #endif
 179 #endif
 180
 181         if (bBackprojectView)
 182       m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
 183
 184 #ifdef HAVE_SGP
 185     if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
 186       EZPlot ezplotProj;
 187
 188       std::ostringstream osXLength;
 189       osXLength << "xlength " << dGraphWidth;
 190
 191       ezplotProj.ezset ("clear");
 192       ezplotProj.ezset ("title Raw Projection");
 193       ezplotProj.ezset ("xticks major 5");
 194       ezplotProj.ezset ("yticks major 5");
 195       ezplotProj.ezset ("xlabel ");
 196       ezplotProj.ezset ("ylabel ");
 197       ezplotProj.ezset ("yporigin 0.55");
 198       ezplotProj.ezset ("ylength 0.45");
 199       ezplotProj.ezset (osXLength.str().c_str());
 200       ezplotProj.ezset ("box.");
 201       ezplotProj.ezset ("grid.");
 202 #if 0  // workaround c++ optimizer bug, now disabled by using /O1 in code
 203       double* pdDetval = new double [m_rProj.nDet()];
 204       for (unsigned int id = 0; id < m_rProj.nDet(); id++) {
 205         pdDetval[id] = detval[id];
 206       }
 207       ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_rProj.nDet());
 208       delete pdDetval;
 209 #else
 210       ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
 211 #endif
 212       pSGP->setTextPointSize (12);
 213       ezplotProj.plot (pSGP);
 214
 215       ezplotProj.ezset ("clear");
 216       ezplotProj.ezset ("title Filtered Projection");
 217       ezplotProj.ezset ("xticks major 5");
 218       ezplotProj.ezset ("xlabel ");
 219       ezplotProj.ezset ("ylabel ");
 220       ezplotProj.ezset ("yticks major 5");
 221       ezplotProj.ezset ("yporigin 0.10");
 222       ezplotProj.ezset ("ylength 0.45");
 223       ezplotProj.ezset (osXLength.str().c_str());
 224       ezplotProj.ezset ("box");
 225       ezplotProj.ezset ("grid");
 226       ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj,  m_nFilteredProjections);
 227       pSGP->setTextPointSize (12);
 228       ezplotProj.plot (pSGP);
 229
 230 }
 231 #endif  //HAVE_SGP
 232   }
 233
 234   delete adFilteredProj;
 235 }
 236