Filter Sequence

The class implements the functionality to specify a list (or just one) of filters/ regularization functionals / denoisers that can be used in any of the LEAP iterative reconstruction algorithms or can be used to create algorithm just for denoising. These filters can be either differentiable (e.g., TV) or non-differentiable (e.g. median filter), but those filters that are not differentiable will be ignored by gradient-based iterative reconstruction algorithms. Currently the only non-gradient-based iterative reconstruction algorithm in LEAP is ASDPOCS.

For example, if one wishes to employ a Total Variation (TV) regularizer, do the following:

filters = filterSequence(1.0e0)
filters.append(TV(leapct, delta=0.02/20.0, p=1.0))

where the argument to filterSequence (in this case 1.0e0) specifies the regularization strength. Larger values enforce a stronger regularization penalty and thus produce smoother results. This parameter should optimized by the user and optimal values vary on the input data and task. It is recommended that one first test various powers of ten, e.g., 1e-3, 1e-2, …, 1e3 and then narrow ones search once they find something that seems to work pretty well.

For TV denoising, it is important to set the delta parameter to an appropriate value. See more information about this below in the TV description.

class leap_filter_sequence.filterSequence(beta=1.0)

This class defines a weighted sum of filters (i.e., regularizers)

Parameters:: beta (float) – the overall strength of the sequence of filters, if zero no filters are applied

append(newFilter)

Append a new filter to the list

Note that in the case of ASDPOCS, the order of the sequence of filters matters because they are applied sequentially. When used in a gradient-based algorithm, the cost of filter sequence is given by self.beta * sum_n filters[n].cost(f) and note that each filter can have its own weight, given by filters[n].weight

Parameters:: newFilter (object whose base class is denoisingFilter) – the denoising filter to add to the sequence of filters

apply(f): Applies all the filters in sequence

clear(): Removes all filters from the list

cost(f): Calculates the accumulated cost of all filters

gradient(f): Calculates the accumulated gradient of all filters

quadForm(f, d): Calculates the accumulated quadratic form of all filters

class leap_filter_sequence.denoisingFilter(leapct)

Parent class for denoising filters

All filters (i.e., regularizers) must be written as a Python class which inherets this class. The apply function must be defined for all filters. This function performs denoising of the given input. If the filter is differentiable, you must define cost, gradient, and quadForm functions. Nondifferentiable filters (e.g., median filter, bilateral filter, etc.) may also be defined, but can only be used by ASDPOCS. On the other hand all differentiable filters can be used in any of the LEAP iterative reconstruction algorithms.

apply(f)

Filters (denoises) the input

If this function is not defined in the child class, then this performs one gradient descent iteration as a means to apply the given filter. This only works for differentiable filters.

Parameters:: f (C contiguous float32 numpy or torch array) – volume data, current estimate
Returns:: The filtered (denoised) input (C contiguous float32 numpy or torch array, same dimensions as the input)

cost(f)

Calculates the cost of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The cost (float) of the cost functional

gradient(f)

Calculates the gradient of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The gradient (C contiguous float32 numpy or torch array; has the same dimensions as the input) of the cost functional

quadForm(f, d)

Calculates the quadratic form of the given volume and descent direction

Parameters:

f (C contiguous float32 numpy or torch array) – volume data, current estimate
d (C contiguous float32 numpy or torch array) – volume data, descent direction

Returns:

The quadratic form (float) of the cost functional

class leap_filter_sequence.TV(leapct, delta=0.0, p=1.2, weight=1.0, f_0=None)

This class defines a filter based on leapct anisotropic Total Variation (TV) regularizer

This regularization functional uses a Huber-like loss function applied to the differences of neighboring samples (in 3D). One can switch between using 6 or 26 neighbors using the "set_numTVneighbors" function. The aTV functional with Huber-like loss function is given by

\[\begin{split}\begin{eqnarray} R(f) &:=& \sum_{\boldsymbol{i}} \sum_{\boldsymbol{j} \in N_{\boldsymbol{i}}} \|\boldsymbol{i} - \boldsymbol{j}\|^{-1} h(f_\boldsymbol{i} - f_\boldsymbol{j}) \\ h(t) &:=& \begin{cases} \frac{1}{2}t^2, & \text{if } |t| \leq delta \\ \frac{delta^{2 - p}}{p}|t|^p + delta^2\left(\frac{1}{2} - \frac{1}{p}\right), & \text{if } |t| > delta \end{cases} \end{eqnarray}\end{split}\]

