sylva.coreg

sylva.coreg

Marker-free coregistration of scan positions in forests.

Registers scans with no targets and no initial alignment, from the trees themselves.

1. Per scan (prepare_scan): a raster terrain model gives height above ground; horizontal slices through the stem band are fitted with circles and linked into a stem map; locally planar points (stems, ground, logs) are kept for ICP and foliage is dropped. Tilted scans are levelled first, with the scanner's own attitude.
2. Per pair (register_pair): reflective targets where both scans saw three or more, otherwise a global stem-map match (stem distances do not depend on the unknown transform, so a sorted pair table makes the search near-exhaustive) whose height comes from the two terrain models, refined by a robust point-to-plane ICP. A pair is accepted only if it fits over all points and above the ground, so the ground plane alone cannot confirm a wrong pair.
3. Whole survey (coregister): accepted pairs are edges of a pose graph, each weighted by the directions its surfaces constrain, solved with a Huber kernel and outlier rejection; scans left over are retried against the combined registered survey; optionally all poses are refined jointly (refine_joint). The report lists, per pair, how far apart the same tree lands from the two scans, a quality measure that needs no ground truth.

Scans can be held fixed at trusted poses, so new or badly registered positions join an existing project, and approximate poses (a RiSCAN SOP, GNSS) can serve as priors.

Examples:

>>> from sylva import coreg, read_riscan_project, riscan
>>> project = read_riscan_project("survey.PROJ")
>>> scans = project.with_scans(require_sop=False)
>>> result = coreg.coregister([p.rxp for p in scans], names=[p.name for p in scans],
...                           levelling=[p.levelling for p in scans],
...                           reflectors=[p.reflectors() for p in scans],
...                           approximate_positions=riscan.gnss_to_local([p.gnss for p in scans]))
>>> print(result.report())
>>> result.save("transforms.json")

sylva.coreg.pipeline

End-to-end coregistration of a survey.

The pipeline is staged, and each stage can run on its own:

1. Per scan (prepare_scan): a terrain model, the stem map, and a planarity-filtered subsample for ICP. The expensive part, done once per scan however many pairs are tried.
2. Per pair (register_pair): reflective targets where both scans saw them, otherwise a global stem-map match with its height taken from the two terrain models, gives a coarse transform that ICP refines; the pair is accepted only if it fits over all points and above the ground, and its terrain agrees.
3. Whole survey (coregister_prepared): accepted pairs become the edges of a pose graph, each weighted by the directions its surfaces constrain, solved with outlier rejection; scans left over are retried against the combined registered survey, and optionally every pose is refined jointly.

Every stage records its own quality, because the useful question is not "did it run" but "which scans can I trust". Stem matches take their height from the shared ground rather than from stem bases (vertical error dominates stem-based registration, as Tremblay & Béland 2018 and GlobalMatch, Wang et al. 2023, both found); edges carry the anisotropic information of their point-to-plane correspondences instead of an isotropic weight; outlier edges are rejected even while some scans are unregistered; and recovered scans are tied in by pairwise measurements. Scans can be held fixed at trusted poses, so that new scans are registered into an existing project; and approximate poses (a RiSCAN SOP, GNSS and compass) can serve as priors, refusing results that move a scanner implausibly far and placing scans that see too few stems.

CoregConfig dataclass

Settings of the whole pipeline.

ground_min_coverage class-attribute instance-attribute

ground_min_coverage: float | None = 0.8

Least fraction of the elevations from -30 to +5 degrees an azimuth of the scan must sample for the terrain in that direction to be fitted from it; see _refit_visible_ground. Tilted scans sample almost none in the azimuths along their tilt axis, upright scans all of them. Needs the scanner's position, so it applies to scans read from file or given origin. None fits the terrain from every return.

icp_min_planarity class-attribute instance-attribute

icp_min_planarity: float = 0.35

Planarity gate on the ICP points of every scan; 0 keeps foliage, which costs accuracy and can trap ICP in a local minimum.

icp_max_height class-attribute instance-attribute

icp_max_height: float = 12.0

Returns higher above ground (m) are left out of ICP: the upper canopy is the least repeatable part of the scene and the most moved by wind.

use_reflectors class-attribute instance-attribute

use_reflectors: bool = True

Try reflective targets before stems where both scans saw at least min_reflector_matches. Far more precise, but they usually cover only part of a survey.

trusted_reflector_matches class-attribute instance-attribute

trusted_reflector_matches: int = 5

A reflector match with at least this many targets, agreeing within trusted_reflector_rmse, is accepted even when ICP then fails its tests. Scans tens of metres apart share targets long before they share enough surface for ICP's fitness test, and so many targets cannot form a congruent pattern by coincidence. The ICP pose is kept if it stayed within reflector_tolerance of the targets', otherwise the targets' own. 0 always requires ICP to agree.

min_match_inliers class-attribute instance-attribute

min_match_inliers: int = 5

Stem correspondences a coarse match needs to reach ICP. Screening only saves ICP time; ICP fitness decides, so this is permissive.

max_match_ambiguity class-attribute instance-attribute

max_match_ambiguity: float = 0.8

A coarse match whose best distinctly different rival has more than this share of its inliers is ambiguous: both are refined by ICP and the pair is accepted only if one fits the above-ground points clearly better.

max_match_rmse class-attribute instance-attribute

max_match_rmse: float = 0.3

Largest scatter (m) of a coarse match's own correspondences.

max_coarse_stem_rmse class-attribute instance-attribute

max_coarse_stem_rmse: float = float('inf')

Optional cap (m) on the median horizontal disagreement of matched stems; off by default.

height_from_ground class-attribute instance-attribute

height_from_ground: bool = True

Take a stem match's vertical offset from the two terrain models where both saw ground, not from the matched stems. A stem's height is its own scan's terrain height under it plus 1.3 m, and under understory or on a slope that terrain is least certain exactly at the stems, while the shared ground between two scanners is wide and well sampled.

ground_radius class-attribute instance-attribute

ground_radius: float = 30.0

Terrain cells further than this (m) from a scanner are left out of height comparisons: far from the scanner, ground returns are sparse and grazing, and the terrain model there is mostly the understory.

min_ground_cells class-attribute instance-attribute

min_ground_cells: int = 50

Shared observed cells a height comparison needs.

max_ground_disagreement class-attribute instance-attribute

max_ground_disagreement: float = 0.25

A pair whose terrain models still differ by more than this (m, median over the shared ground) after ICP is refused: ICP slid vertically, or the match is wrong. None disables the test.

max_pair_distance class-attribute instance-attribute

max_pair_distance: float = 40.0

Pairs whose approximate positions are further apart (m) are not tried. Scale it with the range: with a 30 m gate no pair more than 35 m apart ever registered (0 of 1687).

screen_pairs class-attribute instance-attribute

screen_pairs: bool = True

Match every pair first and run ICP only on the survivors.

max_pairs_per_scan class-attribute instance-attribute

max_pairs_per_scan: int | None = None

After screening keep only each scan's best N pairs; None keeps all.

min_icp_fitness class-attribute instance-attribute

min_icp_fitness: float = 0.04

