如何在嵌入3D的表面上绘制测地曲线? [英] How to plot geodesic curves on a surface embedded in 3D?
问题描述
我想到了此视频或此
I have in mind this video, or this simulation, and I would like to reproduce the geodesic lines on some sort of surface in 3D, given by a function f(x,y), from some starting point.
中点方法似乎在计算上和代码上都很密集，因此我想问问是否有一种方法可以根据法线向量在不同点生成近似测地曲线.每个点都有一个与之相关的切向量空间，因此，似乎知道法线并不能确定向前移动曲线的特定方向.
The midpoint method seems computationally and code intense, and so I'd like to ask if there is a way to generate an approximate geodesic curve based on the normal vector to the surface at different points. Each point has a tangent vector space associated with it, and therefore, it seems like knowing the normal vector does not determine a specific direction to move forward the curve.
我曾尝试过与Geogebra一起工作，但我意识到可能有必要转向其他软件平台，例如Python(或Poser?)，Matlab或其他平台.
I have tried working with Geogebra, but I realize that it may be necessary to shift to other software platforms, such as Python (or Poser?), Matlab, or others.
这个想法有可能吗，我能获得一些有关如何实现它的想法吗?
Is this idea possible, and can I get some ideas as to how to implement it?
如果它提供了有关如何回答该问题的一些想法，以前有一个答案(现在已消除了不幸)，建议使用函数形式z = F(x，y)的地形的中点方法，从端点之间的直线，分成短段[我假定XY平面(?)上的直线]，然后抬起[我假定XY平面(?)上线段之间的节点].接下来，它建议找到一个中点"(我猜想连接曲面上每对连续的投影点对的线段的中点(?))，然后投影它" [我想这些中点中的每个都靠近但不完全在[?(?)]]在表面上正交(在法线方向上)，使用等式Z + t = F(X + t Fx，Y + t Fy)[我想这是一个旨在为零的点积...
In case it provides some ideas as to how to answer the question, there previously was an answer (now unfortunatley erased) suggesting the midpoint method for a terrain with the functional form z = F(x,y), starting with the straight line between the endpoints, splitting in short segments [I presume the straight line on the XY plane (?)], and lifting [I presume the nodes between segments on the XY plane (?)] on the surface. Next it suggested finding "a midpoint" [I guess a midpoint of the segments joining each consecutive pairs of projected points on the surface (?)], and projecting "it" [I guess each one of these midpoints close, but not quite on the surface(?)] orthogonally on the surface (in the direction of the normal), using the equation Z + t = F(X + t Fx, Y + t Fy) [I guess this is a dot product meant to be zero...
(?)]，其中(X，Y，Z)是中点的坐标，Fx，Fy是F的偏导数，t是未知数[这是我理解这一点的主要问题...我是什么?一旦找到它就应该处理这个t吗?是否像(X + t，Y + t，Z + t)一样将其添加到(X，Y，Z)的每个坐标上?接着?].这是t中的非线性方程，可通过牛顿迭代解决.
(?)], where (X,Y,Z) are the coordinates of the midpoint, Fx, Fy the partial derivatives of F, and t the unknown [that is my main issue understanding this... What am I supposed to do with this t once I find it? Add it to each coordinate of (X,Y,Z) as in (X+t, Y+t, Z+t)? And then?]. This is a non-linear equation in t, solved via Newton's iterations.
作为更新/书签，Alvise Vianello在此页面页面在GitHub 上.非常感谢你！
As an update / bookmark, Alvise Vianello has kindly posted a Python computer simulation of geodesic lines inspired on this page on GitHub. Thank you very much!
推荐答案
我有一种方法应该适用于任意3D表面，即使该表面上有孔或有噪声.目前它的运行速度很慢，但它似乎可以正常工作，并且可能会为您提供一些有关如何执行此操作的想法.
I have an approach that should be applicable to an arbitrary 3D surface, even when that surface has holes in it or is noisy. It's pretty slow right now, but it seems to work and may give you some ideas for how to do this.
基本前提是微分几何，其前提是:
The basic premise is a differential geometric one and is to:
1.)生成一个代表您的表面的点集
1.) Generate a pointset representing your surface
2.)从该点集生成一个k最近邻图(我在这里也将尺寸间的距离归一化，因为我觉得它更准确地捕获了邻居"的概念)
2.) Generate a k nearest neighbors proximity graph from this pointset (I also normalized distances across dimensions here as I felt it captured the notion of "neighbors" more accurately)
3.)通过使用该点及其邻居作为矩阵的列来计算与该邻近图中的每个节点相关的切线空间，然后在该矩阵上执行SVD. SVD之后，左奇异矢量为我的切线空间提供了新的基础(前两个列矢量是我的平面矢量，第三个垂直于该平面)
3.) Calculate the tangent spaces associated with each node in this proximity graph by using the point and its neighbors as columns of a matrix that I then perform SVD on. After SVD, the left singular vectors give me a new basis for my tangent space (the first two column vectors are my plane vectors, and the third is normal to the plane)
4.)使用dijkstra的算法从该邻近图上的开始节点移动到结束节点，但不是使用欧几里得距离作为边缘权重，而是使用通过切线空间平行传输的向量之间的距离.
4.) Use dijkstra's algorithm to move from a starting node to an ending node on this proximity graph, but instead of using euclidean distance as edge weights, use the distance between vectors being parallel transported via tangent spaces.
它受到本文的启发(减去所有展开的内容): https://arxiv.org/pdf /1806.09039.pdf
It's inspired by this paper (minus all the unfolding): https://arxiv.org/pdf/1806.09039.pdf
请注意，我留下了一些我正在使用的辅助功能，这些功能可能与您没有直接关系(大部分是平面绘图工作).
Note that I left a few helper functions I was using in that probably aren't relevant to you directly (the plane plotting stuff mostly).
您要查看的函数是get_knn，build_proxy_graph，generate_tangent_spaces和geodesic_single_path_dijkstra.
The functions you'll want to look at are get_knn, build_proxy_graph, generate_tangent_spaces, and geodesic_single_path_dijkstra.
实施方式也可能会得到改善.
The implementation could also probably be improved.
代码如下:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mayavi import mlab
from sklearn.neighbors import NearestNeighbors
from scipy.linalg import svd
import networkx as nx
import heapq
from collections import defaultdict
def surface_squares(x_min, x_max, y_min, y_max, steps):
x = np.linspace(x_min, x_max, steps)
y = np.linspace(y_min, y_max, steps)
xx, yy = np.meshgrid(x, y)
zz = xx**2 + yy**2
return xx, yy, zz
def get_meshgrid_ax(x, y, z):
# fig = plt.figure()
# ax = fig.gca(projection='3d')
# ax.plot_surface(X=x, Y=y, Z=z)
# return ax
fig = mlab.figure()
su = mlab.surf(x.T, y.T, z.T, warp_scale=0.1)
def get_knn(flattened_points, num_neighbors):
# need the +1 because each point is its own nearest neighbor
knn = NearestNeighbors(num_neighbors+1)
# normalize flattened points when finding neighbors
neighbor_flattened = (flattened_points - np.min(flattened_points, axis=0)) / (np.max(flattened_points, axis=0) - np.min(flattened_points, axis=0))
knn.fit(neighbor_flattened)
dist, indices = knn.kneighbors(neighbor_flattened)
return dist, indices
def rotmatrix(axis, costheta):
""" Calculate rotation matrix
Arguments:
- `axis` : Rotation axis
- `costheta` : Rotation angle
"""
x, y, z = axis
c = costheta
s = np.sqrt(1-c*c)
C = 1-c
return np.matrix([[x*x*C+c, x*y*C-z*s, x*z*C+y*s],
[y*x*C+z*s, y*y*C+c, y*z*C-x*s],
[z*x*C-y*s, z*y*C+x*s, z*z*C+c]])
def plane(Lx, Ly, Nx, Ny, n, d):
""" Calculate points of a generic plane
Arguments:
- `Lx` : Plane Length first direction
- `Ly` : Plane Length second direction
- `Nx` : Number of points, first direction
- `Ny` : Number of points, second direction
- `n` : Plane orientation, normal vector
- `d` : distance from the origin
"""
x = np.linspace(-Lx/2, Lx/2, Nx)
y = np.linspace(-Ly/2, Ly/2, Ny)
# Create the mesh grid, of a XY plane sitting on the orgin
X, Y = np.meshgrid(x, y)
Z = np.zeros([Nx, Ny])
n0 = np.array([0, 0, 1])
# Rotate plane to the given normal vector
if any(n0 != n):
costheta = np.dot(n0, n)/(np.linalg.norm(n0)*np.linalg.norm(n))
axis = np.cross(n0, n)/np.linalg.norm(np.cross(n0, n))
rotMatrix = rotmatrix(axis, costheta)
XYZ = np.vstack([X.flatten(), Y.flatten(), Z.flatten()])
X, Y, Z = np.array(rotMatrix*XYZ).reshape(3, Nx, Ny)
eps = 0.000000001
dVec = d #abs((n/np.linalg.norm(n)))*d#np.array([abs(n[i])/np.linalg.norm(n)*val if abs(n[i]) > eps else val for i, val in enumerate(d)]) #
X, Y, Z = X+dVec[0], Y+dVec[1], Z+dVec[2]
return X, Y, Z
def build_proxy_graph(proxy_n_dist, proxy_n_indices):
G = nx.Graph()
for distance_list, neighbor_list in zip(proxy_n_dist, proxy_n_indices):
# first element is always point
current_node = neighbor_list[0]
neighbor_list = neighbor_list[1:]
distance_list = distance_list[1:]
for neighbor, dist in zip(neighbor_list, distance_list):
G.add_edge(current_node, neighbor, weight=dist)
return G
def get_plane_points(normal_vec, initial_point, min_range=-10, max_range=10, steps=1000):
steps_for_plane = np.linspace(min_range, max_range, steps)
xx, yy = np.meshgrid(steps_for_plane, steps_for_plane)
d = -initial_point.dot(normal_vec)
eps = 0.000000001
if abs(normal_vec[2]) < eps and abs(normal_vec[1]) > eps:
zz = (-xx*normal_vec[2] - yy*normal_vec[0] - d)/normal_vec[1]
else:
zz = (-xx*normal_vec[0] - yy*normal_vec[1] - d)/normal_vec[2]
return xx, yy, zz
# def plot_tangent_plane_at_point(pointset, flattened_points, node, normal_vec):
# ax = get_meshgrid_ax(x=pointset[:, :, 0], y=pointset[:, :, 1], z=pointset[:, :, 2])
# node_loc = flattened_points[node]
# print("Node loc: {}".format(node_loc))
# xx, yy, zz = plane(10, 10, 500, 500, normal_vec, node_loc)
# # xx, yy, zz = get_plane_points(normal_vec, node_loc)
# print("Normal Vec: {}".format(normal_vec))
# ax.plot_surface(X=xx, Y=yy, Z=zz)
# ax.plot([node_loc[0]], [node_loc[1]], [node_loc[2]], markerfacecolor='k', markeredgecolor='k', marker='o', markersize=10)
# plt.show()
def generate_tangent_spaces(proxy_graph, flattened_points):
# This depth should gaurantee at least 16 neighbors
tangent_spaces = {}
for node in proxy_graph.nodes():
neighbors = list(nx.neighbors(proxy_graph, node))
node_point = flattened_points[node]
zero_mean_mat = np.zeros((len(neighbors)+1, len(node_point)))
for i, neighbor in enumerate(neighbors):
zero_mean_mat[i] = flattened_points[neighbor]
zero_mean_mat[-1] = node_point
zero_mean_mat = zero_mean_mat - np.mean(zero_mean_mat, axis=0)
u, s, v = svd(zero_mean_mat.T)
# smat = np.zeros(u.shape[0], v.shape[0])
# smat[:s.shape[0], :s.shape[0]] = np.diag(s)
tangent_spaces[node] = u
return tangent_spaces
def geodesic_single_path_dijkstra(flattened_points, proximity_graph, tangent_frames, start, end):
# short circuit
if start == end:
return []
# Create min priority queue
minheap = []
pred = {}
dist = defaultdict(lambda: 1.0e+100)
# for i, point in enumerate(flattened_points):
R = {}
t_dist = {}
geo_dist = {}
R[start] = np.eye(3)
t_dist[start] = np.ones((3,))
dist[start] = 0
start_vector = flattened_points[start]
for neighbor in nx.neighbors(proxy_graph, start):
pred[neighbor] = start
dist[neighbor] = np.linalg.norm(start_vector - flattened_points[neighbor])
heapq.heappush(minheap, (dist[neighbor], neighbor))
while minheap:
r_dist, r_ind = heapq.heappop(minheap)
if r_ind == end:
break
q_ind = pred[r_ind]
u, s, v = svd(tangent_frames[q_ind].T*tangent_frames[r_ind])
R[r_ind] = np.dot(R[q_ind], u * v.T)
t_dist[r_ind] = t_dist[q_ind]+np.dot(R[q_ind], tangent_frames[q_ind].T * (r_dist - dist[q_ind]))
geo_dist[r_ind] = np.linalg.norm(t_dist[r_ind])
for neighbor in nx.neighbors(proxy_graph, r_ind):
temp_dist = dist[r_ind] + np.linalg.norm(flattened_points[neighbor] - flattened_points[r_ind])
if temp_dist < dist[neighbor]:
dist[neighbor] = temp_dist
pred[neighbor] = r_ind
heapq.heappush(minheap, (dist[neighbor], neighbor))
# found ending index, now loop through preds for path
current_ind = end
node_path = [end]
while current_ind != start:
node_path.append(pred[current_ind])
current_ind = pred[current_ind]
return node_path
def plot_path_on_surface(pointset, flattened_points, path):
# ax = get_meshgrid_ax(x=pointset[:, :, 0], y=pointset[:, :, 1], z=pointset[:, :, 2])
# ax.plot(points_in_path[:, 0], points_in_path[:, 1], points_in_path[:, 2], linewidth=10.0)
# plt.show()
get_meshgrid_ax(x=pointset[:, :, 0], y=pointset[:, :, 1], z=pointset[:, :, 2])
points_in_path = flattened_points[path]
mlab.plot3d(points_in_path[:, 0], points_in_path[:, 1], points_in_path[:, 2] *.1)
mlab.show()
"""
True geodesic of graph.
Build proximity graph
Find tangent space using geodisic neighborhood at each point in graph
Parallel transport vectors between tangent space points
Use this as your distance metric
Dijkstra's Algorithm
"""
if __name__ == "__main__":
x, y, z = surface_squares(-5, 5, -5, 5, 500)
# plot_meshgrid(x, y, z)
pointset = np.stack([x, y, z], axis=2)
proxy_graph_num_neighbors = 16
flattened_points = pointset.reshape(pointset.shape[0]*pointset.shape[1], pointset.shape[2])
flattened_points = flattened_points
proxy_n_dist, proxy_n_indices = get_knn(flattened_points, proxy_graph_num_neighbors)
# Generate a proximity graph using proxy_graph_num_neighbors
# Nodes = number of points, max # of edges = number of points * num_neighbors
proxy_graph = build_proxy_graph(proxy_n_dist, proxy_n_indices)
# Now, using the geodesic_num_neighbors, get geodesic neighborshood for tangent space construction
tangent_spaces = generate_tangent_spaces(proxy_graph, flattened_points)
node_to_use = 2968
# 3rd vector of tangent space is normal to plane
# plot_tangent_plane_at_point(pointset, flattened_points, node_to_use, tangent_spaces[node_to_use][:, 2])
path = geodesic_single_path_dijkstra(flattened_points, proxy_graph, tangent_spaces, 250, 249750)
plot_path_on_surface(pointset, flattened_points, path)
请注意，我安装并设置了mayavi以获得不错的输出图像(matplotlib没有真正的3d渲染，因此其图很烂).但是，如果您要使用matplotlib代码，则没有做.如果这样做，只需在路径绘图仪中将缩放比例缩小.1，然后取消注释绘图代码即可.无论如何，这是z = x ^ 2 + y ^ 2的示例图像.白线是测地路径:
Note that I installed and set up mayavi to get a decent output image (matplotlib doesn't have real 3d rendering and consequently, its plots suck). I did however leave the matplotlib code in if you want to use it. If you do, just remove the scaling by .1 in the path plotter and uncomment the plotting code. Anyways, here's an example image for z=x^2+y^2. The white line is the geodesic path:
您也可以很容易地对此进行调整，以返回dijkstra算法中节点之间的所有成对测地距离(请参阅本文的附录，以了解进行此操作所需的细微修改).然后，您可以在表面上绘制任何线条.
You could also fairly easily adjust this to return all the pairwise geodesic distances between nodes from dijkstra's algorithm (look in the appendix of the paper to see the minor modifications you'll need to do this). Then you could draw whatever lines you want on your surface.
这篇关于如何在嵌入3D的表面上绘制测地曲线?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
