Inpainting missing data

Missing data in an image can be an issue, especially when one wants to perform Fourier analysis. This tutorial explains how to fill-up missing pixels with values which looks “realistic” and introduce as little perturbation as possible for subsequent analysis. The user should keep the mask nearby and only consider the values of actual pixels and never the one inpainted.

This tutorial will use fully synthetic data to allow comparison between actual (syntetic) data with inpainted values.

The first part of the tutorial is about the generation of a challenging 2D diffraction image with realistic noise and to describe the metric used, then comes the actual tutorial on how to use the inpainting. Finally a benchmark is used based on the metric determined.

Creation of the image

A realistic challenging image should contain:

• Bragg peak rings. We chose LaB6 as guinea-pig, with very sharp peaks, at the limit of the resolution of the detector
• Some amorphous content
• strong polarization effect
• Poissonian noise

One image will be generated but then multiple ones with different noise to discriminate the effect of the noise from other effects.

%pylab nbagg

Populating the interactive namespace from numpy and matplotlib

import pyFAI
from pyFAI.gui import jupyter
import pyFAI.test.utilstest
from pyFAI.calibrant import get_calibrant
import time
start_time = time.time()

WARNING:test_integrate_widget:pyFAI.integrate_widget tests disabled (DISPLAY env. variable not set)

detector = pyFAI.detector_factory("Pilatus2MCdTe")
mask = detector.mask.copy()
nomask = numpy.zeros_like(mask)
detector.mask=nomask
ai = pyFAI.AzimuthalIntegrator(detector=detector)
ai.setFit2D(200, 200, 200)
ai.wavelength = 3e-11
print(ai)

Detector Pilatus CdTe 2M     PixelSize= 1.720e-04, 1.720e-04 m
Wavelength= 3.000000e-11m
SampleDetDist= 2.000000e-01m        PONI= 3.440000e-02, 3.440000e-02m       rot1=0.000000  rot2= 0.000000  rot3= 0.000000 rad
DirectBeamDist= 200.000mm   Center: x=200.000, y=200.000 pix        Tilt=0.000 deg  tiltPlanRotation= 0.000 deg

LaB6 = get_calibrant("LaB6")
LaB6.wavelength = ai.wavelength
print(LaB6)
r = ai.array_from_unit(unit="q_nm^-1")
decay_b = numpy.exp(-(r-50)**2/2000)
bragg = LaB6.fake_calibration_image(ai, Imax=1e4, W=1e-6) * ai.polarization(factor=1.0) * decay_b
decay_a = numpy.exp(-r/100)
amorphous = 1000*ai.polarization(factor=1.0)*ai.solidAngleArray() * decay_a
img_nomask = bragg + amorphous
#Not the same noise function for all images two images
img_nomask = numpy.random.poisson(img_nomask)
img_nomask2 = numpy.random.poisson(img_nomask)
img = numpy.random.poisson(img_nomask)
img[numpy.where(mask)] = -1
fig,ax = subplots(1,2, figsize=(10,5))
jupyter.display(img=img, label="With mask", ax=ax[0])
jupyter.display(img=img_nomask, label="Without mask", ax=ax[1])

LaB6 Calibrant at wavelength 3e-11

<IPython.core.display.Javascript object>

<matplotlib.axes._subplots.AxesSubplot at 0x7f4cb1a061d0>

Note the aliassing effect on the displayed images.

We will measure now the effect after 1D intergeration. We do not correct for polarization on purpose to highlight the defect one wishes to whipe out. We use a R-factor to describe the quality of the 1D-integrated signal.