where \(N_{\boldsymbol{i}}\) is a neighborhood around the 3D voxel index \(\boldsymbol{i} = (i_1, i_2, i_3)\). It is recommended that one set the value of delta to be smaller than the difference of two materials that one wishes to discriminate between. For example, if the reconstructed values of two materials are 0.02 and 0.04, than one should set delta to be smaller than 0.04-0.02 = 0.02. The guidance of “smaller than” is rather vague; I usually start with dividing this difference by 20, i.e., set delta = (0.04-0.02)/20.0. The delta parameter needs to be optimized for your specific problem, but the above strategy is a good place to start.

Parameters:

leapct (object of the tomographicModels class)
delta (float) – parameter for the Huber-like loss function used in TV
weight (float) – the regularizaion strength of this denoising filter term
f_0 (C contiguous float32 numpy or torch array) – a prior volume; this is optional but if specified this class calculates TV(f-f_0), e.g., PICCS

apply(f)

Filters (denoises) the input

If this function is not defined in the child class, then this performs one gradient descent iteration as a means to apply the given filter. This only works for differentiable filters.

Parameters:: f (C contiguous float32 numpy or torch array) – volume data, current estimate
Returns:: The filtered (denoised) input (C contiguous float32 numpy or torch array, same dimensions as the input)

cost(f)

Calculates the anisotropic Total Variation functional of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The cost (float) of the cost functional

gradient(f)

Calculates the gradient of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The gradient (C contiguous float32 numpy or torch array; has the same dimensions as the input) of the cost functional

quadForm(f, d)

Calculates the quadratic form of the given volume and descent direction

Parameters:

f (C contiguous float32 numpy or torch array) – volume data, current estimate
d (C contiguous float32 numpy or torch array) – volume data, descent direction

Returns:

The quadratic form (float) of the cost functional

class leap_filter_sequence.BlurFilter(leapct, FWHM)

This class defines a filter based on leapct.tomographicModels.BlurFilter which is a basic low pass filter

Parameters:

leapct (object of the tomographicModels class)
FWHM (float) – The Full Width at Half Maximum (given in number of pixels) of the low pass filter

apply(f)

Filters (denoises) the input

If this function is not defined in the child class, then this performs one gradient descent iteration as a means to apply the given filter. This only works for differentiable filters.

Parameters:: f (C contiguous float32 numpy or torch array) – volume data, current estimate
Returns:: The filtered (denoised) input (C contiguous float32 numpy or torch array, same dimensions as the input)

class leap_filter_sequence.BilateralFilter(leapct, spatialFWHM, intensityFWHM, scale=1.0)

This class defines a filter based on leapct.tomographicModels.BilateralFilter

Parameters:

leapct (object of the tomographicModels class)
spatialFWHM (float) – The Full Width at Half Maximum (given in number of pixels) of the filter
intensityFWHM (float) – The FWHM (given in the voxel value domain, which is usually mm^-1) of the filter
scale (float) – The FWHM of a low pass filter applied to the voxel values used in the intensity filter

apply(f)

Filters (denoises) the input

If this function is not defined in the child class, then this performs one gradient descent iteration as a means to apply the given filter. This only works for differentiable filters.

Parameters:: f (C contiguous float32 numpy or torch array) – volume data, current estimate
Returns:: The filtered (denoised) input (C contiguous float32 numpy or torch array, same dimensions as the input)

class leap_filter_sequence.GuidedFilter(leapct, r, epsilon)

This class defines a filter based on leapct.tomographicModels.GuidedFilter

Parameters:

leapct (object of the tomographicModels class)
r (int) – the window radius (in number of pixels)
epsilon (float) – the degree of smoothing

apply(f)

Filters (denoises) the input

If this function is not defined in the child class, then this performs one gradient descent iteration as a means to apply the given filter. This only works for differentiable filters.

Parameters:: f (C contiguous float32 numpy or torch array) – volume data, current estimate
Returns:: The filtered (denoised) input (C contiguous float32 numpy or torch array, same dimensions as the input)

class leap_filter_sequence.MedianFilter(leapct, threshold=0.0, windowSize=3)

This class defines a filter based on leapct.tomographicModels.MedianFilter

Parameters:

leapct (object of the tomographicModels class)
threshold (float) – the threshold of whether to use the filtered value or not
windowSize (int) – the window size; can be 3 or 5

apply(f)

Filters (denoises) the input

