X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=b2d3062c7dbe17934f1e78d4dbf6ceb55d6a0998;hp=30d47185e3cd08120a9c66c826d6946afb105f71;hb=44ba9ce559d2d52cbd7bbea6bcd76242840fd3eb;hpb=595e63c804284d460ce4d032c3848b75bc57186d

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 30d4718..b2d3062 100644
--- a/libctsim/filter.cpp
+++ b/libctsim/filter.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: filter.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**  $Id: filter.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
@@ -46,40 +46,44 @@
  *					ANALYTIC solutions, use numint = 0
  */
 
-double *
-filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint)
+SignalFilter::SignalFilter (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint)
 {
-  double *f = new double [n];
+  m_vecFilter = new double [n];
+  m_filterType = filt_type;
+  m_bw = bw;
+
   double xinc = (xmax - xmin) / (n - 1);
 
-  if (filt_type == FILTER_SHEPP) {
-    double a = 2 * bw;
+  if (m_filterType == FILTER_SHEPP) {
+    double a = 2 * m_bw;
     double c = - 4. / (a * a);
     int center = (n - 1) / 2;
     int sidelen = center;
-    f[center] = 4. / (a * a);
+    m_vecFilter[center] = 4. / (a * a);
     
     for (int i = 1; i <= sidelen; i++ )
-      f [center + i] = f [center - i] = c / (4 * (i * i) - 1);
+      m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
   } else if (domain == D_FREQ) {
     double x;
     int i;
     for (x = xmin, i = 0; i < n; x += xinc, i++)
-      f[i] = filter_frequency_response (filt_type, x, bw, param);
+      m_vecFilter[i] = frequencyResponse (x, param);
   } else if (domain == D_SPATIAL) {
     double x;
     int i;
     for (x = xmin, i = 0; i < n; x += xinc, i++)
       if (numint == 0)
-	f[i] = filter_spatial_response_analytic (filt_type, x, bw, param);
+	m_vecFilter[i] = spatialResponseAnalytic (x, param);
       else
-	f[i] = filter_spatial_response_calc (filt_type, x, bw, param, numint);
+	m_vecFilter[i] = spatialResponseCalc (x, param, numint);
   } else {
     sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain);
-    return (NULL);
   }
-  
-  return (f);
+}
+
+SignalFilter::~SignalFilter (void)
+{
+  delete m_vecFilter;
 }
 
 
