Geolocating a Random Island via Geometry and CUDA Programming

Author: gralhix #004 yassa9 gralhix004 | Date: 16-08-2026 Note: This is a genuine human effort; no LLM generation was used in the original work.

This document serves as a comprehensive write-up for the challenge gralhix 004, created by Sofia Santos | Gralhix.

Task Briefing: The provided image depicts a resort situated on an island. The primary objective is to determine the exact location and answer: In which cardinal direction was the camera facing when the photo was taken?

While using Google Lens would be the easy way out, I felt that would waste a great opportunity for a technical challenge. Instead, I opted to solve this using a combination of mathematics and high-performance programming.

You can access the full source code, local execution files, and the final report on GitHub.

---

🛠 Project Checklist

• Extract Metadata
• Create Geometric Fingerprint
• Filter Global Landmass Dataset
• Implement CUDA Matching Kernel
• Deduplicate Results

---

Step A: Metadata Analysis

The first step in any OSINT challenge is checking the metadata.

Result: As suspected, the metadata contained no GPS coordinates or useful timestamps.

---

Step B: Establishing the Fingerprint

Upon analyzing the image, I identified three distinct landmasses:

• P0: The primary islet containing the resort.
• P1: The island to the right.
• P2: The island in the left foreground (featuring a mountain peak).

Because the photo was taken by a drone, the elevation is unknown and the perspective is skewed. I couldn't build a perfect bird's-eye model, so I relied on intuition to determine the relative distances and the internal angles of the triangle formed by these three points.

I developed a custom GUI tool, 01_triangle_gui.py, to record pixel coordinates and compute the geometry. To account for human error when clicking the center of the islands, I implemented a $\pm 20\%$ tolerance band for the search parameters.

---

Step C: The Global Search

With the geometric fingerprint locked, I needed to compare it against every landmass on Earth.

Dataset: I used the land-polygons-split-4326 set from OpenStreetMap, which contains global coastline vectors in WGS84 (approx. 882 MB).

The Filtering Pipeline

I spent several days tweaking heuristic filters through trial and error to narrow down the candidates.

1. Tropical Latitude Bounding Box Since the vegetation and architecture suggest a tropical climate, I discarded everything outside: $-30^\circ \le \text{latitude} \le 30^\circ$

2. Local Density Filter To avoid crowded archipelagos or dense coastlines, I defined $N_{5\text{km}}(p)$ as the number of centroids within 5km of point $p$. $N_{5\text{km}}(p) \le 10$ If an islet has more than 10 neighbors, it is too crowded to match the isolated nature of the photo.

3. Clustering I searched for points within a 20km radius. A "cluster" is defined as any group containing at least 3 points: $\left|\{q : \text{dist}(p,q) \le 20\,\text{km}\}\right| \ge 3$ Using query_ball_point(f_coords, CLUSTER_RADIUS_KM / 111.0), this reduced the search space to 23,500 clusters.

4. Generating Triplets For each cluster, I calculated all possible combinations of 3 points. The number of combinations is given by: $\binom{n}{3} = \frac{n(n-1)(n-2)}{6}$ To prevent a combinatorial explosion (e.g., 60 points $\approx$ 34,220 triples), I capped each cluster at 60 points using a stratified_sample function. This ensured a mix of small, medium, and large islands.

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.concatenate([
        order[:n_small], 
        order[-n_large:], 
        order[mid_start : mid_start + n_mid]
    ])
    return idx_arr[keep]

5. GPU Matching via CUDA I assigned every single triple to its own CUDA thread for parallel processing.

• Identification: Threads sorted points by area to identify P0 (the smallest).
• Orientation: A 2D cross product determined the assignment of P1 and P2: $\text{cross} = x_a y_b - x_b y_a$ $P1 = \begin{cases} a & \text{cross} > 0 \\ b & \text{cross} \le 0 \end{cases}$
• Geometry: The thread then calculated the angle $\theta_0$ and the distance ratio $r$: $\theta_0 = \arccos\left(\frac{\vec{d_1} \cdot \vec{d_2}}{|\vec{d_1}||\vec{d_2}|}\right), \qquad r = \frac{|\vec{d_1}|}{|\vec{d_2}|}$

If the angle, ratio, and sizes fell within the tolerance windows, the result was saved using an atomicAdd counter to prevent race conditions.

if (hit) {
    unsigned long long slot = atomicAdd(out_count, 1ULL);
    out_p0[slot] = p0idx;
    out_p1[slot] = p1idx;
    out_p2[slot] = p2idx;
}

Performance Metrics:

• Hardware: NVIDIA GeForce RTX 3050 (sm_86)
• VRAM Used: 5169 MB
• Kernel Execution Time: 204.1 ms
• Throughput: 80.7 million triples processed in parallel.

6. Deduplication Since a single physical triple could be processed by multiple threads if it belonged to overlapping clusters, a final deduplication step was performed...