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.1 2000/06/19 02:59:34 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 double *vec_proj = new double [ndet]; // projection data
63 int n_filtered_proj = ndet;
64 double* filtered_proj = new double [n_filtered_proj]; // convolved result
66 #ifdef HAVE_BSPLINE_INTERP
67 int spline_order = 0, zoom_factor = 0;
68 if (interp_type == I_BSPLINE) {
69 zoom_factor = interp_param;
72 n_filtered_proj = (ndet - 1) * (zoom_factor + 1) + 1;
76 int n_vec_filter = 2 * ndet - 1;
77 double filterBW = 1. / det_inc;
78 double filterMin = -detlen;
79 double filterMax = detlen;
80 double* vec_filter = filter_generate (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0);
83 double *plot_xaxis = NULL; /* array for plotting */
85 plot_xaxis = new double [n_vec_filter];
87 if (trace > TRACE_TEXT) {
90 double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
91 for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
94 gid = ezplot (plot_xaxis, vec_filter, n_vec_filter);
95 cio_put_str ("Press any key to continue");
99 if (trace >= TRACE_TEXT) {
100 printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", proj.nView(), proj.nDet(), proj.detStart(), proj.detInc());
104 Backproject* bj = selectBackprojector (backproj_type, proj, im, interp_type);
106 for (int iview = 0; iview < nview; iview++) {
107 if (trace >= TRACE_TEXT)
108 printf ("Reconstructing view %d (last = %d)\n", iview, nview - 1);
110 DetectorArray& darray = proj.getDetectorArray (iview);
111 DetectorValue* detval = darray.detValues();
113 for (int j = 0; j < ndet; j++)
114 vec_proj[j] = detval[j];
116 for (int j = 0; j < ndet; j++)
117 filtered_proj[j] = convolve (vec_proj, vec_filter, det_inc, j, ndet, FUNC_BOTH);
120 if (trace >= TRACE_PLOT) {
122 ezset ("xticks major 5.");
125 ezset ("xlength .5.");
128 ezset ("ufinish yes.");
129 ezplot (vec_proj, plot_xaxis, proj.nDet());
131 ezset ("xticks major 5.");
134 ezset ("ustart yes.");
135 ezset ("xporigin .5.");
136 ezset ("xlength .5.");
139 gid = ezplot (filtered_proj, plot_xaxis, proj.nDet());
143 #ifdef HAVE_BSPLINE_INTERP
144 if (interp_type == I_BSPLINE)
145 bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
148 if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
149 bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
150 ezplot_1d (filtered_proj, n_filtered_proj);
155 bj->BackprojectView (filtered_proj, darray.viewAngle());
158 if (trace >= TRACE_PLOT) {
160 printf ("Do you want to exit with current pic (y/n) -- ");
161 fgets(str, sizeof(str), stdin);
162 sgp2_close (sgp2_get_active_win());
163 if (tolower(str[0]) == 'y') {
172 delete filtered_proj;