Examples

Simple phase

This example illustrates the simple usage of the qpimage.QPImage class for reading and managing quantitative phase data. The attribute QPImage.pha yields the background-corrected phase data and the attribute QPImage.bg_pha yields the background phase image.

_images/simple_phase.jpg

simple_phase.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import matplotlib.pylab as plt
import numpy as np
import qpimage

size = 200
# background phase image with a tilt
bg = np.repeat(np.linspace(0, 1, size), size).reshape(size, size)
# phase image with random noise
phase = np.random.rand(size, size) + bg

# create QPImage instance
qpi = qpimage.QPImage(data=phase, bg_data=bg, which_data="phase")

# plot the properties of `qpi`
plt.figure(figsize=(8, 3))
plot_kw = {"vmin": -1,
           "vmax": 2}

plt.subplot(131, title="fake input phase")
plt.imshow(phase, **plot_kw)

plt.subplot(132, title="fake background phase")
plt.imshow(qpi.bg_pha, **plot_kw)

plt.subplot(133, title="corrected phase")
plt.imshow(qpi.pha, **plot_kw)

plt.tight_layout()
plt.show()

Background image tilt correction

This example illustrates background tilt correction with qpimage. In contrast to the ‘simple_phase.py’ example, the known background data is not given to the qpimage.QPImage class. In this particular example, the background tilt correction achieves an error of about 1% which is sufficient in most quantitative phase imaging applications.

_images/background_tilt.jpg

background_tilt.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
import matplotlib.pylab as plt
import numpy as np
import qpimage

size = 200
# background phase image with a tilt
bg = np.repeat(np.linspace(0, 1, size), size).reshape(size, size)
bg = .6 * bg - .8 * bg.transpose() + .2
# phase image with random noise
rsobj = np.random.RandomState(47)
phase = rsobj.rand(size, size) - .5 + bg

# create QPImage instance
qpi = qpimage.QPImage(data=phase, which_data="phase")
# compute background with 2d tilt approach
qpi.compute_bg(which_data="phase",  # correct phase image
               fit_offset="fit",  # use bg offset from tilt fit
               fit_profile="tilt",  # perform 2D tilt fit
               border_px=5,  # use 5 px border around image
               )

# plot the properties of `qpi`
fig = plt.figure(figsize=(8, 2.5))
plot_kw = {"vmin": -1,
           "vmax": 1}

ax1 = plt.subplot(131, title="input data")
map1 = ax1.imshow(phase, **plot_kw)
plt.colorbar(map1, ax=ax1, fraction=.046, pad=0.04)

ax2 = plt.subplot(132, title="tilt-corrected")
map2 = ax2.imshow(qpi.pha, **plot_kw)
plt.colorbar(map2, ax=ax2, fraction=.046, pad=0.04)

ax3 = plt.subplot(133, title="tilt error")
map3 = ax3.imshow(bg - qpi.bg_pha)
plt.colorbar(map3, ax=ax3, fraction=.046, pad=0.04)

# disable axes
[ax.axis("off") for ax in [ax1, ax2, ax3]]

plt.tight_layout(pad=0, h_pad=0, w_pad=0)
plt.show()

Background image offset correction

This example illustrates the different background offset correction methods implemented in qpimage. The phase image data contains two gaussian noise distributions for which these methods yield different background phase offsets.

_images/background_offset.jpg

background_offset.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import matplotlib.pylab as plt
import numpy as np
import qpimage

size = 200  # the size of the image
bg = 2.5  # the center of the background phase distribution
scale = .1  # the spread of the background phase distribution