@@ -89,21 +93,27 @@ filter_generate (const FilterType filt_type, double bw, double xmin, double xmax
  *  		 		        response
  *
  * SYNOPSIS
- *   y = filter_spatial_response_calc (filt_type, x, bw, param, n)
+ *   y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
  *   double y			Filter's response in spatial domain
  *   int filt_type		Type of filter (definitions in ct.h)
  *   double x			Spatial position to evaluate filter
- *   double bw			Bandwidth of window
+ *   double m_bw			Bandwidth of window
  *   double param		General parameter for various filters
  *   int n			Number of points to calculate integrations
  */
 
 double 
-filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n)
+SignalFilter::spatialResponseCalc (double x, double param, int n) const
+{
+  return (spatialResponseCalc (m_filterType, m_bw, x, param, n));
+}
+
+double 
+SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double param, int n)
 {
   double zmin, zmax;
 
-  if (filt_type == FILTER_TRIANGLE) {
+  if (fType == FILTER_TRIANGLE) {
     zmin = 0;
     zmax = bw;
   } else {
@@ -115,7 +125,7 @@ filter_spatial_response_calc (int filt_type, double x, double bw, double param,
   double z = zmin;
   double q [n];
   for (int i = 0; i < n; i++, z += zinc)
-    q[i] = filter_frequency_response (filt_type, z, bw, param) * cos (TWOPI * z * x);
+    q[i] = frequencyResponse (fType, bw, z, param) * cos (TWOPI * z * x);
   
   double y = 2 * integrateSimpson (zmin, zmax, q, n);
   
@@ -127,21 +137,28 @@ filter_spatial_response_calc (int filt_type, double x, double bw, double param,
  *    filter_frequency_response			Return filter frequency response
  *
  * SYNOPSIS
- *    h = filter_frequency_response (filt_type, u, bw, param)
+ *    h = filter_frequency_response (filt_type, u, m_bw, param)
  *    double h			Filters frequency response at u
  *    int filt_type		Type of filter
  *    double u			Frequency to evaluate filter at
- *    double bw			Bandwidth of filter
+ *    double m_bw			Bandwidth of filter
  *    double param		General input parameter for various filters
  */
 
 double 
-filter_frequency_response (int filt_type, double u, double bw, double param)
+SignalFilter::frequencyResponse (double u, double param) const
+{
+  return frequencyResponse (m_filterType, m_bw, u, param);
+}
+
+
+double 
+SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double param)
 {
   double q;
   double au = fabs (u);
 
-  switch (filt_type) {
+  switch (fType) {
   case FILTER_BANDLIMIT:
     if (au >= bw / 2)
       q = 0.;
@@ -192,9 +209,7 @@ filter_frequency_response (int filt_type, double u, double bw, double param)
     break;
   default:
     q = 0;
-    sys_error (ERR_WARNING,
-	       "Frequency response for filter %d not implemented [filter_frequency_response]",
-	       filt_type);
+    sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", fType);
     break;
   }
   return (q);
@@ -208,16 +223,22 @@ filter_frequency_response (int filt_type, double u, double bw, double param)
  *				response
  *
  * SYNOPSIS
- *   y = filter_spatial_response_analytic (filt_type, x, bw, param)
+ *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
  *   double y			Filter's response in spatial domain
  *   int filt_type		Type of filter (definitions in ct.h)
  *   double x			Spatial position to evaluate filter
- *   double bw			Bandwidth of window
+ *   double m_bw			Bandwidth of window
  *   double param		General parameter for various filters
  */
 
 double 
-filter_spatial_response_analytic (int filt_type, double x, double bw, double param)
+SignalFilter::spatialResponseAnalytic (double x, double param) const
+{
+  return spatialResponseAnalytic (m_filterType, m_bw, x, param);
+}
+
+double 
+SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, double param)
 {
   double q, temp;
   double u = TWOPI * x;
@@ -225,7 +246,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par
   double b = PI / bw;
   double b2 = TWOPI / bw;
 
-  switch (filt_type) {
+  switch (fType) {
   case FILTER_BANDLIMIT:
     q = bw * sinc(u * w, 1.0);
     break;
@@ -237,8 +258,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par
     q = sinc(b-u,w) + sinc(b+u,w);
     break;
   case FILTER_G_HAMMING:
-    q = 2 * param * sin(u*w)/u + (1-param) *
-      (sinc(b2-u, w) + sinc(b2+u, w));
+    q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
     break;
   case FILTER_ABS_BANDLIMIT:
     q = 2 * integral_abscos (u, w);
@@ -264,9 +284,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par
     break;
   case FILTER_ABS_SINC:
   default:
-    sys_error (ERR_WARNING,
-	       "Analytic filter type %d not implemented [filter_spatial_response_analytic]",
-	       filt_type);
+    sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", fType);
     q = 0;
     break;
   }
@@ -287,12 +305,6 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par
  *   v = sin(x * mult) / x;
  */
 
-double 
-sinc (double x, double mult)
-{
-  return (fabs(x) > F_EPSILON ? (sin (x * mult) / x) : 1.0);
-}
-
 
 /* NAME
  *   integral_abscos			Returns integral of u*cos(u)
@@ -308,7 +320,7 @@ sinc (double x, double mult)
  */
 
 double 
-integral_abscos (double u, double w)
+SignalFilter::integral_abscos (double u, double w)
 {
   if (fabs (u) > F_EPSILON)
     return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
@@ -317,3 +329,96 @@ integral_abscos (double u, double w)
 }
 
 
+/* NAME
+ *    convolve			Discrete convolution of two functions
+ *
+ * SYNOPSIS
+ *    r = convolve (f1, f2, dx, n, np, func_type)
+ *    double r			Convolved result
+ *    double f1[], f2[]		Functions to be convolved
+ *    double dx			Difference between successive x values
+ *    int n			Array index to center convolution about
+ *    int np			Number of points in f1 array
+ *    int func_type		EVEN or ODD or EVEN_AND_ODD function f2
+ *
+ * NOTES
+ *    f1 is the projection data, its indices range from 0 to np - 1.
+ *    The index for f2, the filter, ranges from -(np-1) to (np-1).
+ *    There are 3 ways to handle the negative vertices of f2:
+ *	1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
+ *	   All filters used in reconstruction are even.
+ *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n] 
+ *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
+ *	   for negative indices.  Since f2 must range from -(np-1) to (np-1),
+ *	   if we add (np - 1) to f2's array index, then f2's index will
+ *	   range from 0 to 2 * (np - 1), and the origin, x = 0, will be
+ *	   stored at f2[np-1].
+ */
+
+double 
+SignalFilter::convolve (const double func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+{
+  double sum = 0.0;
+
+  if (func_type == FUNC_BOTH) {
+#if UNOPTIMIZED_CONVOLUTION
+    for (int i = 0; i < np; i++)
+      sum += func[i] * m_vecFilter[n - i + (np - 1)];
+#else
+    double* f2 = m_vecFilter + n + (np - 1);
+    for (int i = 0; i < np; i++)
+      sum += *func++ * *f2--;
+#endif
+  }  else if (func_type == FUNC_EVEN) {
+    for (int i = 0; i < np; i++) {
+      int k = abs (n - i);
+      sum += func[i] * m_vecFilter[k];
+    }
+  } else if (func_type == FUNC_ODD) {
+    for (int i = 0; i < np; i++) {
+      int k = n - i;
+      if (k < 0)
+	sum -= func[i] * m_vecFilter[k];
+      else
+	sum += func[i] * m_vecFilter[k];
+    }
+  } else
+    sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
+
+  return (sum * dx);
+}
+
+
+double 
+SignalFilter::convolve (const float func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+{
+  double sum = 0.0;
+
+  if (func_type == FUNC_BOTH) {
+#if UNOPTIMIZED_CONVOLUTION
+    for (int i = 0; i < np; i++)
+      sum += func[i] * m_vecFilter[n - i + (np - 1)];
+#else
+    double* f2 = m_vecFilter + n + (np - 1);
+    for (int i = 0; i < np; i++)
+      sum += *func++ * *f2--;
+#endif
+  }  else if (func_type == FUNC_EVEN) {
+    for (int i = 0; i < np; i++) {
+      int k = abs (n - i);
+      sum += func[i] * m_vecFilter[k];
+    }
+  } else if (func_type == FUNC_ODD) {
+    for (int i = 0; i < np; i++) {
+      int k = n - i;
+      if (k < 0)
+	sum -= func[i] * m_vecFilter[k];
+      else
+	sum += func[i] * m_vecFilter[k];
+    }
+  } else
+    sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
+
+  return (sum * dx);
+}
+