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\$ – Christopher Pearson Jun 18 '15 at 21:48
  • \$\begingroup\$ Sorry, I put it just to clarify some of the code here. But thank you it is useful! \$\endgroup\$ – FacundoGFlores Jun 19 '15 at 13:00

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.
| improve this answer | |
  • \$\begingroup\$ Thank you a lot! Have you any idea about vectorizing ComputePhaseLaplacian? \$\endgroup\$ – FacundoGFlores Jun 18 '15 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\$ – TheBlackCat Jun 18 '15 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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