Het geoloceren van een willekeurige afbeelding van een eilandje met behulp van geometrie en CUDA GPU-programmering

Opdracht

Er is een foto van een resort op een eiland. De volgende vragen moeten beantwoord worden: a) Wat is de naam van het resort? b) Wat zijn de coördinaten van het eiland? c) In welke windrichting keek de camera toen de foto werd genomen?

Om deze uitdaging niet simpelweg met Google Lens op te lossen, is er gekozen voor een aanpak op basis van wiskunde en programmering.

a] Metadata

De eerste stap is het controleren van de metadata. Met exiftool op Linux Void werd het volgende resultaat verkregen:

> exiftool main.png
File Type                       : WEBP (lossless)
MIME Type                       : image/webp
Image Width                     : 736
Image Height                    : 515

Zoals verwacht is er geen bruikbare informatie aanwezig: geen EXIF, geen GPS en geen gegevens over het cameramodel.

b] Het maken van het kenmerkprofiel (fingerprint)

Op de afbeelding zijn drie landmassa's zichtbaar:

• P0: het eilandje zelf.
• P1: het rechter eiland.
• P2: het linker voorste eiland (met een bergtop).

Omdat de foto met een drone is genomen en de hoogte niet in de metadata staat, kon er geen exact perspectiefmodel van bovenaf worden gemaakt. De focus ligt daarom op de relatieve afstanden tussen de drie eilanden en de hoeken van de gevormde driehoek.

Met een kleine GUI (01trianglegui.py) zijn de pixelcoördinaten van elk punt vastgelegd om de geometrie van de driehoek te berekenen. Er is een tolerantieband van ±20% toegevoegd bij het zoeken, omdat het handmatig aanklikken van de exacte centra niet perfect precies is.

c] Zoeken

Vervolgens worden alle landmassa's op aarde vergeleken met dit profiel. Als dataset is het land-polygons-split-4326 set van OpenStreetMap gebruikt (globale kustlijn-vectoren in WGS84, 882 MB).

Om de hoeveelheid data te reduceren, zijn de volgende heuristische filters toegepast:

01] Tropische breedtegraad-begrenzing

$$ -30^\circ \le \text{breedtegraad} \le 30^\circ $$ Het eilandje ziet er tropisch uit, dus alles buiten de tropen wordt direct verwijderd. Hiermee blijven er 141.131 landpolygonen over.

02] Filter voor lokale dichtheid

$$ N{5\text{km}}(p) \le 10 $$ $N{5\text{km}}(p)$ telt hoeveel andere centroids zich binnen 5 km van punt (p) bevinden. Als een eilandje meer dan 10 buren in deze straal heeft, bevindt het zich waarschijnlijk in een dicht rifveld of een archipel, en niet in een kleine, geïsoleerde groep van 3-4 eilanden. Dit reduceert het aantal kandidaten tot 51.576.

03] Clustering

Voor elk overgebleven punt wordt gezocht naar andere punten binnen een straal van 20 km. Alleen punten die deel uitmaken van een cluster van minimaal 3 punten worden behouden, aangezien zij anders geen driehoek kunnen vormen.

tree = cKDTree(f_coords)
neigh = tree.query_ball_point(
    f_coords,
    CLUSTER_RADIUS_KM / 111.0)
clusters = set(tuple(sorted(n)) for n in neigh if len(n) >= 3)

$$ \left|\{q : \text{dist}(p,q) \le 20\,\text{km}\}\right| \ge 3 $$ Dit resulteert in 23.500 clusters.

04] Genereren van tripletten

Elke combinatie van 3 punten binnen een cluster wordt een kandidaat-driehoek. Dit aantal combinaties $C(n, 3)$ groeit snel. Om dit beheersbaar te houden, is elk cluster beperkt tot 60 punten via een gestratificeerde steekproef (een derde kleine eilanden, een derde grote, en een derde gemiddelde eilanden).

$$ \binom{n}{3} = \frac{n(n-1)(n-2)}{6} $$

