如何有效地确定垂直切割面与网格的交点

计算科学 计算几何 网格遍历
2021-12-11 11:36:31

我有一个垂直切割平面列表,并且我有一个多边形网格(它是一个 2D+0.5D 网格,类似于具有额外维度的 2D 网格,附加到每个点)。可以假设网格包含顶点、边和面zVEF

网格是真实的 3D 地形。

现在,给定一个垂直剖切面列表,我想找到剖切面和多边形网格之间的所有交点(形成一条折线),以便为每个剖切面构建地形的垂直视图。

目前,我的方法是使用Aabb 树针对切割部分测试每个网格三角形,然后随后计算切割部分和三角形在彼此的 Aabb 树范围内的交点。但是这种方法并没有利用地形的网格结构,如果我们能够利用网格结构,我认为我们可以显着加快计算速度。

我可以为此目的使用更有效的算法吗?

1个回答

您有一个高度值的二维网格。我将假设您 ,个正方形ij(i,j),(i+1,j),(i+1,j+1)(i,j),(i,j+1),(i+1,j+1)

您的垂直切割平原可以定义为称为的二维线。可以通过使用网格边缘上的 2 个已知点或使用二维线方程来表示。lli=mj+cl=m,c

交点是与沿该 2D 线的三角形的任何边的交点。l

以 2D 线为剖切面的 2D 网格

所有点的高度都是通过在由被相交的三角形边缘连接的两个坐标阵列给出的高度之间使用线性插值来计算的。

沿剖切面的距离由勾股定理给出

#!/usr/bin/env python

import math
from random import choice, uniform

mesh = [[10,5,8],[12,7,5],[4,7,3],[4,6,4],[8,3,1]]
print mesh

iMax = len(mesh)-1
jMax = len(mesh[0])-1

# pick random vertical cut plane
if choice([True,False]):
  p1 = (uniform(0,iMax), choice([0,jMax]))
else:
  p1 = (choice([0,iMax]), uniform(0,jMax))
while True:
  if choice([True,False]):
    p2 = (uniform(0,iMax), choice([0,jMax]))
  else:
    p2 = (choice([0,iMax]), uniform(0,jMax))
  if p1[0] != p2[0] and  p1[1] != p2[1]:
    break

print "P1",p1
print "P2",p2

if p1[1] == p2[1]:
  m = float('Inf') # Not handling this, needs fix
  c = None
else:
  m = float(p2[0]-p1[0])/(p2[1]-p1[1])
  c = p1[0]-m*p1[1]

# i=mj+c formular
print "i={0}j+{1}".format(m,c)

def floor(a):
  return int(math.floor(a))

def ceil(a):
  return int(math.ceil(a))

iL = floor(min(p1[0], p2[0]))
iU = ceil(max(p1[0], p2[0]))
jL = floor(min(p1[1], p2[1]))
jU = ceil(max(p1[1], p2[1]))

print "points to plot unordered, d is distance from p1 (pythagoras), h is height as given by linear interpolation"

def linearInterpolation(h0,h1,r):
  return (h1-h0)*r+h0

def pythag(a,b):
  return math.sqrt(a**2 + b**2)

for i in xrange(iL, iU+1):
  j = (i-c)/m
  if 0 <= j <= jMax:
    j1, j2 = floor(j), ceil(j)
    h = linearInterpolation(mesh[i][j1], mesh[i][j2], j - j1)
    d = pythag(p1[0]-i, p1[1]-j)
    print (d,h)

for j in xrange(jL, jU+1):
  i = m*j+c
  if 0 <= i <= iMax:
    i1, i2 = floor(i), ceil(i)
    h = linearInterpolation(mesh[i1][j], mesh[i2][j], i - i1)
    d = pythag(p1[0]-i, p1[1]-j)
    print (d,h)

for v in xrange(-iU-1, jU):
  i = (c - v)/(1.0-m)
  j = m*i+c
  if 0 <= i <= iMax and 0 <= j <= jMax:
    i1, i2 = floor(i), ceil(i)
    j1, j2 = floor(j), ceil(j)
    r = pythag(i - i1, j - j1) / math.sqrt(2)
    h = linearInterpolation(mesh[i1][j1], mesh[i2][j2], r)
    d = pythag(p1[0]-i, p1[1]-j)
    print (d,h)