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-2000 Kevin Rosenberg
  10 **
  11 **  $Id: reconstruct.cpp,v 1.6 2000/12/29 15:45:06 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, double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int iTrace, SGP* pSGP)
  53   : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace), m_bFail(false), m_adPlotXAxis(0)
  54 {
  55   m_nFilteredProjections = m_rProj.nDet() * interpFactor;
  56
  57 #ifdef HAVE_BSPLINE_INTERP
  58   int spline_order = 0, zoom_factor = 0;
  59   if (interp_type == I_BSPLINE) {
  60     zoom_factor = interpFactor;
  61     spline_order = 3;
  62     zoom_factor = 3;
  63     m_nFilteredProjections = (m_nDet - 1) * (zoom_factor + 1) + 1;
  64   }
  65 #endif
  66
  67   double filterBW = 1. / m_rProj.detInc();
  68   m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_rProj.detInc(), m_rProj.nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace, m_rProj.geometry(), m_rProj.focalLength(), pSGP);
  69
  70   if (m_pProcessSignal->fail()) {
  71     m_bFail = true;
  72     m_strFailMessage = "Error creating ProcessSignal: ";
  73     m_strFailMessage += m_pProcessSignal->failMessage();
  74     delete m_pProcessSignal; m_pProcessSignal = NULL;
  75     return;
  76   }
  77
  78   m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor);
  79   if (m_pBackprojector->fail()) {
  80     m_bFail = true;
  81     m_strFailMessage = "Error creating backprojector: ";
  82     m_strFailMessage += m_pBackprojector->failMessage();
  83     delete m_pBackprojector; m_pBackprojector = NULL;
  84     delete m_pProcessSignal; m_pProcessSignal = NULL;
  85     return;
  86   }
  87
  88 #if HAVE_SGP
  89   m_adPlotXAxis = new double [m_rProj.nDet()];
  90   double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc();
  91   double xInc = m_rProj.detInc();
  92
  93   for (int i = 0; i < m_rProj.nDet(); i++, x += xInc)
  94     m_adPlotXAxis[i] = x;
  95 #endif
  96 }
  97
  98 Reconstructor::~Reconstructor ()
  99 {
 100   delete m_pBackprojector;
 101   delete m_pProcessSignal;
 102   delete m_adPlotXAxis;
 103 }
 104
 105
 106 void
 107 Reconstructor::plotFilter (SGP* pSGP)
 108 {
 109 #if HAVE_SGP
 110   int nVecFilter = m_pProcessSignal->getNFilterPoints();
 111   double* adPlotXAxis = new double [nVecFilter];
 112
 113   if (nVecFilter > 0 && pSGP)  {
 114     double f = m_pProcessSignal->getFilterMin();
 115     double filterInc = m_pProcessSignal->getFilterIncrement();
 116     for (int i = 0; i < nVecFilter; i++, f += filterInc)
 117       adPlotXAxis[i] = f;
 118
 119     if (m_pProcessSignal->getFilter()) {
 120       EZPlot ezplot;
 121
 122       ezplot.ezset ("title Filter Response");
 123       ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter);
 124       ezplot.plot (pSGP);
 125     }
 126   }\r
 127   delete adPlotXAxis;
 128 #endif
 129 }
 130
 131
 132 void
 133 Reconstructor::reconstructAllViews ()
 134 {
 135   reconstructView (0, m_rProj.nView());
 136   delete m_pBackprojector; m_pBackprojector = NULL;
 137 }
 138
 139
 140 void
 141 Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView)
 142 {
 143   double* adFilteredProj = new double [m_nFilteredProjections];   // filtered projections
 144
 145   if (iViewCount <= 0)
 146     iViewCount = m_rProj.nView() - iStartView;
 147
 148   for (int iView = iStartView; iView < (iStartView + iViewCount); iView++)  {
 149     if (m_iTrace == Trace::TRACE_CONSOLE)
 150                 std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
 151
 152     const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
 153     const DetectorValue* detval = rDetArray.detValues();
 154
 155     m_pProcessSignal->filterSignal (detval, adFilteredProj);
 156
 157 #ifdef HAVE_BSPLINE_INTERP
 158     if (interp_type == I_BSPLINE)
 159         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 160
 161 #ifdef HAVE_SGP
 162     if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
 163         bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
 164       ezplot_1d (adFilteredProj, m_nFilteredProjections);
 165     }
 166 #endif
 167 #endif
 168 \r
 169         if (bBackprojectView)
 170       m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
 171
 172 #ifdef HAVE_SGP
 173     if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
 174       EZPlot ezplotProj;
 175
 176       ezplotProj.ezset ("clear");
 177       ezplotProj.ezset ("title Raw Projection");
 178       ezplotProj.ezset ("xticks major 5");
 179       ezplotProj.ezset ("xlabel ");
 180       ezplotProj.ezset ("ylabel ");
 181       ezplotProj.ezset ("yporigin .5");
 182       ezplotProj.ezset ("ylength .5");
 183       ezplotProj.ezset ("box.");
 184       ezplotProj.ezset ("grid.");
 185       ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
 186       ezplotProj.plot (pSGP);
 187       ezplotProj.ezset ("clear");
 188       ezplotProj.ezset ("title Filtered Projection");
 189       ezplotProj.ezset ("xticks major 5");
 190       ezplotProj.ezset ("xlabel ");
 191       ezplotProj.ezset ("ylabel ");
 192       ezplotProj.ezset ("ylength .5");
 193       ezplotProj.ezset ("box");
 194       ezplotProj.ezset ("grid");
 195       ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj,  m_nFilteredProjections);
 196       ezplotProj.plot (pSGP);
 197     }
 198 #endif  //HAVE_SGP
 199   }
 200 \r
 201   delete adFilteredProj;\r
 202 }
 203