Fitness decreases with the distance between scans, since less of each scan is shared. On a VZ-400 survey with positions 30-35 m apart, scored against its reflector-based registration, every pair ICP placed with a fitness of 0.045-0.10 was within 10 cm of the reference, and a gate of 0.10 left 10 of its 14 scans unregistered.

min_icp_fitness_above_ground class-attribute instance-attribute

min_icp_fitness_above_ground: float = 0.03

Fitness over source points above fitness_min_height. Ground is planar and fits under any horizontal shift, so on open sites it can carry a wrong pair past min_icp_fitness alone; stems, branches and logs are what fix the horizontal position. On the same survey, correct pairs scored at least 0.043 and wrong ones at most 0.012.

stem_agreement_tolerance class-attribute instance-attribute

stem_agreement_tolerance: float = 0.25

If ICP leaves the matched stems this much (m) further apart than the coarse match did, the coarse transform is kept: ICP slid somewhere unrelated. Loose on purpose; tightening it rejected good refinements.

max_coarse_to_fine_shift class-attribute instance-attribute

max_coarse_to_fine_shift: float = 2.0

ICP moving further than this (m) from the coarse solution diverged.

recover_unregistered class-attribute instance-attribute

recover_unregistered: bool = True

Retry each unregistered scan against the combined registered survey: a scan can overlap several registered scans partly without overlapping any one of them enough.

refine_multiview class-attribute instance-attribute

refine_multiview: bool = False

Refine every registered pose in one joint Gauss-Newton over the point-to-plane correspondences of every accepted pair plus the matched stems (sylva.coreg.refine_joint).

refinement_max_shift class-attribute instance-attribute

refinement_max_shift: float = 0.3

A refinement moving any scan further (m) is discarded.

information_patch_points class-attribute instance-attribute

information_patch_points: float = 100.0

ICP correspondences counted as one independent observation when weighting pose-graph edges (~sylva.coreg.plane_edge_information): neighbouring residuals share the error of the surface they sample.

information_min_sigma class-attribute instance-attribute

information_min_sigma: float = 0.005

Floor (m) on the point-to-plane residual that weights an edge.

max_prior_shift class-attribute instance-attribute

max_prior_shift: float = 5.0

With priors: results putting a scanner further (m) from its prior position are refused.

max_prior_rotation class-attribute instance-attribute

max_prior_rotation: float | None = None

With priors: if set, also refuse results rotated more (degrees) from the prior. Off, since recorded headings can be wrong by any amount.

workers class-attribute instance-attribute

workers: int = 0

Scans prepared and pairs registered at once; 0 picks from the cores and free memory (memory_per_worker_gb per preparing scan).

riscan_filter class-attribute instance-attribute

riscan_filter: str = 'none'

Apply RiSCAN PRO's RXP import filter when reading (sylva.riscan.riscan_like_mask): "none", "current" (drops echoes within 0.5 m of the scanner, as a current RiSCAN import) or "legacy" (also the weak, isolated echoes the older conversion discarded, a fifth to a third of a scan).

riegl_options class-attribute instance-attribute

riegl_options: dict = field(default_factory=dict)

RXP reading (build it from a RiSCAN export settings file with reading_options): min_range, max_range (m, from the scanner), min_deviation/max_deviation, min_reflectance/max_reflectance, min_amplitude/max_amplitude, echoes, library. None bound by default.

min_points_per_scan class-attribute instance-attribute

min_points_per_scan: int = 1000

Scans with fewer points are set aside rather than processed.

to_dict

to_dict() -> dict

The whole configuration as plain data.

PairResult dataclass

Registration of one pair; transform maps scan i into scan j.

Attributes:

Name	Type	Description
matched_source, matched_target		Positions of the matched stems in each scan's own frame.
coarse_stem_rmse, fine_stem_rmse		Median matched-stem distance under the coarse and the ICP transform.
fitness_above	float	ICP fitness over source points above ground; NaN if unknown.
rival	MatchResult | None	A distinctly different coarse hypothesis nearly as well supported, resolved by ICP.
used_icp	bool	False when the coarse transform was kept over ICP's; icp then scores the coarse transform, not the refinement that was discarded.
ground_offset	float	Median height (m) of the target's terrain over the source's where both saw ground, under transform; NaN if they share too little.

trusted class-attribute instance-attribute

trusted: bool = False

Accepted on the strength of its reflector match although ICP failed (CoregConfig.trusted_reflector_matches).

stem_rmse property

stem_rmse: float

Stem agreement of the transform this pair reports.

summary

summary() -> str

One line for the report.

ScanFeatures dataclass

Everything later stages need from one scan; build with prepare_scan.

Attributes:

