#### NOTE New to DeepInverse? Get started with the basics with the [5 minute quickstart tutorial](https://deepinv.org/auto_examples/basics/demo_quickstart.html.md#sphx-glr-auto-examples-basics-demo-quickstart-py).. # Tour of blur operators This example provides a tour of 2D blur operators in DeepInverse. In particular, we show how to use DiffractionBlurs (Fresnel diffraction), motion blurs and space varying blurs. ```Python import torch import deepinv as dinv from deepinv.utils.plotting import plot from deepinv.utils import load_example ``` --- First, let’s load some test images. ```Python dtype = torch.float32 device = "cuda" if torch.cuda.is_available() else "cpu" img_size = (173, 125) x_rgb = load_example( "CBSD_0010.png", grayscale=False, device=device, dtype=dtype, img_size=img_size ) x_gray = load_example( "barbara.jpeg", grayscale=True, device=device, dtype=dtype, img_size=img_size ) # Next, set the global random seed from pytorch to ensure reproducibility of the example. torch.manual_seed(0) torch.cuda.manual_seed(0) ``` We are now ready to explore the different blur operators. ## Convolution Basics The class [`deepinv.physics.Blur`](https://deepinv.org/api/stubs/deepinv.physics.Blur.html.md#deepinv.physics.Blur) implements convolution operations with kernels. For instance, here is the convolution of a grayscale image with a grayscale filter: ```Python filter_0 = dinv.physics.functional.gaussian_blur(sigma=(2, 0.1), angle=0.0) physics = dinv.physics.Blur(filter_0, device=device) y = physics(x_gray) plot( [x_gray, filter_0, y], titles=["signal", "filter", "measurement"], suptitle="Grayscale convolution", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` When a single channel filter is used, all channels are convolved with the same filter: ```Python physics = dinv.physics.Blur(filter_0, device=device) y = physics(x_rgb) plot( [x_rgb, filter_0, y], titles=["signal", "filter", "measurement"], suptitle="RGB image + grayscale filter convolution", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` By default, the boundary conditions are `'valid'`, but other options among (`'circular'`, `'reflect'`, `'replicate'`) are possible: ```Python physics = dinv.physics.Blur(filter_0, padding="reflect", device=device) y = physics(x_rgb) plot( [x_rgb, filter_0, y], titles=["signal", "filter", "measurement"], suptitle="Reflection boundary conditions", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` For circular boundary conditions, an FFT implementation is also available. It is slower that [`deepinv.physics.Blur`](https://deepinv.org/api/stubs/deepinv.physics.Blur.html.md#deepinv.physics.Blur), but inherits from [`deepinv.physics.DecomposablePhysics`](https://deepinv.org/api/stubs/deepinv.physics.DecomposablePhysics.html.md#deepinv.physics.DecomposablePhysics), so that the pseudo-inverse and regularized inverse are computed faster and more accurately. ```Python physics = dinv.physics.BlurFFT(img_size=x_rgb[0].shape, filter=filter_0, device=device) y = physics(x_rgb) plot( [x_rgb, filter_0, y], titles=["signal", "filter", "measurement"], suptitle="FFT convolution with circular boundary conditions", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` One can also change the blur filter in the forward pass as follows: ```Python filter_90 = dinv.physics.functional.gaussian_blur(sigma=(2, 0.1), angle=90.0).to( device=device, dtype=dtype ) y = physics(x_rgb, filter=filter_90) plot( [x_rgb, filter_90, y], titles=["signal", "filter", "measurement"], suptitle="Changing the filter on the fly", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` When applied to a new image, the last filter is used: ```Python y = physics(x_gray, filter=filter_90) plot( [x_gray, filter_90, y], titles=["signal", "filter", "measurement"], suptitle="Effect of on the fly change is persistent", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` We can also define color filters. In that situation, each channel is convolved with the corresponding channel of the filter: ```Python psf_size = 9 filter_rgb = torch.zeros((1, 3, psf_size, psf_size), device=device, dtype=dtype) filter_rgb[:, 0, :, psf_size // 2 : psf_size // 2 + 1] = 1.0 / psf_size filter_rgb[:, 1, psf_size // 2 : psf_size // 2 + 1, :] = 1.0 / psf_size filter_rgb[:, 2, ...] = ( torch.diag(torch.ones(psf_size, device=device, dtype=dtype)) / psf_size ) y = physics(x_rgb, filter=filter_rgb) plot( [x_rgb, filter_rgb, y], titles=["signal", "Colour filter", "measurement"], suptitle="Color image + color filter convolution", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` ## Blur generators More advanced kernel generation methods are provided with the toolbox thanks to the [`deepinv.physics.generator.PSFGenerator`](https://deepinv.org/api/stubs/deepinv.physics.generator.PSFGenerator.html.md#deepinv.physics.generator.PSFGenerator). In particular, motion blurs generators are implemented. ### Motion blur generators ```Python from deepinv.physics.generator import MotionBlurGenerator ``` In order to generate motion blur kernels, we just need to instantiate a generator with specific the psf size. In turn, motion blurs can be generated on the fly by calling the `step()` method. Let’s illustrate this now and generate 3 motion blurs. First, we instantiate the generator: ```Python psf_size = 31 motion_generator = MotionBlurGenerator((psf_size, psf_size), device=device, dtype=dtype) ``` To generate new filters, we call the step() function: ```Python filters = motion_generator.step(batch_size=3) # the `step()` function returns a dictionary: print(filters.keys()) plot( [f for f in filters["filter"]], suptitle="Examples of randomly generated motion blurs", ) ``` ```none dict_keys(['filter']) /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` Other options, such as the regularity and length of the blur trajectory can also be specified: ```Python motion_generator = MotionBlurGenerator( (psf_size, psf_size), l=0.6, sigma=1, device=device, dtype=dtype ) filters = motion_generator.step(batch_size=3) plot([f for f in filters["filter"]], suptitle="Different length and regularity") ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` ### Diffraction blur generators We also implemented diffraction blurs obtained through Fourier optics and definition of the psf through the pupil plane expanded in Zernike polynomials ```Python from deepinv.physics.generator import DiffractionBlurGenerator diffraction_generator = DiffractionBlurGenerator( (psf_size, psf_size), device=device, dtype=dtype ) ``` Then, to generate new filters, it suffices to call the step() function as follows: ```Python filters = diffraction_generator.step(batch_size=3) ``` In this case, the `step()` function returns a dictionary containing the filters, their pupil function, Zernike coefficients and cutoff frequency: ```Python print(filters.keys()) # Note that we use **0.5 to increase the image dynamics plot( [f for f in filters["filter"] ** 0.5], suptitle="Examples of randomly generated diffraction blurs", ) plot( [ f for f in torch.angle(filters["pupil"][:, None]) * torch.abs(filters["pupil"][:, None]) ], suptitle="Corresponding pupil phases", ) print("Coefficients of the decomposition on Zernike polynomials") print(filters["coeff"]) ``` ```none dict_keys(['filter', 'coeff', 'pupil', 'fc']) /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) Coefficients of the decomposition on Zernike polynomials tensor([[[-0.0160, 0.0017, 0.0578, -0.0089, 0.0600, -0.0567, 0.0046, -0.0250]], [[-0.0392, 0.0006, -0.0610, 0.0743, -0.0712, -0.0380, 0.0690, -0.0289]], [[-0.0091, 0.0396, -0.0102, 0.0660, 0.0314, 0.0264, -0.0037, 0.0203]]], device='cuda:0') ``` We can change the cutoff frequency (below 1/4 to respect Shannon’s sampling theorem) ```Python filters = diffraction_generator.step(batch_size=3, fc=1 / 8) plot( [f for f in filters["filter"] ** 0.5], suptitle="A different cutoff frequency", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` We provide a helper property to get the list of Zernike polynomials used in the decomposition: ```Python zernike_polynomials = diffraction_generator.zernike_polynomials print("Zernike polynomials used:\n", "\n ".join(zernike_polynomials)) # It is also possible to directly specify the Zernike decomposition. # For instance, if the pupil is null, the PSF is the Airy pattern n_zernike = len( zernike_polynomials ) # number of Zernike coefficients in the decomposition filters = diffraction_generator.step(coeff=torch.zeros(3, n_zernike, device=device)) plot( [f for f in filters["filter"][:, None] ** 0.3], suptitle="Airy pattern", ) ``` ```none Zernike polynomials used: Zernike(n = 2, m = 0) -- Defocus Zernike(n = 2, m = -2) -- Oblique Astigmatism Zernike(n = 2, m = 2) -- Vertical Astigmatism Zernike(n = 3, m = -1) -- Vertical Coma Zernike(n = 3, m = 1) -- Horizontal Coma Zernike(n = 3, m = -3) -- Vertical Trefoil Zernike(n = 3, m = 3) -- Oblique Trefoil Zernike(n = 4, m = 0) -- Primary Spherical /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` Notice that you can activate the aberrations you want in the ANSI/Noll nomenclature [https://en.wikipedia.org/wiki/Zernike_polynomials#OSA/ANSI_standard_indices](https://en.wikipedia.org/wiki/Zernike_polynomials#OSA/ANSI_standard_indices) ```Python diffraction_generator = DiffractionBlurGenerator( (psf_size, psf_size), fc=1 / 8, zernike_index=( 5, 6, ), # or equivalently zernike_index = ((2, 2), (3, -3)) for (n,m) indices index_convention="ansi", device=device, dtype=dtype, ) filters = diffraction_generator.step(batch_size=3) plot( [f for f in filters["filter"] ** 0.5], suptitle="PSF obtained with astigmatism only", ) print( "Zernike polynomials used:\n", "\n ".join(diffraction_generator.zernike_polynomials) ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) Zernike polynomials used: Zernike(n = 2, m = 2) -- Vertical Astigmatism Zernike(n = 3, m = -3) -- Vertical Trefoil ``` Generate physically consistent chromatic PSFs by passing a tuple of cutoff frequencies to `step()`. Each channel is assigned its own cutoff frequency `fc = NA * pixel_size / lambda`, which reproduces the wavelength-dependent diffraction limit. The Zernike coefficients represent wavefront errors in units of wavelength. Therefore, keeping the same physical optical path difference (OPD) across channels requires rescaling the coefficients approximately as `w(lambda) = w(lambda0) * lambda0 / lambda`. As a result, shorter wavelengths exhibit both a narrower diffraction pattern and stronger phase aberrations for the same underlying physical wavefront. ```Python diffraction_generator = DiffractionBlurGenerator( (psf_size, psf_size), device=device, dtype=dtype, ) lambdaR = 650 # in nm lambdaG = 550 # in nm lambdaB = 450 # in nm NA = 1.51 # numerical aperture of oil pixel_size = 50 # in nm fc = (NA * pixel_size / lambdaR, NA * pixel_size / lambdaG, NA * pixel_size / lambdaB) filters = diffraction_generator.step(batch_size=3, fc=fc) plot( [f for f in filters["filter"] ** 0.5], suptitle="Chromatic PSFs", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` ### Generator Mixture During training, it’s more robust to train on multiple family of operators. This can be done seamlessly with the [`deepinv.physics.generator.GeneratorMixture`](https://deepinv.org/api/stubs/deepinv.physics.generator.GeneratorMixture.html.md#deepinv.physics.generator.GeneratorMixture). ```Python from deepinv.physics.generator import GeneratorMixture torch.cuda.manual_seed(4) torch.manual_seed(6) generator = GeneratorMixture( ([motion_generator, diffraction_generator]), probs=[0.5, 0.5] ) for i in range(4): filters = generator.step(batch_size=3) plot( [f for f in filters["filter"]], suptitle=f"Random PSF generated at step {i + 1}", ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` ## Space varying blurs with Eigen PSFs by product convolution Space varying blurs are also available using [`deepinv.physics.SpaceVaryingBlur`](https://deepinv.org/api/stubs/deepinv.physics.SpaceVaryingBlur.html.md#deepinv.physics.SpaceVaryingBlur) We plot the impulse responses at different spatial locations by convolving a Dirac comb with the operator. ```Python from deepinv.physics.generator import ( ProductConvolutionBlurGenerator, ) from deepinv.physics.blur import SpaceVaryingBlur img_size = (256, 256) n_eigenpsf = 7 spacing = (32, 32) padding = "valid" batch_size = 1 delta = 16 # Now, scattered random psfs are synthesized and interpolated spatially pc_generator = ProductConvolutionBlurGenerator( psf_generator=MotionBlurGenerator((17, 17), device=device, dtype=dtype), img_size=img_size, n_eigen_psf=n_eigenpsf, spacing=spacing, padding=padding, device=device, ) params_pc = pc_generator.step(batch_size) physics = SpaceVaryingBlur(**params_pc, device=device) dirac_comb = dinv.utils.dirac_comb((1, 1) + img_size, step=delta, device=device) psf_grid = physics(dirac_comb) plot( psf_grid, titles="Space varying impulse responses", rescale_mode="clip", figsize=(5, 5), ) image = dinv.utils.load_example( "celeba_example.jpg", img_size=img_size, resize_mode="resize", device=device ) blurry_image = physics(image) plot( [image, blurry_image], titles=["Original image", "Blurry image"], rescale_mode="clip", figsize=(5, 5), ) ``` ## Space varying blur with tiles While the [`deepinv.physics.blur.SpaceVaryingBlur`](https://deepinv.org/api/stubs/deepinv.physics.SpaceVaryingBlur.html.md#deepinv.physics.SpaceVaryingBlur) physics uses global eigen-PSFs to define the space varying blur, we also provide a local version of space varying blur with the [`deepinv.physics.TiledSpaceVaryingBlur`](https://deepinv.org/api/stubs/deepinv.physics.TiledSpaceVaryingBlur.html.md#deepinv.physics.TiledSpaceVaryingBlur) class. In this case, the image is decomposed into overlapping tiles, and each tile is convolved with a different kernel. The kernels are then blended together with a unity partition to ensure a smooth transition between tiles. ```Python # This physics is particularly useful for large images when one can perform local estimation of the PSFs. from deepinv.physics.blur import TiledSpaceVaryingBlur from deepinv.physics.generator import TiledBlurGenerator img_size = (512, 512) # size of the image to blur patch_size = (128, 128) # size of the tiles on which local convolution is performed stride = (64, 64) # stride between tiles psf_generator = MotionBlurGenerator( (25, 25), device=device, dtype=dtype, ) generator = TiledBlurGenerator( psf_generator=psf_generator, patch_size=patch_size, stride=stride, device=device ) filters = generator.step(batch_size=batch_size, img_size=img_size)["filters"] physics = TiledSpaceVaryingBlur( filters=filters, patch_size=patch_size, stride=stride, device=device, use_fft=True, ) dirac_comb = dinv.utils.dirac_comb((1, 1) + img_size, step=32, device=device) y = physics(dirac_comb) plot( y.abs() ** 0.5, suptitle="Impulse responses of the tiled space varying blur", figsize=(5, 5), ) image = dinv.utils.load_example( "celeba_example.jpg", img_size=img_size, resize_mode="resize", device=device ) blurry_image = physics(image) plot( [image, blurry_image], titles=["Original image", "Blurry image"], rescale_mode="clip", figsize=(5, 5), ) ``` ```none /local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:519: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust. fig.subplots_adjust(top=0.75) ``` **Total running time of the script:** (0 minutes 2.105 seconds) 🏷 Tags: Deblurring, Microscopy, Astronomy [`Download Jupyter notebook: demo_blur_tour.ipynb`](https://deepinv.org/auto_examples/physics/_downloads/effc3c80ec8f6639c1a64d6a3d9b4abe/demo_blur_tour.ipynb) [`Download Python source code: demo_blur_tour.py`](https://deepinv.org/auto_examples/physics/_downloads/d52a422358a79844ce5583cb13687059/demo_blur_tour.py) [`Download zipped: demo_blur_tour.zip`](https://deepinv.org/auto_examples/physics/_downloads/3277dc2a8e87596e02b3badecab36654/demo_blur_tour.zip) [Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)