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