Name	Type	Description
icp_points	ndarray	Planar subsample for ICP (float32: every scan's stays resident).
icp_heights	ndarray	Height above ground of each ICP point.
levelling	ndarray	level_from_scan rotation applied to the raw points first; every pose is in the levelled frame, and SurveyResult.transform_for composes it back out.
origin	ndarray	Scanner position in the levelled frame (the origin for a raw scan).
error	str	Why the scan is unusable; empty if it is fine.

usable property

usable: bool

Can the scan take part in registration at all?

location

location(pose: ndarray) -> np.ndarray

Where the scanner stood in the world, given this scan's world_from_scan.

SurveyResult dataclass

Outcome of coregistering a survey.

Attributes:

Name	Type	Description
poses	list[ndarray]	world_from_levelled_scan per scan; see transform_for for the transform of the raw scan.
registered	list[bool]	Per scan, connected to the reference by accepted pairs.
edge_to_pair	list[int]	Index into pairs of each pose-graph edge.

transform_for

transform_for(index: int) -> np.ndarray

world_from_scan of one scan, applying directly to its raw points.

pair_for_edge

pair_for_edge(edge_index: int) -> PairResult | None

The pair behind a pose-graph edge.

rejected_pairs

rejected_pairs() -> list[PairResult]

Accepted pairs the global solve found inconsistent with the rest.

consistency

consistency(robust: bool = True) -> dict[tuple[int, int], float]

Horizontal stem disagreement per accepted pair under the final poses (m).

The field-usable quality check: no ground truth needed, it measures whether the same tree lands in the same place from two scans. Only the horizontal distance counts: a stem's height comes from its own scan's terrain model, which the pair's ground check covers.

Parameters:

Name	Type	Description	Default
robust	bool	Median rather than RMSE: in a dense stand a few coarse correspondences pair neighbouring trees, and would dominate an RMSE.	True

Returns:

Type	Description
dict	{(i, j): distance}.

report

report() -> str

Plain-text summary for a log or a QC record.

save

save(path: str | Path) -> Path

Write the transforms and quality as JSON (transforms.json).

Returns:

Type	Description
Path

register_pair

register_pair(source: ScanFeatures, target: ScanFeatures, config: CoregConfig | None = None, *, initial: ndarray | None = None, match: MatchResult | None = None, i: int = 0, j: int = 1, target_icp: ICPTarget | None = None) -> PairResult

Coarse-match then ICP-refine one pair of prepared scans.

Parameters:

Name	Type	Description	Default
source	ScanFeatures	Prepared scans; the result maps source into target.	required
target	ScanFeatures	Prepared scans; the result maps source into target.	required
config	CoregConfig | None	Pipeline settings.	None
initial	ndarray | None	Skip the coarse matching and start ICP here (GNSS, a target-based prealignment, a previous run).	None
match	MatchResult | None	A stem match already computed (by screening).	None
i	int	Indices recorded on the result.	0
j	int	Indices recorded on the result.	0
target_icp	ICPTarget | None	target's ICP pyramid, prepared once for all its pairs (ICPTarget, built with config.icp); built here if None.	None

Returns:

Type	Description
PairResult

load_transforms

load_transforms(path: str | Path) -> dict[str, np.ndarray]

Read the world_from_scan of every registered scan from SurveyResult.save.

Parameters:

Name	Type	Description	Default
path	str | Path	A transforms.json.	required

Returns:

Type	Description
dict	{name: (4, 4) matrix}.

prepare_scan

prepare_scan(cloud, config: CoregConfig | None = None, name: str = '', reflectors: list[Reflector] | None = None, levelling: ndarray | None = None, origin=None) -> ScanFeatures

Fit the ground, detect stems and build the ICP subsample of one scan.

Parameters:

Name	Type	Description	Default
cloud		A sylva.PointCloud, an (n, 3) array, or a path (.rxp through RiVLib, anything else through sylva.read).	required
config	CoregConfig | None	Pipeline settings.	None
name	str	Scan name, for reports.	''
reflectors	list[Reflector] | None	Targets the scan saw, in its own frame.	None
levelling	ndarray | None	(4, 4) rotation applied to the points (and targets) first, for a scanner that was not upright; the scanner's own attitude is the usual source, or the rotation of a RiSCAN SOP. Recorded, so poses can be expressed back in the scanner's frame.	None
origin		Scanner position in the (levelled) cloud; the origin for a raw scan.	None

Returns:

Type	Description
ScanFeatures

reading_options

reading_options(settings: str | Path | None = None, **bounds) -> dict

RXP reading options from a RiSCAN PRO export settings file and explicit bounds.

Parameters:

Name	Type	Description	Default
settings	str | Path | None	A RiSCAN export filter settings file (sylva.riscan.read_export_settings).	None
**bounds		min_range, max_range, min_deviation, max_deviation, min_reflectance, max_reflectance, min_amplitude, max_amplitude; None leaves a bound to the file (or open). Explicit bounds override the file.	{}

Returns:

Type	Description
dict	For CoregConfig.riegl_options.

coregister

coregister(clouds: Sequence, config: CoregConfig | None = None, *, names: Sequence[str] | None = None, pairs: Iterable[tuple[int, int]] | None = None, reflectors: Sequence[list] | None = None, approximate_positions: ndarray | None = None, levelling: Sequence[ndarray | None] | None = None, fixed: dict[int, ndarray] | None = None, priors: Sequence[ndarray | None] | None = None, progress: Callable[[str], None] | None = None) -> SurveyResult

Coregister a whole survey.

Parameters:

Name	Type	Description	Default
clouds	Sequence	Scans: paths, sylva.PointCloud or (n, 3) arrays.	required
config	CoregConfig | None	Pipeline settings.	None
names	Sequence[str] | None	Scan names.	None
pairs	Iterable[tuple[int, int]] | None	Pairs to try; all by default, which suits plot-scale surveys since the extra edges make the global solve robust.	None
reflectors	Sequence[list] | None	Per scan, the targets it saw (in its own frame).	None
approximate_positions	ndarray | None	(n, 3) rough scanner positions (GNSS); pairs further apart than max_pair_distance are skipped. NaN rows are unknown.	None
levelling	Sequence[ndarray | None] | None	Per scan, the (4, 4) rotation that levels it.	None
fixed	dict[int, ndarray] | None	{index: world_from_levelled_scan} of scans whose poses are trusted and held; the others are registered into their frame.	None
priors	Sequence[ndarray | None] | None	Per scan, an approximate world_from_levelled_scan (a RiSCAN SOP); see coregister_prepared.	None
progress	Callable[[str], None] | None	Called with progress messages; printed if None and verbose.	None

Returns:

Type	Description
SurveyResult	Poses are world_from_scan with the world fixed to the reference scan (or to the fixed scans).

coregister_prepared

coregister_prepared(scans: Sequence[ScanFeatures], config: CoregConfig | None = None, *, pairs: Iterable[tuple[int, int]] | None = None, approximate_positions: ndarray | None = None, fixed: dict[int, ndarray] | None = None, priors: Sequence[ndarray | None] | None = None, progress: Callable[[str], None] | None = None, started: float | None = None) -> SurveyResult

Pairwise registration and the global solve on prepared scans.

Parameters:

Name	Type	Description	Default
scans	Sequence[ScanFeatures]	From prepare_scan; the expensive part, which depends on no other scan, so it can be cached or computed elsewhere.	required
config	CoregConfig | None	Pipeline settings.	None
pairs	Iterable[tuple[int, int]] | None	Pairs to try; all usable pairs within max_pair_distance if None.	None
approximate_positions	ndarray | None	(n, 3) rough scanner positions; taken from priors if None.	None
fixed	dict[int, ndarray] | None	{index: world_from_levelled_scan} held fixed. Pairs between two fixed scans are not tried.	None
priors	Sequence[ndarray | None] | None	Per scan, an approximate world_from_levelled_scan or None. Stem matching is global and can, rarely, confirm a wrong alignment in a repetitive stand: any pair or placement that puts a scanner more than max_prior_shift from its prior position is refused. Scans whose stems do not match are placed from their prior (place_from_prior). Only the position is checked by default: recorded headings can be wrong by any amount.	None
progress	Callable[[str], None] | None	Called with progress messages.	None
started	float | None	time.perf_counter() at the start, for the reported duration.	None

Returns:

Type	Description
SurveyResult

merge_clouds

merge_clouds(clouds: Sequence, result: SurveyResult, *, voxel: float | None = 0.02, only_registered: bool = True, riegl_options: dict | None = None, riscan_filter: str = 'none')

Every scan moved into the reference frame and concatenated.

Parameters:

Name	Type	Description	Default
clouds	Sequence	The scans as passed to coregister (paths, clouds or arrays).	required
result	SurveyResult	The registration.	required
voxel	float | None	Thin each scan, and the merged cloud, to one point per voxel (m); the overlap otherwise carries several copies of every surface.	0.02
only_registered	bool	Leave out unregistered scans.	True
riegl_options	dict | None	RXP reading options (CoregConfig.riegl_options).	None
riscan_filter	str	RiSCAN import filter applied when reading (CoregConfig.riscan_filter).	'none'

Returns:

Type	Description
PointCloud	With a scan_id attribute.

place_from_prior

place_from_prior(scan: ScanFeatures, survey: Sequence[ScanFeatures], poses: Sequence[ndarray], prior: ndarray, config: CoregConfig | None = None, neighbours: int | None = None) -> tuple[PairResult, list[int]]

Place one scan from an approximate pose instead of from stems.

For scans that see too few stems to match. The prior (a RiSCAN SOP, GNSS and compass) is first corrected in height by the median offset between the scan's ground and its registered neighbours' ground, since a vertical error of metres is common and beyond ICP's reach; ICP then starts with a 1.5 m correspondence distance. Accepted on ICP fitness, above-ground fitness and RMSE, and only within max_prior_shift of the prior.

Parameters:

Name	Type	Description	Default
scan	ScanFeatures	The scan to place.	required
survey	Sequence[ScanFeatures]	Registered scans.	required
poses	Sequence[ndarray]	world_from_scan of each of survey.	required
prior	ndarray	Approximate world_from_scan of scan.	required
config	CoregConfig | None	Pipeline settings.	None
neighbours	int | None	Registered scans nearest the prior that ICP runs against (recovery_neighbours if None).	None

Returns:

Name	Type	Description
result	PairResult	transform is world_from_scan; check success.
used	list of int	Indices into survey of the scans ICP ran against.

sylva.coreg.stems

Stem maps: the features a coarse alignment is found from.

Forests defeat the usual registration cues: there are almost no large planes, and keypoint descriptors are unstable because foliage changes between viewpoints and with wind. Stems are stable, sparse, roughly vertical and described by three numbers (x, y, diameter), so reducing each scan to a stem map turns registration into a small 2-D pattern-matching problem that can be solved with no initial guess.

The detector is the Rust one behind sylva.trees.detect_stems, run in its coregistration mode (0.4 m slice spacing, no prefilter, over-wide clusters skipped, RANSAC (Fischler & Bolles 1981) in blocks of 32), so the stem maps are deterministic.

Stem dataclass

A detected stem.

x, y, z are the axis at reference_height above the ground, in the scan's own frame (z is an elevation in that frame, not a height). axis points up.

quality property

quality: float

0-1 confidence from the fit residual, the arc seen and the slices linked.

StemMap dataclass

The stems of one scan.

positions property

positions: ndarray

(n, 3) stem positions at the reference height.

sorted_by_quality

sorted_by_quality() -> StemMap

Best first (numpy's default sort, so ties always fall the same way).

top

top(n: int) -> StemMap

The n best stems; matching cost is quadratic in the count.

transformed

transformed(T: ndarray) -> StemMap

Every stem position and axis moved by a rigid transform.

save

save(path: str | Path) -> Path

Write the stems as JSON.

load classmethod

load(path: str | Path) -> StemMap

Read stems written by save.

StemDetectionConfig dataclass

Settings of detect_stems.

They suit plot-scale TLS (10-30 m range, stems 5-120 cm). In dense understorey raise min_slices and min_coverage; for buttressed stems raise slice_min_height above the buttresses.

max_cluster_extent class-attribute instance-attribute

max_cluster_extent: float = 2.0

Clusters wider than this (m) are too big to be one stem and are skipped.

detect_stems

detect_stems(points: ndarray, ground: GroundModel | None = None, config: StemDetectionConfig | None = None, *, name: str = '', heights: ndarray | None = None) -> StemMap

Detect the stems of one scan.

Horizontal slices through the stem band are clustered and fitted with RANSAC circles, the circles are linked up the stem, and each chain gives an axis, a diameter at breast height and a quality.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 3) points in the scan's own (levelled) frame.	required
ground	GroundModel | None	Terrain model; fitted if None.	None
config	StemDetectionConfig | None	Detector settings.	None
name	str	Name of the stem map.	''
heights	ndarray | None	Heights above ground of points, if already known.	None

Returns:

Type	Description
StemMap	Stems, best first, with z on the scan's own terrain.

stem_map_from_arrays

stem_map_from_arrays(xy: ndarray, dbh=None, z=None, name: str = '') -> StemMap

A stem map from plain arrays (tests and external stem lists).

Parameters:

Name	Type	Description	Default
xy	ndarray	(n, 2) positions.	required
dbh		Diameters (m); 0.3 if None.	None
z		Elevations; 0 if None.	None
name	str	Name of the map.	''

Returns:

Type	Description
StemMap

sylva.coreg.matching

Global coarse registration by stem-map matching.

Aligning two stem maps is finding the rigid motion that best superimposes two partly overlapping 2-D point patterns whose points carry a diameter. It is searched exhaustively over pairs of stems, not sampled: the distance between two stems does not depend on the unknown transform, so a pair in one scan can only correspond to a pair at the same separation in the other, and a sorted pair table makes that a binary search. Two correspondences give yaw and translation, the right model for levelled scans; ICP recovers any residual tilt. The search runs in the Rust core.

MatchConfig dataclass

Settings of match_stem_maps.

min_pair_distance class-attribute instance-attribute

min_pair_distance: float = 2.0

Stem pairs closer than this (m) are ambiguous and ignored.

max_pair_distance class-attribute instance-attribute

max_pair_distance: float = 35.0

Stem pairs further apart than this (m) are rarely seen together.

inlier_tolerance class-attribute instance-attribute

inlier_tolerance: float = 0.4

Horizontal distance (m) at which two stems count as the same tree.

max_stems class-attribute instance-attribute

max_stems: int = 70

The best stems of each scan taken into matching.

early_exit_inliers class-attribute instance-attribute

early_exit_inliers: int = 40

Stop once a hypothesis has this many inliers. High, so that the search normally runs its budget: in a planted lattice the first convincing hypothesis is often a shifted copy, and only a fuller search shows the true one and how close the two are.

distinct_yaw_deg class-attribute instance-attribute

distinct_yaw_deg: float = 5.0

A rival must differ by distinct_translation m or this much yaw.

MatchResult dataclass

Outcome of match_stem_maps; transform maps source into target.

Attributes:

Name	Type	Description
correspondences	ndarray	(n, 2) indices (source, target) of matched stems, into the stem maps as given.
success	bool	At least min_inliers stems matched consistently. Check it first.
ambiguity	float	Inliers of the best distinctly different alignment over the best's. Near 1 the pattern matches itself elsewhere (a planted lattice) and the winner is a coin toss.
rival	MatchResult | None	That alternative, refined, for ICP to decide between.

inlier_fraction property

inlier_fraction: float

Inliers over the smaller stem map.

match_stem_maps

match_stem_maps(source: StemMap, target: StemMap, config: MatchConfig | None = None) -> MatchResult

Yaw and translation aligning source stems onto target stems, with no initial guess.

Parameters:

Name	Type	Description	Default
source	StemMap	Stem maps.	required
target	StemMap	Stem maps.	required
config	MatchConfig | None	Matcher settings.	None

Returns:

Type	Description
MatchResult	Check success before trusting transform.

sylva.coreg.icp

Fine registration by iterative closest point, tuned for forests.

Two departures from a textbook ICP (Besl & McKay 1992) matter:

• Point-to-plane with a planarity gate. The point-to-plane cost of Chen & Medioni (1992) is kept to planar neighbourhoods (stems, ground, logs), the parts stable between viewpoints: foliage is a view-dependent mess whose normals are noise.
• Robust, trimmed correspondences. Forest scans overlap only partly, so many source points have no true match. A Huber (1964) weight with an adaptive scale and a trim (as in the trimmed ICP of Chetverikov et al. 2002) phased in over the first iterations of each level keep them from biasing the pose without discarding the far points that carry the rotation.

The solver is a damped Gauss-Newton on SE(3) with a capped step, over a coarse-to-fine voxel pyramid. The loop runs in the Rust core.

ICPConfig dataclass

Settings of icp.

The pyramid should span only the error left by the coarse stem match, a few centimetres to decimetres: a 0.4 m level with a 2 m search blurs the stand into a blob and converges confidently to the wrong pose.

max_distances class-attribute instance-attribute

max_distances: tuple[float, ...] | None = (0.8, 0.4, 0.2, 0.12)

Correspondence cut-off per level (m); None for 5 x voxel.

method class-attribute instance-attribute

method: str = 'point_to_plane'

"point_to_plane" or "point_to_point".

robust class-attribute instance-attribute

robust: str = 'huber'

"huber", "tukey" or "none".

robust_scale class-attribute instance-attribute

robust_scale: float = 0.05

Lower bound (m) of the adaptive robust scale.

trim_fraction class-attribute instance-attribute

trim_fraction: float = 0.85

Keep this fraction of the closest correspondences once trimming is in.

trim_ramp class-attribute instance-attribute

trim_ramp: int = 3

Iterations over which trimming is phased in, per level.

min_planarity class-attribute instance-attribute

min_planarity: float = 0.25

Planarity a target point needs (point-to-plane); 0 disables the gate.

fitness_threshold class-attribute instance-attribute

fitness_threshold: float = 0.1

Distance (m) within which a source point counts as fitted.

max_points class-attribute instance-attribute

max_points: int = 120000

Cap per level after thinning.

distances

distances() -> tuple[float, ...]

Correspondence cut-off of each level.

ICPResult dataclass

Outcome of icp; transform maps source into target, initial guess included.

Attributes:

Name	Type	Description
fitness	float	Fraction of (voxel-thinned) source points with a target point within fitness_threshold.
inlier_rmse	float	RMSE over those points (m).
information	PlaneInformation | None	How well the point-to-plane correspondences pin the final transform; None for point-to-point ICP or too little overlap.

PlaneInformation dataclass

How well point-to-plane correspondences pin a transform.

Attributes:

Name	Type	Description
hessian	ndarray	(6, 6) robust-weighted Gauss-Newton matrix sum w a a^T of the residuals n . (T p - q), for an update exp(xi) @ T with xi = [omega, v] in the target frame. Its weak directions are the ones the surfaces leave free: flat ground fixes height, roll and pitch but lets the scans slide horizontally.
sigma	float	Weighted RMS point-to-plane residual (m).
n	int	Correspondences.

ICPTarget

A target's ICP pyramid, built once and reused.

Every ICP rebuilds the target's voxel pyramid, normals and search trees, over half its time on a large scan. The pyramid depends only on the target and the pyramid settings, so registering several scans against a prepared target gives exactly the results of passing the points.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 3) target points.	required
config	ICPConfig | None	ICP settings; the pyramid is built for its voxel_sizes, max_points, method, normal_neighbours and seed, and can only be used with settings that agree on those.	None

icp

icp(source: ndarray, target: ndarray | ICPTarget, initial: ndarray | None = None, config: ICPConfig | None = None) -> ICPResult

Align source onto target.

Parameters:

Name	Type	Description	Default
source	ndarray	(n, 3) points in their own frame.	required
target	ndarray | ICPTarget	(m, 3) points in their own frame, or an ICPTarget prepared with the same pyramid settings.	required
initial	ndarray | None	Initial guess of the source-to-target transform. Forest scans need a good one, from sylva.coreg.match_stem_maps.	None
config	ICPConfig | None	Settings; the defaults suit plot-scale forest scans.	None

Returns:

Type	Description
ICPResult

plane_information

plane_information(source: ndarray, target: ndarray | ICPTarget, transform: ndarray, config: ICPConfig | None = None) -> PlaneInformation | None

Point-to-plane information of a transform, as at the end of icp.

The correspondences, planarity gate, trim and robust weights of the finest pyramid level, without a step: for transforms ICP did not produce (a coarse match kept over ICP's, a placement).

Parameters:

Name	Type	Description	Default
source	ndarray	As for icp.	required
target	ndarray	As for icp.	required
transform	ndarray	Source-to-target transform.	required
config	ICPConfig | None	ICP settings.	None

Returns:

Type	Description
PlaneInformation or None	None if fewer than 10 correspondences survive.

evaluate_registration

evaluate_registration(source: ndarray, target: ndarray, transform: ndarray, threshold: float = 0.1, *, max_points: int = 200000, voxel: float | None = 0.05, seed: int = 0) -> tuple[float, float, int]

Score a registration without changing it.

Fitness, the fraction of source points with a target point within threshold, says more than RMSE in a forest: a low RMSE over a handful of points usually means the clouds slid into a spurious minimum.

Parameters:

Name	Type	Description	Default
source	ndarray	(n, 3) points.	required
target	ndarray	(n, 3) points.	required
transform	ndarray	Source-to-target transform.	required
threshold	float	Inlier distance (m).	0.1
max_points	int	Random cap on each cloud after thinning.	200000
voxel	float | None	Voxel centroids at this size first (m); None to skip.	0.05
seed	int	Seed of the random cap.	0

Returns:

Name	Type	Description
fitness	float
inlier_rmse	float
n_inliers	int

sylva.coreg.ground

The terrain model coregistration measures heights from.

Stem detection works on height above ground, not raw z: on sloping terrain a fixed elevation band cuts through the canopy on one side of a plot and the litter on the other. This is a raster DTM: a low percentile of z per cell, deep pits rejected, gaps filled from the nearest observed cell, a grey opening and smoothing, then a slope limit so a shrub or a log cannot lift the surface. It is deliberately not sylva.ground's CSF: the registration was tuned on the stems and heights this model gives.

GroundModel dataclass

A raster terrain model with bilinear interpolation.

Attributes:

Name	Type	Description
elevation	ndarray	(ny, nx) terrain height, gap-free.
origin	ndarray	(x0, y0), taken as the centre of cell [0, 0] when sampling (a convention a quarter-cell off the binning, kept because the registration was tuned on the heights it gives).
cell_size	float	Grid spacing (m).
observed	ndarray	Cells backed by ground returns rather than filled.

slope_deg property

slope_deg: float

Mean terrain slope (degrees), a sanity check on the fit.

height_at

height_at(xy: ndarray) -> np.ndarray

Terrain height at (n, 2) locations (clamped at the edges).

support

support(xy: ndarray) -> np.ndarray

Were (n, 2) locations backed by ground returns?

normalise

normalise(points: ndarray, dtype=np.float64) -> np.ndarray

Heights above ground of (n, 3) points.

fit_ground

fit_ground(points: ndarray, cell_size: float = 0.5, *, percentile: float = 5.0, max_slope: float = 1.0, smooth_cells: int = 3, opening_cells: int = 5, min_points_per_cell: int = 1, pit_depth: float = 3.0, pit_window: int = 9, max_points: int | None = 8000000, seed: int = 0) -> GroundModel

Fit a raster terrain model to a point cloud.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 3) points, z up.	required
cell_size	float	Resolution (m).	0.5
percentile	float	Per-cell z percentile taken as ground; robust to the few below-ground multipath returns.	5.0
max_slope	float	Largest rise over run between neighbouring cells.	1.0
smooth_cells	int	Windows (cells) of the smoothing and of the grey opening that removes low objects such as logs.	3
opening_cells	int	Windows (cells) of the smoothing and of the grey opening that removes low objects such as logs.	3
min_points_per_cell	int	Points a cell needs to count as observed.	1
pit_depth	float	A cell more than pit_depth m below the median of the pit_window cells around it is a pit, not ground. Tilted VZ-400i scans carry echoes 20-30 m under the terrain, and one such cell otherwise drags the surface down through the slope limit.	3.0
pit_window	float	A cell more than pit_depth m below the median of the pit_window cells around it is a pit, not ground. Tilted VZ-400i scans carry echoes 20-30 m under the terrain, and one such cell otherwise drags the surface down through the slope limit.	3.0
max_points	int | None	Fit from at most this many points, drawn at random (None: all).	8000000
seed	int	Seed of that draw.	0

Returns:

Type	Description
GroundModel

Raises:

Type	Description
ValueError	On an empty cloud, a non-positive cell size or no observed cell.

sylva.coreg.geometry

Resampling, local geometry and nearest-neighbour search for coregistration.

Scans arrive with tens of millions of points at wildly varying density (dense near the scanner, sparse at range), which slows ICP down and biases it towards whichever region is densest; everything here works on voxel centroids.

KdTree

Nearest-neighbour search over a fixed set of 2-D or 3-D points.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 2) or (n, 3) points.	required