If this function is not defined in the child class, then this performs one gradient descent iteration as a means to apply the given filter. This only works for differentiable filters.

Parameters:: f (C contiguous float32 numpy or torch array) – volume data, current estimate
Returns:: The filtered (denoised) input (C contiguous float32 numpy or torch array, same dimensions as the input)

class leap_filter_sequence.LpNorm(leapct, delta=0.0, p=1.0, weight=1.0, f_0=None, FWHM=0.0)

This class defines a filter based on the L_p norm (raised to the p power) of the input

This cost functional of this regularizer is given by

\[\begin{split}\begin{eqnarray} C(f) &:=& \sum h(B(f - f_0)) \\ h(t) &:=& \begin{cases} \frac{1}{2}t^2, & \text{if } |t| \leq delta \\ \frac{delta^{2 - p}}{p}|t|^p + delta^2\left(\frac{1}{2} - \frac{1}{p}\right), & \text{if } |t| > delta \end{cases} \end{eqnarray}\end{split}\]

The B operator is either the identity transform (does nothing), a low-pass, or a high-pass filter. See the FWHM parameter description.

Parameters:

leapct (object of the tomographicModels class)
delta (float) – parameter for the Huber-like loss function used in TV
p (float) – The p-value of the L_p norm
weight (float) – the regularizaion strength of this denoising filter term
f_0 (C contiguous float32 numpy or torch array) – a prior volume; this parameter is optional
FWHM (float) – the units are in number of pixels. if FWHM > 1.0, B is a low-pass filter with the specified FWHM and if FWHM < -1.0, B is a high-pass filter with the specified FWHM, otherwise B is the identity transform

cost(f)

Calculates the cost of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The cost (float) of the cost functional

gradient(f)

Calculates the gradient of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The gradient (C contiguous float32 numpy or torch array; has the same dimensions as the input) of the cost functional

quadForm(f, d)

Calculates the quadratic form of the given volume and descent direction

Parameters:

f (C contiguous float32 numpy or torch array) – volume data, current estimate
d (C contiguous float32 numpy or torch array) – volume data, descent direction

Returns:

The quadratic form (float) of the cost functional

class leap_filter_sequence.histogramSparsity(leapct, mus=None, weight=1.0)

This class defines a filter that encourages sparisty in the histogram domain

Warning: this is a nonconvex regularizer. It is best to use this filter after an initial reconstruction is performed.

Parameters:

leapct (object of the tomographicModels class)
mus (numpy array or list of floats) – list of target values expected in the reconstruction volume
weight (float) – the regularizaion strength of this denoising filter term

cost(f)

Calculates the cost of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The cost (float) of the cost functional

gradient(f)

Calculates the gradient of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The gradient (C contiguous float32 numpy or torch array; has the same dimensions as the input) of the cost functional

quadForm(f, d)

Calculates the quadratic form of the given volume and descent direction

Parameters:

f (C contiguous float32 numpy or torch array) – volume data, current estimate
d (C contiguous float32 numpy or torch array) – volume data, descent direction

Returns:

The quadratic form (float) of the cost functional

class leap_filter_sequence.azimuthalFilter(leapct, FWHM, p, weight=1.0)

This class defines a filter based on leapct AzimuthalBlur filter

This denoising filter applies a high pass filter in the azimuthal direction of the reconstructed volume. If is useful for reconstructing objects that have a sparse azimuthal gradient, such as pipes with delaminations, voids, or high density inclusions.

The functional is given by the L_p norm (raised to the p power) of the azimuthal high pass filter of the input, i.e., ||H(f)||_p^p, where H is the azimuthal high pass filter

Parameters:

leapct (object of the tomographicModels class)
FWHM (float) – The FWHM (measured in degrees) of the high pass filter
p (float) – the p-value of the L_p norm
weight (float) – the regularizaion strength of this denoising filter term

cost(f)

Calculates the cost of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The cost (float) of the cost functional

gradient(f)

Calculates the gradient of the given volume

Parameters:: f (C contiguous float32 numpy or torch array) – volume data
Returns:: The gradient (C contiguous float32 numpy or torch array; has the same dimensions as the input) of the cost functional

quadForm(f, d)

Calculates the quadratic form of the given volume and descent direction

Parameters:

f (C contiguous float32 numpy or torch array) – volume data, current estimate
d (C contiguous float32 numpy or torch array) – volume data, descent direction

Returns:

The quadratic form (float) of the cost functional