libctsim/convolve.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **      Name:           Convolve.c              Convolution functions
   5 **      Date Started:   13 Mar 86
   6 **
   7 **  This is part of the CTSim program
   8 **  Copyright (C) 1983-2000 Kevin Rosenberg
   9 **
  10 **  This program is free software; you can redistribute it and/or modify
  11 **  it under the terms of the GNU General Public License (version 2) as
  12 **  published by the Free Software Foundation.
  13 **
  14 **  This program is distributed in the hope that it will be useful,
  15 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 **  GNU General Public License for more details.
  18 **
  19 **  You should have received a copy of the GNU General Public License
  20 **  along with this program; if not, write to the Free Software
  21 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  22 ******************************************************************************/
  23
  24 #include "ct.h"
  25
  26
  27 /* NAME
  28  *    convolve                  Discrete convolution of two functions
  29  *
  30  * SYNOPSIS
  31  *    r = convolve (f1, f2, dx, n, np, func_type)
  32  *    double r                  Convolved result
  33  *    double f1[], f2[]         Functions to be convolved
  34  *    double dx                 Difference between successive x values
  35  *    int n                     Array index to center convolution about
  36  *    int np                    Number of points in f1 array
  37  *    int func_type             EVEN or ODD or EVEN_AND_ODD function f2
  38  *
  39  * NOTES
  40  *    f1 is the projection data, its indices range from 0 to np - 1.
  41  *    The index for f2, the filter, ranges from -(np-1) to (np-1).
  42  *    There are 3 ways to handle the negative vertices of f2:
  43  *      1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
  44  *         All filters used in reconstruction are even.
  45  *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
  46  *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
  47  *         for negative indices.  Since f2 must range from -(np-1) to (np-1),
  48  *         if we add (np - 1) to f2's array index, then f2's index will
  49  *         range from 0 to 2 * (np - 1), and the origin, x = 0, will be
  50  *         stored at f2[np-1].
  51  */
  52
  53 double
  54 convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type)
  55 {
  56   double sum = 0.0;
  57
  58   if (func_type == FUNC_BOTH) {
  59 #if UNOPTIMIZED_CONVOLUTION
  60     for (int i = 0; i < np; i++)
  61       sum += f1[i] * f2[n - i + (np - 1)];
  62 #else
  63     f2 += n + (np - 1);
  64     for (int i = 0; i < np; i++)
  65       sum += *f1++ * *f2--;
  66 #endif
  67   }  else if (func_type == FUNC_EVEN) {
  68     for (int i = 0; i < np; i++) {
  69       int k = abs (n - i);
  70       sum += f1[i] * f2[k];
  71     }
  72   } else if (func_type == FUNC_ODD) {
  73     for (int i = 0; i < np; i++) {
  74       int k = n - i;
  75       if (k < 0)
  76         sum -= f1[i] * f2[k];
  77       else
  78         sum += f1[i] * f2[k];
  79     }
  80   } else
  81     sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
  82
  83   return (sum * dx);
  84 }