Geolocating a random island using geometry and CUDA programming

Geolocating a random island using geometry and CUDA programming

使用几何学与 CUDA 编程对随机岛屿进行地理定位

gralhix004 | Geolocating Random Islet Image Using Geometry & CUDA GPU Programming 16-08-2026 gralhix004 | 使用几何学与 CUDA GPU 编程对随机小岛图像进行地理定位 2026年8月16日

NOTE: this is a genuine human work, didnt use LLM generation. I’m writing this page as a writeup for this challenge gralhix 004 made by Sofia Santos | Gralhix. You can view, clone and locally try all code files and the final report with all instructions here at github. 注意:这是纯人工完成的工作,未使用大语言模型生成。我撰写此页面是为了记录由 Sofia Santos | Gralhix 发起的 gralhix 004 挑战。你可以在 GitHub 上查看、克隆并本地运行所有代码文件及包含完整说明的最终报告。

Task briefing: This is a photo of a resort located on an island. a) What is the name of the resort? b) What are the coordinates of the island? c) In which cardinal direction was the camera facing when the photo was taken? In my opinion, solving this challenge with google lens is wasting a fun opportunity, so decided to solve it with math and programming. 任务简报:这是一张位于某岛屿上的度假村照片。a) 度假村的名字是什么?b) 该岛屿的坐标是多少?c) 拍摄照片时相机朝向哪个方位?在我看来,用 Google Lens 解决这个挑战是浪费了一次有趣的实践机会,所以我决定用数学和编程来解决它。

a] Metadata

a] 元数据

Of course, first thing u look for is the metadata. Ran that on my linux void: 当然,你首先要看的就是元数据。我在我的 Linux Void 系统上运行了以下命令:

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

As expected, nothing useful here. No EXIF, no GPS, no camera make or model. 不出所料,这里没有任何有用的信息。没有 EXIF,没有 GPS,也没有相机品牌或型号。

b] Building the fingerprint

b] 构建指纹

U can see from the img, there are 3 landmasses: P0: the islet itself, P1: the right island, P2: the left front island ( having mountain peak ) 从图中可以看到有 3 块陆地:P0:小岛本身,P1:右侧岛屿,P2:左前方岛屿(带有山峰)。

I couldnt make a correct perspective model of birdview of this image, as clearly the image is taken by a drone and cant estimate the elevation at all (and not found in the metadata). So I had to estimate that by intuition, I just want the relative distances between the 3 islands and angles of that triangle. I built a small click GUI 01_triangle_gui.py that records pixel coordinates for each point in order and computes the triangle’s geometry. Since clicking exact centers by eye isn’t perfectly precise, I added a ±20% tolerance band around both values when searching. 我无法为这张图像建立准确的鸟瞰透视模型,因为照片显然是由无人机拍摄的,且完全无法估算海拔(元数据中也没有)。所以我只能凭直觉进行估算,我只需要这 3 个岛屿之间的相对距离和三角形的角度。我编写了一个小型点击 GUI 程序 01_triangle_gui.py,按顺序记录每个点的像素坐标并计算三角形的几何形状。由于肉眼点击中心点并不完全精确,我在搜索时为这两个数值添加了 ±20% 的容差范围。

c] 搜索

With the fingerprint locked in, the next step is checking every real landmass on Earth against it ! I used OpenStreetMap’s split land polygon set as the dataset land-polygons-split-4326, full global coastline vectors in WGS84 which has size of 882 MB. I created heuristic filters (all by just intuition and non tangible proofs), spent days (yea full days) tweaking values and tons of trial and error 😭 untill I got this working filters recipe. 指纹确定后,下一步就是将地球上每一块真实的陆地与它进行比对!我使用了 OpenStreetMap 的分割陆地多边形数据集 land-polygons-split-4326,这是 WGS84 坐标系下的全球海岸线矢量数据,大小为 882 MB。我创建了启发式过滤器(完全基于直觉和非实证),花费了数天(真的是整整几天)来调整数值,经过无数次试错 😭,最终得到了这套有效的过滤方案。

01] Tropical latitude bounding box 01] 热带纬度边界框 $$ -30° \le latitude \le 30° $$ the islet in the photo reads as tropical, so I decided that anything outside the tropics is thrown out immediately, before doing any expensive geometry work. Exactly 141,131 land polygons survive that band filter. 照片中的小岛看起来属于热带,所以我决定在进行任何昂贵的几何计算之前,先将热带以外的所有区域直接剔除。最终有 141,131 个陆地多边形通过了该纬度过滤。