def stratified_sample(idx_arr, area_arr, cap):
    order = np.argsort(area_arr[idx_arr])
    n_small = cap // 3
    n_large = cap // 3
    n_mid = cap - n_small - n_large
    mid_start = max(0, (len(idx_arr) - n_large - n_mid) // 2)
    keep = np.unique(np.concatenate([
        order[:n_small],
        order[-n_large:],
        order[mid_start:mid_start + n_mid],
    ]))
    return idx_arr[keep]

def gen_cluster_triples(idx_arr):
    local = np.array(list(
        itertools.combinations(range(len(idx_arr)), 3)),
        dtype=np.int64)
    return idx_arr[local]

In totaal worden er 80.690.777 tripletten gegenereerd.

05] Matching op de GPU

Met CUDA wordt elk triplett toegewezen aan één thread. De thread sorteert de 3 punten op oppervlakte om P0 (het kleinste, het resort-eilandje) te bepalen. De draairichting van de andere twee punten bepaalt P1 en P2 via een 2D-kruisproduct:

$$ \text{cross} = xa yb - xb ya $$ $$ P1 = \begin{cases} a & \text{cross} > 0 \\ b & \text{cross} \le 0 \end{cases} $$

Vervolgens worden de hoek bij P0 en de afstandratio berekend: $$ \theta0 = \arccos\left(\frac{\vec{d1} \cdot \vec{d2}}{|\vec{d1}||\vec{d2}|}\right), \qquad r = \frac{|\vec{d1}|}{|\vec{d_2}|} $$

Een triplett overleeft als de hoek, de ratio, de grootte van P0, de afstand tussen P0 en P1, en beide zijdelengtes binnen de tolerantievensters vallen.

GPU Statistieken:

• GPU: NVIDIA GeForce RTX 3050 (sm_86)
• VRAM gebruik: 5169 MB
• Kernel tijd: 204.1 ms
• Resultaat: 158.784 tripletten passeren het masker.

06] Deduplicatie

Omdat hetzelfde fysieke triplett door meerdere overlappende clusters kan worden geraakt, worden de resultaten gefilterd op unieke identiteit. Dit brengt het aantal terug naar 8.915 unieke tripletten.

seen = set()
uniq = []
for i in range(len(p0_all)):
    key = (p0_all[i], p1_all[i], p2_all[i])
    if key not in seen:
        seen.add(key)
        uniq.append(i)

07] De open rechthoek

Elk overgebleven triplett wordt getest op de aanwezigheid van open water. Er wordt een rechthoek gebouwd langs de rand P0 → P1 (aan de zijde waar P2 zich niet bevindt). Als er andere landmassa's in deze rechthoek liggen, wordt de kandidaat verworpen.

width = np.hypot(x1, y1)
u = np.array([x1, y1]) / width
v = np.array([-u[1], u[0]])
# p2 zit door constructie aan de +v zijde,
# dus de check vindt plaats aan de -v zijde
length = 2 * width
corners_local = [
    (0, 0), (x1, y1),
    (x1 - v[0]*length, y1 - v[1]*length),
    (-v[0]*length, -v[1]*length),
]

Dit reduceert de lijst van 8.915 naar 948 kandidaten.

d] Controle van de vorm van het koraaleiland

In deze fase wordt alleen gekeken naar P0 om te bepalen of de vorm overeenkomt met een koraaleiland (coral cay).

1] Compactheid (Polsby Popper Score): $$ PP = \frac{4\pi \cdot \text{oppervlakte}}{\text{omtrek}^2} $$ Een score van 1.0 is een perfecte cirkel. Koraaleilanden zijn vaak rond door golfafzetting; waarden onder de 0.5 worden verwijderd.

2] Micro Cay Halo Check: Er wordt geteld hoeveel kleine landfragmenten (kleiner dan 0.05 km²) zich binnen 1,5 km van P0 bevinden. Echte rifsystemen hebben vaak kleine zandbanken rond het hoofdeiland. Er moet er minimaal één aanwezig zijn.

Na deze checks blijven er 213 van de 948 kandidaten over.

e] Controle van de ovale vorm

Een geometrisch filter op het polygoon van P0. Er wordt een minimale roterende rechthoek om de vorm geplaatst.

Aspectratio: $$ \text{aspect} = \frac{\text{lange zijde}}{\text{korte zijde}} \in [1.05,\ 2.2] $$ Een waarde te dicht bij 1.0 is een perfecte cirkel; een waarde boven 2.2 is te langgerekt.

