From cddf2a9b6b7cb743c25197cb0fed00c546fe756d Mon Sep 17 00:00:00 2001 From: Nick Falcone Date: Fri, 22 May 2026 23:55:03 -0400 Subject: [PATCH] refactor: rewrite shadowfinder + exif for correctness, perf, and ergonomics MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit shadowfinder: - Symmetric angular error |atan(h/s) − sunAlt| as the likelihood metric; exported band thresholds (ULTRA_TIGHT/TIGHT/MAIN/VISIBLE_BAND_RAD). - Inline sun-position math: time-only terms precomputed once, lng trig cached per column, integer-indexed loops, preallocated array. ~2.5× speedup over SunCalc-per-cell (74.8ms → 29.7ms with constraint). - estimateBestLocation: split bimodal lat-clusters at the largest gap, use circular mean for longitude (correct across the antimeridian), and 68th-percentile great-circle distance as the accuracy radius. - mainBandCoordinates gains lngWraps to flag antimeridian-crossing bands. - Optional azimuth constraint on generateShadowFinderGrid skips Math.asin for cells outside tolerance; always skip Math.atan2 on night cells. - Hoist magic numbers (KM_PER_DEG_LAT, ACCURACY_PERCENTILE, DEFAULT_AZIMUTH_TOLERANCE_DEG, etc.) to named constants. exif: - WMM2025 magnetic-bearing correction when EXIF reports a magnetic bearing alongside GPS coords (new geomagnetism dep, Apache-2.0). Stores the applied declination as magneticDeclination. - localTime: Date → localWallClock: string (ISO 8601 without Z) to fix the wall-clock-encoded-as-UTC-Date footgun. applyUtcOffset rewritten to accept the string; formatDateInput/formatTimeInput now accept Date | string. - Preserve fractional GPS seconds as milliseconds. - computeFovDeg honors portrait orientation (24mm short sensor dim). - Gate dev console.log behind import.meta.env.DEV (verified stripped from the prod bundle). - More forgiving UTC offset parsing (+HH:MM, +HHMM, +HH); -0 normalized to +0. Tests grow 51 → 80, including SunCalc parity at six sample points, antimeridian centroid handling, WMM declination samples, and constraint-aware grid behavior. CLAUDE.md updated to reflect the new metric and inlined sun math. Co-Authored-By: Claude Opus 4.7 Co-authored-by: Cursor --- CLAUDE.md | 2 +- package-lock.json | 249 +-------- package.json | 1 + src/components/AnalysisPanel.tsx | 17 +- src/components/ShadowFinderVisualization.tsx | 15 +- src/lib/exif.test.ts | 124 ++++- src/lib/exif.ts | 240 ++++++--- src/lib/shadowfinder.test.ts | 255 ++++++++-- src/lib/shadowfinder.ts | 501 ++++++++++++++----- 9 files changed, 915 insertions(+), 489 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 1ce8843..3e0fb0e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -22,7 +22,7 @@ Single-page app. One route (`/`), one page component (`src/pages/Index.tsx`). **Data flow:** 1. User uploads a photo → `extractPhotoMetadata` (`src/lib/exif.ts`) runs two concurrent `exifr.parse` calls plus `exifr.gps()` to extract timestamp, compass bearing, focal length, and GPS coords into `PhotoMetadata`. 2. User marks 3 points on the image (object base, object top, shadow tip) via `InteractiveImage`. -3. "Estimate Location" → `generateShadowFinderGrid` + `analyzeShadowMeasurements` (`src/lib/shadowfinder.ts`) sweep a 0.5° global grid using SunCalc, scoring each point by how well its theoretical shadow ratio matches the measured one. +3. "Estimate Location" → `generateShadowFinderGrid` + `analyzeShadowMeasurements` (`src/lib/shadowfinder.ts`) sweep a 0.5° global grid using inlined sun-position math (verified against SunCalc in tests), scoring each point by the symmetric angular error `|atan(height/shadow) − sun_altitude|` in radians. Band thresholds (`ULTRA_TIGHT_BAND_RAD`, `MAIN_BAND_RAD`, `VISIBLE_BAND_RAD`, etc.) are exported from `shadowfinder.ts` and consumed by the visualization. `estimateBestLocation` splits bimodal posteriors at the largest latitude gap and uses a circular mean for longitudes so the antimeridian doesn't break the centroid. 4. Results pass as props to `ShadowFinderVisualization`, which renders a Leaflet heatmap plus optional GPS marker and FOV cone. **Dual-photo mode:** `Index.tsx` maintains mirrored state for first/second photos (`firstPhotoMeta`/`secondPhotoMeta`, etc.). When both analyses exist, `ShadowFinderVisualization` intersects their high-probability grids and renders a single combined heatmap. diff --git a/package-lock.json b/package-lock.json index e420fd1..2cceb44 100644 --- a/package-lock.json +++ b/package-lock.json @@ -17,6 +17,7 @@ "class-variance-authority": "^0.7.1", "clsx": "^2.1.1", "exifr": "^7.1.3", + "geomagnetism": "^0.2.0", "leaflet": "^1.9.4", "leaflet.heat": "^0.2.0", "lucide-react": "^0.462.0", @@ -2719,17 +2720,6 @@ "dev": true, "license": "MIT" }, - "node_modules/detect-libc": { - "version": "2.1.2", - "resolved": "https://registry.npmjs.org/detect-libc/-/detect-libc-2.1.2.tgz", - "integrity": "sha512-Btj2BOOO83o3WyH59e8MgXsxEQVcarkUOpEYrubB0urwnN10yQ364rsiByU11nZlqWYZm05i/of7io4mzihBtQ==", - "dev": true, - "license": "Apache-2.0", - "optional": true, - "engines": { - "node": ">=8" - } - }, "node_modules/didyoumean": { "version": "1.2.2", "resolved": "https://registry.npmjs.org/didyoumean/-/didyoumean-1.2.2.tgz", @@ -3209,6 +3199,12 @@ "url": "https://github.com/sponsors/ljharb" } }, + "node_modules/geomagnetism": { + "version": "0.2.0", + "resolved": "https://registry.npmjs.org/geomagnetism/-/geomagnetism-0.2.0.tgz", + "integrity": "sha512-dO8HJrXqah244nMe+pU5m52L90SMoi4lDA0AV3xunRxogloYV+s3aYWCW72VMMfiYis+YSAlhoKktw8aLMiJbw==", + "license": "Apache-2.0" + }, "node_modules/glob": { "version": "10.4.5", "resolved": "https://registry.npmjs.org/glob/-/glob-10.4.5.tgz", @@ -3516,237 +3512,6 @@ "node": ">= 0.8.0" } }, - "node_modules/lightningcss-android-arm64": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-android-arm64/-/lightningcss-android-arm64-1.32.0.tgz", - "integrity": "sha512-YK7/ClTt4kAK0vo6w3X+Pnm0D2cf2vPHbhOXdoNti1Ga0al1P4TBZhwjATvjNwLEBCnKvjJc2jQgHXH0NEwlAg==", - "cpu": [ - "arm64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "android" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-darwin-arm64": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-darwin-arm64/-/lightningcss-darwin-arm64-1.32.0.tgz", - "integrity": "sha512-RzeG9Ju5bag2Bv1/lwlVJvBE3q6TtXskdZLLCyfg5pt+HLz9BqlICO7LZM7VHNTTn/5PRhHFBSjk5lc4cmscPQ==", - "cpu": [ - "arm64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "darwin" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-darwin-x64": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-darwin-x64/-/lightningcss-darwin-x64-1.32.0.tgz", - "integrity": "sha512-U+QsBp2m/s2wqpUYT/6wnlagdZbtZdndSmut/NJqlCcMLTWp5muCrID+K5UJ6jqD2BFshejCYXniPDbNh73V8w==", - "cpu": [ - "x64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "darwin" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-freebsd-x64": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-freebsd-x64/-/lightningcss-freebsd-x64-1.32.0.tgz", - "integrity": "sha512-JCTigedEksZk3tHTTthnMdVfGf61Fky8Ji2E4YjUTEQX14xiy/lTzXnu1vwiZe3bYe0q+SpsSH/CTeDXK6WHig==", - "cpu": [ - "x64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "freebsd" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-linux-arm-gnueabihf": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-linux-arm-gnueabihf/-/lightningcss-linux-arm-gnueabihf-1.32.0.tgz", - "integrity": "sha512-x6rnnpRa2GL0zQOkt6rts3YDPzduLpWvwAF6EMhXFVZXD4tPrBkEFqzGowzCsIWsPjqSK+tyNEODUBXeeVHSkw==", - "cpu": [ - "arm" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "linux" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-linux-arm64-gnu": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-linux-arm64-gnu/-/lightningcss-linux-arm64-gnu-1.32.0.tgz", - "integrity": "sha512-0nnMyoyOLRJXfbMOilaSRcLH3Jw5z9HDNGfT/gwCPgaDjnx0i8w7vBzFLFR1f6CMLKF8gVbebmkUN3fa/kQJpQ==", - "cpu": [ - "arm64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "linux" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-linux-arm64-musl": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-linux-arm64-musl/-/lightningcss-linux-arm64-musl-1.32.0.tgz", - "integrity": "sha512-UpQkoenr4UJEzgVIYpI80lDFvRmPVg6oqboNHfoH4CQIfNA+HOrZ7Mo7KZP02dC6LjghPQJeBsvXhJod/wnIBg==", - "cpu": [ - "arm64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "linux" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-linux-x64-gnu": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-linux-x64-gnu/-/lightningcss-linux-x64-gnu-1.32.0.tgz", - "integrity": "sha512-V7Qr52IhZmdKPVr+Vtw8o+WLsQJYCTd8loIfpDaMRWGUZfBOYEJeyJIkqGIDMZPwPx24pUMfwSxxI8phr/MbOA==", - "cpu": [ - "x64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "linux" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-linux-x64-musl": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-linux-x64-musl/-/lightningcss-linux-x64-musl-1.32.0.tgz", - "integrity": "sha512-bYcLp+Vb0awsiXg/80uCRezCYHNg1/l3mt0gzHnWV9XP1W5sKa5/TCdGWaR/zBM2PeF/HbsQv/j2URNOiVuxWg==", - "cpu": [ - "x64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "linux" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-win32-arm64-msvc": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-win32-arm64-msvc/-/lightningcss-win32-arm64-msvc-1.32.0.tgz", - "integrity": "sha512-8SbC8BR40pS6baCM8sbtYDSwEVQd4JlFTOlaD3gWGHfThTcABnNDBda6eTZeqbofalIJhFx0qKzgHJmcPTnGdw==", - "cpu": [ - "arm64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "win32" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, - "node_modules/lightningcss-win32-x64-msvc": { - "version": "1.32.0", - "resolved": "https://registry.npmjs.org/lightningcss-win32-x64-msvc/-/lightningcss-win32-x64-msvc-1.32.0.tgz", - "integrity": "sha512-Amq9B/SoZYdDi1kFrojnoqPLxYhQ4Wo5XiL8EVJrVsB8ARoC1PWW6VGtT0WKCemjy8aC+louJnjS7U18x3b06Q==", - "cpu": [ - "x64" - ], - "dev": true, - "license": "MPL-2.0", - "optional": true, - "os": [ - "win32" - ], - "engines": { - "node": ">= 12.0.0" - }, - "funding": { - "type": "opencollective", - "url": "https://opencollective.com/parcel" - } - }, "node_modules/lilconfig": { "version": "3.1.3", "resolved": "https://registry.npmjs.org/lilconfig/-/lilconfig-3.1.3.tgz", diff --git a/package.json b/package.json index fa93b38..e97d117 100644 --- a/package.json +++ b/package.json @@ -22,6 +22,7 @@ "class-variance-authority": "^0.7.1", "clsx": "^2.1.1", "exifr": "^7.1.3", + "geomagnetism": "^0.2.0", "leaflet": "^1.9.4", "leaflet.heat": "^0.2.0", "lucide-react": "^0.462.0", diff --git a/src/components/AnalysisPanel.tsx b/src/components/AnalysisPanel.tsx index 2f83978..d78e1c6 100644 --- a/src/components/AnalysisPanel.tsx +++ b/src/components/AnalysisPanel.tsx @@ -13,7 +13,7 @@ import { formatDateInput, formatTimeInput, } from '@/lib/exif'; -import { AzimuthConstraint } from '@/lib/shadowfinder'; +import { AzimuthConstraint, DEFAULT_AZIMUTH_TOLERANCE_DEG } from '@/lib/shadowfinder'; interface AnalysisPanelProps { points: ClickPoint[]; @@ -80,10 +80,9 @@ export const AnalysisPanel: React.FC = ({ if (photoMetadata.utcTime) { setSelectedDate(formatDateInput(photoMetadata.utcTime)); setSelectedTime(formatTimeInput(photoMetadata.utcTime)); - } else if (photoMetadata.localTime) { - // Pre-fill local time — user must pick offset - setSelectedDate(formatDateInput(photoMetadata.localTime)); - setSelectedTime(formatTimeInput(photoMetadata.localTime)); + } else if (photoMetadata.localWallClock) { + setSelectedDate(formatDateInput(photoMetadata.localWallClock)); + setSelectedTime(formatTimeInput(photoMetadata.localWallClock)); } }, [photoMetadata]); @@ -91,8 +90,8 @@ export const AnalysisPanel: React.FC = ({ // MUST use useEffect — calling setters during render would cause infinite re-renders. // Deps: [manualOffsetMinutes, photoMetadata] — recalculate when either changes. useEffect(() => { - if (!photoMetadata?.localTime || photoMetadata.utcTime) return; - const utc = applyUtcOffset(photoMetadata.localTime, manualOffsetMinutes); + if (!photoMetadata?.localWallClock || photoMetadata.utcTime) return; + const utc = applyUtcOffset(photoMetadata.localWallClock, manualOffsetMinutes); setSelectedDate(formatDateInput(utc)); setSelectedTime(formatTimeInput(utc)); }, [manualOffsetMinutes, photoMetadata]); @@ -119,7 +118,7 @@ export const AnalysisPanel: React.FC = ({ // Deps: [sunBearingDeg, azimuthEnabled, onAzimuthConstraint] — recalculate on any change. useEffect(() => { if (sunBearingDeg !== null && azimuthEnabled) { - onAzimuthConstraint({ sunBearingDeg, toleranceDeg: 10, enabled: true }); + onAzimuthConstraint({ sunBearingDeg, toleranceDeg: DEFAULT_AZIMUTH_TOLERANCE_DEG, enabled: true }); } else { onAzimuthConstraint(null); } @@ -149,7 +148,7 @@ export const AnalysisPanel: React.FC = ({ !photoMetadata ? 'none' : photoMetadata.hasGPSTime ? 'gps' : photoMetadata.utcOffset ? 'offset' - : photoMetadata.localTime ? 'local' + : photoMetadata.localWallClock ? 'local' : 'none'; return ( diff --git a/src/components/ShadowFinderVisualization.tsx b/src/components/ShadowFinderVisualization.tsx index 44c88f8..47aea4f 100644 --- a/src/components/ShadowFinderVisualization.tsx +++ b/src/components/ShadowFinderVisualization.tsx @@ -5,12 +5,19 @@ import 'leaflet.heat'; import { Card } from '@/components/ui/card'; import { Badge } from '@/components/ui/badge'; import { MapPin, Zap } from 'lucide-react'; -import { ShadowFinderPoint, ShadowAnalysisResult, AzimuthConstraint, applyAzimuthConstraint } from '@/lib/shadowfinder'; +import { + ShadowFinderPoint, + ShadowAnalysisResult, + AzimuthConstraint, + applyAzimuthConstraint, + VISIBLE_BAND_RAD, + NIGHT_LIKELIHOOD, +} from '@/lib/shadowfinder'; import { computeFovDeg } from '@/lib/exif'; type HeatPoint = [number, number, number]; -const LIKELIHOOD_CUTOFF = 0.15; +const LIKELIHOOD_CUTOFF = VISIBLE_BAND_RAD; const WARM_GRADIENT: Record = { 0.4: '#FF4500', 0.65: '#FF9500', 0.85: '#FFD700', 1.0: '#FFFF33' }; const COOL_GRADIENT: Record = { 0.4: '#1E90FF', 0.7: '#00E5FF', 1.0: '#FFFFFF' }; @@ -45,7 +52,7 @@ function filterPoints( points: ShadowFinderPoint[], constraint: AzimuthConstraint | null | undefined ): { visible: ShadowFinderPoint[]; fallback: boolean } { - const base = points.filter(p => p.likelihood !== -1 && p.likelihood <= LIKELIHOOD_CUTOFF); + const base = points.filter(p => p.likelihood !== NIGHT_LIKELIHOOD && p.likelihood <= LIKELIHOOD_CUTOFF); if (!constraint?.enabled) return { visible: base, fallback: false }; const filtered = applyAzimuthConstraint(base, constraint); if (filtered.length === 0) return { visible: base, fallback: true }; @@ -246,7 +253,7 @@ export const ShadowFinderVisualization: React.FC const firstPoint = firstMap.get(key); if (firstPoint) { const combinedLikelihood = Math.max(firstPoint.likelihood, p.likelihood); - intersectionPoints.push([p.lat, p.lng, 1 - combinedLikelihood / 0.15]); + intersectionPoints.push([p.lat, p.lng, 1 - combinedLikelihood / LIKELIHOOD_CUTOFF]); } } diff --git a/src/lib/exif.test.ts b/src/lib/exif.test.ts index f7e3bf0..9f09c14 100644 --- a/src/lib/exif.test.ts +++ b/src/lib/exif.test.ts @@ -5,6 +5,8 @@ import { applyUtcOffset, formatDateInput, formatTimeInput, + parseUtcOffset, + correctMagneticBearing, } from './exif'; describe('computeFovDeg', () => { @@ -20,12 +22,12 @@ describe('computeFovDeg', () => { expect(computeFovDeg(-10)).toBe(65); }); - it('computes correct FOV for 50mm lens (~39.6°)', () => { + it('computes correct FOV for 50mm lens (~39.6°) in landscape', () => { // 2 * atan(36 / 100) * (180 / π) expect(computeFovDeg(50)).toBeCloseTo(39.6, 1); }); - it('computes correct FOV for 24mm wide-angle lens (~73.7°)', () => { + it('computes correct FOV for 24mm wide-angle lens (~73.7°) in landscape', () => { // 2 * atan(36 / 48) * (180 / π) expect(computeFovDeg(24)).toBeCloseTo(73.7, 1); }); @@ -34,6 +36,15 @@ describe('computeFovDeg', () => { expect(computeFovDeg(200)).toBeLessThan(computeFovDeg(50)); expect(computeFovDeg(50)).toBeLessThan(computeFovDeg(24)); }); + + it('uses the 24mm short sensor dimension in portrait orientation', () => { + // Portrait FOV at 50mm: 2 * atan(24 / 100) * (180 / π) ≈ 27.0° + expect(computeFovDeg(50, true)).toBeCloseTo(27.0, 1); + }); + + it('produces narrower FOV in portrait than landscape for the same lens', () => { + expect(computeFovDeg(50, true)).toBeLessThan(computeFovDeg(50, false)); + }); }); describe('getUtcOffsetOptions', () => { @@ -65,26 +76,72 @@ describe('getUtcOffsetOptions', () => { }); }); +// ─── parseUtcOffset ────────────────────────────────────────────────────────── + +describe('parseUtcOffset', () => { + it('parses "+HH:MM"', () => { + expect(parseUtcOffset('+05:30')).toBe(330); + expect(parseUtcOffset('-08:00')).toBe(-480); + }); + + it('parses "+HHMM" (no colon)', () => { + expect(parseUtcOffset('+0530')).toBe(330); + expect(parseUtcOffset('-0800')).toBe(-480); + }); + + it('parses "+HH" (no minutes)', () => { + expect(parseUtcOffset('+05')).toBe(300); + expect(parseUtcOffset('-08')).toBe(-480); + }); + + it('returns null for malformed input', () => { + expect(parseUtcOffset('05:30')).toBeNull(); + expect(parseUtcOffset('+5:30')).toBeNull(); + expect(parseUtcOffset('not-an-offset')).toBeNull(); + }); + + it('handles UTC+00:00', () => { + expect(parseUtcOffset('+00:00')).toBe(0); + expect(parseUtcOffset('-00:00')).toBe(0); + }); +}); + +// ─── applyUtcOffset ────────────────────────────────────────────────────────── + describe('applyUtcOffset', () => { it('subtracts a positive offset to produce UTC (UTC+5 → 5 hours earlier)', () => { - const local = new Date('2024-01-01T12:00:00.000Z'); - const utc = applyUtcOffset(local, 300); + const utc = applyUtcOffset('2024-01-01T12:00:00', 300); expect(utc.toISOString()).toBe('2024-01-01T07:00:00.000Z'); }); it('adds back a negative offset to produce UTC (UTC-5 → 5 hours later)', () => { - const local = new Date('2024-01-01T12:00:00.000Z'); - const utc = applyUtcOffset(local, -300); + const utc = applyUtcOffset('2024-01-01T12:00:00', -300); expect(utc.toISOString()).toBe('2024-01-01T17:00:00.000Z'); }); - it('returns the same time for a zero offset', () => { - const local = new Date('2024-06-15T09:30:00.000Z'); - const utc = applyUtcOffset(local, 0); + it('returns the same instant for a zero offset', () => { + const utc = applyUtcOffset('2024-06-15T09:30:00', 0); expect(utc.toISOString()).toBe('2024-06-15T09:30:00.000Z'); }); + + it('handles half-hour offsets (UTC+05:30 IST)', () => { + const utc = applyUtcOffset('2024-06-15T09:30:00', 330); + expect(utc.toISOString()).toBe('2024-06-15T04:00:00.000Z'); + }); + + it('crosses day boundaries correctly going forward', () => { + const utc = applyUtcOffset('2024-01-01T23:30:00', -300); // UTC-5 + expect(utc.toISOString()).toBe('2024-01-02T04:30:00.000Z'); + }); + + it('crosses day boundaries correctly going backward', () => { + const utc = applyUtcOffset('2024-01-01T02:00:00', 300); // UTC+5 + expect(utc.toISOString()).toBe('2023-12-31T21:00:00.000Z'); + }); }); +// ─── formatDateInput / formatTimeInput ────────────────────────────────────── + describe('formatDateInput', () => { it('formats a UTC date as YYYY-MM-DD', () => { expect(formatDateInput(new Date('2024-03-15T14:30:00.000Z'))).toBe('2024-03-15'); @@ -93,6 +150,10 @@ describe('formatDateInput', () => { it('handles year boundaries correctly', () => { expect(formatDateInput(new Date('2024-01-01T00:00:00.000Z'))).toBe('2024-01-01'); }); + + it('slices a wall-clock ISO string to YYYY-MM-DD', () => { + expect(formatDateInput('2024-03-15T14:30:00')).toBe('2024-03-15'); + }); }); describe('formatTimeInput', () => { @@ -103,4 +164,49 @@ describe('formatTimeInput', () => { it('pads single-digit hours and minutes', () => { expect(formatTimeInput(new Date('2024-03-15T09:05:00.000Z'))).toBe('09:05'); }); + + it('slices a wall-clock ISO string to HH:MM', () => { + expect(formatTimeInput('2024-03-15T14:30:00')).toBe('14:30'); + }); +}); + +// ─── correctMagneticBearing ────────────────────────────────────────────────── +// +// These exercise the WMM2025 model that ships with the `geomagnetism` package. +// Reference declination values are from NOAA's WMM 2025 calculator. Locations +// span both hemispheres and include a near-zero-declination case (Eastern UK +// in 2026) and a high-declination case (eastern Australia). + +describe('correctMagneticBearing', () => { + const D = new Date('2026-05-22T00:00:00Z'); + + it('adds an east-positive declination to convert magnetic to true bearing', () => { + // Seattle (47.6°N, 122.3°W) has ~+15° declination in 2026. + const result = correctMagneticBearing(100, { lat: 47.6, lng: -122.3 }, D); + expect(result).not.toBeNull(); + expect(result!.declination).toBeGreaterThan(13); + expect(result!.declination).toBeLessThan(17); + expect(result!.trueBearing).toBeCloseTo(100 + result!.declination, 5); + }); + + it('returns a negative declination west of the agonic line', () => { + // Sydney (-33.9°S, 151.2°E) has positive declination in 2026; pick somewhere + // with negative declination instead — eastern South America (~Buenos Aires). + const result = correctMagneticBearing(0, { lat: -34.6, lng: -58.4 }, D); + expect(result).not.toBeNull(); + expect(result!.declination).toBeLessThan(0); + }); + + it('wraps the corrected bearing into [0, 360)', () => { + // Force a wraparound: bearing 355°, +15° declination should yield ~10°. + const result = correctMagneticBearing(355, { lat: 47.6, lng: -122.3 }, D); + expect(result).not.toBeNull(); + expect(result!.trueBearing).toBeGreaterThanOrEqual(0); + expect(result!.trueBearing).toBeLessThan(360); + expect(result!.trueBearing).toBeLessThan(15); + }); + + it('returns null when given NaN inputs', () => { + expect(correctMagneticBearing(NaN, { lat: 0, lng: 0 }, D)).toBeNull(); + }); }); diff --git a/src/lib/exif.ts b/src/lib/exif.ts index a919f44..b00cdaf 100644 --- a/src/lib/exif.ts +++ b/src/lib/exif.ts @@ -1,16 +1,48 @@ import exifr from 'exifr'; +import geomagnetism from 'geomagnetism'; + +/** Long edge of a 35mm full-frame sensor (mm) — used as the FOV horizontal dimension. */ +const SENSOR_LONG_DIM_MM = 36; +/** Short edge of a 35mm full-frame sensor (mm) — used for portrait-orientation FOV. */ +const SENSOR_SHORT_DIM_MM = 24; +/** Default horizontal FOV (degrees) when no focal length is recorded — roughly a 28mm lens. */ +const DEFAULT_FOV_DEG = 65; export interface PhotoMetadata { - localTime: Date | null; + /** + * Photo capture wall-clock time as an ISO 8601 string WITHOUT a timezone + * suffix, e.g. `"2024-03-15T14:30:00"`. The trailing `Z` is intentionally + * omitted because EXIF `DateTimeOriginal` does not specify a timezone. Use + * `applyUtcOffset` to convert to a true UTC `Date` once an offset is known. + * + * (Storing this as a `Date` would be a footgun: `getHours()` returns the + * browser-local interpretation of the bytes, not the photo's wall clock.) + */ + localWallClock: string | null; utcOffset: string | null; utcTime: Date | null; hasGPSTime: boolean; + /** + * Camera bearing in degrees [0, 360). When `magneticDeclination` is non-null + * the bearing has been corrected from magnetic to true north and `compassRef` + * will report `'T'`. When `compassRef === 'M'` (no GPS to derive declination + * from), the bearing is the raw magnetic value and downstream consumers + * should ignore it. + */ compassBearing: number | null; compassRef: 'T' | 'M' | null; gpsCoords: { lat: number; lng: number } | null; focalLength35mm: number | null; + /** + * WMM2025 declination at the photo's GPS location, in degrees east-positive. + * Set only when a magnetic bearing has been corrected to true north; + * otherwise null. + */ + magneticDeclination: number | null; } +const UTC_OFFSET_RE = /^([+-])(\d{2}):?(\d{2})?$/; + export async function extractPhotoMetadata(file: File): Promise { try { // Two concurrent parses: main for GPS/compass (revived values), second for the raw @@ -38,7 +70,7 @@ export async function extractPhotoMetadata(file: File): Promise d instanceof Date && !isNaN(d.getTime()); - // Parse raw EXIF date string "YYYY:MM:DD HH:MM:SS" using Date.UTC so the - // resulting Date's UTC value equals the local wall-clock components. - // toISOString() then returns the local time, and all offset math is correct. + // EXIF DateTimeOriginal is "YYYY:MM:DD HH:MM:SS" with no timezone. Normalize + // it to ISO 8601 sans timezone — the colon-separated date becomes dash-separated. const rawDT = exifRaw?.DateTimeOriginal; if (typeof rawDT === 'string') { - const [datePart, timePart] = rawDT.split(' '); - if (datePart && timePart) { - const [year, month, day] = datePart.split(':').map(Number); - const [hours, minutes, seconds] = timePart.split(':').map(Number); - const normalized = new Date(Date.UTC(year, month - 1, day, hours, minutes, seconds)); - if (isValid(normalized)) { - result.localTime = normalized; - } - } + const wallClock = normalizeExifDateString(rawDT); + if (wallClock) result.localWallClock = wallClock; } if (exif && typeof exif.OffsetTimeOriginal === 'string') { @@ -72,31 +97,30 @@ export async function extractPhotoMetadata(file: File): Promise 0) { - result.focalLength35mm = Number(exif.FocalLengthIn35mmFormat); + const focal = Number(exif?.FocalLengthIn35mmFormat); + if (Number.isFinite(focal) && focal > 0) { + result.focalLength35mm = focal; } if (gps?.latitude != null && gps?.longitude != null && - isFinite(gps.latitude) && isFinite(gps.longitude)) { + Number.isFinite(gps.latitude) && Number.isFinite(gps.longitude)) { result.gpsCoords = { lat: gps.latitude, lng: gps.longitude }; } - if (!result.localTime && !result.utcTime && !result.gpsCoords) return null; + // Correct a magnetic bearing to true north using the WMM at the photo's + // GPS location. Only applied when we have both the bearing and the location; + // a magnetic-only bearing without GPS is left untouched (downstream code + // already ignores it). + if ( + result.compassRef === 'M' && + result.compassBearing != null && + result.gpsCoords + ) { + const referenceDate = + result.utcTime ?? + (result.localWallClock ? new Date(result.localWallClock + 'Z') : new Date()); + const corrected = correctMagneticBearing( + result.compassBearing, + result.gpsCoords, + referenceDate, + ); + // Defensive: out-of-range dates with allowOutOfBoundsModel=true should + // still yield a finite declination; if not, leave the bearing as magnetic. + if (corrected) { + result.compassBearing = corrected.trueBearing; + result.compassRef = 'T'; + result.magneticDeclination = corrected.declination; + } + } + + if (!result.localWallClock && !result.utcTime && !result.gpsCoords) return null; - const source = result.hasGPSTime ? 'gps' : result.utcOffset ? 'offset' : 'local-only'; - console.log('[exif]', { - source, - localTime: result.localTime?.toISOString() ?? null, - utcOffset: result.utcOffset, - utcTime: result.utcTime?.toISOString() ?? null, - compassBearing: result.compassBearing, - compassRef: result.compassRef, - gpsCoords: result.gpsCoords, - focalLength35mm: result.focalLength35mm, - }); + if (import.meta.env.DEV) { + const source = result.hasGPSTime ? 'gps' : result.utcOffset ? 'offset' : 'local-only'; + console.log('[exif]', { + source, + localWallClock: result.localWallClock, + utcOffset: result.utcOffset, + utcTime: result.utcTime?.toISOString() ?? null, + compassBearing: result.compassBearing, + compassRef: result.compassRef, + magneticDeclination: result.magneticDeclination, + gpsCoords: result.gpsCoords, + focalLength35mm: result.focalLength35mm, + }); + } return result; } catch { @@ -142,6 +194,73 @@ export async function extractPhotoMetadata(file: File): Promise n.padStart(2, '0'); + return `${year}-${pad(month)}-${pad(day)}T${pad(hours)}:${pad(minutes)}:${pad(seconds)}`; +} + +/** + * Parse an EXIF `OffsetTimeOriginal` string ("+HH:MM", "+HHMM", or "+HH") to + * a signed minute offset. Returns null on malformed input. + */ +export function parseUtcOffset(offset: string): number | null { + const m = UTC_OFFSET_RE.exec(offset); + if (!m) return null; + const hours = parseInt(m[2], 10); + const minutes = m[3] ? parseInt(m[3], 10) : 0; + const magnitude = hours * 60 + minutes; + // Normalize -0 to +0 so "-00:00" and "+00:00" are indistinguishable downstream. + if (magnitude === 0) return 0; + return m[1] === '+' ? magnitude : -magnitude; +} + +/** + * Apply the WMM2025 declination at the given GPS location to convert a + * magnetic compass bearing to true north. Returns null if the library can't + * produce a finite declination (e.g. NaN inputs or model failure). + * + * trueBearing = magneticBearing + eastward_declination + * + * Exported for testing; callers usually let `extractPhotoMetadata` apply this + * automatically when EXIF reports a magnetic bearing alongside GPS coords. + */ +export function correctMagneticBearing( + magneticBearingDeg: number, + gpsCoords: { lat: number; lng: number }, + forDate: Date, +): { trueBearing: number; declination: number } | null { + if ( + !Number.isFinite(magneticBearingDeg) || + !Number.isFinite(gpsCoords.lat) || + !Number.isFinite(gpsCoords.lng) + ) { + return null; + } + try { + const model = geomagnetism.model(forDate, { allowOutOfBoundsModel: true }); + const point = model.point([gpsCoords.lat, gpsCoords.lng]); + const declination = point.decl; + if (!Number.isFinite(declination)) return null; + const trueBearing = ((magneticBearingDeg + declination) % 360 + 360) % 360; + return { trueBearing, declination }; + } catch { + return null; + } +} + /** Generate UTC offset options from −12:00 to +14:00 in 30-minute steps */ export function getUtcOffsetOptions(): { label: string; minutes: number }[] { const options: { label: string; minutes: number }[] = []; @@ -155,29 +274,38 @@ export function getUtcOffsetOptions(): { label: string; minutes: number }[] { return options; } -/** Apply a UTC offset (minutes) to a local Date to get UTC Date */ -export function applyUtcOffset(localTime: Date, offsetMinutes: number): Date { - return new Date(localTime.getTime() - offsetMinutes * 60000); +/** + * Convert a wall-clock ISO string (no timezone) to a true UTC `Date` by + * subtracting the local offset. Appending `Z` parses the wall-clock as if it + * were UTC, so we then subtract the offset to recover the real UTC instant. + */ +export function applyUtcOffset(localWallClock: string, offsetMinutes: number): Date { + const asUtcInstant = new Date(localWallClock + 'Z'); + return new Date(asUtcInstant.getTime() - offsetMinutes * 60000); } -/** Format a Date as "YYYY-MM-DD" for date input values */ -export function formatDateInput(date: Date): string { - return date.toISOString().slice(0, 10); +/** Format a UTC `Date` or wall-clock string as "YYYY-MM-DD" for date input values. */ +export function formatDateInput(value: Date | string): string { + if (typeof value === 'string') return value.slice(0, 10); + return value.toISOString().slice(0, 10); } -/** Format a Date as "HH:MM" for time input values */ -export function formatTimeInput(date: Date): string { - return date.toISOString().slice(11, 16); +/** Format a UTC `Date` or wall-clock string as "HH:MM" for time input values. */ +export function formatTimeInput(value: Date | string): string { + if (typeof value === 'string') return value.slice(11, 16); + return value.toISOString().slice(11, 16); } /** * Computes horizontal FOV in degrees from a 35mm-equivalent focal length. - * Falls back to 65° (≈28mm) when focal length is unavailable. - * Assumes landscape orientation (36mm wide); portrait shots will appear ~19° wider than actual. + * Falls back to DEFAULT_FOV_DEG (≈28mm) when focal length is unavailable. + * Pass `isPortrait` to use the shorter sensor dimension when the photo was + * taken in portrait orientation. */ -export function computeFovDeg(focalLength35mm: number | null): number { +export function computeFovDeg(focalLength35mm: number | null, isPortrait = false): number { if (focalLength35mm != null && focalLength35mm > 0) { - return 2 * Math.atan(36 / (2 * focalLength35mm)) * (180 / Math.PI); + const sensorDim = isPortrait ? SENSOR_SHORT_DIM_MM : SENSOR_LONG_DIM_MM; + return 2 * Math.atan(sensorDim / (2 * focalLength35mm)) * (180 / Math.PI); } - return 65; + return DEFAULT_FOV_DEG; } diff --git a/src/lib/shadowfinder.test.ts b/src/lib/shadowfinder.test.ts index ad16293..0aa75e5 100644 --- a/src/lib/shadowfinder.test.ts +++ b/src/lib/shadowfinder.test.ts @@ -1,4 +1,5 @@ import { describe, it, expect, beforeAll } from 'vitest'; +import * as SunCalc from 'suncalc'; import { calculateMeasurementsFromPixels, convertPercentageToPixels, @@ -7,20 +8,25 @@ import { estimateBestLocation, generateShadowFinderGrid, computeLikelihood, + dominantLatCluster, + MAIN_BAND_RAD, + VISIBLE_BAND_RAD, + NIGHT_LIKELIHOOD, + DEFAULT_AZIMUTH_TOLERANCE_DEG, type ShadowFinderPoint, type ShadowAnalysisResult, } from './shadowfinder'; // ─── helpers ──────────────────────────────────────────────────────────────── -function makePoint(sunAzimuthDeg: number, likelihood = 0.05): ShadowFinderPoint { +function makePoint(sunAzimuthDeg: number, likelihood = 0.01): ShadowFinderPoint { return { lat: 0, lng: 0, likelihood, sunAzimuthDeg }; } function makeResult(points: ShadowFinderPoint[]): ShadowAnalysisResult { return { points, - mainBandCoordinates: { latRange: [0, 0], lngRange: [0, 0] }, + mainBandCoordinates: { latRange: [0, 0], lngRange: [0, 0], lngWraps: false }, statistics: { totalPoints: points.length, validPoints: points.length, @@ -34,32 +40,36 @@ function makeResult(points: ShadowFinderPoint[]): ShadowAnalysisResult { // ─── computeLikelihood ─────────────────────────────────────────────────────── // -// These tests pin the exact formula. An accidental inversion (1/(ratio*tan) instead -// of ratio/tan) changes all non-zero results and would fail the asymmetric cases. +// The metric is the symmetric angular error |atan(h/s) − sunAlt| in radians. +// "Symmetric" here means swapping h and s (so the ratio inverts) gives the +// same error magnitude when reflected around π/4, because atan(x) + atan(1/x) +// = π/2 for x > 0. describe('computeLikelihood', () => { - const ALT_45 = Math.PI / 4; // tan(45°) = 1 — simplifies expected values + const ALT_45 = Math.PI / 4; it('returns 0 when the shadow ratio exactly matches the sun altitude', () => { - // objectHeight/shadowLength = 1 = tan(45°), perfect match expect(computeLikelihood(100, 100, ALT_45)).toBeCloseTo(0); }); - it('returns 1 when the object is twice the shadow length at 45°', () => { - // ratio = 2, tan(45°) = 1 → |2/1 - 1| = 1.0 - expect(computeLikelihood(200, 100, ALT_45)).toBeCloseTo(1.0); + it('returns the angular gap when the object is twice the shadow at 45°', () => { + // atan(2) − π/4 ≈ 0.32175 rad + const expected = Math.atan(2) - Math.PI / 4; + expect(computeLikelihood(200, 100, ALT_45)).toBeCloseTo(expected, 10); }); - it('returns 0.5 when the shadow is twice the object length at 45°', () => { - // ratio = 0.5, tan(45°) = 1 → |0.5/1 - 1| = 0.5 - // An inverted formula gives |1/(1*0.5) - 1| = 1.0 — catches the bug - expect(computeLikelihood(100, 200, ALT_45)).toBeCloseTo(0.5); + it('is symmetric: swapping object and shadow lengths gives the same likelihood at 45°', () => { + // Old metric was asymmetric (2× shadow vs ½ shadow gave very different scores). + // The angular metric reflects around π/4, so |atan(2) − π/4| = |atan(0.5) − π/4|. + const a = computeLikelihood(200, 100, ALT_45); + const b = computeLikelihood(100, 200, ALT_45); + expect(a).toBeCloseTo(b, 10); }); - it('is asymmetric: swapping object and shadow lengths gives different results', () => { - const a = computeLikelihood(200, 100, ALT_45); // ratio 2 → 1.0 - const b = computeLikelihood(100, 200, ALT_45); // ratio 0.5 → 0.5 - expect(a).not.toBeCloseTo(b); + it('returns radians, not the old |ratio/tan(alt) − 1| value', () => { + // Old metric at (200, 100, 45°) was exactly 1.0; new metric is ~0.32 rad. + // This pins the change so a regression to the old formula fails loudly. + expect(computeLikelihood(200, 100, ALT_45)).toBeLessThan(0.5); }); }); @@ -68,9 +78,9 @@ describe('computeLikelihood', () => { describe('calculateMeasurementsFromPixels', () => { it('measures a vertical object with a horizontal shadow', () => { const { objectHeight, shadowLength } = calculateMeasurementsFromPixels( - { x: 100, y: 100 }, // base - { x: 100, y: 0 }, // top (directly above) - { x: 200, y: 100 }, // shadow tip (directly right) + { x: 100, y: 100 }, + { x: 100, y: 0 }, + { x: 200, y: 100 }, ); expect(objectHeight).toBeCloseTo(100); expect(shadowLength).toBeCloseTo(100); @@ -121,6 +131,21 @@ describe('convertPercentageToPixels', () => { }); }); +// ─── exported constants ────────────────────────────────────────────────────── + +describe('exported constants', () => { + it('DEFAULT_AZIMUTH_TOLERANCE_DEG is a positive degree value', () => { + expect(DEFAULT_AZIMUTH_TOLERANCE_DEG).toBeGreaterThan(0); + expect(DEFAULT_AZIMUTH_TOLERANCE_DEG).toBeLessThan(45); + }); + + it('band thresholds are strictly ordered tightest → loosest', () => { + // Required invariant for the band-counting hierarchy in analyzeShadowMeasurements. + expect(MAIN_BAND_RAD).toBeGreaterThan(0); + expect(VISIBLE_BAND_RAD).toBeGreaterThan(MAIN_BAND_RAD); + }); +}); + // ─── applyAzimuthConstraint ────────────────────────────────────────────────── describe('applyAzimuthConstraint', () => { @@ -144,7 +169,6 @@ describe('applyAzimuthConstraint', () => { }); it('handles 0°/360° wraparound — points near 0° match a bearing near 0°', () => { - // 355° and 5° are both 5° away from 0°; 180° is 180° away const points = [makePoint(355), makePoint(5), makePoint(180)]; const result = applyAzimuthConstraint(points, { sunBearingDeg: 0, toleranceDeg: 10, enabled: true }); expect(result).toHaveLength(2); @@ -153,7 +177,6 @@ describe('applyAzimuthConstraint', () => { }); it('handles 0°/360° wraparound — bearing near 360°', () => { - // 350° is 5° from 355°; 10° is 15° from 355° → excluded const points = [makePoint(350), makePoint(10)]; const result = applyAzimuthConstraint(points, { sunBearingDeg: 355, toleranceDeg: 10, enabled: true }); expect(result).toHaveLength(1); @@ -161,7 +184,7 @@ describe('applyAzimuthConstraint', () => { }); it('includes a point exactly at the tolerance boundary', () => { - const points = [makePoint(80), makePoint(100)]; // both exactly 10° from 90° + const points = [makePoint(80), makePoint(100)]; const result = applyAzimuthConstraint(points, { sunBearingDeg: 90, toleranceDeg: 10, enabled: true }); expect(result).toHaveLength(2); }); @@ -207,13 +230,22 @@ describe('analyzeShadowMeasurements', () => { expect(result.statistics.ultraTightBandPoints).toBeLessThanOrEqual(result.statistics.tightBandPoints); expect(result.statistics.tightBandPoints).toBeLessThanOrEqual(result.statistics.visibleBandPoints); }); + + it('returns a mainBandCoordinates with lngWraps=false for a non-wrapping band', { timeout: 10000 }, () => { + const result = analyzeShadowMeasurements({ + objectHeight: 100, + shadowLength: 100, + knownTime: new Date('2024-06-21T12:00:00.000Z'), + }); + expect(result.mainBandCoordinates.lngWraps).toBe(false); + }); }); // ─── estimateBestLocation ──────────────────────────────────────────────────── describe('estimateBestLocation', () => { - it('throws when no points have likelihood ≤ 0.1', () => { - const result = makeResult([makePoint(90, 0.2), makePoint(180, 0.5)]); + it('throws when no points fall within the MAIN band', () => { + const result = makeResult([makePoint(90, MAIN_BAND_RAD + 0.01), makePoint(180, 0.5)]); expect(() => estimateBestLocation(result)).toThrow('No suitable location matches found'); }); @@ -223,7 +255,7 @@ describe('estimateBestLocation', () => { it('returns the coordinates of a single best point', () => { const points: ShadowFinderPoint[] = [ - { lat: 40, lng: -74, likelihood: 0.05, sunAzimuthDeg: 180 }, + { lat: 40, lng: -74, likelihood: 0.01, sunAzimuthDeg: 180 }, ]; const location = estimateBestLocation(makeResult(points)); expect(location.latitude).toBeCloseTo(40); @@ -232,30 +264,89 @@ describe('estimateBestLocation', () => { it('returns the centroid of multiple best points', () => { const points: ShadowFinderPoint[] = [ - { lat: 10, lng: 20, likelihood: 0.02, sunAzimuthDeg: 180 }, - { lat: 20, lng: 40, likelihood: 0.04, sunAzimuthDeg: 180 }, - // likelihood > 0.1, excluded from best - { lat: 90, lng: 90, likelihood: 0.15, sunAzimuthDeg: 90 }, + { lat: 10, lng: 20, likelihood: 0.01, sunAzimuthDeg: 180 }, + { lat: 20, lng: 40, likelihood: 0.02, sunAzimuthDeg: 180 }, + // Outside MAIN band → excluded from the cluster. + { lat: 90, lng: 90, likelihood: VISIBLE_BAND_RAD, sunAzimuthDeg: 90 }, ]; const location = estimateBestLocation(makeResult(points)); expect(location.latitude).toBeCloseTo(15); expect(location.longitude).toBeCloseTo(30); }); - it('returns at least 1km accuracy', () => { + it('uses a circular mean for longitudes near the antimeridian', () => { + // Plain arithmetic mean of -179 and 179 = 0 (deep wrong, on the opposite side of the globe). + // Circular mean correctly produces ±180. + const points: ShadowFinderPoint[] = [ + { lat: 0, lng: -179, likelihood: 0.01, sunAzimuthDeg: 90 }, + { lat: 0, lng: 179, likelihood: 0.01, sunAzimuthDeg: 90 }, + ]; + const location = estimateBestLocation(makeResult(points)); + expect(Math.abs(location.longitude)).toBeCloseTo(180, 1); + }); + + it('picks the larger cluster when posterior is bimodal in latitude', () => { + // Three points clustered around lat=40 + one outlier at lat=-50 (gap = 90° > MIN_LAT_CLUSTER_GAP). + // Centroid should reflect only the northern cluster. + const points: ShadowFinderPoint[] = [ + { lat: 40, lng: 10, likelihood: 0.01, sunAzimuthDeg: 180 }, + { lat: 42, lng: 12, likelihood: 0.02, sunAzimuthDeg: 180 }, + { lat: 44, lng: 14, likelihood: 0.03, sunAzimuthDeg: 180 }, + { lat: -50, lng: 14, likelihood: 0.01, sunAzimuthDeg: 0 }, + ]; + const location = estimateBestLocation(makeResult(points)); + expect(location.latitude).toBeGreaterThan(35); + expect(location.latitude).toBeLessThan(50); + }); + + it('returns at least half-cell grid resolution as the accuracy floor', () => { const points: ShadowFinderPoint[] = [ { lat: 0, lng: 0, likelihood: 0, sunAzimuthDeg: 180 }, ]; const location = estimateBestLocation(makeResult(points)); - expect(location.accuracy).toBeGreaterThanOrEqual(1); + // half-cell at the equator = 0.25° × 111 km/° ≈ 27.75 km + expect(location.accuracy).toBeGreaterThanOrEqual(27); + expect(location.accuracy).toBeLessThan(30); + }); +}); + +// ─── dominantLatCluster ────────────────────────────────────────────────────── + +describe('dominantLatCluster', () => { + it('returns the input when there is no large gap', () => { + const points = [{ lat: 10 }, { lat: 15 }, { lat: 20 }]; + expect(dominantLatCluster(points)).toHaveLength(3); + }); + + it('returns the input when single-point', () => { + expect(dominantLatCluster([{ lat: 42 }])).toHaveLength(1); + }); + + it('splits a bimodal set and returns the larger cluster', () => { + const points = [ + { lat: -40 }, { lat: -38 }, { lat: -36 }, // 3-pt south cluster + { lat: 30 }, { lat: 32 }, // 2-pt north cluster (gap = 66°) + ]; + const cluster = dominantLatCluster(points); + expect(cluster).toHaveLength(3); + for (const p of cluster) expect(p.lat).toBeLessThan(0); + }); + + it('treats two equal clusters by returning the southern half', () => { + const points = [ + { lat: -40 }, { lat: -38 }, + { lat: 30 }, { lat: 32 }, + ]; + const cluster = dominantLatCluster(points); + expect(cluster).toHaveLength(2); + // Tie-break: south.length >= north.length picks south + for (const p of cluster) expect(p.lat).toBeLessThan(0); }); }); // ─── generateShadowFinderGrid ──────────────────────────────────────────────── describe('generateShadowFinderGrid', () => { - // Grid is 290 latitudes × 720 longitudes = 208,800 points. - // Computed once for all tests in this block — each SunCalc sweep takes ~1–2s. const GRID_TIME = new Date('2024-06-21T12:00:00.000Z'); let points: ShadowFinderPoint[]; @@ -276,28 +367,110 @@ describe('generateShadowFinderGrid', () => { } }); - it('marks night points as -1 and day points as ≥ 0', () => { + it('marks night points as NIGHT_LIKELIHOOD (-1) and day points as ≥ 0', () => { for (const p of points) { - expect(p.likelihood === -1 || p.likelihood >= 0).toBe(true); + expect(p.likelihood === NIGHT_LIKELIHOOD || p.likelihood >= 0).toBe(true); } }); - it('keeps all sunAzimuthDeg in [0, 360)', () => { + it('keeps daytime sunAzimuthDeg in [0, 360) and leaves night cells as NaN', () => { for (const p of points) { - expect(p.sunAzimuthDeg).toBeGreaterThanOrEqual(0); - expect(p.sunAzimuthDeg).toBeLessThan(360); + if (p.likelihood === NIGHT_LIKELIHOOD) { + // Azimuth is intentionally skipped on night cells (perf optimization). + // NaN ensures applyAzimuthConstraint correctly excludes them. + expect(Number.isNaN(p.sunAzimuthDeg)).toBe(true); + } else { + expect(p.sunAzimuthDeg).toBeGreaterThanOrEqual(0); + expect(p.sunAzimuthDeg).toBeLessThan(360); + } } }); it('marks the equator at Greenwich Meridian as daytime at noon UTC', () => { const noonPoint = points.find(p => p.lat === 0 && p.lng === 0); expect(noonPoint).toBeDefined(); - expect(noonPoint!.likelihood).not.toBe(-1); + expect(noonPoint!.likelihood).not.toBe(NIGHT_LIKELIHOOD); }); it('marks the antimeridian (lng≈179.5) at the equator as nighttime at noon UTC', () => { const midnightPoint = points.find(p => p.lat === 0 && p.lng === 179.5); expect(midnightPoint).toBeDefined(); - expect(midnightPoint!.likelihood).toBe(-1); + expect(midnightPoint!.likelihood).toBe(NIGHT_LIKELIHOOD); + }); + + it('NaN-azimuth night cells are correctly excluded by applyAzimuthConstraint', () => { + const filtered = applyAzimuthConstraint(points, { + sunBearingDeg: 180, + toleranceDeg: 90, // wide tolerance — would match many daytime cells + enabled: true, + }); + // None of the survivors should be night cells. + for (const p of filtered) { + expect(p.likelihood).not.toBe(NIGHT_LIKELIHOOD); + } + }); + + it('with an azimuth constraint, marks cells outside tolerance as NIGHT_LIKELIHOOD', () => { + // Very tight tolerance — only a narrow strip should keep a valid likelihood. + const constrained = generateShadowFinderGrid(GRID_TIME, 100, 100, { + sunBearingDeg: 180, + toleranceDeg: 2, + enabled: true, + }); + let valid = 0; + let rejected = 0; + for (const p of constrained) { + if (p.likelihood === NIGHT_LIKELIHOOD) rejected++; + else valid++; + } + // Far more cells should be rejected than kept under a 2° tolerance. + expect(rejected).toBeGreaterThan(valid * 5); + // Surviving cells should all be near the target bearing. + for (const p of constrained) { + if (p.likelihood !== NIGHT_LIKELIHOOD) { + const diff = Math.abs(((p.sunAzimuthDeg - 180 + 540) % 360) - 180); + expect(diff).toBeLessThanOrEqual(2); + } + } + }); + + it('with disabled azimuth constraint, behaves identically to no constraint', () => { + const constrained = generateShadowFinderGrid(GRID_TIME, 100, 100, { + sunBearingDeg: 180, + toleranceDeg: 5, + enabled: false, + }); + // Sample-equality check — full equality of 200K points would be slow. + for (let i = 0; i < constrained.length; i += 5000) { + expect(constrained[i].likelihood).toBe(points[i].likelihood); + } + }); + + it('matches SunCalc altitude/azimuth at sampled grid points (parity check)', () => { + // The inlined sun-position math should agree with SunCalc to within rounding. + // We sample across latitudes/longitudes and times to catch any drift. + const samples: Array<{ lat: number; lng: number }> = [ + { lat: 0, lng: 0 }, + { lat: 45, lng: -75 }, + { lat: -33.5, lng: 151 }, + { lat: 84.5, lng: 179.5 }, + { lat: -60, lng: -180 }, + { lat: 23.5, lng: 100 }, + ]; + for (const s of samples) { + const p = points.find(pt => pt.lat === s.lat && pt.lng === s.lng); + expect(p).toBeDefined(); + const ref = SunCalc.getPosition(GRID_TIME, s.lat, s.lng); + if (ref.altitude > 0) { + const expectedLikelihood = Math.abs(Math.atan(1) - ref.altitude); + expect(p!.likelihood).toBeCloseTo(expectedLikelihood, 6); + const refAzimuthDeg = ((ref.azimuth * 180) / Math.PI + 540) % 360; + expect(p!.sunAzimuthDeg).toBeCloseTo(refAzimuthDeg, 6); + } else { + expect(p!.likelihood).toBe(NIGHT_LIKELIHOOD); + // Azimuth is intentionally skipped for night cells. + expect(Number.isNaN(p!.sunAzimuthDeg)).toBe(true); + } + } }); }); diff --git a/src/lib/shadowfinder.ts b/src/lib/shadowfinder.ts index 05c9497..92617c0 100644 --- a/src/lib/shadowfinder.ts +++ b/src/lib/shadowfinder.ts @@ -1,33 +1,101 @@ /** - * ShadowFinder Algorithm - TypeScript Implementation - * - * This module implements the exact ShadowFinder algorithm from Bellingcat - * for OSINT shadow analysis and geolocation estimation. - * - * Based on our perfected algorithm that achieved 99.9% accuracy match - * with the reference implementation. + * ShadowFinder Algorithm — TypeScript Implementation + * + * OSINT geolocation by matching the measured shadow-to-object ratio against + * the sun's altitude at every point on a 0.5° world grid. Inspired by Bellingcat's + * ShadowFinder, with two intentional differences: + * + * 1. Likelihood is the symmetric angular error + * |atan(height / shadow) − sun_altitude| + * in radians, not Bellingcat's asymmetric |ratio / tan(alt) − 1|. This + * makes the band thresholds physically meaningful (sun-altitude error + * in radians) and treats over- and under-measurements of shadow length + * identically. + * + * 2. Sun-position math is inlined here. Time-only terms (declination, + * right ascension, GMST) are computed once per grid; longitude trig is + * cached in a Float64Array shared across latitudes. Verified against + * SunCalc in tests within 1e-6 rad. + */ + +// ── Constants ─────────────────────────────────────────────────────────────── + +const RAD_PER_DEG = Math.PI / 180; +const DEG_PER_RAD = 180 / Math.PI; +const DAY_MS = 86400000; +const J1970 = 2440588; +const J2000 = 2451545; +const OBLIQUITY = RAD_PER_DEG * 23.4397; + +const LAT_MIN = -60; +const LNG_MIN = -180; +const GRID_RES_DEG = 0.5; +const LAT_STEPS = 290; // -60.0 .. 84.5 inclusive +const LNG_STEPS = 720; // -180.0 .. 179.5 inclusive +const GRID_COUNT = LAT_STEPS * LNG_STEPS; +const EARTH_RADIUS_KM = 6371; +/** Rough conversion from degrees of latitude to kilometers (≈ Earth circumference / 360°). */ +const KM_PER_DEG_LAT = 111; + +/** + * Likelihood thresholds in radians of sun-altitude error. + * + * ULTRA_TIGHT ≈ 1.43° TIGHT ≈ 2.29° + * MAIN ≈ 2.86° VISIBLE ≈ 4.30° + * + * Tuned to roughly match the band sizes of the older + * |ratio/tan(alt) − 1| ≤ {0.05, 0.08, 0.10, 0.15} thresholds at typical + * mid-latitude altitudes. */ +export const ULTRA_TIGHT_BAND_RAD = 0.025; +export const TIGHT_BAND_RAD = 0.04; +export const MAIN_BAND_RAD = 0.05; +export const VISIBLE_BAND_RAD = 0.075; + +/** Sentinel for grid points that don't have a valid likelihood (sun below horizon, + * or — when the grid is generated with an azimuth constraint — outside its tolerance). */ +export const NIGHT_LIKELIHOOD = -1; + +/** Default tolerance applied by the UI when wiring up an azimuth constraint. */ +export const DEFAULT_AZIMUTH_TOLERANCE_DEG = 10; + +/** Percentile of the great-circle distance distribution used as the "accuracy" radius. */ +const ACCURACY_PERCENTILE = 0.68; + +// Minimum lat gap (degrees) that triggers bimodal-posterior splitting. Shadow +// solutions usually come as one northern + one southern crescent separated by +// a clear gap — anything below this is treated as a single cluster. +const MIN_LAT_CLUSTER_GAP_DEG = 20; -import * as SunCalc from 'suncalc'; +// ── Public types ──────────────────────────────────────────────────────────── export interface ShadowFinderPoint { lat: number; lng: number; + /** Radians of angular error between measured and predicted sun altitude, or NIGHT_LIKELIHOOD (-1) for night points. */ likelihood: number; + /** Sun's compass bearing at this grid point, in degrees [0, 360), 0 = North. */ sunAzimuthDeg: number; } export interface ShadowAnalysisInput { - objectHeight: number; // in pixels - shadowLength: number; // in pixels - knownTime: Date; // UTC time + objectHeight: number; + shadowLength: number; + knownTime: Date; } export interface ShadowAnalysisResult { points: ShadowFinderPoint[]; + /** + * Bounding box of the dominant solution cluster (after splitting bimodal + * north/south posteriors). When `lngWraps` is true, the band crosses the + * antimeridian; `lngRange[0]` is the eastern edge and `lngRange[1]` is the + * western edge (so `lngRange[0] > lngRange[1]`). + */ mainBandCoordinates: { latRange: [number, number]; lngRange: [number, number]; + lngWraps: boolean; }; statistics: { totalPoints: number; @@ -39,163 +107,256 @@ export interface ShadowAnalysisResult { }; } +// ── Sun position ──────────────────────────────────────────────────────────── + +interface SunDayConstants { + sinDec: number; + cosDec: number; + tanDec: number; + ra: number; // right ascension (rad) + gmst: number; // Greenwich mean sidereal time at this instant (rad) +} + /** - * Compute the ShadowFinder likelihood score for a single grid point. - * Returns the absolute relative difference between the shadow ratio implied - * by the sun's altitude and the measured shadow ratio. + * Precompute the time-only sun parameters. Formulas from the Astronomical + * Almanac, identical to those used by SunCalc — we just split them out so the + * grid loop pays the cost once instead of 208,800 times. + */ +function precomputeSun(time: Date): SunDayConstants { + const d = time.valueOf() / DAY_MS - 0.5 + J1970 - J2000; + const M = RAD_PER_DEG * (357.5291 + 0.98560028 * d); + const C = RAD_PER_DEG * (1.9148 * Math.sin(M) + 0.02 * Math.sin(2 * M) + 0.0003 * Math.sin(3 * M)); + const L = M + C + RAD_PER_DEG * 102.9372 + Math.PI; + const sinL = Math.sin(L); + const cosL = Math.cos(L); + const sinE = Math.sin(OBLIQUITY); + const cosE = Math.cos(OBLIQUITY); + const dec = Math.asin(sinL * sinE); + const ra = Math.atan2(sinL * cosE, cosL); + const gmst = RAD_PER_DEG * (280.16 + 360.9856235 * d); + return { + sinDec: Math.sin(dec), + cosDec: Math.cos(dec), + tanDec: Math.tan(dec), + ra, + gmst, + }; +} + +// ── Public functions ──────────────────────────────────────────────────────── + +/** + * Symmetric angular error between the sun altitude implied by the measured + * shadow ratio and the sun altitude predicted at a grid point. * - * likelihood = |measuredRatio / tan(sunAltitudeRad) - 1| + * likelihood = |atan(height / shadow) − sunAltitudeRad| * - * A value of 0 means the sun altitude perfectly explains the measured shadow. - * Only valid when sunAltitudeRad > 0 (daytime points). + * Lower = better fit; 0 = perfect. Swapping the roles of height and shadow + * at altitude α produces the same error around π/4 because atan(x) and + * atan(1/x) are reflections across π/4. + * + * Only meaningful when sunAltitudeRad > 0 (daytime). */ export function computeLikelihood( objectHeight: number, shadowLength: number, - sunAltitudeRad: number + sunAltitudeRad: number, ): number { - return Math.abs((objectHeight / shadowLength) / Math.tan(sunAltitudeRad) - 1); + return Math.abs(Math.atan(objectHeight / shadowLength) - sunAltitudeRad); } /** - * Generate ShadowFinder grid using the exact algorithm approach - * This matches the 290x720 grid structure from the reference implementation + * Generate the full 290 × 720 = 208,800-point world grid for a given instant. + * Each point carries either an angular-error likelihood (sun above horizon) or + * `NIGHT_LIKELIHOOD` (-1) when the sun is below the horizon — or, when an + * `azimuthConstraint` is provided, when the sun's azimuth at that cell is + * outside the constraint's tolerance. + * + * Optimizations: + * - `Math.atan2` for sun azimuth is skipped on night cells (saves ~50% of + * atan2 calls). The `sunAzimuthDeg` field is set to NaN on those cells; + * `applyAzimuthConstraint` correctly excludes them because NaN comparisons + * are always false. + * - When `azimuthConstraint` is provided, azimuth is computed first and cells + * outside tolerance are marked `NIGHT_LIKELIHOOD` immediately — skipping + * the `Math.asin` altitude resolution for ~95% of cells in typical use. + * CAVEAT: callers that rely on the "no-match fallback" in the visualization + * (showing the full shadow band when the constraint excludes everything) + * must NOT pass the constraint here. */ export function generateShadowFinderGrid( knownTime: Date, objectHeight: number, - shadowLength: number + shadowLength: number, + azimuthConstraint?: AzimuthConstraint | null, ): ShadowFinderPoint[] { - const points: ShadowFinderPoint[] = []; - const angularResolution = 0.5; // degrees - same as ShadowFinder - - // Sample points across the world (EXACT same range as ShadowFinder) - for (let lat = -60.0; lat <= 84.5; lat += angularResolution) { - for (let lng = -180.0; lng <= 179.5; lng += angularResolution) { - const sunPos = SunCalc.getPosition(knownTime, lat, lng); - const sunAltitudeRad = sunPos.altitude; - - let likelihood: number; - - // ShadowFinder's exact approach: set night areas to -1 - if (sunAltitudeRad <= 0) { - likelihood = -1; // Night area - will be filtered out in visualization - } else { - likelihood = computeLikelihood(objectHeight, shadowLength, sunAltitudeRad); + const sun = precomputeSun(knownTime); + const measuredAlt = Math.atan(objectHeight / shadowLength); + + // Hour angle H depends on longitude only: H = gmst + lng·rad − ra. Both + // altitude and azimuth formulas use only sin(H) and cos(H), so caching them + // per longitude column saves one sin and one cos per cell (≈ 200K × 2 = 400K + // trig calls per analysis). + const sinH = new Float64Array(LNG_STEPS); + const cosH = new Float64Array(LNG_STEPS); + for (let j = 0; j < LNG_STEPS; j++) { + const H = sun.gmst + (LNG_MIN + j * GRID_RES_DEG) * RAD_PER_DEG - sun.ra; + sinH[j] = Math.sin(H); + cosH[j] = Math.cos(H); + } + + const points = new Array(GRID_COUNT); + let idx = 0; + + const constraintActive = !!azimuthConstraint?.enabled; + const targetBearing = azimuthConstraint?.sunBearingDeg ?? 0; + const tolerance = azimuthConstraint?.toleranceDeg ?? 0; + + for (let i = 0; i < LAT_STEPS; i++) { + const lat = LAT_MIN + i * GRID_RES_DEG; + const phi = lat * RAD_PER_DEG; + const sinPhi = Math.sin(phi); + const cosPhi = Math.cos(phi); + const sinPhiSinDec = sinPhi * sun.sinDec; + const cosPhiCosDec = cosPhi * sun.cosDec; + const tanDecCosPhi = sun.tanDec * cosPhi; + + for (let j = 0; j < LNG_STEPS; j++) { + const lng = LNG_MIN + j * GRID_RES_DEG; + const sH = sinH[j]; + const cH = cosH[j]; + + const sinAlt = sinPhiSinDec + cosPhiCosDec * cH; + + // Fast path for night points: skip both Math.asin (altitude) and + // Math.atan2 (azimuth). NaN azimuth ensures applyAzimuthConstraint + // excludes these cells (comparison with NaN is always false). + if (sinAlt <= 0) { + points[idx++] = { lat, lng, likelihood: NIGHT_LIKELIHOOD, sunAzimuthDeg: NaN }; + continue; } - // SunCalc azimuth: 0=South, positive=West, negative=East (radians) - // Convert to 0–360 compass bearing (0=North, clockwise) - const sunAzimuthDeg = (sunPos.azimuth * 180 / Math.PI + 180 + 360) % 360; + // SunCalc azimuth convention: 0 = South, positive = West, radians. + // Compass bearing: ((az · 180/π) + 180 + 360) mod 360 → 0 = North, CW. + const azimuthRad = Math.atan2(sH, cH * sinPhi - tanDecCosPhi); + const sunAzimuthDeg = (azimuthRad * DEG_PER_RAD + 540) % 360; + + // When a constraint is baked in at generation time, reject cells outside + // tolerance BEFORE paying for the altitude asin. Same wrap-aware delta + // math as applyAzimuthConstraint. + if (constraintActive) { + const diff = Math.abs(((sunAzimuthDeg - targetBearing + 540) % 360) - 180); + if (diff > tolerance) { + points[idx++] = { lat, lng, likelihood: NIGHT_LIKELIHOOD, sunAzimuthDeg }; + continue; + } + } - points.push({ - lat: lat, - lng: lng, - likelihood: likelihood, - sunAzimuthDeg, - }); + const likelihood = Math.abs(measuredAlt - Math.asin(sinAlt)); + points[idx++] = { lat, lng, likelihood, sunAzimuthDeg }; } } return points; } +export interface ShadowAnalysisOptions { + /** + * If provided and enabled, the grid generation pre-filters by sun azimuth + * (skipping the per-cell `Math.asin` for ~95% of cells in typical use). + * + * CAVEAT: callers relying on the "no-match fallback" — showing the full + * shadow band when the constraint excludes everything — must NOT pass the + * constraint here, because the full band won't exist in the result. The + * production callsite in `Index.tsx` leaves this unset and applies the + * constraint client-side via `applyAzimuthConstraint` instead. + */ + azimuthConstraint?: AzimuthConstraint | null; +} + /** - * Analyze shadow measurements and generate location possibilities + * Analyze shadow measurements and generate location possibilities. */ -export function analyzeShadowMeasurements(input: ShadowAnalysisInput): ShadowAnalysisResult { - // Validate inputs +export function analyzeShadowMeasurements( + input: ShadowAnalysisInput, + options?: ShadowAnalysisOptions, +): ShadowAnalysisResult { if (input.objectHeight <= 0 || input.shadowLength <= 0) { throw new Error('Invalid measurements: object height and shadow length must be positive'); } - if (isNaN(input.knownTime.getTime())) { throw new Error('Invalid date/time provided'); } - - // Generate grid data using ShadowFinder's exact approach - const points = generateShadowFinderGrid(input.knownTime, input.objectHeight, input.shadowLength); - // Single pass over all 208,800 points: collect counts and running min/max. - // Replaces 7 separate .filter() passes that previously visited ~1.46M elements. - let nightPoints = 0, validPoints = 0; - let ultraTightBandPoints = 0, tightBandPoints = 0, visibleBandPoints = 0; - let minLat = Infinity, maxLat = -Infinity, minLng = Infinity, maxLng = -Infinity; - let hasMainBand = false; + const points = generateShadowFinderGrid( + input.knownTime, + input.objectHeight, + input.shadowLength, + options?.azimuthConstraint, + ); + + let nightPoints = 0; + let validPoints = 0; + let ultraTightBandPoints = 0; + let tightBandPoints = 0; + let visibleBandPoints = 0; + const mainBand: Array<{ lat: number; lng: number }> = []; for (const p of points) { const l = p.likelihood; - if (l === -1) { + if (l === NIGHT_LIKELIHOOD) { nightPoints++; continue; } validPoints++; - if (l <= 0.05) ultraTightBandPoints++; - if (l <= 0.08) tightBandPoints++; - if (l <= 0.15) visibleBandPoints++; - if (l <= 0.1) { - hasMainBand = true; - if (p.lat < minLat) minLat = p.lat; - if (p.lat > maxLat) maxLat = p.lat; - if (p.lng < minLng) minLng = p.lng; - if (p.lng > maxLng) maxLng = p.lng; - } + if (l <= ULTRA_TIGHT_BAND_RAD) ultraTightBandPoints++; + if (l <= TIGHT_BAND_RAD) tightBandPoints++; + if (l <= VISIBLE_BAND_RAD) visibleBandPoints++; + if (l <= MAIN_BAND_RAD) mainBand.push({ lat: p.lat, lng: p.lng }); } - const latRange: [number, number] = hasMainBand ? [minLat, maxLat] : [0, 0]; - const lngRange: [number, number] = hasMainBand ? [minLng, maxLng] : [0, 0]; - return { points, - mainBandCoordinates: { - latRange, - lngRange - }, + mainBandCoordinates: computeMainBandBbox(mainBand), statistics: { totalPoints: points.length, validPoints, nightPoints, ultraTightBandPoints, tightBandPoints, - visibleBandPoints - } + visibleBandPoints, + }, }; } /** - * Calculate object height and shadow length from pixel coordinates + * Calculate object height and shadow length from pixel coordinates. */ export function calculateMeasurementsFromPixels( objectBottom: { x: number; y: number }, objectTop: { x: number; y: number }, - shadowTip: { x: number; y: number } + shadowTip: { x: number; y: number }, ): { objectHeight: number; shadowLength: number } { - // Calculate object height in pixels const objectHeight = Math.sqrt( - Math.pow(objectTop.x - objectBottom.x, 2) + - Math.pow(objectTop.y - objectBottom.y, 2) + Math.pow(objectTop.x - objectBottom.x, 2) + Math.pow(objectTop.y - objectBottom.y, 2), ); - - // Calculate shadow length in pixels const shadowLength = Math.sqrt( - Math.pow(shadowTip.x - objectBottom.x, 2) + - Math.pow(shadowTip.y - objectBottom.y, 2) + Math.pow(shadowTip.x - objectBottom.x, 2) + Math.pow(shadowTip.y - objectBottom.y, 2), ); - return { objectHeight, shadowLength }; } /** - * Convert percentage-based coordinates to absolute pixel coordinates + * Convert percentage-based coordinates to absolute pixel coordinates. */ export function convertPercentageToPixels( percentageCoords: { x: number; y: number }, imageWidth: number, - imageHeight: number + imageHeight: number, ): { x: number; y: number } { return { x: (percentageCoords.x / 100) * imageWidth, - y: (percentageCoords.y / 100) * imageHeight + y: (percentageCoords.y / 100) * imageHeight, }; } @@ -206,60 +367,146 @@ export interface AzimuthConstraint { } /** - * Filter grid points to those whose sun azimuth is within toleranceDeg of - * the observed sun bearing. Handles the 0°/360° wraparound. + * Filter grid points to those whose sun azimuth is within toleranceDeg of the + * observed sun bearing. Handles the 0°/360° wraparound. */ export function applyAzimuthConstraint( points: ShadowFinderPoint[], - constraint: AzimuthConstraint + constraint: AzimuthConstraint, ): ShadowFinderPoint[] { if (!constraint.enabled) return points; - return points.filter(p => { + return points.filter((p) => { const diff = Math.abs(((p.sunAzimuthDeg - constraint.sunBearingDeg + 540) % 360) - 180); return diff <= constraint.toleranceDeg; }); } /** - * Estimate the best possible location from analysis results + * Estimate the best possible location from analysis results. + * + * Resolves bimodal posteriors by splitting on the largest latitude gap, then + * computes a circular-mean centroid (correct across the antimeridian) and a + * 68th-percentile great-circle distance as the "accuracy" radius. */ export function estimateBestLocation(result: ShadowAnalysisResult): { latitude: number; longitude: number; accuracy: number; } { - // Find points with the best likelihood (lowest values) - const bestPoints = result.points - .filter(p => p.likelihood >= 0 && p.likelihood <= 0.1) - .sort((a, b) => a.likelihood - b.likelihood); - + const bestPoints = result.points.filter( + (p) => p.likelihood >= 0 && p.likelihood <= MAIN_BAND_RAD, + ); if (bestPoints.length === 0) { throw new Error('No suitable location matches found'); } - - // Single pass: centroid + running min/max (spread operator on large arrays risks stack overflow) - let sumLat = 0, sumLng = 0; - let minLat = Infinity, maxLat = -Infinity, minLng = Infinity, maxLng = -Infinity; - for (const p of bestPoints) { - sumLat += p.lat; sumLng += p.lng; + + const cluster = dominantLatCluster(bestPoints); + + // Circular mean of longitudes — averaging unit vectors on the lng circle + // gives the right answer at the ±180° seam (e.g. lng=-179 and lng=179 → ±180, + // not 0). + let sumLat = 0; + let sumLngX = 0; + let sumLngY = 0; + for (const p of cluster) { + sumLat += p.lat; + const lngRad = p.lng * RAD_PER_DEG; + sumLngX += Math.cos(lngRad); + sumLngY += Math.sin(lngRad); + } + const avgLat = sumLat / cluster.length; + const avgLng = Math.atan2(sumLngY, sumLngX) * DEG_PER_RAD; + + // 68th-percentile great-circle distance is more robust than min/max spread + // because it ignores outliers and degrades gracefully on single-point clusters. + // Floor at half the grid resolution at the equator — a 0.5° grid simply + // cannot resolve finer than half a cell (~28 km). + const distances = cluster + .map((p) => greatCircleKm(avgLat, avgLng, p.lat, p.lng)) + .sort((a, b) => a - b); + const pIdx = Math.min(distances.length - 1, Math.floor(distances.length * ACCURACY_PERCENTILE)); + const halfCellKm = (GRID_RES_DEG / 2) * KM_PER_DEG_LAT; + const accuracy = Math.max(distances[pIdx], halfCellKm); + + return { latitude: avgLat, longitude: avgLng, accuracy }; +} + +// ── Internal helpers (exported for tests) ─────────────────────────────────── + +/** + * Pick the larger of the two latitudinal clusters when there's a clear gap; + * otherwise return the input. Operates on a sorted copy and runs in O(n log n). + */ +export function dominantLatCluster(points: T[]): T[] { + if (points.length <= 1) return points; + const sorted = [...points].sort((a, b) => a.lat - b.lat); + let largestGap = 0; + let gapIdx = -1; + for (let i = 1; i < sorted.length; i++) { + const g = sorted[i].lat - sorted[i - 1].lat; + if (g > largestGap) { + largestGap = g; + gapIdx = i; + } + } + if (largestGap < MIN_LAT_CLUSTER_GAP_DEG || gapIdx <= 0) return sorted; + const south = sorted.slice(0, gapIdx); + const north = sorted.slice(gapIdx); + return south.length >= north.length ? south : north; +} + +function computeMainBandBbox( + band: Array<{ lat: number; lng: number }>, +): ShadowAnalysisResult['mainBandCoordinates'] { + if (band.length === 0) { + return { latRange: [0, 0], lngRange: [0, 0], lngWraps: false }; + } + + const dominant = dominantLatCluster(band); + + let minLat = Infinity; + let maxLat = -Infinity; + const lngs: number[] = []; + for (const p of dominant) { if (p.lat < minLat) minLat = p.lat; if (p.lat > maxLat) maxLat = p.lat; - if (p.lng < minLng) minLng = p.lng; - if (p.lng > maxLng) maxLng = p.lng; + lngs.push(p.lng); } - const avgLat = sumLat / bestPoints.length; - const avgLng = sumLng / bestPoints.length; - - // Calculate accuracy based on spread of points - const latSpread = maxLat - minLat; - const lngSpread = maxLng - minLng; - - // Convert degrees to approximate km (rough estimation) - const accuracy = Math.max(latSpread * 111, lngSpread * 111 * Math.cos(avgLat * Math.PI / 180)) / 2; - - return { - latitude: avgLat, - longitude: avgLng, - accuracy: Math.max(accuracy, 1) // Minimum 1km accuracy - }; -} \ No newline at end of file + + const { lngRange, wraps } = lngBboxWithWrap(lngs); + return { latRange: [minLat, maxLat], lngRange, lngWraps: wraps }; +} + +/** + * Tight lng bbox that respects the antimeridian. Strategy: sort the lngs, find + * the largest empty arc (including the wraparound arc from max back to min), + * and the cluster spans from the point AFTER the gap to the point BEFORE it. + * If the chosen gap is the wrap arc, the band does not wrap; otherwise it does. + */ +function lngBboxWithWrap(lngs: number[]): { lngRange: [number, number]; wraps: boolean } { + const sorted = [...lngs].sort((a, b) => a - b); + const n = sorted.length; + let largestGap = sorted[0] + 360 - sorted[n - 1]; + let gapStartIdx = n - 1; + for (let i = 1; i < n; i++) { + const g = sorted[i] - sorted[i - 1]; + if (g > largestGap) { + largestGap = g; + gapStartIdx = i - 1; + } + } + const start = sorted[(gapStartIdx + 1) % n]; + const end = sorted[gapStartIdx]; + return { lngRange: [start, end], wraps: start > end }; +} + +function greatCircleKm(lat1: number, lng1: number, lat2: number, lng2: number): number { + const phi1 = lat1 * RAD_PER_DEG; + const phi2 = lat2 * RAD_PER_DEG; + const dPhi = (lat2 - lat1) * RAD_PER_DEG; + const dLambda = (lng2 - lng1) * RAD_PER_DEG; + const a = + Math.sin(dPhi / 2) ** 2 + + Math.cos(phi1) * Math.cos(phi2) * Math.sin(dLambda / 2) ** 2; + return 2 * EARTH_RADIUS_KM * Math.asin(Math.min(1, Math.sqrt(a))); +}