JP2024066688A - 収差推定方法、収差推定装置、演算装置、プログラム、および記憶媒体 - Google Patents

Abstract

【課題】広い波長幅を有する光源を用いても、高速かつ高精度で被検光学系の収差を推定することが可能な収差推定方法を提供する。【解決手段】収差推定方法は、互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、複数のデフォーカス位置に対応する複数の第１光強度分布を取得する計測ステップ（Ｓ１）と、複数のデフォーカス位置に対応する複数の第２光強度分布を算出する算出ステップ（Ｓ３２）と、複数の第１光強度分布または複数の第２光強度分布を合成する合成ステップ（Ｓ３３、Ｓ３７）とを有し、合成ステップの結果に基づいて被検光学系の収差を推定する。【選択図】図１

Description

本発明は、収差推定方法、収差推定装置、演算装置、プログラム、および記憶媒体に関する。

従来、レーザなどの単色光源を用いて複数のデフォーカス位置で光強度分布を計測し、その計測結果を最も良く再現する収差を最適化演算により推定する収差推定方法が知られている。しかし、レーザなどの単色光源は、他の光源に比べて高価であるとともに管理コストが大きい。また、計測条件によっては、レーザなどの単色光源を利用できない場合がある。このため、レーザなどの単色光源に代えて、ＬＥＤまたは白熱灯などの波長幅を有する光源を用いた収差推定手法が望まれる。

ところで、最適化演算を用いた収差推定手法では、推定された波面収差から像面での光強度分布を算出し、計測した光強度分布と一致するかを判定する。このため、波長幅を有する光源を用いる場合など、波面収差から光強度分布を算出する算出条件が計測条件と異なると、推定結果に誤差が生じてしまう。

非特許文献１には、波面収差から光強度分布を算出する際に、各波長での光強度分布を算出し、それらの光高度分布を足し合わせることで計測条件と同じ計算条件を再現する方法が開示されている。非特許文献２には、単色で計算された光強度分布にボケカーネルを２次元平面内で畳み込むことで、計算条件を計測条件に近づける方法が開示されている。特許文献１には、計測された光強度分布に対してボケカーネルを使った画像の先鋭化処理を加えることで、計測された光強度分布を単色での計算結果に近づける方法が開示されている。

国際公開第２０１９／１４２３１３号

Ｊａｍｅｓ． Ｒ． Ｆｉｅｎｕｐ， "Ｐｈａｓｅ ｒｅｔｒｉｅｖａｌ ｆｏｒｕｎｄｅｒｓａｍｐｌｅｄ ｂｒｏａｄｂａｎｄ ｉｍａｇｅｓ" Ｊｏｕｒｎａｌ ｏｆ ｔｈｅ Ｏｐｔｉｃａｌ Ｓｏｃｉｅｔｙ ｏｆ Ａｍｅｒｉｃａ Ａ Ｖｏｌ．１６， ｐｐ．１８３１， Ｊｕｌｙ １９９９， ＵＳＡ． Ａｌｄｅｎ． Ｓ． Ｊｕｒｌｉｎｇ ａｎｄ Ｊａｍｅｓ． Ｒ．Ｆｉｅｎｕｐ， "Ａ Ｆａｓｔ Ａｐｐｒｏｘｉｍａｔｉｏｎ Ｍｅｔｈｏｄ ｆｏｒ Ｂｒｏａｄｂａｎｄ Ｐｈａｓｅ Ｒｅｔｒｉｅｖａｌ" ｉｎ Ｉｍａｇｉｎｇ ａｎｄ Ａｐｐｌｉｅｄ Ｏｐｔｉｃｓ， ＯＳＡ Ｔｅｃｈｎｉｃａｌ Ｄｉｇｅｓｔ （ＣＤ） （Ｏｐｔｉｃａｌ Ｓｏｃｉｅｔｙ ｏｆ Ａｍｅｒｉｃａ， Ｊｕｌｙ ２０１１）， ｐａｐｅｒ ＳＭＣ３．

非特許文献１に開示されている方法では、波長幅の広がりに合わせて複数の波長で光強度分布を算出する必要があるため、計算時間が増大する。非特許文献２に開示されている方法では、高速な演算が可能であるが、波長幅の広がりにより生じる光強度分布の変化は単なる２次元平面上に広がるボケとは異なるため、高精度で収差を推定することはできない。特に、波長によって焦点位置が異なる軸上色収差は、２次元平面内でのボケでは再現できないため、被検光学系に軸上色収差が存在する場合、高精度で収差を推定することはできない。同様に、特許文献１に開示されている方法では、軸上色収差による影響はボケカーネルよる画像の先鋭化処理では除去できないため、高精度で収差を推定することはできない。

そこで本発明は、広い波長幅を有する光源を用いても、高速かつ高精度で被検光学系の収差を推定することが可能な収差推定方法を提供することを目的とする。

本発明の一側面としての収差推定方法は、互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得する計測ステップと、前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出する算出ステップと、前記複数の第１光強度分布または前記複数の第２光強度分布を合成する合成ステップとを有し、前記合成ステップの結果に基づいて被検光学系の収差を推定する。

本発明の他の目的及び特徴は、以下の実施形態において説明される。

本発明によれば、広い波長幅を有する光源を用いても、高速かつ高精度で被検光学系の収差を推定することが可能な収差推定方法を提供することができる。

