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.13 2001/02/20 04:48:45 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(), 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   delete m_pBackprojector; m_pBackprojector = NULL;
 143 }
 144
 145
 146 void
 147 Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth)
 148 {
 149   double* adFilteredProj = new double [m_nFilteredProjections];   // filtered projections
 150
 151   if (iViewCount <= 0)
 152     iViewCount = m_rProj.nView() - iStartView;
 153
 154   for (int iView = iStartView; iView < (iStartView + iViewCount); iView++)  {
 155     if (m_iTrace == Trace::TRACE_CONSOLE)
 156                 std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
 157
 158     const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
 159     const DetectorValue* detval = rDetArray.detValues();
 160
 161     m_pProcessSignal->filterSignal (detval, adFilteredProj);
 162
 163 #ifdef HAVE_BSPLINE_INTERP
 164     if (interp_type == I_BSPLINE)
 165         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 166
 167 #ifdef HAVE_SGP
 168     if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
 169         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 170       ezplot_1d (adFilteredProj, m_nFilteredProjections);
 171     }
 172 #endif
 173 #endif
 174
 175         if (bBackprojectView)
 176       m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
 177
 178 #ifdef HAVE_SGP
 179     if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
 180       EZPlot ezplotProj;
 181
 182       std::ostringstream osXLength;
 183       osXLength << "xlength " << dGraphWidth;
 184
 185       ezplotProj.ezset ("clear");
 186       ezplotProj.ezset ("title Raw Projection");
 187       ezplotProj.ezset ("xticks major 5");
 188       ezplotProj.ezset ("yticks major 5");
 189       ezplotProj.ezset ("xlabel ");
 190       ezplotProj.ezset ("ylabel ");
 191       ezplotProj.ezset ("yporigin 0.55");
 192       ezplotProj.ezset ("ylength 0.45");
 193       ezplotProj.ezset (osXLength.str().c_str());
 194       ezplotProj.ezset ("box.");
 195       ezplotProj.ezset ("grid.");
 196 #if 0  // workaround c++ optimizer bug, now disabled by using /O1 in code
 197       double* pdDetval = new double [m_rProj.nDet()];
 198       for (unsigned int id = 0; id < m_rProj.nDet(); id++) {
 199         pdDetval[id] = detval[id];
 200       }
 201       ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_rProj.nDet());
 202       delete pdDetval;
 203 #else
 204       ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
 205 #endif
 206       pSGP->setTextPointSize (12);
 207       ezplotProj.plot (pSGP);
 208
 209       ezplotProj.ezset ("clear");
 210       ezplotProj.ezset ("title Filtered Projection");
 211       ezplotProj.ezset ("xticks major 5");
 212       ezplotProj.ezset ("xlabel ");
 213       ezplotProj.ezset ("ylabel ");
 214       ezplotProj.ezset ("yticks major 5");
 215       ezplotProj.ezset ("yporigin 0.10");
 216       ezplotProj.ezset ("ylength 0.45");
 217       ezplotProj.ezset (osXLength.str().c_str());
 218       ezplotProj.ezset ("box");
 219       ezplotProj.ezset ("grid");
 220       ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj,  m_nFilteredProjections);
 221       pSGP->setTextPointSize (12);
 222       ezplotProj.plot (pSGP);
 223
 224 }
 225 #endif  //HAVE_SGP
 226   }
 227
 228   delete adFilteredProj;
 229 }
 230