X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Freconstruct.cpp;h=a0d74167f2e45a914f6411fde3041fe195750e7d;hp=591aed28282888f9441c2e190b1a5a38093740d5;hb=d84febd85035c7985334618cdda43df21824efbe;hpb=6afa21de8aa00b405de47584efe108c71df33e1b diff --git a/libctsim/reconstruct.cpp b/libctsim/reconstruct.cpp index 591aed2..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-2000 Kevin Rosenberg -** -** $Id: reconstruct.cpp,v 1.3 2000/12/16 02:44:26 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,27 +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, 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; @@ -64,8 +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(), 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; @@ -75,7 +82,7 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch return; } - m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor); + m_pBackprojector = new Backprojector (*m_pProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI); if (m_pBackprojector->fail()) { m_bFail = true; m_strFailMessage = "Error creating backprojector: "; @@ -85,18 +92,21 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch return; } -#if HAVE_SGP - m_adPlotXAxis = new double [m_rProj.nDet()]; - double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc(); - double xInc = m_rProj.detInc(); +#ifdef HAVE_SGP + 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; @@ -106,7 +116,7 @@ Reconstructor::~Reconstructor () void Reconstructor::plotFilter (SGP* pSGP) { -#if HAVE_SGP +#ifdef HAVE_SGP int nVecFilter = m_pProcessSignal->getNFilterPoints(); double* adPlotXAxis = new double [nVecFilter]; @@ -117,13 +127,13 @@ Reconstructor::plotFilter (SGP* pSGP) adPlotXAxis[i] = f; if (m_pProcessSignal->getFilter()) { - EZPlot ezplot (*pSGP); + EZPlot ezplot; ezplot.ezset ("title Filter Response"); ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter); - ezplot.plot(); + ezplot.plot (pSGP); } - } + } delete adPlotXAxis; #endif } @@ -132,71 +142,99 @@ Reconstructor::plotFilter (SGP* pSGP) void Reconstructor::reconstructAllViews () { - reconstructView (0, m_rProj.nView()); - delete m_pBackprojector; m_pBackprojector = NULL; + reconstructView (0, m_pProj->nView()); + postProcessing(); +} + +void +Reconstructor::postProcessing() +{ + m_pBackprojector->PostProcessing(); } void -Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP) +Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth) { 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) - 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 - m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle()); + if (bBackprojectView) + m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle()); #ifdef HAVE_SGP if (m_iTrace >= Trace::TRACE_PLOT && pSGP) { - EZPlot ezplotProj (*pSGP); + EZPlot ezplotProj; + + std::ostringstream osXLength; + osXLength << "xlength " << dGraphWidth; ezplotProj.ezset ("clear"); ezplotProj.ezset ("title Raw Projection"); ezplotProj.ezset ("xticks major 5"); + ezplotProj.ezset ("yticks major 5"); ezplotProj.ezset ("xlabel "); ezplotProj.ezset ("ylabel "); - ezplotProj.ezset ("yporigin .5"); - ezplotProj.ezset ("ylength .5"); + ezplotProj.ezset ("yporigin 0.55"); + ezplotProj.ezset ("ylength 0.45"); + ezplotProj.ezset (osXLength.str().c_str()); ezplotProj.ezset ("box."); ezplotProj.ezset ("grid."); - ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet()); - ezplotProj.plot(); +#if 0 // workaround c++ optimizer bug, now disabled by using /O1 in code + 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_pProj->nDet()); + delete pdDetval; +#else + ezplotProj.addCurve (m_adPlotXAxis, detval, m_pProj->nDet()); +#endif + pSGP->setTextPointSize (9); + ezplotProj.plot (pSGP); + ezplotProj.ezset ("clear"); ezplotProj.ezset ("title Filtered Projection"); ezplotProj.ezset ("xticks major 5"); ezplotProj.ezset ("xlabel "); ezplotProj.ezset ("ylabel "); - ezplotProj.ezset ("ylength .5"); + ezplotProj.ezset ("yticks major 5"); + ezplotProj.ezset ("yporigin 0.10"); + ezplotProj.ezset ("ylength 0.45"); + ezplotProj.ezset (osXLength.str().c_str()); ezplotProj.ezset ("box"); ezplotProj.ezset ("grid"); ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj, m_nFilteredProjections); - ezplotProj.plot(); - } + pSGP->setTextPointSize (9); + ezplotProj.plot (pSGP); + +} #endif //HAVE_SGP } - - delete adFilteredProj; + + delete adFilteredProj; }