各実施例における収差測定装置の概略図である。 各実施例における収差推定方法を示すフローチャートである。 実施例１における収差推定方法を示すフローチャートである。 実施例１における推定結果を示す図である。 実施例２における収差推定方法を示すフローチャートである。

以下、本発明の実施形態について、図面を参照しながら詳細に説明する。各図において、同一の部材については同一の参照符号を付し、重複する説明は省略する。

まず、図１を参照して、本実施形態における収差推定装置１００について説明する。図１は、収差推定装置１００の概略図である。光源１０１から発せられた光は、被写体１０２を照明する。被写体１０２を通過した光は、被検光学系１０３によって撮像素子１０４上に結像され、光学像を形成する。撮像素子１０４は、駆動装置（電動ステージ）１０５上に設置されている。駆動装置１０５は、コンピュータ（制御部）１０６によって制御されており、撮像素子１０４を指定の位置まで矢印に示される光軸に沿った方向（光軸方向）に移動させる。撮像素子１０４は、移動した位置において光学像の光強度分布を取得し、取得した光強度分布をコンピュータ１０６または不図示のデータ保持装置に保存する。コンピュータ１０６は、取得された光強度分布に対して後処理を実行することで、被検光学系１０３の収差を推定する。このときコンピュータ１０６は、後述の取得部と算出部と合成部と推定部とを有する演算装置として機能する。後処理は、コンピュータ１０６により実行されるか、または別の演算装置により実行されてもよい。また、ネットワークを通じてクラウド上に存在する演算装置が後処理を実行してもよい。表示部１０７は、推定された収差を表示する。

光源１０１は、例えば、発光スペクトルに幅（波長幅）を有する（広い波長幅を有する）光源である。広い波長幅を有する光源は、例えば、ＬＥＤ、蛍光灯、または白熱灯であるが、これらに限定されるものではない。熱源を用意し、熱源から発せられる赤外光を照明光として用いてもよい。また光源１０１は、バンドパスフィルタ等の波長を制限する素子を備えていてもよい。例えば、光源１０１としてバンドパスフィルタを備えたハロゲンランプを用いることで、所望の波長幅を有する光を照明光とすることができる。また、生体観察等に用いられる蛍光染料から発せられる蛍光を利用してもよい。

被写体１０２は、微小な開口を有する部材であり、例えばピンホールであるが、これに限定されるものではない。また、光ファイバの端部を被写体１０２として設置し、反対側の端部を光源１０１によって照明すれば同様の効果が得られる。被写体１０２は、ピンホールのような透過部によって光を空間的に制限する部材に限定されるものではなく、散乱部または発光部によって光の発生領域を限定する部材でもよい。例えば、被写体１０２として黒体の平面基盤上に描かれた微小な白円や微小な発光塗料を用いてもよい。微小な透過領域または微小な散乱領域、または微小な発光領域を有するものであれば、被写体１０２として機能する。また、光源１０１が被写体１０２を兼ねてもよい。例えば、被検光学系１０３として顕微鏡等の拡大光学系を用いる場合、微小な蛍光色素からの光を撮像素子１０４上に結像してもよい。この場合、微小な蛍光色素が光源１０１となり、被写体１０２を兼ねる。

被検光学系１０３は、計測対象となる光学系（単独の光学系）、または計測対象となる光学系と別の光学系との組み合わせでもよい。例えば、計測対象がカメラレンズまたは望遠鏡である場合、被検光学系１０３は、カメラレンズまたは望遠鏡そのものである。また、計測対象が顕微鏡の対物レンズである場合、被検光学系１０３は、顕微鏡対物レンズと結像レンズとの組み合わせでもよい。計測対象の光学系を含み、被写体１０２の光学像を形成できる光学系であれば、被検光学系１０３として用いることができる。また、天体望遠鏡または顕微鏡等で用いられる補償光学では、大気や生体試料によって収差を計測対象としている。これらの収差も本実施形態により計測可能であるため、被検光学系１０３は、これらの収差を発生させる大気や生体試料を含んでいてもよい。

撮像素子１０４は、ＣＣＤセンサまたはＣＭＯＳセンサ等の撮像センサであり、光電変換素子（画素）がアレイ状に配置されたデバイスである。駆動装置１０５は、電動ステージ等の１次元方向に駆動可能な装置を用いればよい。また、駆動装置１０５が動かす対象は、撮像素子１０４に限定されるものではない。駆動装置１０５は、被検光学系１０３を駆動してもよく、または被写体１０２を駆動してもよい。また、被検光学系１０３がフォーカス調整機構を有する場合、フォーカス調整機構と通信可能な制御部を用いてフォーカス調整機構を駆動させてもよい。いずれの方法においても、光学像に与えられるデフォーカスを変化させる機構を有していればよい。本実施形態では、駆動装置１０５が撮像素子１０４を動かすことによってデフォーカス量を変えるとして説明する。このため、光軸方向の位置とデフォーカス位置は同義である。また、被検光学系１０３によって形成される被写体１０２の結像位置を原点とし、それを基準に撮像素子１０４が置かれる位置を定める。

