Terrain editing on the fly

[ShayBM_Datumate]

Hi all

We are using the latest Cesium version in our website, with a custom terrain provider that we manage. The terrain is used for display purposes, and also for very exact calculations.
We are trying to develop a feature where the user will be able to select draw a polygon on the map and the terrain in that specific area will be flattened instantly.

Does anyone know what is the best way to achieve this functionality?

[omar]

What is the source data for your terrain provider? If you’re using heightmaps you could edit that source and re-generate the terrain.

This forum thread implies using a shader to alter the Z coordinate of the tiles may work but I haven’t seen an implementation of this: How to flat the overlap area between a polygon and a 3D tiles to a specific height?.

I’m curious to hear more about your goals with this feature. Is the idea to flatten it so you can add new data there instead?

[ShayBM_Datumate]

Hi Omar, thank you for replying

Here are answers to your questions:

• We are using quantized mesh, not heightmaps, since we noticed they look better.
• The forum thread you are referring to is indeed what I try to achieve, but not on 3D tiles but on the terrain, because we are not displaying 3D tiles, and also we are using the terrain for high resolution calculations (for example: calculating and displaying cross sections of up to 2 cm accuracy).
• The goal in this feature is to allow the user to fix issues or inaccuracies that he finds in the terrain.

I ended up solving this by extending CesiumTerrainProvider and overriding the requestTileGeometry method. I went over all of the mesh vertexes and if they were inside the polygon the user selected, I interpolated their height using the height of the polygon’s perimeter vertexes.

The results are looking good (notice the curb area between the parking spots):

Before:

After:

However, I have a few things I need to do to improve this solution, and I could use your help:

1. Sometimes the results do not look as good as the images above, probably because I’m using the same mesh with different heights instead of creating a new mesh. I need to find a performance efficient way of recreating it.
2. I want to know what tiles are intersecting the user’s polygon so I won’t have to do this for irrelevant tiles.
3. I want to reload only relevant terrain tiles instead of removing and re-adding the entire terrain for it to update.

Thank you for all your help!

[omar]

That looks pretty good! If you can post what edits you needed to make to the source code (either here or in a new topic in the General category) that could be very helpful for others in the community, and perhaps help push this as a feature in the library if there is interest.

@ShayBM_Datumate

I want to know what tiles are intersecting the user’s polygon so I won’t have to do this for irrelevant tiles.

This isn’t exposed in the public API but this method describes how to get all terrain tiles in view: openlayers 2 - Cesium 3D - Determining the Map Scale of the viewed Globe - Geographic Information Systems Stack Exchange. You can get the bounds of each and do your intersection that way.

@ShayBM_Datumate

I want to reload only relevant terrain tiles instead of removing and re-adding the entire terrain for it to update.

I’m not sure there’s an easy way to do this without making changes to the source.

[ShayBM_Datumate]

@omar

You can get the bounds of each and do your intersection that way.

I found a better way to achieve this. Since I have the level, x and y of the tile in the requestTileGeometry method, and also have access to the terrain provider tiling scheme, I used this._tilingScheme.tileXYToRectangle(x, y, level) to get the tile’s rectangle and checked if it overlaps the user’s polygon.

@omar

I’m not sure there’s an easy way to do this without making changes to the source.

Bummer… but I guess it won’t be so bad to have the terrain reloaded. We’ll see when I’m done

@ShayBM_Datumate

Sometimes the results do not look as good as the images above, probably because I’m using the same mesh with different heights instead of creating a new mesh. I need to find a performance efficient way of recreating it.

Right now I’m focused on this issue. I’m trying to break the mesh triangles that cross the polygon’s perimeter into smaller triangles so that the mesh will look good.

Will update soon with the results and post some of the code here

[ShayBM_Datumate]

Hi,
Sorry for the delay, I didn’t have the time to write this until now.

Here’s the gist of my solution:

import { Request, TerrainData, Rectangle } from 'cesium';
import * as turf from '@turf/turf';
import { Feature, Polygon } from '@turf/turf';

const MAX_SHORT = 32767;

interface ModelEdit {
  polygon: Feature<Polygon>;
  polygonTriangles: Feature<Polygon>[]; // TIN algorithm result on polygon
}

export class TerrainProviderEdit extends Cesium.CesiumTerrainProvider {
  private modelEdits: ModelEdit[] = [];

  constructor({ url, modelEdits }: { url: string; modelEdits: ModelEdit[] }) {
    super({ url });

    this.modelEdits = modelEdits;
  }

