How to efficiently determine the intersection of a vertical cutting plane with a mesh

Question

I have a list of vertical cut planes, and I have a polygonal mesh ( it's a 2D+0.5D mesh, something like a 2D mesh with an extra dimension, $z$ attached to each point). One can assume that the mesh contains vertices $V$ , edges $E$ and face $F$.

The mesh is a real life 3D terrain.

Now, given a list of vertical cut planes, I want to find all the intersection points ( to form a polyline) between the cut planes and the polygonal mesh, in order to build a vertical view of the terrain for each of the cut plane.

Currently, my approach is that I use the Aabb tree to test each of the mesh triangles against the cut section, and then compute the intersection points subsequently for the cut section and the triangle that are within each other's Aabb tree range. But this approach doesn't make use of the mesh structure of the terrain, and I think we can significantly speed up the computation if we can make use of the mesh structure.

Is there a more efficient algorithm I can use for this purpose?

Answer 1

You have an 2D mesh of height values. I'm going to assume you create the faces of the surface for the $i$,$j$th square using two triangles of $\left\langle\left(i,j\right),\left(i+1,j\right),\left(i+1,j+1\right)\right\rangle$ and $\left\langle\left(i,j\right),\left(i,j+1\right),\left(i+1,j+1\right)\right\rangle$.

Your vertical cut plain can be defined as a 2D line called $l$. $l$ can be represented by using 2 known points on the edge of the mesh or by using an 2D line equation $i = mj+c$, $l=\left\langle{m,c}\right\rangle$.

The intersection points are the intersection with any of the edges of triangles along that 2D line $l$.

All the points heights are worked out by using linear interpolation between the heights given by the two array coordinated connected by the triangle edge being intersected.

Distance along the cut plane is given by Pythagorean theorem.

#!/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)