gable_roof_extraction.py

## gable roof extraction
#
# Script for detecting gable roof segment pairs (neighboring segments with
# similar inclination but opposite orientation).
#
# Workflow:
# 1. Reads a point cloud file and stores the intermediate results in a .odm file.
# 2. Computes normals and extracts roof segments (alpha shapes, bounding boxes).
# 3. Identifies potential gable roof pairs based on the angle between segment normals.
# 4. Computes ridge lines by intersecting the corresponding planes.
# 5. Visualizes the segment pairs and ridge lines, either iteratively
#    (argument: --mode iterative) or together in a single figure (--mode all).
# 6. Plots can be shown interactively, saved, or both (--outputMode show|save|both).
#
# For an overview of all parameters, run:
#     python gable_roof_extraction.py --help
 
 
 
import opals
from opals import Histo, Import, Normals, Segmentation, pyDM
import numpy as np
import math
import matplotlib.pyplot as plt
from opals.tools.opalsargparse import ArgumentParser, Boolean
import os
 
#----------setting arguments for parsing----------
 
 
 
def get_args():
    parser = ArgumentParser(description="The program reads a point cloud file, extracts possible gable roof segments and calculates the ridge line. The output is a 3d plot of the roof segments and their ridge lines")
 
    # input file
    parser.add_argument("-inFile",
                        type=str,
                        required=True,
                        help="Point cloud input file to be analyzed"
                        )
 
    # output
    parser.add_argument("-outFile",
                        default=None,
                        type=str,
                        help="Base name for generated output files")
 
    # preprocessing
    parser.add_argument("-searchRadius",
                        type=float,
                        default=1.0,
                        help="Search radius (in meters) for normal computation and segmentation"
                        )
    parser.add_argument("-neighbors",
                        type=int,
                        default=8,
                        help="Number of neighbors used for normal vector computation"
                        )
    parser.add_argument("-skipIfExists",
                        type=Boolean,
                        default=True,
                        help="Skip Import, Normals, and Segmentation steps if outputs already exist"
                        )
    # segmentation
    parser.add_argument("-minSegSize",
                        type=int,
                        default=50,
                        help="Minimum number of points required per segment"
                        )
 
    parser.add_argument("-maxDist",
                        type=float,
                        default=0.05,
                        help="Maximum allowed point-to-plane distance during plane extraction"
                        )
 
    parser.add_argument("-maxSigma",
                        type=float,
                        default=0.15,
                        help="Maximum sigma0 value allowed during plane adjustment"
                        )
 
    parser.add_argument("-normalSigma0",
                        type=float,
                        default=0.02,
                        help="Points with a NormalSigma0 below this threshold are used as seed candidates"
                        )
 
    parser.add_argument("-alphaRadius",
                        type=float,
                        default=2,
                        help="Radius for alpha shape construction"
                        )
 
    parser.add_argument("-radius",
                        type=float,
                        default=1,
                        help="Search radius for local homogeneity checks")
    # analysis
    parser.add_argument("-searchExpand",
                        type=float,
                        default= 2.0,
                        help="Expands the bounding box by this value (in meters) when searching for intersecting segments"
                        )
    # visualization
    parser.add_argument("-mode",
                        type=str,
                        default="all",
                        choices=["all", "iterative"],
                        help="Plot mode: 'all' for a single figure, 'iterative' for one figure per pair"
                        )
 
    parser.add_argument("-outputMode",
                        type=str,
                        default="show",
                        choices=["show", "save", "both"],
                        help="Output mode: show plots, save them as .svg, or both"
                        )
 
    args = parser.parse_args()
 
    if args.outFile is None:
        base = os.path.splitext(os.path.basename(args.inFile))[0]
        args.outFile = base
 
    return args
 
 
#------------functions------------
 
## @package python.demo.gable_roof_extraction
# Creates a search box for a spatial query and expands it to make sure to find intersecting segments
def get_boxes(geometry, expand):
    box = pyDM.Box(geometry)
    box.expand(expand)
 
    return box
 
# [Checking Roof Pairs section] start
 
