Skip to content

Conversation

@ansonl
Copy link
Contributor

@ansonl ansonl commented Aug 19, 2025

I have refined Chris' difference mesh generation code, removed some redundant steps, and added some new steps needed to get more seamless meshes. I have also refactored the code and class organization a bit to play better with code editor autocomplete such as VSCode intellisense.

This is a draft of my changes to Touch Terrain. I'll edit this PR and add my changes for a little while until it is ready to merge.

Difference Mesh

Tweaked Chris' code to get the "difference mesh" to generate at the correct height with the right walls. The difference mesh is the mesh made from the subtraction of a "top" and "bottom" DEM. The walls are created at the right locations by dilating twice.

Two modes are referred to in the code:

Normal Mode

Generates a mesh for the passed in raster with bottom surface flat along the base.

top_elevation_hint should be passed when generating a raster that will be the "bottom" raster of a related Difference Mesh.

Difference Mesh Mode

Generates a mesh that is sandwiched between two rasters.

Config values related to Difference Mesh mode are:

importedDEM_interp = None
"Optional raster file for interpolating at edges"
top_elevation_hint = None
"elevation raster for the future top of the model that would be used for a future difference mesh. Used for Normal mode where Difference Mesh will be created in the future with the same top raster."
bottom_elevation = None
"elevation raster for the bottom of the model."
bottom_floor_elev: None|float = None
"Set bottom raster to an elevation in locations where bottom is NaN but top raster is not NaN. Defaults to min_elev-1. If set to less than min_elev, difference mesh at that point will go thru base."
bottom_thru_base
"if mesh should drop thru to base (as if the bottom was 0)"
dirty_triangles
"allow degenerate triangles for difference mesh"

Mesh Limitations

TouchTerrain's current mesh generation is limited to creating connected volumes with quads as the "top" and "bottom" surfaces of the volume. Real vertical quads connecting a top and bottom quad as "walls" are only allowed at the true X/Y borders of a connected area. This means that a vertical quad (wall) is not generated to connect vertices of 2 quads that are both on the "top" or "bottom" surface.

   (✓ yes)       (✗ no)           (✗ no)
 ̅ ̅ ̅\                ̅ ̅ ̅\                  ̅ ̅ ̅\
    \___              |______              |
                                            \___

↓ This limitation leads to mesh generations where 2 volumes share edges. Because an edge should only have 2 faces, this is technically non-manifold.
tt-non-manifold

image

The 2 volumes are normally closed on their own and the only connection between them are shared edges with no shared volume going through them so this should not cause issues for most usages. It's easy to tell where one volume begins and ends by eye. Hopefully our mesh processing software is smart enough.

The mesh size reduction method in #99 still works despite these technically non-manifold edges because Blender correctly recognizes these edges as part of the "outside edge border".

In the future, we should look into adding a way for top/bottom quad of a cell to connect to the quad along the same surface of another cell by a vertical quad.

  • A mostly well defined case would be when cell's bottom quad center value is at the base and it could connect to bordering bottom quads' edges by a vertical quad to avoid the non-manifold edges seen in the animation above where the 2 width gap straddling the highlighted edges that should be along the base is forced to have a "transition" distance of at least 1 cell to get cells that are completely along the base (or actually gapped like in the picture).
    image
image
  • The code doesn't check connectivity between faces right now since the points between faces are assumed connected. If we have a gap for a vertical wall due to forcing the bottom quad to the base, when generating vertices for a quad in the final step, we can check if quads next to each other have touching corner vertices. If the shared vertices are not connected, generate a quad wall to bridge that gap?

Make the difference mesh manifold at edges

Remove zero height volumes that occur near the edges of the difference mesh (where the top and bottom mesh have the same Z coordinates) so that the mesh is manifold. This cleanup is in cell.remove_zero_height_volumes()

The removal in NW and SE directions works best when splitting edge rotation is set to favor less steep split edges. Thus remove_zero_height_volumes() is only active and needed when split_rotation = 1

↓ Before with zero volume shapes
before_empty_volumes2
↓ After removing zero volume shapes
after_no_empty_volumes2

Quad Triangulation Splitting Edge Rotation

Add feature to rotate the splitting edge when triangulating quads so that the edge orientation can favor less/more steep split edges + flat borders. This occurs in quad.get_triangles()

The original behavior of constant edge orientation from NW-SE is still the default behavior.

↓ Before
single_triangle_direction
↓ After
two_triangle_direction

split_rotation: None | int = None
"""Should quad triangulation rotate the splitting edge based on the slope of the created edge?
None -> NW>SW edges
1 -> Rotate for less steep along split edges > Steeper faces along the split.
2 -> Rotate for more steep along split edges > Less steep faces along the split.
"""

Fixed W/E X coordinate flipped bug