query

query(queries: ndarray, distance_upper_bound: float = np.inf) -> tuple[np.ndarray, np.ndarray]

Nearest point to each query.

Parameters:

Name	Type	Description	Default
queries	ndarray	(q, d) points of the tree's dimension.	required
distance_upper_bound	float	Only neighbours within this distance count.	inf

Returns:

Name	Type	Description
distance	ndarray	(q,), inf where there is none.
index	ndarray	(q,), n where there is none (as scipy's cKDTree).

voxel_downsample

voxel_downsample(points: ndarray, voxel: float, *, centroid: bool = True, return_counts: bool = False)

One point per voxel.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 3) points.	required
voxel	float	Voxel size (m), on a grid anchored at floor(p / voxel).	required
centroid	bool	Average the points of each voxel; otherwise keep the first.	True
return_counts	bool	Also return how many points fed each output point.	False

Returns:

Type	Description
ndarray or (ndarray, ndarray)	(m, 3) points, ordered by voxel, and optionally (m,) counts.

estimate_normals

estimate_normals(points: ndarray, k: int = 20, radius: float | None = None) -> tuple[np.ndarray, np.ndarray]

Unit normals and planarity by local PCA.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 3) points.	required
k	int	Neighbours, the point itself included.	20
radius	float | None	Neighbours further than this (m) are left out; a point with fewer than 3 within it gets a zero normal and planarity 0.	None

