fullseye

光学(レンズ・回折・偏光) — 使い方ガイド

この族は何をする道具箱か

レンズより上、画素より下の層です。産業ビジョン(検査ライン)でも Physical AI(ロボット知覚)でも、画像処理を始める前に誰も撮っていない決定があります — どの焦点距離か、絞りはどこか、被写界深度はどれだけか、回折で潰れる最小欠陥は何 µm か、偏光板でテカりは消えるか。それらは全部閉形式の計算で、この族はそれを第一級の op にしたものです。21 op / 4 カテゴリ(numpy + scipy のみ、台帳は opsoptics.py、実体は optics.py):

データ種は既存語彙の再利用が基本です: table(dict、計測値の束/ABCD 素子リスト)、matrix(ABCD 2x2・Mueller 4x4 の実行列 — mat_svd や mat_cond にそのまま流せる)、image2d(PSF・開口・強度像)、cimage(2-D complex = 複素場と Jones 行列)、pairs((n,2) の曲線 = MTF・cos⁴)。新語は 2 つだけ、jones(長さ 2 固定の complex ベクトル)と stokes(長さ 4 固定 + 偏光度 ≤ 1 の物理制約)で、これは signal / cpoints に相乗りさせると「256 点の正弦波を Stokes 枠に渡せる」という型レベルの嘘になるためです。

既存 op との棲み分け(重複させていないもの)

やりたいこと 使う op 置き場所
面で曲がる実際の光線(反射・Snell 屈折・Fresnel 反射率・全反射) reflect / refract / snell_angle / fresnel_reflectance / normal_from_reflection match3d(この族は近軸・スカラなので、光線と面の相互作用はそちら)
円板画像から Zernike 係数をフィットする fit_zernike match3d(wavefront_stats はその返り dict をそのまま食い、match3d 自身の基底ビルダーを再利用するので規約がずれない)
PSF によるぼかし・逆畳み込み vol_gaussian_psf / vol_richardson_lucy / cx_wiener_deconvolve volrestore / complexops(psf_to_mtf は特性化するだけで復元はしない)
FFT・複素画像・位相アンラップ cx_fft 系 / phase_unwrap complexops
位相シフト干渉法・縞投影(N-step) wrapped_phase / unwrap_phase_2d / phase_to_height fringe(一般 N-step が既にあるので 4-step PSI は置かない)

ファミリ共通の入力契約(fail-closed)

全 op が入力を検証してから計算します。以下は 2026-09-01 の敵対監査で実際に見つかったバグか、それを塞ぐために書いた罠です:

代表的なパイプライン(op の繋がり)

検査機を 1 台、紙の上で設計しきる筋(検証済み examples/optics_imaging.py そのもの)。倍率 → 系の行列 → 深度 → 回折限界 → 収差 → 偏光、とデータ種が table → matrix → table → pairs → table で繋がります。

flowchart LR
    A[必要視野・センサ寸法] -->|倍率 m| B[thin_lens 焦点距離・物体距離]
    B -->|素子リスト table| C[abcd_matrix 系全体 2x2]
    C -->|matrix| D[abcd_trace 結像判定 B=0 / det=1]
    B --> E[depth_of_field 近・遠・過焦点]
    E -->|決めた絞り N| F[mtf_diffraction 回折限界 pairs]
    G[測った点像 image2d] --> H[psf_to_mtf 実測 MTF pairs]
    F --> I{実測は限界を超えていないか}
    H --> I
    J[match3d.fit_zernike 波面係数 table] --> K[wavefront_stats RMS/PV/Strehl]
    K --> I

偏光の筋は独立した 2 本立てで、互いの検算になります(同じ素子を両方の代数で組み、同じ Stokes ベクトルに落ちることを確かめる — 符号規約の取り違えはこの突き合わせでしか捕まりません)。

flowchart LR
    P[jones_element 偏光子/波長板/回転子] -->|cimage 2x2| Q[jones_apply]
    Q -->|jones| R[stokes_from_jones]
    S[mueller_element 同じ素子 + 減偏光子] -->|matrix 4x4| T[mueller_apply]
    R -->|stokes| U[stokes_analyze 偏光度・方位・楕円率]
    T -->|stokes| U
    R -.一致を検算.- T

使い方(最小の 1 本)

import optics as O

# 20 mm の視野を倍率 -0.43 で撮る 50 mm レンズ
lens = O.thin_lens(focal_mm=50.0, object_mm=50.0 * (1 - 1 / -0.43))
system = O.abcd_matrix([("free", lens["object_mm"]), ("lens", 50.0),
                        ("free", lens["image_mm"])])
print(O.abcd_trace(system)["imaging"])          # True = 本当に共役面

