幾何学とCUDAプログラミングでランダムな島を地理位置特定する

gralhix004 | 幾何学とCUDA GPUプログラミングによるランダムな小島の画像地理位置特定 16-08-2026 注: これは本物の人間による作品であり、LLM生成は使用していません。 このページは、Sofia Santos | Gralhix 作成のチャレンジ gralhix 004 の writeup として書いています。すべてのコードファイルと最終レポートは、指示とともにgithubで表示、クローン、ローカルで試すことができます。 タスク概要: これは、島にあるリゾート地の写真です。 a) リゾートの名前は何ですか？ b) その島の座標は何ですか？ c) 写真が撮影されたとき、カメラはどの主要な方向を向いていましたか？ 私の意見では、Googleレンズでこのチャレンジを解決するのは楽しい機会を無駄にすることなので、数学とプログラミングで解決することにしました。 a] メタデータ もちろん、最初に調べるのはメタデータです。 Linux void で実行しました: > exiftool main.png File Type : WEBP (lossless) MIME Type : image/webp Image Width : 736 Image Height : 515 予想通り、有用なものはありませんでした。EXIF、GPS、カメラのメーカーやモデルもありません。 b] フィンガープリントの構築 画像から、3つの陸地が見えます: P0: 小島そのもの P1: 右側の島 P2: 左手前の島（山の頂上がある） この画像は明らかにドローンで撮影されたもので、標高を全く推定できないため、鳥瞰図の正しい透視モデルを構築できませんでした（メタデータにも見つかりませんでした）。そのため、直感で推定する必要がありました。3つの島間の相対的な距離と、その三角形の角度だけが欲しかったのです。 ピクセル座標を順番に記録し、三角形の幾何学を計算する小さなクリックGUI 01_triangle_gui.py を作成しました。 目視で正確な中心をクリックするのは完全に正確ではないため、検索時に両方の値の周りに±20%の許容範囲を追加しました。 c] 検索 フィンガープリントが確定したら、次のステップは地球上のすべての実際の陸地をそれと比較することです！ OpenStreetMap の分割された陸地ポリゴンセットをデータセット land-polygons-split-4326 (WGS84 の完全なグローバル海岸線ベクトル、サイズは 882 MB) として使用しました。 直感と非実証的な証拠のみに基づいて、ヒューリスティックなフィルターを作成しました。値の調整と数え切れないほどの試行錯誤に何日も費やしました😭、ようやくこれが機能するフィルターレシピになりました。 01] 熱帯緯度バウンディングボックス $$ -30° \le latitude \le 30° $$ 写真の小島は熱帯のように見えるので、計算コストの高い幾何学的な作業を行う前に、熱帯の外にあるものはすべてすぐに除外することにしました。このバンドフィルターを通過したのは正確に 141,131 の陸地ポリゴンです。 02] 局所密度フィルター $$ N_{5\text{km}}(p) \le 10 $$ $ N_{5\text{km}}(p) $ は、点 (p) から 5 km 以内にある他の重心の数を数えます。上限は 10 です。小島が近くに 10 を超える隣接点を持つ場合、それは密集したサンゴ礁地帯、混雑した海岸線、または群島の雑然とした場所にあることを意味し、写真が示すような小さく孤立した3〜4島のグループではありません。これにより候補は 51,576 に減りました。 03] クラスタリング 生き残った各点について、20 km 以内にある他のすべての点を見つけます（ヒューリスティック、画像から目視）。少なくとも 2 つの近い隣接点（合計 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 $$ これは 23,500 のクラスターにまで減少します。 04] 三つ組の生成 各クラスターについて、その中の 3 点のすべての組み合わせが候補となる三角形になります。これは $ C(n, 3) $ で、大きなクラスターでは急速に爆発します。例えば、60 点のクラスターだけでも 34,220 の三つ組が生成されます。そのため、各クラスターはまず 60 点に制限され、ランダムではなくサイズでサンプリングされます。 $$ \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] サンプリングは、クラスター全体やランダムな切り取りではなく、小さい島を 1/3、大きい島を 1/3、サイズ分布の中央から 1/3 を取得します。 23,500 のクラスターから合計 80,690,777 の三つ組が生成されます！！ 05] マッチング、GPU上 各三つ組に 1 つの CUDA スレッドを割り当てました。各スレッドは 3 点を陸地面積でソートして P0（最小、リゾートの小島）を選択し、他の 2 つの向きから 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; } 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} $$ P0 から a、次に b へ歩きます。cross > 0 なら左折（反時計回り）です。cross < 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}|} $$ 角度、比率、P0 のサイズ、P0 と P1 の間の分離、および両側の長さがすべてフィンガープリントの許容範囲内に収まる場合、その三つ組は生き残ります。合格したスレッドは、アトミックカウンターを使用して結果を共有出力配列に書き込みます。これにより、同時に終了した 2 つのスレッドが互いを上書きすることはありません。 if (hit) { unsigned long long slot = atomicAdd(out_count, 1ULL); out_p0[slot] = p0idx; out_p1[slot] = p1idx; out_p2[slot] = p2idx; } カーネルから直接 CLI に表示されます: gpu: NVIDIA GeForce RTX 3050 (sm_86) vram used: 5169 MB kernel time: 204.1 ms 8070万の三つ組が並列で 1 スレッドずつ入力され、158,784 がマスクを通過しました。 06] 重複排除 同じ物理的な三つ組が、複数の重複するクラスターに属していたために複数の 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) 重複排除後、8,915 のユニークな三つ組。 07] 開いた長方形 生き残った各三つ組は、もう 1 つのテストを受けます。その隣の空間は、写真が示すように実際に開けた水域であるか？ P0→P1 の辺に沿って、P2 がない方の側に長方形が構築され、その中に他のものが座っていないか陸地データセットに対してチェックされます。 width = np.hypot(x1, y1) u = np.array([x1, y1]) / width v = np.array([-u[1], u[0]]) # p2 は構築により +v 側にあります。 # そのため、チェックは -v で行われます。 length = 2 * width corners_local = [ (0, 0), (x1, y1), (x1 - v[0]*length, y1 - v[1]*length), (-v[0]*length, -v[1]*length), ] 3 つの候補となる島自体以外のものがその長方形に交差する場合、候補はドロップされます。そこに陸地があるということは、写真が実際に示している、遮るもののない開けた水域ではないということです。 8,915 のユニークな三つ組が 948 に減少しました。 そして以下は 948 の候補地の地図です。 d] サンゴ礁の形状チェック この段階では、P0、つまりリゾートの小島のみに注目し、その形状が実際にサンゴ礁の形状のように見えるかどうかをチェックします。 1] コンパクトさ、形状が円にどれだけ近いか: Polsby-Popper スコア: $$ PP = \frac{4\pi \cdot \text{area}}{\text{perimeter}^2} $$ def compactness(row): return (4 * np.pi * row.area_km2) / (row.perim_km ** 2 + 1e-12) 1.0 は完全な円であり、低いほどギザギザまたは細長い輪郭になります。サンゴ礁は波の堆積により円形になる傾向があります。