Last modified: December 2013

URL: http://cxc.harvard.edu/ciao/ahelp/aconvolve.html
AHELP for CIAO 4.11

aconvolve

Context: tools

Synopsis

Convolve an N-dimensional image with a kernel

Syntax

aconvolve  infile outfile kernelspec [writekernel] [kernelfile]
[writefft] [fftroot] [method] [edges] [const] [pad] [center]
[normkernel] [clobber] [verbose]

Description

`aconvolve' convolves an N-dimensional image with a kernel (specified by the parameter 'kernelspec'). The kernel can be specified in one of several flexible ways: as either a file (image), a function, from the library, or specified on the command line. Other than memory restrictions of the machine, there are no limits to the size or dimensionality of either the kernel or image (i.e. the kernel could be larger than the image). One of the most common applications of this tool is smoothing of an image by convolving with a gaussian.

The convolution of a data set with a kernel yields a measure of "how much" the data set looks like the kernel at every location in the data set. Another way to think about the convolution operation is a filtering of the dataset by the kernel.

Two forms of convolution are available:

Frequency-domain convolution is the standard technique: take two functions, Fourier transform them (using the "FFTW" routine from MIT), take their product in Fourier space, and then back transforming the product to the spatial domain to produce the convolved result. When using method="fft", any size image can be used but the algorithm is faster when the data is "padded" in size to take advantage of certain algorithmic optimisations; see the pad parameter.

In spatial-domain convolution `aconvolve' can treat edge-located pixels in one of 5 ways:


Examples

Example 1

aconvolve test_data.fits outfile.fits "lib:gaus(2,5,1,10,10)" method=fft

Convolve (smooth) the image in test_data.fits with a Gaussian that

The convolution is done by the method of FFTing the input arrays, multiplying, and taking the inverse. It will work with any size array but it is much faster if the size is a power of 2.

Example 2

aconvolve test_data.fits outfile.fits "file:conv.fits" edges=wrap pad=no

Convolve the image in test_data.fits with the image in the file conv.fits, wrapping the edges and storing the data in outfile.fits. The input data are not padded.

Example 3

aconvolve test_data.fits outfile.fits "lib:box(2,1,3,3):(1,1)" pad=yes
method=fft writefft=no

Convolve (smooth) the image test_data.fits with a box that is normalized to 1, has two dimensions, and is 3 units long in each direction. The center of the box is a (1,1). The convolution is done by the method of FFTing the input arrays, multiplying, and taking the inverse. The data are padded to the next power of 2**N before doing the convolution. It will work with any size array but it is much faster if the size is a power of 2.

Example 4

aconvolve test_data.fits outfile.fits
"txt:((1,1,1),(1,1,1),(1,1,1)):(1,1)" edges=constant const=0 pad=yes

Same as Example 3. The "box" is implemented from the command line as a text array. The edges are padded with a constant = 0.

Example 5

aconvolve img.fits shifted_img.fits "txt:((1)):(1,1)" meth=slide
edge=wrap

This shifts both the X and Y by +1.

It is also possible to change the offset (kernel origin) of each axis individually, e.g. to shift X by +512 but not change Y (shift=0):

txt:((1)):(512,0)

Parameters

name type ftype def min max reqd autoname
infile file input       yes  
outfile file output       yes yes
kernelspec string         yes  
writekernel boolean   no        
kernelfile file output ./.       yes
writefft boolean   no        
fftroot file output ./.       yes
method string   slide        
edges string   wrap        
const real   0        
pad boolean   no        
center boolean   no        
normkernel string   area        
clobber boolean   no        
verbose integer   0 0 5    

Detailed Parameter Descriptions

Parameter=infile (file required filetype=input)

Input image file.

The input is an image and can have any number of dimensions. The input can have the following data types: "short" (BITPIX=16), "long" (BITPIX=32), "float" (BITPIX=-32), and "double" (BITPIX=-64).

The image can be a virtual image as defined by the datamodel. Thus one could specify a virtual image by using the "bin" syntax like my_file.fits[EVENTS][bin x=1:100:1, y=1:100:1] to specify a 2D image binned on X and Y columns in EVENTS extension of my_file.fits file.

Parameter=outfile (file required filetype=output autoname=yes)

Output image file.

The output file is an image and has the same dimensions as the input image. The output data type is always "float" (BITPIX=-32). The autoname filename suffix is "_conv".

Parameter=kernelspec (string required)

Kernel specification.

The generic syntax for the kernelspec specification is: [key]:[parameters]:[origin]

Where -

Currently supported libraries are:

box

A multi-dimensional array with constant value. The parameters are (D, N, D1, D2,...DD) where:

gaus

A multi-dimensional, non-rotated Gaussian. The parameters are (D, M, N, S1, S2, ... SD) where:

tophat

A two-dimensional non-rotated elliptical top-hat function. The parameters are (D, N, D1, D2) where:

mexhat

A multi-dimensional, mexican hat function.

mexhat = Norm * exp( -0.5 * (x/s1)^2 + (y/s2)^2 + ...)

The parameters are (D, M, N, S1, S2, ... SD) where:

exp

A multi-dimensional exponential function.

exp = Norm * exp(-abs(x/s1)) * exp(-abs(y/s2)) * ...

The parameters are (D, M, N, S1, S2, ... SD) where:

