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/
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
Craig
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.