The W/E locations for X coordinate in createCells() was flipped and noticed in a previous code comment. This made the assignment of the cell's quad's corner vertices confusing because we were flipping the W/E back again during the vertex creation.

9297859

This is fixed in the linked commit of this PR so that the X coordinates assigned for W/E now reflect the right W or E directions.

Code refactor

Partially updated the code to PEP standard. Such as UpperCaseClassName for class names and lower_case_underscore for variable names. It makes Touch Terrain easier to read and debug in code editors with Python linting such as VSCode.

Python type hints (Requires Python 3.10+)

Added Python type hints to make code readability and debugging easier. Types have been assigned to some of the variables but more type hints need to be added in the future to make the code "pass" compile time type checking.

image

Hovering over type hinted variables now tells you what it represents and what it is used for.

Python docstrings

Added new and converted existing docstrings to be formatted correctly.

image

Numpy

Some files imported numpy as np and some had just numpy so I changed the references in TouchTerrainEarthEngine.py to leave numpy as numpy to be consistent.

Touch Terrain Config

Configuration is defined as a class TouchTerrainConfig instead of a dict[str, Any]. Config values centrally managed in a single location in user_config.py as class attributes allow type hinting and hover to view documentation on that value.

The class attribute also makes it less likely to make a typo with IDE autocomplete versus manually typing a text key for dictionary.

Touch Terrain Tile Info

tile_info is is defined as a class TouchTerrainTileInfo instead of a dict[str, Any]. Similar benefits as the transition to TouchTerrainConfig. TouchTerrainTileInfo is defined in tile_info.py.

TouchTerrainConfig is now stored in tile infos under TouchTerrainTileInfo.config so that the config values can be accessed in a single point of truth instead of copying each config value into a dictionary.

Multiple Raster Version Management with RasterVariants

Lots of raster operations done in Touch Terrain affect some or all of the rasters for a single DEM. Keeping track of multiple variants of a DEM with variables like top, top_nan, top_dilated is easy to forget making a change to a single variant.

Raster variants for the same DEM in various stages of processing are kept in RasterVariants which each variant stored as an attirbute such as original, nan_close, and dilated. There is also a new edge_interpolation which uses the DEM from an optional importedDEM_interp config option to interpolate the top vertex values for the bottom_thru_base case.

class RasterVariants:
    """Holds a raster with processed copies of it
    """
    
    original: Union[None, np.ndarray] # Original full raster
    nan_close: Union[None, np.ndarray] # Raster after NaN close values to bottom and before dilation
    dilated: Union[None, np.ndarray] # Raster after dilation
    
    edge_interpolation: Union[None, np.ndarray] # Original full raster with values past edges for interpolation

RasterVariants supports +, -, * operators so all raster variants stored as numpy.ndarray can be modified at once as if RasterVariants was a single numpy.ndarray. In the example below, all bottom DEM raster variants' arrays are increased by a constant value.

# Before without RasterVariant
bottom += 100
bottom_nan += 100
bottom_dilated += 100
bottom_edge_interpolation += 100

# After using RasterVariant
bottom_raster_variants += 100

ProcessingTile

ProcessingTile contains attributes about the processing tile that were previously passed passed as individual parameters. Now it is easier to pass data along for tile logic by adding it to this new class instead of creating new parameters.

class ProcessingTile:
    tile_info: TouchTerrainTileInfo
    top_raster_variants: RasterVariants
    bottom_raster_variants: Union[None, RasterVariants]

Other Changes/Notes

DEM_name config for locally imported DEM and config_path (automatically set in standalone mode)

zip_file_name or DEM_name or config_path (checked in that order) is now used as the exported ZIP filename.

DEM_name or config_path (checked in that order) is now used as the mesh filename.
If only one tile is generated, the exported filename does not include the _tile_1_1 at the end of the filename.

DEM_name: None | str = 'USGS/3DEP/10m'
"name of DEM source used in Google Earth Engine. for all valid sources, see DEM_sources in TouchTerrainEarthEngine.py. Also used if specifying a custom mesh and zip and extracted folder name."
config_path: None | str = None
"The path of the Touch Terrain config file. Set this during runtime. If DEM_name is None or default value, use config filename for default zip and mesh filenames and unzipped folder name."

tileScale config option

The map scale can now be directly specified with tileScale config option. Specified tile scale takes precedence over tilewidth.

Other PRs

Incorporated #109, #106

Test environment

I used conda and virtual environment for testing and running TouchTerrain. For anyone else working with a local directory install of TouchTerrain in a similar environment, I recommend cloning the repo into its own folder.

git clone XXX
cd ./TOUCHTERRAIN_FOR_CAGEO
# With a new virtual environment called touchterrain-dev
conda activate touchterrain-dev
# Update conda env with touch terrain requirements from environment.yml
conda env update --name touchterrain-dev --file environment.yml --prune
# Install touchterrain as a module in "editable" state so it links back to the local development code folder
pip install -e .