power

A multi-dimensional power law function.

power = Norm * abs(x)^s1 * abs(y)^s2 * ...

The parameters are (D, M, N, S1, S2, ... SD) where:

beta / lorentz

A multi-dimensional beta/lorentzian function.

beta = Norm / ((x^2 +y^2+...)/ s1^2)

The parameters are (D, M, N, S1) where:

sinc

A multi-dimensional sinc function.

sinc= Norm * sin( 2*pi * x^2/s1 * y^2/s1 * ...

The parameters are (D, M, N, S1) where:

cone

A multi-dimensional cone function.

cone = Norm * 1 - sqrt(x^2+y^2+...)  / s1

The parameters are (D, N, S1)

pyramid

A multi-dimensional stepped pyramid function.

pyramid = Norm * 1 - MIN(x,y,...)  / s1

The parameters are (D, N, S1)

sphere

A hemisphere/ball kernel

sphere = sqrt( (x - s1)^2 + ( y - s1 ) ^ + ... )

The parameters are (D, R)

hanning

Hann filter

hanning = Norm*(0.5-0.5*cos(2*pi*sqrt(x^2+y^2+...)/(S1-1) ))

The parameters are (D, N, S1)

racos

Generic raised cosine filter. With an alpha=0.5 it is identical to the Hanning filter. With alpha=0.54 it is a Hamming filter.

racos = Norm*(alpha-(1-alpha)*cos(2*pi*sqrt(x^2+y^2+...)/(S1-1)))

The parameters are (D, N, A, S1)

The (optional) origin specification is used to locate the origin of the kernel. Typically 2D kernels have the origin in the "center" of the image. Many 1D kernels have the origin as the first element in the array.

The user can explicitly specify the origin. This alleviates common restrictions in some other tools (e.g. length of data axes must be an odd number).

If no origin is specified, the default is to assign the origin to the center of the array (rounded down). Thus a (5,4) array will have its origin set at (2,2).

Parameter=writekernel (boolean default=no)

Output kernel to an image?

Parameter=kernelfile (file filetype=output default=./. autoname=yes)

File name for kernel image. Autoname suffix is "_kernel".

Parameter=writefft (boolean default=no)

Output FFT of data and kernel?

Parameter=fftroot (file filetype=output default=./. autoname=yes)

Root file name for FFT output files. The autoname suffix for the real image is "_kernel_real". For the imaginary image, "_kernel_imag". The data file suffixes are "_data_real" and "_data_imag". ALWAYS use autonaming in this parameter.

Parameter=method (string default=slide)

Convolution method: (slide|fft).

The convolution of a data set with a kernel yields a measure of "how much" the data set looks like the kernel at every location in the data set. Another way to think about the convolution operation is a filtering of the dataset by the kernel.

Two convolution methods are implemented. Under some simple constraints they will provide the same results.

method=slide

Sliding cell convolution is a convolution from first principles.

O(a1,a2,...aN) = S_x1 S_x2... S_xN

I(x1,x2,..xN)*K(a1-x1,a2-x2,...aN-xN).

[replace S_x by summation signs]

method=fft

The FFT of the data and the kernel is computed. The arrays are multiplied and the inverse FFT is taken. The FFT route used is FFTW.

Internally the kernel and data MUST be the same size. The program will pad either/or both the data and/or kernel to the maximum length along either axis, defined by the larger of the two datasets.

If the edges of the data are "wrapped" (see below) then the FFT and slide methods will yield the same results.

For moderately-large to large kernels, the FFT convolution is much faster ( O(N log N) as opposed to O(N**2) ), however a) it requires much more memory and b) it restricts the edge treatment to "wrap".

Parameter=edges (string default=wrap)

How should the code handle edges? (wrap/nearest/mirror/constant/renorm).

As the kernel moves across the data, the edge of the kernel may fall off the data space. When this happens the program needs to know how the user wants to handle the edges. Five methods for handling the edges are implemented: For the 5 cases described below, consider the following example data space:

  data = { 1, 2, 3}

If you set method="fft" then you are implicitly using the edges="wrap" option.

Parameter=const (real default=0)

Constant value to use at edges.

Parameter=pad (boolean default=no)

Pad the data to the next integer power of 2.

The user can specify that the data be padded such that the length of each data axes is promoted to the next integer power of 2. The data are always padded on the "right" hand side, thus the array 1,2,3,4,5 would be padded to 1,2,3,4,5,0,0,0. This way of padding results in no change in the origin.

Parameter=center (boolean default=no)

Center FFT output?

If center = yes, zero frequency will appear at the center of the FFT images. This option is currently inoperative.

Parameter=normkernel (string default=area)

Normalize the kernel?

Allows aconvolve to automatically normalize the kernel. There are three options.

Parameter=clobber (boolean default=no)

Clobber existing output?

Parameter=verbose (integer default=0 min=0 max=5)

Verbosity level: 0 - no output, 5 - max display.

Changes in CIAO 4.11


Bugs

There are no known bugs for this tool.

See Also

tools
acrosscorr, apowerspectrum, arestore, arfcorr, csmooth, dmcoords, dmfilth, dmregrid, install_marx, make_psf_asymmetry_region, psf_project_ray, psfsize_srcs, simulate_psf, src_psffrac