##
# Checks if normals might correspond to gable roof pair
def check_gable_roof_normals(nx1, ny1, nz1, nx2, ny2, nz2, cogz1, cogz2):
    angle_thres = math.pi / 180 * 5 # tolerance of 5 degrees
    rel_distance_thres = 0.10
    height_diff_thres = 0.1
 
    if None in (nx1, ny1, nz1, nx2, ny2, nz2):
        return False
    if nz1 < 0:
        nx1, ny1, nz1 = -nx1, -ny1, -nz1
    if nz2 < 0:
        nx2, ny2, nz2 = -nx2, -ny2, -nz2
 
    alpha1 = math.atan2(ny1, nx1)
    d1 = math.hypot(nx1, ny1)
    alpha2 = math.atan2(ny2, nx2) + math.pi
    d2 = math.hypot(nx2, ny2)
    d = max(d1, d2)
    dalpha = alpha1 - alpha2
 
    while dalpha < -math.pi:
        dalpha += 2 * math.pi
    while dalpha > math.pi:
        dalpha += 2 * -math.pi
 
    angle_crit = abs(dalpha) < angle_thres
    distance_crit = (1 - d1 / d) < rel_distance_thres and (1 - d2 / d) < rel_distance_thres
    heightdiff_crit = abs(cogz1-cogz2) < height_diff_thres
 
    return angle_crit and distance_crit and heightdiff_crit
 
# [Checking Roof Pairs section] finished
 
##
# Computes the intersection line between two roof segments (returns one point and a direction vector).
# Because of large coordinate values, reduced coordinates are used 
def compute_ridge_line(c1, n1, c2, n2):
    n1 = np.asarray(n1, dtype=float)
    n2 = np.asarray(n2, dtype=float)
    c1 = np.asarray(c1, dtype=float)
    c2 = np.asarray(c2, dtype=float)
 
    # plane 1: NormalX[1]*Xred + NormalY[1]*Yred + NormalZ[1]*Zred = n1 * Xred = 0
    # plane 2: NormalX[2]*Xred + NormalY[2]*Yred + NormalZ[2]*Zred + offset = 0
    # offset = - NormalX[2]*(cogX[2]- cogX[1]) - NormalY[2]*(cogY[2]- cogY[1]) -NormalZ[2]*(cogZ[2]- cogZ[1])
    offset = -np.dot(n2, c2 - c1)
 
    # direction of itersection line
    v = np.cross(n1, n2)
    v_squared = np.dot(v, v)
 
    # base point, d1,d2 are the distance from origin
    d1 = 0
    d2 = -offset
    p0_red = (d1 * np.cross(n2, v) + d2 * np.cross(v, n1)) / v_squared
    p0_global = p0_red + c1
 
    return p0_global, v
 
def compute_ridge_line_segment(p0, v, poa1, poa2):
    lambdas = []
    v_2 = np.dot(v, v)
    for i in range(len(poa1[0])):
        p = np.array([poa1[0][i], poa1[1][i], poa1[2][i]])
        lam = np.dot((p - p0), v) / v_2
        lambdas.append(lam)
 
    for i in range(len(poa2[0])):
        p = np.array([poa2[0][i], poa2[1][i], poa2[2][i]])
        lam = np.dot((p - p0), v) / v_2
        lambdas.append(lam)
 
    p1 = p0 + min(lambdas) * v
    p2 = p0 + max(lambdas) * v
 
    return p1, p2
 
def set_axes_radius(ax, origin, radius):
    ax.set_xlim3d([origin[0] - radius, origin[0] + radius])
    ax.set_ylim3d([origin[1] - radius, origin[1] + radius])
    ax.set_zlim3d([origin[2] - radius, origin[2] + radius])
 
##
# Makes axes of a 3D plot have equal scale so that spheres appear as spheres,
# cubes as cubes, etc.  This is one possible solution to Matplotlib's
# ax.set_aspect('equal') and ax.axis('equal') not working for 3D.
# Input:    ax: a matplotlib axis, e.g., as output from plt.gca().
#     
def set_axes_equal(ax):
    limits = np.array([
        ax.get_xlim3d(),
        ax.get_ylim3d(),
        ax.get_zlim3d(),
    ])
 
    origin = np.mean(limits, axis=1)
    radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
    set_axes_radius(ax, origin, radius)
 
