.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/07-tours/03-bias-field.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_07-tours_03-bias-field.py: ================== Receive bias field ================== When you combine the images of a receive array by root sum of squares, what you get is the object weighted by the root sum of squares of the coil sensitivities. That weight is smooth, multiplicative and largest near the elements, so it shades the image: the same tissue looks brighter near the array than far from it. The shading biases segmentation, intensity-based registration and any quantitative comparison across the field of view. In this tour you shade a T1-weighted head with a head array that has no anterior elements, estimate the field with N4 [#tustison]_ through :func:`bartorch.tools.bias_field_correct`, and compare the corrected image and the estimated field with the object and the true field. Prerequisites: lesson 1, :doc:`../01-basics/01-tensors-and-commands`. **Learning objectives** - Relate the shading of a root-sum-of-squares image to the receive sensitivities. - Estimate and remove the bias field with N4, and judge the result by the uniformity of a tissue class and by the intensity histogram. - Choose the mask and the fitting grid of the estimate. - Recognise what N4 does not recover: the scale of the field, the steep part of a field of large range, and a smooth intensity variation that belongs to the object. .. GENERATED FROM PYTHON SOURCE LINES 33-84 .. code-block:: Python import torch import bartorch import bartorch.tools as bt .. GENERATED FROM PYTHON SOURCE LINES 85-92 The object ---------- Your object is an axial slice of the BrainWeb T1-weighted head [#brainweb]_ through the lateral ventricles, 128 x 128 over a 220 mm field of view. The BrainWeb tissue model also tells you which voxels are at least 90 % white or grey matter, and you measure the uniformity of each class over those. .. GENERATED FROM PYTHON SOURCE LINES 92-98 .. code-block:: Python SIZE, FOV_MM = 128, 220.0 .. GENERATED FROM PYTHON SOURCE LINES 99-108 The shaded image ---------------- The array is BART's analytical eight-element head coil with its two anterior elements removed, as in an open-face head coil: six elements. Each coil image is the object weighted by one element's sensitivity :math:`S_c`, plus independent complex Gaussian noise. Their root sum of squares is the object weighted by the bias field :math:`B = \sqrt{\sum_c |S_c|^2}`, and the noise adds a Rician floor in the background. .. GENERATED FROM PYTHON SOURCE LINES 108-127 .. code-block:: Python def shaded_image(elements): grid = bt.grid(D=(SIZE, SIZE, 1)) sensitivities = bt.coils(t=grid, n=8)[elements, 0] bias = bartorch.rss(sensitivities, axes=(0,)).abs() sensitivities = sensitivities / bias[brain].mean() coil_images = bt.noise(sensitivities * image, n=1e-4, s=3) shaded = bartorch.rss(coil_images, axes=(0,)).abs() return shaded, bias / bias[brain].mean() shaded, bias = shaded_image([0, 1, 2, 5, 6, 7]) .. rst-class:: sphx-glr-script-out .. code-block:: none bias field over the brain: 0.56 to 1.27 (2nd to 98th percentile) .. GENERATED FROM PYTHON SOURCE LINES 128-142 N4 correction ------------- N4 works on the logarithm of the image, where the field is no longer a factor but an additive term: the log of the object plus a smooth field, represented by cubic B-splines. It estimates the field by alternately sharpening the histogram of the log intensities and fitting the splines to what the sharpening removed, over a hierarchy of control-point grids. It knows nothing about the coil; all it needs is an object whose intensities form classes. The field is fitted on the image shrunk by ``shrink_factor``. The default of 4 would leave a 32 x 32 grid of this 128 matrix, so you pass 2, which keeps 64 x 64. Ask for the field as well as the corrected image: .. GENERATED FROM PYTHON SOURCE LINES 142-148 .. code-block:: Python corrected, estimate = bt.bias_field_correct( shaded, shrink_factor=2, return_field=True, ) .. GENERATED FROM PYTHON SOURCE LINES 149-153 You measure the uniformity of a tissue class by its coefficient of variation, the standard deviation over the mean. The object has some too, from partial volume at the class boundaries, so that is the floor to compare with. .. GENERATED FROM PYTHON SOURCE LINES 153-169 .. rst-class:: sphx-glr-script-out .. code-block:: none object coefficient of variation: white matter 0.011, grey matter 0.035; white/grey 1.35 shaded coefficient of variation: white matter 0.135, grey matter 0.193; white/grey 1.38 N4-corrected coefficient of variation: white matter 0.035, grey matter 0.073; white/grey 1.35 .. GENERATED FROM PYTHON SOURCE LINES 170-177 The correction takes the white-matter variation from more than ten times the object's to about three times it, and keeps the white/grey contrast of the object. A multiplicative field is only determined up to a constant factor, and the correction leaves that factor in the image. So you compare the estimate with the true field after scaling both to unit mean over the brain: .. GENERATED FROM PYTHON SOURCE LINES 177-263 .. code-block:: Python estimate = estimate / estimate[brain].mean() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_03-bias-field_001.png :alt: object, shaded, N4-corrected :srcset: /auto_examples/07-tours/images/sphx_glr_03-bias-field_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_03-bias-field_002.png :alt: true field, N4 estimate :srcset: /auto_examples/07-tours/images/sphx_glr_03-bias-field_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_03-bias-field_003.png :alt: image / object on the line, field on the line :srcset: /auto_examples/07-tours/images/sphx_glr_03-bias-field_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_03-bias-field_004.png :alt: intensity histogram of the brain :srcset: /auto_examples/07-tours/images/sphx_glr_03-bias-field_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none estimated / true field over the brain: median 0.997, within 10 % in 93 % of the voxels .. GENERATED FROM PYTHON SOURCE LINES 264-280 The images share one window, each scaled to its white-matter mean. In the shaded image you see the frontal lobes darker than the occipital lobes. Along the dashed line the ratio of the shaded image to the object follows the field, and the histogram of the brain has no separate grey- and white-matter peaks. After correction the ratio stays close to one except at the frontal pole, where the estimated field does not fall as far as the true one, and the two peaks separate again. Over the brain the estimated field is within 10 % of the true one in more than nine voxels out of ten. The mask -------- N4 fits the field only over a mask, by default Otsu's threshold of the image. You can give it a mask of the whole head instead, which adds the scalp and the skull, or a brain mask such as a skull-stripping tool gives you; here the BrainWeb tissue model provides it. .. GENERATED FROM PYTHON SOURCE LINES 280-290 .. code-block:: Python masks = {"Otsu": None, "head": torch.as_tensor(head), "brain": brain} for name, mask in masks.items(): trial = bt.bias_field_correct(shaded, mask=mask, shrink_factor=2) spread = variation(trial, white) print(f"{name:5s} mask: white-matter variation {spread:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Otsu mask: white-matter variation 0.035 head mask: white-matter variation 0.063 brain mask: white-matter variation 0.023 .. GENERATED FROM PYTHON SOURCE LINES 291-301 N4 sharpens one histogram over the whole mask. The scalp and the skull bring intensity classes of their own, so the fit over the head mask is the least uniform. The brain mask, which holds only the classes you judge the field by, gives the most uniform white matter. Residual bias ------------- The estimate gets worse as the range of the field grows. Keep only the four posterior elements, and the field falls steeply towards the frontal pole: .. GENERATED FROM PYTHON SOURCE LINES 301-328 .. code-block:: Python steep, steep_bias = shaded_image([0, 1, 6, 7]) steep_corrected, steep_estimate = bt.bias_field_correct( steep, shrink_factor=2, return_field=True, ) .. image-sg:: /auto_examples/07-tours/images/sphx_glr_03-bias-field_005.png :alt: shaded, N4-corrected, field on the line :srcset: /auto_examples/07-tours/images/sphx_glr_03-bias-field_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none four posterior elements: field 0.35 to 1.55 white-matter variation: shaded 0.262, N4-corrected 0.096; estimate within 10 % in 83 % of the voxels .. GENERATED FROM PYTHON SOURCE LINES 329-350 N4 overestimates the field at the frontal pole: its estimate does not fall as far as the true one, so the corrected frontal cortex stays darker than the occipital cortex. A residual like this is why you still compare a corrected image across regions with care. There is a second limit, and it is not a failure of the fit. N4 removes any smooth intensity variation, whatever its origin: a receive field, a transmit field in a gradient-echo image, or a genuine slow change of the tissue signal. In a SENSE reconstruction the sensitivities are normalized to unit root sum of squares, which removes the receive shading inside the reconstruction. N4 is for images whose sensitivities you do not have, such as a root-sum-of-squares combination, and for the transmit field, which no receive calibration measures. Exercise: B-spline fitting grid ------------------------------- The field is fitted on the shrunk image and evaluated on the full grid. A field that is smooth on the scale of the head needs few grid points, and each N4 iteration gets cheaper as their number falls. Vary ``shrink_factor`` and watch the white-matter uniformity: .. GENERATED FROM PYTHON SOURCE LINES 350-355 .. code-block:: Python for shrink in (1, 2, 4, 8): trial = bt.bias_field_correct(shaded, shrink_factor=shrink) spread = variation(trial, white) print(f"shrink_factor {shrink}: white-matter variation {spread:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none shrink_factor 1: white-matter variation 0.029 shrink_factor 2: white-matter variation 0.035 shrink_factor 4: white-matter variation 0.036 shrink_factor 8: white-matter variation 0.056 .. GENERATED FROM PYTHON SOURCE LINES 356-373 Factors of 1 to 4 do about equally well here. At 8 the shrunk image is 16 x 16 and holds too few voxels of each tissue class, and the variation rises. As a spec --------- What this tour built, stated the way you would ask an agent for it: .. code-block:: text Correct the receive bias field of a 2D root-sum-of-squares magnitude image with bt.bias_field_correct (N4), shrink_factor 2, returning the field. Compare a mask of None (Otsu), the whole head and the brain. Report the coefficient of variation over white and grey matter before and after, the white/grey contrast, and the fraction of brain voxels where the estimated field, scaled to unit brain mean, is within 10 % of the true field. .. GENERATED FROM PYTHON SOURCE LINES 375-386 References ---------- .. [#tustison] Tustison NJ, Avants BB, Cook PA, Zheng Y, Egan A, Yushkevich PA, Gee JC. N4ITK: improved N3 bias correction. *IEEE Trans Med Imaging* 29(6):1310-1320 (2010). https://doi.org/10.1109/TMI.2010.2046908 .. [#brainweb] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, Evans AC. Design and construction of a realistic digital brain phantom. *IEEE Trans Med Imaging* 17(3):463-468 (1998). https://doi.org/10.1109/42.712135 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.720 seconds) .. _sphx_glr_download_auto_examples_07-tours_03-bias-field.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03-bias-field.ipynb <03-bias-field.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03-bias-field.py <03-bias-field.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03-bias-field.zip <03-bias-field.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_