Coverage for lunar_toolkit/selection.py: 93%
325 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Illumination-diverse image selection for shape-from-shading."""
3from __future__ import annotations
5import math
6from dataclasses import dataclass, field
7from typing import TYPE_CHECKING, Any, TypeAlias
9import numpy as np
10import shapely
12from .geometry import AoiLike, footprint_geometry
13from .stereo import angle_between, aoi_plane, azimuth_difference, observation_footprints
15if TYPE_CHECKING: # pragma: no cover
16 from collections.abc import Callable
18 import pandas as pd
20Prefer: TypeAlias = "str | Callable[[pd.DataFrame], pd.Series]"
23def sort_by_sun_azimuth(
24 obs_df: pd.DataFrame,
25 start: float = 0.0,
26 column: str = "sub_solar_azimuth",
27) -> pd.DataFrame:
28 """Sort rows by Sun azimuth, measured clockwise from *start* degrees.
30 Rows with an unknown azimuth go last. The index is preserved.
31 """
32 import pandas as pd
34 az = pd.to_numeric(obs_df[column], errors="coerce")
35 key = (az - start) % 360.0
36 return obs_df.assign(_k=key).sort_values("_k", kind="stable", na_position="last").drop(columns="_k")
39def aoi_coverage(obs_df: pd.DataFrame, aoi: AoiLike) -> pd.Series:
40 """Fraction (0–1) of the AOI covered by each row's footprint.
42 *aoi* is ``{"point": (lat, lon), "radius_km": r}``,
43 ``{"lat": (s, n), "lon": (w, e)}`` or anything
44 [`lunar_toolkit.geometry.load_aoi`][lunar_toolkit.geometry.load_aoi] accepts.
45 """
46 import pandas as pd
48 geoms, target = _aoi_cover(obs_df, aoi)
49 area = target.area
50 vals = [g.intersection(target).area / area if g is not None and area > 0 else 0.0 for g in geoms]
51 return pd.Series(vals, index=obs_df.index, name="aoi_coverage", dtype=float)
54def _aoi_cover(obs_df: pd.DataFrame, aoi: AoiLike) -> tuple[list, Any]:
55 """Per row the footprint in the AOI's tangent plane (``None`` if absent / disjoint), and the AOI."""
56 plane, target, centre, radius = aoi_plane(aoi)
57 rings, centres, radii = observation_footprints(obs_df)
58 near = angle_between(centres, centre[None, :]) < radii + radius + 1e-6
59 geoms = [
60 footprint_geometry(r, plane.ring) if r is not None and near[k] else None for k, r in enumerate(rings)
61 ]
62 return geoms, target
65def _candidates(
66 obs_df: pd.DataFrame,
67 incidence: tuple[float, float] | None,
68 max_emission: float | None,
69 aoi: AoiLike | None,
70 min_aoi_coverage: float,
71) -> pd.DataFrame:
72 import pandas as pd
74 df = obs_df.copy()
75 df["sub_solar_azimuth"] = pd.to_numeric(df["sub_solar_azimuth"], errors="coerce")
76 df = df[df["sub_solar_azimuth"].notna()]
77 if incidence is not None:
78 inc = pd.to_numeric(df["incidence"], errors="coerce")
79 df = df[inc.between(min(incidence), max(incidence))]
80 if max_emission is not None:
81 emi = pd.to_numeric(df["emission"], errors="coerce")
82 df = df[~(emi > max_emission)]
83 if aoi is not None and not df.empty:
84 df["aoi_coverage"] = aoi_coverage(df, aoi)
85 df = df[(df["aoi_coverage"] > 0) & (df["aoi_coverage"] >= min_aoi_coverage)]
86 else:
87 df["aoi_coverage"] = np.nan
88 return df
91def _quality(df: pd.DataFrame, prefer: Prefer) -> pd.Series:
92 """Selection score (higher is better); see
93 [`select_illumination_diverse`][lunar_toolkit.selection.select_illumination_diverse]."""
94 import pandas as pd
96 if callable(prefer):
97 return pd.Series(prefer(df), index=df.index, dtype=float)
98 ps = (
99 pd.to_numeric(df["pixel_scale"], errors="coerce")
100 if "pixel_scale" in df
101 else pd.Series(np.nan, df.index)
102 )
103 res = (ps.min() / ps).fillna(0.5) if ps.notna().any() else pd.Series(0.5, index=df.index)
104 cov = df["aoi_coverage"].fillna(1.0)
105 if prefer == "balanced":
106 return cov * np.sqrt(res)
107 if prefer == "coverage":
108 return cov + 1e-3 * res
109 if prefer == "resolution": 109 ↛ 111line 109 didn't jump to line 111 because the condition on line 109 was always true
110 return res + 1e-3 * cov
111 raise ValueError(f"prefer must be 'balanced', 'coverage', 'resolution' or a callable, not {prefer!r}")
114def _thin_to(df: pd.DataFrame, max_images: int) -> pd.DataFrame:
115 """Drop picks until *max_images* remain, always removing the lower-scored
116 image of the closest pair in azimuth (keeps the azimuth spread)."""
117 df = df.copy()
118 while len(df) > max_images:
119 az = df["sub_solar_azimuth"].to_numpy(float)
120 d = azimuth_difference(az[:, None], az[None, :])
121 np.fill_diagonal(d, np.inf)
122 a, b = np.unravel_index(np.argmin(d), d.shape)
123 sc = df["selection_score"].to_numpy(float)
124 drop = a if sc[a] < sc[b] or (sc[a] == sc[b] and a > b) else b
125 df = df.iloc[np.delete(np.arange(len(df)), drop)]
126 return df
129def select_illumination_diverse(
130 obs_df: pd.DataFrame,
131 azimuth_bin_deg: float = 30.0,
132 incidence: tuple[float, float] | None = None,
133 per_bin: int = 1,
134 max_images: int | None = None,
135 prefer: Prefer = "balanced",
136 aoi: AoiLike | None = None,
137 min_aoi_coverage: float = 0.0,
138 incidence_bin_deg: float | None = None,
139 max_emission: float | None = None,
140 method: str = "bins",
141 azimuth_origin: float = 0.0,
142) -> pd.DataFrame:
143 """Choose images spanning diverse Sun azimuths (for SfS / PSR mapping).
145 Parameters
146 ----------
147 obs_df : DataFrame
148 Observations in the common schema.
149 azimuth_bin_deg : float
150 Bin width (``method="bins"``) or minimum separation
151 (``method="greedy"``) in sub-solar azimuth, degrees.
152 incidence : (min, max), optional
153 Keep only images with incidence in this range (e.g. ``(80, 95)`` at
154 the poles; SfS needs lit but low-Sun images).
155 per_bin : int
156 Images kept per bin. With ``method="greedy"``: at most this many
157 picks may lie within *azimuth_bin_deg* of each other.
158 max_images : int, optional
159 Cap on the result. Excess picks are removed by repeatedly dropping
160 the lower-scored image of the pair closest in azimuth, so the
161 azimuth spread is kept.
162 prefer : {"balanced", "coverage", "resolution"} or callable
163 How to rank images competing for a bin. With ``res`` = finest
164 ``pixel_scale`` among candidates / own ``pixel_scale`` (0.5 if
165 unknown) and ``cov`` = AOI coverage fraction (1 without *aoi*):
166 ``"balanced"`` = ``cov * sqrt(res)``; ``"coverage"`` ranks by
167 coverage then resolution; ``"resolution"`` by resolution then
168 coverage. A callable gets the candidate DataFrame and returns the
169 score (higher is better).
170 aoi : dict, str, Path or Aoi, optional
171 ``{"point": (lat, lon), "radius_km": r}``, ``{"lat": ..., "lon": ...}``
172 or anything [`lunar_toolkit.geometry.load_aoi`][lunar_toolkit.geometry.load_aoi] accepts; images
173 not overlapping it are dropped.
174 min_aoi_coverage : float
175 Minimum fraction of the AOI an image must cover (with *aoi*).
176 incidence_bin_deg : float, optional
177 Also bin by incidence (``method="bins"`` only), giving azimuth ×
178 incidence cells — useful when shadow length diversity matters too.
179 max_emission : float, optional
180 Drop images with larger emission angles (off-nadir views are
181 foreshortened and complicate SfS).
182 method : {"bins", "greedy"}
183 ``"bins"``: fixed bins ``[origin + k*w, origin + (k+1)*w)``, best
184 *per_bin* per bin (two picks in neighbouring bins may still be close
185 in azimuth). ``"greedy"``: candidates in descending score; accept
186 one unless *per_bin* accepted images already lie within
187 *azimuth_bin_deg* of it, so with ``per_bin=1`` all picks are at
188 least *azimuth_bin_deg* apart (the Bertone et al. criterion).
189 azimuth_origin : float
190 Start of the first azimuth bin, degrees.
192 Returns
193 -------
194 DataFrame
195 The selected rows (original index kept) sorted by
196 ``sub_solar_azimuth`` (see [`sort_by_sun_azimuth`][lunar_toolkit.selection.sort_by_sun_azimuth]), with
197 extra columns ``aoi_coverage`` (NaN without *aoi*), ``azimuth_bin`` (bin
198 start, deg; NaN for greedy), ``incidence_bin`` (if binned) and
199 ``selection_score``.
201 Notes
202 -----
203 This uses each image's single scene-centre Sun azimuth; for a large AOI
204 where per-pixel diversity matters, use
205 [`select_for_coverage`][lunar_toolkit.selection.select_for_coverage].
206 """
207 import pandas as pd
209 if azimuth_bin_deg <= 0 or per_bin < 1: 209 ↛ 210line 209 didn't jump to line 210 because the condition on line 209 was never true
210 raise ValueError("azimuth_bin_deg must be > 0 and per_bin >= 1")
211 df = _candidates(obs_df, incidence, max_emission, aoi, min_aoi_coverage)
212 if df.empty: 212 ↛ 213line 212 didn't jump to line 213 because the condition on line 212 was never true
213 return df.assign(azimuth_bin=pd.Series(dtype=float), selection_score=pd.Series(dtype=float))
214 df["selection_score"] = _quality(df, prefer)
215 df = df.sort_values("selection_score", ascending=False, kind="stable")
217 if method == "bins":
218 rel = (df["sub_solar_azimuth"] - azimuth_origin) % 360.0
219 df["azimuth_bin"] = (np.floor(rel / azimuth_bin_deg) * azimuth_bin_deg + azimuth_origin) % 360.0
220 keys = ["azimuth_bin"]
221 if incidence_bin_deg:
222 inc = pd.to_numeric(df["incidence"], errors="coerce")
223 df["incidence_bin"] = np.floor(inc / incidence_bin_deg) * incidence_bin_deg
224 keys.append("incidence_bin")
225 picked = df.groupby(keys, sort=False, dropna=False).head(per_bin)
226 elif method == "greedy": 226 ↛ 237line 226 didn't jump to line 237 because the condition on line 226 was always true
227 df["azimuth_bin"] = np.nan
228 accepted: list[float] = []
229 keep = []
230 for idx, a in zip(df.index, df["sub_solar_azimuth"].to_numpy(float)):
231 near = int(np.sum(azimuth_difference(np.array(accepted), a) < azimuth_bin_deg)) if accepted else 0
232 if near < per_bin:
233 accepted.append(a)
234 keep.append(idx)
235 picked = df.loc[keep]
236 else:
237 raise ValueError(f"method must be 'bins' or 'greedy', not {method!r}")
239 if max_images is not None and len(picked) > max_images:
240 picked = _thin_to(picked, max_images)
241 return sort_by_sun_azimuth(picked)
244def _grid(target, grid_km: float | None, max_cells: int = 20000) -> tuple[np.ndarray, np.ndarray]:
245 """Cell centres ``(x, y)`` of a regular grid over *target* that fall inside it."""
246 x0, y0, x1, y1 = target.bounds
247 if grid_km is None:
248 grid_km = max(max(x1 - x0, y1 - y0) / 100.0, 1e-3)
249 nx = max(int(math.ceil((x1 - x0) / grid_km)), 1)
250 ny = max(int(math.ceil((y1 - y0) / grid_km)), 1)
251 if nx * ny > max_cells: 251 ↛ 252line 251 didn't jump to line 252 because the condition on line 251 was never true
252 f = math.sqrt(nx * ny / max_cells)
253 nx, ny = max(int(nx / f), 1), max(int(ny / f), 1)
254 gx, gy = np.meshgrid(
255 x0 + (np.arange(nx) + 0.5) * (x1 - x0) / nx, y0 + (np.arange(ny) + 0.5) * (y1 - y0) / ny
256 )
257 gx, gy = gx.ravel(), gy.ravel()
258 inside = shapely.intersects_xy(target, gx, gy)
259 return gx[inside], gy[inside]
262def select_for_coverage(
263 obs_df: pd.DataFrame,
264 aoi: AoiLike,
265 azimuth_sep_deg: float = 30.0,
266 views_per_cell: int = 3,
267 incidence: tuple[float, float] | None = None,
268 max_images: int | None = None,
269 prefer: Prefer = "balanced",
270 max_emission: float | None = None,
271 grid_km: float | None = None,
272 min_gain: int = 1,
273) -> pd.DataFrame:
274 """Greedy, per-location illumination-diverse selection over an AOI.
276 The AOI is sampled on a grid (``grid_km`` spacing; default ~100 cells
277 across). Each grid cell keeps the Sun azimuths of the selected images
278 covering it. A candidate *helps* a cell if it covers the cell, the cell
279 has fewer than *views_per_cell* views and the candidate's azimuth is at
280 least *azimuth_sep_deg* from every azimuth already there. At each step
281 the image maximising ``helped_cells * (0.5 + 0.5 * score)`` is taken
282 (score as in [`select_illumination_diverse`][lunar_toolkit.selection.select_illumination_diverse],
283 normalised to 0–1), until no image helps *min_gain* cells or *max_images* is reached.
285 Returns
286 -------
287 DataFrame
288 Selected rows sorted by ``sub_solar_azimuth`` with extra columns
289 ``aoi_coverage``, ``selection_score``, ``selection_order`` (0 = first
290 pick) and ``cells_gained`` (fraction of AOI cells the pick helped).
291 ``attrs["view_fractions"]`` maps a view count *k* to the fraction of
292 AOI cells with exactly *k* non-redundant azimuths
293 (``views_per_cell`` = fully served), and ``attrs["grid_km"]`` /
294 ``attrs["n_cells"]`` describe the grid.
295 """
296 import pandas as pd
298 df = _candidates(obs_df, incidence, max_emission, aoi, 0.0)
299 if df.empty: 299 ↛ 300line 299 didn't jump to line 300 because the condition on line 299 was never true
300 out = df.assign(
301 selection_score=pd.Series(dtype=float),
302 selection_order=pd.Series(dtype=int),
303 cells_gained=pd.Series(dtype=float),
304 )
305 out.attrs["view_fractions"] = {}
306 return out
307 score = _quality(df, prefer).to_numpy(float)
308 smax = np.nanmax(score) if np.isfinite(score).any() else 1.0
309 norm = np.nan_to_num(score / (smax or 1.0), nan=0.0)
310 df["selection_score"] = score
312 geoms, target = _aoi_cover(df, aoi)
313 cx, cy = _grid(target, grid_km)
314 if not len(cx): 314 ↛ 315line 314 didn't jump to line 315 because the condition on line 314 was never true
315 raise ValueError("AOI grid is empty; pass a smaller grid_km")
316 cover = np.stack(
317 [shapely.intersects_xy(g, cx, cy) if g is not None else np.zeros(len(cx), bool) for g in geoms]
318 )
319 n_cells = len(cx)
320 az = df["sub_solar_azimuth"].to_numpy(float)
321 views = np.full((n_cells, views_per_cell), np.nan)
322 count = np.zeros(n_cells, int)
323 available = np.ones(len(df), bool)
324 order: list[int] = []
325 gained: list[float] = []
326 while available.any() and (max_images is None or len(order) < max_images):
327 best, best_val, best_mask = -1, 0.0, None
328 open_cells = count < views_per_cell
329 for k in np.flatnonzero(available):
330 m = cover[k] & open_cells
331 if not m.any():
332 continue
333 d = azimuth_difference(views[m], az[k])
334 ok = np.all(np.isnan(d) | (d >= azimuth_sep_deg), axis=1)
335 g = int(ok.sum())
336 val = g * (0.5 + 0.5 * norm[k])
337 if g >= min_gain and val > best_val:
338 best, best_val = int(k), val
339 best_mask = np.flatnonzero(m)[ok]
340 if best_mask is None:
341 break
342 views[best_mask, count[best_mask]] = az[best]
343 count[best_mask] += 1
344 available[best] = False
345 order.append(best)
346 gained.append(len(best_mask) / n_cells)
348 picked = df.iloc[order].copy()
349 picked["selection_order"] = np.arange(len(order))
350 picked["cells_gained"] = gained
351 out = sort_by_sun_azimuth(picked)
352 out.attrs["view_fractions"] = {int(k): float(np.mean(count == k)) for k in range(views_per_cell + 1)}
353 out.attrs["grid_km"] = float(np.sqrt(target.area / n_cells))
354 out.attrs["n_cells"] = int(n_cells)
355 return out
358def selection_distances(
359 coverage: np.ndarray,
360 azimuth: np.ndarray,
361 weights: tuple[float, float] = (1.0, 1.0),
362 chunk: int = 256,
363) -> np.ndarray:
364 """Image-to-image distance of Bertone et al. (2026, sec. 2) for k-medoids selection.
366 ``d_ij = sqrt((w_a * a_ij)**2 + (w_c * c_ij)**2)`` with ``a_ij`` the
367 circular Sun azimuth (sub-solar longitude) difference and ``c_ij`` the
368 fraction of grid cells covered by exactly one of the two images, each
369 divided by its largest value over all pairs (so both span 0-1).
371 Parameters
372 ----------
373 coverage : (n_images, n_cells) bool array
374 Per-image coverage (or lit) masks on a common grid.
375 azimuth : (n_images,) array
376 Sun azimuths, degrees (finite).
377 weights : (w_a, w_c)
378 Weights of the azimuth and coverage terms.
379 chunk : int
380 Rows per block of the coverage products (memory).
382 Returns
383 -------
384 (n_images, n_images) float array
385 """
386 cov = np.asarray(coverage, bool)
387 az = np.asarray(azimuth, float)
388 n, cells = cov.shape
389 a = azimuth_difference(az[:, None], az[None, :])
390 f = cov.astype(np.float32)
391 counts = f.sum(axis=1)
392 c = np.empty((n, n))
393 for i in range(0, n, chunk):
394 inter = f[i : i + chunk] @ f.T
395 c[i : i + chunk] = (counts[i : i + chunk, None] + counts[None, :] - 2.0 * inter) / max(cells, 1)
396 np.fill_diagonal(c, 0.0)
397 a = a / np.max(a) if np.max(a) > 0 else a
398 c = c / np.max(c) if np.max(c) > 0 else c
399 return np.hypot(weights[0] * a, weights[1] * c)
402def _build(d: np.ndarray, k: int) -> np.ndarray:
403 """PAM BUILD: the most central object, then each object that lowers the total distance most."""
404 med = [int(np.argmin(d.sum(axis=1)))]
405 near = d[:, med[0]].copy()
406 for _ in range(1, k):
407 gain = np.clip(near[:, None] - d, 0.0, None).sum(axis=0)
408 gain[med] = -1.0
409 c = int(np.argmax(gain))
410 med.append(c)
411 near = np.minimum(near, d[:, c])
412 return np.sort(med)
415def kmedoids(
416 dist: np.ndarray,
417 k: int,
418 init: str = "build",
419 seed: int | None = 0,
420 n_init: int = 1,
421 max_iter: int = 100,
422) -> tuple[np.ndarray, np.ndarray]:
423 """Alternating k-medoids (Park & Jun 2009) on a precomputed distance matrix.
425 Assign every object to its nearest medoid, move each medoid to the
426 member with the smallest total distance to its cluster, repeat until the
427 medoids stop changing.
429 Parameters
430 ----------
431 dist : (n, n) array
432 Symmetric distances.
433 k : int
434 Number of clusters (all objects are medoids when ``k >= n``).
435 init : {"build", "park_jun", "random"}
436 ``"build"``: the PAM BUILD step (most central object, then greedy
437 cost reduction; deterministic and robust to duplicate images);
438 ``"park_jun"``: the *k* objects with the smallest
439 ``v_j = sum_i d_ij / sum_l d_il`` (Park & Jun's own start, which
440 favours central objects); ``"random"``: *k* distinct objects drawn
441 with *seed*, best of *n_init* runs.
442 seed : int, optional
443 Random generator seed (``init="random"``; ties are always broken by
444 the lowest index).
445 n_init : int
446 Random restarts (``init="random"``).
447 max_iter : int
448 Iteration cap.
450 Returns
451 -------
452 (medoids, labels)
453 Sorted medoid indices and each object's cluster (index into
454 *medoids*).
455 """
456 d = np.asarray(dist, float)
457 n = len(d)
458 if n == 0 or k <= 0: 458 ↛ 459line 458 didn't jump to line 459 because the condition on line 458 was never true
459 return np.empty(0, int), np.full(n, -1)
460 if k >= n:
461 return np.arange(n), np.arange(n)
462 if init == "build":
463 starts = [_build(d, k)]
464 elif init == "park_jun":
465 row = d.sum(axis=1)
466 v = (d / np.where(row > 0, row, 1.0)[:, None]).sum(axis=0)
467 starts = [np.sort(np.argsort(v, kind="stable")[:k])]
468 elif init == "random":
469 rng = np.random.default_rng(seed)
470 starts = [np.sort(rng.choice(n, k, replace=False)) for _ in range(max(1, n_init))]
471 else:
472 raise ValueError(f"init must be 'build', 'park_jun' or 'random', not {init!r}")
473 best: tuple[float, np.ndarray, np.ndarray] | None = None
474 for med in starts:
475 for _ in range(max_iter): 475 ↛ 486line 475 didn't jump to line 486 because the loop on line 475 didn't complete
476 labels = np.argmin(d[:, med], axis=1)
477 labels[med] = np.arange(k) # duplicates of a medoid stay in their own cluster
478 new = med.copy()
479 for j in range(k):
480 members = np.flatnonzero(labels == j)
481 new[j] = members[np.argmin(d[np.ix_(members, members)].sum(axis=1))]
482 new = np.sort(new)
483 if np.array_equal(new, med):
484 break
485 med = new
486 labels = np.argmin(d[:, med], axis=1)
487 labels[med] = np.arange(k)
488 cost = float(d[np.arange(n), med[labels]].sum())
489 if best is None or cost < best[0] - 1e-12:
490 best = (cost, med, labels)
491 assert best is not None
492 return best[1], best[2]
495@dataclass
496class KMedoidsSelection:
497 """Result of [`select_kmedoids`][lunar_toolkit.selection.select_kmedoids] (indices refer to the input
498 rows).
500 Attributes
501 ----------
502 chosen : list of int
503 Selected images in pick order: medoids first, then the coverage and
504 illumination refinements.
505 stage : list of str
506 Per entry of *chosen*: ``"medoid"``, ``"coverage"`` or ``"illumination"``.
507 medoids, labels : array
508 The k-medoids clustering (``labels[i]`` indexes *medoids*; -1 for
509 images without an azimuth).
510 directions : (n_cells,) int array
511 Distinct illumination directions reached per cell (capped at *views*).
512 uncovered, poorly_lit : float
513 Fractions of the potentially lit cells with no selected image, and
514 with fewer than *views* directions.
515 """
517 chosen: list[int]
518 stage: list[str]
519 medoids: np.ndarray
520 labels: np.ndarray
521 directions: np.ndarray
522 uncovered: float
523 poorly_lit: float
524 history: list[dict[str, Any]] = field(default_factory=list)
527def select_kmedoids(
528 lit: np.ndarray,
529 azimuth: np.ndarray,
530 n_clusters: int = 30,
531 *,
532 coverage: np.ndarray | None = None,
533 views: int = 3,
534 azimuth_sep: float = 30.0,
535 max_uncovered: float = 0.005,
536 max_poorly_lit: float = 0.01,
537 max_images: int | None = None,
538 weights: tuple[float, float] = (1.0, 1.0),
539 init: str = "build",
540 seed: int | None = 0,
541 n_init: int = 1,
542) -> KMedoidsSelection:
543 """k-medoids image selection with coverage and illumination refinement (Bertone et al. 2026, sec. 2).
545 1. Cluster the images with [`kmedoids`][lunar_toolkit.selection.kmedoids] on
546 [`selection_distances`][lunar_toolkit.selection.selection_distances] (Sun azimuth and coverage
547 differences) and take the *n_clusters* medoids.
548 2. Add the image lighting the most still-uncovered cells until fewer
549 than *max_uncovered* of the potentially lit cells (lit in any image)
550 are uncovered.
551 3. Add the image giving the most cells a new direction (Sun azimuth at
552 least *azimuth_sep* from the cell's directions so far, up to *views*)
553 until fewer than *max_poorly_lit* of the potentially lit cells have
554 fewer than *views* directions.
556 Steps 2-3 stop early when no image helps or *max_images* is reached;
557 ties go to the lowest index, so the result is deterministic (for
558 ``init="random"`` given *seed*).
560 Parameters
561 ----------
562 lit : (n_images, n_cells) bool array
563 Lit (valid and above the shadow threshold) cells per image on a
564 common grid, e.g. ``SfsPipeline._lit_masks()``.
565 azimuth : (n_images,) array
566 Sun azimuths in degrees; images with NaN are never selected.
567 n_clusters : int
568 Medoids (the paper's ~30 user-defined images).
569 coverage : (n_images, n_cells) bool array, optional
570 Footprint masks for the distance (default *lit*).
571 views, azimuth_sep : int, float
572 Directions wanted per cell and their minimum separation (3, 30 deg).
573 max_uncovered, max_poorly_lit : float
574 Stopping fractions of steps 2 and 3 (0.5 % and 1 %).
575 max_images : int, optional
576 Cap on the total (medoids are truncated to it first).
577 weights, init, seed, n_init
578 See [`selection_distances`][lunar_toolkit.selection.selection_distances] and
579 [`kmedoids`][lunar_toolkit.selection.kmedoids].
580 """
581 lit = np.asarray(lit, bool)
582 az = np.asarray(azimuth, float)
583 n, cells = lit.shape
584 cov = lit if coverage is None else np.asarray(coverage, bool)
585 valid = np.flatnonzero(np.isfinite(az))
586 labels = np.full(n, -1)
587 if len(valid): 587 ↛ 594line 587 didn't jump to line 594 because the condition on line 587 was always true
588 med, lab = kmedoids(
589 selection_distances(cov[valid], az[valid], weights), n_clusters, init, seed, n_init
590 )
591 medoids = valid[med]
592 labels[valid] = lab
593 else:
594 medoids = np.empty(0, int)
595 cap = n if max_images is None else max(0, max_images)
596 chosen = [int(i) for i in medoids[:cap]]
597 stage = ["medoid"] * len(chosen)
599 potential = lit[valid].any(axis=0) if len(valid) else np.zeros(cells, bool)
600 n_pot = max(int(potential.sum()), 1)
601 dirs = np.full((cells, views), np.nan)
602 count = np.zeros(cells, int)
604 def add(k: int) -> None:
605 idx = np.flatnonzero(lit[k] & (count < views))
606 if len(idx): 606 ↛ exitline 606 didn't return from function 'add' because the condition on line 606 was always true
607 d = azimuth_difference(dirs[idx], az[k])
608 idx = idx[np.all(np.isnan(d) | (d >= azimuth_sep), axis=1)]
609 dirs[idx, count[idx]] = az[k]
610 count[idx] += 1
612 for k in chosen:
613 add(k)
614 covered = lit[chosen].any(axis=0) if chosen else np.zeros(cells, bool)
615 available = np.zeros(n, bool)
616 available[valid] = True
617 available[chosen] = False
618 history: list[dict[str, Any]] = []
620 def fractions() -> tuple[float, float]:
621 return (
622 float(np.sum(potential & ~covered) / n_pot),
623 float(np.sum(potential & (count < views)) / n_pot),
624 )
626 while len(chosen) < cap and fractions()[0] >= max_uncovered and available.any():
627 cand = np.flatnonzero(available)
628 gain = (lit[cand] & ~covered & potential).sum(axis=1)
629 if gain.max() <= 0: 629 ↛ 630line 629 didn't jump to line 630 because the condition on line 629 was never true
630 break
631 k = int(cand[np.argmax(gain)])
632 chosen.append(k)
633 stage.append("coverage")
634 covered |= lit[k]
635 available[k] = False
636 add(k)
637 history.append({"image": k, "stage": "coverage", "gain": int(gain.max())})
639 while len(chosen) < cap and fractions()[1] >= max_poorly_lit and available.any():
640 open_cells = potential & (count < views)
641 best, best_gain = -1, 0
642 for j in np.flatnonzero(available):
643 idx = np.flatnonzero(lit[j] & open_cells)
644 if len(idx) <= best_gain:
645 continue
646 d = azimuth_difference(dirs[idx], az[j])
647 g = int(np.all(np.isnan(d) | (d >= azimuth_sep), axis=1).sum())
648 if g > best_gain:
649 best, best_gain = int(j), g
650 if best < 0: 650 ↛ 651line 650 didn't jump to line 651 because the condition on line 650 was never true
651 break
652 chosen.append(best)
653 stage.append("illumination")
654 covered |= lit[best]
655 available[best] = False
656 add(best)
657 history.append({"image": best, "stage": "illumination", "gain": best_gain})
659 unc, poor = fractions()
660 return KMedoidsSelection(chosen, stage, medoids, labels, count, unc, poor, history)
663def coverage_matrix(
664 obs_df: pd.DataFrame, aoi: AoiLike, grid_km: float | None = None
665) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
666 """Footprint coverage of the AOI grid cells: ``(cover[n_images, n_cells], x, y)``.
668 The grid is the one [`select_for_coverage`][lunar_toolkit.selection.select_for_coverage] uses (cell
669 centres in the AOI's tangent plane, km).
670 """
671 geoms, target = _aoi_cover(obs_df, aoi)
672 cx, cy = _grid(target, grid_km)
673 cover = (
674 np.stack(
675 [shapely.intersects_xy(g, cx, cy) if g is not None else np.zeros(len(cx), bool) for g in geoms]
676 )
677 if geoms
678 else np.zeros((0, len(cx)), bool)
679 )
680 return cover, cx, cy
683def select_kmedoids_footprints(
684 obs_df: pd.DataFrame,
685 aoi: AoiLike,
686 n_clusters: int = 30,
687 incidence: tuple[float, float] | None = None,
688 max_emission: float | None = None,
689 grid_km: float | None = None,
690 **options: Any,
691) -> pd.DataFrame:
692 """[`select_kmedoids`][lunar_toolkit.selection.select_kmedoids] on catalogue footprints over an AOI
693 (footprint = lit, no shadows).
695 *options* go to [`select_kmedoids`][lunar_toolkit.selection.select_kmedoids]. Returns the selected rows
696 sorted by ``sub_solar_azimuth`` with ``selection_order``,
697 ``selection_stage`` and ``cluster`` columns; ``attrs`` holds
698 ``uncovered``, ``poorly_lit`` and ``n_cells``.
699 """
700 import pandas as pd
702 df = _candidates(obs_df, incidence, max_emission, aoi, 0.0)
703 if df.empty: 703 ↛ 704line 703 didn't jump to line 704 because the condition on line 703 was never true
704 out = df.assign(
705 selection_order=pd.Series(dtype=int),
706 selection_stage=pd.Series(dtype=str),
707 cluster=pd.Series(dtype=int),
708 )
709 out.attrs.update(uncovered=1.0, poorly_lit=1.0, n_cells=0)
710 return out
711 cover, cx, _ = coverage_matrix(df, aoi, grid_km)
712 res = select_kmedoids(cover, df["sub_solar_azimuth"].to_numpy(float), n_clusters, **options)
713 picked = df.iloc[res.chosen].copy()
714 picked["selection_order"] = np.arange(len(res.chosen))
715 picked["selection_stage"] = res.stage
716 picked["cluster"] = res.labels[res.chosen]
717 out = sort_by_sun_azimuth(picked)
718 out.attrs.update(uncovered=res.uncovered, poorly_lit=res.poorly_lit, n_cells=int(len(cx)))
719 return out