# 許容錯乱円 2 画素(3.45 µm 画素)での被写界深度と、その絞りの回折限界
dof = O.depth_of_field(50.0, 5.6, lens["object_mm"], coc_mm=2 * 3.45e-3)
mtf = O.mtf_diffraction(f_number=5.6, wavelength_um=0.55)
print(dof["depth_mm"], mtf[-1, 0])              # 深度 [mm], カットオフ [cyc/mm]

# 金属面のテカりを直交偏光で消す(Jones)
blocked = O.jones_apply(O.jones_element("polarizer", 90.0)
                        @ O.jones_element("polarizer", 0.0), [1.0, 0.0])
print(abs(blocked).max())                        # 0.0

設計(design) — 近軸の先を実光線で

上の 4 カテゴリは「設計の出発点」を閉形式で出します。実レンズがそこからどれだけずれるか — 像はどこに結び、どれだけボケ、どの面が原因で、製造ばらつきで歩留まりはどうなるか — は面を 1 枚ずつ実光線で通さないと分かりません。raytrace.py はそのための逐次光線追跡で、台帳では opsoptics の design カテゴリ(12 op)に載ります。全 op の共通入力は lens_system が返す検証済みの処方(table) です:

flowchart LR
    A[面リスト R/t/n/k/ap] --> B[lens_system 処方 table]
    G[glass n_d,V_d] --> A
    B --> C[paraxial_trace EFL/BFL/瞳/f値]
    B --> D[spot_stats / spot_diagram / ray_fan]
    B --> E[opd_map 波面 image2d]
    E --> F[wavefront_from_opd Zernike/Strehl]
    B --> H[seidel_coefficients 面ごとの S_I..S_V, C_L/C_T]
    B --> I[tolerance_analysis Monte-Carlo + 感度]
    C -.thick_lens 閉形式と 1e-9 で一致.- C

1. singlet を近軸 → スポット → Seidel と流す(検証済み examples/lens_design_demo.py の筋):

import raytrace as RT

lens = RT.lens_system()                          # 平凸 BK7 f=100 f/4、絞りは第 1 面
p = RT.paraxial_trace(lens)
print(round(p["efl"], 3), round(p["bfl"], 3), p["fno"])      # 100.0 96.704 4.0
spot = RT.spot_stats(lens)                       # 軸上の RMS スポット半径 [mm]
se = RT.seidel_coefficients(lens, field=5.0)     # 面ごとの三次収差
print(round(spot["rms_radius"], 4), round(se["waves"]["S_I"] / 8, 2))   # 0.1304 11.29
assert abs(p["efl"] - 100.0) < 1e-6 and se["waves"]["S_I"] > 0            # 補正不足(正)
assert abs(sum(r["S_I"] for r in se["per_surface"]) - se["total"]["S_I"]) < 1e-12

2. 放物面鏡は軸上で完全結像(OPD が 0)、同半径の球面鏡はそうでない:

import numpy as np
import raytrace as RT

para = RT.example_system("paraboloid")           # R=-200, k=-1, f=100
sph = RT.example_system("sphere_mirror")
opd = RT.opd_map(para, fill=np.nan)              # 射出瞳基準球に対する OPD [waves]
print(np.nanmax(np.abs(opd)) < 1e-6)             # True: 無収差(stigmatic)
print(RT.spot_stats(para)["rms_radius"] < 1e-9,  # True
      round(RT.spot_stats(sph)["rms_radius"], 3))  # 0.119(球面収差)
assert RT.spot_stats(para, field=1.0)["rms_radius"] > 0.01     # 軸外はコマ

3. 公差表 — どの面のどの誤差が効くか:

import raytrace as RT

tol = RT.tolerance_analysis(RT.lens_system(),
                            {"radius_pct": 0.5, "thickness_mm": 0.05,
                             "index": 0.001, "decenter_mm": 0.02, "tilt_deg": 0.05},
                            trials=20, seed=1)
r = tol["rms_spot"]
print(round(tol["nominal"]["rms_spot"], 4), round(r["p95"], 4), tol["failed"])
top = tol["sensitivity"][0]                      # |d_rms_spot| の大きい順
print(top["surface"], top["parameter"], round(top["d_efl"], 3))   # 0 R 0.5
assert tol["failed"] == 0 and top["parameter"] == "R"

棲み分け: optics = 近軸/波動(閉形式、設計の出発点)、raytrace = 実光線・設計(処方から数値で)。面での反射・屈折のスカラ公式(match3d.refract 等)を 1 本の光線に使うのは今までどおり match3d、系を通すのが raytrace です。

実硝材・非球面・色収差(2026-09-03 追加)

