Image Filtering

A filter makes a response map from an image: it brings out edges, blobs, textures, directions or frequencies. Radiomic features of the response map then describe these patterns. Pictologics has the filters of IBSI 2, and it checks them against the IBSI 2 reference maps. The IBSI 2 reference manual gives the definition and the formula of each filter.

Overview

A filtered configuration usually has these steps:

1. Resample the image to cubic voxels, for example 1 × 1 × 1 mm. IBSI filters a resampled image.
2. Resegment the ROI by intensity (for example a HU range), when you need it. Put this step before the filter: after a filter, a resegment step reads the response values, not the HU values.
3. Filter the image.
4. Discretise the response map, for texture, histogram and IVH features.
5. Extract the features of the response map.

Every filter returns a float32 response map, half the memory of float64. Two rules set the type of the filter passes:

• Mean, Gaussian, Laws and wavelets: the passes run in float32 for every image. A float32 pass moves half the bytes, and the passes still sum each line in double (numba kernels with the arithmetic of scipy). The mean filter and the Laws energy sum each window on its own: the exact mean of the float32 values. So a value does not depend on where the image starts. For a float64 image, the response changes by the float32 rounding of the passes: about 1e-7 of the largest response.
• LoG, Simoncelli, Riesz and Gabor: a float64 image stays float64 in the passes and the FFTs, and the last step writes float32. A float64 image below 131,072 voxels (about 51³) runs its LoG passes in float32, as they are faster there. With a padding_value, the padded image counts. The normalized convolution of a source mask also keeps a float64 image in float64.

Filter Reference

• Direct calls: the functions are in pictologics.filters. A direct call takes spacing_mm=1.0 by default, so sigma_mm then counts voxels. Give the spacing of your image. The pipeline gives it for Gaussian, LoG, Gabor and Riesz-LoG.
• IBSI 2 Phase 1: 28 of 28 compared tests pass. 5 tests of 2D filters are not run: Pictologics filters 3D images (the Gabor filter works on 2D planes, as IBSI defines it). 3 tests have no reference map: 9.b.2 and 10.b.2 need structure-tensor steering, which Pictologics does not have, and the map of 10.a is missing. See the IBSI 2 Phase 1 report.

The Filter Figures

Each filter section below has a figure: the middle slice of a 3D test image of 64 × 64 × 64 voxels of 1 mm (left), and the response of the filter (right). Move the slider, or press Play, to see the response at other settings.

The test image has:

• a large sphere (a blob) and three small spheres (small blobs);
• a cube (edges and corners);
• two striped patches with a period of 6 mm: the left patch changes along x, the right patch along y;
• noise.

The colours of the responses:

• Grey (mean, Gaussian, Laws energy, Gabor): dark is low, and bright is high.
• Blue and orange (LoG, wavelets, Simoncelli, Riesz): blue is below 0, dark grey is 0, and orange is above 0.

Each response has its own scale. The spatial filters use the mirror boundary, and the Simoncelli and Riesz filters the periodic boundary (the defaults of a pipeline step). The script dev/docs_figures.py makes the figures with the Pictologics filters.

Boundaries and Padding

A filter reads voxels past the image edge. The boundary sets their values:

• Defaults: a pipeline filter step uses "mirror", except the Simoncelli and Riesz filters, which use "periodic". A direct call uses "zero" (spatial filters) or "periodic" (FFT filters). Give boundary in every step to make it clear.

• Constant value padding (IBSI 2 Z3VE): with "constant", padding_value sets the constant, for example -1000 to pad a CT image with air:

  {"step": "filter", "params": {
      "type": "log",
      "sigma_mm": 2.0,
      "boundary": "constant",
      "padding_value": -1000.0,
  }}
  

  Every filter takes padding_value. A value other than 0 needs the "constant" (or "zero") boundary; with another boundary, the filter raises ValueError. A spatial filter gives the response of the image padded with the constant as far as the filter reads (no padded copy is made). A Laws energy step then pads the response with 0.

FFT Filter Boundaries

