(57)【要約】【課題】ブロック処理を行うゆらぎ補正において、補正ずれを抑制すること。【解決手段】時間的に連続して撮像された複数の画像データを用いて大気のゆらぎによる影響を補正する画像処理装置を提供する。画像処理装置は、複数の画像データの処理対象の領域の画像に対応する複数の部分画像データを取得する取得手段と、取得された複数の部分画像データのそれぞれをフーリエスペクトルに変換する変換手段と、変換された複数のフーリエスペクトルの低周波成分の位相を用いて高周波成分を含む位相を推定する推定手段と、推定手段によって推定された結果と、取得手段によって取得された部分画像データに基づいて決定される基準データとに基づいて、推定手段によって推定された少なくとも一つの部分画像データに対応する結果を修正する修正手段と、修正手段によって修正された結果に基づいて画像データに対応する出力画像データを生成する生成手段とを備える。【選択図】図３
2019169000.tif 000002
【特許請求の範囲】
【請求項１】
時間的に連続して撮像された複数の画像データを用いて大気のゆらぎによる影響を補正する画像処理装置であって、
前記複数の画像データの処理対象の領域の画像に対応する複数の部分画像データを取得する取得手段と、
前記取得された複数の部分画像データのそれぞれをフーリエスペクトルに変換する変換手段と、
前記変換された複数のフーリエスペクトルの低周波成分の位相を用いて高周波成分を含む位相を推定する推定手段と、
前記推定手段によって推定された結果と、前記取得手段によって取得された部分画像データに基づいて決定される基準データとに基づいて、前記推定手段によって推定された少なくとも一つの部分画像データに対応する結果を修正する修正手段と、
前記修正手段によって修正された結果に基づいて前記画像データに対応する出力画像データを生成する生成手段と
を備えることを特徴とする画像処理装置。
【請求項２】
前記修正手段は、前記推定手段によって推定された結果と、前記取得手段によって取得された部分画像データに基づいて決定される基準データとの変化量が所定の範囲に含まれる場合、前記推定手段によって推定された結果を修正することを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記修正手段は、前記推定手段によって推定された位相を修正し、
前記生成手段は、
前記位相が修正されたフーリエスペクトルを、前記処理対象の領域に対応する部分画像データに逆変換し、
各領域に対応する部分画像データを合成することで前記出力画像データを生成することを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記修正手段は、
前記変換手段によって変換された複数のフーリエスペクトルの位相から前記基準データを決定し、
前記推定された位相と前記基準データとの差を平面で近似し、近似した平面の傾きを前記変化量として算出し、前記傾きが前記所定の範囲に含まれるかを判定することを特徴とする請求項３に記載の画像処理装置。
【請求項５】
前記修正手段は、前記推定された位相から前記近似した平面を引くことで前記修正を行うことを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記推定手段は、前記フーリエスペクトルの振幅をさらに推定し、
前記生成手段は、前記推定された振幅を含み、かつ前記位相が修正されたフーリエスペクトルを逆変換することを特徴とする請求項３から５のいずれか一項に記載の画像処理装置。
【請求項７】
前記推定手段は、前記フーリエスペクトルの振幅をさらに推定し、
前記修正手段は、
前記推定された位相および振幅を含むフーリエスペクトルを前記処理対象の領域に対応する部分画像データに逆変換し、
前記逆変換した部分画像データを修正することを特徴とする請求項２に記載の画像処理装置。
【請求項８】
前記修正手段は、
前記複数の部分画像データから前記基準データを決定し、
前記逆変換した部分画像データと前記基準データとから算出される変化量が前記所定の範囲に含まれるかを判定することを特徴とする請求項７に記載の画像処理装置。
【請求項９】
前記生成手段は、前記修正された部分画像データを含む各領域に対応する部分画像データを合成することで前記出力画像データを生成することを特徴とする請求項７または８に記載の画像処理装置。
【請求項１０】
前記推定手段によって推定された結果と前記基準データとの変化量は、画像の平行移動ずれに相当することを特徴とする請求項２から９のいずれか一項に記載の画像処理装置。
【請求項１１】
前記推定手段によって推定された結果と前記基準データとの変化量は、画像の回転ずれに相当することを特徴とする請求項２から９のいずれか一項に記載の画像処理装置。
【請求項１２】
前記推定手段によって推定された結果と前記基準データとの変化量は、画像の拡大ずれまたは縮小ずれに相当することを特徴とする請求項２から９のいずれか一項に記載の画像処理装置。
【請求項１３】
時間的に連続して撮像された複数の画像データを用いて大気のゆらぎによる影響を補正する画像処理方法であって、
前記複数の画像データの処理対象の領域の画像に対応する複数の部分画像データを取得する取得ステップと、
前記取得された複数の部分画像データのそれぞれをフーリエスペクトルに変換する変換ステップと、
前記変換された複数のフーリエスペクトルの低周波成分の位相を用いて高周波成分を含む位相を推定する推定ステップと、
前記推定ステップにおいて推定された結果と、前記取得ステップにおいて取得された部分画像データに基づいて決定される基準データとに基づいて、前記推定ステップにおいて推定された少なくとも一つの部分画像データに対応する結果を修正する修正ステップと、
前記修正ステップにおいて修正された結果に基づいて前記画像データに対応する出力画像データを生成する生成ステップと
を備えることを特徴とする画像処理方法。
【請求項１４】
コンピュータを、請求項１から１２のいずれか一項に記載の画像処理装置の各手段として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、大気のゆらぎに起因する画像への影響を補正する技術に関する。
【背景技術】
【０００２】
従来、監視カメラシステムなどにおいて、大気のゆらぎに起因する画像の劣化を補正する技術（以下、「ゆらぎ補正」と呼ぶ）が用いられている。ここで大気ゆらぎとは、大気の乱流や温度差により、屈折率が空間的および時間的に変化する自然現象を指す。大気ゆらぎは陽炎とも呼ばれる。大気ゆらぎは光の屈折や位相のばらつきを引き起こし、その結果として、特に望遠レンズで遠方を撮像した場合に、被写体の像に時間的に変化するぼけ及び歪みが発生する。このような画像の劣化は、被写体の視認性を大きく低下させるため、除去されることが望ましい。
【０００３】
複数フレームの画像から、ぼけ及び歪みをまとめて低減した画像を生成（再構成）する技術として、スペックルイメージングが知られている。以降の説明では、動画像を構成する時系列の画像データをフレームと呼ぶ。スペックルイメージングは、元来天文分野で開発された手法であり、短時間露光された画像からオリジナルの画像を再構成する技術である。スペックルイメージングにおいては、再構成される画像のフーリエスペクトルの振幅成分と位相成分とを独立に推定する処理が行われる。そして、推定された振幅成分と位相成分とを用いて逆フーリエ変換をすることで、ぼけ及び歪みを低減した画像が再構成される。特許文献１では、スペックルイメージングをゆらぎ補正に応用する技術が開示されている。
【０００４】
また、スペックルイメージングを行う際には、画像の領域ごとの部分画像に対してスペックルイメージングを行い、領域ごとに得られた部分画像を合成して最終出力画像を生成する処理（以下、「ブロック処理」と呼ぶ）が行われる。特許文献１においてもブロック処理を行う例が示されている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】米国特許出願公開第２００４／０００５０９８号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、特許文献１の技術では、スペックルイメージングによって生成される一部の部分画像に位置ずれなどの補正ずれが生じることがある。すると、それらを合成して得られる最終出力画像において特有のアーティファクトが生じてしまい、視認性が低下する。
【０００７】
本発明は、ブロック処理を行うゆらぎ補正において、補正ずれを抑制することを目的とする。
【課題を解決するための手段】
【０００８】
本発明の一態様に係る画像処理装置は、時間的に連続して撮像された複数の画像データを用いて大気のゆらぎによる影響を補正する画像処理装置であって、前記複数の画像データの処理対象の領域の画像に対応する複数の部分画像データを取得する取得手段と、前記取得された複数の部分画像データのそれぞれをフーリエスペクトルに変換する変換手段と、前記変換された複数のフーリエスペクトルの低周波成分の位相を用いて高周波成分を含む位相を推定する推定手段と、前記推定手段によって推定された結果と、前記取得手段によって取得された部分画像データに基づいて決定される基準データとに基づいて、前記推定手段によって推定された少なくとも一つの部分画像データに対応する結果を修正する修正手段と、前記修正手段によって修正された結果に基づいて前記画像データに対応する出力画像データを生成する生成手段とを備えることを特徴とする。
【発明の効果】
【０００９】
本発明によれば、ブロック処理を行うゆらぎ補正において、補正ずれを抑制することができる。
【図面の簡単な説明】
【００１０】
【図１】画像処理装置の構成図。
【図２】画像処理装置の機能ブロック図。
【図３】ゆらぎ補正部のブロック図。
【図４】ゆらぎ補正を説明する図。
【図５】ゆらぎ補正と位置ずれを説明する図。
【図６】位置ずれ修正を説明する図。
【図７】ゆらぎ補正処理のフローチャート。
【図８】ゆらぎ補正処理のフローチャート。
【図９】位置補償処理のフローチャート。
【図１０】位相差分の例を説明する図。
【図１１】位置ずれ修正を説明する図。
【図１２】ゆらぎ補正部のブロック図。
【図１３】ゆらぎ補正処理のフローチャート。
【図１４】位置補償処理のフローチャート。
【発明を実施するための形態】
【００１１】
以下、本発明の実施形態について、図面を参照して説明する。なお、以下の実施形態は本発明を限定するものではなく、また、本実施形態で説明されている特徴の組み合わせの全てが本発明の解決手段に必須のものとは限らない。なお、同一の構成については、同じ符号を付して説明する。
【００１２】
＜＜実施形態１＞＞
本実施形態では、スペックルイメージングの過程で導出される位相データを用いて、ゆらぎ補正によって補正ずれが生じているかを判定する。本実施形態では、補正ずれの一例として画像位置の平行移動ずれ（位置ずれ）を扱う。そして、判定結果に基づいて、ゆらぎ補正によって生じる位置ずれを修正する例を述べる。
【００１３】
＜画像処理装置のハードウェア構成＞
図１は、本実施形態の画像処理装置のハードウェア構成の一例を示す図である。本実施形態の画像処理装置１００は、ＣＰＵ１０１、ＲＡＭ１０２、ＲＯＭ１０３、二次記憶装置１０４、入力インターフェース１０５、および出力インターフェース１０６を含む。画像処理装置１００の各構成要素は、システムバス１０７によって相互に接続されている。画像処理装置１００は、入力インターフェース１０５を介して外部記憶装置１０８および操作部１１０に接続されている。画像処理装置１００は、出力インターフェース１０６を介して外部記憶装置１０８および表示装置１０９に接続されている。
【００１４】
ＣＰＵ１０１は、ＲＡＭ１０２をワークメモリとして、ＲＯＭ１０３に格納されたプログラムを実行し、システムバス１０７を介して画像処理装置１００の各構成要素を統括的に制御する。これにより、後述する様々な処理が実行される。二次記憶装置１０４は、画像処理装置１００で取り扱われる種々のデータを記憶する記憶装置であり、本実施形態ではＨＤＤが用いられる。ＣＰＵ１０１は、システムバス１０７を介して二次記憶装置１０４へのデータの書き込みおよび二次記憶装置１０４に記憶されたデータの読出しを行う。なお、二次記憶装置１０４にはＨＤＤの他に、光ディスクドライブやフラッシュメモリなど、様々な記憶デバイスを用いることが可能である。
【００１５】
入力インターフェース１０５は、例えばＵＳＢやＩＥＥＥ１３９４等のシリアルバスインターフェースである。画像処理装置１００は、入力インターフェース１０５を介して、外部装置からデータや命令等を入力する。本実施形態では、画像処理装置１００は、入力インターフェース１０５を介して、外部記憶装置１０８（例えば、ハードディスク、メモリカード、ＣＦカード、ＳＤカード、ＵＳＢメモリなどの記憶媒体）からデータを取得する。また本実施形態では、画像処理装置１００は、操作部１１０に入力されたユーザの指示を、入力インターフェース１０５を介して取得する。操作部１１０は、マウスやキーボードなどの入力装置であり、ユーザの指示を入力する。
【００１６】
出力インターフェース１０６は、入力インターフェース１０５と同様にＵＳＢやＩＥＥＥ１３９４等のシリアルバスインターフェースである。なお、出力インターフェース１０６は、例えばＤＶＩやＨＤＭＩ（登録商標）等の映像出力端子であってもよい。画像処理装置１００は、出力インターフェース１０６を介して、外部装置にデータ等を出力する。本実施形態では、画像処理装置１００は、出力インターフェース１０６を介して表示装置１０９（液晶ディスプレイなどの各種画像表示デバイス）に、ＣＰＵ１０１によって処理されたデータ（例えば、画像データ）を出力する。なお、画像処理装置１００の構成要素は上記以外にも存在し得るが、ここでは説明を省略する。
【００１７】
＜画像処理装置の機能ブロック図＞
図２は、画像処理装置１００の機能ブロックの例を示す図である。画像処理装置１００は、画像データ取得部２０１、画像分割部２０２、ゆらぎ補正部２０３、および画像結合部２０４を有する。
【００１８】
図３は、ゆらぎ補正部２０３の詳細を示すブロック図である。ゆらぎ補正部２０３は、フーリエ変換部３０１、振幅推定部３０２、位相推定部３０３、位置補償部３０４、および逆フーリエ変換部３０５を有する。位置補償部３０４は、基準位相決定部３５１、位置ずれ判定部３５２、および位置ずれ修正部３５３を有する。各部で行われる処理の説明については、後述するフローチャートの説明とともに行うこととする。
【００１９】
図２および図３に示す各部の機能は、ＣＰＵ１０１がＲＯＭ１０３内に格納された制御プログラムを読み込み実行することで実現される。なお、各部に相当する専用の処理回路を備えるように画像処理装置１００を構成するようにしてもよい。
【００２０】
＜一般的なゆらぎ補正の説明＞
図４は、一般的なゆらぎ補正の一例を説明する図である。先に説明したように、複数の画像データ（フレーム画像）を用いて、大気のゆらぎによる、画像のぼけ及び歪み等をまとめて低減した画像を生成することが行われている。図４では、複数のフレーム４０１（入力画像データ）が画像処理装置１００に入力された場合に、１枚のゆらぎ補正済みの出力画像データ４４１が得られる例を示している。ゆらぎ補正を行う手法として、「スペックルイメージング処理」が知られている。
【００２１】
ゆらぎ補正処理においては、フレーム４０１の処理対象の領域に対応する部分画像をフレームごとに切り出して、複数の部分画像データ４０２が取得される。そして、複数の部分画像データ４０２をそれぞれフーリエ変換して、周波数空間の成分（位相成分４１１と振幅成分４２１）にそれぞれ変換する。そして、周波数空間において、再構成される画像の位相成分４１２と振幅成分４２２とをそれぞれ独立して推定する処理が行われる。振幅成分４２２は、例えば入力された各部分画像のパワースペクトルの平均を導出し、導出された平均に対し撮像装置の光学系や大気の影響による劣化モデルに基づいて決定された高周波成分を増幅するフィルタをかけることで推定される。位相成分４１２は、例えば入力された各部分画像のバイスペクトラムに基づき、既知の低周波成分の位相を用いて高周波成分の位相を逐次的に推定していくことで求められる。位相成分４１２の推定処理の詳細は後述する。
【００２２】
このように推定された位相成分４１２と振幅成分４２２とを逆フーリエ変換（逆変換）することで、ゆらぎ補正済みの部分画像データ４３１が得られる。このような処理を、フレームの各領域を処理対象の領域として処理することで、各領域に対応するゆらぎ補正済み部分画像データが得られることになる。そして、ゆらぎ補正済み部分画像データ４３１を合成することで、ゆらぎ補正がされた出力画像データ４４１が得られることになる。
【００２３】
このように、部分画像データからぼけ及び歪みを低減した後に、各部分画像データを合成して、出力画像データを生成する処理（ブロック処理）が行われる。このとき一部の部分画像において、スペックルイメージングによって生成される画像に位置ずれなどの補正ずれが生じることがある。すると、それらを合成して得られる出力画像データにおいて特有のアーティファクトが生じてしまい、視認性が低下してしまうことがある。
【００２４】
図５は、ゆらぎ補正と位置ずれを説明する図である。図５（ａ）は、理想的な画像５１０を示している。画像５１０は、周波数成分５１５−１、５１５−２、５１５−３・・・から構成されている。しかしながら、実際には、大気のゆらぎの影響によって、撮像された画像は、図５（ｂ）の画像５２０となってしまっている。大気のゆらぎの影響を受けた画像５２０は、周波数毎に異なる方向で異なる量の位置ずれが生じている。このため、異なる方向で異なる量の位置ずれが生じている周波数成分５２５−１、５２５−２、５２５−３・・・によって構成される画像５２０の境界がボケて見える。図４に示すゆらぎ補正処理（スペックルイメージング処理）は、大気のゆらぎの影響を補正して理想的な画像５１０を得るための処理である。
【００２５】
しかしながら、図４に示すようなゆらぎ補正処理（スペックルイメージング処理）を行った場合、図５（ｃ）の周波数成分５３５−１、５３５−２、５３５−３・・・のように、全ての周波数に同じ方向で同じ量の位置ずれが残存することがある。これらの周波数成分５３５−１、５３５−２、５３５−３・・・によって構成される画像５３０のボケは解消されているが、理想的な画像５１０に比べて位置ずれが生じてしまっている。
【００２６】
本実施形態では、このような位置ずれを修正する処理を行うものである。より詳細には、本実施形態では、図４に示すようなスペックルイメージングの過程で導出される位相成分４１１を用いて位置ずれ判定を行い、判定結果に基づいて位置ずれ修正を行う形態を説明する。
【００２７】
＜本実施形態のゆらぎ補正の説明＞
図６は、本実施形態で説明するゆらぎ補正における位置ずれを修正する処理を説明する図である。部分画像データ４０２をフーリエ変換する処理は、図４の例と同様である。本実施形態では、スペックルイメージングの過程で導出される位相成分４１１を用いて位置ずれ判定を行い、判定結果に基づいて位置ずれ修正を行う形態である。よって、図６において周波数空間の成分は、位相成分のみを図示している。振幅成分については、図４で説明した例と同様であり、図示を省略している。
【００２８】
本実施形態では、ゆらぎ補正部２０３は、部分画像データの位相成分４１１を用いて、位置ずれの判定に用いられる基準データ（基準位相成分６１０）を決定する。ゆらぎ補正部２０３は、推定された位相成分４１２と、基準位相成分６１０とを比較して位置ずれ量を取得する。そして、位置ずれ量が所定の閾値より大きい場合、位置ずれが生じていると判定し、ゆらぎ補正部２０３は、位置ずれを修正した位相成分６１５を導出する。ゆらぎ補正部２０３は、この位置ずれを修正した位相成分６１５と、推定された振幅成分４２２（図６では不図示）とを逆フーリエ変換して、位置ずれを修正した、ゆらぎ補正済みの部分画像データ６２０を生成する。そして、各注目領域に対応する部分画像データ６２０を合成して、出力画像データ６３０が生成される。以下、本実施形態の処理を詳細に説明する。
【００２９】
＜全体のフローチャート＞
図７は、本実施形態における処理の流れを示すフローチャートである。図７に示す処理は、ＣＰＵ１０１が、ＲＯＭ１０３内に格納された制御プログラムを読み込んでＲＡＭ１０２に展開し、これをＣＰＵ１０１が実行することで実現される。あるいはまた、図７におけるステップの一部または全部の機能をＡＳＩＣや電子回路等のハードウェアで実現してもよい。なお、各処理の説明における記号「Ｓ」は、当該フローチャートにおけるステップであることを意味する（本明細書において同様である）。以下、図２と図７とを用いて説明する。
【００３０】
Ｓ７０１において画像データ取得部２０１は、ゆらぎ補正を行う処理対象である画像データを取得する。本実施形態では、動画像を構成する時間的に連続した複数のフレームを取得する。時間的に連続した複数の静止画像のデータを取得する形態でもよい。
【００３１】
Ｓ７０２において画像分割部２０２は、Ｓ７０１で取得した画像データを複数の小領域データ（部分画像データ）に分割する。即ち、各画像データを、複数の部分画像データに分割する。ここで複数の部分画像データは、領域の重複が許容される。部分画像データのサイズは、大気ゆらぎの影響が一定であるとみなせる範囲（アイソプラナティックパッチ）よりも小さいことが好ましい。Ｓ７０２の処理により、Ｓ７０１で取得された複数の画像データのそれぞれが、各小領域に対応する部分画像データに分割される。以降の処理では、小領域ごとに処理が行われる。
【００３２】
Ｓ７０３においてゆらぎ補正部２０３は、全ての小領域に対して処理をしたかを判定する。全ての小領域に対して処理をした場合、Ｓ７０６に進む。未処理の小領域がある場合、Ｓ７０４に進み、ゆらぎ補正部２０３は、未処理の小領域の中から処理対象とする小領域を決定する。
【００３３】
Ｓ７０５においてゆらぎ補正部２０３は、処理対象の小領域の部分画像データを取得し、ゆらぎ補正処理を行う。即ち、ゆらぎ補正部２０３は、処理対象の小領域に対応する、複数フレーム分の複数の部分画像データを用いてゆらぎ補正処理を行う。詳細は後述する。
【００３４】
Ｓ７０４において画像結合部２０４は、Ｓ７０５でゆらぎ補正された各小領域に対応する推定画像（部分画像）を合成して、出力画像を生成する。
【００３５】
＜ゆらぎ補正処理のフローチャート＞
次に、図３と図８とを用いてＳ７０５のゆらぎ補正処理の詳細を説明する。Ｓ８０１においてフーリエ変換部３０１は、各フレームの処理対象の小領域に対応する部分画像データのそれぞれに対してフーリエ変換を行ってフーリエスペクトル（周波数成分）を生成する。Ｓ８０１により、フレーム数分のフーリエスペクトルが得られる。
【００３６】
Ｓ８０２において振幅推定部３０２は、Ｓ８０１で得られた各フレーム（各部分画像）のフーリエスペクトルを用いて、推定画像のフーリエスペクトルＯの振幅を推定する。具体的には、例えば特許文献１に示されているように、式（１）に従いフーリエスペクトルＯの振幅が導出される。
【００３７】
【数１】
2019169000.tif 000003
ここで、Ａは、Ｓ８０１で得られた各フレームのフーリエスペクトルから導出されるパワースペクトルのフレーム平均である。Ｔは、Ｓ７０１で取得された画像データの撮像に用いられた撮像装置の光学伝達関数である。Ｓは、スペックル伝達関数である。Ｔは、予め大気ゆらぎなどの外乱がない状態で点光源を撮像して求めてもよいし、モデルを用いて数値計算によって求めてもよい。Ｓは、大気ゆらぎに起因する画像の劣化を表現する関数であり、予め計測してデータを保持しても良いし、モデルを用いて数値計算によって求めても良い。また、大気ゆらぎの強度に応じて異なる複数の計測データやモデルを保持しておき、選択して用いるようにしても良い。 
Ｓ８０３において位相推定部３０３は、Ｓ８０１で得られた各フレームのフーリエスペクトルを用いて、推定画像のフーリエスペクトルＯの位相（推定位相φ）を推定する。特許文献１に示されている通り、この推定位相φは式（２）によって決定される。
【００３８】
【数２】
2019169000.tif 000004
ここで、＜＞は複数フレーム（ｎ個のフレーム）にわたる平均を表す。またａｒｇ［＊］は複素数＊の偏角を表す。ＩB,n（ｕ，ｖ）は、ｎフレーム目の部分画像データのフーリエスペクトルＩnを用いて式（３）によって定義される（ｕ，ｖ）に対するバイスペクトラムである。
【００３９】
【数３】
2019169000.tif 000005
式（２）によると、周波数平面上の（０，０）、（１，０）、（０，１）の３点における位相を初期位相として予め決めておけば、逐次的に任意の座標における位相を推定することができる。以下、位相を推定する処理を説明する。
【００４０】
式（２）によれば、ｕとｖの位相を用いてｕ＋ｖの位相を推定できる。例えば、ｕ＝（０，１）、ｖ＝（１，０）とすると、これらを用いてｕ＋ｖ＝（１，１）の位相を推定することができる。一方、ｕ＝（０，１）の位相は−ｕ＝（０，−１）のマイナス１倍である。よって、ｕの位相が既知ということであれば−ｕの位相も既知となる。以上の点を考慮すると、（０，０）、（０，１）、（１，０）、（０，−１）、（−１，０）、の５個を初期位相として推定済みの位相とすることができる。
【００４１】
以降のステップでは、直前のステップまでで推定済みの位相から任意に２個を選択して推定を行い、推定済みの位相をさらに用いて以降のステップで用いる。例えば、以下の推定が可能である。
【００４２】
（０，１）、（０，１）から（０，２）の位相を推定
（０，１）、（１，０）から（１，１）の位相を推定
（０，１）、（０，−１）から（０，０）の位相を推定
（０，１）、（−１，０）から（−１，１）の位相を推定
（１，０）、（１，０）から（２，０）の位相を推定
（１，０）、（０，−１）から（１，−１）の位相を推定
（１，０）、（−１，０）から（０，０）の位相を推定
（０，−１）、（０，−１）から（０，−２）の位相を推定
（０，−１）、（−１，０）から（−１，−１）の位相を推定
（−１，０）、（−１，０）から（−２，０）の位相を推定
このような処理を繰り返すことで、全ての周波数座標における位相が推定される。なお、座標ベクトルの絶対値が周波数の高低を示している。即ち、上述の例では、絶対値が小さい（低周波）座標から絶対値の大きい（高周波）座標の位相を推定していることになる。
【００４３】
上述した初期位相として、例えば複数フレームを平均した部分画像データのフーリエスペクトルの位相の値を用いることができる。なお、初期位相は上記と異なる３点に対して定めても良いし、周波数平面上の４点以上に対して定めても良い。
【００４４】
Ｓ８０４において位置補償部３０４は、Ｓ８０３で得られた推定位相φに対して位置ずれを判定し、判定結果に基づいて、位置ずれを修正する処理を行う。位置ずれを判定し、判定結果に基づいて、位置ずれを修正する処理のことを位置補償処理ともいう。位置補償処理の詳細は後述する。
【００４５】
Ｓ８０５において逆フーリエ変換部３０５は、推定画像のフーリエスペクトルＯに対して逆フーリエ変換を行って推定画像を生成する。即ち、逆フーリエ変換部３０５は、振幅推定され、かつ位相推定（必要に応じて位置ずれ修正）がされたフーリエスペクトルＯに対して逆フーリエ変換を行う。以上の処理により、処理対象の小領域に対応する部分画像データのゆらぎ補正が行われる。なお、図８の例では、振幅推定の後に位相推定を行うフローチャートを示しているが、この限りではない。位相推定の後に振幅推定を行っても良いし、振幅推定と位相推定とを並行して行っても良い。
【００４６】
図９は、Ｓ８０４の位置補償処理の詳細を示すフローチャートである。Ｓ９０１において基準位相決定部３５１は、Ｓ８０１で得られた各フレームのフーリエスペクトルを用いて、基準位相を決定する。基準位相は、例えば複数フレームの平均を用いることができる。この平均は、ゆらぎ補正で用いるフレーム全てについてでも良く、一部でも良い。また、任意の１フレームの位相でも良い。あるいは、平均値でなく、中央値など別の統計量で基準位相を求めても良い。また、別途の位置推定処理を行って基準画像を生成しても良い。本実施形態では、ｕ＝（ｕ1、ｕ2）とした場合、基準位相として、式（４）で導出されるｎフレーム分の平均のφaveを用いる。
【００４７】
【数４】
2019169000.tif 000006
Ｓ９０２において位置ずれ判定部３５２は、Ｓ８０３で得られた推定位相φとＳ９０１で導出された基準位相とに基づいて位置ずれ判定を行う。位置ずれ判定の一例を示す。まず、推定位相φと基準位相φaveとの変化量を導出する。例えば、推定位相φと基準位相φaveの差分（を２πで割った余り）である位相差分を算出する。図１０は位相差分の一例である。図１０（ａ）は、位置ずれが無い場合の例であり、図１０（ｂ）は、位置ずれが小さい場合の例であり、図１０（ｃ）は、位置ずれが大きい場合の例である。なお、位相差分は実際には２次元のデータであるが、図１０では、その内１つの次元に対応する要素がゼロのもののみ抜き出した１次元のデータを図示している。横軸は周波数を表し右側ほど高周波である。縦軸は位相を表し、全ての位相を０から２πの範囲でプロットされる。黒丸は、上記位相差分をプロットしたものである。図１０（ａ）から（ｃ）は、いずれも、原点を通る直線（２πになったらゼロに戻る）であり、傾きがずれ量に比例する。傾きゼロが、ずれ量ゼロに相当し、傾きが大きいほどずれ量が大きくなる。一般に、２つの画像が平行移動で重なる関係にある場合、そのフーリエスペクトルの位相差分は、適当なＡを用いて式（５）で近似（図中２重線）できることが知られている。
【００４８】
【数５】
2019169000.tif 000007
このときのＡの絶対値が大きいほど平行移動量が大きいことを意味する。位置ずれ判定部３５２は、Ａの絶対値が予め定める閾値より大きい場合に位置ずれが生じていると判定する。そうでない場合、位置ずれが生じていないと判定する。そして位置ずれ修正部３５３は、位置ずれが生じている場合、位置ずれを修正すべきであると決定し、そうでない場合に修正すべきでないと決定する。
【００４９】
推定位相φと基準位相φaveとの変化量が所定の範囲に含まれる場合、例えば、位相差分が予め定める閾値を超える場合に、位置ずれを修正する理由を説明する。ゆらぎ補正は、元来、局所的な位置ずれを補正するものである。つまり、位置ずれが小さい場合とは、本来的なゆらぎ補正をした結果のものと考えられる。従って、位相差分が所定の閾値を超えていない場合には、位置ずれを修正しない。また、ごく小さい位置ずれであれば、見た目がほとんど変わらないので、処理時間短縮のために、位置ずれを修正しない場合もある。また、１画素以下の位置ずれを修正しようとすると、補間処理によるボケの様な弊害が生じる場合もある。よって、本実施形態では、所定の閾値を超える場合に位置ずれを修正する処理を行っている。
【００５０】
ここまで簡単のために１次元のデータで説明したが、実際には２次元のデータに対して本処理は行われる。すなわち、位置ずれ判定部３５２は、位相差分を平面（を２πで割った余り）で近似し、近似した平面の傾きを変化量として算出し、傾きが所定の範囲に含まれるかを判定する。例えば、位置ずれ判定部３５２は、位相差分を平面（を２πで割った余り）で近似し、その平面を特徴づけるパラメータを閾値と比較して位置ずれが生じているか否かを判定する。具体例として、位相差分を適当なパラメータＡとＢを用いて式（６）で近似し、このときのＡとＢの絶対値の内で最小のものと閾値とを比較して判定を行うなどの方法がある。
【００５１】
【数６】
2019169000.tif 000008
Ｓ９０３において位置ずれ修正部３５３は、Ｓ９０２の判定の結果、変化量（位相差分）が所定の範囲に含まれる場合には、Ｓ９０４に進み、位置ずれ修正処理を行う。Ｓ９０４において位置ずれ修正部３５３は、Ｓ９０２での判定結果に基づいて、Ｓ８０３で得られた推定位相φを修正する。この処理は、推定位相φから位相差分φappを引くことで実現できる。位相差分が閾値よりも大きくない場合には、Ｓ９０４の処理は、スキップされる。
【００５２】
以上説明したように、本実施形態では、スペックルイメージングにおいて推定された位相の位置ずれが生じている場合に、周波数空間において位相成分を修正することで、位置ずれを修正する処理を行う。このような処理によれば、大気のゆらぎを低減した画像を生成し、かつ、位置ずれが生じている場合には位置ずれを修正した画像を生成することができる。
【００５３】
＜＜実施形態２＞＞
実施形態１では、周波数空間において、推定された位相成分と基準位相との差分に基づいて位置ずれが生じているかを判定し、位置ずれが生じている場合に、周波数空間において位置ずれを修正する形態を説明した。
【００５４】
本実施形態では、スペックルイメージングの結果導出される画像データ（即ち、逆フーリエ変換した画像データ）を用いて位置ずれ判定を行う。そして、判定結果に基づいて実空間において位置ずれ修正を行う形態を説明する。
【００５５】
図１１は、本実施形態の位置ずれ修正の概要を説明する図である。複数の部分画像データ４０２からゆらぎ補正済み部分画像データ４３１が得られる処理については図４に示す例と同様である。本実施形態では、部分画像データ４０２から基準データとなる画像データ（基準画像データ１１１０という）が決定される。そして、ゆらぎ補正済み部分画像データ４３１と基準画像データ１１１０とを比較して位置ずれ量が取得される。そして、位置ずれ量が所定の閾値より大きい場合、位置ずれが修正される。このように位置ずれが修正されたゆらぎ補正済みの部分画像データ１１２０を用いて出力画像データ１１３０が生成される。
【００５６】
＜ゆらぎ補正部のブロック図＞
図１２は、本実施形態のゆらぎ補正部２０３のブロック図である。本実施形態のゆらぎ補正部２０３は、実施形態１のゆらぎ補正部と同様の機能を有するが、各機能部の接続が異なる。具体的には、位置補償部１２１０の位置が異なる。また、位置補償部１２１０が有する機能も異なる。位置補償部１２１０は、基準画像決定部１２１１、位置ずれ判定部１２１２、位置ずれ修正部１２１３としての機能を有する。
【００５７】
図１３は、本実施形態における処理の流れを示すフローチャートである。Ｓ１３０１、Ｓ１３０２、Ｓ１３０３は、Ｓ８０１、Ｓ８０２、Ｓ８０３と同等の処理である。Ｓ１３０４において逆フーリエ変換部３０５は、推定画像のフーリエスペクトルに対して逆フーリエ変換を行って推定画像を生成する。Ｓ１３０５において位置補償部１２１０は、推定画像（部分画像）に対して位置補償を行う。
【００５８】
図１４は、Ｓ１３０５の位置補償処理の詳細を説明する図である。Ｓ１４０１において基準画像決定部１２１１は、Ｓ７０２で得られた、処理対象の小領域に対応する部分画像データを用いて、基準画像を決定する。本実施形態では、基準画像として、複数フレーム分の部分画像データのフレーム平均を用いる。なお、基準画像は、例えばゆらぎ補正で用いるフレーム全てについてでも良く、一部でも良い。また、任意の１フレームの画像でも良い。あるいは、平均値でなく、中央値など別の統計量で基準画像を求めても良い。
【００５９】
Ｓ１４０２において位置ずれ判定部１２１２は、Ｓ１３０４で得られた推定画像とＳ１４０１で得られた基準画像とを用いて位置ずれ判定を行う。位置ずれ判定の一例を示す。まず、Ｓ１３０４で得られた推定画像とＳ１４０１で得られた基準画像との相互相関関数を最大にする平行移動量ｄを導出する。そして位置ずれ判定部１２１２は、ｄが予め定める所定の範囲に含まれる場合、例えば予め定める閾値より大きい場合、位置ずれが生じていると判定する。そうでない場合、位置ずれが生じていないと判定する。
【００６０】
Ｓ１４０３において位置ずれ修正部３５３は、Ｓ１４０２の判定の結果、平行移動量が閾値よりも大きい場合には、Ｓ１４０４に進み、位置ずれ修正処理を行う。この処理は、推定画像を平行移動量ｄだけ平行移動することで実現できる。平行移動量が閾値よりも大きくない場合には、Ｓ１４０４の処理はスキップする。
【００６１】
以上説明したように、本実施形態によれば、大気のゆらぎを低減した画像を生成し、かつ、位置ずれが生じている場合には位置ずれを修正した画像を生成することができる。
【００６２】
＜＜その他の実施形態＞＞
上述した実施形態においては、位置ずれの例として平行移動ずれの場合を例に挙げて説明したが、この限りではない。他の位置ずれを修正することもできる。例として、回転（または拡縮）ずれを補正してもよい。この場合、ステップＳ１４０２で推定画像と基準画像の相互相関関数を最大にする回転量ｒ（または拡縮量ｅ）を求め、ステップＳ１４０４で回転量ｒ（または拡縮量ｅ）だけ回転（または拡大もしくは縮小）させれば良い。
【００６３】
本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサーがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。
【符号の説明】
【００６４】
３０１  フーリエ変換部
３０２  振幅推定部
３０３  位相推定部
３０４  位置補償部
３０５  逆フーリエ変換部
【図１】
2019169000.tif 000009
【図２】
2019169000.tif 000010
【図３】
2019169000.tif 000011
【図４】
2019169000.tif 000012
【図５】
2019169000.tif 000013
【図６】
2019169000.tif 000014
【図７】
2019169000.tif 000015
【図８】
2019169000.tif 000016
【図９】
2019169000.tif 000017
【図１０】
2019169000.tif 000018
【図１１】
2019169000.tif 000019
【図１２】
2019169000.tif 000020
【図１３】
2019169000.tif 000021
【図１４】
2019169000.tif 000022