simplest_laplacian_deformation.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm


n = 25

q = n * n
A = np.eye(q)

for i in range(1, n-1):
   for j in range(1, n-1):
       idx = n * i + j
       neigh = [n*(i-1)+j, n*(i+1)+j, n*i+j-1, n*i+j+1]
       for k in neigh:
           A[idx][k] = -1 / len(neigh)

b = np.zeros(q)
b[n*(n//2) + n//2] = 1


z=np.linalg.solve(A, b)
zmat = z.reshape(n, n)
# plt.imshow(zmat)
# plt.colorbar()
# plt.show()


import mpl_toolkits.mplot3d
x = np.arange(n)
y = np.arange(n)
x, y = np.meshgrid(x, y)

fig = plt.figure()
ax = fig.gca(projection='3d')
surf=ax.plot_surface(x, y, zmat, cmap=cm.coolwarm)

plt.show()