.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/1_registration/simpleToyExample_metamorphosis.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_1_registration_simpleToyExample_metamorphosis.py: .. _plot_simpleToyExample_metamorphosis: MetaMorphosis behaviour with different values of rho ======================================================= This example shows how the MetaMorphosis behaves with different values of rho. Let $\rho \in [0,1]$, we define the evolution of the image over time as: .. math:: \dot{I_{t}} = - \sqrt{ \rho } v_{t} \cdot \nabla I_{t} + \sqrt{ 1-\rho } z. Note: The square roots may seem surprising, but the reason will become clear at the end. We define the Hamiltonian: .. math:: H(I,p,v,z) = - (p |\dot{ I}) - \frac{1}{2} (Lv|v)_{2} - \frac{1}{2}(z,z)_{2}. By calculating the optimal trajectories we obtain the system: .. math:: \left\{ \begin{array}{rl} v &= - \sqrt{ \rho } K_{V} (p \nabla I)\\ \dot{p} &= -\sqrt{ \rho } \nabla \cdot (pv) \\ z &= \sqrt{ 1 - \rho } p \\ \dot{I} &= - \sqrt{ \rho } v_{t} \cdot \nabla I_{t} + \sqrt{ 1-\rho } z.\end{array} \right. We notice that we can rewrite: .. math:: \dot{ I_{t}} = - \sqrt{ \rho } v_{t} \cdot \nabla I_{t} + \sqrt{ 1-\rho } z = - \rho K_{V}(p_{t}\nabla I) \cdot \nabla I_{t} + (1-\rho) p_{t}. This means that the optimal dynamics is the weighted average between the creation of the field and the photometric addition controlled by the momentum $p$. .. GENERATED FROM PYTHON SOURCE LINES 37-38 Import the necessary packages .. GENERATED FROM PYTHON SOURCE LINES 38-63 .. code-block:: Python try: import sys, os # add the parent directory to the path base_path = os.path.join(os.path.dirname(os.path.abspath(__file__)),'..') sys.path.insert(0,base_path) import __init__ except NameError: pass import torch import matplotlib.pyplot as plt from demeter import * import demeter.utils.reproducing_kernels as rk import demeter.metamorphosis as mt import demeter.utils.torchbox as tb location = os.getcwd() if 'runner' in location: location = os.path.dirname(os.path.dirname(location)) EXPL_SAVE_FOLDER = os.path.join(location,"saved_optim/") .. GENERATED FROM PYTHON SOURCE LINES 64-69 Open and visualise images before registration. The source and target are 'C' shapes. The source is a 'C' shape that is deformed. The target is a 'C' shape that was cut in half changing its topology and a point was added. The goal is to register the source to the target by deforming the 'C' shape and accounting the cut and the point as intensity additions. .. GENERATED FROM PYTHON SOURCE LINES 69-88 .. code-block:: Python print("ROOT_DIRECTORY : ",ROOT_DIRECTORY) source_name,target_name = 'm0t', 'm1c' # other suggestion of images to try # source_name,target_name = '17','20' # easy, only deformation # source_name,target_name = '08','m1c' # hard, big deformation ! size = (300,300) S = tb.reg_open(source_name,size = size) T = tb.reg_open(target_name,size = size) fig, ax = plt.subplots(1,3,figsize=(10,5)) ax[0].imshow(S[0,0],**DLT_KW_IMAGE) ax[0].set_title('source') ax[1].imshow(T[0,0],**DLT_KW_IMAGE) ax[1].set_title('target') ax[2].imshow(tb.imCmp(S,T,'seg'),origin='lower') ax[2].set_title('superposition of S and T') plt.show() .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_001.png :alt: source, target, superposition of S and T :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ROOT_DIRECTORY : /home/runner/work/Demeter_metamorphosis/Demeter_metamorphosis .. GENERATED FROM PYTHON SOURCE LINES 89-99 Before choosing the optimisation method, we need to define a reproducing kernel and choose a good sigma. The simpler reproducing kernel is the Gaussian kernel. To choose the sigma, we can use the helper functions. `get_sigma_from_img_ratio` and `plot_kernel_on_image`. The first one will compute a good sigma to match the level of details desired. A big sigma will produce a smoother deformation field that will register better big structures. A smaller sigma will register better small details. The subdivisions is basically in how many parts we want to divide the image to get the size of wanted details. The second function will plot the kernel on the image to help us validate our choice of sigma. .. GENERATED FROM PYTHON SOURCE LINES 99-109 .. code-block:: Python image_subdivisions = 10 sigma = rk.get_sigma_from_img_ratio(T.shape,subdiv = image_subdivisions) kernelOperator = rk.GaussianRKHS(sigma,kernel_reach=4) rk.plot_kernel_on_image(kernelOperator,image= T,subdiv=image_subdivisions) plt.show() .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_002.png :alt: sigma = (32.189490394340204, 32.189490394340204), subdiv = 10 :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none kernel shape: torch.Size([1, 257, 257]) kernel shape: 257 257 /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/torch/functional.py:513: UserWarning: torch.meshgrid: in an upcoming release, it will be required to pass the indexing argument. (Triggered internally at ../aten/src/ATen/native/TensorShape.cpp:3609.) return _VF.meshgrid(tensors, **kwargs) # type: ignore[attr-defined] x, y torch.Size([257, 257]) torch.Size([257, 257]) .. GENERATED FROM PYTHON SOURCE LINES 110-111 Perform a first Metamorphosis registration .. GENERATED FROM PYTHON SOURCE LINES 111-137 .. code-block:: Python if torch.cuda.is_available(): device = 'cuda:0' else: device = 'cpu' S = S.to(device) T = T.to(device) dx_convention = 'square' # dx_convention = 'pixel' rho = 0.05 # # data_cost = mt.Ssd_normalized(T) data_cost = mt.Ssd(T) mr = mt.metamorphosis(S,T,0, rho, cost_cst=.001, # If the end result is far from the target, try decreasing the cost constant (reduce regularisation) kernelOperator=kernelOperator, integration_steps=10, # If the deformation is big or complex, try increasing the number of integration steps n_iter=15, # If the optimisation did not converge, try increasing the number of iterations grad_coef=1, # if the optimisation diverged, try decreasing the gradient coefficient dx_convention=dx_convention, data_term=data_cost, hamiltonian_integration=True # Set to true if you want to have control over the intermediate steps of the optimisation ) # mr.save(f'round_to_mot_rho{rho}',light_save = True) .. rst-class:: sphx-glr-script-out .. code-block:: none Progress: [#---------] 13.33% (Ssd : ,802.5087). Progress: [##--------] 20.00% (Ssd : ,163.6530). Progress: [###-------] 26.67% (Ssd : , 89.7214). Progress: [###-------] 33.33% (Ssd : , 71.7436). Progress: [####------] 40.00% (Ssd : , 64.0281). Progress: [#####-----] 46.67% (Ssd : , 56.8949). Progress: [#####-----] 53.33% (Ssd : , 52.6968). Progress: [######----] 60.00% (Ssd : , 49.6450). Progress: [#######---] 66.67% (Ssd : , 46.6622). Progress: [#######---] 73.33% (Ssd : , 44.5396). Progress: [########--] 80.00% (Ssd : , 42.8393). Progress: [#########-] 86.67% (Ssd : , 41.6330). Progress: [#########-] 93.33% (Ssd : , 40.6257). Progress: [##########] 100.00% Done... (Ssd : , 39.3657). Computation of forward done in 0:01:53s and 0.855cents s Computation of metamorphosis done in 0:01:53s and 0.855cents s .. GENERATED FROM PYTHON SOURCE LINES 138-143 .. code-block:: Python _, fig_ax = mr.plot() fig_cmp = fig_ax[0] fig_def = mr.plot_deform() mr.mp.plot() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_003.png :alt: Lambda = 0.001 rho = 0.05 :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_004.png :alt: source, target, Integrated source image, comparaison deformed image with target :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_005.png :alt: diffeo = tensor(True) :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_006.png :alt: t = 0.0, diffeo = tensor(True), t = 0.2, diffeo = tensor(True), t = 0.4, diffeo = tensor(True), t = 0.7, diffeo = tensor(True), t = 1.0, diffeo = tensor(True) :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_006.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-0.24801728129386902..1.3701016902923584]. (
, array([[, , ], [, , ], [, , ], [, , ], [, , ]], dtype=object)) .. GENERATED FROM PYTHON SOURCE LINES 144-150 #################################################################### We will test different values of rho to see how the registration behaves with different values of rho. To save you time I already computed the optimisation for the values of rho in the files listed in list_optim. Feel free to try yourselves If you want to recompute them by setting recompute to True. The number of rho to test is set by n_plot. .. GENERATED FROM PYTHON SOURCE LINES 150-196 .. code-block:: Python list_optim = [ "2D_23_01_2025_simpleToyExample_rho_0.00_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.11_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.22_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.33_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.44_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.56_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.67_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.78_000.pk1", "2D_23_01_2025_simpleToyExample_rho_0.89_000.pk1", "2D_23_01_2025_simpleToyExample_rho_1.00_000.pk1", ] recompute = False n_plot = 10 rho_list = torch.linspace(0,1,n_plot) fig,ax = plt.subplots(2,n_plot,figsize=(20,5)) for i,rho in enumerate(rho_list): print(f'\nrho = {rho}, {i+1}/{n_plot}') if recompute: mr = mt.metamorphosis(S,T,0, rho, cost_cst=.001, kernelOperator=kernelOperator, integration_steps=10, n_iter=30, grad_coef=.1, dx_convention=dx_convention, data_term=data_cost, hamiltonian_integration=True ) mr.save(f'simpleToyExample_rho_{rho:.2f}',light_save = True) else: mr = mt.load_optimize_geodesicShooting(list_optim[i],path =EXPL_SAVE_FOLDER) # mr.plot_cost() ax[0,i].set_title(f'rho = {rho:.2f}') ax[0,i].imshow(mr.mp.image[0,0].detach().cpu(),**DLT_KW_IMAGE) deform = mr.mp.get_deformator() img_deform = tb.imgDeform(S.cpu(),deform,dx_convention=dx_convention) ax[1,i].imshow(img_deform[0,0].detach().cpu(),**DLT_KW_IMAGE) plt.show() # sphinx_gallery_thumbnail_number = 4 .. image-sg:: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_007.png :alt: rho = 0.00, rho = 0.11, rho = 0.22, rho = 0.33, rho = 0.44, rho = 0.56, rho = 0.67, rho = 0.78, rho = 0.89, rho = 1.00 :srcset: /auto_examples/1_registration/images/sphx_glr_simpleToyExample_metamorphosis_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none rho = 0.0, 1/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.00_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.0, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.1111111119389534, 2/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.11_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.1111111119389534, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.2222222238779068, 3/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.22_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.2222222238779068, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.3333333432674408, 4/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.33_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.3333333432674408, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.4444444477558136, 5/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.44_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.4444444477558136, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.5555555820465088, 6/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.56_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.5555555820465088, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.6666666865348816, 7/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.67_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.6666666865348816, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.7777777910232544, 8/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.78_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.7777777910232544, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 0.8888888955116272, 9/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_0.89_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =0.8888888955116272, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) rho = 1.0, 10/10 DT: None New optimiser loaded (2D_23_01_2025_simpleToyExample_rho_1.00_000.pk1) : Metamorphosis_Shooting(cost_parameters : {, rho =1.0, lambda =0.001 }, geodesic integrator : Metamorphosis_integrator( (kernelOperator): GaussianRKHS,2D filter :fft_filter, sigma :(32.189490394340204, 32.189490394340204) kernel_size :(1, 225, 225) kernel_reach :7 normalized :True ) integration method : _step_full_semiLagrangian optimisation method : LBFGS_torch # geodesic steps =10 ) .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 58.686 seconds) .. _sphx_glr_download_auto_examples_1_registration_simpleToyExample_metamorphosis.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: simpleToyExample_metamorphosis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: simpleToyExample_metamorphosis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: simpleToyExample_metamorphosis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_