pirt.gaussfun - Gaussian diffusion

The gaussfun module implements functions for diffusion and Gaussian derivatives for data of any dimension.

Contents of this module:

  • gaussiankernel - Create a Gaussian kernel
  • gaussiankernel2 - Create a 2D Gaussian kernel
  • diffusionkernel - Create a discrete analog to the Gaussian kernel
  • gfilter - Filter data with Gaussian (derivative) kernels
  • diffuse - Filter data with true discrete diffusion kernels
  • gfilter2 - Filter data by specifying sigma in world coordinates
  • diffuse2 - Diffuse data by specifying sigma in world coordinates

pirt.gaussiankernel(sigma, order=0, N=None, returnt=False, warn=True)

Creates a 1D gaussian derivative kernel with the given sigma and the given order. (An order of 0 is a “regular” Gaussian.)

The returned kernel is a column vector, thus working in the first dimension (in images, this often is y).

The returned kernel is odd by default. Using N one can specify the full kernel size (if not int, the ceil operator is applied). By specifying a negative value for N, the tail length (number of elements on both sides of the center element) can be specified. The total kernel size than becomes ceil(-N)*2+1. Though the method to supply it may be a bit obscure, this measure can be handy, since the tail length if often related to the sigma. If not given, the optimal N is determined automatically, depending on sigma and order.

If the given scale is a small for the given order, a warning is produced (unless warn==True).

—– Used Literature:

Koenderink, J. J. The structure of images. Biological Cybernetics 50, 5 (1984), 363-370.

Lindeberg, T. Scale-space for discrete signals. IEEE Transactions on Pattern Analysis and Machine Intelligence 12, 3 (1990), 234-254.

Ter Haar Romeny, B. M., Niessen, W. J., Wilting, J., and Florack, L. M. J. Differential structure of images: Accuracy of representation. In First IEEE International Conference on Image Processing, (Austin, TX) (1994).

pirt.diffusionkernel(sigma, N=4, returnt=False)

A discrete analog to the continuous Gaussian kernel, as proposed by Toni Lindeberg.

N is the tail length factor (relative to sigma).

pirt.gfilter(L, sigma, order=0, mode='constant', warn=True)

Gaussian filterering and Gaussian derivative filters.

Parameters:
L : np.ndarray

The input data to filter

sigma : scalar or list-of-scalars

The smoothing parameter, can be given for each dimension

order : int or list-of-ints

The order of the derivative, can be given for each dimension

mode : {‘reflect’,’constant’,’nearest’,’mirror’, ‘wrap’}

Determines how edge effects are handled. (see scipy.ndimage.convolve1d)

warn : boolean

Whether to show a warning message if the sigma is too small to represent the required derivative.

Notes

Makes use of the seperability property of the Gaussian by convolving 1D kernels in each dimension.

pirt.gfilter2(L, scale, order=0, mode='reflect', warn=True)

Apply Gaussian filtering by specifying a scale in world coordinates rather than a sigma. This function determines the sigmas to apply, based on the sampling of the elements.

See gfilter for more information.

(If L is not an Aarray, this function yields the same result as gfilter.)

pirt.diffuse(L, sigma)

Diffusion using a discrete variant of the diffusion operator.

Parameters:
L : np.ndarray

The input data to filter

sigma : scalar or list-of-scalars

The smoothing parameter, can be given for each dimension

pirt.diffuse2(L, scale, mode='nearest')

Apply diffusion by specifying a scale in world coordinates rather than a sigma. This function determines the sigmas to apply, based on the sampling of the elements.

See diffuse for more information.

(If L is not an Aarray, this function yields the same result as diffuse.)


Technically implemented in pyramid.py, but very much related to Gaussian diffusion:

class pirt.ScaleSpacePyramid(data, min_scale=None, scale_offset=0, use_buffer=False, level_factor=2)

The scale space pyramid class provides a way to manage a scale space pyramid. Given an input image (of arbitrary dimension), it provides two simple methods to obtain the image at the a specified scale or level.

Parameters:
data : numpy array

An array of any dimension. Should preferably be of float type.

min_scale : scalar, optional

The minimum scale to sample from the pyramid. If not given, scale_offset is used. If larger than zero, the image is smoothed to this scale before creating the zeroth level. If the smoothness is sufficient, the data is also downsampled. This makes a registration algorithm much faster, because the image data for the final scales does not have a unnecessary high resolution.

scale_offset : scalar

The scale of the given data. Use this if the data is already smooth. Be careful not to set this value too high, as aliasing artifacts may be introduced. Default zero.

use_buffer : bool

Whether a result obtained with get_scale() is buffered for later use. Only one image is buffered. Default False.

level_factor : scalar

The scale distance between two levels. A larger number means saving a bit of memory in trade of speed. You’re probably fine with 2.0.

Notes

Note that this scale space representation handles anisotropic arrays and that scale is expressed in world units.

Note that images at higher levels do not always have a factor 2 sampling difference with the original! This is because the first and last pixel are kept the same, and the number of pixels is decreased with factors of two (or almost a factor of two if the number is uneven).

The images always have the same offset though.

We adopt the following concepts:
  • level: the level in the pyramid. Each level is a factor two smaller in size (in each dimension) than the previous.
  • scale: the scale in world coordinates
calculate(levels=None, min_shape=None)

Create the image pyramid now. Specify either the amount of levels, or the minimum shape component of the highest level. If neither levels nor min_shape is given, uses min_shape=8.

Returns (max_level, max_sigma) of the current pyramid.

get_level(self, level)

get_level(level):

Get the image at the specified (integer) level, zero being the lowest level (the original image).

Each level is approximately a factor two smaller in size that the previous level. All levels are buffered.

The returned image has two added properties: _pyramid_scale and _pyramid_level, wich specify the image scale and level in the pyramid.

get_scale(scale)

Get the image at the specified scale (expressed in world units). For higher scales, the image has a smaller shape than the original image. If min_scale and scale_offset are not used, a scale of 0 represents the original image.

To calculate the result, the image at the level corresponding to the nearest lower scale is obtained, and diffused an extra bit to obtain the requested scale.

The result is buffered (if the pyramid was instantiated with use_buffer=True), such that calling this function multiple times with the same scale is much faster. Only buffers the last used scale.

The returned image has two added properties: _pyramid_scale and _pyramid_level, wich specify the image scale and level in the pyramid.