jeremybyu / polylidar Goto Github PK

import polylidar
import math
import numpy as np
from polylidar import MatrixDouble, Polylidar3D
import laspy
import time
import matplotlib.pyplot as plt
from polylidar.polylidarutil import (generate_test_points, plot_triangles, get_estimated_lmax,
                                     plot_triangle_meshes, get_triangles_from_list, get_colored_planar_segments, plot_polygons)

from polylidar.polylidarutil import (plot_polygons_3d, generate_3d_plane, set_axes_equal, plot_planes_3d,
                                     scale_points, rotation_matrix, apply_rotation)
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Read LAS file
# las_file_path = 'house1.las'
las_file_path = '/home/tarun/Downloads/house1.las'
las_file = laspy.read(las_file_path)

# Extract x, y, and z coordinates from the LAS points
points = np.vstack((las_file.x, las_file.y, las_file.z)).transpose()

polylidar_kwargs = dict(alpha=0.0, lmax=10, min_triangles=20, z_thresh=0.0, norm_thresh_min=0.94)
polylidar = Polylidar3D(**polylidar_kwargs)

points_mat = MatrixDouble(points, copy=False)
t1 = time.time()
mesh, planes, polygons = polylidar.extract_planes_and_polygons(points_mat)
t2 = time.time()
print("Took {:.2f} milliseconds".format((t2 - t1) * 1000))

triangles = np.asarray(mesh.triangles)
fig, ax = plt.subplots(figsize=(10, 10), nrows=1, ncols=1,
                       subplot_kw=dict(projection='3d'))
# plot all triangles
plot_planes_3d(points, triangles, planes, ax)
plot_polygons_3d(points, polygons, ax)
# plot points
ax.scatter(*scale_points(points), c='k', s=0.1)
set_axes_equal(ax)
ax.view_init(elev=15., azim=-35)
plt.show()
print("")

  File "examples/python/realsense_mesh.py", line 196, in <module>
    main()
  File "examples/python/realsense_mesh.py", line 188, in main
    pcd = pcd.rotate(R_Standard_d400[:3, :3], center=False)
TypeError: rotate(): incompatible function arguments. The following argument types are supported:
    1. (self: open3d.open3d_pybind.geometry.Geometry3D, R: numpy.ndarray[float64[3, 3]], center: numpy.ndarray[float64[3, 1]]) -> open3d.open3d_pybind.geometry.Geometry3D

Invoked with: geometry::PointCloud with 24259 points., array([[ 1.,  0.,  0.],
       [ 0.,  0.,  1.],
       [ 0., -1.,  0.]]); kwargs: center=FalseTraceback (most recent call last):
  File "examples/python/realsense_mesh.py", line 196, in <module>
    main()
  File "examples/python/realsense_mesh.py", line 188, in main
    pcd = pcd.rotate(R_Standard_d400[:3, :3], center=False)
TypeError: rotate(): incompatible function arguments. The following argument types are supported:
    1. (self: open3d.open3d_pybind.geometry.Geometry3D, R: numpy.ndarray[float64[3, 3]], center: numpy.ndarray[float64[3, 1]]) -> open3d.open3d_pybind.geometry.Geometry3D

Invoked with: geometry::PointCloud with 24259 points., array([[ 1.,  0.,  0.],
       [ 0.,  0.,  1.],
       [ 0., -1.,  0.]]); kwargs: center=False

"""Paths and patches"""

from matplotlib.patches import PathPatch
from matplotlib.path import Path
from numpy import asarray, concatenate, ones


class Polygon(object):
    # Adapt Shapely or GeoJSON/geo_interface polygons to a common interface
    def __init__(self, context):
        if isinstance(context, dict):
            self.context = context['coordinates']
        else:
            self.context = context

    @property
    def exterior(self):
        return (getattr(self.context, 'exterior', None)
                or self.context[0])

    @property
    def interiors(self):
        value = getattr(self.context, 'interiors', None)
        if value is None:
            value = self.context[1:]
        return value


def PolygonPath(polygon):
    """Constructs a compound matplotlib path from a Shapely or GeoJSON-like
    geometric object"""

    def coding(ob):
        # The codes will be all "LINETO" commands, except for "MOVETO"s at the
        # beginning of each subpath
        n = len(getattr(ob, 'coords', None) or ob)
        vals = ones(n, dtype=Path.code_type) * Path.LINETO
        vals[0] = Path.MOVETO
        return vals

    if hasattr(polygon, 'geom_type'):  # Shapely
        ptype = polygon.geom_type
        if ptype == 'Polygon':
            polygon = [Polygon(polygon)]
        elif ptype == 'MultiPolygon':
            polygon = [Polygon(p) for p in polygon]
        else:
            raise ValueError(
                "A polygon or multi-polygon representation is required")

    else:  # GeoJSON
        polygon = getattr(polygon, '__geo_interface__', polygon)
        ptype = polygon["type"]
        if ptype == 'Polygon':
            polygon = [Polygon(polygon)]
        elif ptype == 'MultiPolygon':
            polygon = [Polygon(p) for p in polygon['coordinates']]
        else:
            raise ValueError(
                "A polygon or multi-polygon representation is required")

    vertices = concatenate([
        concatenate([asarray(t.exterior.coords)[:, :2]] +
                    [asarray(r)[:, :2] for r in t.interiors])
        for t in polygon])
    codes = concatenate([
        concatenate([coding(t.exterior.coords)] +
                    [coding(r) for r in t.interiors]) for t in polygon])

    return Path(vertices, codes)


def PolygonPatch(polygon, **kwargs):
    """Constructs a matplotlib patch from a geometric object

    The `polygon` may be a Shapely or GeoJSON-like object with or without holes.
    The `kwargs` are those supported by the matplotlib.patches.Polygon class
    constructor. Returns an instance of matplotlib.patches.PathPatch.

    Example (using Shapely Point and a matplotlib axes):

      >>> b = Point(0, 0).buffer(1.0)
      >>> patch = PolygonPatch(b, fc='blue', ec='blue', alpha=0.5)
      >>> axis.add_patch(patch)

    """
    return PathPatch(PolygonPath(polygon), **kwargs)