以下、本実施形態が機能する原理と、その原理に基づいた収差推定方法を順に説明する。まず、本実施形態が機能する原理を説明する。以下、算出される収差は、被検光学系１０３の収差（波面収差）Ｗ（ξ、η）として説明を行う。（ξ，η）は瞳の半径で規格化された瞳空間での直交座標である。本実施形態によって計測が可能な収差はこれに限定されるものではない。波面収差Ｗ（ξ、η）が得られれば、簡単な演算によって横収差量や縦収差量を算出可能である。収差Ｗ（ξ、η）をＺｅｒｎｉｋｅ多項式で展開することで、ザイデル収差に変換可能である。

計測された光強度分布からコンピュータ１０６によって被検光学系１０３の収差Ｗ（ξ、η）を推定する方法として、最適化演算が例示される。最適化演算では、以下の式（１）で表される目的関数Ｆ（Ｗ）を最小化するＷ（ξ、η）が繰り返し演算によって探索される。

ここで、Ｉ_ｓ（ｘ，ｙ，ｚ；Ｗ）は収差Ｗ（ξ、η）から算出される光強度分布、Ｉ_ｍ（ｘ，ｙ，ｚ）は計測された光強度分布、ｘ、ｙは光軸と垂直な平面での２次元直交座標、ｚは光軸方向の座標、ｊは計測を行ったｚの位置番号、Ｊは計測の総数である。目的関数は種々の形が存在する。式（１）に示されるＩ_ｓとＩ_ｍの差分二乗和以外にも、振幅分布（光強度分布の平方根）の差分二乗和等がある。問題に応じて目的関数Ｆ（Ｗ）は適宜選択すればよい。

目的関数Ｆ（Ｗ）を最小化する方法は種々存在し、最急降下法、共役勾配法、準ニュートン法等がある。また目的関数を陽に計算せずとも、特許文献１に記載の方法のようにフーリエ変換を繰り返すだけでも最適化は実現される。本実施形態では、説明を容易にするため、目的関数Ｆ（Ｗ）を用いた最適化を示す。

式（１）で表される目的関数Ｆ（Ｗ）を算出するには、収差Ｗ（ξ、η）から光強度分布Ｉ_ｓ（ｘ，ｙ，ｚ；Ｗ）を算出する必要がある。波長λの単色光に対して収差Ｗ（ξ、η）から光強度分布Ｉ’（ｘ，ｙ，ｚ；Ｗ）を算出する方法は数２となる。

ＦＴはフーリエ変換を表す演算子である。Ｔ（ξ、η）は瞳の透過率分布を表す関数で、通常の光学系では円形開口となる。ＮＡは被検光学系１０３の像側開口数である。一般に、式（２）の左辺と右辺は比例関係であるが、説明を簡単にするため比例係数を省略し、等号で表記している。

式（２）は単色での計算式であるため、光源１０１が波長幅を持つ場合、Ｉ_ｓの計算に式（２）を用いると、計算条件とＩ_ｍを計測した計測条件に不一致が生じてしまう。式（１）で示されるように、最適化演算ではＩ_ｓとＩ_ｍの差分値を小さくするＷが探索されるため、高精度な収差を推定するには、収差Ｗ（ξ、η）から光強度分布Ｉ_ｓを計算する過程に波長幅の影響を補正する処理が必要となる。光源１０１が持つ波長幅によって生じる不一致の大きな要因の１つが被検光学系１０３の色収差である。特に、波面収差は軸上で計測されることが多いため、軸上色収差が不一致の主要因となる。そこで本実施形態では、軸上色収差による影響を高速、高精度に付与もしくは除去する処理を加えることで、波長幅のある光源を用いた場合でも高速、高精度な収差の推定を実現する。

結像理論から波長幅がある光源１０１によって形成される光強度分布Ｉ（ｘ，ｙ，ｚ；Ｗ）は、単色での光強度分布Ｉ’（ｘ，ｙ，ｚ；Ｗ；λ）を足し合わせたものとなり、以下の式（３）で表される。

ここで、Ｓ（λ）は光源１０１の発光スペクトルおよび撮像素子１０４の感度特性によって定まる。場合によっては、被検光学系１０３の透過率スペクトルも考慮される。λ_１及びλ_２はＳ（λ）が値を持つλの下端と上端である。式（３）を用いれば、波長幅のある光源１０１を用いた場合の光強度分布Ｉを算出できるが、波長による積分は演算負荷が高い。そこで、式（３）を精度良く近似する。

軸上色収差とは、波長によって焦点位置が異なる収差である。この収差の影響を、焦点位置の違いのみに注目し、他の影響を無視する。すなわち、波長λにおける光強度分布Ｉ’（ｘ，ｙ，ｚ；Ｗ；λ）は、波長λはλ_０で一定のまま焦点位置のみを変化させた光強度分布とみなし、以下の式（４）で表される。

また、波長の違いによる焦点位置の違いは、以下の式（５）で表される。

ここで、ｚ’は色収差による焦点位置のずれ量、αは焦点位置のずれ量と波長の比例係数である。式（４）を式（３）に代入すると、以下の式（６）が得られる。