Returns:

Name	Type	Description
normals	ndarray	(n, 3).
planarity	ndarray	(l1 - l0) / l2 in [0, 1]: how plane-like each neighbourhood is.

planar_filter

planar_filter(points: ndarray, min_planarity: float = 0.35, voxel: float | None = 0.05, k: int = 20, radius: float | None = 0.15) -> np.ndarray

Keep the locally planar points: stems, ground and logs.

The most useful step before ICP in a forest. Crown returns have no stable geometry between viewpoints (occlusion, leaves, wind) and make thousands of plausible but wrong correspondences; filtering both clouds leaves the stable skeleton, vertical stem surfaces plus the terrain that fixes the height.

Parameters:

Name	Type	Description	Default
points	ndarray	(n, 3) points.	required
min_planarity	float	Planarity a point needs.	0.35
voxel	float | None	Voxel centroids at this size first (m); None to skip.	0.05
k	int	Neighbourhood of the planarity.	20
radius	int	Neighbourhood of the planarity.	20

Returns:

Type	Description
ndarray	(m, 3) planar points.

sylva.coreg.reflectors

Retro-reflective targets: the strongest registration feature, where present.

A target is a point, not a fitted cylinder, so it is located to millimetres, and three shared targets fix all six degrees of freedom where a stem map fixes four. Where scans saw targets they are tried before stems.

