X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Freconstr.cpp;fp=libctsim%2Freconstr.cpp;h=0000000000000000000000000000000000000000;hb=2d39e823ba389fc68e5317c422b55be006094252;hp=9c887800c8f32fcf7a5e5fae1c4982f1f0794836;hpb=a95e41ac40cd2f3a4401d921618604cf33f2a904;p=ctsim.git diff --git a/libctsim/reconstr.cpp b/libctsim/reconstr.cpp deleted file mode 100644 index 9c88780..0000000 --- a/libctsim/reconstr.cpp +++ /dev/null @@ -1,170 +0,0 @@ -/***************************************************************************** -** 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.2 2000/06/20 17:54:51 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; - int n_filtered_proj = ndet; - double filtered_proj [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; - - SignalFilter filter (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0); - -#if HAVE_SGP - SGP_ID gid; - double plot_xaxis [n_vec_filter]; // array for plotting - - 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, filter.getFilter(), 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++) - filtered_proj[j] = filter.convolve (detval, 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 (detval, 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; - - return (im); -} -