X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Freconstr.cpp;fp=libctsim%2Freconstr.cpp;h=11b9fc56bcf415c1a4a17120806d4c69455a5591;hb=99dd1d6ed10db1f669a5fe6af71225a50fc0ddfb;hp=0000000000000000000000000000000000000000;hpb=2c61ff85796550481227f2fbec53506a6b5bd365;p=ctsim.git diff --git a/libctsim/reconstr.cpp b/libctsim/reconstr.cpp new file mode 100644 index 0000000..11b9fc5 --- /dev/null +++ b/libctsim/reconstr.cpp @@ -0,0 +1,181 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: rec_t.c Image Reconstruction Procedures +** Programmer: Kevin Rosenberg +** Date Started: Aug 84 +** +** GLOBAL FUNCTIONS +** proj_reconst() Reconstruct Image from Projections +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: reconstr.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $ +** +** 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 +** published by the Free Software Foundation. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program; if not, write to the Free Software +** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +******************************************************************************/ + +#include "ct.h" + + +/* NAME + * proj_reconst Reconstruct Image from Projections + * + * SYNOPSIS + * im = proj_reconst (im, proj, filt_type, filt_param, interp_type) + * IMAGE *im Output image + * RAYSUM *proj Raysum data if in memory + * 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 + * For each projection view + * Convolve raysum array with filter + * Backproject raysums and add (summate) to image array + * end + */ + +ImageFile& +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) +{ + int ndet = proj.nDet(); + int nview = proj.nView(); + double det_inc = proj.detInc(); + double detlen = (ndet - 1) * det_inc; + double *vec_proj = new double [ndet]; // projection data + int n_filtered_proj = ndet; + double* filtered_proj = new double [n_filtered_proj]; // convolved result + +#ifdef HAVE_BSPLINE_INTERP + int spline_order = 0, zoom_factor = 0; + if (interp_type == I_BSPLINE) { + zoom_factor = interp_param; + spline_order = 3; + zoom_factor = 3; + n_filtered_proj = (ndet - 1) * (zoom_factor + 1) + 1; + } +#endif + + int n_vec_filter = 2 * ndet - 1; + double filterBW = 1. / det_inc; + double filterMin = -detlen; + double filterMax = detlen; + double* vec_filter = filter_generate (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0); + +#if HAVE_SGP + double *plot_xaxis = NULL; /* array for plotting */ + SGP_ID gid; + plot_xaxis = new double [n_vec_filter]; + + if (trace > TRACE_TEXT) { + int i; + double f; + double filterInc = (filterMax - filterMin) / (n_vec_filter - 1); + for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc) + plot_xaxis[i] = f; + + gid = ezplot (plot_xaxis, vec_filter, n_vec_filter); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + sgp2_close (gid); + } + if (trace >= TRACE_TEXT) { + printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", proj.nView(), proj.nDet(), proj.detStart(), proj.detInc()); + } +#endif //HAVE_SGP + + Backproject* bj = selectBackprojector (backproj_type, proj, im, interp_type); + + for (int iview = 0; iview < nview; iview++) { + if (trace >= TRACE_TEXT) + printf ("Reconstructing view %d (last = %d)\n", iview, nview - 1); + + DetectorArray& darray = proj.getDetectorArray (iview); + DetectorValue* detval = darray.detValues(); + + for (int j = 0; j < ndet; j++) + vec_proj[j] = detval[j]; + + for (int j = 0; j < ndet; j++) + filtered_proj[j] = convolve (vec_proj, vec_filter, det_inc, j, ndet, FUNC_BOTH); + +#ifdef HAVE_SGP + if (trace >= TRACE_PLOT) { + ezset ("clear."); + ezset ("xticks major 5."); + ezset ("xlabel "); + ezset ("ylabel "); + ezset ("xlength .5."); + ezset ("box."); + ezset ("grid."); + ezset ("ufinish yes."); + ezplot (vec_proj, plot_xaxis, proj.nDet()); + ezset ("clear."); + ezset ("xticks major 5."); + ezset ("xlabel "); + ezset ("ylabel "); + ezset ("ustart yes."); + ezset ("xporigin .5."); + ezset ("xlength .5."); + ezset ("box"); + ezset ("grid"); + gid = ezplot (filtered_proj, plot_xaxis, proj.nDet()); + } +#endif //HAVE_SGP + +#ifdef HAVE_BSPLINE_INTERP + if (interp_type == I_BSPLINE) + bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj); + +#ifdef HAVE_SGP + if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) { + bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj); + ezplot_1d (filtered_proj, n_filtered_proj); + } +#endif +#endif + + bj->BackprojectView (filtered_proj, darray.viewAngle()); + +#ifdef HAVE_SGP + if (trace >= TRACE_PLOT) { + char str[256]; + printf ("Do you want to exit with current pic (y/n) -- "); + fgets(str, sizeof(str), stdin); + sgp2_close (sgp2_get_active_win()); + if (tolower(str[0]) == 'y') { + break; + } + } +#endif //HAVE_SGP + } + + delete bj; + delete vec_proj; + delete filtered_proj; + delete vec_filter; + +#ifdef HAVE_SGP + delete plot_xaxis; +#endif + + return (im); +} +