# Weighted phase Laplacian formula

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

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]) \

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

• very strange use of docstrings... See PEP257 for the proper usage. – Christopher Pearson Jun 18 '15 at 21:48
• Sorry, I put it just to clarify some of the code here. But thank you it is useful! – FacundoGFlores Jun 19 '15 at 13:00

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.
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.
• 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. – TheBlackCat Jun 18 '15 at 16:48