Python scipy.ndimage.filters.laplace() Examples

def drlse_edge(phi_0, g, lmda, mu, alfa, epsilon, timestep, iters, potentialFunction):  # Updated Level Set Function
    """

    :param phi_0: level set function to be updated by level set evolution
    :param g: edge indicator function
    :param lmda: weight of the weighted length term
    :param mu: weight of distance regularization term
    :param alfa: weight of the weighted area term
    :param epsilon: width of Dirac Delta function
    :param timestep: time step
    :param iters: number of iterations
    :param potentialFunction: choice of potential function in distance regularization term.
%              As mentioned in the above paper, two choices are provided: potentialFunction='single-well' or
%              potentialFunction='double-well', which correspond to the potential functions p1 (single-well)
%              and p2 (double-well), respectively.
    """
    phi = phi_0.copy()
    [vy, vx] = np.gradient(g)
    for k in range(iters):
        phi = NeumannBoundCond(phi)
        [phi_y, phi_x] = np.gradient(phi)
        s = np.sqrt(np.square(phi_x) + np.square(phi_y))
        smallNumber = 1e-10
        Nx = phi_x / (s + smallNumber)  # add a small positive number to avoid division by zero
        Ny = phi_y / (s + smallNumber)
        curvature = div(Nx, Ny)
        if potentialFunction == 'single-well':
            distRegTerm = filters.laplace(phi, mode='wrap') - curvature  # compute distance regularization term in equation (13) with the single-well potential p1.
        elif potentialFunction == 'double-well':
            distRegTerm = distReg_p2(phi)  # compute the distance regularization term in eqaution (13) with the double-well potential p2.
        else:
            print('Error: Wrong choice of potential function. Please input the string "single-well" or "double-well" in the drlse_edge function.')
        diracPhi = Dirac(phi, epsilon)
        areaTerm = diracPhi * g  # balloon/pressure force
        edgeTerm = diracPhi * (vx * Nx + vy * Ny) + diracPhi * g * curvature
        phi = phi + timestep * (mu * distRegTerm + lmda * edgeTerm + alfa * areaTerm)
    return phi 

def distReg_p2(phi):
    """
        compute the distance regularization term with the double-well potential p2 in equation (16)
    """
    [phi_y, phi_x] = np.gradient(phi)
    s = np.sqrt(np.square(phi_x) + np.square(phi_y))
    a = (s >= 0) & (s <= 1)
    b = (s > 1)
    ps = a * np.sin(2 * np.pi * s) / (2 * np.pi) + b * (s - 1)  # compute first order derivative of the double-well potential p2 in equation (16)
    dps = ((ps != 0) * ps + (ps == 0)) / ((s != 0) * s + (s == 0))  # compute d_p(s)=p'(s)/s in equation (10). As s-->0, we have d_p(s)-->1 according to equation (18)
    return div(dps * phi_x - phi_x, dps * phi_y - phi_y) + filters.laplace(phi, mode='wrap') 

def measure_sharpness(img):
    """Measures the sharpeness of an image using the variance of the laplacian

    img: PIL.Image

    Returns the variance of the laplacian. Higher values mean a sharper image
    """
    img_gray = np.array(img.convert('L'), dtype=np.float32)
    return np.var(laplace(img_gray)) 

def test_laplace_comprehensions():
    np.random.seed(0)

    a = np.random.random((3, 12, 14))
    d = da.from_array(a, chunks=(3, 6, 7))

    l2s = [da_ndf.laplace(d[i]) for i in range(len(d))]
    l2c = [da_ndf.laplace(d[i])[None] for i in range(len(d))]

    dau.assert_eq(np.stack(l2s), da.stack(l2s))
    dau.assert_eq(np.concatenate(l2c), da.concatenate(l2c)) 

def test_laplace_compare():
    s = (10, 11, 12)
    a = np.arange(float(np.prod(s))).reshape(s)
    d = da.from_array(a, chunks=(5, 5, 6))

    dau.assert_eq(
        sp_ndf.laplace(a),
        da_ndf.laplace(d)
    ) 

def gvf(f, mu=0.05, iterations=30, anisotropic=False,
        ignore_second_term=False):
    # Gradient vector flow
    # Translated from https://github.com/smistad/3D-Gradient-Vector-Flow-for-Matlab
    f = (f - f.min()) / (f.max() - f.min())
    f = enforce_mirror_boundary(
        f)  # Enforce the mirror conditions on the boundary

    dx, dy, dz = np.gradient(f)  # Initialse with normal gradients
    '''
    Initialise the GVF vectors following S3 in
    Yu, Zeyun, and Chandrajit Bajaj. 
    "A segmentation-free approach for skeletonization of gray-scale images via anisotropic vector diffusion." 
    CVPR, 2004. CVPR 2004. 
    It only uses one of the surronding neighbours with the lowest intensity
    '''
    magsq = dx**2 + dy**2 + dz**2

    # Set up the initial vector field
    u = dx.copy()
    v = dy.copy()
    w = dz.copy()

    for i in tqdm(range(iterations)):
        # The boundary might not matter here
        # u = enforce_mirror_boundary(u)
        # v = enforce_mirror_boundary(v)
        # w = enforce_mirror_boundary(w)

        # Update the vector field
        if anisotropic:
            G = g_all(u, v, w)
            u += mu / 6. * div(np.sum(G * d(u), axis=0))
            v += mu / 6. * div(np.sum(G * d(v), axis=0))
            w += mu / 6. * div(np.sum(G * d(w), axis=0))
        else:
            u += mu * 6 * laplace(u)
            v += mu * 6 * laplace(v)
            w += mu * 6 * laplace(w)

        if not ignore_second_term:
            u -= (u - dx) * magsq
            v -= (v - dy) * magsq
            w -= (w - dz) * magsq

    return np.stack((u, v, w), axis=0)