はじめに
サイバネットシステム株式会社 光学エンジニアの高橋です。
主に光学設計ソフトウェア「Ansys Zemax OpticStudio」(以下、Zemax)の
販売サポートや受託設計・解析を担当しております。
三次元測定機の実測値をZemaxデータに反映したい
以前、レンズメーカーでカメラレンズの設計開発を行っていた頃、
超高精度三次元測定機「UA3P」などに代表される
形状測定機で得られた生の点群データを取り扱う場面が良くあったのですが
- Excelなどで原点オフセットを調整
- 任意断面を抽出（x–z）,（y–z）
- そのデータを ZEMAX の メリットファンクション に入力
- Zemaxの最適化機能で R, k, 非球面係数 を最適化し性能の劣化度合をシミュレーション
…という何とも煩わしい作業を行っていました。
毎回手作業が多く、断面抽出・データ整形・パラメータ入力の繰り返しはかなり煩雑です。
そこで「測定データからZemax非球面係数へのフィットを自動化できないか」と考えて作ってみたものが、今回ご紹介するスクリプトです。
偶数非球面のモデル式
フィッティング対象はZemaxの面タイプ「Even Asphere」です。
通常、レンズ面の子午線方向の高さ z(x)は以下の式で表されます。
z(r) = \frac{c r^2}{1 + \sqrt{\,1 - (1+k)c^2 r^2\,}}
+ \sum_{i=2}^{n} \alpha_{2i} r^{2i}
\ , \ c = \frac{1}{R}
𝑅：曲率半径
𝑘：コーニック定数
α𝑖：非球面係数
上式が偶数次非球面の一般式といえます。
実際にZemaxで面タイプに「Even Asphere」を選択した場合は
16次の項までが有効になるので以下式になります。
z(r) = \frac{c r^2}{1 + \sqrt{\,1 - (1+k)c^2 r^2\,}}
+ \alpha_4 r^4 + \alpha_6 r^6 + \alpha_8 r^8 + \alpha_{10} r^{10}
+ \alpha_{12} r^{12} + \alpha_{14} r^{14} + \alpha_{16} r^{16}
どうフィッティングするか
光学設計に慣れている方なら感覚的に理解されると思いますが、
いきなりR, k, αのすべてを同時に変数として最適化してしまうと
以下の問題が起こりやすくなります。
- 収束性が悪化し、局所解に落ち込む
- パラメータ間の強い相関により、物理的に不自然な値に引きずられる
- 高次項の𝛼𝑖が低次形状のズレを吸収してしまい、本来の半径・コーニックの寄与バランスが失われる（非球面の設計を行う上でこれが一番厄介）
そこで本スクリプトを組む場合も、段階的に変数を開放するフィッティングを採用します。
具体的には、
- まずは半径 𝑅 のみをフィットし、大まかな曲率を合わせる
- 1の結果を初期値にしてフル（R, k, α）でフィッティング
細かいことを言うと、1と2の間に「低次の項までの最適化」を間に挟んで
最終的に高次を含めたフルフィットの方が良さそう、とか
2の最適化時のRは1で決定した数値で固定で良いのでは、とか色々ありますが
一旦これで進めます。
スクリプトの要点だけ紹介
1. Rのみの初期合わせ（マルチスタート）
広い初期値レンジ×繰り返し最適化で安定化を図ります。
def fit_R_only_multistart(x_data, z_data, n_runs=10,
k_fixed=0.0, R_min=-50.0, R_max=50.0):
best_R, best_rms = None, None
for _ in range(n_runs):
R_init = np.random.uniform(R_min, R_max)
def model_for_fit(x, R): # kは固定してRだけ
return asphere_simple(x, R, fixed_k=k_fixed)
try:
popt, _ = curve_fit(model_for_fit, x_data, z_data, p0=[R_init])
except RuntimeError:
continue
R_tmp = popt[0]
rms = np.sqrt(np.mean((z_data - model_for_fit(x_data, R_tmp))**2))
if best_rms is None or rms < best_rms:
best_R, best_rms = R_tmp, rms
return best_R, best_rms
2. R近傍でフルモデル最適化（制約付き）
1の処理で得られたRの周辺に探索域を絞り、かつ各係数にはある程度の数値制約を設け妥当性を担保
def fit_full_multistart_nearR(x_data, z_data, R_base,
n_runs=10, R_spread=0.5,
k_min=-10.0, k_max=10.0):
best_params, best_rms = None, None
for _ in range(n_runs):
R_init = np.random.uniform(R_base - R_spread, R_base + R_spread)
k_init = np.random.uniform(k_min, k_max)
alpha_init = [np.random.uniform(-1e-1, 1e-1) for _ in range(7)]
p0 = [R_init, k_init] + alpha_init
lower = [R_base - R_spread, k_min] + [-1e-1]*7
upper = [R_base + R_spread, k_max] + [ 1e-1]*7
try:
popt, _ = curve_fit(asphere_full, x_data, z_data, p0=p0,
bounds=[lower, upper], maxfev=10000)
except RuntimeError:
continue
rms = np.sqrt(np.mean((z_data - asphere_full(x_data, *popt))**2))
if best_rms is None or rms < best_rms:
best_params, best_rms = popt, rms
return best_params, best_rms
3.R1/R2の符号合わせ（測定系との整合）
R1なら z→−z にして符号規約を合わせます。
測定面がR1面のときは測定機とZemaxで符号が逆転します（よくハマるポイント）。
x_data, z_data = data[:,0], data[:,1]
if surface_type == "R1": # 測定定義に合わせて符号反転
z_data = -z_data
4.ZOS-APIでLDE(レンズデータエディタ)へ自動反映
Radius/Conic/Par2..Par8 に R, k, α4..α16 を流し込み、別名保存。
surf = TheLDE.GetSurfaceAt(surface_number)
surf.Radius = best_popt[0] # R
surf.Conic = best_popt[1] # k
for i, col in enumerate([Par2, Par3, Par4, Par5, Par6, Par7, Par8], start=2):
surf.GetSurfaceCell(col).DoubleValue = best_popt[i]
TheSystem.SaveAs(zemax_file.replace(".zmx", "_fitted.zmx"))
GUIで“誰でも回せる”運用
あとこだわりポイントとしては、Tkinterのモジュールを使って簡易GUIを付けてみました。
実行すると下記のようなポップ画面が現れ、
入力データと測定面（R1/R2）の選択ができます。
適用させたいレンズデータの何番目の面が対象かを間違わずに入力することも重要です。
これであれば現場メンバーにツール配布しても運用しやすいのではないかと。
実行結果イメージ
下図のようなフィッテイング画面がまず現れます。
フィッテイング精度をここで確認します。
Zemaxのレンズデータが自動で生成されます。
面10にフィッテイングされた係数が適用されています（見づらくてすみません）。
問い合わせ
本記事は要点のみの紹介でした。
スクリプト全文／導入支援／カスタマイズ（データ前処理やロバスト化、他面タイプ対応など）をご希望の方は、サイバネットシステム光学エンジニアリングサービスまでお問い合わせください。
本件とは関係なく、光学設計・解析・光学教育などのご依頼もぜひ！
※本記事は筆者個人の見解であり、所属組織の公式見解を示すものではありません。