ここで、ｚ_１＝α（λ_１－λ_０）、ｚ_２＝α（λ_２－λ_０）である。式（６）により波長幅がある場合の光強度分布が単色での演算に変換されるが、ｚ方向の積分を実行するために積分領域におけるＩ’（ｘ，ｙ，ｚ；Ｗ；λ_０）をすべて演算する必要があり、時間がかかる。そこで、積分領域内のすべてのｚに対してＩ’（ｘ，ｙ，ｚ；Ｗ；λ_０）を算出するのではなく、少数の位置ｚ_ｊにおけるＩ’（ｘ，ｙ，ｚ_ｊ；Ｗ；λ_０）を算出し、他の位置でのＩ’（ｘ，ｙ，ｚ；Ｗ；λ_０）は補間演算によって代替すると考える。説明を簡単にするため、まずは２か所の位置ｚ_１及びｚ_２で算出を行い、他のｚでのＩ’（ｘ，ｙ，ｚ；Ｗ；λ_０）は補間演算で代用したとする。最も単純な線形補間を行ったとすると、Ｉ’（ｘ，ｙ，ｚ；Ｗ；λ_０）は、以下の式（７）で表される。

式（７）におけるＩ’（ｘ，ｙ，ｚ；Ｗ；λ_０）は、ｚ方向に対してｚの１次関数となるため、数６の積分は単なるｚの積分となる。これは、解析的に実行可能なため高速な処理が可能となる。

本演算を一般化する。ｚ方向にＮ個の離散的な位置ｚ_ｎ（ｎ＝１，２，，，Ｎ）で光強度分布Ｉ’（ｘ，ｙ，ｚ_ｎ；Ｗ；λ_０）が既知の時、位置ｚ（ｚ_１≦ｚ≦ｚ_Ｎ）における光強度分布Ｉ’（ｘ，ｙ，ｚ；Ｗ；λ_０）は、補間演算によって、以下の式（８）で表される。

ここで、ｆ_ｎ（ｚ）は適当な基底関数、ａ_ｎ（ｘ、ｙ；Ｗ；λ_０）はＩ’（ｘ，ｙ，ｚ_ｎ；Ｗ；λ_０）から線形変換によって定まる（ｘ、ｙ）における展開係数である。ｆ_ｎ（ｚ）としてはｚの多項式やフーリエ級数の基底関数を用いることができる。ｆ_ｎ（ｚ）としてｚの多項式を取り、Ｎ＝２とすれば、式（８）は式（７）に帰着する。また、フーリエ級数の基底関数をｆ_ｎ（ｚ）としてとれば、ａ_ｎ（ｘ、ｙ；Ｗ；λ_０）は、Ｉ’（ｘ，ｙ，ｚ_ｊ；Ｗ；λ_０）を各ｘ、ｙ座標においてｚ方向に逆フーリエ変換した時のフーリエ係数となる。

ここで、光強度分布Ｉ’（ｘ，ｙ，ｚ；Ｗ；λ_０）の各ｘ、ｙ座標におけるｚ方向の振る舞いを考える。Ｉ’（ｘ、ｙ、ｚ；Ｗ；λ_０）は、被検光学系１０３によって形成される像面での光強度分布であるため、ｚ方向の振る舞いは焦点深度λ_０／ＮＡ^２程度の細かさの変化しか存在しない。換言すれば、光強度分布Ｉ’（ｘ，ｙ，ｚ；Ｗ；λ_０）は、ｚ方向に周波数空間で帯域制限された関数であるため、式（８）は以下の式（９）のように表される。

ここで、Δｚはｚ_ｎの間隔である。式（９）を式（６）に代入すると、以下の式（１０）および式（１１）が得られる。

本実施形態の収差推定では、Ｉ（ｘ、ｙ、ｚ；Ｗ）として興味があるのは計測を行った位置ｚ_ｊでの分布である。従って、ｚ＝ｚ_ｊとすると、式（１０）は以下の式（１２）のように表される。

実際の演算においては、一連の計測位置ｚ_ｊ（ｊ＝１，２，，，Ｊ）とｚ_ｎ（ｎ＝１，２，，，Ｎ）を異なる値にする必要はないため、ｚ_ｊとｚ_ｎが同じ値を持つとして構わない。その上で、式（１２）を行列の形式で書き直すと、以下の式（１３）および式（１４）のように表される。

式（１３）および式（１４）より、波長幅がある光源１０１を使って得られる光強度分布Ｉ（ｘ、ｙ、ｚ_ｊ）は単色で離散的な複数の位置で得られた光強度分布Ｉ’（ｘ、ｙ、ｚ_ｊ；Ｗ；λ_０）の線形結合、すなわち重み付き和で算出できることがわかる。より広義には、Ｉ（ｘ、ｙ、ｚ_ｊ）は複数の光強度分布Ｉ’（ｘ、ｙ、ｚ_ｊ；Ｗ；λ_０）から合成処理によって算出できる。このようにして算出された光強度分布Ｉを数１のＩ_ｓとして用いることで、波長幅がある光源１０１を使った場合においても計算条件と計測条件を一致させることができ、高精度な収差推定が可能となる。

特に、式（６）または式（１２）からわかるように、軸上色収差によって生じる光強度分布の変化は２次元平面内での光強度分布のボケ広がりでは無く、光軸方向への光強度分布の重畳と考えることができる。このことから、２次元平面内でのボケカーネルを使った畳み込み処理よりも、本実施形態の式（１３）に示される光軸方向に異なる位置で算出された光強度分布Ｉ’の重み付き和を補正処理として追加した方が高精度な収差推定には有効となる。また、式（１３）による変換は高々Ｎ×Ｎの行列による演算のため、高速な処理が可能となる。加えて、β及びＢはＷに依らない単なる１次元の積分演算で得られるため、一度計算すれば最適化演算の最中に再計算する必要は無い。

