.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/applications/plot_cornea_spot_inpainting.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. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_applications_plot_cornea_spot_inpainting.py: ============================================ Restore spotted cornea image with inpainting ============================================ Optical coherence tomography (OCT) is a non-invasive imaging technique used by ophthalmologists to take pictures of the back of a patient's eye [1]_. When performing OCT, dust may stick to the reference mirror of the equipment, causing dark spots to appear on the images. The problem is that these dirt spots cover areas of in-vivo tissue, hence hiding data of interest. Our goal here is to restore (reconstruct) the hidden areas based on the pixels near their boundaries. This tutorial is adapted from an application shared by Jules Scholler [2]_. The images were acquired by Viacheslav Mazlin (see :func:`skimage.data.palisades_of_vogt`). .. [1] David Turbert, reviewed by Ninel Z Gregori, MD (2023) `What Is Optical Coherence Tomography? `_, American Academy of Ophthalmology. .. [2] Jules Scholler (2019) "Image denoising using inpainting" https://www.jscholler.com/2019-02-28-remove-dots/ .. GENERATED FROM PYTHON SOURCE LINES 25-34 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import plotly.io import plotly.express as px import skimage as ski .. GENERATED FROM PYTHON SOURCE LINES 35-38 The dataset we are using here is an image sequence (a movie!) of human in-vivo tissue. Specifically, it shows the *palisades of Vogt* of a given cornea sample. .. GENERATED FROM PYTHON SOURCE LINES 40-42 Load image data =============== .. GENERATED FROM PYTHON SOURCE LINES 42-49 .. code-block:: Python image_seq = ski.data.palisades_of_vogt() print(f'number of dimensions: {image_seq.ndim}') print(f'shape: {image_seq.shape}') print(f'dtype: {image_seq.dtype}') .. rst-class:: sphx-glr-script-out .. code-block:: none number of dimensions: 3 shape: (60, 1440, 1440) dtype: uint16 .. GENERATED FROM PYTHON SOURCE LINES 50-57 The dataset is an image stack with 60 frames (time points) and 2 spatial dimensions. Let us visualize 10 frames by sampling every six time points: We can see some changes in illumination. We take advantage of the ``animation_frame`` parameter in Plotly's ``imshow`` function. As a side note, when the ``binary_string`` parameter is set to ``True``, the image is represented as grayscale. .. GENERATED FROM PYTHON SOURCE LINES 57-68 .. code-block:: Python fig = px.imshow( image_seq[::6, :, :], animation_frame=0, binary_string=True, labels={'animation_frame': '6-step time point'}, title='Sample of in-vivo human cornea', ) fig.update_layout(autosize=False, minreducedwidth=250, minreducedheight=250) plotly.io.show(fig) .. raw:: html :file: images/sphx_glr_plot_cornea_spot_inpainting_001.html .. GENERATED FROM PYTHON SOURCE LINES 69-80 Aggregate over time =================== First, we want to detect those dirt spots where the data are lost. In technical terms, we want to *segment* the dirt spots (for all frames in the sequence). Unlike the actual data (signal), the dirt spots do not move from one frame to the next; they are still. Therefore, we begin by computing a time aggregate of the image sequence. We shall use the median image to segment the dirt spots, the latter then standing out with respect to the background (blurred signal). Complementarily, to get a feel for the (moving) data, let us compute the variance. .. GENERATED FROM PYTHON SOURCE LINES 80-97 .. code-block:: Python image_med = np.median(image_seq, axis=0) image_var = np.var(image_seq, axis=0) assert image_var.shape == image_med.shape print(f'shape: {image_med.shape}') fig, ax = plt.subplots(ncols=2, figsize=(12, 6)) ax[0].imshow(image_med, cmap='gray') ax[0].set_title('Image median over time') ax[1].imshow(image_var, cmap='gray') ax[1].set_title('Image variance over time') fig.tight_layout() .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_002.png :alt: Image median over time, Image variance over time :srcset: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none shape: (1440, 1440) .. GENERATED FROM PYTHON SOURCE LINES 98-114 Use local thresholding ====================== To segment the dirt spots, we use thresholding. The images we are working with are unevenly illuminated, which causes spatial variations in the (absolute) intensities of the foreground and the background. For example, the average background intensity in one region may be different in another (distant) one. It is therefore more fitting to compute different threshold values across the image, one for each region. This is called adaptive (or local) thresholding, as opposed to the usual thresholding procedure which employs a single (global) threshold for all pixels in the image. When calling the ``threshold_local`` function from the ``filters`` module, we may change the default neighborhood size (``block_size``), i.e., the typical size (number of pixels) over which illumination varies, as well as the ``offset`` (shifting the neighborhood's weighted mean). Let us try two different values for ``block_size``: .. GENERATED FROM PYTHON SOURCE LINES 114-121 .. code-block:: Python thresh_1 = ski.filters.threshold_local(image_med, block_size=21, offset=15) thresh_2 = ski.filters.threshold_local(image_med, block_size=43, offset=15) mask_1 = image_med < thresh_1 mask_2 = image_med < thresh_2 .. GENERATED FROM PYTHON SOURCE LINES 122-124 Let us define a convenience function to display two plots side by side, so it is easier for us to compare them: .. GENERATED FROM PYTHON SOURCE LINES 124-136 .. code-block:: Python def plot_comparison(plot1, plot2, title1, title2): fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6), sharex=True, sharey=True) ax1.imshow(plot1, cmap='gray') ax1.set_title(title1) ax2.imshow(plot2, cmap='gray') ax2.set_title(title2) plot_comparison(mask_1, mask_2, "block_size = 21", "block_size = 43") .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_003.png :alt: block_size = 21, block_size = 43 :srcset: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 137-145 The "dirt spots" appear to be more distinct in the second mask, i.e., the one resulting from using the larger ``block_size`` value. We noticed that increasing the value of the offset parameter from its default zero value would yield a more uniform background, letting the objects of interest stand out more visibly. Note that toggling parameter values can give us a deeper understanding of the method being used, which can typically move us closer to the desired results. .. GENERATED FROM PYTHON SOURCE LINES 145-152 .. code-block:: Python thresh_0 = ski.filters.threshold_local(image_med, block_size=43) mask_0 = image_med < thresh_0 plot_comparison(mask_0, mask_2, "No offset", "Offset = 15") .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_004.png :alt: No offset, Offset = 15 :srcset: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 153-166 Remove fine-grained features ============================ We use morphological filters to sharpen the mask and focus on the dirt spots. The two fundamental morphological operators are *dilation* and *erosion*, where dilation (resp. erosion) sets the pixel to the brightest (resp. darkest) value of the neighborhood defined by a structuring element (footprint). Here, we use the ``diamond`` function from the ``morphology`` module to create a diamond-shaped footprint. An erosion followed by a dilation is called an *opening*. First, we apply an opening filter, in order to remove small objects and thin lines, while preserving the shape and size of larger objects. .. GENERATED FROM PYTHON SOURCE LINES 166-171 .. code-block:: Python footprint = ski.morphology.diamond(3) mask_open = ski.morphology.opening(mask_2, footprint) plot_comparison(mask_2, mask_open, "mask before", "after opening") .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_005.png :alt: mask before, after opening :srcset: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 172-182 Since "opening" an image starts with an erosion operation, bright regions which are smaller than the structuring element have been removed. When applying an opening filter, tweaking the footprint parameter lets us control how fine-grained the removed features are. For example, if we used ``footprint = ski.morphology.diamond(1)`` in the above, we could see that only smaller features would be filtered out, hence retaining more spots in the mask. Conversely, if we used a disk-shaped footprint of same radius, i.e., ``footprint = ski.morphology.disk(3)``, more of the fine-grained features would be filtered out, since the disk's area is larger than the diamond's. .. GENERATED FROM PYTHON SOURCE LINES 184-185 Next, we can make the detected areas wider by applying a dilation filter: .. GENERATED FROM PYTHON SOURCE LINES 185-189 .. code-block:: Python mask_dilate = ski.morphology.dilation(mask_open, footprint) plot_comparison(mask_open, mask_dilate, "Before", "After dilation") .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_006.png :alt: Before, After dilation :srcset: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 190-192 Dilation enlarges bright regions and shrinks dark regions. Notice how, indeed, the white spots have thickened. .. GENERATED FROM PYTHON SOURCE LINES 194-202 Inpaint each frame separately ============================= We are now ready to apply inpainting to each frame. For this we use function ``inpaint_biharmonic`` from the ``restoration`` module. It implements an algorithm based on biharmonic equations. This function takes two arrays as inputs: The image to restore and a mask (with same shape) corresponding to the regions we want to inpaint. .. GENERATED FROM PYTHON SOURCE LINES 202-210 .. code-block:: Python image_seq_inpainted = np.zeros(image_seq.shape) for i in range(image_seq.shape[0]): image_seq_inpainted[i] = ski.restoration.inpaint_biharmonic( image_seq[i], mask_dilate ) .. GENERATED FROM PYTHON SOURCE LINES 211-214 Let us visualize one restored image, where the dirt spots have been inpainted. First, we find the contours of the dirt spots (well, of the mask) so we can draw them on top of the restored image: .. GENERATED FROM PYTHON SOURCE LINES 214-237 .. code-block:: Python contours = ski.measure.find_contours(mask_dilate) # Gather all (row, column) coordinates of the contours x = [] y = [] for contour in contours: x.append(contour[:, 0]) y.append(contour[:, 1]) # Note that the following one-liner is equivalent to the above: # x, y = zip(*((contour[:, 0], contour[:, 1]) for contour in contours)) # Flatten the coordinates x_flat = np.concatenate(x).ravel().round().astype(int) y_flat = np.concatenate(y).ravel().round().astype(int) # Create mask of these contours contour_mask = np.zeros(mask_dilate.shape, dtype=bool) contour_mask[x_flat, y_flat] = 1 # Pick one frame sample_result = image_seq_inpainted[12] # Normalize it (so intensity values range [0, 1]) sample_result /= sample_result.max() .. GENERATED FROM PYTHON SOURCE LINES 238-241 We use function ``label2rgb`` from the ``color`` module to overlay the restored image with the segmented spots, using transparency (alpha parameter). .. GENERATED FROM PYTHON SOURCE LINES 241-253 .. code-block:: Python color_contours = ski.color.label2rgb( contour_mask, image=sample_result, alpha=0.4, bg_color=(1, 1, 1) ) fig, ax = plt.subplots(figsize=(6, 6)) ax.imshow(color_contours) ax.set_title('Segmented spots over restored image') fig.tight_layout() .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_007.png :alt: Segmented spots over restored image :srcset: /auto_examples/applications/images/sphx_glr_plot_cornea_spot_inpainting_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 254-257 Note that the dirt spot located at (x, y) ~ (719, 1237) stands out; ideally, it should have been segmented and inpainted. We can see that we 'lost' it to the opening processing step, when removing fine-grained features. .. GENERATED FROM PYTHON SOURCE LINES 257-259 .. code-block:: Python plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 16.467 seconds) .. _sphx_glr_download_auto_examples_applications_plot_cornea_spot_inpainting.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/scikit-image/scikit-image/v0.24.0?filepath=notebooks/auto_examples/applications/plot_cornea_spot_inpainting.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_cornea_spot_inpainting.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_cornea_spot_inpainting.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_cornea_spot_inpainting.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_