.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/applications/plot_fluorescence_nuclear_envelope.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_applications_plot_fluorescence_nuclear_envelope.py>`
        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_fluorescence_nuclear_envelope.py:


======================================================
Measure fluorescence intensity at the nuclear envelope
======================================================

This example reproduces a well-established workflow in bioimage data analysis
for measuring the fluorescence intensity localized to the nuclear envelope, in
a time sequence of cell images (each with two channels and two spatial
dimensions) which shows a process of protein re-localization from the
cytoplasmic area to the nuclear envelope. This biological application was
first presented by Andrea Boni and collaborators in [1]_; it was used in a
textbook by Kota Miura [2]_ as well as in other works ([3]_, [4]_).
In other words, we port this workflow from ImageJ Macro to Python with
scikit-image.

.. [1] Boni A, Politi AZ, Strnad P, Xiang W, Hossain MJ, Ellenberg J (2015)
       "Live imaging and modeling of inner nuclear membrane targeting reveals
       its molecular requirements in mammalian cells" J Cell Biol
       209(5):705–720. ISSN: 0021-9525.
       :DOI:`10.1083/jcb.201409133`
.. [2] Miura K (2020) "Measurements of Intensity Dynamics at the Periphery of
       the Nucleus" in: Miura K, Sladoje N (eds) Bioimage Data Analysis
       Workflows. Learning Materials in Biosciences. Springer, Cham.
       :DOI:`10.1007/978-3-030-22386-1_2`
.. [3] Klemm A (2020) "ImageJ/Fiji Macro Language" NEUBIAS Academy Online
       Course: https://www.youtube.com/watch?v=o8tfkdcd3DA
.. [4] Vorkel D and Haase R (2020) "GPU-accelerating ImageJ Macro image
       processing workflows using CLIJ" https://arxiv.org/abs/2008.11799

.. GENERATED FROM PYTHON SOURCE LINES 31-42

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    import plotly.io
    import plotly.express as px
    from scipy import ndimage as ndi

    from skimage import filters, measure, morphology, segmentation
    from skimage.data import protein_transport









.. GENERATED FROM PYTHON SOURCE LINES 43-44

We start with a single cell/nucleus to construct the workflow.

.. GENERATED FROM PYTHON SOURCE LINES 44-49

.. code-block:: Python


    image_sequence = protein_transport()

    print(f'shape: {image_sequence.shape}')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    shape: (15, 2, 180, 183)




.. GENERATED FROM PYTHON SOURCE LINES 50-51

The dataset is a 2D image stack with 15 frames (time points) and 2 channels.

.. GENERATED FROM PYTHON SOURCE LINES 51-65

.. code-block:: Python


    vmin, vmax = 0, image_sequence.max()

    fig = px.imshow(
        image_sequence,
        facet_col=1,
        animation_frame=0,
        zmin=vmin,
        zmax=vmax,
        binary_string=True,
        labels={'animation_frame': 'time point', 'facet_col': 'channel'},
    )
    plotly.io.show(fig)




.. raw:: html
    :file: images/sphx_glr_plot_fluorescence_nuclear_envelope_001.html





.. GENERATED FROM PYTHON SOURCE LINES 66-68

To begin with, let us consider the first channel of the first image (step
``a)`` in the figure below).

.. GENERATED FROM PYTHON SOURCE LINES 68-71

.. code-block:: Python


    image_t_0_channel_0 = image_sequence[0, 0, :, :]








.. GENERATED FROM PYTHON SOURCE LINES 72-79

Segment the nucleus rim
=======================
Let us apply a Gaussian low-pass filter to this image in order to smooth it
(step ``b)``).
Next, we segment the nuclei, finding the threshold between the background
and foreground with Otsu's method: We get a binary image (step ``c)``). We
then fill the holes in the objects (step ``c-1)``).

.. GENERATED FROM PYTHON SOURCE LINES 79-87

.. code-block:: Python


    smooth = filters.gaussian(image_t_0_channel_0, sigma=1.5)

    thresh_value = filters.threshold_otsu(smooth)
    thresh = smooth > thresh_value

    fill = ndi.binary_fill_holes(thresh)








.. GENERATED FROM PYTHON SOURCE LINES 88-91

Following the original workflow, let us remove objects which touch the image
border (step ``c-2)``). Here, we can see that part of another nucleus was
touching the bottom right-hand corner.

.. GENERATED FROM PYTHON SOURCE LINES 91-95

.. code-block:: Python


    clear = segmentation.clear_border(fill)
    clear.dtype





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    dtype('bool')



.. GENERATED FROM PYTHON SOURCE LINES 96-98

We compute both the morphological dilation of this binary image
(step ``d)``) and its morphological erosion (step ``e)``).

.. GENERATED FROM PYTHON SOURCE LINES 98-103

.. code-block:: Python


    dilate = morphology.binary_dilation(clear)

    erode = morphology.binary_erosion(clear)








.. GENERATED FROM PYTHON SOURCE LINES 104-107

Finally, we subtract the eroded from the dilated to get the nucleus rim
(step ``f)``). This is equivalent to selecting the pixels which are in
``dilate``, but not in ``erode``:

.. GENERATED FROM PYTHON SOURCE LINES 107-110

.. code-block:: Python


    mask = np.logical_and(dilate, ~erode)








.. GENERATED FROM PYTHON SOURCE LINES 111-112

Let us visualize these processing steps in a sequence of subplots.

.. GENERATED FROM PYTHON SOURCE LINES 112-144

.. code-block:: Python


    fig, ax = plt.subplots(2, 4, figsize=(12, 6), sharey=True)

    ax[0, 0].imshow(image_t_0_channel_0, cmap=plt.cm.gray)
    ax[0, 0].set_title('a) Raw')

    ax[0, 1].imshow(smooth, cmap=plt.cm.gray)
    ax[0, 1].set_title('b) Blur')

    ax[0, 2].imshow(thresh, cmap=plt.cm.gray)
    ax[0, 2].set_title('c) Threshold')

    ax[0, 3].imshow(fill, cmap=plt.cm.gray)
    ax[0, 3].set_title('c-1) Fill in')

    ax[1, 0].imshow(clear, cmap=plt.cm.gray)
    ax[1, 0].set_title('c-2) Keep one nucleus')

    ax[1, 1].imshow(dilate, cmap=plt.cm.gray)
    ax[1, 1].set_title('d) Dilate')

    ax[1, 2].imshow(erode, cmap=plt.cm.gray)
    ax[1, 2].set_title('e) Erode')

    ax[1, 3].imshow(mask, cmap=plt.cm.gray)
    ax[1, 3].set_title('f) Nucleus Rim')

    for a in ax.ravel():
        a.set_axis_off()

    fig.tight_layout()




.. image-sg:: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_002.png
   :alt: a) Raw, b) Blur, c) Threshold, c-1) Fill in, c-2) Keep one nucleus, d) Dilate, e) Erode, f) Nucleus Rim
   :srcset: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 145-149

Apply the segmented rim as a mask
=================================
Now that we have segmented the nuclear membrane in the first channel, we use
it as a mask to measure the intensity in the second channel.

.. GENERATED FROM PYTHON SOURCE LINES 149-165

.. code-block:: Python


    image_t_0_channel_1 = image_sequence[0, 1, :, :]
    selection = np.where(mask, image_t_0_channel_1, 0)

    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 6), sharey=True)

    ax0.imshow(image_t_0_channel_1)
    ax0.set_title('Second channel (raw)')
    ax0.set_axis_off()

    ax1.imshow(selection)
    ax1.set_title('Selection')
    ax1.set_axis_off()

    fig.tight_layout()




.. image-sg:: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_003.png
   :alt: Second channel (raw), Selection
   :srcset: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 166-170

Measure the total intensity
===========================
The mean intensity is readily available as a region property in a labeled
image.

.. GENERATED FROM PYTHON SOURCE LINES 170-177

.. code-block:: Python


    props = measure.regionprops_table(
        mask.astype(np.uint8),
        intensity_image=image_t_0_channel_1,
        properties=('label', 'area', 'intensity_mean'),
    )








.. GENERATED FROM PYTHON SOURCE LINES 178-179

We may want to check that the value for the total intensity

.. GENERATED FROM PYTHON SOURCE LINES 179-182

.. code-block:: Python


    selection.sum()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    np.uint64(28350)



.. GENERATED FROM PYTHON SOURCE LINES 183-184

can be recovered from:

.. GENERATED FROM PYTHON SOURCE LINES 184-187

.. code-block:: Python


    props['area'] * props['intensity_mean']





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    array([28350.])



.. GENERATED FROM PYTHON SOURCE LINES 188-193

Process the entire image sequence
=================================
Instead of iterating the workflow for each time point, we process the
multidimensional dataset directly (except for the thresholding step).
Indeed, most scikit-image functions support nD images.

.. GENERATED FROM PYTHON SOURCE LINES 193-200

.. code-block:: Python


    n_z = image_sequence.shape[0]  # number of frames

    smooth_seq = filters.gaussian(image_sequence[:, 0, :, :], sigma=(0, 1.5, 1.5))
    thresh_values = [filters.threshold_otsu(s) for s in smooth_seq[:]]
    thresh_seq = [smooth_seq[k, ...] > val for k, val in enumerate(thresh_values)]








.. GENERATED FROM PYTHON SOURCE LINES 201-215

Alternatively, we could compute ``thresh_values`` without using a list
comprehension, by reshaping ``smooth_seq`` to make it 2D (where the first
dimension still corresponds to time points, but the second and last
dimension now contains all pixel values), and applying the thresholding
function on the image sequence along its second axis:

.. code-block:: python

    thresh_values = np.apply_along_axis(filters.threshold_otsu,
                                        axis=1,
                                        arr=smooth_seq.reshape(n_z, -1))

We use the following flat structuring element for morphological
computations (``np.newaxis`` is used to prepend an axis of size 1 for time):

.. GENERATED FROM PYTHON SOURCE LINES 215-219

.. code-block:: Python


    footprint = ndi.generate_binary_structure(2, 1)[np.newaxis, ...]
    footprint





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    array([[[False,  True, False],
            [ True,  True,  True],
            [False,  True, False]]])



.. GENERATED FROM PYTHON SOURCE LINES 220-222

This way, each frame is processed independently (pixels from consecutive
frames are never spatial neighbors).

.. GENERATED FROM PYTHON SOURCE LINES 222-225

.. code-block:: Python


    fill_seq = ndi.binary_fill_holes(thresh_seq, structure=footprint)








.. GENERATED FROM PYTHON SOURCE LINES 226-250

When clearing objects which touch the image border, we want to make sure
that the bottom (first) and top (last) frames are not considered as borders.
In this case, the only relevant border is the edge at the greatest (x, y)
values. This can be seen in 3D by running the following code:

.. code-block:: python

    import plotly.graph_objects as go

    sample = fill_seq
    (n_Z, n_Y, n_X) = sample.shape
    Z, Y, X = np.mgrid[:n_Z, :n_Y, :n_X]

    fig = go.Figure(
        data=go.Volume(
            x=X.flatten(),
            y=Y.flatten(),
            z=Z.flatten(),
            value=sample.flatten(),
            opacity=0.5,
            slices_z=dict(show=True, locations=[n_z // 2])
        )
    )
    fig.show()

.. GENERATED FROM PYTHON SOURCE LINES 250-259

.. code-block:: Python


    border_mask = np.ones_like(fill_seq)
    border_mask[n_z // 2, -1, -1] = False
    clear_seq = segmentation.clear_border(fill_seq, mask=border_mask)

    dilate_seq = morphology.binary_dilation(clear_seq, footprint=footprint)
    erode_seq = morphology.binary_erosion(clear_seq, footprint=footprint)
    mask_sequence = np.logical_and(dilate_seq, ~erode_seq)








.. GENERATED FROM PYTHON SOURCE LINES 260-263

Let us give each mask (corresponding to each time point) a different label,
running from 1 to 15. We use ``np.min_scalar_type`` to determine the
minimum-size integer dtype needed to represent the number of time points:

.. GENERATED FROM PYTHON SOURCE LINES 263-269

.. code-block:: Python


    label_dtype = np.min_scalar_type(n_z)
    mask_sequence = mask_sequence.astype(label_dtype)
    labels = np.arange(1, n_z + 1, dtype=label_dtype)
    mask_sequence *= labels[:, np.newaxis, np.newaxis]








.. GENERATED FROM PYTHON SOURCE LINES 270-272

Let us compute the region properties of interest for all these labeled
regions.

.. GENERATED FROM PYTHON SOURCE LINES 272-296

.. code-block:: Python


    props = measure.regionprops_table(
        mask_sequence,
        intensity_image=image_sequence[:, 1, :, :],
        properties=('label', 'area', 'intensity_mean'),
    )
    np.testing.assert_array_equal(props['label'], np.arange(n_z) + 1)

    fluorescence_change = [
        props['area'][i] * props['intensity_mean'][i] for i in range(n_z)
    ]

    fluorescence_change /= fluorescence_change[0]  # normalization

    fig, ax = plt.subplots()
    ax.plot(fluorescence_change, 'rs')
    ax.grid()
    ax.set_xlabel('Time point')
    ax.set_ylabel('Normalized total intensity')
    ax.set_title('Change in fluorescence intensity at the nuclear envelope')
    fig.tight_layout()

    plt.show()




.. image-sg:: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_004.png
   :alt: Change in fluorescence intensity at the nuclear envelope
   :srcset: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 297-300

Reassuringly, we find the expected result: The total fluorescence
intensity at the nuclear envelope increases 1.3-fold in the initial five
time points, and then becomes roughly constant.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.726 seconds)


.. _sphx_glr_download_auto_examples_applications_plot_fluorescence_nuclear_envelope.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_fluorescence_nuclear_envelope.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_fluorescence_nuclear_envelope.ipynb <plot_fluorescence_nuclear_envelope.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_fluorescence_nuclear_envelope.py <plot_fluorescence_nuclear_envelope.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_