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-2001 Kevin Rosenberg
11 ** $Id: reconstruct.cpp,v 1.15 2001/03/01 07:30:49 kevin Exp $
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.
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.
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 ******************************************************************************/
31 * Reconstructor::Reconstructor Reconstruct Image from Projections
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
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
52 Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const char* const filterName,
53 double filt_param, const char* const filterMethodName, const int zeropad,
54 const char* filterGenerationName, const char* const interpName,
55 int interpFactor, const char* const backprojectName, const int iTrace, SGP* pSGP)
56 : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace),
57 m_bFail(false), m_adPlotXAxis(0)
59 m_nFilteredProjections = m_rProj.nDet() * interpFactor;
61 #ifdef HAVE_BSPLINE_INTERP
62 int spline_order = 0, zoom_factor = 0;
63 if (interp_type == I_BSPLINE) {
64 zoom_factor = interpFactor;
67 m_nFilteredProjections = (m_nDet - 1) * (zoom_factor + 1) + 1;
71 double filterBW = 1. / m_rProj.detInc();
72 m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_rProj.detInc(),
73 m_rProj.nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace,
74 m_rProj.geometry(), m_rProj.focalLength(), m_rProj.sourceDetectorLength(), pSGP);
76 if (m_pProcessSignal->fail()) {
78 m_strFailMessage = "Error creating ProcessSignal: ";
79 m_strFailMessage += m_pProcessSignal->failMessage();
80 delete m_pProcessSignal; m_pProcessSignal = NULL;
84 m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor);
85 if (m_pBackprojector->fail()) {
87 m_strFailMessage = "Error creating backprojector: ";
88 m_strFailMessage += m_pBackprojector->failMessage();
89 delete m_pBackprojector; m_pBackprojector = NULL;
90 delete m_pProcessSignal; m_pProcessSignal = NULL;
95 m_adPlotXAxis = new double [m_rProj.nDet()];
96 double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc();
97 double xInc = m_rProj.detInc();
99 for (int i = 0; i < m_rProj.nDet(); i++, x += xInc)
100 m_adPlotXAxis[i] = x;
104 Reconstructor::~Reconstructor ()
106 delete m_pBackprojector;
107 delete m_pProcessSignal;
108 delete m_adPlotXAxis;
113 Reconstructor::plotFilter (SGP* pSGP)
116 int nVecFilter = m_pProcessSignal->getNFilterPoints();
117 double* adPlotXAxis = new double [nVecFilter];
119 if (nVecFilter > 0 && pSGP) {
120 double f = m_pProcessSignal->getFilterMin();
121 double filterInc = m_pProcessSignal->getFilterIncrement();
122 for (int i = 0; i < nVecFilter; i++, f += filterInc)
125 if (m_pProcessSignal->getFilter()) {
128 ezplot.ezset ("title Filter Response");
129 ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter);
139 Reconstructor::reconstructAllViews ()
141 reconstructView (0, m_rProj.nView());
146 Reconstructor::postProcessing()
148 m_pBackprojector->PostProcessing();
153 Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool bBackprojectView, double dGraphWidth)
155 double* adFilteredProj = new double [m_nFilteredProjections]; // filtered projections
158 iViewCount = m_rProj.nView() - iStartView;
160 for (int iView = iStartView; iView < (iStartView + iViewCount); iView++) {
161 if (m_iTrace == Trace::TRACE_CONSOLE)
162 std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n";
164 const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView);
165 const DetectorValue* detval = rDetArray.detValues();
167 m_pProcessSignal->filterSignal (detval, adFilteredProj);
169 #ifdef HAVE_BSPLINE_INTERP
170 if (interp_type == I_BSPLINE)
171 bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
174 if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) {
175 bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj);
176 ezplot_1d (adFilteredProj, m_nFilteredProjections);
181 if (bBackprojectView)
182 m_pBackprojector->BackprojectView (adFilteredProj, rDetArray.viewAngle());
185 if (m_iTrace >= Trace::TRACE_PLOT && pSGP) {
188 std::ostringstream osXLength;
189 osXLength << "xlength " << dGraphWidth;
191 ezplotProj.ezset ("clear");
192 ezplotProj.ezset ("title Raw Projection");
193 ezplotProj.ezset ("xticks major 5");
194 ezplotProj.ezset ("yticks major 5");
195 ezplotProj.ezset ("xlabel ");
196 ezplotProj.ezset ("ylabel ");
197 ezplotProj.ezset ("yporigin 0.55");
198 ezplotProj.ezset ("ylength 0.45");
199 ezplotProj.ezset (osXLength.str().c_str());
200 ezplotProj.ezset ("box.");
201 ezplotProj.ezset ("grid.");
202 #if 0 // workaround c++ optimizer bug, now disabled by using /O1 in code
203 double* pdDetval = new double [m_rProj.nDet()];
204 for (unsigned int id = 0; id < m_rProj.nDet(); id++) {
205 pdDetval[id] = detval[id];
207 ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_rProj.nDet());
210 ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet());
212 pSGP->setTextPointSize (12);
213 ezplotProj.plot (pSGP);
215 ezplotProj.ezset ("clear");
216 ezplotProj.ezset ("title Filtered Projection");
217 ezplotProj.ezset ("xticks major 5");
218 ezplotProj.ezset ("xlabel ");
219 ezplotProj.ezset ("ylabel ");
220 ezplotProj.ezset ("yticks major 5");
221 ezplotProj.ezset ("yporigin 0.10");
222 ezplotProj.ezset ("ylength 0.45");
223 ezplotProj.ezset (osXLength.str().c_str());
224 ezplotProj.ezset ("box");
225 ezplotProj.ezset ("grid");
226 ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj, m_nFilteredProjections);
227 pSGP->setTextPointSize (12);
228 ezplotProj.plot (pSGP);
234 delete adFilteredProj;