feat: add Galaxy-ready CSP inventory CLI

.dockerignore

@@ -1,4 +1,5 @@
 .git
+.github
 .Rproj.user
 .Rhistory
 .RData
@@ -14,4 +15,6 @@
 .vscode
 __pycache__
 data/
-tests/testthat/_snaps/
+tests
+tests/testthat/_snaps/
+inst/extdata

.github/workflows/container.yml

@@ -0,0 +1,37 @@
+name: Container
+
+on:
+  push:
+    branches: [main]
+    tags: ["v*"]
+  workflow_dispatch:
+
+permissions:
+  contents: read
+  packages: write
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: docker/setup-buildx-action@v3
+      - uses: docker/login-action@v3
+        with:
+          registry: ghcr.io
+          username: ${{ github.actor }}
+          password: ${{ secrets.GITHUB_TOKEN }}
+      - uses: docker/metadata-action@v5
+        id: metadata
+        with:
+          images: ghcr.io/3dtrees-earth/3dtrees_csp_standsegmentation
+          tags: |
+            type=raw,value=0.2.0,enable=${{ github.ref == 'refs/heads/main' }}
+            type=ref,event=tag
+            type=sha
+      - uses: docker/build-push-action@v6
+        with:
+          context: .
+          push: true
+          tags: ${{ steps.metadata.outputs.tags }}
+          labels: ${{ steps.metadata.outputs.labels }}

Dockerfile

