warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Extend GiST index to 2-D: altitude + inclination

Browse Source 
The 1-D altitude-band index only pruned ~25% of the 22k satellite
catalog (eliminates MEO/GEO/HEO but 75% is LEO).  Adding inclination
as a second indexed dimension prunes an additional ~40% of remaining
candidates — objects in equatorial or low-inclination orbits that
geometrically cannot pass over the observer's latitude.

Key changes:
- tle_alt_range (16 bytes) → tle_orbital_key (32 bytes) with
  inc_low/inc_high fields
- All 8 GiST support functions updated for 2-D bounding boxes
- Penalty uses margin (half-perimeter) not area to avoid degeneracy
  when leaf entries have zero-width inclination ranges
- Picksplit selects split dimension by normalized spread
- && operator now checks altitude AND inclination overlap
- <-> operator remains altitude-only (conjunction screening is
  altitude-dominant)
- SQL operator comments updated for 2-D semantics
- Test adds Equatorial-LEO satellite at ISS altitude but 5° inclination
  to validate inclination-based pruning
This commit is contained in:
Ryan Malloy 2026-02-15 18:10:19 -07:00
parent 5552bf3280
commit a792e7e083
4 changed files with 284 additions and 194 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
sql/pg_orbit--0.1.0.sql
View File

@ -444,11 +444,11 @@ COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS
-- GiST operator support functions -- GiST operator support functions
-- ============================================================ -- ============================================================
-- Overlap operator: do altitude bands intersect? -- Overlap operator: do orbital keys overlap in altitude AND inclination?
CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
-- Altitude distance operator -- Altitude distance operator (altitude-only, for KNN ordering)
CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8 CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
@ -461,7 +461,7 @@ CREATE OPERATOR && (
JOIN = areajoinsel JOIN = areajoinsel
); );
COMMENT ON OPERATOR && (tle, tle) IS 'Altitude band overlap — necessary condition for conjunction'; COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction';
CREATE OPERATOR <-> ( CREATE OPERATOR <-> (
LEFTARG = tle, LEFTARG = tle,
@ -470,11 +470,11 @@ CREATE OPERATOR <-> (
COMMUTATOR = <-> COMMUTATOR = <->
); );
COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping)'; COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping). Altitude-only — does not account for inclination. Use && for 2-D filtering.';
-- ============================================================ -- ============================================================
-- GiST operator class for altitude-band indexing -- GiST operator class for 2-D orbital indexing (altitude + inclination)
-- ============================================================ -- ============================================================
-- GiST internal support functions -- GiST internal support functions

371
src/gist_tle.c
View File

