Start a new topic

Creator scripting - finding the height of a point on a plane directly beneath a known vertex

Original Post by: chrisell Wed Feb 13 22:00:18 2013


I could use some help (couldn't we all?)


I have a cloud of points in creator and I need to know, for each point, IF it's over a polygon and if so, what the height of that polygon is directly underneath each point.


Essentially I need to project a virtual vertical line down from each point and see if it intersects any polys and if it does, what's the z-value of the intersection.


I'm stuck for how to even begin to solve this little problem - any guidance or code snippets would be appreciated.


Original Post by: OutOfMemory Tue Mar 26 17:03:39 2013


Is the polygon on the XY plane (normal pointing Z+) or is the polygon rotated off axis?


Test for vertice Above or Below a polygon:

You could get the Normal(XYZ) of the polygon. With normal position stored then walk each of the vertices in the point cloud data. For each vertice get its position(XYZ) then test if weather the vertex(z) is > than polygon(z) then that vertice is above the polygon.


Determin if vertice is with in the area of a polygon:

Give this a look see might get you moving http://www.blackpawn.com/texts/pointinpoly/

Original Post by: TedLai Mon Apr 29 21:12:54 2013


You can iterate your collection of points and collection of polygons and invoke the function compute_height_above_ground_in_poly, as defined below:


# reference: http://tog.acm.org/resources/GraphicsGems/gemsiv/ptpoly_haines

def is_pt_in_poly_2d(pt, poly):

inside = False

n = mgCountChild(poly)

if n > 1:

vertex = mgGetChildNth(poly, n)

ok, x, y, z = mgGetVtxCoord(vertex)

prev_pt = mgMakeCoord2d(x, y)

is_prev_pt_above_pt = (prev_pt.y >= pt.y)

vertex = mgGetChild(poly)

while vertex:

ok, x, y, z = mgGetVtxCoord(vertex)

this_pt = mgMakeCoord2d(x, y)

is_this_pt_above_pt = (this_pt.y >= pt.y)

if not is_prev_pt_above_pt == is_this_pt_above_pt:

b = ((this_pt.y - pt.y) * (prev_pt.x - this_pt.x) >= (this_pt.x - pt.x) * (prev_pt.y - this_pt.y))

if(b == is_this_pt_above_pt):

inside = not inside

is_prev_pt_above_pt = is_this_pt_above_pt

prev_pt = this_pt

vertex = mgGetNext(vertex)

return inside


def compute_height_above_ground_in_poly(poly, pt):

height = None

if is_pt_in_poly_2d(pt, poly):

vertex = mgGetChild(poly)

ok, x, y, z = mgGetVtxCoord(vertex)

if not ok:

return height

pt1 = mgMakeCoord3d(x, y, z)

vertex = mgGetNext(vertex)

ok, x, y, z = mgGetVtxCoord(vertex)

if not ok:

return height

pt2 = mgMakeCoord3d(x, y, z)

vertex = mgGetNext(vertex)

ok, x, y, z = mgGetVtxCoord(vertex)

if not ok:

return height

pt3 = mgMakeCoord3d(x, y, z)

vec1 = mgCoord3dSubtract(pt2, pt1)

vec2 = mgCoord3dSubtract(pt3, pt1)

normal = mgCoord3dCross(vec1, vec2)

if not normal.z == 0.0:

height = ((x - pt.x) * normal.x + (y - pt.y) * normal.y) / normal.z + z

return height


Since it is inefficient to iterate your collection of polygons over and over again, you may want to build a bounding volume hierarchy for the polygons to cut down the number of iterations. Here is the code for making a bounding box tree from a list of polygons:


import math

from collections import deque


def mgcoord3d_getitem(self, i):

if i == 0:

return self.x

elif i == 1:

return self.y

elif i == 2:

return self.z

else:

return None


mgcoord3d.__getitem__ = mgcoord3d_getitem


def make_bounding_box(poly_list):

box = None

for poly in poly_list:

vertex = mgGetChild(poly)

while vertex:

ok, x, y, z = mgGetVtxCoord(vertex)

pt = mgMakeCoord3d(x, y, z)

if box == None:

box = mgMakeBox(pt, pt)

else:

mgBoxExpandCoord3d(box, pt)

vertex = mgGetNext(vertex)

return box


def compute_sum(pt_list):

partial_sum = None

if pt_list:

partial_sum = mgCoord3dZero()

for pt in pt_list:

partial_sum = mgCoord3dAdd(partial_sum, pt)

return partial_sum

def compute_mean(pt_list):

mean = None

sum = compute_sum(pt_list)

if sum:

mean = mgCoord3dDivide(sum, len(pt_list))

return mean

def compute_absolute_deviation(pt_list, m):

abs_dev = None

if pt_list:

abs_dev = mgCoord3dZero()

for pt in pt_list:

abs_dev.x += math.fabs(pt.x - m.x)

abs_dev.y += math.fabs(pt.y - m.y)

abs_dev.z += math.fabs(pt.z - m.z)

return abs_dev

def compute_centroid(poly):

pt_list = []

vertex = mgGetChild(poly)

while vertex:

ok, x, y, z = mgGetVtxCoord(vertex)

pt_list.append(mgMakeCoord3d(x, y, z))

vertex = mgGetNext(vertex)

return compute_mean(pt_list)


def get_max_component_index(pt):

if pt.x >= pt.y and pt.x >= pt.z:

return 0

if pt.y >= pt.x and pt.y >= pt.z:

return 1

if pt.z >= pt.x and pt.z >= pt.x:

return 2


class BBT:

root = None

adj_list = {}

leaf_to_poly = {}

ignore_func = None

ignore_func_arg = None

iterator_depth = -1

def split_bounding_box(self, bounding_box):

poly_list = self.adj_list[bounding_box][0]

poly_list_length = len(poly_list)

if (poly_list_length == 1):

self.leaf_to_poly[bounding_box] = poly_list[0]

self.adj_list[bounding_box] = []

elif (poly_list_length > 1):

centroid_list = self.adj_list[bounding_box][1]

centroid_mean = compute_mean(centroid_list)

centroid_abs_dev = compute_absolute_deviation(centroid_list, centroid_mean)

max_abs_dev_index = get_max_component_index(centroid_abs_dev)

poly_sublist_high = []

centroid_sublist_high = []

poly_sublist_low = []

centroid_sublist_low = []

for i in range(poly_list_length):

poly = poly_list

centroid = centroid_list

centroid_comp = centroid[max_abs_dev_index]

centroid_mean_comp = centroid_mean[max_abs_dev_index]

if math.fabs(centroid_comp - centroid_mean_comp) <= 0.0001:

if len(poly_sublist_high) < len(poly_sublist_low):

poly_sublist_high.append(poly)

centroid_sublist_high.append(centroid)

else:

poly_sublist_low.append(poly)

centroid_sublist_low.append(centroid)

elif centroid_comp > centroid_mean_comp:

poly_sublist_high.append(poly)

centroid_sublist_high.append(centroid)

elif centroid_comp < centroid_mean_comp:

poly_sublist_low.append(poly)

centroid_sublist_low.append(centroid)

bounding_box_high = make_bounding_box(poly_sublist_high)

bounding_box_low = make_bounding_box(poly_sublist_low)

self.adj_list[bounding_box_high] = [poly_sublist_high, centroid_sublist_high]

self.adj_list[bounding_box_low] = [poly_sublist_low, centroid_sublist_low]

self.adj_list[bounding_box] = [bounding_box_high, bounding_box_low]

self.split_bounding_box(bounding_box_high)

self.split_bounding_box(bounding_box_low)

__split_bounding_box = split_bounding_box

def __init__(self, poly_list = None):

if poly_list:

centroid_list = []

for poly in poly_list:

centroid = compute_centroid(poly)

if centroid:

centroid_list.append(centroid)

else:

print "error: failed to compute the centroid of " + mgGetName(poly)

return

bounding_box = make_bounding_box(poly_list)

self.root = bounding_box

self.adj_list = {bounding_box : [poly_list, centroid_list]}

self.__split_bounding_box(bounding_box)

def __iter__(self):

if self.root:

if not (self.ignore_func and self.ignore_func(self.root, self.ignore_func_arg)):

# make a queue consisted of ordered pairs of the form (box, depth) to

# enable breadth first search

self.queue = deque([(self.root, 0)])

return self

def next(self):

if self.queue:

head = self.queue[0]

child_list = self.adj_list[head[0]]

self.queue.popleft()

for box in child_list:

if not (self.ignore_func and self.ignore_func(box, self.ignore_func_arg)):

self.queue.append((box, head[1] + 1))

self.iterator_depth = head[1]

return head[0]

else:

self.iterator_depth = -1

raise StopIteration

def get_poly(self, box):

return self.leaf_to_poly.get(box)


Suppose your collection of points is a Python list of mgcoord2d objects, named pt_list, and your collection of polygons is the currently selected polygons in the current database. Then the following code snippet prints out the answers you are looking for:


def compute_height_above_ground(box_tree, pt):

def ignore_func(box, pt):

return not (box.min.x <= pt.x and box.max.x >= pt.x and box.min.y <= pt.y and box.max.y >= pt.y)


height_list = []

box_tree.ignore_func = ignore_func

box_tree.ignore_func_arg = pt

for box in box_tree:

poly = box_tree.get_poly(box)

if poly:

height = compute_height_above_ground_in_poly(poly, pt)

if height:

height_list.append(height)

return height_list


poly_list = []

db = mgGetCurrentDb()

select_list = mgGetSelectList(db)

if mgGetRecListLevel(select_list) == fltPolygon:

for i in range(mgGetRecListCount(select_list)):

poly, mat = mgGetNextRecInList(select_list)

poly_list.append(poly)

box_tree = BBT(poly_list)

if box_tree:

for pt in pt_list:

height_list = compute_height_above_ground(box_tree, pt)

print "Height above ground at (", pt.x, ",", pt.y, "): ", height_list

Login to post a comment