/*****************************************************************************
** FILE IDENTIFICATION
**
-** Name: reconstruct.cpp Reconstruction class
+** Name: reconstruct.cpp Reconstruction class
** Programmer: Kevin Rosenberg
** Date Started: Aug 84
**
** This is part of the CTSim program
-** Copyright (c) 1983-2001 Kevin Rosenberg
-**
-** $Id: reconstruct.cpp,v 1.16 2001/03/11 15:27:30 kevin Exp $
+** Copyright (c) 1983-2009 Kevin Rosenberg
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
*
* SYNOPSIS
* im = proj.reconstruct (im, filt_type, filt_param, interp_type)
- * IMAGE *im Output image
- * int filt_type Type of convolution filter to use
- * double filt_param Filter specific parameter
- * Currently, used only with Hamming filters
- * int interp_type Type of interpolation method to use
+ * IMAGE *im Output image
+ * int filt_type Type of convolution filter to use
+ * double filt_param Filter specific parameter
+ * Currently, used only with Hamming filters
+ * int interp_type Type of interpolation method to use
*
* ALGORITHM
*
- * Calculate one-dimensional filter in spatial domain
- * Allocate & clear (zero) the 2d output image array
+ * Calculate one-dimensional filter in spatial domain
+ * Allocate & clear (zero) the 2d output image array
* For each projection view
- * Convolve raysum array with filter
- * Backproject raysums and add (summate) to image array
+ * Convolve raysum array with filter
+ * Backproject raysums and add (summate) to image array
* end
*/
-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,
- ReconstructionROI* pROI, SGP* pSGP)
- : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace),
- m_bFail(false), m_adPlotXAxis(0)
+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,
+ ReconstructionROI* pROI, bool bRebinToParallel, SGP* pSGP)
+ : m_rOriginalProj(rProj),
+ m_pProj(bRebinToParallel ? m_rOriginalProj.interpolateToParallel() : &m_rOriginalProj),
+ m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0),
+ m_iTrace(iTrace), m_bRebinToParallel(bRebinToParallel), m_bFail(false), m_adPlotXAxis(0)
{
- m_nFilteredProjections = m_rProj.nDet() * interpFactor;
+ m_nFilteredProjections = m_pProj->nDet() * interpFactor;
#ifdef HAVE_BSPLINE_INTERP
int spline_order = 0, zoom_factor = 0;
}
#endif
- double filterBW = 1. / m_rProj.detInc();
- 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(), m_rProj.sourceDetectorLength(), pSGP);
+ double filterBW = 1. / m_pProj->detInc();
+ m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_pProj->detInc(),
+ m_pProj->nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace,
+ m_pProj->geometry(), m_pProj->focalLength(), m_pProj->sourceDetectorLength(), pSGP);
if (m_pProcessSignal->fail()) {
m_bFail = true;
return;
}
- m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI);
+ m_pBackprojector = new Backprojector (*m_pProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI);
if (m_pBackprojector->fail()) {
m_bFail = true;
m_strFailMessage = "Error creating backprojector: ";
}
#ifdef HAVE_SGP
- m_adPlotXAxis = new double [m_rProj.nDet()];
- double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc();
- double xInc = m_rProj.detInc();
+ m_adPlotXAxis = new double [m_pProj->nDet()];
+ double x = - ((m_pProj->nDet() - 1) / 2) * m_pProj->detInc();
+ double xInc = m_pProj->detInc();
- for (int i = 0; i < m_rProj.nDet(); i++, x += xInc)
+ for (int i = 0; i < m_pProj->nDet(); i++, x += xInc)
m_adPlotXAxis[i] = x;
#endif
}
Reconstructor::~Reconstructor ()
{
+ if (m_bRebinToParallel)
+ delete m_pProj;
+
delete m_pBackprojector;
delete m_pProcessSignal;
delete m_adPlotXAxis;
void
Reconstructor::reconstructAllViews ()
{
- reconstructView (0, m_rProj.nView());
+ reconstructView (0, m_pProj->nView());
postProcessing();
}
double* adFilteredProj = new double [m_nFilteredProjections]; // filtered projections
if (iViewCount <= 0)
- iViewCount = m_rProj.nView() - iStartView;
-
+ iViewCount = m_pProj->nView() - iStartView;
+
for (int iView = iStartView; iView < (iStartView + iViewCount); iView++) {
- if (m_iTrace == Trace::TRACE_CONSOLE)
- std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
-
- const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
+ if (m_iTrace == Trace::TRACE_CONSOLE)
+ std::cout <<"Reconstructing view " << iView << " (last = " << m_pProj->nView() - 1 << ")\n";
+
+ const DetectorArray& rDetArray = m_pProj->getDetectorArray (iView);
const DetectorValue* detval = rDetArray.detValues();
m_pProcessSignal->filterSignal (detval, adFilteredProj);
#ifdef HAVE_BSPLINE_INTERP
- if (interp_type == I_BSPLINE)
- bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
-
+ if (interp_type == I_BSPLINE)
+ bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
+
#ifdef HAVE_SGP
if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
- bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
+ bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
ezplot_1d (adFilteredProj, m_nFilteredProjections);
}
#endif
#endif
- if (bBackprojectView)
+ if (bBackprojectView)
m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
#ifdef HAVE_SGP
ezplotProj.ezset ("box.");
ezplotProj.ezset ("grid.");
#if 0 // workaround c++ optimizer bug, now disabled by using /O1 in code
- double* pdDetval = new double [m_rProj.nDet()];
- for (unsigned int id = 0; id < m_rProj.nDet(); id++) {
+ double* pdDetval = new double [m_pProj->nDet()];
+ for (unsigned int id = 0; id < m_pProj->nDet(); id++) {
pdDetval[id] = detval[id];
}
- ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_rProj.nDet());
+ ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_pProj->nDet());
delete pdDetval;
#else
- ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
+ ezplotProj.addCurve (m_adPlotXAxis, detval, m_pProj->nDet());
#endif
- pSGP->setTextPointSize (12);
+ pSGP->setTextPointSize (9);
ezplotProj.plot (pSGP);
ezplotProj.ezset ("clear");
ezplotProj.ezset ("box");
ezplotProj.ezset ("grid");
ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj, m_nFilteredProjections);
- pSGP->setTextPointSize (12);
+ pSGP->setTextPointSize (9);
ezplotProj.plot (pSGP);
-}
+}
#endif //HAVE_SGP
}
projects / ctsim.git / blobdiff