|
|
@@ -0,0 +1,2073 @@
|
|
|
+"""Detect a family of convex upper-left to lower-right quadratic curves.
|
|
|
+
|
|
|
+The detector starts from a high-recall curve mask. It splits the skeleton at
|
|
|
+junctions, lets local line fragments vote in a vertex-form parabola parameter
|
|
|
+space, scores complete hypotheses against the original pixels, and finally
|
|
|
+selects a mutually non-crossing family. No endpoint labels or language model
|
|
|
+predictions are used.
|
|
|
+"""
|
|
|
+
|
|
|
+from __future__ import annotations
|
|
|
+
|
|
|
+from dataclasses import asdict, dataclass, replace
|
|
|
+from itertools import combinations
|
|
|
+from contextvars import ContextVar
|
|
|
+from functools import wraps, lru_cache
|
|
|
+import math
|
|
|
+from typing import Any
|
|
|
+
|
|
|
+import cv2
|
|
|
+import numpy as np
|
|
|
+
|
|
|
+from .grid import PlotBounds
|
|
|
+from .models import Point
|
|
|
+
|
|
|
+
|
|
|
+_IMAGE_CACHE: ContextVar[dict | None] = ContextVar("curve_image_cache", default=None)
|
|
|
+
|
|
|
+
|
|
|
+def _image_cached(function):
|
|
|
+ """Reuse immutable image evidence only inside one extraction call."""
|
|
|
+ @wraps(function)
|
|
|
+ def call(*args, **kwargs):
|
|
|
+ cache = _IMAGE_CACHE.get()
|
|
|
+ if cache is None:
|
|
|
+ return function(*args, **kwargs)
|
|
|
+ def key(value):
|
|
|
+ if isinstance(value, (str, int, float, bool, type(None))):
|
|
|
+ return value
|
|
|
+ if isinstance(value, (tuple, list)):
|
|
|
+ return tuple(key(item) for item in value)
|
|
|
+ # Keep referenced objects alive; recycled object ids must not alias.
|
|
|
+ cache.setdefault("references", {})[id(value)] = value
|
|
|
+ return ("object", id(value))
|
|
|
+ token = (function.__name__, key(args), tuple((name, key(value)) for name, value in kwargs.items()))
|
|
|
+ if token not in cache:
|
|
|
+ cache[token] = function(*args, **kwargs)
|
|
|
+ return cache[token]
|
|
|
+ return call
|
|
|
+
|
|
|
+
|
|
|
+def _image_scope(function):
|
|
|
+ @wraps(function)
|
|
|
+ def call(*args, **kwargs):
|
|
|
+ token = _IMAGE_CACHE.set({})
|
|
|
+ try:
|
|
|
+ return function(*args, **kwargs)
|
|
|
+ finally:
|
|
|
+ _IMAGE_CACHE.reset(token)
|
|
|
+ return call
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class QuadraticCurveFamilyConfig:
|
|
|
+ """Scale-independent controls for quadratic-family extraction."""
|
|
|
+
|
|
|
+ junction_radius_ratio: float = 0.004
|
|
|
+ junction_radius_line_widths: float = 0.65
|
|
|
+ minimum_segment_span_ratio: float = 0.012
|
|
|
+ minimum_segment_line_widths: float = 1.75
|
|
|
+ sharp_turn_angle_degrees: float = 60.0
|
|
|
+ sharp_turn_simplification_ratio: float = 0.002
|
|
|
+ sharp_turn_simplification_line_widths: float = 1.25
|
|
|
+ sharp_turn_maximum_arm_line_widths: float = 3.5
|
|
|
+ sharp_turn_maximum_arm_fraction: float = 0.20
|
|
|
+ minimum_tangent_slope: float = 0.06
|
|
|
+ maximum_tangent_slope: float = 18.0
|
|
|
+ minimum_curvature: float = 0.30
|
|
|
+ maximum_curvature: float = 14.0
|
|
|
+ curvature_sample_count: int = 56
|
|
|
+ vertex_bin_size: float = 0.022
|
|
|
+ minimum_vote_segments: int = 1
|
|
|
+ maximum_seed_x_overlap_ratio: float = 0.50
|
|
|
+ maximum_seed_y_overlap_ratio: float = 0.50
|
|
|
+ segment_vote_beam_width: int = 64
|
|
|
+ # Gap is only a ranking scale, never a hard search radius.
|
|
|
+ association_gap_scale_ratio: float = 0.14
|
|
|
+ association_proposal_line_widths: float = 4.0
|
|
|
+ association_distance_weight: float = 0.08
|
|
|
+ association_curvature_weight: float = 1.20
|
|
|
+ association_tangent_absolute_tolerance: float = 0.75
|
|
|
+ association_tangent_relative_tolerance: float = 0.45
|
|
|
+ association_joint_median_line_widths: float = 2.5
|
|
|
+ association_joint_p80_line_widths: float = 4.0
|
|
|
+ association_joint_overlap_line_widths: float = 2.0
|
|
|
+ association_joint_minimum_overlap: float = 0.65
|
|
|
+ association_seed_drift_line_widths: float = 0.75
|
|
|
+ association_seed_overlap_drop: float = 0.10
|
|
|
+ association_joint_curvature_prior_weight: float = 0.002
|
|
|
+ association_path_beam_width: int = 32
|
|
|
+ association_path_branch_limit: int = 6
|
|
|
+ maximum_vote_hypotheses: int = 800
|
|
|
+ curve_sample_count: int = 260
|
|
|
+ support_distance_ratio: float = 0.0045
|
|
|
+ support_distance_line_widths: float = 1.0
|
|
|
+ direct_support_distance_ratio: float = 0.0025
|
|
|
+ direct_support_distance_line_widths: float = 0.55
|
|
|
+ minimum_x_span_ratio: float = 0.07
|
|
|
+ minimum_y_span_ratio: float = 0.24
|
|
|
+ minimum_pixel_support: float = 0.48
|
|
|
+ minimum_direct_support: float = 0.32
|
|
|
+ refinement_corridor_ratio: float = 0.025
|
|
|
+ refinement_corridor_line_widths: float = 2.5
|
|
|
+ column_half_window_ratio: float = 0.0015
|
|
|
+ column_half_window_line_widths: float = 0.4
|
|
|
+ maximum_column_band_line_widths: float = 5.0
|
|
|
+ maximum_gap_ratio: float = 0.17
|
|
|
+ maximum_candidate_count: int = 100
|
|
|
+ duplicate_distance_ratio: float = 0.014
|
|
|
+ minimum_family_separation_ratio: float = 0.010
|
|
|
+ # Reward evidence-backed members without assuming a fixed family size.
|
|
|
+ # Keeping this configurable makes the legacy selector reproducible while
|
|
|
+ # allowing dense/low-resolution charts to retain weaker outer members.
|
|
|
+ family_size_reward: float = 0.55
|
|
|
+ family_gap_irregularity_penalty: float = 1.40
|
|
|
+ # Optional semantic evidence. This is a requested family cardinality, not
|
|
|
+ # a filename-specific constant; candidates must still pass all geometric
|
|
|
+ # compatibility and pixel-support checks.
|
|
|
+ target_family_size: int | None = None
|
|
|
+ maximum_family_size: int = 15
|
|
|
+ family_beam_width: int = 900
|
|
|
+ endpoint_frontier_radius: float = 0.035
|
|
|
+ endpoint_frontier_tolerance: float = 0.095
|
|
|
+ endpoint_suppression_radius: float = 0.042
|
|
|
+ endpoint_locus_initial_tolerance: float = 0.080
|
|
|
+ endpoint_locus_final_tolerance: float = 0.070
|
|
|
+ endpoint_locus_minimum_curvature: float = -4.5
|
|
|
+ endpoint_locus_maximum_curvature: float = 0.5
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class ShortLineSegment:
|
|
|
+ """One junction-separated fragment used as quadratic voting evidence."""
|
|
|
+
|
|
|
+ segment_id: str
|
|
|
+ midpoint: Point
|
|
|
+ normalized_midpoint: Point
|
|
|
+ tangent_slope_normalized: float
|
|
|
+ normalized_x_range: tuple[float, float]
|
|
|
+ normalized_y_range: tuple[float, float]
|
|
|
+ length_px: float
|
|
|
+ span_px: float
|
|
|
+ fit_rmse_px: float
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class QuadraticCurve:
|
|
|
+ """One selected unlabeled quadratic characteristic curve."""
|
|
|
+
|
|
|
+ curve_id: str
|
|
|
+ coefficients_normalized: tuple[float, float, float]
|
|
|
+ vertex_normalized: Point
|
|
|
+ x_domain_normalized: tuple[float, float]
|
|
|
+ pixel_points: tuple[Point, ...]
|
|
|
+ # Segments that produced the original parameter-space vote.
|
|
|
+ supporting_segment_ids: tuple[str, ...]
|
|
|
+ # Segments associated again after the final global pixel refit.
|
|
|
+ associated_segment_ids: tuple[str, ...]
|
|
|
+ x_span_ratio: float
|
|
|
+ y_span_ratio: float
|
|
|
+ pixel_support: float
|
|
|
+ direct_pixel_support: float
|
|
|
+ maximum_gap_ratio: float
|
|
|
+ tangent_agreement: float
|
|
|
+ vote_score: float
|
|
|
+ global_score: float
|
|
|
+ confidence: float
|
|
|
+ proposed_segment_ids: tuple[str, ...] = ()
|
|
|
+ rejected_seed_ids: tuple[str, ...] = ()
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class QuadraticCurveFamilyResult:
|
|
|
+ """Serializable output of the complete model-free pipeline."""
|
|
|
+
|
|
|
+ coordinate_space: str
|
|
|
+ plot_bounds: PlotBounds
|
|
|
+ parameters: dict[str, Any]
|
|
|
+ adaptive_scale: dict[str, float | int]
|
|
|
+ skeleton_pixel_count: int
|
|
|
+ junction_count: int
|
|
|
+ segment_count: int
|
|
|
+ vote_hypothesis_count: int
|
|
|
+ scored_candidate_count: int
|
|
|
+ curves: tuple[QuadraticCurve, ...]
|
|
|
+ family_score: float
|
|
|
+
|
|
|
+ def to_dict(self) -> dict[str, Any]:
|
|
|
+ return asdict(self)
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class QuadraticCurveFamilyDiagnostics:
|
|
|
+ """Intermediate masks and labels used by the visualization script."""
|
|
|
+
|
|
|
+ skeleton: np.ndarray
|
|
|
+ junction_zone: np.ndarray
|
|
|
+ segment_labels: np.ndarray
|
|
|
+ segments: tuple[ShortLineSegment, ...]
|
|
|
+ adaptive_scale: dict[str, float | int]
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class _AdaptiveScale:
|
|
|
+ line_width_px: float
|
|
|
+ junction_radius_px: int
|
|
|
+ minimum_segment_span_px: float
|
|
|
+ minimum_segment_pixels: int
|
|
|
+ support_distance_px: float
|
|
|
+ direct_support_distance_px: float
|
|
|
+ refinement_corridor_px: float
|
|
|
+ column_half_window_px: int
|
|
|
+ maximum_column_band_px: int
|
|
|
+
|
|
|
+ def to_dict(self) -> dict[str, float | int]:
|
|
|
+ return asdict(self)
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class _VoteHypothesis:
|
|
|
+ curvature: float
|
|
|
+ vertex_x: float
|
|
|
+ vertex_y: float
|
|
|
+ vote_score: float
|
|
|
+ segment_ids: tuple[str, ...]
|
|
|
+ proposed_segment_ids: tuple[str, ...] = ()
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class _ScoredCandidate:
|
|
|
+ coefficients: np.ndarray
|
|
|
+ vertex: tuple[float, float]
|
|
|
+ domain: tuple[float, float]
|
|
|
+ points: np.ndarray
|
|
|
+ segment_ids: tuple[str, ...]
|
|
|
+ x_span: float
|
|
|
+ y_span: float
|
|
|
+ support: float
|
|
|
+ direct_support: float
|
|
|
+ maximum_gap: float
|
|
|
+ tangent_agreement: float
|
|
|
+ vote_score: float
|
|
|
+ score: float
|
|
|
+ family_position: float
|
|
|
+ proposed_segment_ids: tuple[str, ...] = ()
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class _SegmentFitEvidence:
|
|
|
+ median_distance_px: float
|
|
|
+ p80_distance_px: float
|
|
|
+ overlap_fraction: float
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class _SegmentGroupFit:
|
|
|
+ coefficients: np.ndarray
|
|
|
+ evidence: dict[str, _SegmentFitEvidence]
|
|
|
+ loss: float
|
|
|
+
|
|
|
+
|
|
|
+@dataclass(frozen=True)
|
|
|
+class _AssociationPath:
|
|
|
+ # Canonical x order includes internal insertions, not just two outer arms.
|
|
|
+ added: tuple[str, ...]
|
|
|
+ fit: _SegmentGroupFit
|
|
|
+
|
|
|
+
|
|
|
+def _segment_fit_evidence(
|
|
|
+ coefficients: np.ndarray,
|
|
|
+ points: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ overlap_band_px: float,
|
|
|
+) -> _SegmentFitEvidence:
|
|
|
+ """Measure all skeleton pixels against the untrimmed equation."""
|
|
|
+
|
|
|
+ width = max(bounds.right - bounds.left, 1.0)
|
|
|
+ height = max(bounds.bottom - bounds.top, 1.0)
|
|
|
+ pixel_slopes = np.polyval(np.polyder(coefficients), points[:, 0]) * height / width
|
|
|
+ distances = (
|
|
|
+ np.abs(points[:, 1] - np.polyval(coefficients, points[:, 0]))
|
|
|
+ * height / np.sqrt(1.0 + pixel_slopes**2)
|
|
|
+ )
|
|
|
+ return _SegmentFitEvidence(
|
|
|
+ median_distance_px=float(np.median(distances)),
|
|
|
+ p80_distance_px=float(np.percentile(distances, 80)),
|
|
|
+ overlap_fraction=float(np.mean(distances <= overlap_band_px)),
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _fit_segment_group(
|
|
|
+ segment_ids: tuple[str, ...],
|
|
|
+ normalized_pixels: dict[str, np.ndarray],
|
|
|
+ seed_ids: frozenset[str],
|
|
|
+ reference_coefficients: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ line_width_px: float,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> _SegmentGroupFit:
|
|
|
+ """Fit all fragments with balanced weights; never trim away a fragment.
|
|
|
+
|
|
|
+ Each seed receives twice the weight of a later fragment, independent of
|
|
|
+ pixel count. Iterative slope weights approximate normal-distance fitting.
|
|
|
+ A weak curvature prior stabilizes a single short seed's quadratic fit.
|
|
|
+ Every original skeleton point participates in the subsequent evidence
|
|
|
+ check, even when a fragment contains many pixels.
|
|
|
+ """
|
|
|
+
|
|
|
+ groups = [normalized_pixels[item] for item in segment_ids]
|
|
|
+ points = np.concatenate(groups)
|
|
|
+ base_weights = np.concatenate(
|
|
|
+ [
|
|
|
+ np.full(
|
|
|
+ len(group), np.sqrt((2.0 if item in seed_ids else 1.0) / len(group))
|
|
|
+ )
|
|
|
+ for item, group in zip(segment_ids, groups, strict=True)
|
|
|
+ ]
|
|
|
+ )
|
|
|
+ width = max(bounds.right - bounds.left, 1.0)
|
|
|
+ height = max(bounds.bottom - bounds.top, 1.0)
|
|
|
+ center_x = float(np.average(points[:, 0], weights=base_weights**2))
|
|
|
+ centered_x = points[:, 0] - center_x
|
|
|
+ design = np.column_stack((centered_x**2, centered_x, np.ones(len(points))))
|
|
|
+ coefficients = reference_coefficients.copy()
|
|
|
+ prior = config.association_joint_curvature_prior_weight
|
|
|
+ for _ in range(4):
|
|
|
+ pixel_slopes = (
|
|
|
+ np.polyval(np.polyder(coefficients), points[:, 0]) * height / width
|
|
|
+ )
|
|
|
+ weights = base_weights / np.sqrt(1.0 + pixel_slopes**2)
|
|
|
+ matrix = np.vstack((design * weights[:, None], [prior, 0.0, 0.0]))
|
|
|
+ target = np.append(points[:, 1] * weights, prior * reference_coefficients[0])
|
|
|
+ a, b, c = np.linalg.lstsq(matrix, target, rcond=None)[0]
|
|
|
+ # A short, almost straight fragment cannot determine curvature well.
|
|
|
+ # Solve at the convexity boundary rather than discarding a tight fit.
|
|
|
+ minimum_curvature = config.minimum_curvature * 0.35
|
|
|
+ if a < minimum_curvature:
|
|
|
+ a = minimum_curvature
|
|
|
+ b, c = np.linalg.lstsq(
|
|
|
+ design[:, 1:] * weights[:, None],
|
|
|
+ (points[:, 1] - a * centered_x**2) * weights,
|
|
|
+ rcond=None,
|
|
|
+ )[0]
|
|
|
+ coefficients = np.asarray(
|
|
|
+ (a, b - 2.0 * a * center_x, c - b * center_x + a * center_x**2)
|
|
|
+ )
|
|
|
+
|
|
|
+ evidence: dict[str, _SegmentFitEvidence] = {}
|
|
|
+ for item, group in zip(segment_ids, groups, strict=True):
|
|
|
+ evidence[item] = _segment_fit_evidence(
|
|
|
+ coefficients, group, bounds,
|
|
|
+ line_width_px * config.association_joint_overlap_line_widths,
|
|
|
+ )
|
|
|
+ loss = float(np.mean([
|
|
|
+ value.median_distance_px / max(line_width_px, 1e-9)
|
|
|
+ + 0.5 * value.p80_distance_px / max(line_width_px, 1e-9)
|
|
|
+ + 1.0 - value.overlap_fraction
|
|
|
+ for value in evidence.values()
|
|
|
+ ]))
|
|
|
+ return _SegmentGroupFit(coefficients, evidence, loss)
|
|
|
+
|
|
|
+
|
|
|
+def _segment_group_is_consistent(
|
|
|
+ trial: _SegmentGroupFit,
|
|
|
+ seed_baseline: _SegmentGroupFit,
|
|
|
+ line_width_px: float,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> bool:
|
|
|
+ """Require each added fragment to fit and protect every original seed."""
|
|
|
+
|
|
|
+ if (
|
|
|
+ not np.all(np.isfinite(trial.coefficients))
|
|
|
+ or trial.coefficients[0] < config.minimum_curvature * 0.35
|
|
|
+ ):
|
|
|
+ return False
|
|
|
+ for item, value in trial.evidence.items():
|
|
|
+ median_limit = line_width_px * config.association_joint_median_line_widths
|
|
|
+ p80_limit = line_width_px * config.association_joint_p80_line_widths
|
|
|
+ overlap_limit = config.association_joint_minimum_overlap
|
|
|
+ if item in seed_baseline.evidence:
|
|
|
+ baseline = seed_baseline.evidence[item]
|
|
|
+ drift = line_width_px * config.association_seed_drift_line_widths
|
|
|
+ median_limit = min(median_limit, baseline.median_distance_px + drift)
|
|
|
+ p80_limit = min(p80_limit, baseline.p80_distance_px + drift)
|
|
|
+ overlap_limit = max(
|
|
|
+ overlap_limit,
|
|
|
+ max(0.0, baseline.overlap_fraction - config.association_seed_overlap_drop),
|
|
|
+ )
|
|
|
+ if (
|
|
|
+ value.median_distance_px > median_limit
|
|
|
+ or value.p80_distance_px > p80_limit
|
|
|
+ or value.overlap_fraction < overlap_limit
|
|
|
+ ):
|
|
|
+ return False
|
|
|
+ return True
|
|
|
+
|
|
|
+
|
|
|
+def _thin_numpy(mask: np.ndarray) -> np.ndarray:
|
|
|
+ """Return a Zhang-Suen skeleton without the optional contrib module."""
|
|
|
+
|
|
|
+ skeleton = (mask > 0).astype(np.uint8)
|
|
|
+ for _ in range(120):
|
|
|
+ removed_total = 0
|
|
|
+ for phase in (0, 1):
|
|
|
+ padded = np.pad(skeleton, 1)
|
|
|
+ neighbors = (
|
|
|
+ padded[:-2, 1:-1],
|
|
|
+ padded[:-2, 2:],
|
|
|
+ padded[1:-1, 2:],
|
|
|
+ padded[2:, 2:],
|
|
|
+ padded[2:, 1:-1],
|
|
|
+ padded[2:, :-2],
|
|
|
+ padded[1:-1, :-2],
|
|
|
+ padded[:-2, :-2],
|
|
|
+ )
|
|
|
+ count = sum(neighbors)
|
|
|
+ transitions = sum(
|
|
|
+ (neighbors[index] == 0)
|
|
|
+ & (neighbors[(index + 1) % len(neighbors)] == 1)
|
|
|
+ for index in range(len(neighbors))
|
|
|
+ )
|
|
|
+ removable = (
|
|
|
+ (skeleton == 1)
|
|
|
+ & (count >= 2)
|
|
|
+ & (count <= 6)
|
|
|
+ & (transitions == 1)
|
|
|
+ )
|
|
|
+ north, east, south, west = (
|
|
|
+ neighbors[0],
|
|
|
+ neighbors[2],
|
|
|
+ neighbors[4],
|
|
|
+ neighbors[6],
|
|
|
+ )
|
|
|
+ if phase == 0:
|
|
|
+ removable &= (north * east * south == 0) & (
|
|
|
+ east * south * west == 0
|
|
|
+ )
|
|
|
+ else:
|
|
|
+ removable &= (north * east * west == 0) & (
|
|
|
+ north * south * west == 0
|
|
|
+ )
|
|
|
+ removed = int(np.count_nonzero(removable))
|
|
|
+ skeleton[removable] = 0
|
|
|
+ removed_total += removed
|
|
|
+ if removed_total == 0:
|
|
|
+ break
|
|
|
+ return skeleton * 255
|
|
|
+
|
|
|
+
|
|
|
+def _thin(mask: np.ndarray) -> np.ndarray:
|
|
|
+ """Return a one-pixel Zhang-Suen skeleton using OpenCV when available."""
|
|
|
+
|
|
|
+ binary = np.where(mask > 0, 255, 0).astype(np.uint8)
|
|
|
+ ximgproc = getattr(cv2, "ximgproc", None)
|
|
|
+ if ximgproc is not None and hasattr(ximgproc, "thinning"):
|
|
|
+ return ximgproc.thinning(
|
|
|
+ binary,
|
|
|
+ thinningType=ximgproc.THINNING_ZHANGSUEN,
|
|
|
+ )
|
|
|
+ return _thin_numpy(binary)
|
|
|
+
|
|
|
+
|
|
|
+def _crossing_number(skeleton: np.ndarray) -> np.ndarray:
|
|
|
+ foreground = (skeleton > 0).astype(np.uint8)
|
|
|
+ padded = np.pad(foreground, 1)
|
|
|
+ neighbors = (
|
|
|
+ padded[:-2, 1:-1],
|
|
|
+ padded[:-2, 2:],
|
|
|
+ padded[1:-1, 2:],
|
|
|
+ padded[2:, 2:],
|
|
|
+ padded[2:, 1:-1],
|
|
|
+ padded[2:, :-2],
|
|
|
+ padded[1:-1, :-2],
|
|
|
+ padded[:-2, :-2],
|
|
|
+ )
|
|
|
+ return sum(
|
|
|
+ (neighbors[index] == 0)
|
|
|
+ & (neighbors[(index + 1) % len(neighbors)] == 1)
|
|
|
+ for index in range(len(neighbors))
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _junction_seeds(skeleton: np.ndarray) -> np.ndarray:
|
|
|
+ """Include multi-pixel junctions invisible to a per-pixel crossing number."""
|
|
|
+ foreground = skeleton > 0
|
|
|
+ seeds = foreground & (_crossing_number(skeleton) >= 3)
|
|
|
+ squares = (foreground[:-1, :-1] & foreground[:-1, 1:]
|
|
|
+ & foreground[1:, :-1] & foreground[1:, 1:])
|
|
|
+ # Zhang–Suen can leave a 2x2 core at an X junction. Each core pixel
|
|
|
+ # individually has crossing number two, despite the four outgoing arms.
|
|
|
+ for dy in (0, 1):
|
|
|
+ for dx in (0, 1):
|
|
|
+ seeds[dy:dy + squares.shape[0], dx:dx + squares.shape[1]] |= squares
|
|
|
+ return seeds.astype(np.uint8) * 255
|
|
|
+
|
|
|
+
|
|
|
+def _estimate_line_width(
|
|
|
+ mask: np.ndarray,
|
|
|
+ skeleton: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+) -> float:
|
|
|
+ """Estimate the typical ink width at skeleton centerlines."""
|
|
|
+
|
|
|
+ left = max(0, math.floor(bounds.left))
|
|
|
+ right = min(mask.shape[1], math.ceil(bounds.right) + 1)
|
|
|
+ top = max(0, math.floor(bounds.top))
|
|
|
+ bottom = min(mask.shape[0], math.ceil(bounds.bottom) + 1)
|
|
|
+ plot_mask = mask[top:bottom, left:right]
|
|
|
+ plot_skeleton = skeleton[top:bottom, left:right]
|
|
|
+ inside_distance = cv2.distanceTransform(plot_mask, cv2.DIST_L2, 3)
|
|
|
+ center_distances = inside_distance[plot_skeleton > 0]
|
|
|
+ center_distances = center_distances[center_distances > 0]
|
|
|
+ if not len(center_distances):
|
|
|
+ return 1.0
|
|
|
+ diameter = 2.0 * float(np.median(center_distances))
|
|
|
+ maximum = max(1.0, min(right - left, bottom - top) * 0.03)
|
|
|
+ return float(np.clip(diameter, 1.0, maximum))
|
|
|
+
|
|
|
+
|
|
|
+def _derive_adaptive_scale(
|
|
|
+ mask: np.ndarray,
|
|
|
+ skeleton: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> _AdaptiveScale:
|
|
|
+ width = max(1.0, bounds.right - bounds.left)
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ minimum_dimension = min(width, height)
|
|
|
+ line_width = _estimate_line_width(mask, skeleton, bounds)
|
|
|
+ minimum_span = max(
|
|
|
+ minimum_dimension * config.minimum_segment_span_ratio,
|
|
|
+ line_width * config.minimum_segment_line_widths,
|
|
|
+ )
|
|
|
+ return _AdaptiveScale(
|
|
|
+ line_width_px=line_width,
|
|
|
+ junction_radius_px=max(
|
|
|
+ 1,
|
|
|
+ round(
|
|
|
+ max(
|
|
|
+ minimum_dimension * config.junction_radius_ratio,
|
|
|
+ line_width * config.junction_radius_line_widths,
|
|
|
+ )
|
|
|
+ ),
|
|
|
+ ),
|
|
|
+ minimum_segment_span_px=minimum_span,
|
|
|
+ minimum_segment_pixels=max(1, math.ceil(minimum_span * 0.65)),
|
|
|
+ support_distance_px=max(
|
|
|
+ width * config.support_distance_ratio,
|
|
|
+ line_width * config.support_distance_line_widths,
|
|
|
+ ),
|
|
|
+ direct_support_distance_px=max(
|
|
|
+ width * config.direct_support_distance_ratio,
|
|
|
+ line_width * config.direct_support_distance_line_widths,
|
|
|
+ ),
|
|
|
+ refinement_corridor_px=max(
|
|
|
+ height * config.refinement_corridor_ratio,
|
|
|
+ line_width * config.refinement_corridor_line_widths,
|
|
|
+ ),
|
|
|
+ column_half_window_px=max(
|
|
|
+ 1,
|
|
|
+ round(
|
|
|
+ max(
|
|
|
+ width * config.column_half_window_ratio,
|
|
|
+ line_width * config.column_half_window_line_widths,
|
|
|
+ )
|
|
|
+ ),
|
|
|
+ ),
|
|
|
+ maximum_column_band_px=max(
|
|
|
+ 1, math.ceil(line_width * config.maximum_column_band_line_widths)
|
|
|
+ ),
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _fit_segment(
|
|
|
+ segment_id: str,
|
|
|
+ xs: np.ndarray,
|
|
|
+ ys: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+) -> ShortLineSegment | None:
|
|
|
+ points = np.column_stack((xs, ys)).astype(float)
|
|
|
+ center = points.mean(axis=0)
|
|
|
+ centered = points - center
|
|
|
+ covariance = centered.T @ centered
|
|
|
+ _, eigenvectors = np.linalg.eigh(covariance)
|
|
|
+ direction = eigenvectors[:, -1]
|
|
|
+ if direction[0] < 0.0:
|
|
|
+ direction = -direction
|
|
|
+ projections = centered @ direction
|
|
|
+ fitted = center + np.outer(projections, direction)
|
|
|
+ residuals = np.linalg.norm(points - fitted, axis=1)
|
|
|
+ length = float(np.ptp(projections))
|
|
|
+ if length < 1e-6 or abs(float(direction[0])) < 0.045:
|
|
|
+ return None
|
|
|
+ width = max(1.0, bounds.right - bounds.left)
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ normalized_slope = float(direction[1] / direction[0] * width / height)
|
|
|
+ return ShortLineSegment(
|
|
|
+ segment_id=segment_id,
|
|
|
+ midpoint=Point(float(center[0]), float(center[1])),
|
|
|
+ normalized_midpoint=Point(
|
|
|
+ float((center[0] - bounds.left) / width),
|
|
|
+ float((center[1] - bounds.top) / height),
|
|
|
+ ),
|
|
|
+ tangent_slope_normalized=normalized_slope,
|
|
|
+ normalized_x_range=(
|
|
|
+ float((np.min(xs) - bounds.left) / width),
|
|
|
+ float((np.max(xs) - bounds.left) / width),
|
|
|
+ ),
|
|
|
+ normalized_y_range=(
|
|
|
+ float((np.min(ys) - bounds.top) / height),
|
|
|
+ float((np.max(ys) - bounds.top) / height),
|
|
|
+ ),
|
|
|
+ length_px=length,
|
|
|
+ span_px=float(max(np.ptp(xs), np.ptp(ys))),
|
|
|
+ fit_rmse_px=float(np.sqrt(np.mean(np.square(residuals)))),
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _trace_component_path(xs: np.ndarray, ys: np.ndarray) -> np.ndarray:
|
|
|
+ """Approximate the geodesic diameter of a one-pixel component."""
|
|
|
+
|
|
|
+ coordinates = [(int(x), int(y)) for x, y in zip(xs, ys)]
|
|
|
+ coordinate_set = set(coordinates)
|
|
|
+
|
|
|
+ def neighbors(point: tuple[int, int]) -> list[tuple[int, int]]:
|
|
|
+ x, y = point
|
|
|
+ return [
|
|
|
+ (x + dx, y + dy)
|
|
|
+ for dy in (-1, 0, 1)
|
|
|
+ for dx in (-1, 0, 1)
|
|
|
+ if (dx or dy) and (x + dx, y + dy) in coordinate_set
|
|
|
+ ]
|
|
|
+
|
|
|
+ def farthest_from(
|
|
|
+ start: tuple[int, int],
|
|
|
+ ) -> tuple[tuple[int, int], dict[tuple[int, int], tuple[int, int]]]:
|
|
|
+ queue = [start]
|
|
|
+ parent: dict[tuple[int, int], tuple[int, int]] = {}
|
|
|
+ distance = {start: 0}
|
|
|
+ for point in queue:
|
|
|
+ for neighbor in neighbors(point):
|
|
|
+ if neighbor in distance:
|
|
|
+ continue
|
|
|
+ distance[neighbor] = distance[point] + 1
|
|
|
+ parent[neighbor] = point
|
|
|
+ queue.append(neighbor)
|
|
|
+ farthest = max(distance, key=distance.get) # type: ignore[arg-type]
|
|
|
+ return farthest, parent
|
|
|
+
|
|
|
+ if not coordinates:
|
|
|
+ return np.empty((0, 2), dtype=int)
|
|
|
+ first, _ = farthest_from(coordinates[0])
|
|
|
+ second, parent = farthest_from(first)
|
|
|
+ path = [second]
|
|
|
+ while path[-1] != first:
|
|
|
+ path.append(parent[path[-1]])
|
|
|
+ path.reverse()
|
|
|
+ return np.asarray(path, dtype=int)
|
|
|
+
|
|
|
+
|
|
|
+def _rdp_indices(points: np.ndarray, epsilon: float) -> list[int]:
|
|
|
+ """Return retained indices from Ramer-Douglas-Peucker simplification."""
|
|
|
+
|
|
|
+ if len(points) <= 2:
|
|
|
+ return list(range(len(points)))
|
|
|
+ start = points[0].astype(float)
|
|
|
+ end = points[-1].astype(float)
|
|
|
+ vector = end - start
|
|
|
+ length = float(np.linalg.norm(vector))
|
|
|
+ offsets = points.astype(float) - start
|
|
|
+ if length <= 1e-9:
|
|
|
+ distances = np.linalg.norm(offsets, axis=1)
|
|
|
+ else:
|
|
|
+ distances = np.abs(
|
|
|
+ vector[0] * offsets[:, 1] - vector[1] * offsets[:, 0]
|
|
|
+ ) / length
|
|
|
+ split = int(np.argmax(distances))
|
|
|
+ if distances[split] <= epsilon:
|
|
|
+ return [0, len(points) - 1]
|
|
|
+ left = _rdp_indices(points[: split + 1], epsilon)
|
|
|
+ right = _rdp_indices(points[split:], epsilon)
|
|
|
+ return left[:-1] + [split + index for index in right]
|
|
|
+
|
|
|
+
|
|
|
+def _split_component_at_sharp_turns(
|
|
|
+ xs: np.ndarray,
|
|
|
+ ys: np.ndarray,
|
|
|
+ epsilon: float,
|
|
|
+ angle_threshold_degrees: float,
|
|
|
+ maximum_short_arm_px: float = math.inf,
|
|
|
+ maximum_short_arm_fraction: float = 1.0,
|
|
|
+) -> list[tuple[np.ndarray, np.ndarray]]:
|
|
|
+ """Split a skeleton path only at concentrated, non-smooth turns."""
|
|
|
+
|
|
|
+ path = _trace_component_path(xs, ys)
|
|
|
+ if len(path) < 5:
|
|
|
+ return [(xs, ys)]
|
|
|
+ simplified_indices = _rdp_indices(path, epsilon)
|
|
|
+ split_indices: list[int] = []
|
|
|
+ for previous_index, index, next_index in zip(
|
|
|
+ simplified_indices,
|
|
|
+ simplified_indices[1:],
|
|
|
+ simplified_indices[2:],
|
|
|
+ ):
|
|
|
+ incoming = path[index].astype(float) - path[previous_index]
|
|
|
+ outgoing = path[next_index].astype(float) - path[index]
|
|
|
+ denominator = float(np.linalg.norm(incoming) * np.linalg.norm(outgoing))
|
|
|
+ if denominator <= 1e-9:
|
|
|
+ continue
|
|
|
+ cosine = float(
|
|
|
+ np.clip(np.dot(incoming, outgoing) / denominator, -1.0, 1.0)
|
|
|
+ )
|
|
|
+ angle = math.degrees(math.acos(cosine))
|
|
|
+ short_arm = min(
|
|
|
+ float(np.linalg.norm(incoming)), float(np.linalg.norm(outgoing))
|
|
|
+ )
|
|
|
+ concentrated = (
|
|
|
+ short_arm <= maximum_short_arm_px
|
|
|
+ and short_arm / max(len(path) - 1, 1) <= maximum_short_arm_fraction
|
|
|
+ )
|
|
|
+ if angle >= angle_threshold_degrees and concentrated:
|
|
|
+ split_indices.append(index)
|
|
|
+ if not split_indices:
|
|
|
+ return [(xs, ys)]
|
|
|
+ chunks: list[tuple[np.ndarray, np.ndarray]] = []
|
|
|
+ start = 0
|
|
|
+ for end in split_indices + [len(path) - 1]:
|
|
|
+ chunk = path[start : end + 1]
|
|
|
+ if len(chunk) >= 2:
|
|
|
+ chunks.append((chunk[:, 0], chunk[:, 1]))
|
|
|
+ start = end + 1
|
|
|
+ return chunks or [(xs, ys)]
|
|
|
+
|
|
|
+
|
|
|
+def _extract_segments(
|
|
|
+ mask: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> tuple[
|
|
|
+ np.ndarray,
|
|
|
+ np.ndarray,
|
|
|
+ np.ndarray,
|
|
|
+ list[ShortLineSegment],
|
|
|
+ int,
|
|
|
+ _AdaptiveScale,
|
|
|
+]:
|
|
|
+ skeleton = _thin(mask)
|
|
|
+ adaptive_scale = _derive_adaptive_scale(mask, skeleton, bounds, config)
|
|
|
+ junction_seed = _junction_seeds(skeleton)
|
|
|
+ radius = adaptive_scale.junction_radius_px
|
|
|
+ kernel = cv2.getStructuringElement(
|
|
|
+ cv2.MORPH_ELLIPSE, (radius * 2 + 1, radius * 2 + 1)
|
|
|
+ )
|
|
|
+ junction_zone = cv2.dilate(junction_seed, kernel)
|
|
|
+ broken = skeleton.copy()
|
|
|
+ broken[junction_zone > 0] = 0
|
|
|
+ count, labels, stats, _ = cv2.connectedComponentsWithStats(broken, 8)
|
|
|
+ retained_labels = np.zeros_like(labels, dtype=np.int32)
|
|
|
+ segments: list[ShortLineSegment] = []
|
|
|
+ sharp_turn_epsilon = max(
|
|
|
+ min(bounds.right - bounds.left, bounds.bottom - bounds.top)
|
|
|
+ * config.sharp_turn_simplification_ratio,
|
|
|
+ adaptive_scale.line_width_px
|
|
|
+ * config.sharp_turn_simplification_line_widths,
|
|
|
+ )
|
|
|
+ for label in range(1, count):
|
|
|
+ area = int(stats[label, cv2.CC_STAT_AREA])
|
|
|
+ xs = np.flatnonzero(np.any(labels == label, axis=0))
|
|
|
+ ys = np.flatnonzero(np.any(labels == label, axis=1))
|
|
|
+ if (
|
|
|
+ area < adaptive_scale.minimum_segment_pixels
|
|
|
+ or not len(xs)
|
|
|
+ or not len(ys)
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ pixel_y, pixel_x = np.nonzero(labels == label)
|
|
|
+ local = (labels[stats[label, cv2.CC_STAT_TOP]:stats[label, cv2.CC_STAT_TOP] + stats[label, cv2.CC_STAT_HEIGHT],
|
|
|
+ stats[label, cv2.CC_STAT_LEFT]:stats[label, cv2.CC_STAT_LEFT] + stats[label, cv2.CC_STAT_WIDTH]] == label)
|
|
|
+ # A branched or closed component is not a single voting fragment.
|
|
|
+ # Never replace it with its diameter and silently discard the arms.
|
|
|
+ if np.count_nonzero(local & (_crossing_number(local) == 1)) != 2:
|
|
|
+ continue
|
|
|
+ fragments = _split_component_at_sharp_turns(
|
|
|
+ pixel_x,
|
|
|
+ pixel_y,
|
|
|
+ sharp_turn_epsilon,
|
|
|
+ config.sharp_turn_angle_degrees,
|
|
|
+ adaptive_scale.line_width_px
|
|
|
+ * config.sharp_turn_maximum_arm_line_widths,
|
|
|
+ config.sharp_turn_maximum_arm_fraction,
|
|
|
+ )
|
|
|
+ for fragment_x, fragment_y in fragments:
|
|
|
+ if len(fragment_x) < adaptive_scale.minimum_segment_pixels:
|
|
|
+ continue
|
|
|
+ span = max(float(np.ptp(fragment_x)), float(np.ptp(fragment_y)))
|
|
|
+ if span < adaptive_scale.minimum_segment_span_px:
|
|
|
+ continue
|
|
|
+ segment = _fit_segment(
|
|
|
+ f"L{len(segments) + 1:04d}",
|
|
|
+ fragment_x,
|
|
|
+ fragment_y,
|
|
|
+ bounds,
|
|
|
+ )
|
|
|
+ if segment is None:
|
|
|
+ continue
|
|
|
+ slope = segment.tangent_slope_normalized
|
|
|
+ if not (
|
|
|
+ config.minimum_tangent_slope
|
|
|
+ <= slope
|
|
|
+ <= config.maximum_tangent_slope
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ retained_labels[fragment_y, fragment_x] = len(segments) + 1
|
|
|
+ segments.append(segment)
|
|
|
+ junction_count = cv2.connectedComponents(junction_zone, 8)[0] - 1
|
|
|
+ return (
|
|
|
+ skeleton,
|
|
|
+ junction_zone,
|
|
|
+ retained_labels,
|
|
|
+ segments,
|
|
|
+ junction_count,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _curvature_samples(config: QuadraticCurveFamilyConfig) -> np.ndarray:
|
|
|
+ return np.geomspace(
|
|
|
+ config.minimum_curvature,
|
|
|
+ config.maximum_curvature,
|
|
|
+ config.curvature_sample_count,
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _range_overlap_ratio(
|
|
|
+ first: tuple[float, float], second: tuple[float, float]
|
|
|
+) -> float:
|
|
|
+ overlap = max(0.0, min(first[1], second[1]) - max(first[0], second[0]))
|
|
|
+ smaller_span = min(first[1] - first[0], second[1] - second[0])
|
|
|
+ return float(overlap / max(smaller_span, 1e-9))
|
|
|
+
|
|
|
+
|
|
|
+def _segments_can_share_curve(
|
|
|
+ first: ShortLineSegment,
|
|
|
+ second: ShortLineSegment,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> bool:
|
|
|
+ """Reject two independent observations occupying the same x/y region.
|
|
|
+
|
|
|
+ A characteristic curve is modeled as a single-valued function ``v=f(u)``.
|
|
|
+ Two disconnected skeleton fragments whose x and y ranges both overlap
|
|
|
+ substantially are competing observations in the same part of that
|
|
|
+ function, not complementary pieces that may be merged.
|
|
|
+ """
|
|
|
+
|
|
|
+ x_overlap = _range_overlap_ratio(
|
|
|
+ first.normalized_x_range, second.normalized_x_range
|
|
|
+ )
|
|
|
+ y_overlap = _range_overlap_ratio(
|
|
|
+ first.normalized_y_range, second.normalized_y_range
|
|
|
+ )
|
|
|
+ return not (
|
|
|
+ x_overlap >= config.maximum_seed_x_overlap_ratio
|
|
|
+ and y_overlap >= config.maximum_seed_y_overlap_ratio
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+def _compatible_voting_segments(
|
|
|
+ weights: dict[str, float],
|
|
|
+ segment_by_id: dict[str, ShortLineSegment],
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> tuple[tuple[str, ...], float]:
|
|
|
+ """Find a high-weight mutually compatible subset in one vote bin."""
|
|
|
+
|
|
|
+ ordered_ids = sorted(
|
|
|
+ weights, key=lambda segment_id: weights[segment_id], reverse=True
|
|
|
+ )
|
|
|
+ beam: list[tuple[tuple[str, ...], float]] = [((), 0.0)]
|
|
|
+ for segment_id in ordered_ids:
|
|
|
+ expanded = list(beam)
|
|
|
+ segment = segment_by_id[segment_id]
|
|
|
+ for selected_ids, score in beam:
|
|
|
+ if all(
|
|
|
+ _segments_can_share_curve(
|
|
|
+ segment, segment_by_id[previous_id], config
|
|
|
+ )
|
|
|
+ for previous_id in selected_ids
|
|
|
+ ):
|
|
|
+ expanded.append(
|
|
|
+ (selected_ids + (segment_id,), score + weights[segment_id])
|
|
|
+ )
|
|
|
+ expanded.sort(key=lambda item: (item[1], len(item[0])), reverse=True)
|
|
|
+ beam = expanded[: config.segment_vote_beam_width]
|
|
|
+ selected_ids, score = beam[0]
|
|
|
+ return tuple(sorted(selected_ids)), float(score)
|
|
|
+
|
|
|
+
|
|
|
+def _vote_hypotheses(
|
|
|
+ segments: list[ShortLineSegment],
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> list[_VoteHypothesis]:
|
|
|
+ # Keyed by (curvature index, quantized vertex x, quantized vertex y).
|
|
|
+ votes: dict[tuple[int, int, int], dict[str, Any]] = {}
|
|
|
+ curvature_values = _curvature_samples(config)
|
|
|
+ bin_size = config.vertex_bin_size
|
|
|
+ for segment in segments:
|
|
|
+ u = segment.normalized_midpoint.x
|
|
|
+ v = segment.normalized_midpoint.y
|
|
|
+ slope = segment.tangent_slope_normalized
|
|
|
+ segment_weight = math.sqrt(max(segment.length_px, 1.0)) / (
|
|
|
+ 1.0 + 0.20 * segment.fit_rmse_px
|
|
|
+ )
|
|
|
+ for curvature_index, curvature in enumerate(curvature_values):
|
|
|
+ vertex_x = u - slope / (2.0 * curvature)
|
|
|
+ vertex_y = v - curvature * (u - vertex_x) ** 2
|
|
|
+ if not (-0.25 <= vertex_x <= 0.88 and -0.25 <= vertex_y <= 0.82):
|
|
|
+ continue
|
|
|
+ center_x = round(vertex_x / bin_size)
|
|
|
+ center_y = round(vertex_y / bin_size)
|
|
|
+ # Soft 3x3 voting makes the accumulator tolerant of skeleton and
|
|
|
+ # local-tangent noise without blurring the curvature dimension.
|
|
|
+ for offset_x in (-1, 0, 1):
|
|
|
+ for offset_y in (-1, 0, 1):
|
|
|
+ distance = math.hypot(offset_x, offset_y)
|
|
|
+ weight = segment_weight * math.exp(-0.65 * distance * distance)
|
|
|
+ key = (
|
|
|
+ curvature_index,
|
|
|
+ center_x + offset_x,
|
|
|
+ center_y + offset_y,
|
|
|
+ )
|
|
|
+ bucket = votes.setdefault(key, {"weights": {}})
|
|
|
+ bucket["weights"][segment.segment_id] = max(
|
|
|
+ weight,
|
|
|
+ bucket["weights"].get(segment.segment_id, 0.0),
|
|
|
+ )
|
|
|
+ candidates: list[_VoteHypothesis] = []
|
|
|
+ segment_by_id = {segment.segment_id: segment for segment in segments}
|
|
|
+ for (curvature_index, vertex_x_bin, vertex_y_bin), bucket in votes.items():
|
|
|
+ segment_ids, compatible_score = _compatible_voting_segments(
|
|
|
+ bucket["weights"], segment_by_id, config
|
|
|
+ )
|
|
|
+ if len(segment_ids) < config.minimum_vote_segments:
|
|
|
+ continue
|
|
|
+ candidates.append(
|
|
|
+ _VoteHypothesis(
|
|
|
+ curvature=float(curvature_values[curvature_index]),
|
|
|
+ vertex_x=float(vertex_x_bin * bin_size),
|
|
|
+ vertex_y=float(vertex_y_bin * bin_size),
|
|
|
+ vote_score=compatible_score,
|
|
|
+ segment_ids=segment_ids,
|
|
|
+ )
|
|
|
+ )
|
|
|
+ candidates.sort(
|
|
|
+ key=lambda item: (len(item.segment_ids), item.vote_score), reverse=True
|
|
|
+ )
|
|
|
+ # Accumulator peaks occupy neighboring bins. Suppress them here so the
|
|
|
+ # expensive full-image score sees diverse parabola hypotheses.
|
|
|
+ retained: list[_VoteHypothesis] = []
|
|
|
+ for candidate in candidates:
|
|
|
+ duplicate = any(
|
|
|
+ abs(math.log(candidate.curvature / previous.curvature)) < 0.09
|
|
|
+ and abs(candidate.vertex_x - previous.vertex_x) < bin_size * 1.6
|
|
|
+ and abs(candidate.vertex_y - previous.vertex_y) < bin_size * 1.6
|
|
|
+ for previous in retained
|
|
|
+ )
|
|
|
+ if duplicate:
|
|
|
+ continue
|
|
|
+ retained.append(candidate)
|
|
|
+ if len(retained) >= config.maximum_vote_hypotheses:
|
|
|
+ break
|
|
|
+ return retained
|
|
|
+
|
|
|
+
|
|
|
+def _maximum_gap(evidence: np.ndarray, threshold: float = 0.30) -> float:
|
|
|
+ unsupported = evidence < threshold
|
|
|
+ longest = current = 0
|
|
|
+ for value in unsupported:
|
|
|
+ current = current + 1 if value else 0
|
|
|
+ longest = max(longest, current)
|
|
|
+ return float(longest / max(len(evidence), 1))
|
|
|
+
|
|
|
+
|
|
|
+def _best_supported_interval(
|
|
|
+ normalized_x: np.ndarray,
|
|
|
+ evidence: np.ndarray,
|
|
|
+ minimum_span: float,
|
|
|
+) -> tuple[int, int] | None:
|
|
|
+ reward = evidence - 0.31
|
|
|
+ prefix = np.concatenate(([0.0], np.cumsum(reward)))
|
|
|
+ best: tuple[float, int, int] | None = None
|
|
|
+ for start in range(len(normalized_x) - 1):
|
|
|
+ valid = np.flatnonzero(normalized_x[start + 1 :] - normalized_x[start] >= minimum_span)
|
|
|
+ if not len(valid):
|
|
|
+ continue
|
|
|
+ first_end = start + 1 + int(valid[0])
|
|
|
+ end = first_end + int(np.argmax(prefix[first_end + 1 :] - prefix[first_end]))
|
|
|
+ value = float(prefix[end + 1] - prefix[start])
|
|
|
+ span_bonus = 0.40 * float(normalized_x[end] - normalized_x[start])
|
|
|
+ item = (value + span_bonus, start, end)
|
|
|
+ if best is None or item[0] > best[0]:
|
|
|
+ best = item
|
|
|
+ return None if best is None else (best[1], best[2])
|
|
|
+
|
|
|
+
|
|
|
+def _distance_evidence(
|
|
|
+ distance_to_ink: np.ndarray,
|
|
|
+ points: np.ndarray,
|
|
|
+ support_scale: float,
|
|
|
+) -> tuple[np.ndarray, np.ndarray]:
|
|
|
+ xs = np.clip(np.round(points[:, 0]).astype(int), 0, distance_to_ink.shape[1] - 1)
|
|
|
+ ys = np.clip(np.round(points[:, 1]).astype(int), 0, distance_to_ink.shape[0] - 1)
|
|
|
+ distances = distance_to_ink[ys, xs]
|
|
|
+ evidence = np.exp(-np.square(distances / support_scale))
|
|
|
+ return distances, evidence
|
|
|
+
|
|
|
+
|
|
|
+def _candidate_from_equation(
|
|
|
+ coefficients: np.ndarray,
|
|
|
+ vote: _VoteHypothesis,
|
|
|
+ distance_to_ink: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+ adaptive_scale: _AdaptiveScale,
|
|
|
+) -> _ScoredCandidate | None:
|
|
|
+ width = max(1.0, bounds.right - bounds.left)
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ normalized_x = np.linspace(0.0, 1.0, config.curve_sample_count)
|
|
|
+ normalized_y = np.polyval(coefficients, normalized_x)
|
|
|
+ derivative = np.polyval(np.polyder(coefficients), normalized_x)
|
|
|
+ valid = (
|
|
|
+ (normalized_y >= -0.02)
|
|
|
+ & (normalized_y <= 1.02)
|
|
|
+ & (derivative >= config.minimum_tangent_slope * 0.45)
|
|
|
+ )
|
|
|
+ normalized_x = normalized_x[valid]
|
|
|
+ normalized_y = normalized_y[valid]
|
|
|
+ if len(normalized_x) < 20:
|
|
|
+ return None
|
|
|
+ points = np.column_stack(
|
|
|
+ (bounds.left + normalized_x * width, bounds.top + normalized_y * height)
|
|
|
+ )
|
|
|
+ support_scale = adaptive_scale.support_distance_px
|
|
|
+ distances, evidence = _distance_evidence(distance_to_ink, points, support_scale)
|
|
|
+ interval = _best_supported_interval(
|
|
|
+ normalized_x, evidence, config.minimum_x_span_ratio
|
|
|
+ )
|
|
|
+ if interval is None:
|
|
|
+ return None
|
|
|
+ start, end = interval
|
|
|
+ normalized_x = normalized_x[start : end + 1]
|
|
|
+ normalized_y = normalized_y[start : end + 1]
|
|
|
+ points = points[start : end + 1]
|
|
|
+ distances = distances[start : end + 1]
|
|
|
+ evidence = evidence[start : end + 1]
|
|
|
+ x_span = float(np.ptp(normalized_x))
|
|
|
+ y_span = float(np.ptp(normalized_y))
|
|
|
+ if x_span < config.minimum_x_span_ratio or y_span < config.minimum_y_span_ratio:
|
|
|
+ return None
|
|
|
+ support = float(np.mean(evidence))
|
|
|
+ direct_distance = adaptive_scale.direct_support_distance_px
|
|
|
+ direct_support = float(np.mean(distances <= direct_distance))
|
|
|
+ maximum_gap = _maximum_gap(evidence)
|
|
|
+ if (
|
|
|
+ support < config.minimum_pixel_support
|
|
|
+ or direct_support < config.minimum_direct_support
|
|
|
+ or maximum_gap > config.maximum_gap_ratio
|
|
|
+ ):
|
|
|
+ return None
|
|
|
+ derivative = np.polyval(np.polyder(coefficients), normalized_x)
|
|
|
+ tangent_agreement = float(np.mean(derivative > config.minimum_tangent_slope))
|
|
|
+ vote_strength = math.log1p(vote.vote_score)
|
|
|
+ score = (
|
|
|
+ 3.8 * support
|
|
|
+ + 2.2 * direct_support
|
|
|
+ + 2.0 * min(y_span / 0.55, 1.3)
|
|
|
+ + 1.0 * min(x_span / 0.35, 1.3)
|
|
|
+ + 0.55 * tangent_agreement
|
|
|
+ + 0.22 * vote_strength
|
|
|
+ - 2.4 * maximum_gap
|
|
|
+ )
|
|
|
+ return _ScoredCandidate(
|
|
|
+ coefficients=coefficients,
|
|
|
+ vertex=(vote.vertex_x, vote.vertex_y),
|
|
|
+ domain=(float(normalized_x[0]), float(normalized_x[-1])),
|
|
|
+ points=points,
|
|
|
+ segment_ids=vote.segment_ids,
|
|
|
+ x_span=x_span,
|
|
|
+ y_span=y_span,
|
|
|
+ support=support,
|
|
|
+ direct_support=direct_support,
|
|
|
+ maximum_gap=maximum_gap,
|
|
|
+ tangent_agreement=tangent_agreement,
|
|
|
+ vote_score=vote.vote_score,
|
|
|
+ score=float(score),
|
|
|
+ family_position=float(np.median(normalized_x)),
|
|
|
+ proposed_segment_ids=vote.proposed_segment_ids or vote.segment_ids,
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+@_image_cached
|
|
|
+def _column_observations(
|
|
|
+ mask: np.ndarray,
|
|
|
+ x: int,
|
|
|
+ top: int,
|
|
|
+ bottom: int,
|
|
|
+ half_window: int,
|
|
|
+ maximum_band_height: int,
|
|
|
+) -> list[float]:
|
|
|
+ left = max(0, x - half_window)
|
|
|
+ right = min(mask.shape[1], x + half_window + 1)
|
|
|
+ active = np.any(mask[top : bottom + 1, left:right] > 0, axis=1)
|
|
|
+ indices = np.flatnonzero(active)
|
|
|
+ if not len(indices):
|
|
|
+ return []
|
|
|
+ breaks = np.flatnonzero(np.diff(indices) > 1)
|
|
|
+ return [
|
|
|
+ float(top + np.median(group))
|
|
|
+ for group in np.split(indices, breaks + 1)
|
|
|
+ if 0 < len(group) <= maximum_band_height
|
|
|
+ ]
|
|
|
+
|
|
|
+
|
|
|
+def _refine_candidate(
|
|
|
+ candidate: _ScoredCandidate,
|
|
|
+ mask: np.ndarray,
|
|
|
+ distance_to_ink: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+ adaptive_scale: _AdaptiveScale,
|
|
|
+) -> _ScoredCandidate | None:
|
|
|
+ width = max(1.0, bounds.right - bounds.left)
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ coefficients = candidate.coefficients.copy()
|
|
|
+ corridor = adaptive_scale.refinement_corridor_px
|
|
|
+ top = max(0, math.floor(bounds.top))
|
|
|
+ bottom = min(mask.shape[0] - 1, math.ceil(bounds.bottom))
|
|
|
+ half_window = adaptive_scale.column_half_window_px
|
|
|
+ for _ in range(4):
|
|
|
+ xs = np.arange(
|
|
|
+ max(0, math.ceil(bounds.left + candidate.domain[0] * width)),
|
|
|
+ min(mask.shape[1] - 1, math.floor(bounds.left + candidate.domain[1] * width)) + 1,
|
|
|
+ max(1, round(width * 0.002)),
|
|
|
+ )
|
|
|
+ u = (xs - bounds.left) / width
|
|
|
+ predicted = bounds.top + np.polyval(coefficients, u) * height
|
|
|
+ observations: list[tuple[float, float]] = []
|
|
|
+ for x, expected in zip(xs, predicted):
|
|
|
+ values = _column_observations(
|
|
|
+ mask,
|
|
|
+ int(x),
|
|
|
+ top,
|
|
|
+ bottom,
|
|
|
+ half_window,
|
|
|
+ adaptive_scale.maximum_column_band_px,
|
|
|
+ )
|
|
|
+ if not values:
|
|
|
+ continue
|
|
|
+ nearest = min(values, key=lambda value: abs(value - expected))
|
|
|
+ if abs(nearest - expected) <= corridor:
|
|
|
+ observations.append((float(x), nearest))
|
|
|
+ if len(observations) < 12:
|
|
|
+ break
|
|
|
+ array = np.asarray(observations)
|
|
|
+ fit_x = (array[:, 0] - bounds.left) / width
|
|
|
+ fit_y = (array[:, 1] - bounds.top) / height
|
|
|
+ selected = np.ones(len(array), dtype=bool)
|
|
|
+ for _ in range(4):
|
|
|
+ fitted = np.polyfit(fit_x[selected], fit_y[selected], 2)
|
|
|
+ residuals = fit_y - np.polyval(fitted, fit_x)
|
|
|
+ median = float(np.median(residuals[selected]))
|
|
|
+ mad = float(np.median(np.abs(residuals[selected] - median)))
|
|
|
+ tolerance = max(0.003, 3.0 * 1.4826 * mad)
|
|
|
+ updated = np.abs(residuals - median) <= tolerance
|
|
|
+ if np.count_nonzero(updated) < 8 or np.array_equal(updated, selected):
|
|
|
+ break
|
|
|
+ selected = updated
|
|
|
+ if fitted[0] <= config.minimum_curvature * 0.35:
|
|
|
+ break
|
|
|
+ coefficients = fitted
|
|
|
+ curvature = float(coefficients[0])
|
|
|
+ vertex_x = float(-coefficients[1] / (2.0 * curvature))
|
|
|
+ vertex_y = float(np.polyval(coefficients, vertex_x))
|
|
|
+ refined_vote = _VoteHypothesis(
|
|
|
+ curvature=curvature,
|
|
|
+ vertex_x=vertex_x,
|
|
|
+ vertex_y=vertex_y,
|
|
|
+ vote_score=candidate.vote_score,
|
|
|
+ segment_ids=candidate.segment_ids,
|
|
|
+ proposed_segment_ids=candidate.proposed_segment_ids,
|
|
|
+ )
|
|
|
+ return _candidate_from_equation(
|
|
|
+ coefficients,
|
|
|
+ refined_vote,
|
|
|
+ distance_to_ink,
|
|
|
+ bounds,
|
|
|
+ config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+@_image_cached
|
|
|
+def _candidate_distance(
|
|
|
+ first: _ScoredCandidate,
|
|
|
+ second: _ScoredCandidate,
|
|
|
+ bounds: PlotBounds,
|
|
|
+) -> tuple[float, float, float]:
|
|
|
+ left = max(first.domain[0], second.domain[0])
|
|
|
+ right = min(first.domain[1], second.domain[1])
|
|
|
+ if right <= left:
|
|
|
+ return math.inf, 0.0, 0.0
|
|
|
+ xs = np.linspace(left, right, 90)
|
|
|
+ first_y = np.polyval(first.coefficients, xs)
|
|
|
+ second_y = np.polyval(second.coefficients, xs)
|
|
|
+ difference = first_y - second_y
|
|
|
+ stable = difference[np.abs(difference) > 0.003]
|
|
|
+ crossings = float(np.count_nonzero(np.sign(stable[1:]) != np.sign(stable[:-1]))) if len(stable) > 1 else 0.0
|
|
|
+ overlap = (right - left) / max(
|
|
|
+ min(first.x_span, second.x_span), 1e-6
|
|
|
+ )
|
|
|
+ distance_px = float(np.median(np.abs(difference))) * max(
|
|
|
+ 1.0, bounds.bottom - bounds.top
|
|
|
+ )
|
|
|
+ return distance_px, float(overlap), crossings
|
|
|
+
|
|
|
+
|
|
|
+def _deduplicate_candidates(
|
|
|
+ candidates: list[_ScoredCandidate],
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> list[_ScoredCandidate]:
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ retained: list[_ScoredCandidate] = []
|
|
|
+ for candidate in sorted(candidates, key=lambda item: item.score, reverse=True):
|
|
|
+ if any(
|
|
|
+ (lambda values: values[1] > 0.48 and values[0] < height * config.duplicate_distance_ratio)(
|
|
|
+ _candidate_distance(candidate, previous, bounds)
|
|
|
+ )
|
|
|
+ for previous in retained
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ retained.append(candidate)
|
|
|
+ if len(retained) >= config.maximum_candidate_count:
|
|
|
+ break
|
|
|
+ return retained
|
|
|
+
|
|
|
+
|
|
|
+def _compatible_with_family(
|
|
|
+ candidate: _ScoredCandidate,
|
|
|
+ selected: list[_ScoredCandidate],
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+ distinct_endpoints: bool = False,
|
|
|
+) -> bool:
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ for previous in selected:
|
|
|
+ if distinct_endpoints and (
|
|
|
+ abs(candidate.domain[1] - previous.domain[1])
|
|
|
+ < config.endpoint_suppression_radius
|
|
|
+ ):
|
|
|
+ return False
|
|
|
+ distance, overlap, crossings = _candidate_distance(candidate, previous, bounds)
|
|
|
+ if crossings > 0.0 and overlap > 0.18:
|
|
|
+ return False
|
|
|
+ if overlap > 0.50 and distance < height * config.minimum_family_separation_ratio:
|
|
|
+ return False
|
|
|
+ return True
|
|
|
+
|
|
|
+
|
|
|
+def _lower_endpoint(candidate: _ScoredCandidate) -> tuple[float, float]:
|
|
|
+ """Return the normalized right/lower endpoint of a candidate."""
|
|
|
+
|
|
|
+ endpoint_x = candidate.domain[1]
|
|
|
+ endpoint_y = float(np.polyval(candidate.coefficients, endpoint_x))
|
|
|
+ return endpoint_x, endpoint_y
|
|
|
+
|
|
|
+
|
|
|
+def _endpoint_frontier_candidates(
|
|
|
+ candidates: list[_ScoredCandidate],
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> list[_ScoredCandidate]:
|
|
|
+ """Keep one smooth outer-envelope candidate per endpoint position.
|
|
|
+
|
|
|
+ Characteristic curves in this chart family terminate along a smooth lower
|
|
|
+ envelope. Efficiency contours may individually fit a parabola, but their
|
|
|
+ endpoints usually lie inside that envelope or duplicate the horizontal
|
|
|
+ position of a stronger characteristic curve.
|
|
|
+ """
|
|
|
+
|
|
|
+ if len(candidates) <= 3:
|
|
|
+ return candidates
|
|
|
+ endpoints = [_lower_endpoint(candidate) for candidate in candidates]
|
|
|
+ frontier_indices: list[int] = []
|
|
|
+ for index, (endpoint_x, endpoint_y) in enumerate(endpoints):
|
|
|
+ local_bottom = max(
|
|
|
+ other_y
|
|
|
+ for other_x, other_y in endpoints
|
|
|
+ if abs(other_x - endpoint_x) <= config.endpoint_frontier_radius
|
|
|
+ )
|
|
|
+ if local_bottom - endpoint_y <= config.endpoint_frontier_tolerance:
|
|
|
+ frontier_indices.append(index)
|
|
|
+
|
|
|
+ ranked = sorted(
|
|
|
+ frontier_indices,
|
|
|
+ key=lambda index: (
|
|
|
+ candidates[index].score + 0.8 * candidates[index].x_span,
|
|
|
+ candidates[index].direct_support,
|
|
|
+ ),
|
|
|
+ reverse=True,
|
|
|
+ )
|
|
|
+ suppressed: list[int] = []
|
|
|
+ for index in ranked:
|
|
|
+ endpoint_x, _ = endpoints[index]
|
|
|
+ if any(
|
|
|
+ abs(endpoint_x - endpoints[previous][0])
|
|
|
+ < config.endpoint_suppression_radius
|
|
|
+ for previous in suppressed
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ suppressed.append(index)
|
|
|
+ suppressed.sort(key=lambda index: endpoints[index][0])
|
|
|
+ if len(suppressed) < 5:
|
|
|
+ return [candidates[index] for index in suppressed]
|
|
|
+
|
|
|
+ endpoint_array = np.asarray([endpoints[index] for index in suppressed])
|
|
|
+ best: tuple[tuple[int, float, float], np.ndarray] | None = None
|
|
|
+ for sample_indices in combinations(range(len(suppressed)), 3):
|
|
|
+ sample = endpoint_array[list(sample_indices)]
|
|
|
+ coefficients = np.polyfit(sample[:, 0], sample[:, 1], 2)
|
|
|
+ if not (
|
|
|
+ config.endpoint_locus_minimum_curvature
|
|
|
+ <= coefficients[0]
|
|
|
+ <= config.endpoint_locus_maximum_curvature
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ residuals = np.abs(
|
|
|
+ endpoint_array[:, 1]
|
|
|
+ - np.polyval(coefficients, endpoint_array[:, 0])
|
|
|
+ )
|
|
|
+ inliers = residuals <= config.endpoint_locus_initial_tolerance
|
|
|
+ if np.count_nonzero(inliers) < 3:
|
|
|
+ continue
|
|
|
+ refitted = np.polyfit(
|
|
|
+ endpoint_array[inliers, 0], endpoint_array[inliers, 1], 2
|
|
|
+ )
|
|
|
+ if not (
|
|
|
+ config.endpoint_locus_minimum_curvature
|
|
|
+ <= refitted[0]
|
|
|
+ <= config.endpoint_locus_maximum_curvature
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ residuals = np.abs(
|
|
|
+ endpoint_array[:, 1]
|
|
|
+ - np.polyval(refitted, endpoint_array[:, 0])
|
|
|
+ )
|
|
|
+ inliers = residuals <= config.endpoint_locus_final_tolerance
|
|
|
+ inlier_count = int(np.count_nonzero(inliers))
|
|
|
+ if inlier_count < 3:
|
|
|
+ continue
|
|
|
+ median_residual = float(np.median(residuals[inliers]))
|
|
|
+ covered_span = float(np.ptp(endpoint_array[inliers, 0]))
|
|
|
+ key = (inlier_count, covered_span, -median_residual)
|
|
|
+ if best is None or key > best[0]:
|
|
|
+ best = (key, inliers)
|
|
|
+ if best is None:
|
|
|
+ accepted_endpoint_x = [endpoints[index][0] for index in suppressed]
|
|
|
+ else:
|
|
|
+ accepted_endpoint_x = [
|
|
|
+ endpoints[index][0]
|
|
|
+ for index, is_inlier in zip(suppressed, best[1])
|
|
|
+ if is_inlier
|
|
|
+ ]
|
|
|
+
|
|
|
+ # Restore close alternatives inside every accepted endpoint cluster. The
|
|
|
+ # highest-scoring local fit can cross a neighbor; the global beam must be
|
|
|
+ # allowed to choose a compatible runner-up from the same cluster.
|
|
|
+ alternative_radius = config.endpoint_suppression_radius * 0.52
|
|
|
+ return [
|
|
|
+ candidates[index]
|
|
|
+ for index in frontier_indices
|
|
|
+ if any(
|
|
|
+ abs(endpoints[index][0] - accepted_x) <= alternative_radius
|
|
|
+ for accepted_x in accepted_endpoint_x
|
|
|
+ )
|
|
|
+ ]
|
|
|
+
|
|
|
+
|
|
|
+@_image_cached
|
|
|
+def _family_objective(
|
|
|
+ selected: list[_ScoredCandidate],
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+) -> float:
|
|
|
+ if not selected:
|
|
|
+ return 0.0
|
|
|
+ ordered = sorted(selected, key=lambda item: item.family_position)
|
|
|
+ quality = sum(max(item.score - 5.5, 0.0) for item in ordered)
|
|
|
+ positions = np.asarray([item.family_position for item in ordered])
|
|
|
+ curvatures = np.asarray([item.coefficients[0] for item in ordered])
|
|
|
+ objective = (
|
|
|
+ quality
|
|
|
+ + config.family_size_reward * len(ordered)
|
|
|
+ + 1.5 * float(np.ptp(positions))
|
|
|
+ )
|
|
|
+ if len(ordered) >= 3:
|
|
|
+ gaps = np.diff(positions)
|
|
|
+ median_gap = max(float(np.median(gaps)), 1e-4)
|
|
|
+ objective -= config.family_gap_irregularity_penalty * float(
|
|
|
+ np.median(np.abs(gaps - median_gap)) / median_gap
|
|
|
+ )
|
|
|
+ log_curvature = np.log(np.maximum(curvatures, 1e-6))
|
|
|
+ objective -= 0.65 * float(np.median(np.abs(np.diff(log_curvature))))
|
|
|
+ return float(objective)
|
|
|
+
|
|
|
+
|
|
|
+def _select_family(
|
|
|
+ candidates: list[_ScoredCandidate],
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+ distinct_endpoints: bool = False,
|
|
|
+) -> tuple[list[_ScoredCandidate], float]:
|
|
|
+ ordered = sorted(candidates, key=lambda item: item.family_position)
|
|
|
+ beam: list[list[_ScoredCandidate]] = [[]]
|
|
|
+ for candidate in ordered:
|
|
|
+ expanded = list(beam)
|
|
|
+ for selected in beam:
|
|
|
+ if len(selected) >= config.maximum_family_size:
|
|
|
+ continue
|
|
|
+ if _compatible_with_family(
|
|
|
+ candidate,
|
|
|
+ selected,
|
|
|
+ bounds,
|
|
|
+ config,
|
|
|
+ distinct_endpoints=distinct_endpoints,
|
|
|
+ ):
|
|
|
+ expanded.append(selected + [candidate])
|
|
|
+ expanded.sort(
|
|
|
+ key=lambda items: _family_objective(items, bounds, config), reverse=True
|
|
|
+ )
|
|
|
+ if config.target_family_size is None:
|
|
|
+ beam = expanded[: config.family_beam_width]
|
|
|
+ else:
|
|
|
+ # A high-scoring larger partial family must not evict every path to
|
|
|
+ # the requested cardinality before later compatible candidates are
|
|
|
+ # visited. Preserve a bounded beam independently for each size.
|
|
|
+ quota = max(
|
|
|
+ 64,
|
|
|
+ config.family_beam_width // (config.target_family_size + 1),
|
|
|
+ )
|
|
|
+ counts: dict[int, int] = {}
|
|
|
+ beam = []
|
|
|
+ for item in expanded:
|
|
|
+ size = len(item)
|
|
|
+ if size > config.target_family_size:
|
|
|
+ continue
|
|
|
+ if counts.get(size, 0) >= quota:
|
|
|
+ continue
|
|
|
+ beam.append(item)
|
|
|
+ counts[size] = counts.get(size, 0) + 1
|
|
|
+ if config.target_family_size is not None:
|
|
|
+ eligible = [
|
|
|
+ items for items in beam if len(items) == config.target_family_size
|
|
|
+ ]
|
|
|
+ if not eligible:
|
|
|
+ raise ValueError(
|
|
|
+ "像素证据无法组成语义要求的"
|
|
|
+ f"{config.target_family_size}条互相兼容曲线"
|
|
|
+ )
|
|
|
+ else:
|
|
|
+ eligible = [items for items in beam if len(items) >= 2]
|
|
|
+ selected = max(
|
|
|
+ eligible or beam,
|
|
|
+ key=lambda items: _family_objective(items, bounds, config),
|
|
|
+ )
|
|
|
+ selected.sort(key=lambda item: item.family_position)
|
|
|
+ return selected, _family_objective(selected, bounds, config)
|
|
|
+
|
|
|
+
|
|
|
+def _revalidate_vote_seeds(
|
|
|
+ candidates: list[_ScoredCandidate], segments: list[ShortLineSegment],
|
|
|
+ labels: np.ndarray, bounds: PlotBounds, config: QuadraticCurveFamilyConfig,
|
|
|
+ adaptive: _AdaptiveScale,
|
|
|
+) -> list[_ScoredCandidate]:
|
|
|
+ """A vote proposes membership; absolute pixel evidence can revoke it."""
|
|
|
+ pixels = {}
|
|
|
+ weights = {}
|
|
|
+ for label, segment in enumerate(segments, 1):
|
|
|
+ y, x = np.nonzero(labels == label)
|
|
|
+ pixels[segment.segment_id] = np.column_stack((
|
|
|
+ (x - bounds.left) / max(bounds.right - bounds.left, 1),
|
|
|
+ (y - bounds.top) / max(bounds.bottom - bounds.top, 1),
|
|
|
+ ))
|
|
|
+ weights[segment.segment_id] = math.sqrt(max(segment.length_px, 1)) / (1 + .2 * segment.fit_rmse_px)
|
|
|
+ accepted = []
|
|
|
+ for candidate in candidates:
|
|
|
+ kept = []
|
|
|
+ for sid in candidate.segment_ids:
|
|
|
+ if sid not in pixels or not len(pixels[sid]):
|
|
|
+ continue
|
|
|
+ ev = _segment_fit_evidence(candidate.coefficients, pixels[sid], bounds,
|
|
|
+ adaptive.line_width_px * config.association_joint_overlap_line_widths)
|
|
|
+ if (ev.overlap_fraction >= config.association_joint_minimum_overlap
|
|
|
+ and ev.median_distance_px <= adaptive.line_width_px * config.association_joint_median_line_widths
|
|
|
+ and ev.p80_distance_px <= adaptive.line_width_px * config.association_joint_p80_line_widths):
|
|
|
+ kept.append(sid)
|
|
|
+ while kept:
|
|
|
+ joint = _fit_segment_group(tuple(kept), pixels, frozenset(kept),
|
|
|
+ candidate.coefficients, bounds, adaptive.line_width_px, config)
|
|
|
+ if _segment_group_is_consistent(joint, joint, adaptive.line_width_px, config):
|
|
|
+ break
|
|
|
+ worst = max(kept, key=lambda sid: (
|
|
|
+ joint.evidence[sid].median_distance_px + joint.evidence[sid].p80_distance_px,
|
|
|
+ sid,
|
|
|
+ ))
|
|
|
+ kept.remove(worst)
|
|
|
+ if not kept:
|
|
|
+ continue
|
|
|
+ ratio = sum(weights[sid] for sid in kept) / max(sum(weights.get(sid, 0) for sid in candidate.segment_ids), 1e-9)
|
|
|
+ vote_score = candidate.vote_score * ratio
|
|
|
+ accepted.append(replace(candidate, segment_ids=tuple(kept), vote_score=vote_score,
|
|
|
+ score=candidate.score + .22 * (math.log1p(vote_score) - math.log1p(candidate.vote_score))))
|
|
|
+ return accepted
|
|
|
+
|
|
|
+
|
|
|
+def _associate_segments_to_final_curves(
|
|
|
+ candidates: list[_ScoredCandidate],
|
|
|
+ segments: list[ShortLineSegment],
|
|
|
+ segment_labels: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+ adaptive_scale: _AdaptiveScale,
|
|
|
+) -> list[tuple[str, ...]]:
|
|
|
+ """Fill internal holes and extend a non-branching, jointly fitted path.
|
|
|
+
|
|
|
+ Candidate discovery covers the whole plot, not just the outermost seeds.
|
|
|
+ A substantial part of a fragment must lie along the extended equation
|
|
|
+ before it can be proposed. Every addition must then preserve all red and
|
|
|
+ orange fragments in a joint fit. Gap and local curvature are soft ranking
|
|
|
+ terms; overlapping alternatives compete in a canonical ordered beam.
|
|
|
+ """
|
|
|
+
|
|
|
+ width = max(1.0, bounds.right - bounds.left)
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ gap_scale = max(min(width, height) * config.association_gap_scale_ratio, 1e-9)
|
|
|
+ line_width = max(adaptive_scale.line_width_px, 1e-9)
|
|
|
+ segment_by_id = {segment.segment_id: segment for segment in segments}
|
|
|
+ endpoints: dict[str, tuple[np.ndarray, np.ndarray]] = {}
|
|
|
+ normalized_pixels: dict[str, np.ndarray] = {}
|
|
|
+ for label, segment in enumerate(segments, start=1):
|
|
|
+ pixel_y, pixel_x = np.nonzero(segment_labels == label)
|
|
|
+ points = np.column_stack((pixel_x, pixel_y)).astype(float)
|
|
|
+ if not len(points):
|
|
|
+ continue
|
|
|
+ center = points.mean(axis=0)
|
|
|
+ centered = points - center
|
|
|
+ _, eigenvectors = np.linalg.eigh(centered.T @ centered)
|
|
|
+ direction = eigenvectors[:, -1]
|
|
|
+ if direction[0] < 0.0:
|
|
|
+ direction = -direction
|
|
|
+ projections = centered @ direction
|
|
|
+ endpoints[segment.segment_id] = (
|
|
|
+ center + float(np.min(projections)) * direction,
|
|
|
+ center + float(np.max(projections)) * direction,
|
|
|
+ )
|
|
|
+ normalized_pixels[segment.segment_id] = np.column_stack(
|
|
|
+ (
|
|
|
+ (pixel_x.astype(float) - bounds.left) / width,
|
|
|
+ (pixel_y.astype(float) - bounds.top) / height,
|
|
|
+ )
|
|
|
+ )
|
|
|
+
|
|
|
+ def ordered(items: tuple[str, ...] | frozenset[str]) -> tuple[str, ...]:
|
|
|
+ return tuple(sorted(
|
|
|
+ items, key=lambda item: (segment_by_id[item].normalized_midpoint.x, item)
|
|
|
+ ))
|
|
|
+
|
|
|
+ @lru_cache(maxsize=32768)
|
|
|
+ def can_insert(segment_id: str, occupied: tuple[str, ...]) -> bool:
|
|
|
+ segment = segment_by_id[segment_id]
|
|
|
+ left, right = segment.normalized_x_range
|
|
|
+ for item in occupied:
|
|
|
+ previous = segment_by_id[item]
|
|
|
+ if not _segments_can_share_curve(segment, previous, config):
|
|
|
+ return False
|
|
|
+ other_left, other_right = previous.normalized_x_range
|
|
|
+ overlap = max(0.0, min(right, other_right) - max(left, other_left))
|
|
|
+ shorter_span = max(min(right - left, other_right - other_left), 1e-9)
|
|
|
+ # A single-valued path cannot take two parallel branches occupying
|
|
|
+ # the same x slot, even when their y ranges do not overlap.
|
|
|
+ if overlap / shorter_span >= config.maximum_seed_x_overlap_ratio:
|
|
|
+ return False
|
|
|
+ return True
|
|
|
+
|
|
|
+ def fragment_score(
|
|
|
+ segment_id: str, path_ids: tuple[str, ...], fit: _SegmentGroupFit,
|
|
|
+ ) -> float:
|
|
|
+ segment = segment_by_id[segment_id]
|
|
|
+ index = path_ids.index(segment_id)
|
|
|
+ neighbors = path_ids[max(0, index - 1):index] + path_ids[index + 1:index + 2]
|
|
|
+ slope = segment.tangent_slope_normalized
|
|
|
+ predicted_slope = float(np.polyval(
|
|
|
+ np.polyder(fit.coefficients), segment.normalized_midpoint.x
|
|
|
+ ))
|
|
|
+ tangent_limit = max(
|
|
|
+ config.association_tangent_absolute_tolerance,
|
|
|
+ config.association_tangent_relative_tolerance
|
|
|
+ * max(abs(slope), abs(predicted_slope)),
|
|
|
+ )
|
|
|
+ gaps: list[float] = []
|
|
|
+ curvature_errors: list[float] = []
|
|
|
+ for item in neighbors:
|
|
|
+ neighbor = segment_by_id[item]
|
|
|
+ delta_x = segment.normalized_midpoint.x - neighbor.normalized_midpoint.x
|
|
|
+ if delta_x >= 0:
|
|
|
+ left_end, right_start = endpoints[item][1], endpoints[segment_id][0]
|
|
|
+ else:
|
|
|
+ left_end, right_start = endpoints[segment_id][1], endpoints[item][0]
|
|
|
+ gaps.append(float(np.linalg.norm(right_start - left_end)))
|
|
|
+ expected_change = 2.0 * float(fit.coefficients[0]) * delta_x
|
|
|
+ observed_change = slope - neighbor.tangent_slope_normalized
|
|
|
+ tolerance = max(
|
|
|
+ tangent_limit,
|
|
|
+ config.association_tangent_relative_tolerance
|
|
|
+ * abs(neighbor.tangent_slope_normalized),
|
|
|
+ )
|
|
|
+ # Do not divide by delta_x: noisy slopes on adjacent short pieces
|
|
|
+ # otherwise imply arbitrarily large or nearly zero curvature.
|
|
|
+ curvature_errors.append(abs(observed_change - expected_change) / tolerance)
|
|
|
+ evidence = fit.evidence[segment_id]
|
|
|
+ fit_error = (
|
|
|
+ evidence.median_distance_px / line_width
|
|
|
+ + 0.5 * evidence.p80_distance_px / line_width
|
|
|
+ + 1.0 - evidence.overlap_fraction
|
|
|
+ )
|
|
|
+ return float(
|
|
|
+ config.association_distance_weight * math.log1p(min(gaps, default=0) / gap_scale)
|
|
|
+ + 0.55 * abs(slope - predicted_slope) / tangent_limit
|
|
|
+ + config.association_curvature_weight * np.mean(curvature_errors or [0.0])
|
|
|
+ + 0.25 * fit_error
|
|
|
+ )
|
|
|
+
|
|
|
+ all_seed_ids = {
|
|
|
+ segment_id for candidate in candidates for segment_id in candidate.segment_ids
|
|
|
+ }
|
|
|
+ proposal_scores: list[dict[str, float]] = []
|
|
|
+ best_paths: list[_AssociationPath | None] = []
|
|
|
+ baselines: list[_SegmentGroupFit | None] = []
|
|
|
+ for candidate in candidates:
|
|
|
+ seed_ids = ordered(tuple(item for item in candidate.segment_ids if item in endpoints))
|
|
|
+ if not seed_ids:
|
|
|
+ proposal_scores.append({})
|
|
|
+ best_paths.append(None)
|
|
|
+ baselines.append(None)
|
|
|
+ continue
|
|
|
+ seeds = frozenset(seed_ids)
|
|
|
+ fit_cache: dict[tuple[str, ...], _SegmentGroupFit] = {}
|
|
|
+
|
|
|
+ def fit_group(items: tuple[str, ...]) -> _SegmentGroupFit:
|
|
|
+ key = tuple(sorted(items))
|
|
|
+ if key not in fit_cache:
|
|
|
+ fit_cache[key] = _fit_segment_group(
|
|
|
+ key, normalized_pixels, seeds, candidate.coefficients,
|
|
|
+ bounds, adaptive_scale.line_width_px, config,
|
|
|
+ )
|
|
|
+ return fit_cache[key]
|
|
|
+
|
|
|
+ seed_baseline = fit_group(seed_ids)
|
|
|
+ baselines.append(seed_baseline)
|
|
|
+
|
|
|
+ objective_cache: dict[tuple[str, ...], float] = {}
|
|
|
+
|
|
|
+ def objective(state: _AssociationPath) -> float:
|
|
|
+ if state.added in objective_cache:
|
|
|
+ return objective_cache[state.added]
|
|
|
+ # Reward supported ink, not the empty distance jumped over.
|
|
|
+ support = sum(segment_by_id[item].length_px for item in state.added)
|
|
|
+ path_ids = ordered(seed_ids + state.added)
|
|
|
+ penalty = sum(fragment_score(item, path_ids, state.fit) for item in state.added)
|
|
|
+ value = float(
|
|
|
+ support / gap_scale
|
|
|
+ + 0.15 * len(state.added)
|
|
|
+ - 0.35 * penalty
|
|
|
+ - 0.35 * max(0.0, state.fit.loss - seed_baseline.loss)
|
|
|
+ )
|
|
|
+ objective_cache[state.added] = value
|
|
|
+ return value
|
|
|
+
|
|
|
+ initial = _AssociationPath((), seed_baseline)
|
|
|
+ beam = [initial]
|
|
|
+ best = initial
|
|
|
+ for _ in range(len(segments)):
|
|
|
+ expanded: dict[tuple[str, ...], _AssociationPath] = {}
|
|
|
+ for state in beam:
|
|
|
+ occupied = ordered(seed_ids + state.added)
|
|
|
+ options: list[tuple[float, _AssociationPath]] = []
|
|
|
+ for segment in segments:
|
|
|
+ segment_id = segment.segment_id
|
|
|
+ if segment_id in occupied or (
|
|
|
+ segment_id in all_seed_ids
|
|
|
+ and segment_id not in candidate.segment_ids
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ if segment_id not in endpoints or not can_insert(segment_id, occupied):
|
|
|
+ continue
|
|
|
+ proposal = _segment_fit_evidence(
|
|
|
+ state.fit.coefficients, normalized_pixels[segment_id], bounds,
|
|
|
+ line_width * config.association_proposal_line_widths,
|
|
|
+ )
|
|
|
+ if proposal.overlap_fraction < config.association_joint_minimum_overlap:
|
|
|
+ continue
|
|
|
+ trial_fit = fit_group(occupied + (segment_id,))
|
|
|
+ if not _segment_group_is_consistent(
|
|
|
+ trial_fit, seed_baseline, adaptive_scale.line_width_px, config
|
|
|
+ ):
|
|
|
+ continue
|
|
|
+ new_state = _AssociationPath(
|
|
|
+ ordered(state.added + (segment_id,)), trial_fit,
|
|
|
+ )
|
|
|
+ options.append((objective(new_state), new_state))
|
|
|
+ branches = sorted(
|
|
|
+ options, key=lambda item: (-item[0], item[1].added),
|
|
|
+ )[: config.association_path_branch_limit]
|
|
|
+ for _, new_state in branches:
|
|
|
+ # The same subset reached in a different insertion order
|
|
|
+ # has exactly the same fit and score; keep it only once.
|
|
|
+ expanded[new_state.added] = new_state
|
|
|
+ if not expanded:
|
|
|
+ break
|
|
|
+ beam = sorted(expanded.values(), key=objective, reverse=True)[
|
|
|
+ : config.association_path_beam_width
|
|
|
+ ]
|
|
|
+ if objective(beam[0]) > objective(best):
|
|
|
+ best = beam[0]
|
|
|
+ best_paths.append(best)
|
|
|
+ accepted = {segment_id: fragment_score(segment_id, ordered(seed_ids + best.added), best.fit)
|
|
|
+ for segment_id in seed_ids}
|
|
|
+ path_ids = ordered(seed_ids + best.added)
|
|
|
+ accepted.update({
|
|
|
+ item: fragment_score(item, path_ids, best.fit) for item in best.added
|
|
|
+ })
|
|
|
+ proposal_scores.append(accepted)
|
|
|
+
|
|
|
+ owners: dict[str, tuple[float, int]] = {}
|
|
|
+ for curve_index, scores in enumerate(proposal_scores):
|
|
|
+ for segment_id, score in scores.items():
|
|
|
+ ownership = (score, curve_index)
|
|
|
+ if segment_id not in owners or ownership < owners[segment_id]:
|
|
|
+ owners[segment_id] = ownership
|
|
|
+ associations: list[tuple[str, ...]] = []
|
|
|
+ for curve_index, (candidate, state, baseline) in enumerate(zip(
|
|
|
+ candidates, best_paths, baselines, strict=True
|
|
|
+ )):
|
|
|
+ owned_seeds = tuple(
|
|
|
+ item for item in candidate.segment_ids
|
|
|
+ if owners.get(item, (0, -1))[1] == curve_index
|
|
|
+ )
|
|
|
+ if not owned_seeds or state is None or baseline is None:
|
|
|
+ associations.append(owned_seeds)
|
|
|
+ continue
|
|
|
+ added = [item for item in state.added if owners[item][1] == curve_index]
|
|
|
+ while added:
|
|
|
+ ids = owned_seeds + tuple(added)
|
|
|
+ fit = _fit_segment_group(
|
|
|
+ ids, normalized_pixels, frozenset(owned_seeds), candidate.coefficients,
|
|
|
+ bounds, adaptive_scale.line_width_px, config,
|
|
|
+ )
|
|
|
+ if _segment_group_is_consistent(
|
|
|
+ fit, baseline, adaptive_scale.line_width_px, config
|
|
|
+ ):
|
|
|
+ break
|
|
|
+ # Ownership can change the fit. Revalidate the whole remaining
|
|
|
+ # path, including its interior, without deleting unrelated tails.
|
|
|
+ worst = max(
|
|
|
+ added, key=lambda item: (
|
|
|
+ fit.evidence[item].median_distance_px / line_width
|
|
|
+ + fit.evidence[item].p80_distance_px / line_width
|
|
|
+ + 1.0 - fit.evidence[item].overlap_fraction,
|
|
|
+ item,
|
|
|
+ ),
|
|
|
+ )
|
|
|
+ added.remove(worst)
|
|
|
+ associations.append(tuple(sorted(owned_seeds + tuple(added))))
|
|
|
+ return associations
|
|
|
+
|
|
|
+
|
|
|
+def _refit_candidate_from_associated_segments(
|
|
|
+ candidate: _ScoredCandidate,
|
|
|
+ associated_segment_ids: tuple[str, ...],
|
|
|
+ segments: list[ShortLineSegment],
|
|
|
+ segment_labels: np.ndarray,
|
|
|
+ distance_to_ink: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ config: QuadraticCurveFamilyConfig,
|
|
|
+ adaptive_scale: _AdaptiveScale,
|
|
|
+) -> _ScoredCandidate:
|
|
|
+ """Refit one equation from all selected interior and outer fragments."""
|
|
|
+
|
|
|
+ width = max(1.0, bounds.right - bounds.left)
|
|
|
+ height = max(1.0, bounds.bottom - bounds.top)
|
|
|
+ label_by_id = {
|
|
|
+ segment.segment_id: label for label, segment in enumerate(segments, 1)
|
|
|
+ }
|
|
|
+ normalized_pixels: dict[str, np.ndarray] = {}
|
|
|
+ for segment_id in associated_segment_ids:
|
|
|
+ label = label_by_id.get(segment_id)
|
|
|
+ if label is None:
|
|
|
+ continue
|
|
|
+ pixel_y, pixel_x = np.nonzero(segment_labels == label)
|
|
|
+ normalized_pixels[segment_id] = np.column_stack(
|
|
|
+ (
|
|
|
+ (pixel_x.astype(float) - bounds.left) / width,
|
|
|
+ (pixel_y.astype(float) - bounds.top) / height,
|
|
|
+ )
|
|
|
+ )
|
|
|
+ if not normalized_pixels:
|
|
|
+ return candidate
|
|
|
+ fitted_group = _fit_segment_group(
|
|
|
+ tuple(sorted(normalized_pixels)), normalized_pixels,
|
|
|
+ frozenset(candidate.segment_ids), candidate.coefficients,
|
|
|
+ bounds, adaptive_scale.line_width_px, config,
|
|
|
+ )
|
|
|
+ coefficients = fitted_group.coefficients
|
|
|
+ curvature = float(coefficients[0])
|
|
|
+ if curvature < config.minimum_curvature * 0.35:
|
|
|
+ return candidate
|
|
|
+ vertex_x = float(-coefficients[1] / (2.0 * curvature))
|
|
|
+ vertex_y = float(np.polyval(coefficients, vertex_x))
|
|
|
+ vote = _VoteHypothesis(
|
|
|
+ curvature=curvature,
|
|
|
+ vertex_x=vertex_x,
|
|
|
+ vertex_y=vertex_y,
|
|
|
+ vote_score=candidate.vote_score,
|
|
|
+ segment_ids=candidate.segment_ids,
|
|
|
+ proposed_segment_ids=candidate.proposed_segment_ids,
|
|
|
+ )
|
|
|
+ return (
|
|
|
+ _candidate_from_equation(
|
|
|
+ coefficients,
|
|
|
+ vote,
|
|
|
+ distance_to_ink,
|
|
|
+ bounds,
|
|
|
+ config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ or candidate
|
|
|
+ )
|
|
|
+
|
|
|
+
|
|
|
+class QuadraticCurveFamilyExtractor:
|
|
|
+ """Run segment voting, global scoring and non-crossing family selection."""
|
|
|
+
|
|
|
+ def __init__(self, config: QuadraticCurveFamilyConfig | None = None) -> None:
|
|
|
+ self.config = config or QuadraticCurveFamilyConfig()
|
|
|
+
|
|
|
+ @_image_scope
|
|
|
+ def extract(
|
|
|
+ self,
|
|
|
+ candidate_mask: np.ndarray,
|
|
|
+ bounds: PlotBounds,
|
|
|
+ ) -> tuple[QuadraticCurveFamilyResult, QuadraticCurveFamilyDiagnostics]:
|
|
|
+ if candidate_mask.ndim != 2:
|
|
|
+ raise ValueError("candidate_mask必须是单通道图像")
|
|
|
+ mask = np.where(candidate_mask > 0, 255, 0).astype(np.uint8)
|
|
|
+ (
|
|
|
+ skeleton,
|
|
|
+ junction_zone,
|
|
|
+ labels,
|
|
|
+ segments,
|
|
|
+ junction_count,
|
|
|
+ adaptive_scale,
|
|
|
+ ) = _extract_segments(mask, bounds, self.config)
|
|
|
+ if not segments:
|
|
|
+ raise ValueError("没有可用于二次曲线投票的左上至右下短线段")
|
|
|
+ votes = _vote_hypotheses(segments, self.config)
|
|
|
+ if not votes:
|
|
|
+ raise ValueError("短线段没有形成二次曲线参数投票峰值")
|
|
|
+ inverse = np.where(mask > 0, 0, 255).astype(np.uint8)
|
|
|
+ distance_to_ink = cv2.distanceTransform(inverse, cv2.DIST_L2, 3)
|
|
|
+ scored: list[_ScoredCandidate] = []
|
|
|
+ for vote in votes:
|
|
|
+ coefficients = np.asarray(
|
|
|
+ [
|
|
|
+ vote.curvature,
|
|
|
+ -2.0 * vote.curvature * vote.vertex_x,
|
|
|
+ vote.curvature * vote.vertex_x**2 + vote.vertex_y,
|
|
|
+ ],
|
|
|
+ dtype=float,
|
|
|
+ )
|
|
|
+ candidate = _candidate_from_equation(
|
|
|
+ coefficients,
|
|
|
+ vote,
|
|
|
+ distance_to_ink,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ if candidate is None:
|
|
|
+ continue
|
|
|
+ refined = _refine_candidate(
|
|
|
+ candidate,
|
|
|
+ mask,
|
|
|
+ distance_to_ink,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ if refined is not None:
|
|
|
+ scored.append(refined)
|
|
|
+ retained = _deduplicate_candidates(scored, bounds, self.config)
|
|
|
+ if not retained:
|
|
|
+ raise ValueError("没有二次曲线通过全局像素证据评分")
|
|
|
+ selected, family_score = _select_family(retained, bounds, self.config)
|
|
|
+ # Apply endpoint-envelope regularization only when the first global
|
|
|
+ # selection is ambiguous. Clean five-curve charts remain untouched.
|
|
|
+ if len(selected) > 6:
|
|
|
+ family_candidates = _endpoint_frontier_candidates(
|
|
|
+ retained, self.config
|
|
|
+ )
|
|
|
+ if family_candidates:
|
|
|
+ try:
|
|
|
+ selected, family_score = _select_family(
|
|
|
+ family_candidates,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ distinct_endpoints=True,
|
|
|
+ )
|
|
|
+ except ValueError:
|
|
|
+ if self.config.target_family_size is None:
|
|
|
+ raise
|
|
|
+ # Semantic cardinality may exceed the conservative endpoint
|
|
|
+ # envelope on a dense chart. Keep the full candidate result
|
|
|
+ # and let the mandatory post-fit QC report any conflicts.
|
|
|
+ pass
|
|
|
+ # The last stage refits each chosen equation against its own pixel
|
|
|
+ # corridor, after family membership has already been decided.
|
|
|
+ finally_refined: list[_ScoredCandidate] = []
|
|
|
+ for candidate in selected:
|
|
|
+ refined = _refine_candidate(
|
|
|
+ candidate,
|
|
|
+ mask,
|
|
|
+ distance_to_ink,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ finally_refined.append(refined or candidate)
|
|
|
+ selected = sorted(
|
|
|
+ finally_refined, key=lambda item: _lower_endpoint(item)[0]
|
|
|
+ )
|
|
|
+ # Alternate full-plot association and a refit from the chosen fragments.
|
|
|
+ # Refine twice, then reconcile the final membership below, without
|
|
|
+ # letting the full mask pull the curve onto a neighbor.
|
|
|
+ for _ in range(2):
|
|
|
+ selected = _revalidate_vote_seeds(selected, segments, labels, bounds, self.config, adaptive_scale)
|
|
|
+ associated_segment_ids = _associate_segments_to_final_curves(
|
|
|
+ selected,
|
|
|
+ segments,
|
|
|
+ labels,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ selected = [
|
|
|
+ _refit_candidate_from_associated_segments(
|
|
|
+ candidate,
|
|
|
+ associated_ids,
|
|
|
+ segments,
|
|
|
+ labels,
|
|
|
+ distance_to_ink,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ for candidate, associated_ids in zip(
|
|
|
+ selected, associated_segment_ids, strict=True
|
|
|
+ )
|
|
|
+ ]
|
|
|
+ selected = _revalidate_vote_seeds(selected, segments, labels, bounds, self.config, adaptive_scale)
|
|
|
+ associated_segment_ids = _associate_segments_to_final_curves(
|
|
|
+ selected,
|
|
|
+ segments,
|
|
|
+ labels,
|
|
|
+ bounds,
|
|
|
+ self.config,
|
|
|
+ adaptive_scale,
|
|
|
+ )
|
|
|
+ # The exported equation must use the final membership, not the previous
|
|
|
+ # pass's membership (which may have gained/lost a fragment).
|
|
|
+ selected = [
|
|
|
+ _refit_candidate_from_associated_segments(
|
|
|
+ candidate, associated_ids, segments, labels, distance_to_ink,
|
|
|
+ bounds, self.config, adaptive_scale,
|
|
|
+ )
|
|
|
+ for candidate, associated_ids in zip(
|
|
|
+ selected, associated_segment_ids, strict=True
|
|
|
+ )
|
|
|
+ ]
|
|
|
+ family_score = _family_objective(selected, bounds, self.config)
|
|
|
+ curves: list[QuadraticCurve] = []
|
|
|
+ for index, (candidate, associated_ids) in enumerate(
|
|
|
+ zip(selected, associated_segment_ids, strict=True), start=1
|
|
|
+ ):
|
|
|
+ confidence = float(
|
|
|
+ np.clip(
|
|
|
+ 0.18
|
|
|
+ + 0.42 * candidate.support
|
|
|
+ + 0.24 * candidate.direct_support
|
|
|
+ + 0.12 * min(candidate.y_span / 0.55, 1.0)
|
|
|
+ - 0.22 * candidate.maximum_gap,
|
|
|
+ 0.05,
|
|
|
+ 0.99,
|
|
|
+ )
|
|
|
+ )
|
|
|
+ curves.append(
|
|
|
+ QuadraticCurve(
|
|
|
+ curve_id=f"C{index:03d}",
|
|
|
+ coefficients_normalized=tuple(
|
|
|
+ round(float(value), 10) for value in candidate.coefficients
|
|
|
+ ),
|
|
|
+ vertex_normalized=Point(
|
|
|
+ round(float(candidate.vertex[0]), 6),
|
|
|
+ round(float(candidate.vertex[1]), 6),
|
|
|
+ ),
|
|
|
+ x_domain_normalized=tuple(
|
|
|
+ round(float(value), 6) for value in candidate.domain
|
|
|
+ ),
|
|
|
+ pixel_points=tuple(
|
|
|
+ Point(round(float(x), 3), round(float(y), 3))
|
|
|
+ for x, y in candidate.points
|
|
|
+ ),
|
|
|
+ supporting_segment_ids=candidate.segment_ids,
|
|
|
+ associated_segment_ids=associated_ids,
|
|
|
+ x_span_ratio=round(candidate.x_span, 4),
|
|
|
+ y_span_ratio=round(candidate.y_span, 4),
|
|
|
+ pixel_support=round(candidate.support, 4),
|
|
|
+ direct_pixel_support=round(candidate.direct_support, 4),
|
|
|
+ maximum_gap_ratio=round(candidate.maximum_gap, 4),
|
|
|
+ tangent_agreement=round(candidate.tangent_agreement, 4),
|
|
|
+ vote_score=round(candidate.vote_score, 4),
|
|
|
+ global_score=round(candidate.score, 4),
|
|
|
+ confidence=round(confidence, 3),
|
|
|
+ proposed_segment_ids=candidate.proposed_segment_ids,
|
|
|
+ rejected_seed_ids=tuple(sorted(set(candidate.proposed_segment_ids) - set(candidate.segment_ids))),
|
|
|
+ )
|
|
|
+ )
|
|
|
+ result = QuadraticCurveFamilyResult(
|
|
|
+ coordinate_space="rectified_image",
|
|
|
+ plot_bounds=bounds,
|
|
|
+ parameters=asdict(self.config),
|
|
|
+ adaptive_scale=adaptive_scale.to_dict(),
|
|
|
+ skeleton_pixel_count=int(np.count_nonzero(skeleton)),
|
|
|
+ junction_count=junction_count,
|
|
|
+ segment_count=len(segments),
|
|
|
+ vote_hypothesis_count=len(votes),
|
|
|
+ scored_candidate_count=len(retained),
|
|
|
+ curves=tuple(curves),
|
|
|
+ family_score=round(family_score, 4),
|
|
|
+ )
|
|
|
+ diagnostics = QuadraticCurveFamilyDiagnostics(
|
|
|
+ skeleton=skeleton,
|
|
|
+ junction_zone=junction_zone,
|
|
|
+ segment_labels=labels,
|
|
|
+ segments=tuple(segments),
|
|
|
+ adaptive_scale=adaptive_scale.to_dict(),
|
|
|
+ )
|
|
|
+ return result, diagnostics
|