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_hainesdef 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 insidedef 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 mathfrom collections import dequedef mgcoord3d_getitem(self, i):    if i == 0:        return self.x    elif i == 1:        return self.y    elif i == 2:        return self.z    else:        return Nonemgcoord3d.__getitem__ = mgcoord3d_getitemdef 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 boxdef 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 2class 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_listpoly_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 `