1 /*****************************************************************************
4 ** Name: reconstruct.cpp Reconstruction class
5 ** Programmer: Kevin Rosenberg
6 ** Date Started: Aug 84
8 ** This is part of the CTSim program
9 ** Copyright (c) 1983-2009 Kevin Rosenberg
11 ** This program is free software; you can redistribute it and/or modify
12 ** it under the terms of the GNU General Public License (version 2) as
13 ** published by the Free Software Foundation.
15 ** This program is distributed in the hope that it will be useful,
16 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ** GNU General Public License for more details.
20 ** You should have received a copy of the GNU General Public License
21 ** along with this program; if not, write to the Free Software
22 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 ******************************************************************************/
29 * Reconstructor::Reconstructor Reconstruct Image from Projections
32 * im = proj.reconstruct (im, filt_type, filt_param, interp_type)
33 * IMAGE *im Output image
34 * int filt_type Type of convolution filter to use
35 * double filt_param Filter specific parameter
36 * Currently, used only with Hamming filters
37 * int interp_type Type of interpolation method to use
41 * Calculate one-dimensional filter in spatial domain
42 * Allocate & clear (zero) the 2d output image array
43 * For each projection view
44 * Convolve raysum array with filter
45 * Backproject raysums and add (summate) to image array
50 Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const char* const filterName,
51 double filt_param, const char* const filterMethodName, const int zeropad,
52 const char* filterGenerationName, const char* const interpName,
53 int interpFactor, const char* const backprojectName, const int iTrace,
54 ReconstructionROI* pROI, bool bRebinToParallel, SGP* pSGP)
55 : m_rOriginalProj(rProj),
56 m_pProj(bRebinToParallel ? m_rOriginalProj.interpolateToParallel() : &m_rOriginalProj),
57 m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0),
58 m_iTrace(iTrace), m_bRebinToParallel(bRebinToParallel), m_bFail(false), m_adPlotXAxis(0)
60 m_nFilteredProjections = m_pProj->nDet() * interpFactor;
62 #ifdef HAVE_BSPLINE_INTERP
63 int spline_order = 0, zoom_factor = 0;
64 if (interp_type == I_BSPLINE) {
65 zoom_factor = interpFactor;
68 m_nFilteredProjections = (m_nDet - 1) * (zoom_factor + 1) + 1;
72 double filterBW = 1. / m_pProj->detInc();
73 m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_pProj->detInc(),
74 m_pProj->nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace,
75 m_pProj->geometry(), m_pProj->focalLength(), m_pProj->sourceDetectorLength(), pSGP);
77 if (m_pProcessSignal->fail()) {
79 m_strFailMessage = "Error creating ProcessSignal: ";
80 m_strFailMessage += m_pProcessSignal->failMessage();
81 delete m_pProcessSignal; m_pProcessSignal = NULL;
85 m_pBackprojector = new Backprojector (*m_pProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI);
86 if (m_pBackprojector->fail()) {
88 m_strFailMessage = "Error creating backprojector: ";
89 m_strFailMessage += m_pBackprojector->failMessage();
90 delete m_pBackprojector; m_pBackprojector = NULL;
91 delete m_pProcessSignal; m_pProcessSignal = NULL;
96 m_adPlotXAxis = new double [m_pProj->nDet()];
97 double x = - ((m_pProj->nDet() - 1) / 2) * m_pProj->detInc();
98 double xInc = m_pProj->detInc();
100 for (int i = 0; i < m_pProj->nDet(); i++, x += xInc)
101 m_adPlotXAxis[i] = x;
105 Reconstructor::~Reconstructor ()
107 if (m_bRebinToParallel)
110 delete m_pBackprojector;
111 delete m_pProcessSignal;
112 delete m_adPlotXAxis;
117 Reconstructor::plotFilter (SGP* pSGP)
120 int nVecFilter = m_pProcessSignal->getNFilterPoints();
121 double* adPlotXAxis = new double [nVecFilter];
123 if (nVecFilter > 0 && pSGP) {
124 double f = m_pProcessSignal->getFilterMin();
125 double filterInc = m_pProcessSignal->getFilterIncrement();
126 for (int i = 0; i < nVecFilter; i++, f += filterInc)
129 if (m_pProcessSignal->getFilter()) {
132 ezplot.ezset ("title Filter Response");
133 ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter);
143 Reconstructor::reconstructAllViews ()
145 reconstructView (0, m_pProj->nView());
150 Reconstructor::postProcessing()
152 m_pBackprojector->PostProcessing();
157 Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth)
159 double* adFilteredProj = new double [m_nFilteredProjections]; // filtered projections
162 iViewCount = m_pProj->nView() - iStartView;
164 for (int iView = iStartView; iView < (iStartView + iViewCount); iView++) {
165 if (m_iTrace == Trace::TRACE_CONSOLE)
166 std::cout <<"Reconstructing view " << iView << " (last = " << m_pProj->nView() - 1 << ")\n";
168 const DetectorArray& rDetArray = m_pProj->getDetectorArray (iView);
169 const DetectorValue* detval = rDetArray.detValues();
171 m_pProcessSignal->filterSignal (detval, adFilteredProj);
173 #ifdef HAVE_BSPLINE_INTERP
174 if (interp_type == I_BSPLINE)
175 bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
178 if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
179 bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
180 ezplot_1d (adFilteredProj, m_nFilteredProjections);
185 if (bBackprojectView)
186 m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
189 if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
192 std::ostringstream osXLength;
193 osXLength << "xlength " << dGraphWidth;
195 ezplotProj.ezset ("clear");
196 ezplotProj.ezset ("title Raw Projection");
197 ezplotProj.ezset ("xticks major 5");
198 ezplotProj.ezset ("yticks major 5");
199 ezplotProj.ezset ("xlabel ");
200 ezplotProj.ezset ("ylabel ");
201 ezplotProj.ezset ("yporigin 0.55");
202 ezplotProj.ezset ("ylength 0.45");
203 ezplotProj.ezset (osXLength.str().c_str());
204 ezplotProj.ezset ("box.");
205 ezplotProj.ezset ("grid.");
206 #if 0 // workaround c++ optimizer bug, now disabled by using /O1 in code
207 double* pdDetval = new double [m_pProj->nDet()];
208 for (unsigned int id = 0; id < m_pProj->nDet(); id++) {
209 pdDetval[id] = detval[id];
211 ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_pProj->nDet());
214 ezplotProj.addCurve (m_adPlotXAxis, detval, m_pProj->nDet());
216 pSGP->setTextPointSize (9);
217 ezplotProj.plot (pSGP);
219 ezplotProj.ezset ("clear");
220 ezplotProj.ezset ("title Filtered Projection");
221 ezplotProj.ezset ("xticks major 5");
222 ezplotProj.ezset ("xlabel ");
223 ezplotProj.ezset ("ylabel ");
224 ezplotProj.ezset ("yticks major 5");
225 ezplotProj.ezset ("yporigin 0.10");
226 ezplotProj.ezset ("ylength 0.45");
227 ezplotProj.ezset (osXLength.str().c_str());
228 ezplotProj.ezset ("box");
229 ezplotProj.ezset ("grid");
230 ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj, m_nFilteredProjections);
231 pSGP->setTextPointSize (9);
232 ezplotProj.plot (pSGP);
238 delete adFilteredProj;