Targets come from what the scanner or RiSCAN PRO already found (a RIEGL .tpl tie-point list beside each scan, or a RiSCAN .rfl reflector list), or are detected in the cloud by their return strength: retro-reflectors return far more energy than anything natural. Reading, detection and matching run in the Rust core.

Reflector dataclass

One retro-reflective target, in the scan's own frame.

ReflectorMatch dataclass

Alignment of two target sets; transform maps source into target.

read_tiepoint_list

read_tiepoint_list(path: str | Path) -> list[Reflector]

Read a RIEGL .tpl tie-point list (JSON).

Parameters:

Name	Type	Description	Default
path	str | Path	The .tpl file.	required

Returns:

Type	Description
list of Reflector	Empty for a missing or unreadable file: a position that found no targets is normal.

read_reflector_list

read_reflector_list(path: str | Path) -> list[Reflector]

Read a RiSCAN PRO .rfl reflector list (RieglRflID).

Each ReflectorN= line holds the fields named by the ReflectorIdx= line (name, index, status, x, y, z, ..., reflectance, diameter, points), with x, y, z in the scanner's own frame. The position-level ScanPosNNN.rfl holds the fine-scanned targets; a scan's own <timestamp>.rfl those found while scanning.

Parameters:

Name	Type	Description	Default
path	str | Path	The .rfl file.	required

Returns:

Type	Description
list of Reflector	Empty for a missing or unreadable file.

detect_reflectors

detect_reflectors(xyz: ndarray, reflectance: ndarray | None, *, min_reflectance: float = 5.0, cluster_radius: float = 0.15, min_points: int = 8, max_extent: float = 0.5) -> list[Reflector]

Find retro-reflective targets in a scan by their return strength.

min_reflectance is the setting that matters and it is not universal: in a scan with targets the reflectance histogram is bimodal, with a wide gap between the scene and the targets; check it first.

Parameters:

Name	Type	Description	Default
xyz	ndarray	(n, 3) points.	required
reflectance	ndarray | None	(n,) reflectance (dB).	required
min_reflectance	float	Returns at least this bright are target candidates.	5.0
cluster_radius	float	Single-link clustering distance (m).	0.15
min_points	int	Returns a target needs.	8
max_extent	float	Clusters larger than this (m) are not targets (a wet surface, a sign).	0.5

Returns:

Type	Description
list of Reflector

match_reflectors

match_reflectors(source: list[Reflector], target: list[Reflector], *, tolerance: float = 0.05, min_inliers: int = 3, distance_tolerance: float = 0.03) -> ReflectorMatch

Align two target sets with no initial guess.

Triangles are searched exhaustively (targets are few), and three correspondences give a full 6-DoF transform. distance_tolerance is tight because target positions are good to millimetres, which is what makes the match unambiguous.

Parameters:

Name	Type	Description	Default
source	list[Reflector]	Targets of the two scans.	required
target	list[Reflector]	Targets of the two scans.	required
tolerance	float	How far (m) a matched target may sit from its partner.	0.05
min_inliers	int	Targets that must correspond.	3
distance_tolerance	float	Agreement (m) of triangle sides.	0.03

Returns:

Type	Description
ReflectorMatch	Its correspondences are ordered by source index.

sylva.coreg.posegraph

Multi-scan global optimisation.

Pairwise registration leaves a survey inconsistent: going A->B->C->A does not return to the start, and the drift shows up as doubled stems. A pose graph finds the poses that best explain every pairwise measurement at once, spreading the closure error over the loop (Lu & Milios 1997).