# compute random phase data
rsobj = np.random.RandomState(42)
data = rsobj.normal(loc=bg, scale=scale, size=size**2)
# Add a second distribution `data2` at random positions `idx`,
# such that there is no pure gaussian distribution.
# (otherwise 'mean' and 'gaussian' cannot be distinguished)
data2 = rsobj.normal(loc=bg*1.1, scale=scale, size=size**2//2)
idx = rsobj.choice(data.size, data.size//2)
data[idx] = data2
# reshape `data` to get a 2D array
data = data.reshape(size, size)

qpi = qpimage.QPImage(data=data, which_data="phase")

cpkw = {"which_data": "phase",  # correct the input phase data
        "fit_profile": "offset",  # perform offset correction only
        "border_px": 5,  # use a border of 5px of the input phase
        "ret_mask": True,  # return the mask image for visualization
        }

mask = qpi.compute_bg(fit_offset="mode", **cpkw)
bg_mode = np.mean(qpi.bg_pha[mask])

qpi.compute_bg(fit_offset="mean", **cpkw)
bg_mean = np.mean(qpi.bg_pha[mask])

qpi.compute_bg(fit_offset="gauss", **cpkw)
bg_gauss = np.mean(qpi.bg_pha[mask])

bg_data = (qpi.pha + qpi.bg_pha)[mask]
# compute histogram
nbins = int(np.ceil(np.sqrt(bg_data.size)))
mind, maxd = bg_data.min(), bg_data.max()
histo = np.histogram(bg_data, nbins, density=True, range=(mind, maxd))
dx = abs(histo[1][1] - histo[1][2]) / 2
hx = histo[1][1:] - dx
hy = histo[0]

# plot the properties of `qpi`
plt.figure(figsize=(8, 4))

ax1 = plt.subplot(121, title="input phase")
map1 = plt.imshow(data)
plt.colorbar(map1, ax=ax1, fraction=.046, pad=0.04)


t2 = "{}px border histogram with {} bins".format(cpkw["border_px"], nbins)
plt.subplot(122, title=t2)
plt.plot(hx, hy, label="histogram", color="gray")
plt.axvline(bg_mode, 0, 1, label="mode", color="red")
plt.axvline(bg_mean, 0, 1, label="mean", color="green")
plt.axvline(bg_gauss, 0, 1, label="gauss", color="orange")
plt.legend()

plt.tight_layout()
plt.show()

Masked background image correction

This example illustrates background correction with qpimage using a mask to exclude regions that do not contain background information.

The phase image of a microgel bead (top left) has two artifacts; there is a tilt-like phase profile added along the vertical axis and there is a second microgel bead in close proximity to the center bead. A regular phase tilt background correction using the image values around a frame of five pixels (see “background_tilt.py” example) does not yield a flat background, because the second bead is fitted into the background which leads to a horizontal background phase profile (top right). By defining a mask (bottom left image), the phase values of the second bead can be excluded from the background tilt fit and a flat background phase is achieved (bottom right).

_images/background_mask.jpg

background_mask.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
import matplotlib.pylab as plt
import numpy as np
import qpimage


# load the experimental data
input_phase = np.load("./data/phase_beads_close.npz")["phase"].astype(float)

# create QPImage instance
qpi = qpimage.QPImage(data=input_phase,
                      which_data="phase")

# background correction without mask
qpi.compute_bg(which_data="phase",
               fit_offset="fit",
               fit_profile="tilt",
               border_px=5,
               )
pha_nomask = qpi.pha

# educated guess for mask
mask = input_phase < input_phase.max() / 10

# background correction with mask
# (the intersection of `mask` and the 5px border is used for fitting)
qpi.compute_bg(which_data="phase",
               fit_offset="fit",
               fit_profile="tilt",
               border_px=5,
               from_mask=mask
               )
pha_mask = qpi.pha

# plot
fig = plt.figure(figsize=(8, 7))
plot_kw = {"vmin": -.1,
           "vmax": 1.5}

ax1 = plt.subplot(221, title="input phase")
map1 = ax1.imshow(input_phase, **plot_kw)
plt.colorbar(map1, ax=ax1, fraction=.044, pad=0.04)

ax2 = plt.subplot(222, title="tilt-corrected (no mask)")
map2 = ax2.imshow(pha_nomask, **plot_kw)
plt.colorbar(map2, ax=ax2, fraction=.044, pad=0.04)

ax3 = plt.subplot(223, title="mask")
map3 = ax3.imshow(1.*mask, cmap="gray_r")
plt.colorbar(map3, ax=ax3, fraction=.044, pad=0.04)

ax4 = plt.subplot(224, title="tilt-corrected (with mask)")
map4 = ax4.imshow(pha_mask, **plot_kw)
plt.colorbar(map4, ax=ax4, fraction=.044, pad=0.04)

# disable axes
[ax.axis("off") for ax in [ax1, ax2, ax3, ax3, ax4]]

plt.tight_layout(h_pad=0, w_pad=0)
plt.show()

Background image 2nd order polynomial correction

This example extends the tilt correction (‘background_tilt.py’) to a second order polynomial correction for samples that exhibit more sophisticated phase aberrations. The phase background correction is computed from a ten pixel wide frame around the image. The phase data shown are computed from a hologram of a single myeloid leukemia cell (HL60) recorded using digital holographic microscopy (DHM) (see [SSM+15]).

_images/background_poly2o.jpg

background_poly2o.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
import matplotlib.pylab as plt
import numpy as np
# The data are stored in a .jpg file (lossy compression).
# If `PIL` is not found, try installing the `pillow` package.
from PIL import Image
import qpimage

edata = np.array(Image.open("./data/hologram_cell_curved_bg.jpg"))

# create QPImage instance
qpi = qpimage.QPImage(data=edata, which_data="hologram")
pha0 = qpi.pha

# background correction using tilt only
qpi.compute_bg(which_data=["phase"],
               fit_offset="fit",
               fit_profile="tilt",
               border_px=10,
               )
pha_tilt = qpi.pha

# background correction using polynomial
qpi.compute_bg(which_data=["phase"],
               fit_offset="fit",
               fit_profile="poly2o",
               border_px=10,
               )
pha_poly2o = qpi.pha

# plot phase data
fig = plt.figure(figsize=(8, 3.3))

phakw = {"cmap": "viridis",
         "interpolation": "bicubic",
         "vmin": pha_poly2o.min(),
         "vmax": pha_poly2o.max()}

ax1 = plt.subplot(131, title="input phase")
map1 = ax1.imshow(pha0, **phakw)
plt.colorbar(map1, ax=ax1, fraction=.067, pad=0.04)

ax2 = plt.subplot(132, title="tilt correction only")
map2 = ax2.imshow(pha_tilt, **phakw)
plt.colorbar(map2, ax=ax2, fraction=.067, pad=0.04)

ax3 = plt.subplot(133, title="polynomial correction")
map3 = ax3.imshow(pha_poly2o, **phakw)
plt.colorbar(map3, ax=ax3, fraction=.067, pad=0.04)

# disable axes
[ax.axis("off") for ax in [ax1, ax2, ax3]]

plt.tight_layout(w_pad=0)
plt.show()

Object-mask background image correction

In some cases, using only the border of the phase image for background correction might not be enough. To increase the area of the background image, it is possible to mask only the cell area. The qpsphere package provides the convenience method qpsphere.cnvnc.bg_phase_mask_for_qpi() which computes the background phase mask based on the position and radius of an automatically detected spherical phase object. The sized of the mask can be tuned with the radial_clearance parameter.

Note that the various methods used in the examples for determining such a phase mask can be combined. Also note that before applying the method discussed here, an initial background correction might be necessary.

_images/background_mask_sphere.jpg

background_mask_sphere.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import matplotlib.pylab as plt
import numpy as np
# The data are stored in a .jpg file (lossy compression).
# If `PIL` is not found, try installing the `pillow` package.
from PIL import Image
import qpimage
import qpsphere

edata = np.array(Image.open("./data/hologram_cell_curved_bg.jpg"))

# create QPImage instance
qpi = qpimage.QPImage(data=edata,
                      which_data="hologram",
                      meta_data={"medium index": 1.335,
                                 "wavelength": 550e-9,
                                 "pixel size": 0.107e-6})
pha0 = qpi.pha

# determine the position of the cell (takes a while)
mask = qpsphere.cnvnc.bg_phase_mask_for_qpi(qpi=qpi,
                                            r0=7e-6,
                                            method="edge",
                                            model="projection",
                                            radial_clearance=1.15)

# background correction using polynomial and mask
qpi.compute_bg(which_data=["phase"],
               fit_offset="fit",
               fit_profile="poly2o",
               from_mask=mask,
               )
pha_corr = qpi.pha

# plot phase data
fig = plt.figure(figsize=(8, 3.3))

phakw = {"cmap": "viridis",
         "interpolation": "bicubic",
         "vmin": pha_corr.min(),
         "vmax": pha_corr.max()}

ax1 = plt.subplot(131, title="input phase")
map1 = ax1.imshow(pha0, **phakw)
plt.colorbar(map1, ax=ax1, fraction=.067, pad=0.04)

ax2 = plt.subplot(132, title="background phase mask")
map2 = ax2.imshow(1.*mask, cmap="gray_r")
plt.colorbar(map2, ax=ax2, fraction=.067, pad=0.04)

ax3 = plt.subplot(133, title="polynomial correction")
map3 = ax3.imshow(pha_corr, **phakw)
plt.colorbar(map3, ax=ax3, fraction=.067, pad=0.04)

# disable axes
[ax.axis("off") for ax in [ax1, ax2, ax3]]

plt.tight_layout(w_pad=0)
plt.show()

Digital hologram of a single cell

This example illustrates how qpimage can be used to analyze digital holograms. The hologram of a single myeloid leukemia cell (HL60) shown was recorded using digital holographic microscopy (DHM). Because the phase-retrieval method used in DHM is based on the discrete Fourier transform, there always is a residual background phase tilt which must be removed for further image analysis. The setup used for recording these data is described in reference [SSM+15], which also contains a description of the hologram-to-phase conversion and phase background correction algorithms which qpimage is based on.

_images/hologram_cell.jpg

hologram_cell.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import matplotlib
import matplotlib.pylab as plt
import numpy as np
import qpimage

# load the experimental data
edata = np.load("./data/hologram_cell.npz")

# create QPImage instance
qpi = qpimage.QPImage(data=edata["data"],
                      bg_data=edata["bg_data"],
                      which_data="hologram",
                      # This parameter allows to pass arguments to the
                      # hologram-analysis algorithm of qpimage.
                      # (see qpimage.holo.get_field)
                      holo_kw={
                               # For this hologram, the "smooth disk"
                               # filter yields the best trade-off
                               # between interference from the central
                               # band and image resolution.
                               "filter_name": "smooth disk",
                               # As can be seen in the hologram image,
                               # the sidebands are not positioned at
                               # an angle of 45° from the central band.
                               # If the filter size is 1/3 (default),
                               # the central band introduces line-
                               # artifacts to the reconstructed image.
                               "filter_size": 1/4
                               }
                      )

amp0 = qpi.amp
pha0 = qpi.pha

# background correction
qpi.compute_bg(which_data=["amplitude", "phase"],
               fit_offset="fit",
               fit_profile="tilt",
               border_px=5,
               )

# plot the properties of `qpi`
fig = plt.figure(figsize=(8, 10))

matplotlib.rcParams["image.interpolation"] = "bicubic"
holkw = {"cmap": "gray",
         "vmin": 0,
         "vmax": 200}

ax1 = plt.subplot(321, title="cell hologram")
map1 = ax1.imshow(edata["data"], **holkw)
plt.colorbar(map1, ax=ax1, fraction=.046, pad=0.04)

ax2 = plt.subplot(322, title="bg hologram")
map2 = ax2.imshow(edata["bg_data"], **holkw)
plt.colorbar(map2, ax=ax2, fraction=.046, pad=0.04)

ax3 = plt.subplot(323, title="input phase [rad]")
map3 = ax3.imshow(pha0)
plt.colorbar(map3, ax=ax3, fraction=.046, pad=0.04)

ax4 = plt.subplot(324, title="input amplitude")
map4 = ax4.imshow(amp0, cmap="gray")
plt.colorbar(map4, ax=ax4, fraction=.046, pad=0.04)

ax5 = plt.subplot(325, title="corrected phase [rad]")
map5 = ax5.imshow(qpi.pha)
plt.colorbar(map5, ax=ax5, fraction=.046, pad=0.04)

ax6 = plt.subplot(326, title="corrected amplitude")
map6 = ax6.imshow(qpi.amp, cmap="gray")
plt.colorbar(map6, ax=ax6, fraction=.046, pad=0.04)

# disable axes
[ax.axis("off") for ax in [ax1, ax2, ax3, ax4, ax5, ax6]]

plt.tight_layout()
plt.show()

Hologram filter choice

There are several parameters that influence the quality of phase and amplitude data retrieved from holograms. This example demonstrates the advantages and disadvantages of a three hologram filters in qpimage. For more information, please have a look at qpimage.holo.get_field().

Several observations can be made:

  • There appears to be a “bleed-through” of phase data into the amplitude data.

  • A (sharp) disk filter introduces ringing artifacts in the amplitude and phase images.

  • A smooth disk filter does not lead to such artifacts, but a dark halo is introduced around the coins in the amplitude image.

  • The amplitude reconstruction with the gaussian filter does not exhibit the dark halo but, due to blurring, reveals less details.

To correctly interpret the data shown, please note that:

  • This is a simulated hologram with no central band. For real data, the “filter_size” parameter also affects the reconstruction quality. Contributions from the central band can lead to strong artifacts. A balance between high resolution (large filter size) and small contributions from the central band (small filter size) usually has be found.

  • It is not trivial to compare a gaussian filter with a disk filter in terms of filter size (sigma vs. radius). The gaussian filter takes into account larger frequencies and suppresses low frequencies. In qpimage, the actual gaussian filter size is chosen such that the resolution approximately matches that of the disk filter with a corresponding radius. In general however, the filter size parameter has to be examined when comparing the two.

  • There is an inherent loss of information (resolution) in the holographic reconstruction process. The side band is isolated with a low-pass filter in Fourier space. The size and shape of this filter determine the resolution of the phase and amplitude images. As a result, the level of detail of all reconstructions shown is lower than that of the original images.

_images/hologram_filters.jpg

hologram_filters.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import matplotlib.pylab as plt
import numpy as np
import qpimage
from skimage import color, data

# image of a galaxy recorded with the Hubble telescope
img1 = color.rgb2grey(data.hubble_deep_field())[354:504, 70:220]
# image of a coin
img2 = color.rgb2grey(data.coins())[150:300, 70:220]

pha = img1/img1.max() * 2 * np.pi
amp = img2/img2.mean()

# create a hologram
x, y = np.mgrid[0:150, 0:150]
hologram = 2 * amp * np.cos(-2 * (x + y) + pha)

filters = ["disk", "smooth disk", "gauss"]
qpis = []

for filter_name in filters:
    qpi = qpimage.QPImage(data=hologram,
                          which_data="hologram",
                          holo_kw={"filter_size": .5,
                                   "filter_name": filter_name})
    qpis.append(qpi)

fig = plt.figure(figsize=(8, 16))

phakw = {"interpolation": "bicubic",
         "cmap": "viridis",
         "vmin": pha.min(),
         "vmax": pha.max(),
         }

ampkw = {"interpolation": "bicubic",
         "cmap": "gray",
         "vmin": amp.min(),
         "vmax": amp.max()
         }

numrows = len(filters) + 1

plt.subplot(numrows, 2, 1, title="original phase")
plt.imshow(pha, **phakw)

ax2 = plt.subplot(numrows, 2, 2, title="original amplitude")
plt.imshow(amp, **ampkw)

for ii in range(len(filters)):
    # phase
    plt.subplot(numrows, 2, 2*ii+3, title=filters[ii]+" filter")
    plt.imshow(qpis[ii].pha, **phakw)
    # amplitude
    plt.subplot(numrows, 2, 2*ii+4, title=filters[ii]+" filter")
    plt.imshow(qpis[ii].amp, **ampkw)

# disable axes
for ax in fig.get_axes():
    ax.axis("off")

plt.tight_layout()
plt.show()