@@ -0,0 +1,56 @@
+FROM rocker/geospatial:4.4.3
+
+ENV DEBIAN_FRONTEND=noninteractive \
+    RGL_USE_NULL=TRUE
+
+ARG RLAS_COMMIT=82cbba42f158d1dfc91efda3207923260a052564
+
+RUN apt-get update \
+    && apt-get install --yes --no-install-recommends \
+        libgl1-mesa-dev \
+        libglu1-mesa-dev \
+        libx11-dev \
+        libxt-dev \
+        curl \
+        patch \
+        python3-pip \
+    && rm -rf /var/lib/apt/lists/*
+
+RUN python3 -m pip install --break-system-packages --no-cache-dir --no-deps \
+        laspy==2.6.1 \
+        lazrs==0.8.1
+
+RUN install2.r --error --skipinstalled --ncpus -1 \
+        BH \
+        RANN \
+        RCSF \
+        Rcpp \
+        RcppArmadillo \
+        colorspace \
+        conicfit \
+        data.table \
+        dbscan \
+        doParallel \
+        foreach \
+        geometry \
+        igraph \
+        jsonlite \
+        lidR \
+        magrittr \
+        rgl \
+        testthat
+
+COPY patches/rlas-read-all-extrabytes.patch /tmp/rlas-read-all-extrabytes.patch
+RUN curl -L -sS "https://github.com/r-lidar/rlas/archive/${RLAS_COMMIT}.tar.gz" -o /tmp/rlas.tar.gz \
+    && mkdir -p /tmp/rlas-src \
+    && tar -xzf /tmp/rlas.tar.gz -C /tmp/rlas-src --strip-components=1 \
+    && patch -d /tmp/rlas-src -p1 < /tmp/rlas-read-all-extrabytes.patch \
+    && R CMD INSTALL /tmp/rlas-src \
+    && rm -rf /tmp/rlas-src /tmp/rlas.tar.gz /tmp/rlas-read-all-extrabytes.patch
+
+WORKDIR /opt/CspStandSegmentation
+COPY . /opt/CspStandSegmentation
+RUN chmod -R a+rX /opt/CspStandSegmentation \
+    && R CMD INSTALL --no-multiarch --with-keep.source /opt/CspStandSegmentation
+
+CMD ["Rscript", "/opt/CspStandSegmentation/exec/run.R", "--help"]

R/forest_inventory.R

@@ -37,6 +37,62 @@ suppress_cat <- function(f, ...) {
   f(...)                        # Call the function and capture its return value
 }
 
+# Quiet, vectorized equivalent of conicfit::CircleFitByPratt(). The upstream
+# implementation emits diagnostics with cat(), which previously required
+# opening and sinking /dev/null for every RANSAC iteration. Forest inventory
+# calls this thousands of times, making output suppression more expensive than
+# the fit itself.
+.circle_fit_by_pratt_quiet <- function(XY) {
+  XY <- as.matrix(XY)
+  centroid <- colMeans(XY)
+  centered <- sweep(XY, 2L, centroid, "-")
+  Xi <- centered[, 1L]
+  Yi <- centered[, 2L]
+  Zi <- Xi * Xi + Yi * Yi
+
+  Mxx <- mean(Xi * Xi)
+  Myy <- mean(Yi * Yi)
+  Mxy <- mean(Xi * Yi)
+  Mxz <- mean(Xi * Zi)
+  Myz <- mean(Yi * Zi)
+  Mzz <- mean(Zi * Zi)
+  Mz <- Mxx + Myy
+  Cov_xy <- Mxx * Myy - Mxy * Mxy
+  Mxz2 <- Mxz * Mxz
+  Myz2 <- Myz * Myz
+  A2 <- 4 * Cov_xy - 3 * Mz * Mz - Mzz
+  A1 <- Mzz * Mz + 4 * Cov_xy * Mz - Mxz2 - Myz2 - Mz * Mz * Mz
+  A0 <- Mxz2 * Myy + Myz2 * Mxx - Mzz * Cov_xy - 2 * Mxz * Myz * Mxy + Mz * Mz * Cov_xy
+  A22 <- A2 + A2
+  epsilon <- 1e-12
+  ynew <- 1e20
+  xnew <- 0
+
+  for (iter in seq_len(20L)) {
+    yold <- ynew
+    ynew <- A0 + xnew * (A1 + xnew * (A2 + xnew * xnew * 4))
+    if (abs(ynew) > abs(yold)) {
+      xnew <- 0
+      break
+    }
+    Dy <- A1 + xnew * (A22 + 16 * xnew * xnew)
+    xold <- xnew
+    xnew <- xold - ynew / Dy
+    if (abs((xnew - xold) / xnew) < epsilon) break
+    if (iter >= 20L || xnew < 0) {
+      xnew <- 0
+      break
+    }
+  }
+
+  DET <- xnew * xnew - xnew * Mz + Cov_xy
+  center <- c(
+    Mxz * (Myy - xnew) - Myz * Mxy,
+    Myz * (Mxx - xnew) - Mxz * Mxy
+  ) / DET / 2
+  matrix(c(center + centroid, sqrt(sum(center * center) + Mz + 2 * xnew)), nrow = 1L)
+}
+
 #' RANSAC circle fitting algorithm specially adapted for tree DBH estimation
 #'
 #' This function fits a circle to a set of points using the RANSAC algorithm it maximizes the points that are in the circle and the number of filled 36 degree angle segments
@@ -87,9 +143,7 @@ ransac_circle_fit <- function(data,n_iterations = 1000L,distance_threshold = 0.0
 
     # Fit circle; keep tryCatch very tight and avoid pipe
     circle <- tryCatch(
-      {
-        CspStandSegmentation::suppress_cat(conicfit::CircleFitByPratt, sample_points)
-      },
+      .circle_fit_by_pratt_quiet(sample_points),
       warning = function(w) NULL,
       error   = function(e) NULL
     )
@@ -536,4 +590,3 @@ plot_inventory <- function(plot, inventory, col = NA, cex = 1.5, label_col = "wh
     rgl::lines3d(c(inventory$X[i] - plot[1], inventory$X[i] - plot[1]), c(inventory$Y[i] - plot[2], inventory$Y[i] - plot[2]), c(inventory$Z[i], inventory$Height[i]), col = ifelse(length(col) >= i, col[i], col), lwd = 2)
   }
 }
-

README.md

@@ -3,7 +3,114 @@
 Authors: Julian Frey and Zoe Schindler, University of Freiburg, Chair of Forest Growth and Dendroecology
 
 
-[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.17294732.svg)](https://doi.org/10.5281/zenodo.17294732)  [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![R-CMD-check](https://github.com/JulFrey/CspStandSegmentation/actions/workflows/r.yml/badge.svg)](https://github.com/JulFrey/CspStandSegmentation/actions/workflows/r.yml)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.17294732.svg)](https://doi.org/10.5281/zenodo.17294732)  [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
+
+## 3Dtrees command-line workflow
+
+This fork retains the upstream segmentation implementation and adds a headless
+CLI and container for the 3Dtrees Galaxy tool. It always creates a DTM and can
+inventory any number of existing instance dimensions. CSP segmentation is
+optional; when enabled, the original points and dimensions are preserved and
+the output gains exactly one `PredInstance_CSP` dimension.
+
+Inventory-only runs avoid loading dimensions that cannot affect the DTM or the
+requested inventories. They retain XYZ, Classification, requested
+instance/species fields, and optional ForestMamba score/semantic fields. CSP
+runs continue to load every dimension because the emitted point cloud must
+preserve them. LASlib can select only the first nine extra-byte records by
+position; when a requested field occurs later, the reader safely falls back to
+all extra bytes, then immediately projects the in-memory cloud back to the
+required fields. The container uses the same pinned `rlas` patch as 3Dtrees
+standardization, so that fallback correctly loads every declared extra byte
+rather than stopping after nine.
+
+```bash
+Rscript exec/run.R \
+  --input input.laz \
+  --output-dir results \
+  --segmentation-spec PredInstance_SAT,species_id_SAT,species_prob_SAT \
+  --segmentation-spec PredInstance_FM,species_id_FM,species_prob_FM \
+  --enable-csp false \
+  --dtm-resolution 0.2 \
+  --random-seed 42
+```
+
+`--segmentation-spec INSTANCE[,SPECIES,SPECIES_PROB]` is repeatable. Species
+dimensions are optional per segmentation but must be supplied as a pair. When
+none are supplied, the combined inventory omits species columns and no species
+composition file is created. Run `Rscript exec/run.R --help` for all controls.
+
+Common controls are the input, repeatable segmentation specs, non-tree IDs,
+optional CSP and seed source, optional native-CRS AOI GeoJSON, DTM resolution
+(default 0.2 m), and random seed. Fine-tuning controls retain upstream defaults,
+including a 0.3 m CSP voxel and one routing worker. CSP geometry features are
+computed only when a non-zero geometry weight requires them.
+
+Outputs include:
+
+- `dtm_full.tif` and optional `dtm_aoi.tif`;
+- one TSV per instance dimension plus `inventory_combined.tsv`;
+- `stand_summary.tsv` and optional `species_composition.tsv`;
+- `effective_seeds.tsv` and `segmented_csp.laz` only when CSP is enabled;
+- `run_metadata.json` and `resource_summary.json`.
+
+Inventory rows include point count, position, height, DBH, crown convex-hull
+area, and an explicit measurement-quality field. ForestMamba inventories also
+include median `PredScore_FM`, a mixed-score flag, and wood/leaf point counts and
+shares when those source dimensions exist. Processing is fail-atomic: the final
+output directory is published only after all requested products succeed.
+
+Existing-instance inventory keeps upstream's 500 RANSAC iterations but uses a
+quiet, vectorized implementation of the same Pratt circle equations. It also
+projects only the point attributes required by each requested segmentation and
+skips the full preservation copy when CSP output is disabled. On the local
+5.9-million-point GFZ benchmark these changes reduced tool time from 59.4 to
+41.4 seconds and process peak RSS from 4.85 to 3.19 GiB.
+
+With selective reading and the standardization `rlas` patch, the same GFZ file
+successfully loaded 14 extra-byte attributes and inventoried SAT and FM together
+in 53.7 seconds at 2.14 GiB process peak RSS. Sampled CPU averaged 100.1%,
+confirming that the one-thread default consumes approximately one core. The run
+produced 64 tree rows plus DTM, stand, and species products without creating a
+point-cloud output.
+
+Inventory-only execution uses two bounded passes. DTM generation defaults to
+300 m tiles with a 5 m buffer and up to 10 workers. Below 50 million points the
+tiles read the source directly. Larger inputs are scanned once in parallel by
+point range, retaining the minimum Z at each deterministic 0.1 m cell centre;
+CSF and TIN rasterization then run on that reduced surface in parallel tiles.
+Each spatial worker uses one lidR thread, avoiding nested oversubscription. The
+result raster retains the aligned input extent, with unsupported edge areas as
+NoData.
+
+The second pass reads only selected inventory fields into disk-backed hash
+partitions by instance ID. Each tree remains complete even when its points span
+spatial chunks. CSP continues to use the full-cloud path to preserve upstream
+global voxel routing and the optional point-cloud output.
+
+On dataset 2056 (1,141,911,324 points, 3.7 GB compact LAZ), the memory-focused
+streaming candidate stage completed in 56.8 seconds with 10 workers, averaged
+9.51 CPU cores, and retained 1,267,661 candidates. Parent peak RSS was 0.78 GiB
+and the conservative sum of all worker peaks was 6.72 GiB under a 50 GB Docker
+limit. The 300 m + 5 m CSF/TIN stage completed in 76.4 seconds; its conservative
+aggregate worker peak was 2.98 GiB. An earlier exact-coordinate prototype
+scaled from 543.4 seconds with one worker to 306.6 seconds with two, 206.6
+seconds with four, and 127.9 seconds with ten. Replacing per-worker XY grids
+with deterministic cell centres produced the final 56.8-second result and cut
+the conservative 10-worker peak sum from 16.6 GiB to 6.72 GiB.
+
+On the 5.9-million-point GFZ reference, reducing to deterministic 0.1 m cell
+centres before CSF/TIN changed the DTM relative to direct full-cloud CSF/TIN by
+about 9.1 cm RMSE (4.0 cm median absolute difference). Automatic mode therefore
+keeps the direct spatial method for clouds below the 50-million-point threshold.
+
+Build and run the pinned container with:
+
+```bash
+docker build -t 3dtrees-csp .
+docker run --rm -v "$PWD:/work" -w /work 3dtrees-csp \
+  Rscript /opt/CspStandSegmentation/exec/run.R --help
+```
 
 
 
@@ -161,4 +268,3 @@ BibTex:
 	file = {Full Text PDF:O\:\\Research\\Projects\\Confobi_IWW\\Literatur\\lit_database\\storage\\R7Q8BFU5\\Larysch et al. - 2025 - Quantifying and mapping the ready-to-use veneer volume of European beech trees based on terrestrial.pdf:application/pdf},
 }
 ```
-