X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Freconstruct.cpp;h=a0d74167f2e45a914f6411fde3041fe195750e7d;hp=c3b6e8b65eb6777cace756dfe7c95544d52a9c70;hb=d84febd85035c7985334618cdda43df21824efbe;hpb=de411914da8b157958e9caae917bf1edeafbb713 diff --git a/libctsim/reconstruct.cpp b/libctsim/reconstruct.cpp index c3b6e8b..a0d7416 100644 --- a/libctsim/reconstruct.cpp +++ b/libctsim/reconstruct.cpp @@ -1,14 +1,12 @@ /***************************************************************************** ** 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 @@ -32,32 +30,34 @@ * * 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; @@ -69,10 +69,10 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch } #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; @@ -82,7 +82,7 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch 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: "; @@ -93,17 +93,20 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch } #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; @@ -139,7 +142,7 @@ Reconstructor::plotFilter (SGP* pSGP) void Reconstructor::reconstructAllViews () { - reconstructView (0, m_rProj.nView()); + reconstructView (0, m_pProj->nView()); postProcessing(); } @@ -156,30 +159,30 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool 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 @@ -201,16 +204,16 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool 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"); @@ -225,10 +228,10 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool 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 }