def main():
    args = get_args()
 
    if args.skipIfExists == False or not os.path.exists(f"{args.outFile}_segments.odm"):
        # import file
        imp = Import.Import()
        imp.inFile = args.inFile
        imp.outFile = f"{args.outFile}.odm"
        imp.run()
        imp.reset()
 
        # calculating normal vectors for segmentation
        norm = Normals.Normals()
        norm.inFile = f"{args.outFile}.odm"
        norm.neighbours = int(args.neighbors)
        norm.searchRadius = float(args.searchRadius)
        norm.run()
        norm.reset()
 
        # creates roofsegments and their alpha shapes as well as bounding boxes
        seg = Segmentation.Segmentation()
        seg.inFile = f"{args.outFile}.odm"
        seg.searchRadius = float(args.radius)
        seg.minSegSize = int(args.minSegSize)
        seg.method = opals.Types.SegmentationMethod.planeExtraction
        seg.planeExtraction.maxDist = float(args.maxDist)
        seg.planeExtraction.maxSigma = float(args.maxSigma)
        seg.alphaRadius = int(args.alphaRadius)
        seg.planeExtraction.seedCalculator = f"NormalSigma0<{args.normalSigma0}  ? NormalSigma0 : invalid"
        seg.byproduct = f"{args.outFile}_segments.odm"
        seg.run()
        seg.reset()
 
    # define layout for iterating over odm file
    f = pyDM.AddInfoLayoutFactory()
    f.addColumn(pyDM.ColumnSemantic.SegmentID)
    f.addColumn(pyDM.ColumnType.bool_, "_AlphaShape")
    f.addColumn(pyDM.ColumnType.double_, "_cogx")
    f.addColumn(pyDM.ColumnType.double_, "_cogy")
    f.addColumn(pyDM.ColumnType.double_, "_cogz")
    f.addColumn(pyDM.ColumnSemantic.NormalX)
    f.addColumn(pyDM.ColumnSemantic.NormalY)
    f.addColumn(pyDM.ColumnSemantic.NormalZ)
    f.addColumn(pyDM.ColumnType.uint32, "_PointCount")
    layout = f.getLayout()
 
    # load roofsegment.odm and Polylines
    odm = pyDM.Datamanager.load(f"{args.outFile}_segments.odm", True, False)
    pi = odm.getPolylineIndex()
 
    # [Searching section] start
 
    #searching for gable roof pairs
    pairs = {} #uses a dictionary to remove redundant pairs via next()
    for geom in odm.geometries(layout):
        box = get_boxes(geom, expand= float(args.searchExpand))
        info = geom.info()
        # makes sure to just use bounding boxes, alpha shapes do not have normal vectors or cog
        if info.get(1):
            continue
        boxID = info.get(0)
        nx1, ny1, nz1 = info.get(5), info.get(6), info.get(7)
        cogz1 = info.get(4)
        pointCountA = info.get(8)
 
        window = pyDM.Window(box)
        segs = pi.searchGeometry(window, pyDM.SpatialQueryMode.intersect) #gets all other boxes intercepting with the given box
 
        for seg in segs:
            seg.cloneAddInfoView(layout, True)
            sinfo = seg.info()
            if sinfo.get(1):
                continue
            segID = sinfo.get(0)
            nx2, ny2, nz2 = sinfo.get(5), sinfo.get(6), sinfo.get(7)
            cogz2 = info.get(4)
            pointCountB = sinfo.get(8)
 
            if check_gable_roof_normals(nx1, ny1, nz1, nx2, ny2, nz2, cogz1, cogz2):
                total_points = pointCountA + pointCountB
                pair = tuple(sorted((boxID,segID)))
                existing_pair = next((p for p in pairs.keys() if boxID in p or segID in p), None)
                if existing_pair:
                    if pairs[existing_pair] < total_points:
                        del pairs[existing_pair]
                        pairs[pair] = total_points
                else:
                    pairs[pair] = total_points
    print(f"found {len(pairs)} potential roof pairs")
 
    # [Searching section] finished
 
    # collect bounding boxes and necessary attributes for ridgeline computation
    boxes = {}
    for geom in odm.geometries(layout):
        info = geom.info()
        if info.get(1):
            continue
        segID = info.get(0)
        n = (info.get(5), info.get(6), info.get(7))
        c = (info.get(2), info.get(3), info.get(4))
        b = pyDM.Box(geom)
        boxes[segID] = dict(n=n, c=c, xmin=b.xmin, ymin=b.ymin, xmax=b.xmax, ymax=b.ymax)
 
    #collect points of alpha shapes
    points_of_alphashape = dict()
    for geom in odm.geometries(layout):
        info = geom.info()
        if not info.get(1):
            continue
        segID = info.get(0)
 
        all_x = [pt.x for pt in geom.points()]
        all_y = [pt.y for pt in geom.points()]
        all_z = [pt.z for pt in geom.points()]
        points_of_alphashape[segID] = [all_x, all_y, all_z]
 
 
    # [Computation section] start
    # computing ridge lines
    ridge_lines = []
    for segA, segB in pairs.keys():
        A, B = boxes[segA], boxes[segB]
        p0, v = compute_ridge_line(A["c"], A["n"], B["c"], B["n"])
 
        p1, p2 = compute_ridge_line_segment(p0, v, points_of_alphashape[segA], points_of_alphashape[segB])
 
        ridge_lines.append((segA, segB, tuple(p1.tolist()), tuple(p2.tolist())))
 
    print(f"{len(ridge_lines)} ridge lines found")
 
    # [Computation section] finished
 
    # function for visualising
 
    def plot_results(pairs, ridge_lines, layout, odm, title=""):
        fig = plt.figure(figsize=(10, 8))
        ax = fig.add_subplot(111, projection= "3d")
        # get segment IDs
        paired_ids = set()
        colors_of_pair = dict()
        for pair in pairs:
            color = np.random.rand(3,)
            for segID in pair:
                paired_ids.add(segID)
                colors_of_pair[segID] = color
 
        #plot gable roof surface
        for geom in odm.geometries(layout):
            info = geom.info()
            segID = info.get(0)
            isAlpha = info.get(1)
            if not isAlpha or segID not in paired_ids:
                continue
 
            #color = np.random.rand(3,)
 
            # gets points from sub-segments
            for part in geom.parts():
                xs = [pt.x for pt in part.points()]
                ys = [pt.y for pt in part.points()]
                zs = [pt.z for pt in part.points()]
                ax.plot(xs, ys, zs, color=colors_of_pair[segID], linewidth=1.0, alpha=0.7)
 
        # drawing ridge lines
        for (a, b, p1, p2) in ridge_lines:
            if (a, b) in pairs or (b, a) in pairs:
                xs = [p1[0], p2[0]]
                ys = [p1[1], p2[1]]
                zs = [p1[2], p2[2]]
                ax.plot(xs, ys, zs, color="k", linewidth=1.0)
 
        set_axes_equal(ax)
        plt.title(title)
        plt.axis("off")
        return fig, ax
 
    # visualization
    mode = args.mode
    output_mode = args.outputMode
    figures = []
 
    if mode == "all":
        # creates one plot containing all roof pairs
        print("creating one plot to show all roof pairs")
        fig, ax = plot_results(pairs.keys(), ridge_lines, layout, odm, "all roof pairs")
        filename=f"{args.outFile}_all_roof_pairs.svg"
        if output_mode in ("save", "both"):
            plt.savefig(filename)
        if output_mode == "save":
            plt.close(fig)
        if output_mode in("show", "both"):
            plt.show()
 
    elif mode == "iterative":
        # creates a .svg file for each pair
        print("creating one plot for each of the pairs")
        for i, pair in enumerate(pairs.keys()):
            fig, ax = plot_results([pair], ridge_lines, layout, odm, f"Roof pair {i + 1}: {pair}")
            filename = f"{args.outFile}_roof_pair_{i + 1}.svg"
            if output_mode in ("save", "both"):
                plt.savefig(filename)
            if output_mode == "save":
                plt.close(fig)
            else:
                figures.append(fig)
        if output_mode in ("show", "both"):
            plt.show()
 
        plt.close(fig)
 
if __name__ == "__main__":
    main()