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}\$


\$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.


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
  • 2
    \$\begingroup\$ very strange use of docstrings... See PEP257 for the proper usage. \$\endgroup\$ Commented Jun 18, 2015 at 21:48
  • \$\begingroup\$ Sorry, I put it just to clarify some of the code here. But thank you it is useful! \$\endgroup\$ Commented Jun 19, 2015 at 13:00

1 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.
  • \$\begingroup\$ Thank you a lot! Have you any idea about vectorizing ComputePhaseLaplacian? \$\endgroup\$ Commented Jun 18, 2015 at 16:45
  • 1
    \$\begingroup\$ I would first generate four copies of phi, each one shifted by one element to handle the k1, k2, k3, and k4 cases. w1 and w3 are just U and V. w2 and w4 are again just a shifted U and V. Gradient can be converted to logic indexing on the resulting arrays. \$\endgroup\$ Commented Jun 18, 2015 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.