search

mfbonfigli / gocesiumtiler

A Cesium.js point cloud 3D tiles generator from LAS files written in Golang
Mozilla Public License 2.0
197 stars 40 forks source link

readme

gocesiumtiler

Build Status

                                             _   _ _
  __ _  ___   ___ ___  ___(_)_   _ _ __ ___ | |_(_) | ___ _ __
 / _  |/ _ \ / __/ _ \/ __| | | | | '_   _ \| __| | |/ _ \ '__|
| (_| | (_) | (_|  __/\__ \ | |_| | | | | | | |_| | |  __/ |
 \__, |\___/ \___\___||___/_|\__,_|_| |_| |_|\__|_|_|\___|_|
  __| | A Cesium Point Cloud tile generator written in golang
 |___/

gocesiumtiler is a tool to convert point cloud stored as LAS files to Cesium.js 3D tiles ready to be streamed, automatically generating the appropriate level of details and including additional information for each point such as color, laser intensity and classification.

What's new: V2.0.0

gocesiumtiler V2.0.0 has been released introducing several important improvements over V1. Please refer to the changelog for a full list of changes.

This release is backward incompatible compared to V1 and also deprecates several options, most of which were not of much interest. Some of these might be added in future minor updates of V2.

Features

gocesiumtiler V2 offers the following features:

New from v2.0.0-gamma: Convertion between EGM and WGS84 elevations is now delegated to Proj. Just make sure to specify the correct input vertical datum and make sure the share folder contains the required vertical grids. These are not included but can be downloaded from the Proj CDN

The currently tool uses the version 9.5.0 of the well-known Proj.4 library to handle coordinate conversion. The input CRS specified by just providing the relative EPSG code, a Proj4 string or WKT text, or letting the library figure it out automatically from the LAS metadata.

Speed is a major concern for this tool, thus it has been chosen to store the data completely in memory. If you don't have enough memory the tool will fail, so if you have really big LAS files and not enough RAM it is advised to split the LAS in smaller chunks to be processed separately.

Information on point intensity and classification is stored in the output tileset Batch Table under the propeties named INTENSITY and CLASSIFICATION.

Demo

You can preview a couple of tilesets generated from this tool at this website

Changelog

Version 2.0.1
Version 2.0.0

Note: The final V2.0.0 version significantly differs from "pre-release" versions V2.0.0 alpha, beta and gamma. The V2.0.0 final changelog consolidates all changes from V1.

Precompiled Binaries

Along with the source code a prebuilt binary for Windows x64 and Linux x64 is provided for each release of the tool in the github page. Binaries for other systems at the moment are not provided.

Environment setup and compiling from sources

Please note that the development setup has changed considerably from V2.0.0 gamma and now leverages Docker for reproducible builds.

In general local development requires:

Please refer to the DEVELOPMENT.md file for further info on how to setup a local development environment in both Windows and Linux.

Installation instructions

  1. Download the latest version of the executable from the Releases section in github and unzip it
  2. Make sure the executable is in the same folder where the share folder is.
  3. (Optional) If you need special grids to convert your data (e.g. in case of EGM to WGS84 elevation conversion or some less common projections), please download the Proj Data grids from the Proj CDN and unpack them in the share folder.
  4. Execute the binary tool with the appropriate flags.

CLI Usage

To run just execute the binary tool with the appropriate flags.

To show the help run:

gocesiumtiler -help

Commands

There are two commands, file and folder:

Flags

Common flags

These flags are applicable to both the file and the folder commands

   --out value, -o value                  full path of the output folder where to save the resulting Cesium tilesets
   --crs value, --epsg value, -e value    String representing the input CRS. For example, and EPSG code like EPSG:4326 or EPSG:28355+5773 or a generic Proj4 or WKT string. Bare numbers will be interpreted as EPSG codes. If empty the system will attempt to autodetect the CRS from the LAS metadata. In case of multiple LAS files, the CRS must be consistent else an error will be thrown.
   --resolution value, -r value           minimum resolution of the 3d tiles, in meters. approximately represets the maximum sampling distance between any two points at the lowest level of detail (default: 20)
   --z-offset value, -z value             z offset to apply to the point, in meters. only use it if the input elevation is referred to the WGS84 ellipsoid or geoid (default: 0)
   --depth value, -d value                maximum depth of the output tree. (default: 10)
   --min-points-per-tile value, -m value  minimum number of points to enforce in each 3D tile (default: 5000)
   --8-bit                                set to interpret the input points color as part of a 8bit color space (default: false)  
   --subsample value                      Approximate percent of points to keep in the final point cloud, between 0.01 (1%) and 1 (100%) (default: 1)
   --help, -h                             show help

Folder command flags

These commands are specific to the folder command:

   --join, -j                             merge the input LAS files in the folder into a single cloud. The LAS files must have the same properties (CRS etc) (default: false)

A note on vertical coordinate conversion

Previous releases of gocesiumtiler had a dedicated flag for EGM to WGS84 ellipsoid elevation conversion. This has been deprecated and now the vertical datum conversion is fully delegated to Proj. This means that to convert the vertical coordinates in case they are not referred to the WGS84 ellipsoid the input CRS definition needs to include the definition for the vertical datum.

If the LAS file includes a 3D CRS definition or separate Horizontal and Vertical CRS definitions gocesiumtiler will attempt to use those info automatically. Else if a vertical definition is not included, or the tool fails to detect the correct CRS, the right definition can be provided manually via the --crs flag.

This can be achieved for example by providing a composite CRS EPSG code.

For example, if a LAS file contains coordinates in the EPSG:32633 CRS with elevations referred to the EGM2008 (EPSG:3855), run the command with:

--crs EPSG:32633+3855

NOTE: This might require to install extra grid files in the share folder. In case of issues please make sure please to download the Proj Data grids from the Proj CDN and save them in the share folder.

Usage examples:

Example 1

Convert all LAS files in folder C:\las, write output tilesets in folder C:\out, assume LAS input coordinates expressed in EPSG:32633, with elevations referred to the EGM2008 (EPSG 3855). Use a resolution of 10 meters, create tiles with minimum 1000 points and enforce a maximum tree depth of 12:

gocesiumtiler folder -out C:\out -crs EPSG:32633+3855 -resolution 10 -min-points-per-tile 1000 -depth 12 C:\las

or, using the shorthand notation:

gocesiumtiler folder -o C:\out -e EPSG:32633+3855 -r 10 -m 1000 -d 12 C:\las

Example 2

Like Example 1 but merge all LAS files in C:\las into a single 3D Tileset

gocesiumtiler folder -out C:\out -crs EPSG:32633+3855 -resolution 10 -min-points-per-tile 1000 -depth 12 -join C:\las

or, using the shorthand notation:

gocesiumtiler folder -o C:\out -e EPSG:32633+3855 -r 10 -m 1000 -d 12 -j C:\las

Example 3

Convert a single LAS file at C:\las\file.las, write the output tileset in the folder C:\out, use the system defaults and let gocesiumtiler extract the CRS metadata from the LAS.

gocesiumtiler file -out C:\out C:\las\file.las

or, using the shorthand notation:

gocesiumtiler file -o C:\out -e 32633 C:\las\file.las

Library Usage in other GO programs

To use the tiler in other go programs just:

go get github.com/mfbonfigli/gocesiumtiler/v2

Then instantiate a tiler object and launch it either via ProcessFiles or ProcessFolder passing in the desired processing options.

A mnimal example is:

package main

import (
    "context"
    "log"

    "github.com/mfbonfigli/gocesiumtiler/v2/tiler"
)

func main() {
    t, err := tiler.NewGoCesiumTiler()
    if err != nil {
        log.Fatal(err)
    }
    ctx := context.TODO()
    err = t.ProcessFiles([]string{"myinput.las"}, "/tmp/myoutput", "EPSG:32632", tiler.NewTilerOptions(
        tiler.WithEightBitColors(true),
        tiler.WithWorkerNumber(2),
        tiler.WithMaxDepth(5),
    ), ctx)
    if err != nil {
        log.Fatal(err)
    }
}

Note that you will require to use cgo for the compilation, for how to setup the build environment please refer to the DEVELOPMENT.md.

Mutators

gocesiumtiler from version 2.0.0 final offers the concept of mutators. Mutators are implementations of the mutator.Mutator interface and can be used to manipulate or discard input points.

The library vends a ZOffset mutator to perform vertical traslation of point clouds and a Subsampler mutator to thin down the points in the output.

Other possible uses of mutators (not yet built in into the library) could be, for example:

To use the mutators just pass them to the tiler options:

mut := []mutator.Mutator{
   mutator.NewZOffset(10),
   mutator.NewSubsampler(0.5),
}

err = t.ProcessFiles([]string{"myinput.las"}, "/tmp/myoutput", "EPSG:32632", tiler.NewTilerOptions(
   tiler.WithEightBitColors(true),
   tiler.WithWorkerNumber(2),
   tiler.WithMaxDepth(5),
   tiler.WithMutators(mut),
), ctx)

To implement a custom mutator, implement the mutator interface:

Mutate(p model.Point, t model.Transform) (model.Point, bool)

The Mutate function receives as input the original point and a Transform object that can be used to trasform from the local CRS to the global EPSG 4978 CRS and back. The output of the Mutate function is the, eventually manipulated, point and a boolean which should be true if the point should appear in the final point cloud, false otherwise.

Note: mutators must be goroutine safe.

Example of a custom mutator: coloring points by class

type ClassBasedColorMutator struct {}

func (m *ClassBasedColorMutator) Mutate(pt model.Point, localToGlobal model.Transform) (model.Point, bool) {
    switch pt.Classification {
    case 2:
        // ground > brown
        pt.R, pt.G, pt.B = 128, 0, 0
    case 3:
        // low vegetation > dark green
        pt.R, pt.G, pt.B = 0, 60, 0
    case 4:
        // medium vegetation > green
        pt.R, pt.G, pt.B = 0, 130, 0
    case 5:
        // high vegetation > light green
        pt.R, pt.G, pt.B = 30, 170, 30
    case 6:
        // building > orange
        pt.R, pt.G, pt.B = 200, 100, 0
    case 9:
        // water > blue
        pt.R, pt.G, pt.B = 5, 170, 200
    default:
        // all other classes > grey
        pt.R, pt.G, pt.B = 30, 30, 30
    }
    return pt, true
}

Algorithms

The sampling occurs using a hybrid, lazy octree data structure. The algorithm works as follows:

  1. All points are stored in a linked tree backed by an underlying array: this provides efficient list manipulations operations (splitting, adding, removing) and avoids dynamic allocations of slices. The coordinates are internally converted to EPSG 4978 and then to a local CRS with Z-axis normal to the WGS84 ellipsoid. The backing array helps with CPU cache friendliness and helps achieving measurable speed improvements of up to 20% compared to a standard linked list.
  2. An octree cell is created. Every cell stores N points (variable). A cell also has a grid spacing property. The root node has a spacing set to the provided resolution.
  3. The points are traversed. Each point will fall into one of the cells the space has been divided by the given grid spacing. If the point is the closest one the the center of the cell, it's taken as new closest for that cell, else it's discarded and parked into a list of the Node octant it belongs to.
  4. When all points are traversed the points closes to the cells the space has been partitioned in will be the points for the current tree node. All others are parked. If an octant however has a number of parked points that is less than the min-points-per-node, the parked points for that octant are rolled up to the current node. Similarly if the current node has a depth equal to the configured max depth.
  5. Whenever the children are retrieved, the previously parked points are used to create child nodes on demand using the same algorithm, lazily.
  6. The points are written in the final artifacts with coordinates relative to the local CRS, however a global transform is applied at the tileset root to convert points back to the EPSG 4978 CRS required by Cesium.

Precompiled Binaries

Along with the source code, a prebuilt binary for both Linux and Windows x64 is provided for each release of the tool in the github page.

Future work and support

Further work needs to be done, such as:

Contributors and their ideas are welcome.

If you have questions you can contact me at m.federico.bonfigli@gmail.com

Versioning

This library uses SemVer for versioning. For the versions available, see the tags on this repository.

Credits

Massimo Federico Bonfigli - Github

License

This project is licensed under the Mozilla Public License 2.0 - see the LICENSE.md file for details.

The software uses third party code and libraries. Their licenses can be found in LICENSE-3RD-PARTIES.md file.

Acknowledgments

License

Lhe library is released under the Mozilla Public License Version 2.0. See the LICENSE.md file for further info.