The Simoncelli and Riesz filters work in the Fourier domain, which is periodic:

• "periodic" (default): the FFT works on the image as it is. This is the fastest path, and the one that IBSI specifies for the Simoncelli tests.
• "zero", "nearest", "mirror": the filter pads the image with the requested mode, filters, and cuts the result back. A composite filter (riesz_log, riesz_simoncelli) pads once around the whole chain. riesz_log then does the LoG and the Riesz transform in one FFT round trip on the padded image.

Filters in a Pipeline

• Checks: add_config checks every filter step. An unknown parameter gets a hint ("kernels" gives "did you mean 'kernel'?"). A missing sigma_mm, lambda_mm or order is an error, and so is a bad wavelet level or decomposition, a bad Riesz order, an unknown Gabor response, or a padding_value without the constant boundary.
• FBS after a filter: the response has other units than the image, so a filter step cancels the FBS start of an earlier resegment step. An FBS discretise step after a filter needs min_val, or a resegment step after the filter (IBSI asks for a fixed start).
• The region around the ROI: often no later step reads the image away from the ROI. Then the mean, Gaussian, LoG, Laws, wavelet and Gabor filters compute only the region around the ROI box. The region grows by the reach of the filter. The pipeline then keeps only that region: the image and the masks are cut to it, as after a discretise step. The features are the same, bit for bit, and a small ROI in a large image costs much less time. A later filter, resample, binarize_mask that keeps the value 0, or whole-image normalise step needs the whole grid. See Tips for Speed and Memory.
• The log: the log entry of each filter step records boundary_requested, boundary_effective, params_requested and params_effective (with the spacing and the boundary that the pipeline gave).

Mean Filter

The mean filter replaces each voxel with the mean of the cube of support voxels around it. It smooths the image and reduces noise.

{"step": "filter", "params": {"type": "mean", "support": 3, "boundary": "mirror"}}

from pictologics.filters import mean_filter

response = mean_filter(image_array, support=3, boundary="mirror")

Gaussian Filter

The Gaussian filter smooths the image with a Gaussian kernel of scale σ (in mm). It keeps the coarse structure and removes fine detail and noise.

{"step": "filter", "params": {"type": "gaussian", "sigma_mm": 2.0, "boundary": "mirror"}}

from pictologics.filters import gaussian_filter

response = gaussian_filter(image_array, sigma_mm=2.0, spacing_mm=(0.8, 0.8, 2.0), boundary="mirror")

Laplacian of Gaussian (LoG)

The LoG filter smooths the image with a Gaussian of scale σ, and then takes the Laplacian (the second derivative). It brings out edges and blobs of a size near σ: a larger σ finds larger blobs.

{"step": "filter", "params": {"type": "log", "sigma_mm": 3.0, "boundary": "mirror"}}

from pictologics.filters import laplacian_of_gaussian

response = laplacian_of_gaussian(image_array, sigma_mm=3.0, spacing_mm=(1.0, 1.0, 1.0), boundary="mirror")

Laws Texture Energy

Laws filters combine three 1D kernels (L: level, E: edge, S: spot, W: wave, R: ripple) into one 3D kernel, for example "L5E5S5". Each kernel brings out one type of texture.

{"step": "filter", "params": {
    "type": "laws",
    "kernel": "E5L5S5",
    "rotation_invariant": True,
    "pooling": "max",
    "compute_energy": True,
}}

from pictologics.filters import laws_filter

response = laws_filter(image_array, kernels="E5L5S5", rotation_invariant=True, compute_energy=True)

Gabor Filter

A Gabor filter is a complex sine wave in a Gaussian envelope. It brings out texture of one wavelength and one direction. Pictologics applies 2D Gabor kernels in the axial plane, or in all three orthogonal planes, as IBSI defines them.

{"step": "filter", "params": {
    "type": "gabor",
    "sigma_mm": 5.0,
    "lambda_mm": 2.0,
    "rotation_invariant": True,
    "delta_theta": 0.7853981633974483,  # pi/4: 8 directions
}}

