1 /*****************************************************************************
5 ** Purpose: Master Image Reconstruction Header
6 ** Programmer: Kevin Rosenberg
7 ** Date Started: July 1, 1984
9 ** This is part of the CTSim program
10 ** Copyright (C) 1983-2000 Kevin Rosenberg
12 ** $Id: ir.h,v 1.27 2000/06/15 19:07:10 kevin Exp $
15 ** This program is free software; you can redistribute it and/or modify
16 ** it under the terms of the GNU General Public License (version 2) as
17 ** published by the Free Software Foundation.
19 ** This program is distributed in the hope that it will be useful,
20 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
21 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 ** GNU General Public License for more details.
24 ** You should have received a copy of the GNU General Public License
25 ** along with this program; if not, write to the Free Software
26 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 ******************************************************************************/
36 /*----------------------------------------------------------------------*/
38 /*----------------------------------------------------------------------*/
40 /* Ray sums are collected along an array of ndet detectors. The data
41 * for these detectors is stored in the structure DETECTARRAY
44 typedef float DETECT_TYPE;
48 DETECT_TYPE *detval; /* Pointer to array of values recorded by detector */
49 int ndet; /* Number of detectors in array */
50 double view_angle; /* View angle in radians */
52 typedef struct DetectorArray DETARRAY;
65 ScannerGeometry geometry; /* Geometry of detectory */
66 int ndet; /* Number of detectors in array */
67 int nview; /* Number of rotated views */
68 int nsample; /* Number of rays per detector */
69 double detlen; /* Total length of detector array */
70 double rotlen; /* Rotation angle length in radians (norm 2PI) */
71 double det_inc; /* Increment between centers of detectors */
72 double rot_inc; /* Increment in rotation angle between views */
73 double radius; /* Radius of rotation. Distance from */
74 /* center of phm to center of det */
75 double phmlen; /* Maximum Length of phantom or area of interest */
77 double xd1,yd1,xd2,yd2; /* Coordinates of detector endpoints */
78 double xs1,ys1,xs2,ys2; /* Coordinates of source endpoints */
79 double angle; /* Starting angle */
82 Detector (const Phantom& phm, const ScannerGeometry geometry, int ndet, int nview, int nsample, const double rot_anglen);
91 int fd; /* Disk file descriptor */
92 int file_mode; /* Current file mode (read or write) */
93 int header_size; /* Size of disk file header */
94 int geometry; /* Geometry of scanner */
95 struct DetectorArray **view; /* Pointer to array of detarray_st pointers */
97 string remark; /* description of raysum data */
98 double calctime; /* time required to calculate raysums */
100 int ndet; /* number of detectors in array */
101 int nview; /* number of rotated views */
102 double rot_start; /* starting view rotation */
103 double rot_inc; /* angle between rotations */
104 double det_start; /* distance of beginning detector to center */
106 double det_inc; /* increment between detectors */
107 double phmlen; /* Length of PHANTOM edge (phm is square) */
109 typedef class Projections RAYSUM;
111 /*----------------------------------------------------------------------*/
113 /*----------------------------------------------------------------------*/
115 /* Codes for Filter types */
117 typedef enum { /* filter types for filter_generate() */
123 FILTER_ABS_BANDLIMIT, /* filter times |x| */
125 FILTER_ABS_G_HAMMING,
130 /* function domains */
132 static const char D_FREQ_STR[]= "freq";
133 static const char D_SPATIAL_STR[]= "spatial";
141 FUNC_EVEN = 1, /* function types, f[-n] = f[n] */
142 FUNC_ODD, /* f[-n] = -f[n] */
143 FUNC_BOTH /* function has both odd & even components */
146 /* interpolation methods */
147 #undef HAVE_BSPLINE_INTERP
148 typedef enum { /* Interpolation methods */
149 I_NEAREST = 1, /* Nearest neighbor */
150 #if HAVE_BSPLINE_INTERP
152 I_1BSPLINE, /* 1st order B-Spline */
156 I_LINEAR /* Linear interpolation */
159 static const int N_EXTRA_DETECTORS=4; /* Number of extra detectors widths when calculating detlen */
161 static const char O_TRACE_NONE_STR[]= "none";
162 static const char O_TRACE_TEXT_STR[]= "text";
163 static const char O_TRACE_PHM_STR[]= "phm";
164 static const char O_TRACE_RAYS_STR[]= "rays";
165 static const char O_TRACE_PLOT_STR[]= "plot";
166 static const char O_TRACE_CLIPPING_STR[]= "clipping";
169 TRACE_NONE, /* No tracing */
170 TRACE_TEXT, /* Minimal status */
171 TRACE_PHM, /* Show phantom */
172 TRACE_RAYS, /* Show all rays */
173 TRACE_PLOT, /* Plot raysums */
174 TRACE_CLIPPING /* Plot clipping */
178 O_PHM_HERMAN, /* Herman head phantom */
179 O_PHM_ROWLAND, /* Rowland head phantom */
180 O_PHM_BROWLAND, /* Bordered Rowland head phantom */
181 O_PHM_UNITPULSE /* Unit pulse phantom */
184 static const char O_PHM_HERMAN_STR[]= "herman";
185 static const char O_PHM_ROWLAND_STR[]= "rowland";
186 static const char O_PHM_BROWLAND_STR[]= "browland";
187 static const char O_PHM_UNITPULSE_STR[]= "unitpulse";
189 static const char O_INTERP_NEAREST_STR[]= "nearest";
190 static const char O_INTERP_LINEAR_STR[]= "linear";
191 static const char O_INTERP_BSPLINE_STR[]= "bspline";
193 static const char O_FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit";
194 static const char O_FILTER_ABS_SINC_STR[]= "abs_sinc";
195 static const char O_FILTER_ABS_COS_STR[]= "abs_cos";
196 static const char O_FILTER_ABS_HAMMING_STR[]= "abs_hamming";
197 static const char O_FILTER_SHEPP_STR[]= "shepp";
198 static const char O_FILTER_BANDLIMIT_STR[]= "bandlimit";
199 static const char O_FILTER_SINC_STR[]= "sinc";
200 static const char O_FILTER_COS_STR[]= "cos";
201 static const char O_FILTER_HAMMING_STR[]= "hamming";
202 static const char O_FILTER_TRIANGLE_STR[]= "triangle";
212 static const char O_BPROJ_TRIG_STR[]= "trig";
213 static const char O_BPROJ_TABLE_STR[]= "table";
214 static const char O_BPROJ_DIFF_STR[]= "diff";
215 static const char O_BPROJ_DIFF2_STR[]= "diff2";
216 static const char O_BPROJ_IDIFF2_STR[]= "idiff2";
218 const static int RAYSUM_TRACE_ROW_TITLE=1;
219 const static int RAYSUM_TRACE_ROW_TITLE2=2;
220 const static int RAYSUM_TRACE_ROW_PHANT_ID=4;
221 const static int RAYSUM_TRACE_ROW_CHROMATIC=7;
222 const static int RAYSUM_TRACE_ROW_SCATTER=8;
223 const static int RAYSUM_TRACE_ROW_PHOT_STAT=9;
224 const static int RAYSUM_TRACE_ROW_NDET=12;
225 const static int RAYSUM_TRACE_ROW_NVIEW=13;
226 const static int RAYSUM_TRACE_ROW_SAMPLES=14;
227 const static int RAYSUM_TRACE_ROW_CURR_VIEW=17;
228 const static int RAYSUM_TRACE_ROW_ATTEN=18;
232 /*************************************************************************
233 * FUNCTION DECLARATIONS
234 ************************************************************************/
237 double convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type);
240 bool phm_add_pelem_kb (Phantom& phm);
241 const Phantom& phm_select (Phantom& phm);
242 int interpolation_select (void);
243 int filter_select (double *filter_param);
246 double *filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint);
247 double filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n);
248 double filter_spatial_response_analytic (int filt_type, double x, double bw, double param);
249 double filter_frequency_response (int filt_type, double u, double bw, double param);
250 double sinc (double x, double mult);
251 double integral_abscos(double u, double w);
254 int opt_set_trace(const char *optarg);
255 const char *name_of_phantom(const int phmid);
256 int opt_set_phantom(const char *optarg);
257 InterpolationType opt_set_interpolation(const char *optarg);
258 const char *name_of_interpolation(int interp_type);
259 FilterType opt_set_filter(const char *optarg);
260 const char *name_of_filter(const int filter);
261 DomainType opt_set_filter_domain(const char *optarg);
262 const char *name_of_filter_domain(const DomainType domain);
263 BackprojType opt_set_backproj(const char *optarg);
264 const char *name_of_backproj(const BackprojType backproj);
267 int raysum_collect (RAYSUM *rs, const Detector& det, const Phantom& phm, const int start_view, const int trace, const int unit_pulse);
268 void rayview (const Phantom& phm, DETARRAY *darray, const Detector& det, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, const int unit_pulse);
269 double phm_ray_attenuation (const Phantom& phm, const double x1, const double y1, const double x2, const double y2);
270 double pelem_ray_attenuation (PhantomElement& pelem, const double x1, const double y1, const double x2, const double y2);
271 void raysum_trace_show_param (const char *label, const char *fmt, int row, int color, ...);
274 RAYSUM *raysum_create(const char *fname, const int nview, const int ndet);
275 RAYSUM *raysum_create_from_det(const char *fname, const Detector& det);
276 RAYSUM *raysum_open(const char *filename);
277 void raysum_alloc_views(RAYSUM *rs);
278 void raysum_free(RAYSUM *rs);
279 int raysum_is_open(RAYSUM *rs);
280 int raysum_close(RAYSUM *rs);
281 int raysum_read_header(RAYSUM *rs);
282 int raysum_write_header(RAYSUM *rs);
283 int raysum_read(RAYSUM *rs);
284 int raysum_write(RAYSUM *rs);
285 void raysum_print(const RAYSUM *rs);
286 void raysum_print_info(const RAYSUM *rs);
287 DETARRAY *detarray_alloc(const int n);
288 void detarray_free(DETARRAY *darray);
289 int detarray_read(RAYSUM *rs, DETARRAY *darray, const int view_num);
290 int detarray_write(RAYSUM *rs, const DETARRAY *darray, const int view_num);
292 /* From phm2image.cpp */
293 void phm_to_imagefile (const Phantom& phm, ImageFile& im, const int col_start, const int col_count, const int nsample, const int trace);
296 void image_filter_response(ImageFile& im, const DomainType domain, double bw, const FilterType filt_type, double filt_param, const int opt_trace);
297 int image_display (const ImageFile& im);
298 int image_display_scale (const ImageFile& im, const int scale, const double pmin, const double pmax);
300 /* From reconstr.cpp */
301 ImageFile& proj_reconst (ImageFile& im, RAYSUM *rs, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, int const ir_trace);
304 void backproj_init (const RAYSUM *rs, ImageFile& im, const BackprojType bproj_method);
305 int backproj_calc (const RAYSUM *rs, ImageFile& im, const double *t, const double view_angle, const int interp_type, const int bproj_method);
306 void backproj_term (const RAYSUM *rs, ImageFile& im, const int bproj_method);
308 void backproj_init_trig (const RAYSUM *rs, ImageFile& im);
309 int backproj_calc_trig (const RAYSUM *rs, ImageFile& im, const double *t,
310 const double view_angle, const int interp_type);
311 void backproj_term_trig (const RAYSUM *rs, ImageFile& im);
312 void backproj_init_table (const RAYSUM *rs, ImageFile& im);
313 int backproj_calc_table (const RAYSUM *rs, ImageFile& im, const double *t,
314 const double view_angle, const int interp_type);
315 void backproj_term_table (const RAYSUM *rs, ImageFile& im);
316 void backproj_init_d (const RAYSUM *rs, ImageFile& im);
317 int backproj_calc_d (const RAYSUM *rs, ImageFile& im, const double *t,
318 const double view_angle, const int interp_type);
319 void backproj_term_d (const RAYSUM *rs, ImageFile& im);
320 void backproj_init_d2 (const RAYSUM *rs, ImageFile& im);
321 int backproj_calc_d2 (const RAYSUM *rs, ImageFile& im, const double *t,
322 const double view_angle, const int interp_type);
323 void backproj_term_d2 (const RAYSUM *rs, ImageFile& im);
324 void backproj_init_id (const RAYSUM *rs, ImageFile& im);
325 int backproj_calc_id (const RAYSUM *rs, ImageFile& im, const double *t,
326 const double view_angle, const int interp_type);
327 void backproj_term_id (const RAYSUM *rs, ImageFile& im);
328 void backproj_init_id2 (const RAYSUM *rs, ImageFile& im);
329 int backproj_calc_id2 (const RAYSUM *rs, ImageFile& im, const double *t,
330 const double view_angle, const int interp_type);
331 void backproj_term_id2 (const RAYSUM *rs, ImageFile& im);