Nodes are world_from_scan poses; edges are measured relative transforms with an information matrix saying how far, and in which directions, each is trusted: an ICP edge's comes from its point-to-plane correspondences (plane_edge_information), so a pair matched mostly on flat ground holds height, roll and pitch firmly and the horizontal position loosely. The solver is Levenberg-Marquardt (Levenberg 1944; Marquardt 1963) on SE(3) with a Huber (1964) kernel and an explicit outlier pass, because in a forest a pairwise match can be confidently and completely wrong when two parts of a stand have similar stem patterns. The graph here is a container; the solve runs in the Rust core.

PoseGraphEdge dataclass

A measured relative transform.

transform maps scan i into scan j (j_from_i).

weight property

weight: float

Scalar trust of the edge, for the spanning-tree initialisation.

OptimisationResult dataclass

Diagnostics of PoseGraph.optimise.

PoseGraph

A pose graph over scan poses.

poses[i] is world_from_scan_i. The world frame is that of the reference node, held fixed. Other nodes can be held fixed too, at poses given in fixed (scans already trusted, which new scans are registered into).

Parameters:

Name	Type	Description	Default
n_nodes	int	Number of scans.	required
reference	int	Node whose pose defines the world frame.	0
fixed	dict[int, ndarray] | None	{node: world_from_scan} of further nodes held fixed.	None

add_edge

add_edge(i: int, j: int, transform: ndarray, *, information: ndarray | None = None, fitness: float = 1.0, rmse: float = 0.01, n_correspondences: int = 0, label: str = '') -> PoseGraphEdge

Add a measurement that scan i maps into scan j by transform.

Returns:

Type	Description
PoseGraphEdge

components

components() -> list[list[int]]

Connected components, as sorted lists of nodes.

initialise

initialise(reference: int | None = None) -> list[np.ndarray]

Initialise the poses by walking a maximum-weight spanning tree.

Chaining the strongest edges first starts the solve close to the optimum, which matters because the SE(3) cost is not convex. Nodes in components without an anchor keep the identity.

Returns:

Type	Description
list of numpy.ndarray

residual

residual(edge: PoseGraphEdge, poses: list[ndarray] | None = None) -> np.ndarray

6-vector error of one edge under poses (default: the current ones).

total_error

total_error(poses: list[ndarray] | None = None, edges: list[PoseGraphEdge] | None = None) -> float

Sum of squared Mahalanobis edge errors over edges (default: all).

optimise

optimise(max_iterations: int = 200, *, tolerance: float = 1e-06, huber_delta: float = 3.0, reject_outliers: bool = True, outlier_sigma: float = 5.0, max_rejection_passes: int = 2) -> OptimisationResult

Levenberg-Marquardt over all edges, with outlier rejection.

Parameters:

Name	Type	Description	Default
max_iterations	int	Cap per rejection pass; a 62-scan survey needed more than 60.	200
tolerance	float	Stop when an iteration improves the error by less than this fraction of it.	1e-06
huber_delta	float	Mahalanobis distance beyond which an edge is down-weighted.	3.0
reject_outliers	bool	Remove edges whose error stays far above the median, then solve again. An edge is never removed if that would cut a node off from the anchors it reached.	True
outlier_sigma	float	How many robust standard deviations above the median is an outlier.	5.0
max_rejection_passes	int	Solve-and-reject passes.	2

Returns:

Type	Description
OptimisationResult

relative

relative(i: int, j: int) -> np.ndarray

Optimised transform mapping scan i into scan j's frame.

default_information

default_information(rmse: float, fitness: float, n_correspondences: int, extent: float = 15.0) -> np.ndarray

Diagonal information matrix from registration quality.

Translation precision scales as the residual over the square root of the correspondence count; rotation precision is that positional uncertainty spread over the scan's extent (m).

Parameters:

Name	Type	Description	Default
rmse	float	ICP quality of the edge.	required
fitness	float	ICP quality of the edge.	required
n_correspondences	float	ICP quality of the edge.	required
extent	float	Spatial extent of a scan (m).	15.0

Returns:

Type	Description
ndarray	(6, 6) information, rotation block first.

adjoint

adjoint(T: ndarray) -> np.ndarray

(6, 6) adjoint of a transform on [omega, v] twists: T @ se3_exp(xi) @ invert(T) == se3_exp(adjoint(T) @ xi).

plane_edge_information

plane_edge_information(hessian: ndarray, sigma: float, n: int, transform: ndarray, *, patch_points: float = 100.0, min_sigma: float = 0.005) -> np.ndarray

Edge information from the point-to-plane correspondences of an ICP.

The Gauss-Newton matrix of the correspondences (~sylva.coreg.icp.PlaneInformation) says in which directions the surfaces pin the transform; it is scaled to the residual and carried from the target frame, where ICP perturbs the transform, into the frame of the edge residual. Neighbouring residuals share the error of the surface they sample, so n correspondences count as n / patch_points independent ones: the information of an edge grows with its overlap, but not to the point where millimetres of disagreement between two edges look like outliers.

Parameters:

Name	Type	Description	Default
hessian	ndarray	(6, 6) sum w a a^T for updates se3_exp(xi) @ transform.	required
sigma	float	Weighted RMS point-to-plane residual (m); floored at min_sigma.	required
n	int	Correspondences behind hessian.	required
transform	ndarray	The measured target_from_source transform.	required
patch_points	float	Correspondences per independent observation.	100.0
min_sigma	float	Floor on sigma (m): a residual of a millimetre says the surfaces agree, not that the pose is known to a millimetre.	0.005

Returns:

Type	Description
ndarray	(6, 6) information of the edge residual se3_log(invert(transform) @ target_from_source), rotation block first.

sylva.coreg.refine

Joint multi-view refinement: every pose solved at once from raw correspondences.

The pose graph sees each accepted pair as one relative transform, which throws away the geometry behind it. This keeps every point-to-plane correspondence of every accepted pair, adds the matched stems as horizontal constraints, and solves all poses together with a robust Gauss-Newton, the reference fixed, so that errors are spread by the data rather than handed from scan to scan. The solve runs in the Rust core.

JointRefinement

Result of refine_joint: the poses and how far each moved.

Attributes:

Name	Type	Description
poses		Refined world_from_scan.
shifts, rotations		How far each scan moved (m, degrees).
residual_before, residual_after		Median absolute point-to-plane residual (m).
correspondences		Largest number of correspondences in one association.

refine_joint

refine_joint(points: list[ndarray], poses: list[ndarray], edges: list[tuple[int, int]], stems: list[ndarray], reference: int, *, voxel_sizes: tuple[float, ...] = (0.1, 0.05, 0.03), max_distances: tuple[float, ...] = (0.3, 0.15, 0.08), rounds: int = 3, iterations: int = 3, points_per_scan: int = 400000, correspondences_per_pair: int = 20000, min_planarity: float = 0.25, normal_neighbours: int = 20, max_normal_angle_deg: float = 45.0, stem_weight: float = 0.05, stem_radius: float = 0.15, stem_scale: float = 0.03, min_voxel_points: int = 1, robust_scale: float = 0.02, prior_translation: float = 0.03, prior_rotation_deg: float = 0.3, max_step_translation: float = 0.05, max_step_rotation_deg: float = 0.5, seed: int = 0, log: Callable[[str], None] | None = None) -> JointRefinement