02] Local density filter 02] 局部密度过滤器 $$ N_{5\text{km}}(p) \le 10 $$ $ N_{5\text{km}}(p) $ counts how many other centroids fall within 5km of point (p). Cap is 10: if an islet has more than 10 neighbors that close, it’s sitting in a dense reef field, a crowded coastline or a archipelago clutter, not a small isolated 3-4 island group like the photo shows. This dropped candidates down to 51,576. $ N_{5\text{km}}(p) $ 用于计算点 (p) 周围 5 公里内有多少个其他质心。上限设为 10:如果一个小岛在如此近的范围内有超过 10 个邻居,说明它位于密集的礁石区、拥挤的海岸线或群岛中,而不是照片中那种 3-4 个岛屿的小型孤立群。这使候选数量降至 51,576 个。

03] Clustering 03] 聚类 For every surviving point, find every other point within 20km (heuristic, by eye from the image). If it has at least 2 neighbors that close (3 points total), it’s a cluster. Points with no cluster of 3+ nearby are dropped, they can’t form a triangle at all. 对于每一个幸存的点,找出 20 公里内(基于图像的肉眼直觉)的所有其他点。如果它至少有 2 个邻居(总共 3 个点),它就是一个聚类。没有 3 个以上点组成的聚类的点会被剔除,因为它们无法构成三角形。

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 $$ That collapses down to 23,500 clusters. 这使聚类数量缩减至 23,500 个。

04] Generating Triplets 04] 生成三元组 For every cluster, every combination of 3 points inside it becomes a candidate triangle. That’s $ C(n, 3) $, which explodes fast for big clusters, for example: a cluster of 60 points already gives 34,220 triples on its own. So each cluster gets capped at 60 points first, sampled by size, not randomly. 对于每个聚类,其中任意 3 个点的组合都成为候选三角形。即 $ C(n, 3) $,对于大型聚类,这个数字会迅速膨胀。例如:一个包含 60 个点的聚类本身就能产生 34,220 个三元组。因此,每个聚类首先被限制在 60 个点以内,并按大小而非随机进行采样。

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]

The sampling takes a third small islands, a third large, a third from the middle of the size distribution, instead of the full cluster or a random cut. 23,500 clusters produce 80,690,777 triples total !! 采样过程从大小分布中分别取三分之一的小岛、三分之一的大岛和三分之一的中等岛屿,而不是使用整个聚类或随机截取。23,500 个聚类总共产生了 80,690,777 个三元组!!

05] Matching, on the GPU 05] 在 GPU 上进行匹配 I gave every triple one CUDA thread. Each thread sorts its 3 points by land area to pick out P0 (smallest, the resort islet), then uses the winding direction of the other two to assign P1 and P2: 我为每个三元组分配了一个 CUDA 线程。每个线程按陆地面积对 3 个点进行排序,以选出 P0(最小的,即度假村小岛),然后利用另外两个点的绕行方向来指定 P1 和 P2:

long long i = blockIdx.x * (long long)blockDim.x + threadIdx.x;
if (i >= n_triples) return;
int pos[3] = {0, 1, 2};
// ... (sorting logic)

P1 vs P2 comes from a 2D cross product, no branching on which cluster the triple came from, just the sign: P1 和 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} $$ Walk from P0 to a, then to b. If cross > 0, that’s a left turn (counterclockwise). If cross < 0, it’s a right turn (clockwise). It’s the same sign trick used to tell if 3 points curve one way or the other. then angle at P0 and the distance ratio, same formulas as the fingerprint step, computed independently by every thread: 从 P0 走到 a,再走到 b。如果叉积 > 0,则是左转(逆时针)。如果叉积 < 0,则是右转(顺时针)。这与判断 3 个点是向哪一侧弯曲的符号技巧相同。然后计算 P0 处的角度和距离比,公式与指纹步骤相同,由每个线程独立计算:

$$ \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}|} $$ A triple survives if angle, ratio, P0’s size, the separation between P0 and P1, and both side lengths all land inside the fingerprint’s tolerance windows. Threads that pass write their result into a shared output array using an atomic counter, so two threads finishing at the same time never overwrite each other: 如果角度、比例、P0 的大小、P0 与 P1 之间的距离以及两条边长都落在指纹的容差范围内,则该三元组通过筛选。通过筛选的线程使用原子计数器将结果写入共享输出数组,这样两个同时完成的线程就不会相互覆盖:

if (hit) { unsigned ...