Untitled
unknown
python
a year ago
2.4 kB
6
Indexable
def shape_operator(m,v): pos = m.positions() # Insert code here which computes the shape operator matrix S that can be used to compute # the principal curvatures kmin and kmax and the mean curvature H # Your code here # ------------------ N = m.vertex_normal(v) F = [] U = [] a_tilde = N a_tilde[1] *= -1 projectionOnA = a_tilde - N*(np.dot(a_tilde,N)) a = (projectionOnA)/np.linalg.norm(projectionOnA) b = np.cross(N,a) pi = pos[v] T = np.array([a,b]) neighbours = list(m.circulate_vertex(v,"v")) U = np.zeros((len(neighbours), 3)) F = np.zeros((len(neighbours), 1)) for j, neighbor in enumerate(neighbours): pj = pos[neighbor] uv = T @ (pj - pi) hj = np.dot(N, (pj - pi)) U[j, :] = [uv[0]**2, uv[0]*uv[1], uv[1]**2] F[j] = hj # Solve the least squares problem coeffs, residuals, rank, s = np.linalg.lstsq(U, F, rcond=None) #Get 2D Shape operator S = -1*np.array([[coeffs[0][0], coeffs[1][0]], [coeffs[1][0], coeffs[2][0]]]) #Convert to 3D S = T.T@S@T eigVal,eigVec = np.linalg.eig(S) #The three vectors #1. Eigenvector closest to 0 -> normal #2. Lowest -> min #3. Highest -> max #First we get the normal vector normalIdx = np.argmin(np.abs(eigVal)) normalVector = eigVec[normalIdx] normalValue = eigVal[normalIdx] #Then we check if any of the other vectors are as close to 0 numericalError = max(1**(-1000),normalValue) zeroValueIdx = [] for x in range(len(eigVal)): absEigenValue = abs(eigVal[x]) if(absEigenValue <= numericalError): zeroValueIdx += [x] nrZeroValues = len(zeroValueIdx) print(zeroValueIdx) #Error if(nrZeroValues == 0): return 0/0 if(nrZeroValues == 1): np.delete(eigVal,normalIdx) np.delete(eigVec,normalIdx) kmax = max(eigVal) kmin = min(eigVal) if(nrZeroValues == 2): print("Two 0's - but we're just gonna ignore it instead of figuring out the principal curvature :D") print(eigVal) np.delete(eigVal,normalIdx) np.delete(eigVec,normalIdx) kmax = max(eigVal) kmin = min(eigVal) # Mean curvature H = np.linalg.det(S) return kmin, kmax, H
Editor is loading...
Leave a Comment