wo = ai.integrate1d(img_nomask, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
wo2 = ai.integrate1d(img_nomask2, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
wm = ai.integrate1d(img, 2000, unit="q_nm^-1", method="splitpixel", mask=mask, radial_range=(0,210))
ax = jupyter.plot1d(wm , label="with_mask")
ax.plot(*wo, label="without_mask")
ax.plot(*wo2, label="without_mask2")
ax.plot(wo.radial, wo.intensity-wm.intensity, label="delta")
ax.plot(wo.radial, wo.intensity-wo2.intensity, label="relative-error")
ax.legend()
print("Between masked and non masked image R= %s"%pyFAI.test.utilstest.Rwp(wm,wo))
print("Between two different non-masked images R'= %s"%pyFAI.test.utilstest.Rwp(wo2,wo))

<IPython.core.display.Javascript object>

Between masked and non masked image R= 5.67315744311
Between two different non-masked images R'= 0.175106399221

# Effect of the noise on the delta image
fig, ax = subplots()
jupyter.display(img=img_nomask-img_nomask2, label="Delta due to noise", ax=ax)
ax.figure.colorbar(ax.images[0])

<IPython.core.display.Javascript object>

<matplotlib.colorbar.Colorbar at 0x7f4cb187ad30>

Inpainting

This part describes how to paint the missing pixels for having a “natural-looking image”. The delta image contains the difference with the original image

#Inpainting:
inpainted = ai.inpainting(img, mask=mask, poissonian=True)
fig, ax = subplots(1, 2, figsize=(12,5))
jupyter.display(img=inpainted, label="Inpainted", ax=ax[0])
jupyter.display(img=img_nomask-inpainted, label="Delta", ax=ax[1])
ax[1].figure.colorbar(ax[1].images[0])

<IPython.core.display.Javascript object>

<matplotlib.colorbar.Colorbar at 0x7f4ca5031668>

# Comparison of the inpained image with the original one:
wm = ai.integrate1d(inpainted, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
wo = ai.integrate1d(img_nomask, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
ax = jupyter.plot1d(wm , label="inpainted")
ax.plot(*wo, label="without_mask")
ax.plot(wo.radial, wo.intensity-wm.intensity, label="delta")
ax.legend()
print("R= %s"%pyFAI.test.utilstest.Rwp(wm,wo))

<IPython.core.display.Javascript object>

R= 1.29109163117

One can see by zooming in that the main effect on inpainting is a broadening of the signal in the inpainted region. This could (partially) be adressed by increasing the number of radial bins used in the inpainting.

Benchmarking and optimization of the parameters

The parameter set depends on the detector, the experiment geometry and the type of signal on the detector. Finer detail require finer slicing.

#Basic benchmarking of execution time for default options:
% timeit inpainted = ai.inpainting(img, mask=mask)

1 loop, best of 3: 1.04 s per loop

wo = ai.integrate1d(img_nomask, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
for j in ("csr", "csr_nosplit"):
    for k in (512,1024,2048, 4096):
        ai.reset()
        for i in range(10):
            inpainted = ai.inpainting(img, mask=mask, poissonian=True, method=j, npt_rad=k, grow_mask=i)
            wm = ai.integrate1d(inpainted, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
            print("method: %s npt_rad=%i grow=%i; R= %s"%(j, k, i,pyFAI.test.utilstest.Rwp(wm,wo)))

method: csr npt_rad=512 grow=0; R= 3.16527395132
method: csr npt_rad=512 grow=1; R= 2.90373554002
method: csr npt_rad=512 grow=2; R= 2.7075280075
method: csr npt_rad=512 grow=3; R= 2.65327421322
method: csr npt_rad=512 grow=4; R= 2.61776326717
method: csr npt_rad=512 grow=5; R= 2.57439873144
method: csr npt_rad=512 grow=6; R= 2.54164759702
method: csr npt_rad=512 grow=7; R= 2.54679849484
method: csr npt_rad=512 grow=8; R= 2.53635706674
method: csr npt_rad=512 grow=9; R= 2.52623386268
method: csr npt_rad=1024 grow=0; R= 1.69183876823
method: csr npt_rad=1024 grow=1; R= 1.31855606592
method: csr npt_rad=1024 grow=2; R= 1.31247003017
method: csr npt_rad=1024 grow=3; R= 1.29097440452
method: csr npt_rad=1024 grow=4; R= 1.27074644319
method: csr npt_rad=1024 grow=5; R= 1.27238306777
method: csr npt_rad=1024 grow=6; R= 1.27167179838
method: csr npt_rad=1024 grow=7; R= 1.26518530245
method: csr npt_rad=1024 grow=8; R= 1.26533812218
method: csr npt_rad=1024 grow=9; R= 1.26587569844
method: csr npt_rad=2048 grow=0; R= 0.895867687266
method: csr npt_rad=2048 grow=1; R= 0.645103411948
method: csr npt_rad=2048 grow=2; R= 0.630137579053
method: csr npt_rad=2048 grow=3; R= 0.624736974026
method: csr npt_rad=2048 grow=4; R= 0.62869452797
method: csr npt_rad=2048 grow=5; R= 0.628601472572
method: csr npt_rad=2048 grow=6; R= 0.631538002653
method: csr npt_rad=2048 grow=7; R= 0.629878287879
method: csr npt_rad=2048 grow=8; R= 0.622824815248
method: csr npt_rad=2048 grow=9; R= 0.631178231084
method: csr npt_rad=4096 grow=0; R= 0.564353619636
method: csr npt_rad=4096 grow=1; R= 0.399846196196
method: csr npt_rad=4096 grow=2; R= 0.410612097165
method: csr npt_rad=4096 grow=3; R= 0.410295214907
method: csr npt_rad=4096 grow=4; R= 0.410335613233
method: csr npt_rad=4096 grow=5; R= 0.406933039484
method: csr npt_rad=4096 grow=6; R= 0.406117616454
method: csr npt_rad=4096 grow=7; R= 0.404649982091
method: csr npt_rad=4096 grow=8; R= 0.407870646779
method: csr npt_rad=4096 grow=9; R= 0.41036433231

WARNING:pyFAI.splitBBoxCSR:Pixel splitting desactivated !

method: csr_nosplit npt_rad=512 grow=0; R= 3.84746010082
method: csr_nosplit npt_rad=512 grow=1; R= 2.95385230424
method: csr_nosplit npt_rad=512 grow=2; R= 2.67580240631
method: csr_nosplit npt_rad=512 grow=3; R= 2.6126805269
method: csr_nosplit npt_rad=512 grow=4; R= 2.58362811691
method: csr_nosplit npt_rad=512 grow=5; R= 2.60564448801
method: csr_nosplit npt_rad=512 grow=6; R= 2.58621310627
method: csr_nosplit npt_rad=512 grow=7; R= 2.541826826
method: csr_nosplit npt_rad=512 grow=8; R= 2.5206857231
method: csr_nosplit npt_rad=512 grow=9; R= 2.52451588041

WARNING:pyFAI.splitBBoxCSR:Pixel splitting desactivated !

method: csr_nosplit npt_rad=1024 grow=0; R= 2.75213259282
method: csr_nosplit npt_rad=1024 grow=1; R= 1.4384241414
method: csr_nosplit npt_rad=1024 grow=2; R= 1.33379422417
method: csr_nosplit npt_rad=1024 grow=3; R= 1.28940297956
method: csr_nosplit npt_rad=1024 grow=4; R= 1.27720031756
method: csr_nosplit npt_rad=1024 grow=5; R= 1.28777067352
method: csr_nosplit npt_rad=1024 grow=6; R= 1.16986522429
method: csr_nosplit npt_rad=1024 grow=7; R= 1.13872560834
method: csr_nosplit npt_rad=1024 grow=8; R= 1.15033670027
method: csr_nosplit npt_rad=1024 grow=9; R= 1.15867431745

WARNING:pyFAI.splitBBoxCSR:Pixel splitting desactivated !

method: csr_nosplit npt_rad=2048 grow=0; R= 2.79811477948
method: csr_nosplit npt_rad=2048 grow=1; R= 1.22533548525
method: csr_nosplit npt_rad=2048 grow=2; R= 1.12880254637
method: csr_nosplit npt_rad=2048 grow=3; R= 1.02906086137
method: csr_nosplit npt_rad=2048 grow=4; R= 0.991830558658
method: csr_nosplit npt_rad=2048 grow=5; R= 0.907918584422
method: csr_nosplit npt_rad=2048 grow=6; R= 0.841481199476
method: csr_nosplit npt_rad=2048 grow=7; R= 0.805653872599
method: csr_nosplit npt_rad=2048 grow=8; R= 0.786059385048
method: csr_nosplit npt_rad=2048 grow=9; R= 0.765824492478

WARNING:pyFAI.splitBBoxCSR:Pixel splitting desactivated !

method: csr_nosplit npt_rad=4096 grow=0; R= 2.66329347592
method: csr_nosplit npt_rad=4096 grow=1; R= 1.19262995518
method: csr_nosplit npt_rad=4096 grow=2; R= 1.15682781524
method: csr_nosplit npt_rad=4096 grow=3; R= 1.15083356666
method: csr_nosplit npt_rad=4096 grow=4; R= 1.12555533128
method: csr_nosplit npt_rad=4096 grow=5; R= 1.10571482782
method: csr_nosplit npt_rad=4096 grow=6; R= 1.06182288301
method: csr_nosplit npt_rad=4096 grow=7; R= 1.0336445245
method: csr_nosplit npt_rad=4096 grow=8; R= 0.992978393303
method: csr_nosplit npt_rad=4096 grow=9; R= 0.962156509678

#Inpainting, best solution found:
ai.reset()
%time inpainted = ai.inpainting(img, mask=mask, poissonian=True, method="csr", npt_rad=4096, grow_mask=5)
fig, ax = subplots(1, 2, figsize=(12,5))
jupyter.display(img=inpainted, label="Inpainted", ax=ax[0])
jupyter.display(img=img_nomask-inpainted, label="Delta", ax=ax[1])
ax[1].figure.colorbar(ax[1].images[0])

CPU times: user 3.84 s, sys: 740 ms, total: 4.58 s
Wall time: 1.49 s

<IPython.core.display.Javascript object>

<matplotlib.colorbar.Colorbar at 0x7f4ca4f16be0>

# Comparison of the inpained image with the original one:
wm = ai.integrate1d(inpainted, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
wo = ai.integrate1d(img_nomask, 2000, unit="q_nm^-1", method="splitpixel", radial_range=(0,210))
ax = jupyter.plot1d(wm , label="inpainted")
ax.plot(*wo, label="without_mask")
ax.plot(wo.radial, wo.intensity-wm.intensity, label="delta")
ax.legend()
print("R= %s"%pyFAI.test.utilstest.Rwp(wm,wo))

<IPython.core.display.Javascript object>

R= 0.40783662849

Conclusion

Inpainting is one of the only solution to fill up the gaps in detector when Fourier analysis is needed. This tutorial explains basically how this is possible using the pyFAI library and how to optimize the parameter set for inpainting. The result may greatly vary with detector position and tilt and the kind of signal (amorphous or more spotty).

print("Execution time: %.3fs"%(time.time()-start_time))

Execution time: 106.425s