glass(nd, vd) は 2 項 Cauchy のモデルでした。glass_catalog("N-BK7") は Schott の公開 Sellmeier 定数(N-BK7 / N-K5 / N-BAK4 / N-SK16 / N-SSK5 / N-BAF10 / N-LAK22 / N-LAK9 / N-LASF9 / N-FK51A / N-F2 / N-SF2 / N-SF5 / N-SF10 / N-SF11 / N-SF6 / N-SF57)と Malitson の溶融石英・CaF2・サファイアの 20 種で、面の n に名前をそのまま書けます。sellmeier(B1, B2, B3, C1, C2, C3) は任意メーカーの定数から同じ形の table を作り、評価した nd / vd を返す(データシートと突き合わせられる)。20 種すべて、データシートの nd を 1e-4・vd を 0.3 以内で再現することをテストが固定しています。

面には k(円錐)に加えて asph=(A4, A6, A8, …)(偶数次非球面係数、mm⁻³, mm⁻⁵, …)を持たせられます。実光線は円錐面交点からの Newton 反復で非球面に載せ(残差 < 1e-8 mm)、法線はサグの勾配、Seidel には 4 次項 G = k c³/8 + A4 が入ります。example_system("asphere") は Descartes の無収差レンズ(物体側が平面、射出面が k = −n² の双曲面)で、スポット RMS が 3e-15 mm — 法線か交点が 1 桁でも違えば通らない検証です。

chromatic_shift(system, wavelengths) は処方を波長ごとに評価し直し、EFL / BFL / 主光線像高 / 基準像面でのスポットを返します(axial_color = F−C の BFL 差、lateral_color = 像高差、rms_polychromatic = 全波長プールの RMS)。example_system("catalog_doublet") は N-BK7 / N-SF2 の実硝材ダブレットです。

import raytrace as RT

g = RT.glass_catalog("N-SF2")
print(round(g["nd"], 5), round(g["vd"], 2))                 # 1.64769 33.82 (データシート)
single = RT.chromatic_shift(RT.lens_system())
doub = RT.chromatic_shift(RT.example_system("catalog_doublet"))
print(round(single["axial_color"], 3), round(doub["axial_color"], 3))   # -1.531 -0.299
assert abs(doub["axial_color"]) < abs(single["axial_color"]) / 4
print(round(RT.spot_stats(RT.example_system("asphere"))["rms_radius"], 12))   # 0.0 (stigmatic)

最適化(optimization) — 処方を減衰最小二乗で動かす

raytrace は処方を評価します。lensopt.optimize_lens(system, variables, fields, wavelengths, efl_target, …) はそれを変える側 — 設計者が毎日回している Levenberg–Marquardt(減衰最小二乗、Wynne & Wormell 1963 / Meiron 1965)です。残差ベクトルは各視野・各波長の瞳充填光線束の主光線からの横ずれ(mm、その RMS がスポット RMS)と、EFL を目標に固定する重い残差。変数は "R0"(曲率として動く — 平面を通り越せる)/ "t1" / "k1" / "A4_1"… の文字列。ヤコビアンは前進差分、正規方程式は λ·diag(JᵀJ) で減衰し、メリットが下がる歩幅だけ採用(λ を適応)。毎歩 lens_system で再検証するので、追跡器が拒否する処方(口径が円錐の範囲外、負の厚み)は返せません。merit_function は同じ残差の評価のみ、bend_singlet は Coddington 形状因子の閉形式(既定 = 球面収差最小の曲げ q = 2(n²−1)/(n+2))。

グラウンドトゥルース(tests/test_lensopt.py): 等凸 (q=0) から始めて EFL=100 固定で 2 半径を動かすと q=0.732(厚肉の総当たり最適 0.730、薄肉閉形式 0.740)に着地。平凸の射出面の k だけ動かすと Descartes の −n² を 1e-3 で再発見し RMS 0.57 mm → 1e-15 mm。A4 は円錐の 4 次項 k c³/8 に 1 % で収束し、A6/A8 を足すと単調に改善。

import lensopt
import raytrace as RT

eq = lensopt.bend_singlet(100.0, 1.5168, 2.0, 5.0, shape_factor=0.0)   # 等凸から
r = lensopt.optimize_lens(eq["system"], variables=["R0", "R1"], efl_target=100.0)
R1, R2 = r["variables"][0]["final"], r["variables"][1]["final"]
print(round((R2 + R1) / (R2 - R1), 3), round(r["efl_final"], 4), r["iterations"])   # 0.732 100.0 …
d = RT.example_system("catalog_doublet"); d["field"] = 3.0
m = lensopt.merit_function(d, wavelengths=[RT.WL_F, RT.WL_D, RT.WL_C])
print(sorted(m["rms_by_field"]), round(m["efl"], 3))
assert r["rms_final"] < r["rms_initial"] and abs(r["efl_final"] - 100.0) < 1e-3

照明設計(illumination) — 欠陥にコントラストがあるか、撮る前に決める