次に、図２を参照して、以上の原理に基づく収差推定方法を説明する。図２は、収差推定方法を示すフローチャートである。

まずステップＳ１において、収差推定装置１００は、複数のデフォーカス位置における光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）を計測する（計測ステップ、第１取得ステップ）。すなわち収差推定装置１００は、互いに異なる複数のデフォーカス位置で被検光学系１０３を介して物体を撮像する（デフォーカス位置に撮像素子の撮像面を配置（移動）しつつ被検光学系を介して物体を撮像する）。そしてコンピュータ１０６などの演算装置（取得部）は、複数のデフォーカス位置に対応する複数の第１光強度分布を取得する。この計測は、駆動装置１０５による駆動および撮像素子１０４による光強度分布の取得をコンピュータ１０６による制御によって繰り返すことで実行される。取得された光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）は、データ保持装置に保存されるか、またはコンピュータ１０６の一時メモリ（記憶部）に保存される。続いてステップＳ２において、コンピュータ１０６は、取得された光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）を読み込む。なお、ステップＳ１にて光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）がコンピュータ１０６の一時メモリに保存されている場合、本ステップはスキップできる。

続いてステップＳ３において、コンピュータ１０６は、光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）に基づいて後処理を実行し、被検光学系１０３の収差を推定する。推定された収差は表示部１０７に送られて表示される。ここで、図３を参照して、ステップＳ３にてコンピュータ１０６が実行する収差推定方法を説明する。図３は、収差推定方法を示すフローチャートである。

まずステップＳ３１において、コンピュータ１０６は、式（１１）に基づいて重み係数βを算出する。式（１１）において、Δｚは、ステップＳ１にて計測を行ったデフォーカス位置の間隔を用いればよい。Ｓ（ｚ）は、光源１０１の発光スペクトルおよび撮像素子１０４の感度の波長依存性の積を式（５）によって変数変換した関数を用いればよい。また、式（１１）における積分は、被検光学系１０３の色収差によって定まる範囲で実行されていればよい。

続いてステップＳ３２において、コンピュータ（算出部）１０６は、収差（推定収差）Ｗから計測を行った位置ｚ_ｊにおける中心波長λ_０に対する単色での光強度分布Ｉ’（ｘ、ｙ、ｚ_ｊ；Ｗ；λ_０）を算出する（算出ステップ、第２取得ステップ）。すなわちコンピュータ１０６は、複数のデフォーカス位置に対応する複数の第２光強度分布を算出する。繰り返し演算の最初はＷとして初期値を設定する。初期値は、無収差を設定してもよいし、別の演算手段で近似的に求めた収差を利用してもよい。また、設計値等から想定される収差を利用することもできる。収差Ｗおよび中心波長λ_０から光強度分布Ｉ’（ｘ、ｙ、ｚ_ｊ；Ｗ；λ_０）を算出する方法としては、式（２）を用いることができる。撮像素子１０４の画素によって光強度分布が離散化される影響を加味した演算を加えることもできる。

続いてステップＳ３３において、コンピュータ（合成部）１０６は、ステップＳ３１にて算出された重みβおよびステップＳ３２にて算出された光強度分布Ｉ’（ｘ、ｙ、ｚ_ｊ；Ｗ；λ_０）から、式（１３）または式（１４）に基づいて合成処理を行う。すなわちコンピュータ１０６は、重みβを用いて複数の第２光強度分布を合成する（合成ステップ）。これにより、色収差の影響が考慮された光強度分布Ｉ（ｘ、ｙ、ｚ_ｊ；Ｗ）を算出することができる。

続いてステップＳ３４において、コンピュータ１０６は、ステップＳ３３にて算出されたＩ（ｘ、ｙ、ｚ_ｊ；Ｗ）をＩ_ｓ（ｘ、ｙ、ｚ_ｊ；Ｗ）として、式（１）の目的関数Ｆ（Ｗ）を算出する。続いてステップＳ３５において、コンピュータ１０６は、算出された目的関数Ｆ（Ｗ）を所定値と比較する。コンピュータ１０６は、目的関数Ｆ（Ｗ）が所定値以下であると判定した場合、演算を終了して収差Ｗを出力する。一方、コンピュータ１０６は、目的関数Ｆ（Ｗ）が所定値を上回っていた場合、ステップＳ３６に進む。

ステップＳ３６において、コンピュータ１０６は、収差Ｗを更新する。このときの更新方法は種々存在するが、目的関数Ｆ（Ｗ）を収差Ｗに対して微分した式に基づいて更新する勾配法を用いればよい。コンピュータ（推定部）１０６は、目的関数Ｆ（Ｗ）が所定値を下回るまで、または所定の繰り返し回数が実行されるまで、ステップＳ３２～Ｓ３６までを繰り返し実行する。以上のフローによって、被検光学系１０３の収差Ｗが推定される。

本実施形態は、数学的にモデル化することができるため、コンピュータ・システムのソフトウェア機能として実装可能である。ここで、コンピュータ・システムのソフトウェア機能は、実行可能なコードを含んだプログラミング（プログラム）を含む。ソフトウェア・コードは、汎用コンピュータで実行可能である。ソフトウェア・コード動作中に、コード、または関連データ記録は、汎用コンピュータ・プラットフォーム内に格納される。しかしながら、その他の場合、ソフトウェアは他の場所に格納される、または適切な汎用コンピュータ・システムにロードされる。したがって、ソフトウェア・コードは、１つまたは複数のモジュールとして、少なくとも１つの機械可読媒体（記憶媒体）で保持可能である。

