1 /*****************************************************************************
4 ** Name: rec_t.c Image Reconstruction Procedures
5 ** Programmer: Kevin Rosenberg
6 ** Date Started: Aug 84
9 ** proj_reconst() Reconstruct Image from Projections
11 ** This is part of the CTSim program
12 ** Copyright (C) 1983-2000 Kevin Rosenberg
14 ** $Id: reconstr.cpp,v 1.2 2000/06/20 17:54:51 kevin Exp $
16 ** This program is free software; you can redistribute it and/or modify
17 ** it under the terms of the GNU General Public License (version 2) as
18 ** published by the Free Software Foundation.
20 ** This program is distributed in the hope that it will be useful,
21 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
22 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 ** GNU General Public License for more details.
25 ** You should have received a copy of the GNU General Public License
26 ** along with this program; if not, write to the Free Software
27 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 ******************************************************************************/
34 * proj_reconst Reconstruct Image from Projections
37 * im = proj_reconst (im, proj, filt_type, filt_param, interp_type)
38 * IMAGE *im Output image
39 * RAYSUM *proj Raysum data if in memory
40 * int filt_type Type of convolution filter to use
41 * double filt_param Filter specific parameter
42 * Currently, used only with Hamming filters
43 * int interp_type Type of interpolation method to use
47 * Calculate one-dimensional filter in spatial domain
48 * Allocate & clear (zero) the 2d output image array
49 * For each projection view
50 * Convolve raysum array with filter
51 * Backproject raysums and add (summate) to image array
56 proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, const int trace)
58 int ndet = proj.nDet();
59 int nview = proj.nView();
60 double det_inc = proj.detInc();
61 double detlen = (ndet - 1) * det_inc;
62 int n_filtered_proj = ndet;
63 double filtered_proj [n_filtered_proj]; // convolved result
65 #ifdef HAVE_BSPLINE_INTERP
66 int spline_order = 0, zoom_factor = 0;
67 if (interp_type == I_BSPLINE) {
68 zoom_factor = interp_param;
71 n_filtered_proj = (ndet - 1) * (zoom_factor + 1) + 1;
75 int n_vec_filter = 2 * ndet - 1;
76 double filterBW = 1. / det_inc;
77 double filterMin = -detlen;
78 double filterMax = detlen;
80 SignalFilter filter (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0);
84 double plot_xaxis [n_vec_filter]; // array for plotting
86 if (trace > TRACE_TEXT) {
89 double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
90 for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
93 gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
94 cio_put_str ("Press any key to continue");
98 if (trace >= TRACE_TEXT) {
99 printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", proj.nView(), proj.nDet(), proj.detStart(), proj.detInc());
103 Backproject* bj = selectBackprojector (backproj_type, proj, im, interp_type);
105 for (int iview = 0; iview < nview; iview++) {
106 if (trace >= TRACE_TEXT)
107 printf ("Reconstructing view %d (last = %d)\n", iview, nview - 1);
109 DetectorArray& darray = proj.getDetectorArray (iview);
110 DetectorValue* detval = darray.detValues();
112 for (int j = 0; j < ndet; j++)
113 filtered_proj[j] = filter.convolve (detval, det_inc, j, ndet, FUNC_BOTH);
116 if (trace >= TRACE_PLOT) {
118 ezset ("xticks major 5.");
121 ezset ("xlength .5.");
124 ezset ("ufinish yes.");
125 ezplot (detval, plot_xaxis, proj.nDet());
127 ezset ("xticks major 5.");
130 ezset ("ustart yes.");
131 ezset ("xporigin .5.");
132 ezset ("xlength .5.");
135 gid = ezplot (filtered_proj, plot_xaxis, proj.nDet());
139 #ifdef HAVE_BSPLINE_INTERP
140 if (interp_type == I_BSPLINE)
141 bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
144 if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
145 bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
146 ezplot_1d (filtered_proj, n_filtered_proj);
151 bj->BackprojectView (filtered_proj, darray.viewAngle());
154 if (trace >= TRACE_PLOT) {
156 printf ("Do you want to exit with current pic (y/n) -- ");
157 fgets(str, sizeof(str), stdin);
158 sgp2_close (sgp2_get_active_win());
159 if (tolower(str[0]) == 'y') {