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.26 2000/06/13 16:20:31 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 ******************************************************************************/
33 static const int POINTS_PER_CIRCLE=90;
48 PElmType type; /* pelm type (box, ellipse, etc) */
49 double atten; /* X-ray attenuation coefficient */
50 double cx,cy; /* center of pelm */
51 double u,v; /* size of pelm */
52 double rot; /* pelm rotation angle (in radians) */
53 double *x, *y; /* ptr to array of points in obj world coord */
54 int pts; /* number of points in outline arrays */
55 double xmin, xmax, ymin, ymax; /* pelm limits */
56 double radius; /* " */
57 GRFMTX_2D p_to_o; /* map from phantom to standard pelm coords */
58 GRFMTX_2D o_to_p; /* map from std pelm coords to phantom coords */
59 class PhmElement *next; /* pointer to next pelm in phantom */
61 typedef class PhmElement PELM;
64 P_PELMS, /* Phantom made of Pelms */
65 P_UNIT_PULSE, /* Special phantom, not made of pelms */
66 P_FILTER /* defined only by this type */
74 PELM *pelm_list; /* pelm linked-list */
76 int n_pelm; /* number of pelms in phantom */
77 double xmin, xmax, ymin, ymax; /* extent of pelms in pelm coordinates */
78 double radius; /* " " */
80 typedef class Phantom PHANTOM;
83 /*----------------------------------------------------------------------*/
85 /*----------------------------------------------------------------------*/
87 /* Ray sums are collected along an array of ndet detectors. The data
88 * for these detectors is stored in the structure DETECTARRAY
91 typedef float DETECT_TYPE;
95 DETECT_TYPE *detval; /* Pointer to array of values recorded by detector */
96 int ndet; /* Number of detectors in array */
97 double view_angle; /* View angle in radians */
99 typedef struct DetectorArray DETARRAY;
112 ScannerGeometry geometry; /* Geometry of detectory */
113 int ndet; /* Number of detectors in array */
114 int nview; /* Number of rotated views */
115 int nsample; /* Number of rays per detector */
116 double detlen; /* Total length of detector array */
117 double rotlen; /* Rotation angle length in radians (norm 2PI) */
118 double det_inc; /* Increment between centers of detectors */
119 double rot_inc; /* Increment in rotation angle between views */
120 double radius; /* Radius of rotation. Distance from */
121 /* center of phm to center of det */
122 double phmlen; /* Maximum Length of phantom or area of interest */
124 double xd1,yd1,xd2,yd2; /* Coordinates of detector endpoints */
125 double xs1,ys1,xs2,ys2; /* Coordinates of source endpoints */
126 double angle; /* Starting angle */
129 typedef class Detector DETECTOR;
135 int fd; /* Disk file descriptor */
136 int file_mode; /* Current file mode (read or write) */
137 int header_size; /* Size of disk file header */
138 int geometry; /* Geometry of scanner */
139 struct DetectorArray **view; /* Pointer to array of detarray_st pointers */
141 char remark[MAXREMARK+1]; /* description of raysum data */
142 double calctime; /* time required to calculate raysums */
144 int ndet; /* number of detectors in array */
145 int nview; /* number of rotated views */
146 double rot_start; /* starting view rotation */
147 double rot_inc; /* angle between rotations */
148 double det_start; /* distance of beginning detector to center */
150 double det_inc; /* increment between detectors */
151 double phmlen; /* Length of PHANTOM edge (phm is square) */
153 typedef class Projections RAYSUM;
156 /*----------------------------------------------------------------------*/
158 /*----------------------------------------------------------------------*/
160 /* Codes for Coordinate Types */
161 /* Defines coords for pelm_is_point_inside() */
164 PELM_COORD, /* Normalized Pelm Coordinates */
165 PHM_COORD /* User's phantom Coordinates */
168 /* Codes for Filter types */
170 typedef enum { /* filter types for filter_generate() */
176 FILTER_ABS_BANDLIMIT, /* filter times |x| */
178 FILTER_ABS_G_HAMMING,
183 /* function domains */
185 static const char D_FREQ_STR[]= "freq";
186 static const char D_SPATIAL_STR[]= "spatial";
194 FUNC_EVEN = 1, /* function types, f[-n] = f[n] */
195 FUNC_ODD, /* f[-n] = -f[n] */
196 FUNC_BOTH /* function has both odd & even components */
199 /* interpolation methods */
200 #undef HAVE_BSPLINE_INTERP
201 typedef enum { /* Interpolation methods */
202 I_NEAREST = 1, /* Nearest neighbor */
203 #if HAVE_BSPLINE_INTERP
205 I_1BSPLINE, /* 1st order B-Spline */
209 I_LINEAR /* Linear interpolation */
212 /* Constants for sizing PHANTOM */
214 static const double PERCENT_PHM_SIZE_INCR=0.0; /* Fractional increase in phantom limits compared to pelm size */
215 static const int N_EXTRA_DETECTORS=4; /* Number of extra detectors widths when calculating detlen */
217 static const char O_TRACE_NONE_STR[]= "none";
218 static const char O_TRACE_TEXT_STR[]= "text";
219 static const char O_TRACE_PHM_STR[]= "phm";
220 static const char O_TRACE_RAYS_STR[]= "rays";
221 static const char O_TRACE_PLOT_STR[]= "plot";
222 static const char O_TRACE_CLIPPING_STR[]= "clipping";
225 TRACE_NONE, /* No tracing */
226 TRACE_TEXT, /* Minimal status */
227 TRACE_PHM, /* Show phantom */
228 TRACE_RAYS, /* Show all rays */
229 TRACE_PLOT, /* Plot raysums */
230 TRACE_CLIPPING /* Plot clipping */
234 O_PHM_HERMAN, /* Herman head phantom */
235 O_PHM_ROWLAND, /* Rowland head phantom */
236 O_PHM_BROWLAND, /* Bordered Rowland head phantom */
237 O_PHM_UNITPULSE /* Unit pulse phantom */
240 static const char O_PHM_HERMAN_STR[]= "herman";
241 static const char O_PHM_ROWLAND_STR[]= "rowland";
242 static const char O_PHM_BROWLAND_STR[]= "browland";
243 static const char O_PHM_UNITPULSE_STR[]= "unitpulse";
245 static const char O_INTERP_NEAREST_STR[]= "nearest";
246 static const char O_INTERP_LINEAR_STR[]= "linear";
247 static const char O_INTERP_BSPLINE_STR[]= "bspline";
249 static const char O_FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit";
250 static const char O_FILTER_ABS_SINC_STR[]= "abs_sinc";
251 static const char O_FILTER_ABS_COS_STR[]= "abs_cos";
252 static const char O_FILTER_ABS_HAMMING_STR[]= "abs_hamming";
253 static const char O_FILTER_SHEPP_STR[]= "shepp";
254 static const char O_FILTER_BANDLIMIT_STR[]= "bandlimit";
255 static const char O_FILTER_SINC_STR[]= "sinc";
256 static const char O_FILTER_COS_STR[]= "cos";
257 static const char O_FILTER_HAMMING_STR[]= "hamming";
258 static const char O_FILTER_TRIANGLE_STR[]= "triangle";
268 static const char O_BPROJ_TRIG_STR[]= "trig";
269 static const char O_BPROJ_TABLE_STR[]= "table";
270 static const char O_BPROJ_DIFF_STR[]= "diff";
271 static const char O_BPROJ_DIFF2_STR[]= "diff2";
272 static const char O_BPROJ_IDIFF2_STR[]= "idiff2";
274 const static int RAYSUM_TRACE_ROW_TITLE=1;
275 const static int RAYSUM_TRACE_ROW_TITLE2=2;
276 const static int RAYSUM_TRACE_ROW_PHANT_ID=4;
277 const static int RAYSUM_TRACE_ROW_CHROMATIC=7;
278 const static int RAYSUM_TRACE_ROW_SCATTER=8;
279 const static int RAYSUM_TRACE_ROW_PHOT_STAT=9;
280 const static int RAYSUM_TRACE_ROW_NDET=12;
281 const static int RAYSUM_TRACE_ROW_NVIEW=13;
282 const static int RAYSUM_TRACE_ROW_SAMPLES=14;
283 const static int RAYSUM_TRACE_ROW_CURR_VIEW=17;
284 const static int RAYSUM_TRACE_ROW_ATTEN=18;
288 /*************************************************************************
289 * FUNCTION DECLARATIONS
290 ************************************************************************/
293 double convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type);
296 int phm_add_pelm_kb (PHANTOM *phm);
297 PHANTOM *phm_select (void);
298 int interpolation_select (void);
299 int filter_select (double *filter_param);
302 double *filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint);
303 double filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n);
304 double filter_spatial_response_analytic (int filt_type, double x, double bw, double param);
305 double filter_frequency_response (int filt_type, double u, double bw, double param);
306 double sinc (double x, double mult);
307 double integral_abscos(double u, double w);
310 int opt_set_trace(const char *optarg);
311 const char *name_of_phantom(const int phmid);
312 int opt_set_phantom(const char *optarg);
313 InterpolationType opt_set_interpolation(const char *optarg);
314 const char *name_of_interpolation(int interp_type);
315 FilterType opt_set_filter(const char *optarg);
316 const char *name_of_filter(const int filter);
317 DomainType opt_set_filter_domain(const char *optarg);
318 const char *name_of_filter_domain(const DomainType domain);
319 BackprojType opt_set_backproj(const char *optarg);
320 const char *name_of_backproj(const BackprojType backproj);
323 PHANTOM *phm_create(const int phmid);
324 PHANTOM *phm_create_from_file(const char *fname);
325 PHANTOM *phm_init(void);
326 void phm_free (PHANTOM *phm);
327 int phm_add_pelm_file(PHANTOM *phm, const char *fname);
328 void phm_add_pelm (PHANTOM *phm, const char *type, const double cx, const double cy,
329 const double u, const double v, const double rot, const double atten);
330 int pelm_make_points(PELM *obj);
331 void pelm_make_xform (PELM *obj);
332 PELM *pelm_alloc(void);
333 void calc_arc(double x[], double y[], const int pts, const double xcent, const double ycent,
334 const double r, const double start, const double stop);
335 void calc_ellipse(double x[], double y[], const int pts, const double u, const double v);
336 int circle_pts(double theta);
337 void phm_print(PHANTOM *phm);
339 void phm_show(const PHANTOM *phm);
340 void phm_draw(const PHANTOM *phm);
344 void phm_std_herman (PHANTOM *phm);
345 void phm_std_rowland (PHANTOM *phm);
346 void phm_std_rowland_bordered (PHANTOM *phm);
349 int raysum_collect(RAYSUM *rs, const DETECTOR *det, const PHANTOM *phm, const int start_view, const int trace, const int unit_pulse);
350 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);
351 double phm_ray_attenuation (const PHANTOM *phm, const double x1, const double y1, const double x2, const double y2);
352 double pelm_ray_attenuation (PELM *pelm, const double x1, const double y1, const double x2, const double y2);
353 int pelm_clip_line (const PELM *pelm, double& x1, double& y1, double& x2, double& y2);
354 void raysum_trace_show_param (const char *label, const char *fmt, int row, int color, ...);
357 DETECTOR *detector_create(const PHANTOM *phm, const ScannerGeometry geometry, int ndet, int nview, int nsample, const double rot_anglen);
358 void detector_free(DETECTOR *det);
361 RAYSUM *raysum_create(const char *fname, const int nview, const int ndet);
362 RAYSUM *raysum_create_from_det(const char *fname, const DETECTOR *det);
363 RAYSUM *raysum_open(const char *filename);
364 void raysum_alloc_views(RAYSUM *rs);
365 void raysum_free(RAYSUM *rs);
366 int raysum_is_open(RAYSUM *rs);
367 int raysum_close(RAYSUM *rs);
368 int raysum_read_header(RAYSUM *rs);
369 int raysum_write_header(RAYSUM *rs);
370 int raysum_read(RAYSUM *rs);
371 int raysum_write(RAYSUM *rs);
372 void raysum_print(const RAYSUM *rs);
373 void raysum_print_info(const RAYSUM *rs);
374 DETARRAY *detarray_alloc(const int n);
375 void detarray_free(DETARRAY *darray);
376 int detarray_read(RAYSUM *rs, DETARRAY *darray, const int view_num);
377 int detarray_write(RAYSUM *rs, const DETARRAY *darray, const int view_num);
379 /* From phm2image.cpp */
380 void phm_to_imagefile (const PHANTOM *phm, ImageFile& im, const int col_start, const int col_count, const int nsample, const int trace);
381 int pelm_is_point_inside(PELM *obj, const double x, const double y, const CoordType coord_type);
384 void image_filter_response(ImageFile& im, const DomainType domain, double bw, const FilterType filt_type, double filt_param, const int opt_trace);
385 int image_display (const ImageFile& im);
386 int image_display_scale (const ImageFile& im, const int scale, const double pmin, const double pmax);
388 /* From reconstr.cpp */
389 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);
392 void backproj_init (const RAYSUM *rs, ImageFile& im, const BackprojType bproj_method);
393 int backproj_calc (const RAYSUM *rs, ImageFile& im, const double *t, const double view_angle, const int interp_type, const int bproj_method);
394 void backproj_term (const RAYSUM *rs, ImageFile& im, const int bproj_method);
396 void backproj_init_trig (const RAYSUM *rs, ImageFile& im);
397 int backproj_calc_trig (const RAYSUM *rs, ImageFile& im, const double *t,
398 const double view_angle, const int interp_type);
399 void backproj_term_trig (const RAYSUM *rs, ImageFile& im);
400 void backproj_init_table (const RAYSUM *rs, ImageFile& im);
401 int backproj_calc_table (const RAYSUM *rs, ImageFile& im, const double *t,
402 const double view_angle, const int interp_type);
403 void backproj_term_table (const RAYSUM *rs, ImageFile& im);
404 void backproj_init_d (const RAYSUM *rs, ImageFile& im);
405 int backproj_calc_d (const RAYSUM *rs, ImageFile& im, const double *t,
406 const double view_angle, const int interp_type);
407 void backproj_term_d (const RAYSUM *rs, ImageFile& im);
408 void backproj_init_d2 (const RAYSUM *rs, ImageFile& im);
409 int backproj_calc_d2 (const RAYSUM *rs, ImageFile& im, const double *t,
410 const double view_angle, const int interp_type);
411 void backproj_term_d2 (const RAYSUM *rs, ImageFile& im);
412 void backproj_init_id (const RAYSUM *rs, ImageFile& im);
413 int backproj_calc_id (const RAYSUM *rs, ImageFile& im, const double *t,
414 const double view_angle, const int interp_type);
415 void backproj_term_id (const RAYSUM *rs, ImageFile& im);
416 void backproj_init_id2 (const RAYSUM *rs, ImageFile& im);
417 int backproj_calc_id2 (const RAYSUM *rs, ImageFile& im, const double *t,
418 const double view_angle, const int interp_type);
419 void backproj_term_id2 (const RAYSUM *rs, ImageFile& im);