以下、本実施形態の好ましい実施例を詳細に説明する。

まず、実施例１における収差推定方法について、シミュレーションによる解析結果を用いて説明する。本実施例は、図１の収差推定装置１００によって実現される。光源１０１は、白色照明部とバンドパスフィルタとを有する。バンドパスフィルタは、５００ｎｍから５５０ｎｍの光を透過する。被写体１０２は、直径３μｍの穴が開いたピンホールである。被検光学系１０３は、等倍でＮＡが０．１７９のレンズである。レンズは波面収差を有し、Ｆｒｉｎｇｅ Ｚｅｒｎｉｋｅ多項式で展開した時の係数で第５項に１０ｎｍ、第８項に５０ｎｍ、第９項に８０ｎｍとなる収差量を有する。また軸上色収差を有し、波長が１ｎｍ異なると焦点位置が６６２ｎｍだけ光軸方向に変化する。駆動装置１０５はコンピュータ１０６によって指示された光軸方向の位置ｚに移動し、撮像素子１０４によって光強度分布を取得する。

本実施例における収差の推定は、図２および図３に示されるフローチャートに沿って行われる。ステップＳ１において、コンピュータ１０６は、駆動装置１０５と撮像素子１０４をそれぞれ制御することで、光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）を異なるデフォーカス位置で取得する。計測を行った位置ｚ_ｊは、－４９．４μｍ、－３２．９μｍ、－１６．５μｍ、０μｍ、１６．５μｍ、３２．９μｍ、４９．４μｍの計７か所である。取得された光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）はコンピュータ１０６のハードディスクに保存される。ステップＳ２においてコンピュータ１０６は保存された光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）を読み込む。

続いてステップＳ２において、コンピュータ１０６は、ハードディスクに保存された光強度分布Ｉ_ｍ（ｘ，ｙ，ｚ_ｊ）を読み出し、一時メモリに格納する。続いてステップＳ３において、コンピュータ１０６は、被検光学系１０３の収差を演算によって推定する。本実施例では演算をより簡単にするため、ステップＳ３にて以下の演算式を用いる。

ｚ_ｊは等間隔で定められており（ｚ_ｊ＝（ｊ－１）Δｚ+ｚ_０）、更にｓｉｎｃ関数の特性からｚ－ｚ_ｊが大きい時にβ（ｚ－ｚ_ｊ）は無視できるほど小さくなるため、式（１４）は以下の式（１５）とみなすことができる。

すなわち、計算対象となる位置ｚ_ｊの前後のデフォーカス位置（ｚ_ｊ－Δｚ及びｚ_ｊ+Δｚ）で計算された光強度分布に重み（β（Δｚ）およびβ（－Δｚ））を掛けて足し合わせるだけで、色収差による影響が考慮された光強度分布Ｉを算出可能である。この式に基づいて、ステップＳ３を実行する。

まずステップＳ３１において、コンピュータ１０６は、式（１６ａ）～（１６ｃ）に基づいて重み係数βを算出する。本実施例では、ｚ_１＝－１６．５μｍ、ｚ_２＝１６．５μｍ、Δｚ＝１６．５μｍとして積分演算を実行する。続いてステップＳ３２において、コンピュータ１０６は、計測を行った各ｚ_ｊにおける単色での光強度分布Ｉ’（ｘ，ｙ，ｚ_ｊ；Ｗ；λ_０）を中心波長λ_０＝５２５ｎｍで算出した。Ｗの初期値としては無収差を用いる。

続いてステップＳ３３において、コンピュータ１０６は、βおよびＩ’（ｘ，ｙ，ｚ_ｊ；Ｗ；λ_０）から、対象となる位置ｚ_ｊの前後の位置（ｚ_ｊ－Δｚ及びｚ_ｊ+Δｚ）で計算された光強度分布Ｉ’に重み（β（Δｚ）およびβ（－Δｚ））を掛けて足し合わせる。これにより、Ｉ（ｘ，ｙ，ｚ_ｊ；Ｗ）を得ることができる。

続いてステップＳ３４において、コンピュータ１０６は、色収差が考慮された光強度分布Ｉ（ｘ，ｙ，ｚ_ｊ；Ｗ）をＩ_ｓ（ｘ，ｙ，ｚ_ｊ；Ｗ）として、式（１）に基づき目的関数Ｆ（Ｗ）を算出する。続いてステップＳ３５において、コンピュータ１０６は、目的関数Ｆ（Ｗ）を評価する。そしてステップＳ３６において、コンピュータ１０６は、Ｗを勾配法によって更新する。

図４は、本実施例の以上の処理により推定された収差Ｗの説明図であり、真値、従来例、および本実施例を比較して示す。収差Ｗは、Ｆｒｉｎｇｅ Ｚｅｒｎｉｋｅ多項式で展開した係数で示している。図４において、横軸はＺｅｒｎｉｋｅ項、縦軸は係数をそれぞれ示す。本実施例の推定結果である破線は、従来例の点線に比べて真値の実線に近い結果を得られている。すなわち、本実施例によって高精度な推定が可能となる。

次に、実施例２における収差推定方法について説明する。本実施例は、実施例１とは異なり、計測された複数Ｉ_ｍの光強度分布に対して補正処理を加えることで、高速かつ高精度な収差推定を実現する。以下、本実施例の原理を述べ、続いて原理に基づいた収差推定方法を述べる。

