.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/2_kernels/Gaussian_RK.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. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_2_kernels_Gaussian_RK.py: .. _gaussianRK: A Gaussian Reproducing Kernel ================================================ This page is not only meant to show you how to use the Gaussian Reproducing Kernel (GRK) but also, to provide you with everything you need to know about kernel operators in the Demeter library: how they are used and the available options. If you are familiar with reproducing kernels on point cloud data, you might be surprised. Since we are working with images, it is as if we had a point cloud where each pixel is a Dirac with its weight being the pixel intensity. As a result, there is a large number of closely spaced points. This is why, instead of using a classical kernel that computes the distance between points, we perform the same operation using a convolution. In metamorphosis, the RK is used to parameterize the space of acceptable vector fields $V$ through its norm, and is used to compute the vector field from the momentum. .. math:: v = - K_{\sigma} ( p \nabla I) where $p$ is the momentum, $I$ the image and $K_{\sigma}$ the Gaussian RK operator. we also call $p \nabla I$ the field momentum. .. GENERATED FROM PYTHON SOURCE LINES 27-28 Importing the necessary libraries .. GENERATED FROM PYTHON SOURCE LINES 28-38 .. code-block:: Python import numpy as np from demeter.constants import DLT_KW_IMAGE import matplotlib.pyplot as plt import torch import demeter.utils.reproducing_kernels as rk import demeter.utils.torchbox as tb import demeter.utils.bspline as bs import demeter.metamorphosis as mt .. GENERATED FROM PYTHON SOURCE LINES 39-46 To use the Gaussian RK, we need to choose a sigma. The larger the sigma, the more the RK will smooth the vector fields. The smaller the sigma, the more the RK will keep the high frequency information of the vector fields. We choose sigmas in pixel units, meaning that the gaussian will reach a value of 0.01 at about a distance of 3 sigma pixels from the center. We define the sigma as a tuple $(\sigma_x,\sigma_y)$ (or $(\sigma_x,\sigma_y,\sigma_y)$ in 3d) to allow for anisotropic smoothing and to give the dimensions of the image. .. GENERATED FROM PYTHON SOURCE LINES 46-53 .. code-block:: Python sigma = (4,2) kernelOp = rk.GaussianRKHS(sigma) print(kernelOp) kernelOp.plot() .. image-sg:: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_001.png :alt: kernel sigma: 2, Gaussian Kernel $\sigma$=(4, 2), kernel sigma: 4 :srcset: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none GaussianRKHS,2D filter :fft_filter, sigma :(4, 2) kernel_size :(1, 49, 25) kernel_reach :6 normalized :True .. GENERATED FROM PYTHON SOURCE LINES 54-65 In the cell above, we initialized a Gaussian RK operator with sigma = (10,5). You can print the operator to see its parameters. You can see the kernel size and its reach. The kernel size is the size of the kernel in pixel units. The reach is the distance at which we want to truncate the kernel. A reach of 3 means will be truncated at 3 sigma from center to edge. This operator in meant to smooth the vector field momentum in the Metamorphosis geodesic equation. Before, to use it in a Metamorphosis registration, let's visualize its effect on a synthetic vector field. .. GENERATED FROM PYTHON SOURCE LINES 65-81 .. code-block:: Python # Create a synthetic vector field H,W = 30,30 s = .5 # Random control matrix cms = torch.rand((2,int(H*s),int(W*s)),dtype=torch.float) * 2 - 1 cms *= 5 field = bs.field2D_bspline(cms, n_pts = (H,W), degree=(1,1), dim_stack=-1)[None] smoothed_field = tb.im2grid(kernelOp( tb.grid2im(field) )) ic(field.shape, smoothed_field.shape) fig,ax = plt.subplots(1,2,figsize=(10,5)) tb.quiver_plot(field, title='Original field', step=1,ax=ax[0],check_diffeo=True) tb.quiver_plot(smoothed_field, title='Smoothed field', step=1,ax=ax[1],check_diffeo=True) plt.show() .. image-sg:: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_002.png :alt: Gaussian RK :srcset: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 82-102 ##################################################################### Note: that we used the `im2grid` and `grid2im` functions to convert the field from image to grid representation and vice versa. The grid representation is a tensor of shape (1,H,W,2) where H and W are the height and width of the image. You might find it strange that we pass the field in image representation to the kernel operator. This mostly for efficiency reasons, see `pytorch stride ` for more information. In the plot above, we can see the effect of the Gaussian RK on a synthetic vector field. We show the "diffeomorphismness" of the field by plotting the sign of the Jacobian of the field. Green means that the Jacobian is positive, red means that the Jacobian is negative, and the image must be all green for the deformation to be a diffeomorphism. We can see that the smoothed field is a diffemorphism ! Registation ^^^^^^^^^^^^^^^^^^^^^ Good, now lets try on a real use case. First let's load a source and image target: Try to change the size and the proposed sigma to see their effects on the registration results. Feel free to try different values that the one proposed. .. GENERATED FROM PYTHON SOURCE LINES 102-121 .. code-block:: Python device = 'cpu' if torch.cuda.is_available(): device = 'cuda:0' print('device used :',device) size = (300,300) # size = (300,150) source = tb.reg_open('01',size = size) target = tb.reg_open('star_sr5_r10',size = size) fig, ax = plt.subplots(1,3,figsize=(10,5)) ax[0].imshow(source[0,0],**DLT_KW_IMAGE) ax[0].set_title('source') ax[1].imshow(target[0,0],**DLT_KW_IMAGE) ax[1].set_title('target') ax[2].imshow(tb.imCmp(source,target,'compose'),origin='lower') ax[2].set_title('superposition of source and target') plt.show() .. image-sg:: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_003.png :alt: source, target, superposition of source and target :srcset: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none device used : cpu .. GENERATED FROM PYTHON SOURCE LINES 122-124 Now let's try to register the source to the target using LDDMM with differents sigma values. .. GENERATED FROM PYTHON SOURCE LINES 124-143 .. code-block:: Python source = source.to(device) target = target.to(device) # sigma = (10,5) # wonky result sigma = (10,10) # good result # sigma = (20,20) # cannot get details. kernelOp = rk.GaussianRKHS(sigma) print(kernelOp) mr = mt.lddmm(source,target,0, kernelOperator = kernelOp, n_iter = 10, cost_cst=1, grad_coef=1, integration_steps=15, dx_convention='square', ) .. rst-class:: sphx-glr-script-out .. code-block:: none GaussianRKHS,2D filter :fft_filter, sigma :(10, 10) kernel_size :(1, 121, 121) kernel_reach :6 normalized :True Progress: [##--------] 20.00% (Ssd : ,655.1315). Progress: [###-------] 30.00% (Ssd : ,285.0678). Progress: [####------] 40.00% (Ssd : ,157.9997). Progress: [#####-----] 50.00% (Ssd : ,126.9168). Progress: [######----] 60.00% (Ssd : ,121.4015). Progress: [#######---] 70.00% (Ssd : ,117.6779). Progress: [########--] 80.00% (Ssd : ,115.2220). Progress: [#########-] 90.00% (Ssd : ,104.2477). Progress: [##########] 100.00% Done... (Ssd : , 97.9167). Computation of forward done in 0:00:27s and 0.860cents s Computation of lddmm done in 0:00:27s and 0.860cents s .. GENERATED FROM PYTHON SOURCE LINES 144-146 Integration and convergence plot. Then we compare the source, target and registered source images .. GENERATED FROM PYTHON SOURCE LINES 146-147 .. code-block:: Python mr.plot() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_004.png :alt: Lambda = 1 rho = 1.0 :srcset: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_005.png :alt: source, target, Integrated source image, comparaison deformed image with target :srcset: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none ((
, array([, ], dtype=object)), (
, array([[, ], [, ]], dtype=object))) .. GENERATED FROM PYTHON SOURCE LINES 148-150 This plot shows the deformation as a grid and arrows. along with the deformed image. .. GENERATED FROM PYTHON SOURCE LINES 150-153 .. code-block:: Python mr.plot_deform() plt.show() .. image-sg:: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_006.png :alt: diffeo = tensor(True) :srcset: /auto_examples/2_kernels/images/sphx_glr_Gaussian_RK_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 30.114 seconds) .. _sphx_glr_download_auto_examples_2_kernels_Gaussian_RK.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Gaussian_RK.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Gaussian_RK.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Gaussian_RK.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_