1 /*****************************************************************************
4 ** Name: Convolve.c Convolution functions
5 ** Date Started: 13 Mar 86
7 ** This is part of the CTSim program
8 ** Copyright (C) 1983-2000 Kevin Rosenberg
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.
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.
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 ******************************************************************************/
28 * convolve Discrete convolution of two functions
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
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
54 convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type)
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)];
64 for (int i = 0; i < np; i++)
67 } else if (func_type == FUNC_EVEN) {
68 for (int i = 0; i < np; i++) {
72 } else if (func_type == FUNC_ODD) {
73 for (int i = 0; i < np; i++) {
81 sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);