from pictologics.filters import gabor_filter

response = gabor_filter(
    image_array,
    sigma_mm=5.0,
    lambda_mm=2.0,
    gamma=1.0,
    spacing_mm=(1.0, 1.0, 1.0),
    rotation_invariant=True,
    delta_theta=0.7853981633974483,
)

Separable Wavelets

A separable wavelet transform splits each axis into a low (L) and a high (H) frequency part. A 3D image so gives 8 sub-bands (LLL, LLH, ..., HHH). Pictologics uses the undecimated (stationary) transform, so the response map has the size of the image.

{"step": "filter", "params": {
    "type": "wavelet",
    "wavelet": "coif1",
    "level": 1,
    "decomposition": "LLH",
    "rotation_invariant": True,
}}

from pictologics.filters import wavelet_transform

response = wavelet_transform(image_array, wavelet="coif1", decomposition="LLH", level=1, rotation_invariant=True)

Simoncelli Wavelet

The Simoncelli wavelet is a non-separable, isotropic wavelet: it has no preferred direction. Use it when the direction of a texture does not matter.

{"step": "filter", "params": {"type": "simoncelli", "level": 1}}

from pictologics.filters import simoncelli_wavelet

response = simoncelli_wavelet(image_array, level=1)

Riesz Transform

The Riesz transform is a higher-order derivative in all directions at once. Its order (l1, l2, l3) gives the derivative order along each axis. Pictologics combines it with a LoG or a Simoncelli filter. It does not steer the filters (no structure-tensor alignment).

With the periodic boundary, the LoG stage of riesz_log uses the zero boundary. With another boundary, the whole chain uses the requested boundary.

{"step": "filter", "params": {"type": "riesz", "variant": "log", "order": [2, 0, 0], "sigma_mm": 3.0}}

from pictologics.filters import riesz_log, riesz_simoncelli, riesz_transform

response = riesz_transform(image_array, order=(2, 0, 0))
response = riesz_log(image_array, sigma_mm=3.0, spacing_mm=(1.0, 1.0, 1.0), order=(1, 0, 0))
response = riesz_simoncelli(image_array, level=1, order=(1, 0, 0))

Source Masks for Sentinel Values

Some images hold a sentinel value where they have no data, for example -2048 HU outside the field of view. These values would spoil a filter response. Every filter takes a source_mask (True: a voxel with data):

In a Pipeline

The pipeline makes the source mask from the source_mode of the configuration, and gives it to the filter. In "auto" mode with sentinel_value, the valid voxels are those without the sentinel value; in "roi_only" mode, the voxels of the ROI.

pipeline.add_config(
    "filtered_features",
    [
        {"step": "resample", "params": {"new_spacing": (1.0, 1.0, 1.0)}},
        {"step": "resegment", "params": {"range_min": -1000, "range_max": 3000}},  # HU, before the filter
        {"step": "filter", "params": {"type": "log", "sigma_mm": 3.0, "boundary": "mirror"}},
        {"step": "discretise", "params": {"method": "FBN", "n_bins": 32}},
        {"step": "extract_features", "params": {"families": ["intensity", "texture"]}},
    ],
    source_mode="auto",
    sentinel_value=-2048,
)

A resample step also removes the voxels without data from the masks. Without a resample step, a resegment step before the filter keeps the sentinel voxels out of the intensity mask (the voxels of the intensity and texture features).

In a Direct Call

from pictologics.filters import gabor_filter, laplacian_of_gaussian, mean_filter

valid = image_array > -2000  # True: a voxel with data

# Mean, Gaussian, LoG and Laws return (response, valid mask) with a source mask
response, response_valid = mean_filter(image_array, support=15, source_mask=valid)
response, response_valid = laplacian_of_gaussian(image_array, sigma_mm=3.0, source_mask=valid)

# The other filters return the response only
response = gabor_filter(image_array, sigma_mm=5.0, lambda_mm=2.0, source_mask=valid)

See also: Source Modes and Sentinel Values, the filter API and the IBSI 2 reports.