Fill ratio (Vulgraad): Een perfecte ellips vult precies $\pi / 4$ van zijn omschrijvende rechthoek. $$ \frac{\text{oppervlakte}{\text{ellipse}}}{\text{oppervlakte}{\text{box}}} = \frac{\pi}{4} \approx 0.785 $$ De drempelwaarde is vastgesteld op 75% van dit maximum: $$ \text{FILL\RATIO\MIN} = 0.75 \times \frac{\pi}{4} \approx 0.589 $$ Vormen zoals halve manen of ringen vallen hieronder, solide afgeronde eilandjes niet.

137 van de 213 kandidaten overleven deze check.

f] NDVI Vegetatie Check

Via de publieke STAC API van Earth Search (Element84) wordt Sentinel-2 satellietbeeldmateriaal gebruikt om te controleren of P0 begroeid is met palmbomen en niet uit kaal zand of rots bestaat.

$$ \text{NDVI} = \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}} $$

Gezonde vegetatie reflecteert sterk in het nabij-infrarood (NIR) en absorbeert rood licht. Een NDVI-drempelwaarde van 0.6 is gehanteerd.

66 van de 137 kandidaten overleven de NDVI-check.

g] Controle van hoogte en bergen

De laatste check hanteert twee condities:

1. P0 zelf moet laag en vlak zijn (consistent met een rifeiland).
2. P2 moet werkelijke hoogte hebben in de richting waar de camera naar keek.

De kijkrichting ($\theta{\text{front}}$) is de bissectrice tussen de richting naar P1 en de richting naar P2: $$ \theta(P0, Pi) = \text{atan2}\Big(\sin(\Delta\lambda)\cos\phii,\ \cos\phi0\sin\phii - \sin\phi0\cos\phii\cos(\Delta\lambda)\Big) $$ $$ \theta{\text{front}} = \theta(P0, P2) + \frac{\big((\theta(P0,P1) - \theta(P0,P_2) + 180) \bmod 360\big) - 180}{2} $$

Er wordt een waaierscan van ±50° rond deze richting uitgevoerd tot een straal van 20 km, waarbij gebruik wordt gemaakt van de Copernicus DEM GLO-30 (30m resolutie).

De overlevingscriteria zijn: $$ \text{hoogte}(P0) \le 50\text{m} $$ $$ 100\text{m} \le \max{\text{boog}}(\text{hoogte}) \le 500\text{m} $$

26 van de 66 kandidaten overleven deze check. De meeste bevinden zich in Zuid-Azië, Australië en Oceanië, met één uitzondering nabij Brazilië.

h] Eindverslag

De overgebleven kandidaten worden in een HTML-tabel gezet met landnaam en Google Maps satellietlinks naar P0, P1 en P2.

Na visuele inspectie van de lijst bleek de achtste kandidaat in de tabel, gelegen in de staat Micronesië, de juiste oplossing te zijn.

i] Antwoorden

a) Wat is de naam van het resort? Oan

b) Wat zijn de coördinaten van het eiland? $7^\circ\,21^\prime\,48.4^{\prime\prime}\,\text{N} \qquad 151^\circ\,45^\prime\,20.7^{\prime\prime}\,\text{E}$ (of $7.363444^\circ,\ 151.755750^\circ$)

c) In welke windrichting keek de camera toen de foto werd genomen? Gebaseerd op de berekening: $P0 = (7.3633,\ 151.755983), \quad P1 = (7.386573,\ 151.739534)$ $\theta = 324.97^\circ \implies$ Noordwest (NW)

j] Data & Licenties

• Kustlijn-polygonen: land-polygons-split-4326 © OpenStreetMap contributors (ODbL 1.0).
• Hoogtegegevens: Copernicus DEM GLO-30 © DLR e.V. en Airbus Defence and Space GmbH.
• Satellietbeelden: Copernicus Sentinel data 2025-2026 via Earth Search (Element 84) on AWS Open Data.
• Landgrenzen: Natural Earth 10m admin-0 (Public Domain).
• Uitdaging & bronfoto: OSINT Exercise #004 door Sofia Santos (gralhix).
• Screenshots: Google Maps / Google Earth.