-
Notifications
You must be signed in to change notification settings - Fork 2
/
curvature_term.py
77 lines (66 loc) · 2.81 KB
/
curvature_term.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
import numpy as np
import math
class meancurvature:
_eps = 2.2204e-16
def __init__(self, rows, cols):
self.nrows = rows
self.ncols = cols
def get_mean_curvature(self, narrow_band, narrow_band_point_loc,phi):
num_pts = len(narrow_band_point_loc)
K = np.zeros(num_pts)
for i in range(len(narrow_band_point_loc)):
pt = narrow_band_point_loc[i]
nr = pt[0]
nc = pt[1]
# Boundary conditions
if all((nr + 1) >= narrow_band[0]):
nr = self.nrows - 2
if (nr - 1) <= 0:
nr = 2
if all((nc + 1) >= narrow_band[1]):
nc = self.ncols - 2
if (nc - 1) <= 0:
nc = 2
# % derivatives
phi_y = phi[(nr, nc + 1)] - phi[(nr, nc - 1)]
phi_x = phi[(nr + 1, nc)] - phi[(nr - 1, nc)]
phi_yy = phi[(nr, nc + 1)] - 2 * phi[(nr, nc)] + phi[(nr, nc - 1)]
phi_xx = phi[(nr + 1, nc)] - 2 * phi[(nr, nc)] + phi[(nr - 1, nc)]
phi_xy = 0.25 * (-phi[(nr - 1, nc - 1)] - phi[(nr + 1, nc + 1)] + phi[(nr - 1, nc + 1)] + phi[(nr + 1, nc - 1)])
# norm
norm = math.sqrt(phi_x ** 2 + phi_y ** 2)
# K[i]
K[i] = ((phi_x ** 2 * phi_yy + phi_y ** 2 * phi_xx - 2 * phi_x * phi_y * phi_xy) /
(phi_x ** 2 + phi_y ** 2 + self._eps) ** (3 / 2)) * norm
return K
def get_mean_curvature_matrix(self,narrow_band, narrow_band_point_loc,phi):
num_pts = len(narrow_band_point_loc)
base = np.zeros(phi.shape)
K = np.zeros(num_pts)
for i in range(len(narrow_band_point_loc)):
pt = narrow_band_point_loc[i]
nr = pt[0]
nc = pt[1]
# Boundary conditions
if all((nr + 1) >= narrow_band[0]):
nr = self.nrows - 2
if (nr - 1) <= 0:
nr = 2
if all((nc + 1) >= narrow_band[1]):
nc = self.ncols - 2
if (nc - 1) <= 0:
nc = 2
# % derivatives
phi_y = phi[(nr, nc + 1)] - phi[(nr, nc - 1)]
phi_x = phi[(nr + 1, nc)] - phi[(nr - 1, nc)]
phi_yy = phi[(nr, nc + 1)] - 2 * phi[(nr, nc)] + phi[(nr, nc - 1)]
phi_xx = phi[(nr + 1, nc)] - 2 * phi[(nr, nc)] + phi[(nr - 1, nc)]
phi_xy = 0.25 * (
-phi[(nr - 1, nc - 1)] - phi[(nr + 1, nc + 1)] + phi[(nr - 1, nc + 1)] + phi[(nr + 1, nc - 1)])
# norm
norm = math.sqrt(phi_x ** 2 + phi_y ** 2)
# K[i]
K[i] = ((phi_x ** 2 * phi_yy + phi_y ** 2 * phi_xx - 2 * phi_x * phi_y * phi_xy) /
(phi_x ** 2 + phi_y ** 2 + self._eps) ** (3 / 2)) * norm
base[narrow_band] = K
return base