Solve all poses together from point-to-plane and stem correspondences.

Levels run coarse to fine. At each, correspondences are associated rounds times under the current poses and each association is followed by iterations Gauss-Newton steps with Huber (1964) reweighting.

Stems are paired afresh at every association, mutually nearest within stem_radius horizontally, with a fixed Huber scale stem_scale, and carry stem_weight of each pair's weight whatever their count, so a few trunks can steer the horizontal directions that a hundred thousand ground points leave open.

Two things keep the solve honest: a weak prior ties every pose to its pose-graph solution (prior_translation m, prior_rotation_deg at one sigma), so a scan whose correspondences are all ground cannot slide; and a step is kept only if it lowers the robust cost, with the damping raised and the step retried otherwise, capped at max_step_translation and max_step_rotation_deg.

Parameters:

Name	Type	Description	Default
points	list[ndarray]	Per scan, its planar ICP points in its own frame.	required
poses	list[ndarray]	Current world_from_scan.	required
edges	list[tuple[int, int]]	Accepted pairs.	required
stems	list[ndarray]	Per scan, its stem positions in its own frame.	required
reference	int	Scan held fixed.	required
log	Callable[[str], None] | None	Called with progress messages.	None

Returns:

Type	Description
JointRefinement

sylva.coreg.simulate

Synthetic plots and virtual scans with known poses, for testing registration.

Field data rarely comes with ground truth. These plots have known stems, and the scans reproduce what breaks forest registration: point density and noise that fall off with range, incidence drop-out, and above all mutual occlusion, which is why two scans of one plot may share only a fraction of their stems.

Tree dataclass

A leaning, tapering stem carrying a crown.

radius_at

radius_at(h)

Stem radius at height h above the base.

axis_point

axis_point(h) -> np.ndarray

Axis position at height h above the base.

Plot dataclass

A plot on an analytic terrain surface.

stem_positions property

stem_positions: ndarray

True stem axis positions 1.3 m above the base.

terrain

terrain(xy: ndarray) -> np.ndarray

Terrain elevation at (n, 2) locations.

Survey dataclass

Simulated scans, each in its own scanner frame, with their true poses.

Attributes:

Name	Type	Description
true_transforms	list[ndarray]	world_from_scan of each scan: what a registration should recover.

overlap_stems

overlap_stems(i: int, j: int) -> int

Trees detectable in both scans: an upper bound on matches.

simulate_plot

simulate_plot(size: float = 40.0, n_trees: int = 45, *, seed: int = 0, dbh_mean: float = 0.28, dbh_sigma: float = 0.45, min_spacing: float = 1.8, slope: float = 0.06, relief: float = 1.2, max_lean_deg: float = 6.0) -> Plot

A plot of n_trees on undulating terrain.

Positions by dart throwing with a minimum spacing, DBH log-normal.

Returns:

Type	Description
Plot

scan_plot

scan_plot(plot: Plot, scanner_xy: ndarray, *, seed: int = 0, max_range: float = 30.0, scanner_height: float = 1.6, angular_step: float = 0.0012, range_noise: float = 0.004, ground_density: float = 150.0, foliage_density: float = 2500.0, understorey_density: float = 70.0, occlusion: bool = True, world_frame: bool = True, min_stem_returns: int = 150) -> tuple[PointCloud, np.ndarray, np.ndarray]

A virtual scan of plot from scanner_xy.

Returns:

Name	Type	Description
cloud	PointCloud	With label (tree index, -1 ground, -2 foliage, -3 understorey) and range attributes.
visible	ndarray	Trees that kept at least min_stem_returns stem points after occlusion: the fair recall denominator, since a stem grazed by a few beams cannot be fitted across several slices.
scanner	ndarray	(3,) scanner position in the plot frame.

simulate_survey

simulate_survey(n_scans: int = 4, *, plot: Plot | None = None, seed: int = 0, scan_radius: float = 12.0, max_range: float = 30.0, levelled: bool = True, tilt_deg: float = 0.5, **scan_kwargs) -> Survey

A multi-scan survey of one plot.

Each cloud is in its own scanner frame with an unknown yaw, as a levelled but unoriented TLS survey. levelled=False adds a little roll and pitch.

Returns:

Type	Description
Survey

sylva.coreg.transforms

Rigid transforms (SO(3) and SE(3)) used throughout coregistration.

Transforms are 4x4 homogeneous matrices in column-vector convention, q = T @ [p, 1]. Twists are ordered xi = [wx, wy, wz, tx, ty, tz], rotation first, the convention of the pose-graph solver. The maps and fits run in the Rust core.

identity

identity() -> np.ndarray

The 4x4 identity transform.

skew

skew(v: ndarray) -> np.ndarray

The 3x3 skew-symmetric matrix of a 3-vector.

so3_exp

so3_exp(w: ndarray) -> np.ndarray

Rotation matrix of a rotation vector.

so3_log

so3_log(R: ndarray) -> np.ndarray

Rotation vector of a rotation matrix.

se3_exp

se3_exp(xi: ndarray) -> np.ndarray

Transform of a twist [w, t].

Parameters:

Name	Type	Description	Default
xi	ndarray	(6,) twist, rotation first.	required

Returns:

Type	Description
ndarray	(4, 4) transform.

se3_log

se3_log(T: ndarray) -> np.ndarray

Twist [w, t] of a transform; the inverse of se3_exp.

invert

invert(T: ndarray) -> np.ndarray

Inverse of a rigid transform.

transform_points

transform_points(T: ndarray, points: ndarray) -> np.ndarray

Apply a transform to (n, 3) points.

transform_vectors

transform_vectors(T: ndarray, vectors: ndarray) -> np.ndarray

Rotate (n, 3) directions (the translation is ignored).

yaw_transform

yaw_transform(yaw: float, tx: float = 0.0, ty: float = 0.0, tz: float = 0.0) -> np.ndarray

A rotation of yaw radians about +z followed by a translation.

kabsch

kabsch(source: ndarray, target: ndarray, weights: ndarray | None = None) -> np.ndarray

Least-squares rigid transform mapping source onto target.

Parameters:

Name	Type	Description	Default
source	ndarray	(n, 3) corresponding points, n >= 3.	required
target	ndarray	(n, 3) corresponding points, n >= 3.	required
weights	ndarray | None	Optional (n,) weights.	None

Returns:

Type	Description
ndarray	(4, 4) proper rigid transform (reflections suppressed).

Raises:

Type	Description
ValueError	On mismatched shapes, fewer than 3 points or zero total weight.

kabsch_2d_yaw

kabsch_2d_yaw(source: ndarray, target: ndarray, weights: ndarray | None = None) -> np.ndarray

Least-squares yaw and translation (4 degrees of freedom) between point sets.

The right estimator for levelled scans, whose roll and pitch are known.

rotation_angle

rotation_angle(T: ndarray) -> float

Rotation magnitude of a transform, in radians (exact near zero).

transform_difference

transform_difference(A: ndarray, B: ndarray) -> tuple[float, float]

(rotation_rad, translation_m) of the discrepancy A^-1 B.