前述のとおり、最適化演算を用いた収差推定では、Ｉ_ｓを算出する計算条件がＩ_ｍを計測する計測条件と一致している必要がある。実施例１においては、Ｉ_ｓの算出過程で式（１３）に基づきＩ’をＩに変換することで不一致を解消した。不一致を解消する他の１つの方法が、計測された光強度分布Ｉ_ｍに補正処理を加えることで、光源１０１が単色の場合に取得される光強度分布Ｉ_ｍ’を算出することである。

式（１３）の両辺に行列Ｂの逆行列を左から掛けると、以下の式（１７）が得られる。

波長幅がある光源１０１によって得られた光強度分布Ｉは、線形変換、すなわち重み付き和によって単色での光強度分布Ｉ’に変換される。よって、式（１７）においてＩをＩ_ｍ、Ｉ’をＩ_ｍ’とし、式（１）のＩ_ｍをＩ_ｍ’に置き換えれば、Ｉ_ｓとの不一致が解消され、高精度な波面収差の推定が可能となる。

図５を参照して、この原理に基づいてコンピュータ１０６が実行する収差推定方法を説明する。図５は、本実施例における収差推定方法を示すフローチャートである。なお、ステップＳ３１は実施例１と同じである。本実施例では、ステップＳ３１の次にステップＳ３７が実行される。ステップＳ３７において、コンピュータ（合成部）１０６は、計測された複数の光強度分布Ｉ_ｍから式（１７）に基づいて合成処理を行い、単色での光強度分布Ｉ_ｍ’を算出する。すなわちコンピュータ１０６は、複数の第１光強度分布を合成する。このときに用いられるＢ^－１は、式（１１）および式（１４）から定まるβ及びＢから算出される。

続くステップＳ３２は、実施例１と同じである。ステップＳ３４において、コンピュータ１０６は、ステップＳ３２にて算出されたＩ’を式（１）におけるＩ_ｓとし、ステップＳ３７にて算出されたＩ_ｍ’を式（１）におけるＩ_ｍとして用いることで、目的関数Ｆ（Ｗ）を算出する。なお、ステップＳ３５、Ｓ３６は実施例１と同じである。

式（１７）ではＢの逆行列を用いているが、これに限定されるものではない。Ｂの疑似逆行列を用いてもよく、または、画像処理で用いられるウィーナーフィルタの形式を取った行列を用いることもできる。いずれの行列を用いたとしても、計測された複数の光強度分布のうち少なくとも２つ以上から合成処理を行うことによって計測条件を式（２）の計算条件に一致させることができる。そして、条件が一致することにより、高精度な収差推定を実現することができる。

前述の各実施例では、説明を容易にするため光軸上近傍で議論を行ったが、光軸外に結像する光に対しても機能する。光軸外では主光線に傾きがあるため、計測位置ｚが変わることによって画像が横にシフトする。従って、式（１３）や式（１７）の演算を行う前に、主光線の傾きに合わせた画像のシフトを加える必要がある。

各実施例では、式（１１）または式（１６）の積分演算を行う範囲（ｚ_１およびｚ_２）を被検光学系１０３の色収差に基づいて定める。しかし、被検光学系１０３の色収差は事前に知ることができない場合がある。その場合、予想される色収差より少し広い範囲で積分範囲を定めておけばよい。または、光強度分布を取得した間隔Δｚに基づいて、積分範囲をΔｚの１倍から３倍程度の範囲としておけば、各実施例は機能する。また、前述のとおり、Ｉ’はｚ方向にλ_０／ＮＡ^２程度の細かさの変化しかないため、Δｚ～λ_０／ＮＡ^２としておけばよい。

各実施例では、式（１２）におけるｚ_ｊとｚ_ｎは同じ値を持つとして、ステップＳ３２での演算を行うが、ｚ_ｊとｚ_ｎは異なっていてもよい。２つ以上の異なる位置ｚで算出された光強度分布Ｉ’に基づいてＩが算出されていれば、各実施例は機能する。

各実施例では、ステップＳ３２におけるＩ’の算出において、すべてのｚ_ｊにおけるＩ’対して同じλ_０を用いるが、これに限定されるものではない。また各実施例では、説明を簡略化するため、ｆ_ｎ（ｚ）としてｓｉｎｃ関数を用いて、式（１３）の線形変換や式（１１）の重み係数βの導出方法を説明した。より一般には、式（８）を式（６）に代入して得られる式（１８）および式（１９）となる。

ａ_ｎ（ｘ，ｙ；Ｗ；λ_０）は、Ｉ’（ｘ、ｙ、ｚ_ｊ；Ｗ；λ_０）の線形変換で求まるため、ＩとＩ’は線形変換、すなわち重み付き和で互いに変換される。各実施例は、計算された光強度分布または計測された光強度分布のうち少なくとも一方の光強度分布に対して、少なくとも２つ以上のデフォーカス位置における光強度分布を被検光学系１０３の色収差に基づいた重み付き和によって変換する。これにより、計算条件と計測条件とを一致させることができる。そして条件が一致したことにより、高精度な収差推定が実現される。より広義には、各実施例は、少なくとも２つのデフォーカス位置における光強度分布に対する合成処理を有する収差推定方法に関する。