@ -1,14 +1,18 @@
/* /*
* gist_tle.c -- GiST operator class for altitude-band indexing on TLE * gist_tle.c -- GiST operator class for 2-D orbital indexing on TLE
* *
* Every TLE defines an orbit whose perigee and apogee altitudes form a * Every TLE defines an orbit whose perigee/apogee altitudes and
* 1-D range (the "altitude band"). Two orbits can only be in proximity * inclination form a 2-D bounding box in (altitude, inclination) space.
* if their altitude bands overlap -- a necessary but not sufficient * Two orbits can only be in proximity if their altitude bands overlap
* condition for conjunction. * AND their inclination ranges overlap -- both necessary but not
* sufficient conditions for conjunction.
* *
* The GiST index stores [perigee_km, apogee_km] ranges as internal * The GiST index stores [alt_low, alt_high] x [inc_low, inc_high] keys,
* keys, enabling fast coarse filtering. The && (overlap) operator is * enabling fast coarse filtering. For leaf entries inc_low == inc_high
* always rechecked: real conjunction screening requires propagation. * (exact value); internal nodes widen to the bounding box of children.
*
* The && (overlap) operator is always rechecked: real conjunction
* screening requires propagation.
* *
* Semi-major axis from Kepler's third law using WGS-72 constants: * Semi-major axis from Kepler's third law using WGS-72 constants:
* a = (KE / n)^(2/3) [earth radii] * a = (KE / n)^(2/3) [earth radii]
@ -37,30 +41,37 @@ PG_FUNCTION_INFO_V1(gist_tle_picksplit);
PG_FUNCTION_INFO_V1(gist_tle_same); PG_FUNCTION_INFO_V1(gist_tle_same);
PG_FUNCTION_INFO_V1(gist_tle_distance); PG_FUNCTION_INFO_V1(gist_tle_distance);
/* Floating-point comparison tolerance (km) */ /* Floating-point comparison tolerance (km and radians) */
#define ALT_EPSILON 1.0e-9 #define KEY_EPSILON 1.0e-9
/* /*
* Altitude band extracted from a TLE's mean elements. * 2-D orbital key extracted from a TLE's mean elements.
* Altitude band (perigee/apogee) plus inclination range.
* This is the GiST internal key -- much cheaper to compare * This is the GiST internal key -- much cheaper to compare
* than propagating two full state vectors. * than propagating two full state vectors.
*
* For leaf entries: inc_low == inc_high (exact inclination).
* For internal nodes: bounding box over all children.
*/ */
typedef struct tle_alt_range typedef struct tle_orbital_key
{ {
double low; /* perigee altitude, km */ double alt_low; /* perigee altitude, km */
double high; /* apogee altitude, km */ double alt_high; /* apogee altitude, km */
} tle_alt_range; double inc_low; /* inclination lower bound, radians */
double inc_high; /* inclination upper bound, radians */
} tle_orbital_key;
/* ---------------------------------------------------------------- /* ----------------------------------------------------------------
* tle_to_alt_range -- compute [perigee, apogee] from mean elements * tle_to_orbital_key -- compute [perigee, apogee] x [inc, inc]
* *
* Uses WGS-72 KE and AE (the only constants valid for SGP4 elements). * Uses WGS-72 KE and AE (the only constants valid for SGP4 elements).
* Degenerate TLEs with n <= 0 map to a zero-width range at 0 km. * Degenerate TLEs with n <= 0 map to zero-width ranges at 0.
* Inclination is stored in radians (same as pg_tle).
* ---------------------------------------------------------------- * ----------------------------------------------------------------
*/ */
static void static void
tle_to_alt_range(const pg_tle *tle, tle_alt_range *range) tle_to_orbital_key(const pg_tle *tle, tle_orbital_key *key)
{ {
double n = tle->mean_motion; /* rad/min */ double n = tle->mean_motion; /* rad/min */
double e = tle->eccentricity; double e = tle->eccentricity;
@ -68,74 +79,91 @@ tle_to_alt_range(const pg_tle *tle, tle_alt_range *range)
if (n <= 0.0) if (n <= 0.0)
{ {
range->low = 0.0; key->alt_low = 0.0;
range->high = 0.0; key->alt_high = 0.0;
key->inc_low = 0.0;
key->inc_high = 0.0;
return; return;
} }
a_er = pow(WGS72_KE / n, 2.0 / 3.0); a_er = pow(WGS72_KE / n, 2.0 / 3.0);
range->low = a_er * (1.0 - e) * WGS72_AE - WGS72_AE; key->alt_low = a_er * (1.0 - e) * WGS72_AE - WGS72_AE;
range->high = a_er * (1.0 + e) * WGS72_AE - WGS72_AE; key->alt_high = a_er * (1.0 + e) * WGS72_AE - WGS72_AE;
/* Guard against numerical inversion from near-zero eccentricity */ /* Guard against numerical inversion from near-zero eccentricity */
if (range->low > range->high) if (key->alt_low > key->alt_high)
{ {
double tmp = range->low; double tmp = key->alt_low;
range->low = range->high; key->alt_low = key->alt_high;
range->high = tmp; key->alt_high = tmp;
} }
/* Leaf entry: exact inclination (radians) */
key->inc_low = tle->inclination;
key->inc_high = tle->inclination;
} }
/* ---------------------------------------------------------------- /* ----------------------------------------------------------------
* range_overlaps -- do two altitude bands share any interval? * key_overlaps -- do two orbital keys overlap in BOTH dimensions?
*
* Altitude bands AND inclination ranges must both overlap.
* ---------------------------------------------------------------- * ----------------------------------------------------------------
*/ */
static inline bool static inline bool
range_overlaps(const tle_alt_range *a, const tle_alt_range *b) key_overlaps(const tle_orbital_key *a, const tle_orbital_key *b)
{ {
return (a->low <= b->high) && (b->low <= a->high); return (a->alt_low <= b->alt_high) && (b->alt_low <= a->alt_high)
&& (a->inc_low <= b->inc_high) && (b->inc_low <= a->inc_high);
} }
/* ---------------------------------------------------------------- /* ----------------------------------------------------------------
* range_contains -- does outer fully contain inner? * key_contains -- does outer fully contain inner in both dimensions?
* ---------------------------------------------------------------- * ----------------------------------------------------------------
*/ */
static inline bool static inline bool
range_contains(const tle_alt_range *outer, const tle_alt_range *inner) key_contains(const tle_orbital_key *outer, const tle_orbital_key *inner)
{ {
return (outer->low <= inner->low) && (inner->high <= outer->high); return (outer->alt_low <= inner->alt_low)
&& (inner->alt_high <= outer->alt_high)
&& (outer->inc_low <= inner->inc_low)
&& (inner->inc_high <= outer->inc_high);
} }
/* ---------------------------------------------------------------- /* ----------------------------------------------------------------
* range_merge -- expand dst to encompass src * key_merge -- expand dst bounding box to encompass src
* ---------------------------------------------------------------- * ----------------------------------------------------------------
*/ */
static inline void static inline void
range_merge(tle_alt_range *dst, const tle_alt_range *src) key_merge(tle_orbital_key *dst, const tle_orbital_key *src)
{ {
if (src->low < dst->low) if (src->alt_low < dst->alt_low)
dst->low = src->low; dst->alt_low = src->alt_low;
if (src->high > dst->high) if (src->alt_high > dst->alt_high)
dst->high = src->high; dst->alt_high = src->alt_high;
if (src->inc_low < dst->inc_low)
dst->inc_low = src->inc_low;
if (src->inc_high > dst->inc_high)
dst->inc_high = src->inc_high;
} }
/* ---------------------------------------------------------------- /* ----------------------------------------------------------------
* range_separation -- minimum gap between two non-overlapping ranges * alt_separation -- minimum altitude gap between two keys
* *
* Returns 0 if the ranges overlap. * Returns 0 if the altitude bands overlap.
* Used for KNN distance (altitude-dominant ordering).
* ---------------------------------------------------------------- * ----------------------------------------------------------------
*/ */
static inline double static inline double
range_separation(const tle_alt_range *a, const tle_alt_range *b) alt_separation(const tle_orbital_key *a, const tle_orbital_key *b)
{ {
if (a->high < b->low) if (a->alt_high < b->alt_low)
return b->low - a->high; return b->alt_low - a->alt_high;
if (b->high < a->low) if (b->alt_high < a->alt_low)
return a->low - b->high; return a->alt_low - b->alt_high;
return 0.0; return 0.0;
} }
@ -149,8 +177,8 @@ range_separation(const tle_alt_range *a, const tle_alt_range *b)
/* /*
* tle_overlap(tle, tle) -> bool [the && operator] * tle_overlap(tle, tle) -> bool [the && operator]
* *
* True if the altitude bands of two TLEs share any interval. * True if two TLEs overlap in both altitude AND inclination.
* This is a fast pre-filter: overlapping bands are necessary * This is a fast pre-filter: overlapping keys are necessary
* but not sufficient for actual conjunction. * but not sufficient for actual conjunction.
*/ */
Datum Datum
@ -158,12 +186,12 @@ tle_overlap(PG_FUNCTION_ARGS)
{ {
pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0); pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0);
pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1); pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1);
tle_alt_range ra, rb; tle_orbital_key ka, kb;
tle_to_alt_range(a, &ra); tle_to_orbital_key(a, &ka);
tle_to_alt_range(b, &rb); tle_to_orbital_key(b, &kb);
PG_RETURN_BOOL(range_overlaps(&ra, &rb)); PG_RETURN_BOOL(key_overlaps(&ka, &kb));
} }
@ -174,18 +202,21 @@ tle_overlap(PG_FUNCTION_ARGS)
* overlap. This is not the physical distance between the objects -- * overlap. This is not the physical distance between the objects --
* it is the gap between their orbital shells, useful for ordering * it is the gap between their orbital shells, useful for ordering
* nearest-neighbor queries without propagation. * nearest-neighbor queries without propagation.
*
* Altitude-only: inclination weighting adds complexity without
* meaningful benefit for KNN conjunction screening.
*/ */
Datum Datum
tle_alt_distance(PG_FUNCTION_ARGS) tle_alt_distance(PG_FUNCTION_ARGS)
{ {
pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0); pg_tle *a = (pg_tle *) PG_GETARG_POINTER(0);
pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1); pg_tle *b = (pg_tle *) PG_GETARG_POINTER(1);
tle_alt_range ra, rb; tle_orbital_key ka, kb;
tle_to_alt_range(a, &ra); tle_to_orbital_key(a, &ka);
tle_to_alt_range(b, &rb); tle_to_orbital_key(b, &kb);
PG_RETURN_FLOAT8(range_separation(&ra, &rb)); PG_RETURN_FLOAT8(alt_separation(&ka, &kb));
} }
@ -196,10 +227,10 @@ tle_alt_distance(PG_FUNCTION_ARGS)
/* /*
* gist_tle_compress -- extract altitude range from a leaf TLE * gist_tle_compress -- extract orbital key from a leaf TLE
* *
* Leaf entries carry the full pg_tle; we compress to tle_alt_range. * Leaf entries carry the full pg_tle; we compress to tle_orbital_key.
* Internal entries are already tle_alt_range from union operations. * Internal entries are already tle_orbital_key from union operations.
*/ */
Datum Datum
gist_tle_compress(PG_FUNCTION_ARGS) gist_tle_compress(PG_FUNCTION_ARGS)
@ -209,18 +240,18 @@ gist_tle_compress(PG_FUNCTION_ARGS)
if (entry->leafkey) if (entry->leafkey)
{ {
pg_tle *tle = (pg_tle *) DatumGetPointer(entry->key); pg_tle *tle = (pg_tle *) DatumGetPointer(entry->key);
tle_alt_range *range = (tle_alt_range *) palloc(sizeof(tle_alt_range)); tle_orbital_key *key = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
tle_to_alt_range(tle, range); tle_to_orbital_key(tle, key);
retval = (GISTENTRY *) palloc(sizeof(GISTENTRY)); retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
gistentryinit(*retval, PointerGetDatum(range), gistentryinit(*retval, PointerGetDatum(key),
entry->rel, entry->page, entry->offset, false); entry->rel, entry->page, entry->offset, false);
} }
else else
{ {
/* Internal node: already a tle_alt_range */ /* Internal node: already a tle_orbital_key */
retval = entry; retval = entry;
} }
@ -241,48 +272,44 @@ gist_tle_decompress(PG_FUNCTION_ARGS)
/* /*
* gist_tle_consistent -- can this subtree contain matches for the query? * gist_tle_consistent -- can this subtree contain matches for the query?
* *
* Strategy RTOverlapStrategyNumber (&&): altitude bands must overlap. * Checks overlap in both altitude AND inclination dimensions.
* Always sets recheck = true because altitude overlap is only a necessary * Always sets recheck = true because 2-D overlap is only a necessary
* condition -- the real conjunction test requires propagation. * condition -- the real conjunction test requires propagation.
*/ */
Datum Datum
gist_tle_consistent(PG_FUNCTION_ARGS) gist_tle_consistent(PG_FUNCTION_ARGS)
{ {
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1); pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1);
StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2); StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
/* arg 3 is the query type OID, unused */ /* arg 3 is the query type OID, unused */
bool *recheck = (bool *) PG_GETARG_POINTER(4); bool *recheck = (bool *) PG_GETARG_POINTER(4);
tle_alt_range *key = (tle_alt_range *) DatumGetPointer(entry->key); tle_orbital_key *key = (tle_orbital_key *) DatumGetPointer(entry->key);
tle_alt_range query_range; tle_orbital_key query_key;
bool result; bool result;
tle_to_alt_range(query, &query_range); tle_to_orbital_key(query, &query_key);
/*
* Altitude overlap is necessary, not sufficient.
* The actual operator (if exact) would need propagation, so always recheck.
*/
*recheck = true; *recheck = true;
switch (strategy) switch (strategy)
{ {
case RTOverlapStrategyNumber: /* && */ case RTOverlapStrategyNumber: /* && */
result = range_overlaps(key, &query_range); result = key_overlaps(key, &query_key);
break; break;
case RTContainedByStrategyNumber: /* <@ */ case RTContainedByStrategyNumber: /* <@ */
if (GIST_LEAF(entry)) if (GIST_LEAF(entry))
result = range_contains(&query_range, key); result = key_contains(&query_key, key);
else else
result = range_overlaps(key, &query_range); result = key_overlaps(key, &query_key);
break; break;
case RTContainsStrategyNumber: /* @> */ case RTContainsStrategyNumber: /* @> */
if (GIST_LEAF(entry)) if (GIST_LEAF(entry))
result = range_contains(key, &query_range); result = key_contains(key, &query_key);
else else
result = range_overlaps(key, &query_range); result = key_overlaps(key, &query_key);
break; break;
default: default:
@ -297,31 +324,30 @@ gist_tle_consistent(PG_FUNCTION_ARGS)
/* /*
* gist_tle_union -- compute bounding altitude range for a set of entries * gist_tle_union -- compute 2-D bounding box for a set of entries
* *
* The union of N altitude ranges is simply [min(low), max(high)]. * The union is [min(alt_low), max(alt_high)] x [min(inc_low), max(inc_high)].
*/ */
Datum Datum
gist_tle_union(PG_FUNCTION_ARGS) gist_tle_union(PG_FUNCTION_ARGS)
{ {
GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0); GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
int *sizep = (int *) PG_GETARG_POINTER(1); int *sizep = (int *) PG_GETARG_POINTER(1);
int i; int i;
tle_alt_range *result; tle_orbital_key *result;
tle_alt_range *cur; tle_orbital_key *cur;
result = (tle_alt_range *) palloc(sizeof(tle_alt_range)); result = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[0].key); cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[0].key);
result->low = cur->low; *result = *cur;
result->high = cur->high;
for (i = 1; i < entryvec->n; i++) for (i = 1; i < entryvec->n; i++)
{ {
cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[i].key); cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
range_merge(result, cur); key_merge(result, cur);
} }
*sizep = sizeof(tle_alt_range); *sizep = sizeof(tle_orbital_key);
PG_RETURN_POINTER(result); PG_RETURN_POINTER(result);
} }
@ -330,26 +356,28 @@ gist_tle_union(PG_FUNCTION_ARGS)
/* /*
* gist_tle_penalty -- cost of inserting a new entry into an existing subtree * gist_tle_penalty -- cost of inserting a new entry into an existing subtree
* *
* Penalty is the increase in altitude span (km) that results from * Uses margin (half-perimeter) instead of area. Leaf entries have
* expanding the subtree's bounding range to include the new entry. * inc_low == inc_high, giving zero area -- area-based penalty would
* A good penalty function keeps the tree balanced and minimizes * degenerate to 0 for every insertion, making the tree random.
* unnecessary range expansion. * Margin remains non-zero for degenerate (zero-width) bounding boxes.
*/ */
Datum Datum
gist_tle_penalty(PG_FUNCTION_ARGS) gist_tle_penalty(PG_FUNCTION_ARGS)
{ {
GISTENTRY *origentry = (GISTENTRY *) PG_GETARG_POINTER(0); GISTENTRY *origentry = (GISTENTRY *) PG_GETARG_POINTER(0);
GISTENTRY *newentry = (GISTENTRY *) PG_GETARG_POINTER(1); GISTENTRY *newentry = (GISTENTRY *) PG_GETARG_POINTER(1);
float *penalty = (float *) PG_GETARG_POINTER(2); float *penalty = (float *) PG_GETARG_POINTER(2);
tle_alt_range *orig = (tle_alt_range *) DatumGetPointer(origentry->key); tle_orbital_key *orig = (tle_orbital_key *) DatumGetPointer(origentry->key);
tle_alt_range *add = (tle_alt_range *) DatumGetPointer(newentry->key); tle_orbital_key *add = (tle_orbital_key *) DatumGetPointer(newentry->key);
double orig_span; double orig_margin;
double merged_span; double merged_margin;
orig_span = orig->high - orig->low; orig_margin = (orig->alt_high - orig->alt_low)
merged_span = fmax(orig->high, add->high) - fmin(orig->low, add->low); + (orig->inc_high - orig->inc_low);
merged_margin = (fmax(orig->alt_high, add->alt_high) - fmin(orig->alt_low, add->alt_low))
+ (fmax(orig->inc_high, add->inc_high) - fmin(orig->inc_low, add->inc_low));
*penalty = (float)(merged_span - orig_span); *penalty = (float)(merged_margin - orig_margin);
PG_RETURN_POINTER(penalty); PG_RETURN_POINTER(penalty);
} }
@ -357,19 +385,19 @@ gist_tle_penalty(PG_FUNCTION_ARGS)
/* /*
* Comparison callback for qsort in picksplit. * Comparison callback for qsort in picksplit.
* Sorts entries by the midpoint of their altitude range. * Sorts entries by a sort key chosen at call time.
*/ */
typedef struct typedef struct
{ {
int index; /* position in the original entry vector */ int index; /* position in the original entry vector */
double midpoint; /* (low + high) / 2 */ double sortval; /* midpoint in the chosen dimension */
} picksplit_item; } picksplit_item;
static int static int
picksplit_cmp(const void *a, const void *b) picksplit_cmp(const void *a, const void *b)
{ {
double ma = ((const picksplit_item *) a)->midpoint; double ma = ((const picksplit_item *) a)->sortval;
double mb = ((const picksplit_item *) b)->midpoint; double mb = ((const picksplit_item *) b)->sortval;
if (ma < mb) if (ma < mb)
return -1; return -1;
@ -382,29 +410,66 @@ picksplit_cmp(const void *a, const void *b)
/* /*
* gist_tle_picksplit -- split an overfull page into two groups * gist_tle_picksplit -- split an overfull page into two groups
* *
* Strategy: sort entries by altitude midpoint, split at the median. * Standard R-tree approach: compute spread in both dimensions, split
* This keeps nearby altitude bands in the same subtree, which is * along whichever dimension has the greater spread. This prevents
* optimal for a 1-D range index. * the tree from becoming biased toward one dimension.
*/ */
Datum Datum
gist_tle_picksplit(PG_FUNCTION_ARGS) gist_tle_picksplit(PG_FUNCTION_ARGS)
{ {
GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0); GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
GIST_SPLITVEC *splitvec = (GIST_SPLITVEC *) PG_GETARG_POINTER(1); GIST_SPLITVEC *splitvec = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
int nentries = entryvec->n; int nentries = entryvec->n;
picksplit_item *items; picksplit_item *items;
tle_alt_range *left_union; tle_orbital_key *left_union;
tle_alt_range *right_union; tle_orbital_key *right_union;
tle_alt_range *cur; tle_orbital_key *cur;
int split_at; int split_at;
int i; int i;
double alt_min, alt_max, inc_min, inc_max;
double alt_spread, inc_spread;
bool split_on_alt;
/* First pass: compute spread in both dimensions */
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[0].key);
alt_min = (cur->alt_low + cur->alt_high) / 2.0;
alt_max = alt_min;
inc_min = (cur->inc_low + cur->inc_high) / 2.0;
inc_max = inc_min;
for (i = 1; i < nentries; i++)
{
double alt_mid, inc_mid;
cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
alt_mid = (cur->alt_low + cur->alt_high) / 2.0;
inc_mid = (cur->inc_low + cur->inc_high) / 2.0;
if (alt_mid < alt_min) alt_min = alt_mid;
if (alt_mid > alt_max) alt_max = alt_mid;
if (inc_mid < inc_min) inc_min = inc_mid;
if (inc_mid > inc_max) inc_max = inc_mid;
}
/*
* Normalize spreads so they're comparable. Altitude is km above
* surface (range ~0-35786), inclination in radians (range 0-pi).
* Divide each by its domain width to get a 0-1 fraction.
*/
alt_spread = (alt_max - alt_min) / 35786.0; /* GEO altitude above surface */
inc_spread = (inc_max - inc_min) / M_PI;
split_on_alt = (alt_spread >= inc_spread);
/* Second pass: compute sort values in the chosen dimension */
items = (picksplit_item *) palloc(sizeof(picksplit_item) * nentries); items = (picksplit_item *) palloc(sizeof(picksplit_item) * nentries);
for (i = 0; i < nentries; i++) for (i = 0; i < nentries; i++)
{ {
cur = (tle_alt_range *) DatumGetPointer(entryvec->vector[i].key); cur = (tle_orbital_key *) DatumGetPointer(entryvec->vector[i].key);
items[i].index = i; items[i].index = i;
items[i].midpoint = (cur->low + cur->high) / 2.0; if (split_on_alt)
items[i].sortval = (cur->alt_low + cur->alt_high) / 2.0;
else
items[i].sortval = (cur->inc_low + cur->inc_high) / 2.0;
} }
qsort(items, nentries, sizeof(picksplit_item), picksplit_cmp); qsort(items, nentries, sizeof(picksplit_item), picksplit_cmp);
@ -417,39 +482,37 @@ gist_tle_picksplit(PG_FUNCTION_ARGS)
splitvec->spl_nleft = 0; splitvec->spl_nleft = 0;
splitvec->spl_nright = 0; splitvec->spl_nright = 0;
/* Compute union ranges and assign entries */ /* Compute union keys and assign entries */
left_union = (tle_alt_range *) palloc(sizeof(tle_alt_range)); left_union = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
right_union = (tle_alt_range *) palloc(sizeof(tle_alt_range)); right_union = (tle_orbital_key *) palloc(sizeof(tle_orbital_key));
/* Seed the unions from the first entry in each half */ /* Seed the unions from the first entry in each half */
cur = (tle_alt_range *) DatumGetPointer( cur = (tle_orbital_key *) DatumGetPointer(
entryvec->vector[items[0].index].key); entryvec->vector[items[0].index].key);
left_union->low = cur->low; *left_union = *cur;
left_union->high = cur->high;
cur = (tle_alt_range *) DatumGetPointer( cur = (tle_orbital_key *) DatumGetPointer(
entryvec->vector[items[split_at].index].key); entryvec->vector[items[split_at].index].key);
right_union->low = cur->low; *right_union = *cur;
right_union->high = cur->high;
for (i = 0; i < nentries; i++) for (i = 0; i < nentries; i++)
{ {
int idx = items[i].index; int idx = items[i].index;
cur = (tle_alt_range *) DatumGetPointer( cur = (tle_orbital_key *) DatumGetPointer(
entryvec->vector[idx].key); entryvec->vector[idx].key);
if (i < split_at) if (i < split_at)
{ {
splitvec->spl_left[splitvec->spl_nleft++] = splitvec->spl_left[splitvec->spl_nleft++] =
(OffsetNumber)(idx + 1); /* 1-based */ (OffsetNumber)(idx + 1); /* 1-based */
range_merge(left_union, cur); key_merge(left_union, cur);
} }
else else
{ {
splitvec->spl_right[splitvec->spl_nright++] = splitvec->spl_right[splitvec->spl_nright++] =
(OffsetNumber)(idx + 1); (OffsetNumber)(idx + 1);
range_merge(right_union, cur); key_merge(right_union, cur);
} }
} }
@ -465,19 +528,20 @@ gist_tle_picksplit(PG_FUNCTION_ARGS)
/* /*
* gist_tle_same -- equality test on compressed keys * gist_tle_same -- equality test on compressed keys
* *
* Two altitude ranges are "same" if both endpoints match within * Two orbital keys are "same" if all four bounds match within
* a small tolerance. This lets the GiST machinery detect * tolerance. This lets the GiST machinery detect duplicate keys.
* duplicate keys.
*/ */
Datum Datum
gist_tle_same(PG_FUNCTION_ARGS) gist_tle_same(PG_FUNCTION_ARGS)
{ {
tle_alt_range *a = (tle_alt_range *) PG_GETARG_POINTER(0); tle_orbital_key *a = (tle_orbital_key *) PG_GETARG_POINTER(0);
tle_alt_range *b = (tle_alt_range *) PG_GETARG_POINTER(1); tle_orbital_key *b = (tle_orbital_key *) PG_GETARG_POINTER(1);
bool *result = (bool *) PG_GETARG_POINTER(2); bool *result = (bool *) PG_GETARG_POINTER(2);
*result = (fabs(a->low - b->low) < ALT_EPSILON && *result = (fabs(a->alt_low - b->alt_low) < KEY_EPSILON
fabs(a->high - b->high) < ALT_EPSILON); && fabs(a->alt_high - b->alt_high) < KEY_EPSILON
&& fabs(a->inc_low - b->inc_low) < KEY_EPSILON
&& fabs(a->inc_high - b->inc_high) < KEY_EPSILON);
PG_RETURN_POINTER(result); PG_RETURN_POINTER(result);
} }
@ -489,17 +553,20 @@ gist_tle_same(PG_FUNCTION_ARGS)
* Returns the minimum altitude-band separation in km. * Returns the minimum altitude-band separation in km.
* For overlapping ranges this is 0, making the entry a candidate. * For overlapping ranges this is 0, making the entry a candidate.
* The planner uses this to drive ORDER BY <-> queries. * The planner uses this to drive ORDER BY <-> queries.
*
* Altitude-only: conjunction screening is altitude-dominant.
* Inclination weighting can be added later if needed.
*/ */
Datum Datum
gist_tle_distance(PG_FUNCTION_ARGS) gist_tle_distance(PG_FUNCTION_ARGS)
{ {
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1); pg_tle *query = (pg_tle *) PG_GETARG_POINTER(1);
/* strategy number at arg 2, unused for single-distance class */ /* strategy number at arg 2, unused for single-distance class */
tle_alt_range *key = (tle_alt_range *) DatumGetPointer(entry->key); tle_orbital_key *key = (tle_orbital_key *) DatumGetPointer(entry->key);
tle_alt_range query_range; tle_orbital_key query_key;
tle_to_alt_range(query, &query_range); tle_to_orbital_key(query, &query_key);
PG_RETURN_FLOAT8(range_separation(key, &query_range)); PG_RETURN_FLOAT8(alt_separation(key, &query_key));
} }

75
test/expected/gist_index.out
View File

@ -1,4 +1,4 @@
-- Test GiST index and operators -- Test GiST index and operators (2-D: altitude + inclination)
CREATE EXTENSION IF NOT EXISTS pg_orbit; CREATE EXTENSION IF NOT EXISTS pg_orbit;
NOTICE: extension "pg_orbit" already exists, skipping NOTICE: extension "pg_orbit" already exists, skipping
-- Create test table with mixed orbit types -- Create test table with mixed orbit types
@ -7,46 +7,59 @@ CREATE TABLE test_orbits (
name text, name text,
tle tle tle tle
); );
-- ISS (LEO, ~400km) -- ISS (LEO, ~400km, inclination 51.64 deg)
INSERT INTO test_orbits (name, tle) VALUES ('ISS', INSERT INTO test_orbits (name, tle) VALUES ('ISS',
'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025 '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'); 2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Hubble (LEO, ~540km) -- Hubble (LEO, ~540km, inclination 28.47 deg)
INSERT INTO test_orbits (name, tle) VALUES ('Hubble', INSERT INTO test_orbits (name, tle) VALUES ('Hubble',
'1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992 '1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992
2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008'); 2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008');
-- GPS IIR-M (MEO, ~20200km) -- GPS IIR-M (MEO, ~20200km, inclination 55.44 deg)
INSERT INTO test_orbits (name, tle) VALUES ('GPS-IIR', INSERT INTO test_orbits (name, tle) VALUES ('GPS-IIR',
'1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993 '1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993
2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006'); 2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006');
-- Equatorial-LEO: same altitude as ISS but near-equatorial inclination (5 deg).
-- Same mean motion (15.501 rev/day) and eccentricity as ISS → same altitude band.
-- Different inclination → should NOT overlap with ISS under 2-D key.
INSERT INTO test_orbits (name, tle) VALUES ('Equatorial-LEO',
'1 99901U 24999A 24001.50000000 .00016717 00000-0 10270-3 0 9990
2 99901 5.0000 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Create GiST index -- Create GiST index
CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle); CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle);
-- Overlap: ISS && Hubble = false (different altitude bands: ~400km vs ~540km) -- 2-D overlap: ISS && Equatorial-LEO should be false (same altitude, different inclination)
SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps
FROM test_orbits a, test_orbits b FROM test_orbits a, test_orbits b
WHERE a.id < b.id WHERE a.id < b.id
ORDER BY a.name, b.name; ORDER BY a.name, b.name;
sat_a | sat_b | overlaps sat_a | sat_b | overlaps
--------+---------+---------- ---------+----------------+----------
Hubble | GPS-IIR | f GPS-IIR | Equatorial-LEO | f
ISS | GPS-IIR | f Hubble | Equatorial-LEO | f
ISS | Hubble | f Hubble | GPS-IIR | f
(3 rows) ISS | Equatorial-LEO | f
ISS | GPS-IIR | f
ISS | Hubble | f
(6 rows)
-- Altitude distance between different orbit regimes -- Altitude distance: ISS <-> Equatorial-LEO should be ~0 (same altitude shell)
SELECT a.name AS sat_a, b.name AS sat_b, SELECT a.name AS sat_a, b.name AS sat_b,
round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km
FROM test_orbits a, test_orbits b FROM test_orbits a, test_orbits b
WHERE a.id < b.id WHERE a.id < b.id
ORDER BY a.name, b.name; ORDER BY a.name, b.name;
sat_a | sat_b | alt_dist_km sat_a | sat_b | alt_dist_km
--------+---------+------------- ---------+----------------+-------------
Hubble | GPS-IIR | 19332 GPS-IIR | Equatorial-LEO | 19451
ISS | GPS-IIR | 19451 Hubble | Equatorial-LEO | 115
ISS | Hubble | 115 Hubble | GPS-IIR | 19332
(3 rows) ISS | Equatorial-LEO | 0
ISS | GPS-IIR | 19451
ISS | Hubble | 115
(6 rows)
-- GiST index scan: find all LEO sats (overlap with ISS) -- GiST index scan: find all sats overlapping ISS (altitude AND inclination)
-- Should return only ISS itself (Equatorial-LEO is at same altitude but wrong inclination)
SET enable_seqscan = off; SET enable_seqscan = off;
SELECT name FROM test_orbits SELECT name FROM test_orbits
WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS') WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS')
@ -63,21 +76,23 @@ SELECT name, round((tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS'))::n
FROM test_orbits FROM test_orbits
WHERE name != 'ISS' WHERE name != 'ISS'
ORDER BY tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS'); ORDER BY tle <-> (SELECT tle FROM test_orbits WHERE name = 'ISS');
name | dist name | dist
---------+------- ----------------+-------
Hubble | 115 Equatorial-LEO | 0
GPS-IIR | 19451 Hubble | 115
(2 rows) GPS-IIR | 19451
(3 rows)
RESET enable_seqscan; RESET enable_seqscan;
-- Self-overlap is always true -- Self-overlap is always true
SELECT name, tle && tle AS self_overlap FROM test_orbits ORDER BY name; SELECT name, tle && tle AS self_overlap FROM test_orbits ORDER BY name;
name | self_overlap name | self_overlap
---------+-------------- ----------------+--------------
GPS-IIR | t Equatorial-LEO | t
Hubble | t GPS-IIR | t
ISS | t Hubble | t
(3 rows) ISS | t
(4 rows)
-- Cleanup -- Cleanup
DROP TABLE test_orbits; DROP TABLE test_orbits;

22
test/sql/gist_index.sql
View File

@ -1,4 +1,4 @@
-- Test GiST index and operators -- Test GiST index and operators (2-D: altitude + inclination)
CREATE EXTENSION IF NOT EXISTS pg_orbit; CREATE EXTENSION IF NOT EXISTS pg_orbit;
-- Create test table with mixed orbit types -- Create test table with mixed orbit types
@ -8,38 +8,46 @@ CREATE TABLE test_orbits (
tle tle tle tle
); );
-- ISS (LEO, ~400km) -- ISS (LEO, ~400km, inclination 51.64 deg)
INSERT INTO test_orbits (name, tle) VALUES ('ISS', INSERT INTO test_orbits (name, tle) VALUES ('ISS',
'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025 '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'); 2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Hubble (LEO, ~540km) -- Hubble (LEO, ~540km, inclination 28.47 deg)
INSERT INTO test_orbits (name, tle) VALUES ('Hubble', INSERT INTO test_orbits (name, tle) VALUES ('Hubble',
'1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992 '1 20580U 90037B 24001.50000000 .00000790 00000+0 39573-4 0 9992
2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008'); 2 20580 28.4705 61.4398 0002797 317.3115 42.7577 15.09395228 00008');
-- GPS IIR-M (MEO, ~20200km) -- GPS IIR-M (MEO, ~20200km, inclination 55.44 deg)
INSERT INTO test_orbits (name, tle) VALUES ('GPS-IIR', INSERT INTO test_orbits (name, tle) VALUES ('GPS-IIR',
'1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993 '1 28874U 05038A 24001.50000000 .00000012 00000+0 00000+0 0 9993
2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006'); 2 28874 55.4408 300.3467 0117034 51.6543 309.5420 2.00557079 00006');
-- Equatorial-LEO: same altitude as ISS but near-equatorial inclination (5 deg).
-- Same mean motion (15.501 rev/day) and eccentricity as ISS → same altitude band.
-- Different inclination → should NOT overlap with ISS under 2-D key.
INSERT INTO test_orbits (name, tle) VALUES ('Equatorial-LEO',
'1 99901U 24999A 24001.50000000 .00016717 00000-0 10270-3 0 9990
2 99901 5.0000 208.9163 0006703 30.1694 61.7520 15.50100486 00001');
-- Create GiST index -- Create GiST index
CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle); CREATE INDEX test_orbits_gist ON test_orbits USING gist (tle);
-- Overlap: ISS && Hubble = false (different altitude bands: ~400km vs ~540km) -- 2-D overlap: ISS && Equatorial-LEO should be false (same altitude, different inclination)
SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps SELECT a.name AS sat_a, b.name AS sat_b, a.tle && b.tle AS overlaps
FROM test_orbits a, test_orbits b FROM test_orbits a, test_orbits b
WHERE a.id < b.id WHERE a.id < b.id
ORDER BY a.name, b.name; ORDER BY a.name, b.name;
-- Altitude distance between different orbit regimes -- Altitude distance: ISS <-> Equatorial-LEO should be ~0 (same altitude shell)
SELECT a.name AS sat_a, b.name AS sat_b, SELECT a.name AS sat_a, b.name AS sat_b,
round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km round((a.tle <-> b.tle)::numeric, 0) AS alt_dist_km
FROM test_orbits a, test_orbits b FROM test_orbits a, test_orbits b
WHERE a.id < b.id WHERE a.id < b.id
ORDER BY a.name, b.name; ORDER BY a.name, b.name;
-- GiST index scan: find all LEO sats (overlap with ISS) -- GiST index scan: find all sats overlapping ISS (altitude AND inclination)
-- Should return only ISS itself (Equatorial-LEO is at same altitude but wrong inclination)
SET enable_seqscan = off; SET enable_seqscan = off;
SELECT name FROM test_orbits SELECT name FROM test_orbits
WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS') WHERE tle && (SELECT tle FROM test_orbits WHERE name = 'ISS')