Scattered in the fog (TSG CTF 2024) WriteUp

# はじめに
この記事は [1日1CTF Advent Calendar 2024](https://qiita.com/advent-calendar/2024/1-day-1-ctf) の 14 日目の記事です。

# 問題
Scattered in the fog (問題出典: TSG CTF 2024)
> 覆水盆に返らず。
> フラグの形式は TSGCTF{[A-Z_]+} です

## 問題概要

```server.py
import cv2
import numpy as np
import open3d as o3d
from scipy.stats import special_ortho_group, norm
import string

alphabet = string.ascii_uppercase + "{}_"
print(len(alphabet))

patchsize = 64
shift = 64
patches = np.zeros((len(alphabet)+1, patchsize, patchsize), np.uint8)

offset = 8
for i in range(len(alphabet)):
    cv2.putText(patches[i], alphabet[i], (offset, patchsize-offset), cv2.FONT_HERSHEY_SIMPLEX, 2.0, 255, 4)
cv2.imwrite("sample.png", patches.reshape(5, 6, patchsize, patchsize).transpose(0,2,1,3).reshape(patchsize * 5, patchsize * 6))

flag = "TSGCTF{_______________________________________________________}"  # secret!
assert len(flag) == 63
radii = (len(flag) - 1) / 2

coords = []
for i in range(len(flag)):
    index = alphabet.index(flag[i])
    img = patches[index]
    xmap, ymap = np.meshgrid(np.arange(patchsize), np.arange(patchsize))
    xs = xmap[np.where(img == 255)].astype(np.float32)[:, None]
    ys = ymap[np.where(img == 255)].astype(np.float32)[:, None]
    zs = np.full_like(xs, shift * (i - radii))

    s = 16 * norm.pdf((i - radii), loc=0, scale=11) / norm.pdf(0, loc=0, scale=11)
    noise_xs = np.round(np.random.normal(0, scale=s, size=xs.shape)) * shift
    noise_ys = np.round(np.random.normal(0, scale=s, size=ys.shape)) * shift

    coords.append(np.concatenate([xs + noise_xs, ys + noise_ys, zs], axis=1))

coords = np.concatenate(coords, axis=0)

# random rotate
rot = special_ortho_group.rvs(3)
coords = coords @ rot

# random translate
trl = np.random.normal(0, scale=100, size=3)
coords += trl

# shuffle order
coords = coords[np.random.permutation(len(coords))]
np.save("problem.npy", coords)

# for visualization
# pcd = o3d.t.geometry.PointCloud(coords)
# o3d.visualization.draw_geometries([pcd.to_legacy()])
```

## 考察

とりあえず `for visualization` のコメントアウトを消してビジュアライズしてみる。コマみたいな図形が観測できた。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/942569/de1d1636-a78d-74e8-3039-e9ba48b2b903.png)

もう少し見てみると、尖った方から見ると「TSGCTF」の文字が確認できた。ただ、中央部はめちゃくちゃに散らばっているので、そのまま読むのは難しそう。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/942569/d5872073-7ba2-5f62-0140-1c506b962db3.png)

ちゃんとコードを読む。

各文字について次の処理を行っている。
```py
for i in range(len(flag)):
    index = alphabet.index(flag[i])
    img = patches[index]
    xmap, ymap = np.meshgrid(np.arange(patchsize), np.arange(patchsize))
    xs = xmap[np.where(img == 255)].astype(np.float32)[:, None]
    ys = ymap[np.where(img == 255)].astype(np.float32)[:, None]
    zs = np.full_like(xs, shift * (i - radii))

    s = 16 * norm.pdf((i - radii), loc=0, scale=11) / norm.pdf(0, loc=0, scale=11)
    noise_xs = np.round(np.random.normal(0, scale=s, size=xs.shape)) * shift
    noise_ys = np.round(np.random.normal(0, scale=s, size=ys.shape)) * shift

    coords.append(np.concatenate([xs + noise_xs, ys + noise_ys, zs], axis=1))
```
各文字の画像から真っ白な格子点を取り出して、フラグの真ん中に行くほど分散が大きくなるような正規分布(から取った値を四捨五入して `shift` 倍したもの)のノイズを加えている。

```py
coords = np.concatenate(coords, axis=0)

# random rotate
rot = special_ortho_group.rvs(3)
coords = coords @ rot

# random translate
trl = np.random.normal(0, scale=100, size=3)
coords += trl

# shuffle order
coords = coords[np.random.permutation(len(coords))]
np.save("problem.npy", coords)
```

その後、適当に回転/並行移動してから順番をぐちゃぐちゃにしている。

まず 回転/並行 移動をもとに戻す方法を考えよう。

フラグの最初/最後の `T` / `}` に対応する点はほとんどノイズが乗らない。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/942569/65bd3eb8-3ee0-ba9a-307f-430c501fd912.png)

よって、この中から 4 つの点を適当に選んで、それに対応する点を頑張って当てればいい。この 4 つの点は端の方にあるため、変換後の座標でも端の方にあることを期待 + それぞれの点の間の距離は変換後の座標でも変わらないことを利用して枝狩りすれば、十分高速。あとは正しく復元できるような 4 点を最初に選べていたことをお祈りする。ダメだったら変えてやり直す。

```py
import cv2
import numpy as np
import open3d as o3d
from scipy.stats import special_ortho_group, norm
import string

alphabet = string.ascii_uppercase + "{}_"
print(len(alphabet))

patchsize = 64
shift = 64
patches = np.zeros((len(alphabet) + 1, patchsize, patchsize), np.uint8)

offset = 8
for i in range(len(alphabet)):
    cv2.putText(
        patches[i],
        alphabet[i],
        (offset, patchsize - offset),
        cv2.FONT_HERSHEY_SIMPLEX,
        2.0,
        255,
        4,
    )
cv2.imwrite(
    "sample.png",
    patches.reshape(5, 6, patchsize, patchsize)
    .transpose(0, 2, 1, 3)
    .reshape(patchsize * 5, patchsize * 6),
)

flag = "TSGCTF{_______________________________________________________}"  # secret!
assert len(flag) == 63
radii = (len(flag) - 1) / 2

coords = []
for i in range(len(flag)):
    index = alphabet.index(flag[i])
    img = patches[index]
    xmap, ymap = np.meshgrid(np.arange(patchsize), np.arange(patchsize))
    xs = xmap[np.where(img == 255)].astype(np.float32)[:, None]
    ys = ymap[np.where(img == 255)].astype(np.float32)[:, None]
    zs = np.full_like(xs, shift * (i - radii))

    coords.append(np.concatenate([xs, ys, zs], axis=1))

coords = np.concatenate(coords, axis=0)

p0, p1, p2, p3 = coords[4], coords[-4], coords[30], coords[-2] # お祈りパート

coords = np.load("problem.npy")
dis_p0_p1 = (p0[0] - p1[0]) ** 2 + (p0[1] - p1[1]) ** 2 + (p0[2] - p1[2]) ** 2
dis_p0_p2 = (p0[0] - p2[0]) ** 2 + (p0[1] - p2[1]) ** 2 + (p0[2] - p2[2]) ** 2
dis_p0_p3 = (p0[0] - p3[0]) ** 2 + (p0[1] - p3[1]) ** 2 + (p0[2] - p3[2]) ** 2
dis_p1_p2 = (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2
dis_p1_p3 = (p1[0] - p3[0]) ** 2 + (p1[1] - p3[1]) ** 2 + (p1[2] - p3[2]) ** 2
dis_p2_p3 = (p2[0] - p3[0]) ** 2 + (p2[1] - p3[1]) ** 2 + (p2[2] - p3[2]) ** 2

coords = sorted(coords, key=lambda x: x[0])
coords = np.array(coords)
ok = False
# p は全て端っこの方の点であることを利用して枝狩り
for q0 in list(range(600)) + list(range(len(coords) - 600, len(coords))):
    for q1 in list(range(600)) + list(range(len(coords) - 600, len(coords))):
        dis_q0_q1 = (
            (coords[q0][0] - coords[q1][0]) ** 2
            + (coords[q0][1] - coords[q1][1]) ** 2
            + (coords[q0][2] - coords[q1][2]) ** 2
        )
        if abs(dis_p0_p1 - dis_q0_q1) < 1e-8:
            for q2 in list(range(600)) + list(range(len(coords) - 600, len(coords))):
                dis_q0_q2 = (
                    (coords[q0][0] - coords[q2][0]) ** 2
                    + (coords[q0][1] - coords[q2][1]) ** 2
                    + (coords[q0][2] - coords[q2][2]) ** 2
                )
                dis_q1_q2 = (
                    (coords[q1][0] - coords[q2][0]) ** 2
                    + (coords[q1][1] - coords[q2][1]) ** 2
                    + (coords[q1][2] - coords[q2][2]) ** 2
                )
                if (
                    abs(dis_p0_p2 - dis_q0_q2) < 1e-8
                    and abs(dis_p1_p2 - dis_q1_q2) < 1e-8
                ):
                    for q3 in list(range(600)) + list(
                        range(len(coords) - 600, len(coords))
                    ):
                        dis_q0_q3 = (
                            (coords[q0][0] - coords[q3][0]) ** 2
                            + (coords[q0][1] - coords[q3][1]) ** 2
                            + (coords[q0][2] - coords[q3][2]) ** 2
                        )
                        dis_q1_q3 = (
                            (coords[q1][0] - coords[q3][0]) ** 2
                            + (coords[q1][1] - coords[q3][1]) ** 2
                            + (coords[q1][2] - coords[q3][2]) ** 2
                        )
                        dis_q2_q3 = (
                            (coords[q2][0] - coords[q3][0]) ** 2
                            + (coords[q2][1] - coords[q3][1]) ** 2
                            + (coords[q2][2] - coords[q3][2]) ** 2
                        )
                        if (
                            abs(dis_p0_p3 - dis_q0_q3) < 1e-8
                            and abs(dis_p1_p3 - dis_q1_q3) < 1e-8
                            and abs(dis_p2_p3 - dis_q2_q3) < 1e-8
                        ):
                            print(q0, q1, q2, q3)
                            ok = True
                            break
                if ok:
                    break
        if ok:
            break
    if ok:
        break
        

import sympy

x0, y0, z0 = sympy.symbols("x0 y0 z0")
rxx, rxy, rxz = sympy.symbols("rxx rxy rxz")
ryx, ryy, ryz = sympy.symbols("ryx ryy ryz")
rzx, rzy, rzz = sympy.symbols("rzx rzy rzz")
s = []
print(q0, q1, q2, q3)
for p,q in zip([p0, p1, p2, p3], [coords[q0], coords[q1], coords[q2], coords[q3]]):
    print(p, q)
    x = q[0] * rxx + q[1] * ryx + q[2] * rzx + x0
    y = q[0] * rxy + q[1] * ryy + q[2] * rzy + y0
    z = q[0] * rxz + q[1] * ryz + q[2] * rzz + z0
    s.append(x - p[0])
    s.append(y - p[1])
    s.append(z - p[2])

sol = sympy.solve(s, [x0, y0, z0, rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz])

x0, y0, z0 = sol[x0], sol[y0], sol[z0]
rxx, rxy, rxz = sol[rxx], sol[rxy], sol[rxz]
ryx, ryy, ryz = sol[ryx], sol[ryy], sol[ryz]
rzx, rzy, rzz = sol[rzx], sol[rzy], sol[rzz]

rot_rev = np.array([[rxx, rxy, rxz], [ryx, ryy, ryz], [rzx, rzy, rzz]], dtype=np.float64)
trl_rev = np.array([x0, y0, z0], dtype=np.float64)
new_coords = coords @ rot_rev + trl_rev

# for visualization
pcd = o3d.t.geometry.PointCloud(new_coords)
o3d.visualization.draw_geometries([pcd.to_legacy()])
```

いい感じに復元できてそう。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/942569/0fbdaa48-2572-2413-9377-b6694ef2d95a.png)

あとは各文字について、全ての点の z 座標は同じであることを利用。
各 z 座標の点の個数から、もとの文字を予測する。

```python
import cv2
import numpy as np
import open3d as o3d
from scipy.stats import special_ortho_group, norm
import string

alphabet = string.ascii_uppercase + "{}_"
print(len(alphabet))

patchsize = 64
shift = 64
patches = np.zeros((len(alphabet) + 1, patchsize, patchsize), np.uint8)

offset = 8
for i in range(len(alphabet)):
    cv2.putText(
        patches[i],
        alphabet[i],
        (offset, patchsize - offset),
        cv2.FONT_HERSHEY_SIMPLEX,
        2.0,
        255,
        4,
    )
cv2.imwrite(
    "sample.png",
    patches.reshape(5, 6, patchsize, patchsize)
    .transpose(0, 2, 1, 3)
    .reshape(patchsize * 5, patchsize * 6),
)

flag = "TSGCTF{_______________________________________________________}"  # secret!
assert len(flag) == 63
radii = (len(flag) - 1) / 2

coords = []
for i in range(len(flag)):
    index = alphabet.index(flag[i])
    img = patches[index]
    xmap, ymap = np.meshgrid(np.arange(patchsize), np.arange(patchsize))
    xs = xmap[np.where(img == 255)].astype(np.float32)[:, None]
    ys = ymap[np.where(img == 255)].astype(np.float32)[:, None]
    zs = np.full_like(xs, shift * (i - radii))
    coords.append(np.concatenate([xs, ys, zs], axis=1))

coords = np.concatenate(coords, axis=0)

p0, p1, p2, p3 = coords[4], coords[-4], coords[30], coords[-2]  # お祈りパート

coords = np.load("problem.npy")
dis_p0_p1 = (p0[0] - p1[0]) ** 2 + (p0[1] - p1[1]) ** 2 + (p0[2] - p1[2]) ** 2
dis_p0_p2 = (p0[0] - p2[0]) ** 2 + (p0[1] - p2[1]) ** 2 + (p0[2] - p2[2]) ** 2
dis_p0_p3 = (p0[0] - p3[0]) ** 2 + (p0[1] - p3[1]) ** 2 + (p0[2] - p3[2]) ** 2
dis_p1_p2 = (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2
dis_p1_p3 = (p1[0] - p3[0]) ** 2 + (p1[1] - p3[1]) ** 2 + (p1[2] - p3[2]) ** 2
dis_p2_p3 = (p2[0] - p3[0]) ** 2 + (p2[1] - p3[1]) ** 2 + (p2[2] - p3[2]) ** 2

coords = sorted(coords, key=lambda x: x[0])
coords = np.array(coords)
ok = False
# p は全て端っこの方の点であることを利用して枝狩り
for q0 in list(range(600)) + list(range(len(coords) - 600, len(coords))):
    for q1 in list(range(600)) + list(range(len(coords) - 600, len(coords))):
        dis_q0_q1 = (
            (coords[q0][0] - coords[q1][0]) ** 2
            + (coords[q0][1] - coords[q1][1]) ** 2
            + (coords[q0][2] - coords[q1][2]) ** 2
        )
        if abs(dis_p0_p1 - dis_q0_q1) < 1e-8:
            for q2 in list(range(600)) + list(range(len(coords) - 600, len(coords))):
                dis_q0_q2 = (
                    (coords[q0][0] - coords[q2][0]) ** 2
                    + (coords[q0][1] - coords[q2][1]) ** 2
                    + (coords[q0][2] - coords[q2][2]) ** 2
                )
                dis_q1_q2 = (
                    (coords[q1][0] - coords[q2][0]) ** 2
                    + (coords[q1][1] - coords[q2][1]) ** 2
                    + (coords[q1][2] - coords[q2][2]) ** 2
                )
                if (
                    abs(dis_p0_p2 - dis_q0_q2) < 1e-8
                    and abs(dis_p1_p2 - dis_q1_q2) < 1e-8
                ):
                    for q3 in list(range(600)) + list(
                        range(len(coords) - 600, len(coords))
                    ):
                        dis_q0_q3 = (
                            (coords[q0][0] - coords[q3][0]) ** 2
                            + (coords[q0][1] - coords[q3][1]) ** 2
                            + (coords[q0][2] - coords[q3][2]) ** 2
                        )
                        dis_q1_q3 = (
                            (coords[q1][0] - coords[q3][0]) ** 2
                            + (coords[q1][1] - coords[q3][1]) ** 2
                            + (coords[q1][2] - coords[q3][2]) ** 2
                        )
                        dis_q2_q3 = (
                            (coords[q2][0] - coords[q3][0]) ** 2
                            + (coords[q2][1] - coords[q3][1]) ** 2
                            + (coords[q2][2] - coords[q3][2]) ** 2
                        )
                        if (
                            abs(dis_p0_p3 - dis_q0_q3) < 1e-8
                            and abs(dis_p1_p3 - dis_q1_q3) < 1e-8
                            and abs(dis_p2_p3 - dis_q2_q3) < 1e-8
                        ):
                            print(q0, q1, q2, q3)
                            ok = True
                            break
                if ok:
                    break
        if ok:
            break
    if ok:
        break


import sympy

x0, y0, z0 = sympy.symbols("x0 y0 z0")
rxx, rxy, rxz = sympy.symbols("rxx rxy rxz")
ryx, ryy, ryz = sympy.symbols("ryx ryy ryz")
rzx, rzy, rzz = sympy.symbols("rzx rzy rzz")
s = []
print(q0, q1, q2, q3)
for p, q in zip([p0, p1, p2, p3], [coords[q0], coords[q1], coords[q2], coords[q3]]):
    print(p, q)
    x = q[0] * rxx + q[1] * ryx + q[2] * rzx + x0
    y = q[0] * rxy + q[1] * ryy + q[2] * rzy + y0
    z = q[0] * rxz + q[1] * ryz + q[2] * rzz + z0
    s.append(x - p[0])
    s.append(y - p[1])
    s.append(z - p[2])

sol = sympy.solve(s, [x0, y0, z0, rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz])

x0, y0, z0 = sol[x0], sol[y0], sol[z0]
rxx, rxy, rxz = sol[rxx], sol[rxy], sol[rxz]
ryx, ryy, ryz = sol[ryx], sol[ryy], sol[ryz]
rzx, rzy, rzz = sol[rzx], sol[rzy], sol[rzz]

rot_rev = np.array(
    [[rxx, rxy, rxz], [ryx, ryy, ryz], [rzx, rzy, rzz]], dtype=np.float64
)
trl_rev = np.array([x0, y0, z0], dtype=np.float64)
new_coords = coords @ rot_rev + trl_rev

dic = {}

for i in range(len(new_coords)):
    z = round(new_coords[i][2])
    if z not in dic:
        dic[z] = 1
    else:
        dic[z] += 1

print(dic)

dic2 = {}
for i in range(len(alphabet)):
    img = patches[i]
    xmap, ymap = np.meshgrid(np.arange(patchsize), np.arange(patchsize))
    xs = xmap[np.where(img == 255)].astype(np.float32)[:, None]
    if len(xs) in dic2:
        dic2[len(xs)] += alphabet[i]
    else:
        dic2[len(xs)] = alphabet[i]
print(dic2)

ans = []
for k, v in dic.items():
    ans.append((k, dic2[v]))
ans = sorted(ans, key=lambda x: x[0])
print(ans)

flag = ""
for a in ans:
    if len(a[1]) == 1:
        flag += a[1]
    else:
        flag += "[" + a[1] + "]"
print(flag)

# for visualization
# pcd = o3d.t.geometry.PointCloud(new_coords)
# o3d.visualization.draw_geometries([pcd.to_legacy()])
```

いくつか候補があるときは、`[]` で囲うようにした。

```:result
T[AS]GCTF[Y{}][AS][AS]K_TH[AS]T_DEMON_TO_[AS]IMUL[AS]TE_THE_REVER[AS]E_OF_[AS]PILL_PROCE[AS][AS][Y{}]
```

英文っぽいので意味が通るように文字を選ぶ。

flag: `TSGCTF{ASK_THAT_DEMON_TO_SIMULATE_THE_REVERSE_OF_SPILL_PROCESS}`