また各実施例では、収差の推定演算として最適化演算を例示しているが、これに限定されるものではない。光強度分布から収差を推定する演算方法は、強度輸送方程式を解く方法や機械学習を用いた方法などもある。多くの手法が単色での演算を仮定しているため、実施例２で示したＩ_ｍに対する合成処理を加えることで、他の演算方法に対しても効果がある。

（その他の実施例）
本発明は、上述の実施例の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。

各実施例によれば、広い波長幅を有する光源を用いても、高速かつ高精度で被検光学系の収差を推定することが可能な収差推定方法、収差推定装置、演算装置、プログラム、および記憶媒体を提供することができる。

各実施例の開示は、以下の構成および方法を含む。

（方法１）
互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得する計測ステップと、
前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出する算出ステップと、
前記複数の第１光強度分布または前記複数の第２光強度分布を合成する合成ステップとを有し、
前記合成ステップの結果に基づいて被検光学系の収差を推定することを特徴とする収差推定方法。
（方法２）
前記合成ステップは、前記複数の第２光強度分布、または前記複数の第１光強度分布に対して、前記複数の光強度分布の重み付き和を算出することを特徴とする方法１に記載の収差推定方法。
（方法３）
前記合成ステップは、前記被検光学系の色収差に基づいて実行されることを特徴とする方法１または２に記載の収差推定方法。
（方法４）
前記合成ステップは、前記被検光学系の色収差に基づいて算出される重み係数を用いることを特徴とする方法１乃至３のいずれかに記載の収差推定方法。
（方法５）
前記計測ステップにおいて、撮像素子を前記被検光学系の光軸方向に沿って移動させることで、前記複数の第１光強度分布を計測することを特徴とする方法１乃至４のいずれかに記載の収差推定方法。
（構成１）
光源と、
被検光学系により形成される光学像を光電変換して画像データを出力する撮像素子と、
前記画像データに基づいて被検光学系の収差を推定する制御部とを有し、
前記制御部は、
互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得し、
前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出し、
前記複数の第１光強度分布または前記複数の第２光強度分布を合成して、前記被検光学系の収差を推定することを特徴とする収差推定装置。
（構成２）
前記光源は、ＬＥＤ、蛍光灯、または白熱灯であることを特徴とする構成６に記載の収差推定装置。
（構成３）
互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得する取得部と、
前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出する算出部と、
前記複数の第１光強度分布または前記複数の第２光強度分布を合成する合成部と、
前記合成部の結果に基づいて被検光学系の収差を推定する推定部と、を有することを特徴とする演算装置。
（構成４）
方法１乃至５のいずれか一項に記載の収差推定方法をコンピュータに実行させることを特徴とするプログラム。
（構成５）
構成４に記載のプログラムを記憶していることを特徴とするコンピュータ読み取り可能な記憶媒体。

以上、本発明の好ましい実施形態について説明したが、本発明はこれらの実施形態に限定されず、その要旨の範囲内で種々の変形及び変更が可能である。

１００ 収差推定装置
１０１ 光源
１０３ 被検光学系
１０４ 撮像素子
１０６ コンピュータ（制御部）

Claims (10)

1. 互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得する計測ステップと、
   前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出する算出ステップと、
   前記複数の第１光強度分布または前記複数の第２光強度分布を合成する合成ステップとを有し、
   前記合成ステップの結果に基づいて被検光学系の収差を推定することを特徴とする収差推定方法。
2. 前記合成ステップは、前記複数の第２光強度分布、または前記複数の第１光強度分布に対して、前記複数の光強度分布の重み付き和を算出することを特徴とする請求項１に記載の収差推定方法。
3. 前記合成ステップは、前記被検光学系の色収差に基づいて実行されることを特徴とする請求項１に記載の収差推定方法。
4. 前記合成ステップは、前記被検光学系の色収差に基づいて算出される重み係数を用いることを特徴とする請求項１に記載の収差推定方法。
5. 前記計測ステップにおいて、撮像素子を前記被検光学系の光軸方向に沿って移動させることで、前記複数の第１光強度分布を計測することを特徴とする請求項１に記載の収差推定方法。
6. 光源と、
   被検光学系により形成される光学像を光電変換して画像データを出力する撮像素子と、
   前記画像データに基づいて被検光学系の収差を推定する制御部とを有し、
   前記制御部は、
   互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得し、
   前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出し、
   前記複数の第１光強度分布または前記複数の第２光強度分布を合成して、前記被検光学系の収差を推定することを特徴とする収差推定装置。
7. 前記光源は、ＬＥＤ、蛍光灯、または白熱灯であることを特徴とする請求項６に記載の収差推定装置。
8. 互いに異なる複数のデフォーカス位置で被検光学系を介して物体を撮像することにより、該複数のデフォーカス位置に対応する複数の第１光強度分布を取得する取得部と、
   前記複数のデフォーカス位置に対応する複数の第２光強度分布を算出する算出部と、
   前記複数の第１光強度分布または前記複数の第２光強度分布を合成する合成部と、
   前記合成部の結果に基づいて被検光学系の収差を推定する推定部と、を有することを特徴とする演算装置。
9. 請求項１乃至５のいずれか一項に記載の収差推定方法をコンピュータに実行させることを特徴とするプログラム。
10. 請求項９に記載のプログラムを記憶していることを特徴とするコンピュータ読み取り可能な記憶媒体。