Especially for MRI, the pixel values can vary between scanner types. This will lead to a very differnt range of pixel intensities across images in the same batch and the algorithms will have a hard time converging. Various techniques for rescaling exist, of which some are implemented here.

@patch
def show_L(l:L, **kwargs): 
    "shows all Tensors in L"
    show_images_3d(torch.stack(tuple(l)), **kwargs)
    
@patch
def apply(l:L, fun):
    "applies a function to each element of L"
    return L(fun(item) for item in l)

fns = Path('../data/series').ls()
images = L(TensorDicom3D.create(fn) for fn in fns)
images.show_L(nrow = 15)

Size correction

If working with data from different scanners and even different sequences, pixel spacing can differ. E.g. T2w images are usually in a higher resolution as DWI images, thus the size of each individual pixel is larger in DWI images than in T2W. Normalizing pixel size to a standartized value can help achiving better performance. However, only the H x W dimension is rescaled, not the D dimension.

from faimed3d.basics import TensorDicom3D, TensorMask3D # for compatibility with show_docs

TensorDicom3D.size_correction[source]

TensorDicom3D.size_correction(im:TensorMask3D'>), new_spacing=1)

TensorMask3D.size_correction[source]

TensorMask3D.size_correction(im:TensorMask3D'>), new_spacing=1)

[im.get_spacing() for im in images]

[(0.6835939884185791, 0.6835939884185791, 2.0),
 (0.6835939884185791, 0.6835939884185791, 6.0),
 (0.6835939884185791, 0.6835939884185791, 6.0),
 (0.6835939884185791, 0.6835939884185791, 6.0),
 (0.6835939884185791, 0.6835939884185791, 1.0)]

Pixel rescaling

rescale intercept: 0028|1052
rescale slope: 0028|1053

from faimed3d.basics import TensorDicom3D # for compatibility with show_docs

TensorDicom3D.rescale_pixeldata[source]

TensorDicom3D.rescale_pixeldata(t:TensorDicom3D)

im = images[0]
print(im.mean())
# the example MRI images have no rescale slope or intercept
im.set_metadata('0028|1052', 5.)
im.set_metadata('0028|1053', 0.5)
im = im.rescale_pixeldata()
im.mean()

TensorDicom3D(173.5362)

TensorDicom3D(91.7681)

Pixel intensity normalization

Mean / Max / Median scaling

TensorDicom3D.mean_scale[source]

TensorDicom3D.mean_scale(t:TensorDicom3D)

Scales pixels by subtracting the mean and dividing by std. 0 pixels are ignored

class MeanScale[source]

MeanScale(p=1.0, **kwargs) :: RandTransform

A transform that before_call its state at each `__call__`

images.apply(MeanScale()).show_L(nrow = 15)

TensorDicom3D.median_scale[source]

TensorDicom3D.median_scale(t:TensorDicom3D)

Scales pixels by subtracting the median and dividing by the IQR. 0 pixels are ignored

class MedianScale[source]

MedianScale(p=1.0, **kwargs) :: RandTransform

A transform that before_call its state at each `__call__`

TensorDicom3D.max_scale[source]

TensorDicom3D.max_scale(t:TensorDicom3D)

class MaxScale[source]

MaxScale(p=1.0, **kwargs) :: RandTransform

A transform that before_call its state at each `__call__`

Histogram scaling

However, just scaling with one mean (or median) and std (or IQR) might not be the optimal solution. Histogram scaling can be used for a more even distribution of pixel values. Different functions are provided for histogram scaling:

Method adapted from Jeremy Howard

See the excellent Kaggle kernel from Jeremy Howard for an explanation.

from torch import Tensor # for compatibility with show_docs

TensorDicom3D.freqhist_bins[source]

TensorDicom3D.freqhist_bins(t:Tensor'>), n_bins=100)

A function to split the range of pixel values into groups, so that each group has around the same number of pixels.
taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

Tensor.freqhist_bins[source]

Tensor.freqhist_bins(t:Tensor'>), n_bins=100)

A function to split the range of pixel values into groups, so that each group has around the same number of pixels.
taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

TensorDicom3D.hist_scaled[source]

TensorDicom3D.hist_scaled(t:Tensor'>), brks=None)

Scales a tensor using `freqhist_bins` to values between 0 and 1
taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

Tensor.hist_scaled[source]

Tensor.hist_scaled(t:Tensor'>), brks=None)

Scales a tensor using `freqhist_bins` to values between 0 and 1
taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

TensorDicom3D.hist_scaled_pt[source]

TensorDicom3D.hist_scaled_pt(t:Tensor'>), brks=None)

same as fastai fucntion for PILDicom

Tensor.hist_scaled_pt[source]

Tensor.hist_scaled_pt(t:Tensor'>), brks=None)

same as fastai fucntion for PILDicom

The sample data from radiopaedia is already windowed and scaled, so the results presented here are not very impressive. However, if applied to MRI data, changes will be stronger.

def f(x): return x.hist_scaled()
images.apply(f).show_L(nrow = 15)

Comparison to SimpleITK adaptive histogram equalization

The method from Jeremy Howard is more simple that methods provided in SimpleITK, but faster.

im1 = TensorDicom3D(images[0][10:20])

image = sitk.Cast(im1.as_sitk(), sitk.sitkFloat32)
output = sitk.AdaptiveHistogramEqualization(image, radius=[20]*3, alpha=0.3, beta=0.3) # ~ 30-60 seconds
TensorDicom3D(sitk.GetArrayFromImage(output)).show() 
im1.hist_scaled().show()

N4 bias field correction

from official SimpleITK docs

The N4 bias field correction algorithm is a popular method for correcting low frequency intensity non-uniformity present in MRI image data known as a bias or gain field. The method has also been successfully applied as flat-field correction in microscopy data. This method assumes a simple parametric model and does not require tissue classification.

corrector = sitk.N4BiasFieldCorrectionImageFilter()
corrector.SetMaximumNumberOfIterations([3]*3)
output = corrector.Execute(sitk.Cast(im1.as_sitk(), sitk.sitkFloat32)) # needs float
TensorDicom3D(sitk.GetArrayFromImage(output)).show() #

faimed3d provides a wrapper class to apply both N4 bias field correction and adaptive histogram equalization, which basically just simplifies reading the data.

class ImageCorrectionWrapper[source]

ImageCorrectionWrapper(n4_max_num_it=3, hist_radius=[5, 5, 5], hist_alpha=0.3, hist_beta=0.3, do_n4=True, do_hist=True, verbose=True)

Piecewise linear histogram matching

[1] N. Laszlo G and J. K. Udupa, “On Standardizing the MR Image Intensity Scale,” Magn. Reson. Med., vol. 42, pp. 1072–1081, 1999.

[2] M. Shah, Y. Xiao, N. Subbanna, S. Francis, D. L. Arnold, D. L. Collins, and T. Arbel, “Evaluating intensity normalization on MRIs of human brain with multiple sclerosis,” Med. Image Anal., vol. 15, no. 2, pp. 267–282, 2011.

Implementation adapted from: https://github.com/jcreinhold/intensity-normalization, ported to pytorch (no use of numpy works in cuda).

In contrast to hist_scaled, the piecewise linear histogram matching need pre-specified values for new scale and landmarks. It should be used to normalize a whole dataset.

get_percentile[source]

get_percentile(t:Tensor, q:float)

Return the ``q``-th percentile of the flattened input tensor's data.

CAUTION:
 * Needs PyTorch >= 1.1.0, as ``torch.kthvalue()`` is used.
 * Values are not interpolated, which corresponds to
   ``numpy.percentile(..., interpolation="nearest")``.

:param t: Input tensor.
:param q: Percentile to compute, which must be between 0 and 100 inclusive.
:return: Resulting value (float).

This function is twice as fast as torch.quantile and has no size limitations

get_landmarks[source]

get_landmarks(t:Tensor, percentiles:Tensor)

Returns the input's landmarks.

:param t (torch.Tensor): Input tensor.
:param percentiles (torch.Tensor): Peraentiles to calculate landmarks for.
:return: Resulting landmarks (torch.tensor).

find_standard_scale[source]

find_standard_scale(inputs, i_min=1, i_max=99, i_s_min=1, i_s_max=100, l_percentile=10, u_percentile=90, step=10)

determine the standard scale for the set of images
Args:
    inputs (list or L): set of TensorDicom3D objects which are to be normalized
    i_min (float): minimum percentile to consider in the images
    i_max (float): maximum percentile to consider in the images
    i_s_min (float): minimum percentile on the standard scale
    i_s_max (float): maximum percentile on the standard scale
    l_percentile (int): middle percentile lower bound (e.g., for deciles 10)
    u_percentile (int): middle percentile upper bound (e.g., for deciles 90)
    step (int): step for middle percentiles (e.g., for deciles 10)
Returns:
    standard_scale (np.ndarray): average landmark intensity for images
    percs (np.ndarray): array of all percentiles used

Tensor.piecewise_hist[source]

Tensor.piecewise_hist(image:Tensor, landmark_percs, standard_scale)

Do the Nyul and Udupa histogram normalization routine with a given set of learned landmarks

Args:
    input_image (TensorDicom3D): image on which to find landmarks
    landmark_percs (torch.tensor): corresponding landmark points of standard scale
    standard_scale (torch.tensor): landmarks on the standard scale
Returns:
    normalized (TensorDicom3D): normalized image

class PiecewiseHistScaling[source]

PiecewiseHistScaling(landmark_percs=None, standard_scale=None, final_scale=None, p=1.0, slicewise=False, slice_dim=None, **kwargs) :: RandTransform

Applies theNyul and Udupa histogram nomalization and rescales the pixel values.

Args:
    input_image (TensorDicom3D): image on which to find landmarks
    landmark_percs (torch.tensor): corresponding landmark points of standard scale
    final_scale (function): final rescaling of values, if none is provided values are
                            scaled to a mean of 0 and a std of 1.
    slicewise (bool): if the scaling should be applied to each slice individually. Slower but leads to more homogeneous images.

Returns:
    If input is TensorMask3D returns input unchanged
    If input is TensorDicom3D returns normalized and rescaled Tensor

standard_scale, percs = find_standard_scale(images)
standard_scale, percs

(tensor([313.1477, 343.0235, 348.7227, 353.8698, 359.4490, 364.1021, 367.5030,
         370.4590, 373.5393, 377.9435, 385.5532]),
 tensor([ 1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 99]))

images.apply(PiecewiseHistScaling(landmark_percs=percs, 
                                  standard_scale=standard_scale)).show_L(nrow = 15)

standard_scale_from_filelist[source]

standard_scale_from_filelist(fns:Series'>))

calculates standard scale from images in a filelist

from fastai.data.core import DataLoaders # for compatibility with show_docs

DataLoaders.standard_scale_from_dls[source]

DataLoaders.standard_scale_from_dls(dls:DataLoaders)

calculates standard scale from images in a dataloader