レンズは欠陥をどれだけ鮮明に写すかを決め、照明は欠陥にコントラストがあるかを決めます。illumdesign.py はその照明側を第一原理で数値化します:

import numpy as np
import illumdesign as ID

low = ID.light_source("ring", radius_mm=60.0, height_mm=15.0, n=24)     # 暗視野
u = ID.illumination_uniformity(ID.irradiance_map(low, size_mm=(50, 50), shape=(48, 48)))
dc = ID.defect_contrast(low, surface="glossy", slopes_deg=[5.0, 10.0])
print(round(u["uniformity"], 3), dc["regime"], [round(r["max_abs"], 3) for r in dc["per_slope"]])
sw = ID.lighting_sweep(surface="mirror", slope_deg=10.0, camera_height_mm=1000.0)
print(sw[int(np.argmax(sw[:, 1])), 0])                                  # 70.0 = 90 - 2*10
d = ID.illumination_design(surface="mirror", defect="topographic", slope_deg=10.0)
print(d["recommended"], d["rule_of_thumb"], d["agrees_with_rule"])   # coaxial coaxial True
assert d["ranking"][0]["score"] >= d["ranking"][-1]["score"]

正直な限界: 発光点は点光源の和(面光源は n を増やして近似)、多重反射・相互反射なし、BRDF は Lambert + 等方 GGX(異方性ヘアラインは未対応)、偏光なし、カメラは 1 点(視野内の視線変化は camera を変えて呼ぶ)。

結像シミュレーション(imaging_sim) — 設計したレンズで撮る

「擬似物理空間に光学系を組み、AI 学習用の欠陥画像を生成したい」という要望の出口が lensimage.py(台帳では opsoptics の imaging_sim カテゴリ、4 op)です。処方(lens_system の table)からセンサが記録する画像まで通します:

グラウンドトゥルース(tests/test_lensimage.py): 無収差瞳(singlet を 1 mm 絞り)は Airy — 第 1 暗環 1.22·λ·F# の 0.1 % 以内、暗環内エネルギー 83.8 %。f/4 singlet(球面収差 11 波)の Strehl 0.011。放物面鏡の歪曲 < 1e-7 %、singlet は樽型 −0.065 %(15 deg)。δ 画像を放物面鏡で描くと画素積分 Airy PSF と 1e-12 で一致。傷マスクの歪曲後 IoU 0.77。

1. PSF・Strehl・歪曲を singlet と doublet で見比べる(検証済み examples/lens_defect_dataset_demo.py の筋):

import raytrace as RT
import lensimage

sg, db = RT.example_system("singlet"), RT.example_system("doublet")
psf = lensimage.psf_from_opd(sg, pixel_pitch_um=5.5)      # 画素積分 PSF、和 = 1
print(psf.shape, round(float(psf.sum()), 6))
g = lensimage.psf_field_grid(db, fields=(0.0, 4.0))
print([round(s, 4) for s in g["strehl"]], [round(r, 4) for r in g["rms_spot_mm"]])
d = lensimage.distortion_map(sg, fields=[0.0, 5.0, 10.0, 15.0])
print([round(v, 4) for v in d["distortion_pct"]])          # 0, ..., -0.0649 (樽型)
assert d["distortion_pct"][-1] < 0.0

2. 欠陥画像を doublet 越しに 2 枚作り、注釈が像に揃っていることを確かめる:

import numpy as np
import raytrace as RT
import lensimage

recs = lensimage.defect_dataset(2, system=RT.example_system("doublet"),
                                size=(96, 96), seed=1)  # out_dir=... で PNG + annotations.json
r = recs[0]
print(r["image"].shape, r["mask"].dtype, [(d["kind"], d["bbox"]) for d in r["defects"]])
print({k: round(v, 4) for k, v in r["lens"].items()})   # efl / fno / rms_spot / 歪曲
assert r["image"].shape == r["mask"].shape and r["mask"].any()
same = lensimage.defect_dataset(2, system=RT.example_system("doublet"), size=(96, 96), seed=1)
assert np.array_equal(r["image"], same[0]["image"])       # seed で再現

正直な限界: 単色・インコヒーレント結像(色ごとの PSF が要るなら波長を変えて呼び足し合わせる)。視野依存 PSF は zones² 点の線形補間で、軸外 PSF は +y 視野の PSF をタイル方位へ回した近似(回転対称な処方が前提 — 偏心・傾きのある処方は +y の PSF のみ)。周辺光量は光線の到達数 × cos⁴ で、瞳収差の重みは入っていません。センサは線形・一様の雑音モデルのみ(クロストーク・PRNU・カラーフィルタなし)。

アルゴリズムの正典(著者・年)

正直な限界

© 2026 Kazufumi Furuse — Fullseye operator documentation. Licensed under Apache-2.0.

This site is open source. Improve this page.