evolution.py

import numpy as np
from tqdm import tqdm

from v1diffusion import operators


def parallel_diffusion(u0, T, dt=0.1, beta=0.5, model=0, progress=True):
    """
    Hypoelliptic evolution in PTR2 according to Boscain-2016

    This method takes an input an image lifted in PTR2 and applies the
    hypoelliptic heat operator \Delta_H at time intervals of dt until final time
    T is reached

    u0: image in PTR2
    T: time to run the evolution
    dt: time interval
    beta: weight coefficient
    progress: to either show or not tqdm animated progress bar

    @returns: an image in PTR2
    """

    u1 = u0 #Initialization of u1 in case the loop is not ran

    t = 0
    for _step in tqdm(range(0, int(T/dt)), disable=(not progress)):
        u1 = u1+dt*operators.hypoelliptic_heat_operator(u1, beta, model)
        t += dt

    return np.asarray(u1)


def orthogonal_diffusion(u0, T, dt=0.1, beta=0.5, model=0, progress=True):
    u1 = u0.copy() #Initialization of u1 in case the loop is not ran

    t = 0
    for _step in tqdm(range(0, int(T/dt)), disable=(not progress)):
        u1 = u1+dt*operators.orthogonal_operator(u1, beta, model)
        t += dt

    return u1


def parallel_concentration(u0, T, dt=0.01, beta=0.5, model=0, progress=True):
    """
    Handler function to regress the hypoelliptic evolution with negative dt
    """
    u1 = u0 #Initialization of u1 in case the loop is not ran

    t = 0
    for _step in tqdm(range(0, int(T/dt)), disable=(not progress)):
        u1 = u1-dt*operators.hypoelliptic_heat_operator(u1, beta, model)
        t += dt

    return np.asarray(u1)


def orthogonal_concentration(u0, T, dt=0.01, beta=0.5, model=0, progress=True):
    """
    Handler function to regress the hypoelliptic evolution with negative dt
    """
    u1 = u0 #Initialization of u1 in case the loop is not ran

    t = 0
    for _step in tqdm(range(0, int(T/dt)), disable=(not progress)):
        u1 = u1-dt*operators.orthogonal_operator(u1, beta, model)
        t += dt

    return np.asarray(u1)


def v1_unsharp_filter(u, C=1, T=1, beta=1, model=0):
    """
    SE(2) unsharp filter as described in Ballerin-Grong-25, Section 4.3.

    The image is blurred along X3 (orthogonal to level lines) by running
    the forward orthogonal diffusion for time T.  The blurred version is
    then used as the "mask" in classical unsharp masking:

        output = u + C * (u - u_blurred)

    u: image in PTR2/SE2
    C: sharpening factor (C=0 → no change, C=1 → standard unsharp)
    T: diffusion time (controls blur radius)
    beta: weight coefficient for the X2 component
    """
    u_blurred = orthogonal_diffusion(u, T, beta=beta, model=model)
    return u + C * (u - u_blurred)


def evolve_vc_hypoelliptic(u0, img, T, a0=1, b0=1, a1=1, b1=1, sigma=1,
 epss=0.1, dt=0.1, progress=True):
    r"""
    Hypoelliptic evolution with varying coefficients in PTR2

    This method takes an input an image lifted in PTR2 and applies
    the hypoelliptic heat operator \Delta_H with varying coefficients
    at time intervals of dt until final time T is reached

    u0: image in PTR2
    img: the original 2D greyscale image
    T: time to run the evolution
    a0: starting weight
    b0: starting weight
    a1: starting weight
    a2: starting weight
    sigma: weight
    epss: weight
    dt: time interval
    progress: to either show or not tqdm animated progress bar

    @returns: an image in PTR2
    """

    u1 = u0.copy() #Initialization of u1 in case the loop is not ran

    eps = np.exp(-img**2/sigma)
    t = 0
    for _step in tqdm(range(0, int(T/dt)), disable=(not progress)):
        a = np.zeros(img.shape)
        b = np.zeros(img.shape)
        a_mask = (a0+a1*eps)/(a0+a1)>epss
        b_mask = (b0+b1*eps)/(b0+b1)>epss
        a[a_mask] = a0 + a1*eps[a_mask]
        b[b_mask] = b0 + b1*eps[b_mask]

        u1 = u0+dt*operators.vc_hypoelliptic_heat_operator(u0, a, b)
        u0 = u1.copy()
        t+=dt

    return u1