I keep the data files in a separate directory at the top level named like touchterrain-dev so the the repo code folder and the data folder are in the same top level directory.

To run TouchTerrain with from TouchTerrain_standalone.py in the code folder but use the data folder as the working directory, you can reference the python file in the code folder while in data folder. Like python ../TOUCHTERRAIN_FOR_CAGEO/TouchTerrain_standalone.py

Or open VSCode in the code folder (workspace folder) and use a VSCode launch configuration setup to debug in the data folder (cwd):

{
      "name": "Python Debugger: TT Standalone config.json",
      "type": "debugpy",
      "request": "launch",
      "program": "${workspaceFolder}/TouchTerrain_standalone.py",
      "console": "integratedTerminal",
      "cwd": "${workspaceFolder}/../touchterrain-dev/",
      "args": "config.json"
}
image

I have added a VSCode launch.json to the repo that I use. It has some left over test cases right now and I will clean up the file and add some test configs in the future. My launch.json is setup to run with the JSON config files and DEMs in the ../touchterrain-dev (data) folder relative to the repo (code) folder as described above.

ansonl added 3 commits August 19, 2025 23:22
…d raster add/minus/scaling in grid.__init__(). Removed tile_info.user_offset because it was redundant to minus the minimum height of the raster then add user_offset (which is the minimum height of raster - user defined min_elev).
…zip_file_name takes precedence for ZIP and folder name
@ansonl ansonl marked this pull request as ready for review August 22, 2025 05:59
nirabo added a commit to nirabo/TouchTerrain_for_CAGEO that referenced this pull request Oct 14, 2025
Testing and validation complete:
- All tests passing (7 passed, 20 skipped)
- Dependencies verified
- Pre-commit formatting applied
- Integration validated

Deferred 138 linting errors in PR ChHarding#111 core files to post-integration cleanup.

Ready for pull request creation.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>
@nirabo nirabo mentioned this pull request Oct 14, 2025
93 tasks
ansonl added 28 commits October 17, 2025 22:56
Find the intersection polygon between each raster cell and the clipping polygon. Sort all individual edges of intersection polygons into buckets stored in RasterVariants based on if the edge lies on a cardinal direction edge of the cell quad. Marks all interior edges as needing walls created.

Add shapely_utils with functions for flattening shapely geometries
…ke_wall"

Created python test case for the inner most wall determination function
…k at horizontal and vertical neighbor cell and were looking at diagonal neighbor by applying direction YX at same time.

Worked on test case for intersection function
…de if direction is out of range which happens when the cell is on the very edge of the raster.
…ll creation.

Moved quad and vertex into their own files for import because my new functions in shapely_utils.py need to import quad but they can't import existing quad from grid_tesselate due to a circular reference since the new functions are called from grid_tesselate. Simplified the code for marking walls and handle case where cell1 or cell2 are disjoint or contained properly in boundary.
…. Changed wall marking to mark walls on the Polygon side of a match instead of cell 1 so that we know which vertex order to create the eventual wall edge for outward normal direction.
…erties for surfacePolygons explicitly as None so that the property exists to check for existence
… the cell is outside the clipping boundary.

The top_hint is used in raster_preparation and affects the top.original raster so it's cells that are out of bounds need to match the top.original cell boundary when they are out of bounds.
…ll_for_walls to handle case where cell1 has no buckets when cell1 is outside the boundary (or otherwise NaN from source) or cell1 is contained properly
…tion borders of the vertex that is deleted because if an entire tri is removed due to having the same Z, there is no longer those 2 borders touching it and there is no new diagonal "wall/border" because the the diagonal edge of the top and bottom have the same Z
…ce mode bottom raster interpolation values match the corresponding "top" raster of the normal mode
…irements.txt. Remove requirements.txt from active use by renaming it. Update pyproject.toml metadata.
…on. RasterVariants polygon_intersection_geometry now only represents disjoint and intersection (partial+complete) cases. Wall marking handles partial/complete intersection case by checking dilated elevation raster to see if we are at the edge of the generated mesh to see if we need a wall.
…raster_preparation relies on a clipped raster and wall marking relies on dilated raster
This occured because we stopped differentiating between partial intersection and full enclosure (contains_properly) for cells with the existence of polygon_intersection_geometry. We can't differentiate between partial and full enclosure that way anymore so we create a new 2d array under RasterVariants to store if a cell is fully contained by the clipping polygon.
This way create_cell() can use that new array to check if a cell is fully contained and use the "quad" for the top/bot of the generated mesh instead of the surface polygon which is now present for all cells that intersect the clipping polygon in some way.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant