title: '数学の絵は、きれいなので誰も確かめない ―― 定理を門にして描く PoC シリーズ'
tags:

• Python
• NumPy
• 数学
• アルゴリズム
• 画像処理
  public_private: true
  public_id: fe04f6eef40119913894

言語 / Language: 日本語 · English

数学の絵は、きれいなので誰も確かめない ―― 定理を門にして描く PoC シリーズ

↑ 腕が 1 → 4 → 16 → 90 → 600 → 4000 本と増えるにつれ絵が正体を現す。どの段でどこまで似るかは描く前に分かる。

図について: この記事が扱う 7 回ぶんの図は、下に全部出してあります(86 枚)。シリーズ外も含めた全 151 回・1165 枚は紙面の計測館に並べてあり、元のファイルはこちらにあります。

この記事について

位相彩色、ニュートンの吸引域、マンデルブロ集合、アポロニウスの円、葉序の螺旋、翼まわりの流れ。数学の図はきれいです。そしてきれいなものは、合っているかを誰も確かめません。

実装が少し間違っていても、出てくる絵はやはりきれいです。色は回るし、境界はフラクタルになるし、円は詰まります。間違いが絵に出ない——これがこの分野を書くときの難所です。

そこで描画コードを足すときの基準を 2 つに絞りました。

(i) 定理が門になるか、(ii) 既存の別の実装が真値になるか。
どちらも満たさないものは書かない。使った式で答え合わせをするのは禁止。

このシリーズは、その基準で作った描画を 1 本ずつ並べます。毎回、絵と一緒に成績表が出ます。コードは fullseye にあり、図はすべてスクリプト自身の出力です。

シリーズの現在地

第 1 回: 写真を一本の線にして、回る円に描かせる

手順

写真 → 点描 → 一本の閉じた線 → 等弧長に打ち直し → フーリエ係数 → 回る円で描き直し

写真の濃淡を点の密度に写し(重みつき Lloyd)、全部の点をちょうど 1 回ずつ通って戻る経路を作り(TSP)、等弧長に打ち直して複素フーリエ級数に載せます。閉じた曲線のフーリエ級数は回る腕の足し算なので、そのまま振り子の連鎖になります。

筋は TSP art(Kaplan & Bosch 2005)です。巡回路が閉じていることがフーリエ級数に載る条件なので、この 2 つはもともと噛み合っています。

門と成績

北斎「神奈川沖浪裏」(202×300 画素)を 9,000 点で描いた実測です。真値に使ったのは下界・対照群・閉形式の 3 種類だけで、特別な道具は要りません。

↑ 左が一筆書き(MST の 1.146 倍)、右が同じ点を座標順に繋いだ対照群(41.4 倍)。同じ点でも繋ぎ方でこれだけ違います。

★ 実装すると、ここで間違える

この筋道を自分で書く人が確実に踏む穴が 4 つあります。どれも例外は出ず、絵はそれらしく出るので、テストを書かないと気づけません。

1. 相関が高くても、密度は暗さに比例していない

「点の密度が濃淡を追う」を相関で測ると 0.99 を超えます。しかし相関は単調性しか見ていません。

重みつき Lloyd が作るのは重心ボロノイ分割です。Gersho の予想により、最適な点密度は重み ρ に対し √ρ に比例します(2 次元で ρ^(d/(d+2)))。つまり重みをそのまま渡すと、濃淡の半分しか濃さに出ません。

実際に指数(log 密度 を log 暗さ に回帰した傾き)を測ると 0.61 です。相関 0.9946 と指数 0.61 は矛盾しません——単調ではあるので相関は高く出ます。

濃淡に比例させたいなら、重みを二乗して渡してください。 実測で指数は 0.77 まで上がります(床と有限の点数のぶん、1 には届きません)。

2. フーリエに載せる前に、線が短くなっている

巡回路を等弧長に打ち直して標本の列にするとき、標本間隔が線分より粗いと角が切られ、線そのものが短くなります。載せる前の段階で情報が落ちるので、そのあとの「予言が当たった」は土台から崩れます。

線分が平均 2.08・中央 1.57 px の巡回路で測ると:

間隔 1.14 px では密な領域を解像できず、線が 6.5 % 縮みます。誤差はきれいに 1/N で落ちるので、標本間隔が線分の中央値を下回るまで増やすのが目安です。

3. 点数が偶数だと、フーリエ係数を二重に数える

複素フーリエ係数を k = −N/2 … +N/2 で取り出すとき、N が偶数だと両端が同じ係数を指します。両方足すとパーセバルの等式が破れます(N=8 で 1.17e-01、N=512 で 7.5e-03)。

厄介なのは、対称性の高い形では気づけないことです。正方形の輪郭で検算すると、最高周波数の係数がたまたま 0 に近いので通ってしまいます。N を 8, 9, 16, 17, 64, 255, 256 と振って、係数の本数とパーセバルの両方を固定してください。

4. ペン幅は調整するものではなく、解くもの

線が重ならない範囲では インク率 ≈ 線長 × ペン幅 / 面積 が成り立つので、目標の濃さからペン幅を逆に解けます。

ただし 2 つ注意があります。ひとつは、線を整数画素の円板で塗るとペン幅が階段になること(0.5 / 1.0 / 1.5 px が同じ絵になります)。被覆率で塗れば連続なノブになります。

もうひとつは、解いた幅でも目標には届かないこと。実測でインク率は目標の 0.890 倍でした。合併は和より小さいので実測は必ず予言以下になり、その差 11 % が線の重なりの量そのものです。ここを誤差と見ずに「重なりの測定値」と読むのが正解です。

向かないこと

墨の線画には向きません。 重心ボロノイは必ず画面全体を埋めるので、暗い画素が数 % しかない絵では点の大半が白地に配られます。鳥獣戯画(暗い画素 5.5 %)で試すと元の絵が読めません。線画は細線化して線そのものを辿るべきで、それは別の道具です。

エッジに沿わせる項は入れませんでした。 試すとエッジ近傍を通る線の割合は 0.421 → 0.530 に上がりますが、濃淡の相関が 0.911 → 0.772、疎密比が 6.5 → 2.4 と両方落ちます。この絵ではエッジが既に暗いので、濃度の対比を上げるほうが両方良くなります(エッジの面積比は 0.235 なので、偶然に線が乗る割合が 0.235。濃淡だけで 0.475 なら既に 2 倍です)。

動かす

pip install -U fullseye
py -3.11 examples/poc_one_stroke_epicycles.py \
    --image my_photo.jpg --points 12000 --out out/mine

点の数は「画素数 × 点数」で効きます。まず --points 3000 で形を見てから増やすのが早いです。

import fullseye as fs

img  = fs.read_image("photo.png")
pts  = fs.ledger.stipple_points_from_image(img, 9000, gamma=1.6)  # 濃いところに点
tour = fs.ledger.stroke_tour_closed(pts)                           # 一本の閉じた線
err  = fs.ledger.stroke_tone_error(img, tour)                      # 何を保ったかを数で
print(err["corr"], err["length_px"])

この回が作った図は全部で 9 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 6 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ Lloyd's energy sum w|x-c(x)|^2 never goes up — that is a property of the algorithm that can be checked exactly, which is why it is the gate rather than the look…

↑ the tone of a line drawing exists at a scale coarser than the line, so both pictures are compared after a Gaussian of sigma 6 px: correlation 0.984 against +0.0…

↑ while the strokes do not overlap the ink fraction is length x width / area, so the width that reproduces the mean tone is solved in closed form (0.90 px here)

↑ Parseval makes this a prediction, not a measurement: the squared error of the order-K reconstruction equals the energy of the coefficients above K, exactly

↑ every row is measured by an operator in this run; the harmonic rows are predicted in closed form before the drawing is made

↑ Parseval: the error above is the energy that was left out, exactly

第 2 回: 複素平面を「面」で見る

手順

有理関数 f(z) → 面の上で値を出す → 位相彩色 / 吸引域 / 脱出時間 → 絵とは別の真値で採点

複素解析の可視化は、曲線の上(積分・巻き数・ローラン係数・等角写像)と領域の上(値の場・位相彩色・吸引域)で別の道具が要ります。この回は後者です。

↑ どれも「きれいだから正しい」が効かない絵。採点は全部この絵の外から来る。

門と成績

おまけ: 揚力係数と薄翼理論の比は、迎角にも速さにもよらず厳密に a/b(円の半径と写像の特異点の距離の比)になります。実測のずれ 2.2e-16、厚みを 0 に近づけると薄翼理論に戻ります。厚みの効果が 1 つの数に落ちる、気持ちのよい結果です。

実装で踏む穴

既存 op が壊れて見えたが、壊れていたのは呼び方でした。 cplx_cr_residual が正則な場で 0.95、乱数で 1.04 を返し「判別していない」と結論しかけたところ、docstring に「行は虚部の増える向き」と書いてありました。こちらの格子は画像の並び(行が下向き)なので、op は共役の場を正しく測っていたのです。

規約違反は例外を出しません。もっともらしく別の量を返します。 なので食い違いは罠にせず門にしました —— 「そのまま渡すと 2」と「反転すると ≈ 0」を同じ検査に両方書いてあります。

なお「ちょうど 2」が出るのは中心差分が厳密になる 2 次以下の多項式のときだけで、一般の正則場では 2 から丸めの床の分だけ引かれます(実測 1.999951)。

向かないこと

脱出時間は反復回数の上限で切るので、境界の近くは常に「まだ分からない」が残ります(上限を上げると境界画素が増え続ける)。閉形式の内部判定は下界で、マンデルブロ集合の内部を全部は覆いません(実測 90.1 %)。位相彩色は位数を読めますが、零点の位置は読めません —— そちらは巻き数を数える op の仕事です。

最初の 1 本

import fullseye as fs

z = fs.cplx_plane_grid((-2.0, 1.0), (-1.5, 1.5), (512, 512))
basins = fs.cplx_newton_basins((1.0, 0.0, -1.0), shape=(512, 512))  # z^2 - 1
print("未収束:", int((basins < 0).sum()))          # Cayley の厳密解なら 0

この回が作った図は全部で 10 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 9 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ 有理関数 (z−z₁)(z−z₂)/(z−p) の位相彩色。

↑ ★絵が定理を証明する側に回る。

↑ 2 次の境界は虚軸 1 本(513 画素の格子で 1026)。

↑ ★★右の明るい部分は反復せずに内側と分かる領域 —— 固定点 z = (1−√(1−4c))/2 が吸引的(|2z*| < 1)なら軌道は決して 脱出しない、という閉形式。*

↑ z → z² + c は |z| を 2 乗するだけなので、c = 0 では「|z| < 1 は 0 へ、|z| > 1 は無限へ」が厳密に言える。

↑ 迎角 8 度のジューコフスキー翼まわりの非粘性流。

↑ 後縁は写像の特異点(dz/dζ = 0)なので、循環が正しくないと1/√距離 で発散する。

↑ ★CL = 2π(a/b)sin(α+β)。

↑ どの行も「使った式」ではない真値で採点している: 既存 op (cplx_winding_number / cplx_cr_residual)、Cayley の定理、閉形式の内部判定、c=0 のジュリア集合、クッタ条件の対照群、薄翼理論との比。

第 3 回: 定理が絵になる

手順

曲率の四つ組 / ファレイの分母 / 黄金角 / IFS の写像 / 曲線の位数 → 絵
                                ↓
             絵を見ずに、定理の側から答え合わせをする

アポロニウスの円、フォードの円、測地ドーム、葉序の螺旋、IFS のアトラクタ、空間充填曲線。実装が少し間違っていても、円は詰まるし螺旋は回るしフラクタルはフラクタルに見えます。

門と成績

★ここがこの回の要点です。 測地ドームの次数は別実装の op(conngraph.graph_degree_table)に数えさせています。自分で数え直して自分と一致しても、何も確かめたことになりません。

デカルトの円定理も同じ構えです。生成は反射 k' = 2(k₁+k₂+k₃) − k₄ で行うので、接触を距離から探し直して定理に入れます —— そうすると中心の計算が正しいかまで一度に効きます。

そして整数充填。種 (−1, 2, 2, 3) から始めた充填はどこまで行っても曲率が整数です。実装が少しずれれば整数から外れます —— 絵では絶対に見えない種類の誤りがここで捕まります。

実装で踏む穴

検査そのものの欠陥を 2 件見つけました。

1. pytest.raises(match="similarit") は op 名 ifs_similarity_dimension に当たって必ず通ります。 しかもプリセット名を間違えていて("fern"、正しくは "barnsley_fern")、「未知プリセットの拒否」を見て合格していました。「拒否した」だけを見て「なぜ拒否したか」を見ないとこうなります。
2. 型を points で登録して型契約の門に捕まりました。points は (N,3) の 3 次元点で、2 次元の点列は pairs です。名前が似ているので取り違えやすい。

そして図が数値検査の盲点を暴きました。 測地ドームの図を三角形 1 枚ずつ焼いていたため、呼び出しごとに自分の外接矩形へ正規化され「各三角形が画面いっぱいに伸びた線の山」になっていました。検査は 20 件とも通っていたので、図を開くまで気づけません。

向かないこと

アポロニウスの充填は depth に対し円が 3 のべきで増えます(depth 10 で 118,100 個)。フォード円は分母の二乗で半径が縮むので、分母 30 を超えると多くが 1 画素を割ります。IFS のカオスゲームは点を打つので、細かい隙間は点数を増やさないと埋まりません(箱数えの次元が閉形式より上に出るのはそのため)。モランの式は相似写像にしか使えません —— アフィンだが相似でないバーンズリーのシダは拒否します。

最初の 1 本

import fullseye as fs

t = fs.ledger.circle_packing_apollonian(curvatures=(-1, 2, 2, 3), depth=4)
print(len(t["x"]), "円 / 曲率は整数のまま:",
      abs(t["curvature"] - t["curvature"].round()).max())

pts  = fs.ledger.phyllotaxis_pattern(n_points=800)
gaps = fs.ledger.neighbour_index_gaps(pts, k=6)
print("番号差の山:", gaps.argsort()[::-1][:5])      # 黄金角なら 8, 13, 21, ...

この回が作った図は全部で 11 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 10 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ 合計は 2·3^depth + 2。

↑ フォードの円(分母 12 まで、47 円)。

↑ 分割 f を 1, 2, 3, 4, 6 と上げても、次数 5 の頂点はちょうど 12 個から動かない(次数 6 だけが増える)。

↑ オイラーの公式 V − E + F = 2 の帰結。

↑ 3 枚とも螺旋に見える。

↑ 黄金角の山は [34, 55, 89, 21, 13, 8, 144] —— すべてフィボナッチ数。

↑ モランの式 Σrᵢᵈ = 1 は写像の縮小率だけから d を出す。

↑ order 5 で 4^5 = 1,024 点をちょうど 1 回ずつ通り、隣り合う点は必ず距離 1。

↑ ヒルベルトは √k の近くを通り(k=32 で 6.38)、走査線はほぼ k に比例する(16.07)。

↑ どの行も「使った式」ではない真値で採点している: デカルトの円定理、整数充填、|ps − qr| = 1、オイラーの関数の和、オイラーの公式(既存 graph_degree_table が数えた)、フィボナッチ、モランの式と既存 fractal_dimension、置換であること。

第 4 回: うなりは一つ ―― 干渉縞と印刷のモアレは同じ数学

手順

膜の固有モード / 二重スリット / 回折格子      ┐
                                          ├→ どちらも「2 つの周波数ベクトルの差」
網点 2 版の重ね / 彫版線 / ハッチ / モザイク  ┘

二重スリットの縞と、2 版を重ねた網点のモアレは教科書では別の章に載っています。でも式は 1 本しかありません。 この回はその 1 本を両側から確かめます。

↑ モードの形は符号つき([-1, 1])。節線の両側で位相が反転するので、image2d を名乗ると腹の位相が消える。

門と成績

気をつけること(2 つとも一度読み違えました)

うなりは低周波側にあります。 FFT の探す範囲を切らずに argmax を取ると、どの角度でも 16.7 px = スクリーンの周期 1/f そのものが返り、「op の予言と全然合わない」と読めます。

ハッチの線は高周波側にあります。 切らずに取ると、線ではなく「どこに墨を置いたか」の濃淡の包絡(元の絵と同じ 24 px)を拾い、向きが 90° ずれて見えます。

どちらも op は正しく、欠陥は測り方にありました。スペクトルには必ず複数の山があり、argmax は「一番強い山」であって「探している山」ではありません。しかも返ってくるのはもっともらしい数なので、例外も NaN も出ずに静かに別の量を答えます。

クラドニ図形についての訂正

よく見る cos·cos − cos·cos の式は、板ではなく膜の解です。実際のクラドニ板は重調和方程式に従う別物で、砂が節線に集まる絵は同じでも周波数比は合いません。ここでは膜と明記し、膜の真値(閉形式の固有値・ベッセル零点)だけで採点しています。

もう 1 つ。円膜の節円は J₀ の零点にあります。J′₀ の零点は腹の位置です —— 取り違えると「op が間違っている」と読み違えます(一度やりました)。

向かないこと

wave_two_slit はフラウンホーファー領域の遠視野で、近接場やフレネル領域は別の op の仕事です。膜のモードは矩形と円の解析解だけで、任意形状は解きません。wave_fringe_period は1 本の支配的な周期を返すので、複数の周期が重なった像では最大の山しか見ません。モアレの予言は同じ形のスクリーン 2 枚の場合で、3 版以上や点形状の違いによる 2 次のうなりは入りません。hatch_field は画素ごとに 1 つの向きしか持てないので、交差や分岐では向きが任意になります —— そこを信じてよいかを言うのは向きではなくコヒーレンスです。

最初の 1 本

import fullseye as fs

# 刷る前にモアレの周期が分かる
p = fs.ledger.halftone_moire_period(lpi_a=60.0, angle_a_deg=45.0,
                                    lpi_b=60.0, angle_b_deg=75.0, pixel_um=25.4)
print("モアレの周期 %.2f px / 向き %.1f 度" % (p["period_px"][0], p["angle_deg"][0]))

# 縞間隔は、作った op とは別の op が測り返す
img = fs.wave_two_slit(wavelength_nm=550.0, slit_sep_um=200.0,
                       distance_mm=200.0, shape=(64, 1024), pixel_um=5.0)
print("実測 %.2f px / lambda*D/d = %.2f px"
      % (fs.wave_fringe_period(img), 0.55 * 200e3 / (200.0 * 5.0)))

この回が作った図は全部で 13 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 12 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ (m, n) モードの節線は縦 m−1 本・横 n−1 本。

↑ 縞間隔は λD/d。

↑ 重なっていることが主張。

↑ ★デューティ 50 % の矩形格子は偶数次が消える(矩形波のフーリエ係数が偶数調波で 0)。

↑ 新聞の 15° / 45° / 75° は差が 30° で、うなりが最も細かくなる組み合わせ。

↑ ★FFT の探す範囲を切らないとスクリーン自身の山(どの角度でも 16.7 px = 1/f)を拾って「予言と全然合わない」と読める。

↑ 平均絶対差 0.0517。

↑ 被覆率は線幅 / 間隔 = w/d。

↑ 向きは既存 structure_tensor_orientation が決める(この族は構造テンソルを再実装しない)。

↑ Lloyd 反復でセルが等エネルギーに近づく。

↑ 下がり続けることが主張。

↑ どの行も「使った式」ではない真値で採点している: 膜の閉形式とベッセルの零点、λD/d を測り返す別の op、既存 grating_wavelengths との往復、描く前のモアレ予言、w/d の閉形式、既存 structure_tensor_orientation、既存 stipple_energy。

第 5 回: 絵では確かめられないもの ―― 力学系と極小曲面

手順

名前で呼んだベクトル場 → 積む → 断面・指数・分岐・次元
極小曲面の厳密なパラメータ表示 → メッシュ → 平均曲率を別の op に測らせる

ローレンツ・アトラクタの図は、積分器が 1 次でも 4 次でも蝶に見えます。 極小曲面の図は、平均曲率が 0 でなくてもきれいな曲面に見えます。この回は「見て分かる」が一切効きません。

門と成績

★採点者が別の族から来ます。 ode_vector_field_grid は flow2d を返すので、PIV(粒子画像流速計)族の既存 op —— piv_divergence / piv_vorticity —— がそのまま食います。あちらはこの族の実装を何も知らないので、「場が正しい」という主張の独立した採点者になります。線形系では発散が厳密に tr(A)、渦度が A₁₀ − A₀₁ です。

極小曲面は定義そのものが門になります。 「平均曲率 H が至るところ 0」という定義を、既存 vertex_curvature がそのまま測ります —— つまり「これは極小曲面だ」という主張を、作り方を知らない op が採点します。 対照群の単位球が 1.00004、半径 1 の円柱が 0.50000 に出るので、門が素通しでないことも同時に分かります。

残差を隠さない

2 か所とも、理由を書く前に正体を実測で特定しました。

(1) 分岐点は必ず少し手前に見えます。 分岐点での収束は幾何的でなく代数的なので、有限の burn-in では届きません(δ も 4.7485 対 4.6692 と 1.7 % 偏る)。だから門は「誤差が小さい」ではなく 「burn-in を伸ばすと置いていった分だけ縮む」 で置きました —— 2,000 → 20,000 で誤差が 10.9 倍縮みます。残差の正体が写像でなく測り方であることの証拠は、そちらです。

(2) 相関次元の偏りは「点数不足」ではありませんでした。 円 1.0061・カントール 0.6408(真値 log2/log3 = 0.6309)は合うのに、平面は 1.8789(真値 2)。最初に疑った点数は外れでした —— 400 → 1,500 → 3,000 点で 1.8825 / 1.8789 / 1.8709 と、むしろ下がります。

正体はべき乗則を見る半径の窓です。既定は対距離の 1〜25 パーセンタイルで、有界な集合ではその上端が境界に当たり、相関和が飽和して傾きが寝ます。窓を小さい r 側に寄せると 1.887 → 1.947 → 2.050 と真値に寄ります。

「有限標本だから」「数値誤差だから」は何にでも当てはまるので、書いた瞬間に調べるのをやめてしまいます。 当てはまる説明と正しい説明は別物で、正しい説明だけが「どのノブを回せば直るか」を教えます。この件では r_lo / r_hi であって max_points ではありませんでした。

ジャイロイドも同じ構えです。level = 0 で場は体心反転に対し奇なので体積比は対称性からちょうど 1/2(格子に依らない)。ですが Schoen のジャイロイドは H = 0 で、節面近似には残差が残ります(|H| / 主曲率スケール = 0.1250)。足場と絵として使うもので、極小性の主張ではありません。

実装で踏む穴

引数で戻り値の型が変わる関数は、型付きの台帳に載せられません。 gyroid_isosurface(..., thickness=...) は thickness を渡すと voxel、渡さないと mesh を返す設計でしたが、宣言できる出力型は 1 つなのでどちらを書いても嘘になります。2 つの op に分けました。

分けた副作用が利得になりました —— 固体形が scikit-image の依存から外れ、numpy だけの環境でも使えるようになったのです。

もう 1 つ。ode_vector_field_grid を最初 (H, W, 2) で返していて型検査に捕まりました。そちらのほうが自然に読めますが、flow2d は既存の型で既存の述語と既存の消費側があるので、転置を返すのは型の嘘です。読みやすさより既存の語彙に合わせた結果、下流の既存 op がそのまま繋がりました。

向かないこと

系は族名か係数配列で受け、callable は取りません —— 型付き台帳は入力を型で登録するので、関数は載らないのです。任意の常微分方程式を渡す口ではありません。リアプノフ指数は接流の解ける系だけで、遅延埋め込みからの推定は入りません。分岐図は 1 次元写像だけです。相関次元は答えが要るときは半径の窓を明示してください(上のとおり既定の窓は有界集合で低く出ます)。極小曲面は閉形式の 4 族 + ジャイロイドで、任意の境界に張るプラトー問題は解きません。管メッシュは自己交差を検出しません(曲率半径より太い管は自分を貫きます)。

最初の 1 本

import numpy as np
import fullseye as fs

# 積分器の次数を実証する(厳密解は expm(At)x0)
A, x0, T = np.array([0.0, 1.0, -1.0, 0.0]), np.array([1.0, 0.0]), 4.0
exact = np.array([np.cos(T), -np.sin(T)])
for dt in (0.04, 0.02, 0.01):
    st = fs.ode_flow_states("linear", A, x0, T, dt, "rk4")
    end = np.array([st["x0"][-1], st["x1"][-1]])
    print("dt = %.3f  誤差 %.3e" % (dt, np.linalg.norm(end - exact)))   # 半分で 1/16

# 「これは極小曲面だ」を、作り方を知らない op が採点する
V, F = fs.ledger.minimal_surface("catenoid", nu=90, nv=140, extent=1.2)
H = np.abs(np.asarray(fs.ledger.vertex_curvature((V, F))))
print("|H| 中央値 %.5f(対照群の単位球は 1.0)" % np.median(H[np.isfinite(H)]))

この回が作った図は全部で 11 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 10 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ 厳密解は expm(At)x₀ なので、誤差は積分器の次数をそのまま出す。

↑ ローレンツ・アトラクタの xz 断面。

↑ Σλ = -13.66666、閉形式 −(σ+1+β) = -13.66667(差 2.6e-06)。

↑ 円 1・平面 2・カントール log2/log3 = 0.6309。

↑ ode_vector_field_grid は flow2d((2, H, W) の (dy, dx))を返すので、piv_divergence / piv_vorticity / piv_flow_magnitude がそのまま使える。

↑ 同じ op で両方を測っている。

↑ ローレンツの軌道を管メッシュに(MATLAB の tubeplot 相当、12510 頂点 / 25000 面)。

↑ 「これは極小曲面だ」という主張を、作り方を知らない op が採点する。

↑ level = 0 で場は体心反転に対し奇なので、両側の体積はちょうど半分ずつ(実測 0.499928 / 0.499981 / 0.499991、格子に依らない)。

↑ どの行も絵からは読めない: expm との誤差の比、トレース恒等式、厳密な分岐点、既知の 3 つの次元、PIV 族の既存 op が測った tr(A)、トーラスの解析解、既存 vertex_curvature が測った H —— そして節面近似の残差は隠さずここに出してある。

第 6 回: 目が嘘をつく絵を作って、測る側を採点する

手順

錯視図を作る → その図が「否定している不変量」を数で返す → 既存の op に測らせる
無限に描き続ける系を回す → 絵ではなく恒等式で採点する
時間依存の量をすべて θ の関数にする → 継ぎ目が最初から存在しない動画になる

ここまでの 5 回は「描いた絵を、絵の外から採点する」話でした。第 6 回は向きが変わります —— 採点される側の絵を、こちらから作ります。

錯視を置くつもりでした。カフェウォールの目地は厳密な水平線で、ミュラー・リヤーの 2 本は厳密に等長で、チェッカーシャドウの 2 マスは浮動小数として同じ値 —— それを op に測らせれば「測る側は騙されない」と言える、と考えたからです。

錯視の章は落としました —— 測っても何も言えなかったから

言えませんでした。 検査が 7 件とも自分で置いた定数を読み返すだけだったからです。厳密に水平に引いた目地が水平だ(傾き 0.0e+00)、平坦に塗ったマスが平坦だ(ptp 0.0e+00)、同値に置いた 2 パッチが同値だ。間に変換が 1 つもありません。 錯視は人間の知覚の性質なので、測る側が騙される理由がそもそも無く、「測る側は騙されない」は発見ではなく同義反復です。絵の側も直線と矩形を引いているだけで、この箱の力を何も示していません。

残したのは 1 件だけ、測る側の偏りです。

↑ 2 本の軸は矢羽根を外すと画素単位で同一(差 0)。ところが「行の黒い区間」を素朴に測ると 223 画素と 225 画素になります。しかも矢羽根が無くても 223 画素(真値 220)—— 反エイリアスの裾を「黒い画素」として数えているからで、これは錯視より前の素の欠陥です。

↑ ★ここでも外しました。「矢羽根が長いほど偏る」と読みましたが、20 px で 2 画素に達したあと 80 px まで動きません。偏りを作るのは長さではなく、端に裾が乗るかどうかだけでした。

止めるまで描き続けるものを、どう採点するか

8 ビットの一行プログラムから続く系譜 —— 10 PRINT のランダム迷路、基本セルオートマトン、ラングトンの蟻、カオスゲーム。これらは「いつ止めても途中」なので、絵では採点できません。

そこで、絵とは独立に成り立つ主張を持つものだけを入れました。

↑ 規則 90 を中央 1 点から 200 行。第 n 行の第 k セルは二項係数 C(n,k) の偶奇に厳密に一致します —— パスカルの三角形を 2 で割った余り。絵を一切見ずに採点できる、という意味でこの族の代表例です。

↑ 始まりはどれも中央の 1 セルだけ。変えたのは 8 ビットの規則番号 1 つです。規則 30 は乱数生成に使われた歴史があり、規則 110 は万能計算ができることが証明されています。

↑ 規則は 2 つだけ(白なら右折して黒く塗る、黒なら左折して白く塗る)。最初の 1 万歩ほどは対称な模様、次に無秩序、そして突然「高速道路」に入り、以後は周期 104 で斜めに進み続けます。

↑ 黒マスの数を歩数に対して。前半は揺れますが、高速道路に入ると 12/104 = 0.1154 マス/歩の直線になります。絵を見なくても、傾きで「入った」が分かります。

↑ 隙間に接する円を入れ続けるので終わりがありません。互いに接する 4 円の曲率はデカルトの円定理 (Σk)² = 2Σk² を厳密に満たします(残差 1.2e-16)—— 円の位置まで複素数の曲率中心で解けるので、当てはめではなく代数で描いています。 ★塗りの明るさは曲率(= 1/半径)—— 定理が試している量をそのまま色にしてあるので、絵の濃淡が主張と同じものを指しています。

↑ 賽を振って頂点へ半分ずつ寄るだけで、シェルピンスキーの三角形が出ます。明るさは点の密度(双一次で積み、asinh で暗部を伸ばしてあります)。描き方は飾りではありません —— 同じ点列でも、切り捨てて積んだ絵から測った箱数え次元は真値 log3/log2 から −0.032 ずれるのに、こう描くと −0.0030 で 10 倍正確になります。

↑ バーンズリーのシダ —— 4 本のアフィン写像を確率で選んで回すだけ。明るさは密度で、正規化は最大値でなく分位点で行っています(根元の 1 点が飛び抜けるので、最大値で割ると他が全部潰れます)。

↑ 場はポテンシャル ψ の回転として作ってあるので、発散が恒等的に 0(非圧縮、残差 7.0e-17)。流線が渦のまわりで閉じているのがその絵です —— 「それらしい雑音」で作った流れ場では、ここが湧き出しと吸い込みになります。 ★線の色と太さは流速 |∇ψ| —— 太さも濃さも一定だと、同じ場から描いても落書きにしかなりません。

↑ グレイ–スコット反応拡散。種は散らして 12 個置いてあります ——中央に 1 個だけだと、同じ歩数でも模様の被覆が 0.086 にしかならず、画面の大半が余白になります(散らすと 0.379)。★格子状に並べると格子そのものが絵に出るので、ずらして置きます。

↑ 左は 10 PRINT CHR$(205.5+RND(1)); : GOTO 10 —— Commodore 64 の一行プログラム。右はトルシェ・タイルで、弧の端点が必ず辺の中点に来るので、どう組んでも曲線が途切れません。

時間軸で循環する画像

↑ 最後のコマから最初のコマへ戻るところに、切れ目がありません。

作り方が要点です。「最後に頭へ戻す」のではなく、時間依存の量をすべて θ の関数にして θ を 0→2π 回すと、t = T は t = 0 と同じ式になります —— 継ぎ目は編集で消すものではなく、最初から存在しません。

↑ 継ぎ目の正解は 0 ではなく 1 です。 比 =(最後のまたぎの差)÷(ふつうのコマ間の差)。1 なら「最後のまたぎが、ほかのまたぎと見分けが付かない」。0 だと継ぎ目が無いのではなく、動きが止まっているという意味になります —— ここを取り違えると「完璧な結果」に見えます。

門と成績

検査が見つけたもの

4 件出ました。どれも絵を見ても気づけない型です。

1. 二項係数を int64 で積んでいて、C(62,31)×31 で黙って溢れた。 リュカの定理の帰結「C(n,k) が奇 ⟺ (n & k) == k」に置き換えて、桁あふれの余地を消しました。
2. チェッカーシャドウで対にするマスを間違えていた。 同じ明暗のマスを 2 つ選んでいて、差が 0.31 出ました。対にすべきは「影の外の暗マス」と「影の中の明マス」。絵は"それらしく"見えるので、数で測るまで気づけません。
3. 「104 歩で 52 マス」と見当で書いた成長率が、実測 0.114/歩 と合わなかった。 測り直すと正味ちょうど 12 マス(整数)。周期の中で塗っては消すので、正味はずっと少ない。推測でなく実測が正しかったという、それだけの話です。
4. ★アポロニウスで平方根の枝を選び損ねていた。 複素曲率中心には枝が 2 つあり、片方は「親 3 円に接しない円」を返します。主値だけを使っていたため、外円の外へ逃げる円の鎖が出ていました —— しかも絵としては「フラクタルっぽい」ので、画像を開いて見るまで分かりません。接することを数で確かめて枝を選ぶようにし、「外へ逃げた円が無い」を門にしました。

そのほかに、定数の配列でも std() は厳密に 0 にならない(平均を引く途中で 1 ulp 残る、実測 5.6e-17)という小さな発見がありました。「全部同じ値」を厳密に言いたいときは ptp(最大 − 最小)を使います。

速度も 1 件。素朴に全画面の距離場を primitive ごとに作っていたので、錯視 12 枚で 17.7 秒、アポロニウス 1 枚で 4.96 秒かかっていました。外接箱だけを触る方式にして 0.07 秒 / 0.03 秒(250 倍・165 倍)。PoC の所要時間は CI の合否そのものなので、これは配管ではなく門の問題です。

向かないこと

作品としての完成度には向きません。この族には主張が検算できる図しか入れていないので、美しいが検算できない図は、描けても入れていません。逆に言えば、ここにある絵は全部「なぜそれで正しいと言えるのか」に答えられます。

また、錯視の強さは測っていません(「どれくらい傾いて見えるか」は心理物理の話で、この族の守備範囲の外です)。測っているのは「傾いていないこと」だけです。

この回が作った図は全部で 16 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 1 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ この PoC の主張は全部この表に入っています。

第 7 回: 測定器のほうを健診する ―― キャリパーが外すのはどこで、なぜか

手順

答えの分かっている絵(錯視)を作る → 産業用のサブピクセル測定器に測らせる
→ 真値と突き合わせる → 外れたら「なぜ外れたか」を突き止め、直せるか試す
→ 直らなかったものは、直らなかったと書く

第 6 回で「真値の付いた絵」が作れるようになりました。第 7 回はその絵の使い道です ―― 測られる側ではなく、測る側を採点します。

Fullseye には HALCON 流のサブピクセル測定器が入っています(measure_pos / measure_pairs / fuzzy_measure_pairing / apply_metrology_model)。工場で寸法を測るための道具です。この道具に、目が確実に外すと分かっている図を測らせるとどうなるか。

★この回で新しい op は 1 つも作っていません。作る側と測る側が同じ箱に入っているので、答え合わせが 1 本のスクリプトで閉じます。

まず、測定器は目より正しい

↑ 青が与えた参照線、橙が測定器が法線方向に探して見つけたエッジ点 41 個。タイルは傾いて見えたままなのに、点は一直線に乗り(rms = 0.0)、当て直した直線の角度は厳密に 0.0 度。

↑ 左の円は小さく、右の円は大きく見えます。円の計測オブジェクトが返した半径は 29.996525 px と 29.996525 px ―― 差は 2.5e-12 px。目が主張する差は、画像の側には存在しません。

↑ 縦軸を 0.2 px 幅まで拡大しても 2 点は重なったまま。どちらも指定の 30 px より 0.0035 px 小さいのですが、そのずれは両方に同じだけ乗るので、差を取ると消えます。

ところが、探針 1 枚では足りなかった

ここが、この回でいちばん学びになったところです。

最初、カフェウォールを 1 通りの設定(ずらし量 0.25)だけで測って、「測定器は錯視では決して動かない」と書きました。きれいな結論です。

そのあと、ずらし量 5 通り × 探索半幅 4 通りの 20 通り全部を測りました。

↑ 探索半幅が目地の帯幅 8 px に届く 10 px と 16 px のとき、ずらし量 0.125 と 0.375 でだけ 0.1429 度傾きます。0.25 は、たまたま外れない側でした。

↑ 何が起きているか。全体図では 8 px の飛びが潰れて見えないので、目地の周り 34 行だけを等倍で切り出しました。下の段では、点が帯の反対側にも乗っています ―― 探索範囲が向こう側の境界に届いているのです。

既定に見える 1 点だけで試すと、この欠陥は構造的に見つかりません。

外したことは、答えではなく残差に出る

そして、ここがいちばん実務に効くところです。

↑ 同じ 20 通りについて、答え(角度)と残差(rms)を重ねたもの。

外れたときの角度は 0.1429 度です。数字だけ見れば「ほぼ 0」で、そのまま通ってしまいます。ところが残差は 0.00 から 1.70 へ桁で動きます。しかも「角度 0 ⇔ rms 0」が 20 通り全部で一致しました。

門は答えではなく残差に置く。「当てた形が本当にそこにあったか」は、答えの側からは見えません。

測定器も外す。ただし目とは別の理由で

↑ 上(矢羽根なし)は測定線上のエッジが 2 本、下(矢羽根あり)は 4 本。

矢羽根は「長く見せる」だけではありません。測定線の上に余計なエッジを置きます。測定器は隣り合う逆極性のエッジを対にするので、軸ではなく矢羽根のストロークを測ってしまいます(幅 5.701 px)。

目は「長く見える」と外し、測定器は「別の物を測る」と外す ―― 外れ方が違います。

↑ 同じ 1 本の軸が、測り方で 4 通りの答えになります。どれもバグではありません。

測り方を宣言すると、静かな誤りが見える拒否に変わる

↑ fuzzy_measure_pairing に「幅は 220 px のはず」と宣言して測らせたもの。矢羽根なしでは適合度 0.9988、矢羽根ありでは 0.0225。

これは「間違った値を返した」のではありません。「想定した幅の構造は見つからない」と数で申告している状態です。隣り合う逆極性の対しか候補にしないので、軸をまたぐ対はそもそも提案されません。

測り方を宣言した分だけ、失敗が静かな誤りから見える拒否に変わります。

3.777 px は錯視ではなく、端の定義だった

矢羽根が無い図でも、答えは 223.777 px で真値 220 と合いません。錯視のせいでしょうか。

↑ 軸の長さを 140〜300 px に振ると、ずれは +3.777 px 一定(振れ幅 0.0000 px)。

比例していれば倍率の誤り、一定なら端の定義のずれです。軸は太さを持って描かれているので、「外側エッジ間の距離」は「線分の長さ」より線幅ぶん長い。引けば真値に戻ります。

動くのは答えではなく、証拠の数のほうだった

↑ 収束して見える 7 本すべてが、測ると水平(最大 |傾き| 4.3e-04 度 = 幅 520 px を渡って 0.0039 px)。

↑ 傾きは 7 本とも 0 に重なって見分けが付きません。動いたのは採用されたエッジ点の数で、ハッチの向きが 1 本おきに反転するため 61 点中 13 と 19 を交互に取ります。

「測れた」と「よく測れた」は別の量です。前者だけ見ていると、この差は見えません。

直らなかったものは、直らなかったと書く

↑ 帯の左右の 2 本は厳密に 1 本の直線です(帯で隠れているだけ)。橙が素朴なキャリパーの点、青が両エッジを対にして中心を取り直した線。

ここでは測定器が外します。共線のはずの 2 区間で、角度が 0.450 度違いました。

原因は突き止められました。極性が塊で交互しているのです(41 点が 9 個の塊、+18 / −23)。細い線には上側と下側の 2 つのエッジがあり、「線の位置」はどちらか(あるいは中心か)を言うまで定義されません。測定器はその 2 つを行き来していました。

*↑ 両エッジを対にして中心を取ると、角度差は 0.450 → 0.115 度、残差は 2.07 → 0.10/0.39 まで下がります。**それでも 0 ではありません。*帯の縁の細い縦線と画像の端が測定線に掛かっており、そこまでは追い込んでいない ―― というのが正直な状態です。

絵を作り直したら、測れる真値が変わった

濃淡と色を付け直しただけで、絵から測った数が真値に近づきました。同じ点列・
同じ op・同じしきい値で、箱数え次元のずれが −0.032 → −0.0030(10 倍)。

原因は積み方でした。前の実装は .astype(np.int64) で切り捨てて整数画素に
積んでいた —— 四捨五入ですらないので、絵全体がちょうど半画素ずれます
(実測の重心差 −0.5009 / −0.4998 画素)。総量は 39.9 万点で相対差
0.00e+00 と合うので、0 次だけ見る検査はこの欠陥に構造的に盲目です。
双一次で積むと重心差は 0.000e+00 —— 1 次元で i(1−f) + (i+1)f = i + f に
戻るからで、2 次元は軸ごとに分かれるので両方とも厳密です。

トーンは asinh にしました。狭義単調なので画素の大小は 1 組も入れ替わりません
(4,998 組で検査)。ただし正規化を最大値でやると、バーンズリーのシダは根元の
1 点が飛び抜けるので他が全部潰れて真っ暗になりました。分位点に替えて
std 0.052 → 0.0996。反応拡散は種が中央に 1 個だけで、模様の被覆が
0.086 しかありませんでした —— 散らして 12 個置くと同じ計算量で 0.379。
★格子状に並べると格子そのものが絵に出るので、ずらして置きます。

★外したこと: アポロニウスの色が大きい円で潰れていたので「順位で
正規化すれば階調が配り直される」と読みました。逆効果でした ——
極小円が多数派なので、大きい円が全部下端に潰れて前より悪くなりました。
対数に戻してガンマ 0.45 で低域を持ち上げるのが正解でした。
「分布を均す」は、見せたい対象が少数派のときは逆に働きます。

正反対の主張の図が、同じ絵だった

超立方体を 2 枚の面で同時に回す章で、「比 1 : 2 は 60 コマで閉じる」と
「比 1 : φ は閉じない」の GIF が バイト単位で同一でした。

原因は投影です。図は落とした 3 次元の x, y の 2 成分だけを描いていました。
rot4(0, t1, 1, t2) の t2 は zw 面の回転なので、x と y に一切効きません
—— 60 コマ通して x,y の差は厳密に 0.000e+00、x,y,z なら 2.494。
章の主題が、その投影では原理的に見えなかったわけです。3 次元を一般の向きから
見る 2×3 の正射影(正規直交性 1.11e-16)に替えて差 1.322 になりました。

★この嘘が通ったのは、「2 枚の図が同じ中身か」を見る門がどこにも無かった
から。枚数も名前も説明も存在も検査していたのに、そこだけ穴でした。書いた図の
sha256 を記録して、同一バイト列が 2 枚出たら欠陥にする門を入れました。
全 151 本を走査して該当は 2 組だけ(もう 1 組は図の枚数が変わって位置番号が
ずれた置き去りファイル)でした。

錯視の章は、測っても何も言えなかった

第 6 回で落とした話です —— 詳しくはその節に書き直してあります。要点だけ: 検査 7 件が自分で置いた定数を読み返すだけで、間に変換が 1 つもありませんでした。残したのは測る側の偏り 1 件だけです。

作った図は 1,155 枚、記事に出ていたのは 317 枚だった

数え直して分かりました。展示館が参照していたのは 27 % で、848 枚が
記事から辿れない状態でした。1 展示あたり 2.1 枚(150 件のうち 139 件が
ちょうど 2 枚)で、PoC は中央値 6 枚・最大 25 枚作っているのに、
残りへの道が記事に無かった。

出す枚数を増やして 318 → 767 枚にし、それでも出ない分のために、
展示ごとに「この回が作った図は全部で N 枚 —— 全部見る」を置きました
(案内 1,155 / 1,155 枚)。追加分の説明は最初の 1 文に切ってあります ——
枚数を増やすのが目的で、詳しい説明は主図と 2 枚目が持っているからです。

健診結果

↑ この表に「絵を見て判断した」行は 1 つもありません。真値は生成器が保証し、答えは既存の測定 op が返したものです。

真値をそのまま回復したのは 3 行、外れたのは 5 行。外れた内訳は 3 通りでした ―― 決めれば戻るもの(探索半幅・線幅の下駄)、拒否に変えられるもの(想定幅の宣言)、そして 戻らなかったもの(ポッゲンドルフ)。

検査が見つけたもの

この回で自分の書いたものが 2 か所、データと食い違っていました。どちらも「生成された図説を読み返して」見つけたものです。

1. カフェウォールの掃引の図説に「半幅を変えると採用点の数が変わる」と書きましたが、4 通りとも 41 点でした。証拠の数が実際に動くのはツェルナーのほうです。
2. 健診表の図説に「外した行は 2 行・決めれば戻る」と書きましたが、実際は 5 行で、ポッゲンドルフは戻っていません。

図説は、書いた時点ではなく、生成されたものを読んで確かめる。 数字を % で埋め込んで自動生成していても、その周りの文章は手で書いているので、同じようにずれます。

向かないこと

この PoC は「測定器は錯視に強い/弱い」を言うものではありません。測り方を決めていない分だけ答えが割れることを、真値つきの図で数として見せるものです。工場で使う測定線の張り方(探索半幅・想定幅・エッジの極性)を決めるときの、机上の確認には使えます。実写のノイズ・照明むら・被写界深度は入っていません。

この回が作った図は全部で 16 枚あります —— 全部見る

この回の残りの図

以下は同じ実行が作った残り 1 枚です。説明は 1 文目だけに切ってあります —— 詳しい読み方は上の主図が持っています。

↑ measure_length は「参照線の法線方向に何 px 探すか」。

出典

• 画像: 葛飾北斎「神奈川沖浪裏」(富嶽三十六景、1830–32 年頃)。メトロポリタン美術館 Open Access が isPublicDomain: true(CC0)で公開している版(object 45434 / 画像 DP130155)を輝度化して縮小。
• 手法: Craig S. Kaplan and Robert Bosch, "TSP Art", Computational Aesthetics (2005)。
• 点密度の指数: A. Gersho, "Asymptotically optimal block quantization", IEEE Trans. Inform. Theory 25(4), 1979。
• ニュートン吸引域: A. Cayley, "The Newton–Fourier imaginary problem", Amer. J. Math. 2, 1879。
• 円の充填: R. Descartes (1643) と F. Soddy, "The kiss precise", Nature 137, 1936; J. Lagarias, C. Mallows, A. Wilks, "Beyond the Descartes circle theorem", Amer. Math. Monthly 109, 2002。
• フォード円: L. R. Ford, "Fractions", Amer. Math. Monthly 45, 1938。
• 葉序: H. Vogel, "A better way to construct the sunflower head", Math. Biosci. 44, 1979。
• 相似次元: P. A. P. Moran, "Additive functions of intervals and Hausdorff measure", Proc. Camb. Phil. Soc. 42, 1946。
• 膜の固有値: Lord Rayleigh, The Theory of Sound, 1877(クラドニの砂の図は板であって、ここで解くのは膜)。
• 二重スリット: T. Young, 1804。
• モアレのベクトル解析: B. Oztan, G. Sharma, R. P. Loce, "Misregistration sensitivity in clustered-dot color halftones", J. Electronic Imaging 17, 2008。
• 点描と量子化: A. Secord, "Weighted Voronoi stippling", NPAR 2002; S. Lloyd, "Least squares quantization in PCM", IEEE Trans. Inf. Theory 28, 1982。
• 力学系: E. N. Lorenz, "Deterministic nonperiodic flow", J. Atmos. Sci. 20, 1963; O. E. Rössler, "An equation for continuous chaos", Phys. Lett. A 57, 1976; M. J. Feigenbaum, J. Stat. Phys. 19, 1978; P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50, 1983; G. Benettin et al., Meccanica 15, 1980。
• 極小曲面: H. A. Schwarz, Gesammelte Mathematische Abhandlungen, 1890; A. H. Schoen, "Infinite periodic minimal surfaces without self-intersections", NASA TN D-5541, 1970。
• 管の枠: R. L. Bishop, "There is more than one way to frame a curve", Amer. Math. Monthly 82, 1975。