1 /*****************************************************************************
5 ** Purpose: Routines for signal-procesing filters
6 ** Progammer: Kevin Rosenberg
7 ** Date Started: Aug 1984
9 ** This is part of the CTSim program
10 ** Copyright (C) 1983-2000 Kevin Rosenberg
12 ** $Id: filter.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $
14 ** This program is free software; you can redistribute it and/or modify
15 ** it under the terms of the GNU General Public License (version 2) as
16 ** published by the Free Software Foundation.
18 ** This program is distributed in the hope that it will be useful,
19 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 ** GNU General Public License for more details.
23 ** You should have received a copy of the GNU General Public License
24 ** along with this program; if not, write to the Free Software
25 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 ******************************************************************************/
32 * SignalFilter::SignalFilter Construct a signal
35 * f = SignalFilter (filt_type, bw, xmin, xmax, n, param, domain, analytic)
36 * double f Generated filter vector
37 * int filt_type Type of filter wanted
38 * double bw Bandwidth of filter
39 * double xmin, xmax Filter limits
40 * int n Number of points in filter
41 * double param General input parameter to filters
42 * int domain FREQ or SPATIAL domain wanted
43 * int numint Number if intervals for calculating discrete inverse fourier xform
44 * for spatial domain filters. For ANALYTIC solutions, use numint = 0
47 SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numint)
49 m_idFilter = convertFilterNameToID (filterName);
50 m_idDomain = convertDomainNameToID (domainName);
51 init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numint);
54 SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
56 init (filterID, bw, xmin, xmax, n, param, domainID, numint);
60 SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
63 m_idFilter = filterID;
64 m_idDomain = domainID;
65 if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
69 m_nameFilter = convertFilterIDToName (m_idFilter);
70 m_nameDomain = convertDomainIDToName (m_idDomain);
75 m_vecFilter = new double[n];
77 double xinc = (m_xmax - m_xmin) / (m_nPoints - 1);
79 if (m_idFilter == FILTER_SHEPP) {
81 double c = - 4. / (a * a);
82 int center = (m_nPoints - 1) / 2;
84 m_vecFilter[center] = 4. / (a * a);
86 for (int i = 1; i <= sidelen; i++ )
87 m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
88 } else if (m_idDomain == DOMAIN_FREQ) {
91 for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
92 m_vecFilter[i] = frequencyResponse (x, param);
93 } else if (m_idDomain == DOMAIN_SPATIAL) {
96 for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
98 m_vecFilter[i] = spatialResponseAnalytic (x, param);
100 m_vecFilter[i] = spatialResponseCalc (x, param, numint);
102 sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", m_idDomain);
107 SignalFilter::~SignalFilter (void)
113 SignalFilter::FilterID
114 SignalFilter::convertFilterNameToID (const char *filterName)
118 if (strcasecmp (filterName, FILTER_BANDLIMIT_STR) == 0)
119 filterID = FILTER_BANDLIMIT;
120 else if (strcasecmp (filterName, FILTER_HAMMING_STR) == 0)
121 filterID = FILTER_G_HAMMING;
122 else if (strcasecmp (filterName, FILTER_SINC_STR) == 0)
123 filterID = FILTER_SINC;
124 else if (strcasecmp (filterName, FILTER_COS_STR) == 0)
125 filterID = FILTER_COSINE;
126 else if (strcasecmp (filterName, FILTER_TRIANGLE_STR) == 0)
127 filterID = FILTER_TRIANGLE;
128 else if (strcasecmp (filterName, FILTER_ABS_BANDLIMIT_STR) == 0)
129 filterID = FILTER_ABS_BANDLIMIT;
130 else if (strcasecmp (filterName, FILTER_ABS_HAMMING_STR) == 0)
131 filterID = FILTER_ABS_G_HAMMING;
132 else if (strcasecmp (filterName, FILTER_ABS_SINC_STR) == 0)
133 filterID = FILTER_ABS_SINC;
134 else if (strcasecmp (filterName, FILTER_ABS_COS_STR) == 0)
135 filterID = FILTER_ABS_COSINE;
136 else if (strcasecmp (filterName, FILTER_SHEPP_STR) == 0)
137 filterID = FILTER_SHEPP;
139 sys_error(ERR_WARNING, "Invalid filter type %s\n", filterName);
140 filterID = FILTER_INVALID;
147 SignalFilter::convertFilterIDToName (const FilterID filterID)
149 const char *name = "";
151 if (filterID == FILTER_SHEPP)
152 name = FILTER_SHEPP_STR;
153 else if (filterID == FILTER_ABS_COSINE)
154 name = FILTER_ABS_COS_STR;
155 else if (filterID == FILTER_ABS_SINC)
156 name = FILTER_ABS_SINC_STR;
157 else if (filterID == FILTER_ABS_G_HAMMING)
158 name = FILTER_ABS_HAMMING_STR;
159 else if (filterID == FILTER_ABS_BANDLIMIT)
160 name = FILTER_ABS_BANDLIMIT_STR;
161 else if (filterID == FILTER_COSINE)
162 name = FILTER_COS_STR;
163 else if (filterID == FILTER_SINC)
164 name = FILTER_SINC_STR;
165 else if (filterID == FILTER_G_HAMMING)
166 name = FILTER_HAMMING_STR;
167 else if (filterID == FILTER_BANDLIMIT)
168 name = FILTER_BANDLIMIT_STR;
169 else if (filterID == FILTER_TRIANGLE)
170 name = FILTER_TRIANGLE_STR;
175 const SignalFilter::DomainID
176 SignalFilter::convertDomainNameToID (const char* const domainName)
180 if (strcasecmp (domainName, DOMAIN_SPATIAL_STR) == 0)
181 dID = DOMAIN_SPATIAL;
182 else if (strcasecmp (domainName, DOMAIN_FREQ_STR) == 0)
185 dID = DOMAIN_INVALID;
191 SignalFilter::convertDomainIDToName (const DomainID domain)
193 const char *name = "";
195 if (domain == DOMAIN_SPATIAL)
196 return (DOMAIN_SPATIAL_STR);
197 else if (domain == DOMAIN_FREQ)
198 return (DOMAIN_FREQ_STR);
205 SignalFilter::response (const char* filterName, const char* domainName, double bw, double x, double filt_param)
208 FilterID filterID = convertFilterNameToID (filterName);
209 DomainID domainID = convertDomainNameToID (domainName);
211 if (domainID == DOMAIN_SPATIAL)
212 response = spatialResponseAnalytic (filterID, bw, x, filt_param);
213 else if (domainID == DOMAIN_FREQ)
214 response = frequencyResponse (filterID, bw, x, filt_param);
220 * filter_spatial_response_calc Calculate filter by discrete inverse fourier
221 * transform of filters's frequency
225 * y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
226 * double y Filter's response in spatial domain
227 * int filt_type Type of filter (definitions in ct.h)
228 * double x Spatial position to evaluate filter
229 * double m_bw Bandwidth of window
230 * double param General parameter for various filters
231 * int n Number of points to calculate integrations
235 SignalFilter::spatialResponseCalc (double x, double param, int n) const
237 return (spatialResponseCalc (m_idFilter, m_bw, x, param, n));
241 SignalFilter::spatialResponseCalc (FilterID filterID, double bw, double x, double param, int n)
245 if (filterID == FILTER_TRIANGLE) {
252 double zinc = (zmax - zmin) / (n - 1);
256 for (int i = 0; i < n; i++, z += zinc)
257 q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
259 double y = 2 * integrateSimpson (zmin, zmax, q, n);
266 * filter_frequency_response Return filter frequency response
269 * h = filter_frequency_response (filt_type, u, m_bw, param)
270 * double h Filters frequency response at u
271 * int filt_type Type of filter
272 * double u Frequency to evaluate filter at
273 * double m_bw Bandwidth of filter
274 * double param General input parameter for various filters
278 SignalFilter::frequencyResponse (double u, double param) const
280 return frequencyResponse (m_idFilter, m_bw, u, param);
285 SignalFilter::frequencyResponse (FilterID filterID, double bw, double u, double param)
288 double au = fabs (u);
291 case FILTER_BANDLIMIT:
297 case FILTER_ABS_BANDLIMIT:
303 case FILTER_TRIANGLE:
313 q = cos(PI * u / bw);
315 case FILTER_ABS_COSINE:
319 q = au * cos(PI * u / bw);
322 q = bw * sinc (PI * bw * u, 1.);
324 case FILTER_ABS_SINC:
325 q = au * bw * sinc (PI * bw * u, 1.);
327 case FILTER_G_HAMMING:
331 q = param + (1 - param) * cos (TWOPI * u / bw);
333 case FILTER_ABS_G_HAMMING:
337 q = au * (param + (1 - param) * cos(TWOPI * u / bw));
341 sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
350 * filter_spatial_response_analytic Calculate filter by analytic inverse fourier
351 * transform of filters's frequency
355 * y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
356 * double y Filter's response in spatial domain
357 * int filt_type Type of filter (definitions in ct.h)
358 * double x Spatial position to evaluate filter
359 * double m_bw Bandwidth of window
360 * double param General parameter for various filters
364 SignalFilter::spatialResponseAnalytic (double x, double param) const
366 return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
370 SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, double param)
373 double u = TWOPI * x;
376 double b2 = TWOPI / bw;
379 case FILTER_BANDLIMIT:
380 q = bw * sinc(u * w, 1.0);
382 case FILTER_TRIANGLE:
383 temp = sinc (u * w, 1.0);
384 q = bw * temp * temp;
387 q = sinc(b-u,w) + sinc(b+u,w);
389 case FILTER_G_HAMMING:
390 q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
392 case FILTER_ABS_BANDLIMIT:
393 q = 2 * integral_abscos (u, w);
395 case FILTER_ABS_COSINE:
396 q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
398 case FILTER_ABS_G_HAMMING:
399 q = 2 * param * integral_abscos(u,w) +
400 (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
404 q = 4. / (PI * bw * bw);
406 q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
409 if (fabs (x) < bw / 2)
414 case FILTER_ABS_SINC:
416 sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
426 * sinc Return sin(x)/x function
430 * double v sinc value
434 * v = sin(x * mult) / x;
439 * integral_abscos Returns integral of u*cos(u)
442 * q = integral_abscos (u, w)
443 * double q Integral value
444 * double u Integration variable
445 * double w Upper integration boundary
448 * Returns the value of integral of u*cos(u)*dV for V = 0 to w
452 SignalFilter::integral_abscos (double u, double w)
454 if (fabs (u) > F_EPSILON)
455 return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
462 * convolve Discrete convolution of two functions
465 * r = convolve (f1, f2, dx, n, np, func_type)
466 * double r Convolved result
467 * double f1[], f2[] Functions to be convolved
468 * double dx Difference between successive x values
469 * int n Array index to center convolution about
470 * int np Number of points in f1 array
471 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
474 * f1 is the projection data, its indices range from 0 to np - 1.
475 * The index for f2, the filter, ranges from -(np-1) to (np-1).
476 * There are 3 ways to handle the negative vertices of f2:
477 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
478 * All filters used in reconstruction are even.
479 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
480 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
481 * for negative indices. Since f2 must range from -(np-1) to (np-1),
482 * if we add (np - 1) to f2's array index, then f2's index will
483 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
484 * stored at f2[np-1].
488 SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
492 #if UNOPTIMIZED_CONVOLUTION
493 for (int i = 0; i < np; i++)
494 sum += func[i] * m_vecFilter[n - i + (np - 1)];
496 double* f2 = m_vecFilter + n + (np - 1);
497 for (int i = 0; i < np; i++)
498 sum += *func++ * *f2--;
506 SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
510 #if UNOPTIMIZED_CONVOLUTION
511 for (int i = 0; i < np; i++)
512 sum += func[i] * m_vecFilter[n - i + (np - 1)];
514 double* f2 = m_vecFilter + n + (np - 1);
515 for (int i = 0; i < np; i++)
516 sum += *func++ * *f2--;