  requestTileGeometry(x: number, y: number, level: number, request: Request): Promise<TerrainData> | undefined {
    const promise = super.requestTileGeometry(x, y, level, request);
    if (!promise || this.modelEdits.length === 0) {
      return promise;
    }

    const tileRectangle: Rectangle = this.tilingScheme.tileXYToRectangle(x, y, level);
    const tilePolygon = GeoUtils.rectangleToPolygon(tileRectangle);  // Create turf polygon from tile rectangle
    const relevantEdits = this.modelEdits.filter(
      edit => turf.booleanOverlap(edit.polygon, tilePolygon) || turf.booleanContains(edit.polygon, tilePolygon)
    );
    if (relevantEdits.length === 0) {
      return promise;
    }

    return promise.then((data: TerrainData) => this.modifyTerrainTile(data, tileRectangle, relevantEdits));
  }

  private modifyTerrainTile(terrainData: TerrainData, tileRectangle: Rectangle, modelEdits: ModelEdit[]) {
    const data = terrainData as any;
    const minimumHeight = data._minimumHeight;
    const maximumHeight = data._maximumHeight;

    const quantizedVertices: Uint16Array = data._quantizedVertices;
    const vertexCount = quantizedVertices.length / 3;

    const positions: number[][] = [];
    for (let i = 0; i < vertexCount; i++) {
      const rawU = quantizedVertices[i];
      const rawV = quantizedVertices[i + vertexCount];
      const rawH = quantizedVertices[i + vertexCount * 2];

      const u = rawU / MAX_SHORT;
      const v = rawV / MAX_SHORT;
      const longitude = Cesium.Math.toDegrees(Cesium.Math.lerp(tileRectangle.west, tileRectangle.east, u));
      const latitude = Cesium.Math.toDegrees(Cesium.Math.lerp(tileRectangle.south, tileRectangle.north, v));

      let height = Cesium.Math.lerp(minimumHeight, maximumHeight, rawH / MAX_SHORT);
      const currentPoint = turf.point([longitude, latitude]);
      const relevantEdit = modelEdits.find(edit => turf.booleanPointInPolygon(currentPoint, edit.polygon));
      if (relevantEdit) {
        const relevantTriangle = relevantEdit.polygonTriangles.find(triangle =>
          turf.booleanPointInPolygon(currentPoint, triangle)
        );
        if (relevantTriangle) {
          height = turf.planepoint(currentPoint, relevantTriangle);
        }
      }
      positions.push([longitude, latitude, height]);
    }

    // TODO: split mesh triangles that are crossing the user polygon's perimiter and recreate mesh

    const heights = positions.map(p => p[2]);
    const newMinHeight = Math.min(...heights);
    const newMaxHeight = Math.max(...heights);

    const newQuantizedVertices = new Uint16Array(positions.length * 3);
    positions.forEach((p, i) => {
      const lonRad = Cesium.Math.toRadians(p[0]);
      newQuantizedVertices[i] = Math.round(
        Cesium.Math.lerp(MAX_SHORT, 0, (lonRad - tileRectangle.east) / (tileRectangle.west - tileRectangle.east))
      );

      const latRad = Cesium.Math.toRadians(p[1]);
      newQuantizedVertices[i + positions.length] = Math.round(
        Cesium.Math.lerp(MAX_SHORT, 0, (latRad - tileRectangle.north) / (tileRectangle.south - tileRectangle.north))
      );

      const relativeHeight = Math.round(
        Cesium.Math.lerp(0, MAX_SHORT, (p[2] - newMinHeight) / (newMaxHeight - newMinHeight))
      );
      newQuantizedVertices[i + positions.length * 2] = relativeHeight;
    });

    data._minimumHeight = newMinHeight;
    data._maximumHeight = newMaxHeight;
    data._quantizedVertices = newQuantizedVertices;

    return data as TerrainData;
  }
}

A few notes:

• I removed some of the code because it was specific to my app, not the general solution to the problem I described in this post (So copy-pasting it probably won’t work ).
• An additional argument is passed to the constructor - a list of ModelEdits which are polygons the user created and wants to be edited. In addition, every polygon is passed with a list of triangles that is the result of the TIN algorithm on that polygon.
• To make this work, every time the user creates a new model edit you should remove and destroy the previous terrain and add a new one, with the new model edit (in addition to any previous ones).
• I added a TODO in the middle of the code to note that this solution is not ideal. It doesn’t break the mesh’s triangles into smaller ones, just changes the elevation of the existing triangles’ points. This causes issues if the user’s polygon intersects any mesh triangles. This is mostly noticeable when trying to edit elements such as bridges, but for our use cases it was fine since the user mostly wants to get rid of holes or small mounds.

If anyone has any questions feel free to ask
Shay

[cesiumTs]

Hello, I see that the flattening effect of your terrain is very good. I am developing this function, but the code you provided cannot be used normally at present. Could you please provide a more detailed demo?Thank you very much.