用几何与CUDA定位随机岛屿：完整技术复盘
用几何与CUDA定位随机岛屿
做地理定位或空间计算的同学必看，这套从指纹构建到GPU并行匹配的完整流程可直接复用到类似任务，代码和参数都开源了，赶紧拿去跑一遍。
gralhix004 | Geolocating Random Islet Image Using Geometry & CUDA GPU Programming
gralhix004 | 使用几何与CUDA GPU编程定位随机小岛图像
16-08-2026
16-08-2026
NOTE: this is a genuine human work, didnt use LLM generation.
注意：这是真实的人类作品，未使用LLM生成。
I'm writing this page as a writeup for this challenge gralhix 004 made by Sofia Santos | Gralhix.
我写这个页面是为了记录由Sofia Santos | Gralhix制作的gralhix 004挑战的解题过程。
You can view, clone and locally try all code files and the final report with all instructions here at github.
你可以在github上查看、克隆并本地尝试所有代码文件和包含所有说明的最终报告。
Task briefing:
任务简报：
This is a photo of a resort located on an island.
这是一张位于岛屿上的度假村的照片。
a) What is the name of the resort?
a) 度假村的名称是什么？
b) What are the coordinates of the island?
b) 岛屿的坐标是什么？
c) In which cardinal direction was the camera facing when the photo was taken?
c) 拍摄照片时相机朝向哪个基本方向？
In my opinion, solving this challenge with google lens is wasting a fun opportunity, so decided to solve it with math and programming.
在我看来，用谷歌镜头解决这个挑战会浪费一个有趣的机会，所以我决定用数学和编程来解决它。
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> exiftool main.png
File Type : WEBP (lossless)
MIME Type : image/webp
Image Width : 736
Image Height : 515As 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:
从图片中可以看到，有3块陆地：
- P0: the islet itself,
- P1: the right island,
- P2: the left front island ( having mountain peak )
- 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.
所以我不得不凭直觉估计，我只想要3个岛屿之间的相对距离和那个三角形的角度。
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.
我构建了一个小的点击GUI 01_triangle_gui.py，它按顺序记录每个点的像素坐标并计算三角形的几何形状。
Since clicking exact centers by eye isn't perfectly precise, I added a ±20% tolerance band around both values when searching.
由于用眼睛点击精确中心并不完美精确，我在搜索时对两个值都添加了±20%的容差范围。
c] SEARCH
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.
我使用了OpenStreetMap的拆分陆地多边形数据集land-polygons-split-4326，这是WGS84坐标系下的全球海岸线矢量数据，大小为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.
我创建了启发式过滤器（全部仅凭直觉和非具体证据），花了数天时间（是的，整天）调整数值，并进行了大量试错😭，直到我得到了这个有效的过滤器配方。
01] Tropical latitude bounding box
01] 热带纬度边界框
$$ -30° \le latitude \le 30° $$
$$ -30° \le 纬度 \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) \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.
$ N_{5\text{km}}(p) $ 计算在点(p)的5公里范围内有多少其他质心。上限为10：如果一个小岛有超过10个如此近的邻居，那么它位于密集的礁石区、拥挤的海岸线或群岛杂乱区，而不是像照片中那样的小型孤立3-4岛群。
This dropped candidates down to 51,576.
这将候选数量减少到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)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 $$
$$ \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个点，按大小采样，而不是随机。
$$ \binom{n}{3} = \frac{n(n-1)(n-2)}{6} $$
$$ \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]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};
for (int a1 = 1; a1 < 3; a1++)
{
int key = pos[a1];
double keyval = a[key];
int j = a1 - 1;
while (j >= 0 && a[pos[j]] > keyval)
{
pos[j + 1] = pos[j];
j--;
}
pos[j + 1] = key;
}long long i = blockIdx.x * (long long)blockDim.x + threadIdx.x;
if (i >= n_triples) return;
int pos[3] = {0, 1, 2};
for (int a1 = 1; a1 < 3; a1++)
{
int key = pos[a1];
double keyval = a[key];
int j = a1 - 1;
while (j >= 0 && a[pos[j]] > keyval)
{
pos[j + 1] = pos[j];
j--;
}
pos[j + 1] = key;
}P1 vs P2 comes from a 2D cross product, no branching on which cluster the triple came from, just the sign:
P1与P2来自2D叉积，不分支于三元组来自哪个聚类，只看符号：
$$ \text{cross} = x_a y_b - x_b y_a $$ $$ P1 = \begin{cases} a & \text{cross} > 0 \\ b & \text{cross} \le 0 \end{cases} $$
$$ \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.
从P0走到a，再到b。如果叉积>0，那是左转（逆时针）。如果叉积<0，那是右转（顺时针）。这和用来判断3个点是否朝某一方向弯曲的符号技巧相同。
then angle at P0 and the distance ratio, same formulas as the fingerprint step, computed independently by every thread:
然后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}|} $$
$$ \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 long long slot = atomicAdd(out_count, 1ULL);
out_p0[slot] = p0idx;
out_p1[slot] = p1idx;
out_p2[slot] = p2idx;
}if (hit)
{
unsigned long long slot = atomicAdd(out_count, 1ULL);
out_p0[slot] = p0idx;
out_p1[slot] = p1idx;
out_p2[slot] = p2idx;
}Now printed in the CLI directly from the kernel:
现在直接从内核在CLI中打印：
gpu: NVIDIA GeForce RTX 3050 (sm_86)
vram used: 5169 MB
kernel time: 204.1 msgpu: NVIDIA GeForce RTX 3050 (sm_86)
vram used: 5169 MB
kernel time: 204.1 ms80.7 million triples go in, one thread each, in parallel. 158,784 pass the mask.
8070万个三元组进入，每个线程一个，并行处理。158,784个通过掩码。
06] Dedup
06] 去重
Since same physical triple can get hit by multiple GPU threads if it belonged to more than one overlapping cluster, so raw matches get collapsed by identity first:
由于同一个物理三元组如果属于多个重叠的簇，可能会被多个GPU线程命中，因此原始匹配首先按身份进行折叠：
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)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)8,915 unique triples after dedup.
去重后得到8,915个唯一三元组。
07] The Open Rectangle
07] 开放矩形
Every surviving triple gets one more test: is the space next to it actually open water, like the photo shows ? A rectangle gets built along the P0→P1 edge, on whichever side P2 is not on, then checked against the land dataset for anything else sitting inside it.
每个幸存的三元组还要经过一次测试：它旁边的空间是否实际上是开阔水域，就像照片显示的那样？沿着P0→P1边构建一个矩形，位于P2不在的那一侧，然后对照陆地数据集检查是否有其他东西位于其中。
width = np.hypot(x1, y1)
u = np.array([x1, y1]) / width
v = np.array([-u[1], u[0]])
# p2 sits on the +v side by construction,
# so the check goes on -v
length = 2 * width
corners_local = [
(0, 0), (x1, y1),
(x1 - v[0]*length, y1 - v[1]*length),
(-v[0]*length, -v[1]*length),
]width = np.hypot(x1, y1)
u = np.array([x1, y1]) / width
v = np.array([-u[1], u[0]])
# p2 sits on the +v side by construction,
# so the check goes on -v
length = 2 * width
corners_local = [
(0, 0), (x1, y1),
(x1 - v[0]*length, y1 - v[1]*length),
(-v[0]*length, -v[1]*length),
]If anything other than the 3 candidate islands themselves intersects that rectangle, the candidate is dropped. Land sitting there means it's not the open, unobstructed water the photo actually shows.
如果除了3个候选岛屿本身之外的任何东西与该矩形相交，则丢弃该候选。那里有陆地意味着它不是照片实际显示的开放、无障碍的水域。
8,915 unique triples down to 948.
8,915个唯一三元组减少到948个。
更进一步：量化金融体系
看懂新闻只是起点——沿量化金融路径，把它变成能交付的工程能力