Weighted phase Laplacian formula

Question

We have the following equation (weighted phase Laplacian):

\$c_{i,j} = U(i,j)\Delta^{^x}_{i,j} - U(i-l, j)\Delta^{x}_{i-1,j} + V(i,j)\Delta^{y}_{i,j} - V(i,j-1)\Delta^{y}_{i,j-1}\$

Where

\$U(i,j) = min(w^2_{i+1}, w^2_{i+1})V(i,j) = min(w^2_{j+1}, w^2_{j+1})\$

with \$U, V, \Delta \$ (N,N) matrices

I coded the following code:

def compute_Weights(w):
    U = np.zeros_like(w)
    V = np.zeros_like(w)
    U = np.minimum(w[1:, :]**2, w[:-1, :]**2)
    V = np.minimum(w[:, 1:]**2, w[:, :-1]**2)

    """add the last ones rows and cols
    """
    U = np.vstack((U, U[-1]))
    V = np.vstack((V.T, V.T[-1])).transpose()

    return U, V

def Gradient(p1, p2):
    r = p1 - p2
    if r > np.pi:
        return r - 2 * np.pi
    if r < -np.pi:
        return r + 2 * np.pi
    return r
"""
    phase is not flattened
    wts are flattened
"""
def ComputePhaseLaplacian(phase, U, V):
    rows, cols = phase.shape
    phi = phase.flatten()
    rho = np.zeros_like(phi)

    for j in range(rows):
        for i in range(cols):
            k = j * rows + i
            k1 = k + 1 if i < cols - 1 else k - 1
            k2 = k - 1 if i > 0 else k + 1
            k3 = k + cols if j < rows - 1 else k - cols
            k4 = k - cols if j > 0 else k + cols


            w1 = U[k]
            w2 = U[k-1] if i > 0 else U[k]
            w3 = V[k]
            w4 = V[k-cols] if j > 0 else V[k]

            rho[k] = w1 * Gradient(phi[k], phi[k1]) \
                + w2 * Gradient(phi[k], phi[k2]) \
                + w3 * Gradient(phi[k], phi[k3]) \
                + w4 * Gradient(phi[k], phi[k4])

    return rho

I want to know if there is a "pythonic way" to write this. I'm a C/C++ programmer so this was my perspective. Maybe something to avoid to work with flattened arrays, and to work with 2D arrays.

I.e:

phase = np.load('phase.npy') #phase 512x512 float array
U,V = ComputeWeigths(W) #W 512x512 float array with values in (0,1)
uphase = ComputePhaseLaplacian(phase, U, V) #uphase 512x512 float array

Answer 1

A few points. In your code, some of these render others irrelevant, but they may be relevant in other code you write:

1. Try to follow PEP8
2. It should be possible to vectorize ComputePhaseLaplacian.
3. It is better to loop over rows and columns directly rather than looping over indexes. You can use enumerate to keep track of the rows and columns.
4. You shouldn't flatten, you can index in two dimensions.
5. Instead of vstack, you should use pad in this case.
6. You create U and V, then immediately overwrite them. This is redundant.
7. You should use ravel instead of flatten since it doesn't make a copy.