(57)【要約】      （修正有）【課題】手振れやピントのずれなどで劣化した撮影画像を修復する方法の提供。【解決手段】手振れやぼやけによって劣化を受けた劣化画像は、有限領域の離散データであるので、従来から提案されている画像修復アルゴリズムを適用するだけでは、所望の修復画像を得ることはできない。そこで、定義された劣化画像の領域よりも外側に、数学的画像モデルに適合するようなデータを埋め込み、さらにその部分にフィルタリングを行った後に画像修復アルゴリズムを適用することによって、所望の修復画像を得ることができる。【選択図】図１
2007102477.tif 000002
【特許請求の範囲】
【請求項１】
数学的モデルに基づいて、劣化した原画像を修復する画像修復装置であって、
撮像領域の境界部を補正する境界値補正部と、
上記境界値補正部によって境界部を補正された画像に基づいて劣化画像補正パラメータを推定するパラメータ推定部と、
上記境界値補正部によって境界部を補正された画像または上記パラメータ推定部によって推定された劣化画像補正パラメータに基づいて、原画像を推定し、修復画像を出力する原画像推定部と
を有し、
上記境界値補正部は、上記撮像領域の境界の外側をパディングした後、上記撮像領域の境界付近を上記パラメータ推定部によって推定された劣化画像補正パラメータまたは上記修復画像に基づいてフィルタリングすることを特徴とする画像修復装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、手振れやピントのずれ、大気の揺らぎ等による影響で、劣化した静止画像を撮像した時、あるいは動画像の撮像において、動画のある１フレームが同様の影響を受け、劣化した静止画像として切り出された時に、前記劣化した画像に対して、画像を修復するための画像修復装置及び画像修復方法に関するものである。
特に、デジタルスチルカメラやデジタルビデオカメラ、携帯電話などに搭載されているカメラの撮像装置に関するものであるが、これらに限定されるものでは無い。
【背景技術】
【０００２】
手に撮像機器を持ったときに起きる手振れや、揺れた車内から被写体を撮影した場合におきる手振れは、撮像機器が揺れており、撮像された写真が一定の方向に流れて劣化している。また、動いている被写体の流し撮りや、被写体の中の一部に高速で移動するものがある場合、画像の一部が静止、あるいは画像の一部が流れて撮像される。このような劣化を受けた撮像画像は、撮像画像中の劣化している部分を検出し、劣化画像を数学的にモデリングできれば、画像修復アルゴリズムを適用することで劣化を修復することができる。
【０００３】
特別な場合を除いて、劣化画像の数学的画像モデルは、２次元の線形モデルとしてモデル化でき、劣化に関する情報はインパルス応答関数で表現される。つまり画像修復とは、この２次元線形モデルにおけるインパルス応答関数の逆関数による逆フィルタを劣化画像に施すことである。（たとえば、参考文献1）
【０００４】
【非特許文献１】Ｈ．Ｃ．Ａｎｄｒｅｗｓ  ａｎｄ  Ｂ．Ｒ．Ｈｕｎｔ，“Ｄｉｇｉｔａｌ  Ｉｍａｇｅ  Ｒｅｓｔｏｒａｔｉｏｎ”，Ｐｒｅｎｔｉｃｅ  Ｈａｌｌ，１９７７。
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、参考文献１による手法では、画像の境界付近において、実際に劣化を受け撮像される画像と、数学的画像モデルとの間にモデル化誤差があり、精度の高い修復画像が得られないばかりか、モデル化誤差が前記逆フィルタによって増幅され、さらに劣化した画像となる。また、処理に離散フーリエ変換（ＤＦＴ）行列を用いた手法では、定義された領域の境界付近に近似が生じ、近似によるモデル化誤差が生じることで修復画像が大幅に劣化してしまう。
また、ＤＦＴを用いる手法と同様のアプローチに、離散コサイン変換（ＤＣＴ）行列や、離散サイン変換（ＤＳＴ）を用いる手法がある。これらも同様に近似によるモデル化誤差が生じることから修復画像が大幅に劣化してしまう。
【課題を解決するための手段】
【０００６】
本発明に係る画像修復装置は、数学的モデルに基づいて、劣化した原画像を修復する画像修復装置であって、撮像領域の境界部を補正する境界値補正部と、上記境界値補正部によって境界部を補正された画像に基づいて劣化画像補正パラメータを推定するパラメータ推定部と、上記境界値補正部によって境界部を補正された画像または上記パラメータ推定部によって推定された劣化画像補正パラメータに基づいて、原画像を推定し、修復画像を出力する原画像推定部とを有し、上記境界値補正部は、上記撮像領域の境界の外側をパディングした後、上記撮像領域の境界付近を上記パラメータ推定部によって推定された劣化画像補正パラメータまたは上記修復画像に基づいてフィルタリングするものである。
【発明の効果】
【０００７】
本発明に係る画像修復装置は上記のように構成されることで、手振れやぼやけにより劣化して撮像された画像を、数学的画像モデルに適合するように、定義されている領域の境界付近の画素値を補正し画像修復アルゴリズムを適用するので、劣化画像に依存することなく、安定した精度の高い画像修復が行える。
また、ブラインド画像修復にも組み込むことができ、本発明により修復画像の精度向上と共に、数学的画像モデルに含まれるパラメータの推定精度が向上する。
【発明を実施するための最良の形態】
【０００８】
以下、図面を参照して本発明の実施の形態を説明する。以下に説明する実施の形態は、主に、デジタルスチルカメラを用いた撮像において、手振れやぼやけにより劣化して撮像された画像の修復方法として用いることに適したものであるが、他の撮像装置にも適用可能な技術である。また、撮像画像をコンピュータに取り込んで処理することも可能な技術である。
【０００９】
実施の形態１．
図１は本発明における図４の境界値補正部１に含まれるパディング部１１の処理を示した図であり、図２は本発明における図４の境界値補正部１に含まれるフィルタリング部１２の処理を示した図であり、図３は本発明における図４の境界値補正部１に含まれるマージ部１３の処理を示した図である。また、図４は本発明における画像修復アルゴリズムを示すブロック図であり、図５は図４の境界値補正部１のブロック図である。
【００１０】
本発明は、ディジタル信号処理における画像修復アルゴリズムおよびブラインド画像修復アルゴリズムに係りなされたものであり、画像修復アルゴリズムは、得られた画像が手振れやぼやけにより劣化している場合、その劣化に関するインパルス応答関数と付加的雑音の性質を与え、原画像を推定（画像修復、画像復元）するものである。また、ブラインド画像修復アルゴリズムは、得られた劣化画像のみを用い、原画像の推定のほか、劣化に関するインパルス応答関数や付加的雑音の性質をも同時に推定するものである。これは、数学的画像モデルに含まれるパラメータを推定することと同じである。
【００１１】
さらに、原画像をも数学的画像モデルでモデル化し、ブラインド画像修復アルゴリズムを適用することもあり、この場合原画像を示す数学的画像モデルに含まれるパラメータをも推定する。このとき、パラメータの推定精度および原画像の推定精度は、アルゴリズムによって大きく依存すると考えられている。しかし、アルゴリズムが大きな問題ではない。数学的画像モデルに含まれるモデル化誤差が大きな要因となり、パラメータおよび原画像の推定精度を大きく下げてしまう。
以下に、劣化画像の数学的画像モデルを記述し、モデル化誤差と本発明の処理を明確にしたのち、ブロック図と画像の模式図を用いて本発明を詳細に説明する。
【００１２】
【数１】
2007102477.tif 000003
（式１）
【００１３】
式（１）においてｉ、ｊは画素の位置を表しており、画面垂直方向にＭ１画素、画面水平方向にＮ１画素となるサイズの画像であれば、ｉ＝１，・・・，Ｍ１、ｊ＝１，・・・、Ｎ１となる整数値を取る。このとき、画像はＭ１×Ｎ１の領域で定義される。また、この定義された領域は撮像される領域と同じである。ｕｉ，ｊは劣化の無い原画像、ｚｉ，ｊは手振れやぼやけと言った撮像時の劣化に起因して撮像された劣化画像、ηｉ，ｊは撮像画像中に含まれる付加的雑音、ｄｋ，ｌは手振れやぼやけを示すインパルス応答関数、Ｋより計算される（２Ｋ＋１）はインパルス応答関数の画面垂直方向の次数、Ｌより計算される（２Ｌ＋１）はインパルス応答関数の画面水平方向の次数を表すとする。ただし、原画像および劣化画像の平均値は０であるとする。これはそれぞれの撮像画像の画素信号からそれぞれの平均値を、各画素信号から引くことにより計算できる。このとき、原画像および劣化画像の統計的性質として、これらの画像の期待値は平均値と等しい（エルゴード性）と仮定すると、Ｅ［・］を・の期待値を表すものとして、Ｅ［ｕｉ，ｊ］＝０、Ｅ［ｚｉ，ｊ］＝０と表すことができる。つまり、ここでは期待値は平均値と同じである。さらに、付加的雑音ηｉ，ｊにはＥ［ｕｉ，ｊηｉa，ｊa］＝０（ｉaおよびｊaは任意）の性質を仮定する。この仮定は、付加的雑音が原画像にほとんど依存しない事実から導かれた仮定である。画像修復アルゴリズムの中には、付加的雑音の性質としてＥ［ηｉ，ｊ］＝０、Ｅ［ηｉ，ｊηｉa，ｊa］＝σ２δｉ−ｉａ，ｊ−ｊａ（δ：クロネッカーのデルタ関数）を仮定するものや、劣化画像中に含まれる付加的雑音の自乗和をある値で仮定するもの、雑音を無視できるものもある。しかしながら、本発明は雑音の仮定によって限定されない。
【００１４】
式（１）を基にした画像修復アルゴリズムあるいはブラインド画像修復アルゴリズムは、莫大な数が提案されており、枚挙にいとまが無い。しかし、それらのアルゴリズムでは、定義された領域の外側の領域に対して、原画像が０であると仮定されていたり、巡回信号が続くと仮定してアルゴリズムに離散フーリエ変換行列（ＤＦＴ）を用いていたり、定義された領域の外側の領域に対して、鏡像の関係で原画像が続いていると仮定して離散コサイン変換（ＤＣＴ）行列もしくは離散サイン変換行列（ＤＳＴ）を用いたりしており、領域の外側に存在する原画像はアルゴリズムに適合するように仮定されている。しかし、領域の外側に存在している原画像は０値を取るわけではなく、領域内に存在する撮像画像が領域外に存在するわけでも無いことから、上記アプローチを適用することでモデル化誤差が生じ、良い修復画像を得ることができない。
また、領域内の劣化画像は、上記領域の外側にある原画像が、劣化の際に領域内に写り込むことになるが、式（１）では領域の外側にある原画像の値は定義されておらず、０値を取るものとして扱われることから、やはりモデル化誤差が生じる。
【００１５】
これまで提案されてきた画像処理アルゴリズムにおける数値実験では、上記モデル化誤差は式（１）における付加的雑音ηｉ，ｊに含まれて扱われており、付加的雑音の性質に誤差を生じている。また、画像修復アルゴリズムはインパルス応答関数の逆関数を用いた逆フィルタリングであるので、画像修復時に付加的雑音を増幅させてしまう欠点がある。つまり、画像修復時にモデル化誤差によって付加的雑音が増幅されて、修復画像中に大きな雑音を含んでしまう結果となる。
上記の境界付近の扱いと修復画像の精度を解析した研究がほとんどなされておらず、これまで発表された画像修復アルゴリズムおよびブラインド画像修復アルゴリズムは、アルゴリズムの性能が正しく評価されず、劣化画像に依存した結果が得られているに過ぎない。そこで本発明では、上記モデル化誤差による影響を抑えるため、劣化画像の境界値補正処理を組み込んだ。
以下に、模式図を用い、モデル化誤差を抑えるために組み込んだ劣化画像の境界値補正処理について詳細に述べる。
【００１６】
まず、図１において、内側の枠で示される領域は撮像される領域であり、数学的画像モデルが定義される領域でもある。撮像画像が劣化した場合、手振れやぼやけによりフィルタリングされた結果が撮像されるので、外側の枠で示した領域までの原画像が撮像される領域まで写り込んで来ることとなる。このとき、写り込まれる可能性のある領域を示したのが図２の斜線部分である。
しかし、数学的画像モデルでは、定義された領域にのみ原画像が定義されており、その領域の外側は便宜的に零値や平均値を取ると仮定されることから、図２の斜線で示した領域に対しては、実際に撮像した画像と異なる値を持つこととなる。これがモデル化誤差の要因である。モデル化誤差を抑えるためには、数学的画像モデルの境界付近における記述を正確にする必要があるが、定義される領域よりも外側から写り込まれるということを、定義される領域のみで正確に記述することは不可能である。したがて、定義される領域よりも外側において、単に便宜的に零値等を取ると仮定したのでは、この誤差を抑えることができない。そこで、定義される領域を拡張し、撮像画像を数学的画像モデルに適合するよう境界値補正処理を施す。
【００１７】
数学的画像モデルにおいて、定義された領域よりも外側にある原画像は、零値や平均値を取るものと仮定されることから、この仮定に適合するようにできるだけ劣化画像の境界値を補正する。いま、原画像の推定値が得られていない場合（アルゴリズムの初期段階）、劣化画像を原画像の初期の推定値（反復アルゴリズムを実行する際の初期値であり、以下原画像の初期推定値と呼ぶこととする。）とする。フィルタリング領域のサイズが（２Ｋ＋１）×（２Ｌ＋１）であることを考慮して、原画像の推定値の上端と下端、左端と右端を図３で示す外側の枠で示される領域へ拡張し（以下、拡張領域と呼ぶ）、劣化画像の上端および下端に２Ｋ行の零値を、左端と右端に２Ｌ列の零値を埋める。以下、零値を埋めることで説明を展開するが、劣化画像の平均値など、得られている情報から定数を定め、その値を埋めても処理は変わらない。拡張領域を２Ｋ行、２Ｌ列としたのは、フィルタリング領域の外側の境界での値が零とならなければならないという数学的画像モデルからの要請を満足するためである。
【００１８】
このとき、原画像の推定値とパラメータの推定値が真値であれば、定義される領域内をインパルス応答関数でフィルタリングして劣化画像を計算すれば、モデル化誤差を無くした補正値を得ることができる。
しかし、原画像の推定値およびパラメータの推定値が真値である保障が無いことから、拡張領域においてモデル化誤差の無い劣化画像を作成すべく、図４の斜線部分のみフィルタリングを行い、図５の斜線部分へ劣化画像を埋める。これにより、できるだけモデル化誤差を無くした境界値補正劣化画像を得ることができる。
【００１９】
以上の処理により、定義される領域の外側に零値が存在すると仮定されたもとでの劣化画像を作成することができた。本処理をブラインド画像修復アルゴリズムに組み込むことによって、画像修復アルゴリズムにおける境界値に対する問題が解決され、パラメータの推定値および原画像の推定値の精度を向上させることができる。
以下に、図６のブロック図を用いて、本発明を詳細に述べる。
【００２０】
図６は、この発明を実施するための実施の形態１における画像修復装置１０１の構成を示すブロック図、図７は画像修復装置の境界値補正部１の構成を示すブロック図である。
【００２１】
まず、図６の画像修復装置１０１において、入力端子１１１から劣化画像を、入力端子１１２からは数学的画像モデルに含まれる未知パラメータを推定するための初期推定値を入力し、出力端子１１３から原画像の推定値（修復画像）が、出力端子１１４からは数学的画像モデルに含まれる未知パラメータの推定値が出力される。
【００２２】
画像修復装置１０１には、数学的画像モデルに含まれる未知パラメータを推定するパラメータ推定部３と、原画像の推定を行う原画像推定部４があり、これらの推定にはディジタル信号処理を用いたブラインド画像修復アルゴリズムを用いる。従来から提案されているアルゴリズムでは、パラメータの推定と原画像の推定のどちらにも劣化画像を入力しているが、本発明の主目的のひとつは、画像修復装置に数学的画像モデルのモデル化誤差の補正処理を組み込むことであるので、境界値補正部１で劣化画像の境界付近の値を補正して、パラメータ推定部３あるいは原画像推定部４へ入力する。スイッチ２は、パラメータの推定を行う場合と、原画像の推定を行う場合で切り替える。
【００２３】
スイッチ２は、採用するブラインド画像修復アルゴリズムによって動作が異なる。例えば、パラメータを推定に逐次予測誤差法を用い、原画像の推定にカルマン・固定区間スムーザを用いる場合、スイッチ２は左側が先にＯＮとなりパラメータ推定部３でパラメータを推定し、その後右側をＯＮにして、推定されたパラメータを用い原画像推定部４にて画像修復を行う。逐次予測誤差法を用いたパラメータの推定方法とカルマン・固定区間スムーザによる原画像の推定方法についてはのちに詳細を述べる。
パラメータ推定アルゴリズムおよび原画像の推定アルゴリズムは、莫大な数の研究がなされてきているが、境界値補正部１を用いた本発明では、これらのアルゴリズムによって限定されるものでは無い。
【００２４】
また、上記と異なるアルゴリズムとしては、先にスイッチ２の右側をＯＮとし、原画像推定部４において原画像の推定を行った後、スイッチ２の左側をＯＮとし、パラメータ推定部３にてパラメータの推定を行うアルゴリズムがある。この場合、原画像の推定とパラメータの推定は、どちらとも収束するまで交互に繰り返される。つまりスイッチ２は左右交互にＯＮとなる。このような解法の一つとしてＥＭアルゴリズム（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ）と呼ばれる最尤推定法が挙げられるが、その他にもパラメータ推定を行う信号処理アルゴリズム、システム同定法がたく存在している。
【００２５】
また、原画像推定部４において採用される画像修復アルゴリズムも枚挙にいとまが無い。それらの手法の大部分は式（１）を基礎としており、式（１）を用いた手法および式（１）を変形して導いた手法による画像復元アルゴリズムであり、いずれも画像復元を行うには式（１）におけるインパルス応答関数の推定が必要となる。インパルス応答関数が推定された下での画像復元アルゴリズムの代表的なものを以下に述べる。
【００２６】
まず、インパルス応答関数の逆関数を用いた逆フィルタリング手法がある。また、逆関数を持たず求めることができない場合、擬似逆または一般逆と呼ばれる逆関数を求める逆フィルタリング手法がある。さらに、実施の形態１で用いた拘束条件付最小自乗法における特徴抽出オペレータを変更したものや、重み付き拘束条件最小自乗法として導出されたアルゴリズムがある。これらの手法は一般的に正規直交変換を用いて空間周波数領域で実行される。その正規直交変換にはＤＦＴ、ＤＣＴ、ＤＳＴが用いられることが多い（例えばＤＣＴを用いる場合は、ＤＣＴ空間周波数領域と呼ぶことが多い。）。また、これらのアルゴリズムは全て、空間領域で反復処理として実行することもできる。
【００２７】
次に、撮像された画像の原画像に統計的性質を仮定し、その性質に基づき画像復元を行うアルゴリズムがある。その代表的なものにウィーナ・フィルタがあり、空間領域、空間周波数領域どちらでも実行できる。また、ウィーナ・フィルタの逐次推定手法としてカルマン・フィルタによる画像復元アルゴリズムがある。本手法は式（１）の下に空間領域、周波数領域、空間周波数領域で行うアルゴリズムを導出できる。また、カルマン・フィルタは逐次推定手法であることから、その初期値に依存する場合もあり、その初期値を推定するために、固定区間スムーザと呼ばれるカルマン・フィルタ実行後に行うアルゴリズムが存在する。これにより、適切な初期値を推定し、カルマン・フィルタを実行し、再び固定区間スムーザを実行することで最良の復元画像を得ることもできる。また、固定区間スムーザの変わりに、固定ラグスムーザ、固定点スムーザを用いても、固定区間スムーザとほぼ同等の画像復元が行える。また、これらのアルゴリズムは事後確率最大化（Ｍａｘｉｍｕｍ  Ａ  Ｐｏｓｔｅｒｉｏｒｉ：ＭＡＰ）推定を用いた画像復元法の一つとも考えられる。また、ＭＡＰ推定には擬似焼き鈍し（Ｓｉｍｕｌａｔｅｄ  Ａｎｅａｌｉｎｇ）法による確率緩和法を用いた手法なども用いられる。
上記以外にも射影フィルタを中心としたアルゴリズムや、エントロピー最大化手法に基づいたアルゴリズムなど、多くの手法が存在する。本発明においては、画像復元アルゴリズムとして上記に示した様々な手法を用いることが可能であり、応用目的や、使用する電子回路の性能やメモリ量等に応じて最適な手法を選択可能である。
【００２８】
さらに、本発明における図６の画像修復装置１０１において、パラメータ制御部５を用いることで、より複雑な反復処理を組み込んだパラメータの推定および原画像の推定を可能とした。パラメータ制御部５は境界値補正部１、パラメータ推定部３、原画像推定部４へのパラメータの受け渡しのタイミングを制御するものであり、これについて以下に述べる。
【００２９】
パラメータ推定アルゴリズムとして、パラメータの初期推定値（初期値）を必要とするアルゴリズムを採用している場合、パラメータ推定部３で推定されるパラメータの推定値は、パラメータ制御部５へ入力されるが、すぐにパラメータ推定部３へ戻すことでこの値をパラメータの初期推定値とし、再びパラメータ推定を行うことができる。本処理を一定回数行えば、パラメータの推定精度が向上する。尚、パラメータの初期推定値の与え方は、パラメータの推定方法に逐次予測誤差法を用いた場合のアルゴリズムを詳細に後述するためそこで同時に示す。
【００３０】
また、パラメータ推定部３で推定されるパラメータの推定値をパラメータ制御部５へ入力し、パラメータ制御部５から境界値補正部１へ出力すれば、精度が高くなった新しいパラメータを用いて境界値補正劣化画像を算出することができる。従って、再びパラメータ推定部３においてパラメータ推定を行えば、パラメータの推定精度が向上する。本処理を一定回数行えば、パラメータの推定精度が向上する。
【００３１】
また、原画像推定部４にて画像修復を実行する際には、数学的画像モデルのパラメータが必要となることから、パラメータ制御部５より精度の高いパラメータの推定値が原画像推定部４に出力される。原画像推定部４からは修復画像が出力端子１１３へ出力されるだけでなく、境界値補正部１へも出力される。境界値補正部１では、入力端子１１１から入力される劣化画像と原画像推定部４から入力される修復画像、パラメータ制御部５から入力されるパラメータの推定値を用い、新たな境界値補正劣化画像を出力する。
【００３２】
上記のように、パラメータの推定値や原画像の推定値が得られる度に、境界値補正部１で境界値補正劣化画像を計算すれば、さらにパラメータの推定値や原画像の推定値の精度が向上する。このフィードバックの処理回数は特に決まっているわけでなく、前もって決定しておくことでも構わない。また、パラメータ制御部５においてパラメータの精度をモニタリングすることや、原画像推定部４から出力される修復画像の精度をモニタリングすることなどをアルゴリズムに組み込むこともできる。様々な手順が考えられるが、これらはどれも本発明に含まれる。
次に、図７を用いて境界値補正部１における劣化画像の境界付近の補正について詳細に述べる。
【００３３】
図７に示す境界値補正部１において入力端子２１からは画像修復装置１０１における入力端子１１１から入力される劣化画像を、入力端子２２からは原画像推定部４から出力される修復画像を、入力端子２３からはパラメータ制御部５から出力されるパラメータの推定値を入力する。そして出力端子２４から境界値補正劣化画像を出力する。
境界値補正劣化画像を算出する際には、入力端子２２から入力される修復画像と、入力端子２１から入力される劣化画像（原画像の初期推定値）のどちらかをスイッチ１４で選択し、パディング部１１で選択された画像に対して、図３の斜線で示されている領域に零埋めを行う（前述したが、零値ではなく劣化画像の平均値や、その他一定の値を埋めても良い）。零埋めされた画像はフィルタリング部１２に出力され、フィルタリング部１２では、図４の斜線で示された領域のみ、入力端子２３から入力されたパラメータの推定値におけるインパルス応答関数でフィルタリングする。フィルタリング部１２から出力された画像はマージ部１３へ入力され、図５の斜線で示された領域に入力端子２１から入力した劣化画像を埋め、出力端子２４へ出力される。最後に、出力端子２４から境界値補正劣化画像が出力される。
【００３４】
ここで、スイッチ１４は、修復画像の精度をモニタリングすることで、アルゴリズム実行中に下側をＯＮとする処理も考えられる。ここで修復画像の精度は、例えば推定されたパラメータｈｋ，ｌを用いて原画像の推定値ｘｉ，ｊを畳み込み演算で劣化させ、中心部（Ｍ−２Ｋ）×（Ｎ−２Ｌ）に対して劣化画像との誤差の自乗和を計算し、この誤差が原画像の初期推定値を用いた場合に計算される誤差の自乗和との大小から決めることができる。しかし、本発明はスイッチ１４の処理方法に限定されるものではない。
以下に、図３乃至図７を用い、逐次予測誤差法によるパラメータ推定を用いた場合の信号の流れから発明を詳細に説明する。アルゴリズムには反復処理を含むことができ、反復回数は自由に設定して構わない。いま説明のために、パラメータ推定部３のパラメータ制御部５との反復回数をｐ１回、パラメータ推定部３とパラメータ制御部５、境界値補正部１との反復回数をｐ２回、原画像推定部４内で行われる原画像の反復推定回数をｑ１回、原画像推定部４と境界値補正部１との反復回数をｑ２回、パラメータ推定を含んだ原画像推定部４と境界値補正部１との反復をｒ回とする。尚、反復の順序に関しては以下に詳述する。
【００３５】
まず、図６で示される画像修復装置１０１において、入力端子１１１から劣化画像を入力し、境界値補正部１へ入力する。入力端子１１２からはパラメータの初期推定値（初期値）を入力し、パラメータ制御部５はパラメータの初期推定値を境界値補正部１へ出力する。
【００３６】
図６で示される境界値補正部１は、この段階では入力端子２２から修復画像が得られないので、スイッチ１４は下側がＯＮとなり、入力端子２１から入力された劣化画像を原画像の初期推定値とし、パディング部１１へ入力される。パディング部１１では図３に示される斜線部分へ零値を埋め、フィルタリング部１２へ出力する。フィルタリング部１２では、入力端子２３から入力されたパラメータの推定値（この段階では初期推定値）を入力し、パラメータの推定値に含まれているインパルス応答関数の推定値で、図４で示される斜線部分のみフィルタリングを行う。さらにマージ部１３では、図５で示される斜線部分に、入力端子２１から入力された劣化画像を埋め込み、出力端子２４から境界値補正画像を出力する。
【００３７】
次に、逐次予測誤差法を実行するので、スイッチ２は左側がＯＮとなり、境界値補正劣化画像がパラメータ推定部３へ入力する。パラメータ推定部３では、逐次予測誤差法を実行し、パラメータの推定値をパラメータ制御部５へ出力する。この時点で１回目のパラメータ推定値が得られたこととなる。パラメータ制御部５は、１回目のパラメータ推定値をパラメータ推定部へ出力して、パラメータ推定部３では１回目のパラメータ推定値を初期推定値として逐次予測誤差法を実行し、パラメータ推定が行われる。このような反復処理をｐ１回行ったのちに得られたｐ１回目のパラメータ推定値を、パラメータ制御部５から境界値補正部１とパラメータ推定部３へ出力する。
【００３８】
境界値補正部１では、ｐ１回目のパラメータの推定値を用いて先に示した処理を行い、境界値補正劣化画像を出力する。境界値補正部１を含めたパラメータ推定の反復はｐ２回であるので、パラメータ制御部５にはｐ１×ｐ２回目のパラメータ推定値が最新となる。この推定値を境界値補正部１と原画像推定部４へ出力し、スイッチ２の右側をＯＮにして、境界値補正部１から出力された境界値補正劣化画像から画像修復を行う。
【００３９】
原画像推定部４では、スイッチ２を介して境界値補正劣化画像を入力し、ｐ１×ｐ２回目のパラメータ推定値を用い、先に述べた画像修復アルゴリズムを用いて、画像修復を行う。画像修復アルゴリズムによっては反復推定が可能なものもあり、実施の形態１．ではｑ１回の反復推定を行うとしている。ｑ１回目に修復された画像は境界値補正部１と出力端子１１３へ出力される。境界値補正部１では、図７の入力端子２２から１回目の修復画像が入力され、スイッチ１４の上側をＯＮとすることにより、前述した手順で境界値補正劣化画像を出力端子２４から出力する。出力された境界値補正劣化画像はスイッチ２を介して原画像推定部４に入力され、画像修復が行われる。このような画像修復をｑ２回繰り返す。画像修復をｑ１×ｑ２回繰り返した時点で、再びパラメータ推定をｐ１×ｐ２回繰り返し、この反復で得られたパラメータの推定値を用いて再びｑ１×ｑ２回の画像修復を行う。これらの処理をｒ回繰り返して得られたｑ１×ｑ２×ｒ回目の修復画像を最終的に得られる原画像の推定値とし、出力端子１１３から出力する。同時にｐ１×ｐ２×ｒ回目のパラメータの推定値をパラメータの最終推定値として出力端子１１４から出力する。
以上の逐次予測誤差法を用いたブラインド画像修復アルゴリズムを、フロー図として図８を参照し、詳細に述べる。
【００４０】
まず、劣化画像と数学的画像モデルに含まれるパラメータの初期推定値を入力する（Ｓ１、Ｓ２）。劣化画像（原画像の初期推定値）とパラメータの初期推定値を用い、境界値補正処理（Ｓ３）を行う。境界値補正劣化画像を用いて逐次予測誤差法によりパラメータの推定（Ｓ４）を行う。得られたパラメータの推定値を初期値として、再び逐次予測誤差法によりパラメータの推定（Ｓ４）を行う。このようなパラメータの反復推定をｐ１回行う。
【００４１】
次に、ｐ１回目のパラメータの推定値を用いて境界値補正処理（Ｓ３）を行う。境界値補正劣化画像が計算されるので、パラメータの推定（Ｓ４）をｐ１回反復し、再び境界値補正処理（Ｓ３）を行う。このような境界値補正処理をｐ２回行う。
【００４２】
ｐ１×ｐ２回目に推定されたパラメータを用いて境界値補正処理（Ｓ５）を行い、境界値補正劣化画像を入力して、原画像の推定（Ｓ６）を行う。推定された原画像を再び初期値として用い、原画像の推定（Ｓ６）を行う。このような原画像の反復推定をｑ１回行う。ｑ１回目に得られた原画像の推定値を用いて境界値補正処理（Ｓ５）を行い、原画像の推定（Ｓ６）を行う。境界値補正処理を含めた原画像の反復推定をｑ２回行う。
【００４３】
このように得られたｑ１×ｑ２回目の原画像の推定値を用いて境界値補正処理（Ｓ３）を行い、パラメータの推定を再び行う。
以上のような反復を行うことにより、最終的に得られたｐ１×ｐ２×ｒ回目のパラメータの推定値と、ｑ１×ｑ２×ｒ回目の原画像の推定値を出力し、修復画像を得る。
次に、逐次予測誤差法によるパラメータ推定方法と、カルマン・固定区間スムーザを用いた原画像の推定方法について詳細に述べる。
【００４４】
逐次予測誤差法およびカルマン・固定区間スムーザは、式（１）中の原画像ｕｉ，ｊに対して統計的数学モデルと呼ばれる数学的画像モデルを導入する必要がある。いま、境界値補正部１により画像サイズが拡張されたことから、その縦サイズをＭ、横サイズをＮと表記し、原画像を示す統計的数学モデルに半因果的確率モデルを導入すると、原画像ｕｉ，ｊは、
【００４５】
【数２】
2007102477.tif 000004
（式２）
ただし，ξi,jはシステム雑音と呼ばれ、平均
【００４６】
【数３】
2007102477.tif 000005
および分散
【００４７】
【数４】
2007102477.tif 000006
の性質を仮定している。このとき、式（２）および式（１）をベクトル・行列表現すれば、
【００４８】
【数５】
2007102477.tif 000007
【００４９】
【数６】
2007102477.tif 000008
で表現できる。ただし、
【００５０】
【数７】
2007102477.tif 000009
【００５１】
【数８】
2007102477.tif 000010
【００５２】
【数９】
2007102477.tif 000011
【００５３】
【数１０】
2007102477.tif 000012
【００５４】
【数１１】
2007102477.tif 000013
【００５５】
【数１２】
2007102477.tif 000014
【００５６】
【数１３】
2007102477.tif 000015
でそれぞれ定義する。
これらの行列はサイズが大きいことから計算機上で演算することは困難であるが、巡回型の行列として近似できれば離散フーリエ変換行列で対角化でき、画像をｉ行ごとに独立して扱うことができる。上式の巡回行列は、それぞれ
【００５７】
【数１４】
2007102477.tif 000016
【００５８】
【数１５】
2007102477.tif 000017
【００５９】
【数１６】
2007102477.tif 000018
として表現でき、さらにこれらを離散フーリエ変換行列で対角化すれば、式（２）および式（１）は、
【００６０】
【数１７】
2007102477.tif 000019
【００６１】
【数１８】
2007102477.tif 000020
（式３）
で記述することができる．ただし、離散フーリエ変換行列を
【００６２】
【数１９】
2007102477.tif 000021
で定義し、
【００６３】
【数２０】
2007102477.tif 000022
【００６４】
【数２１】
2007102477.tif 000023
【００６５】
【数２２】
2007102477.tif 000024
【００６６】
【数２３】
2007102477.tif 000025
を定義した。右肩に示す添え字（ｉ）は、行を示しており、
【００６７】
【数２４】
2007102477.tif 000026
（式４）
【００６８】
【数２５】
2007102477.tif 000027
(式５)
【００６９】
【数２６】
2007102477.tif 000028
(式６)
として表現される。
【００７０】
ここで、原画像の推定のために必要な未知パラメータは式（４）乃至式（６）に含まれており、逐次予測誤差法では
【００７１】
【数２７】
2007102477.tif 000029
と
【００７２】
【数２８】
2007102477.tif 000030
を推定することになる。
以下に、この逐次予測誤差法について詳述する。
【００７３】
【数２９】
2007102477.tif 000031
【００７４】
【数３０】
2007102477.tif 000032
を定義し、ｉ行ごと画像モデルをＡＲＭＡ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデルによる表現をすると、
【００７５】
【数３１】
2007102477.tif 000033
に対して、
【００７６】
【数３２】
2007102477.tif 000034
と表現できる。このとき、本来の未知パラメータ
【００７７】
【数３３】
2007102477.tif 000035
を直接求めるのではなく、
【００７８】
【数３４】
2007102477.tif 000036
を推定した後に、最適化問題を定式化し、解くことにより求める。
（６）式において、行を示す変数の右肩の添え字（ｉ）を省略し、表記の便宜性から、
【００７９】
【数３５】
2007102477.tif 000037
【００８０】
【数３６】
2007102477.tif 000038
【００８１】
【数３７】
2007102477.tif 000039
を定義すると、
【００８２】
【数３８】
2007102477.tif 000040
（式７）
を得る。ここで，z-1は遅延演算子を表すとし、
【００８３】
【数３９】
2007102477.tif 000041
【００８４】
【数４０】
2007102477.tif 000042
を定義すると、
【００８５】
【数４１】
2007102477.tif 000043
の表現を得ることができる。また、次数
【００８６】
【数４２】
2007102477.tif 000044
と２Ｌのうち大きい方を
【００８７】
【数４３】
2007102477.tif 000045
と定義すると、観測画像yjは次数をＲとし、平均が
【００８８】
【数４４】
2007102477.tif 000046
、分散が
【００８９】
【数４５】
2007102477.tif 000047
である正規性白色雑音ejを用いて
【００９０】
【数４６】
2007102477.tif 000048
として表現できる。ejはyjの一段予測誤差（Ｉｎｎｏｖａｔｉｏｎ過程）と呼ばれる。      
【００９１】
【数４７】
2007102477.tif 000049
は安定な多項式であり、
【００９２】
【数４８】
2007102477.tif 000050
を定義することにより、
【００９３】
【数４９】
2007102477.tif 000051
として表すこともできる。さらに、yjのスペクトル密度の関係から、多項式Ｃ(z-1)は、
【００９４】
【数５０】
2007102477.tif 000052
により一意に決定される。このとき、（７）式の対数尤度関数は、未知パラメータを
【００９５】
【数５１】
2007102477.tif 000053
として、
【００９６】
【数５２】
2007102477.tif 000054
（式８）
で表現することができ、θに対する最大化を行うことで、パラメータ推定を行うことができる。ここで、次のような変形を用いて未知パラメータ数を減らし、パラメータの推定を簡潔に行う。
（８）式の両辺を
【００９７】
【数５３】
2007102477.tif 000055
で割り、
【００９８】
【数５４】
2007102477.tif 000056
【００９９】
【数５５】
2007102477.tif 000057
を定義すると、
【０１００】
【数５６】
2007102477.tif 000058
（式９）
を得ることができる。ここで、未知パラメータを、
【０１０１】
【数５７】
2007102477.tif 000059
と再定義すれば、予測誤差規範（Ｐｒｅｄｉｃｔｉｏｎ  Ｅｒｒｏｒ  Ｃｒｉｔｅｒｉｏｎ）
【０１０２】
【数５８】
2007102477.tif 000060
の最小化により推定値
【０１０３】
【数５９】
2007102477.tif 000061
が得られ、
【０１０４】
【数６０】
2007102477.tif 000062
により一段予測誤差の分散が求まる。さらにこれらが求まれば、（９）式にスペクトル因子分解法を適用して一段予測誤差係数
【０１０５】
【数６１】
2007102477.tif 000063
および
【０１０６】
【数６２】
2007102477.tif 000064
を求め、状態空間モデルのシステム雑音と付加的雑音それぞれの分散を
【０１０７】
【数６３】
2007102477.tif 000065
【０１０８】
【数６４】
2007102477.tif 000066
として求めることができる。スペクトル因子分解法については後述する。
次に、未知パラメータの推定値
【０１０９】
【数６５】
2007102477.tif 000067
の推定誤差共分散を
【０１１０】
【数６６】
2007102477.tif 000068
と定義し、初期値を
【０１１１】
【数６７】
2007102477.tif 000069
、
【０１１２】
【数６８】
2007102477.tif 000070
として与えることにより、未知パラメータは逐次的に、
【０１１３】
【数６９】
2007102477.tif 000071
【０１１４】
【数７０】
2007102477.tif 000072
【０１１５】
【数７１】
2007102477.tif 000073
として求めることができる。ただし、
【０１１６】
【数７２】
2007102477.tif 000074
で定義されており、逐次予測誤差法を実行するには、ej、Ψjを推定することが必要となる。これについて以下に述べる。
ejは多項式{Cr}が得られていれば、
【０１１７】
【数７３】
2007102477.tif 000075
を用いて得ることができる。この両辺をθで微分すると、
【０１１８】
【数７４】
2007102477.tif 000076
となるので、
【０１１９】
【数７５】
2007102477.tif 000077
【０１２０】
【数７６】
2007102477.tif 000078
を定義すると、
【０１２１】
【数７７】
2007102477.tif 000079
で表現することができる。すなわちΨjは、
【０１２２】
【数７８】
2007102477.tif 000080
（式１０）
として表現できる。ただし、
【０１２３】
【数７９】
2007102477.tif 000081
を定義した。
次に、
【０１２４】
【数８０】
2007102477.tif 000082
は、
【０１２５】
【数８１】
2007102477.tif 000083
を定義し、
【０１２６】
【数８２】
2007102477.tif 000084
と表現すると、
【０１２７】
【数８３】
2007102477.tif 000085
（式１１）
なる関係があるので、
【０１２８】
【数８４】
2007102477.tif 000086
【０１２９】
【数８５】
2007102477.tif 000087
を定義して、
【０１３０】
【数８６】
2007102477.tif 000088
を求める。いま、
【０１３１】
【数８７】
2007102477.tif 000089
の係数を
【０１３２】
【数８８】
2007102477.tif 000090
として表すこととし、
【０１３３】
【数８９】
2007102477.tif 000091
を定義すると、
【０１３４】
【数９０】
2007102477.tif 000092
で表すことができる。ここに、
【０１３５】
【数９１】
2007102477.tif 000093
なる関係があることから
【０１３６】
【数９２】
2007102477.tif 000094
を得ることができ、
【０１３７】
【数９３】
2007102477.tif 000095
と
【０１３８】
【数９４】
2007102477.tif 000096
の定義から
【０１３９】
【数９５】
2007102477.tif 000097
を得ることができる。また、
【０１４０】
【数９６】
2007102477.tif 000098
として、
【０１４１】
【数９７】
2007102477.tif 000099
【０１４２】
【数９８】
2007102477.tif 000100
、
【０１４３】
【数９９】
2007102477.tif 000101
、
【０１４４】
【数１００】
2007102477.tif 000102
と同様な形で
【０１４５】
【数１０１】
2007102477.tif 000103
、
【０１４６】
【数１０２】
2007102477.tif 000104
、
【０１４７】
【数１０３】
2007102477.tif 000105
、
【０１４８】
【数１０４】
2007102477.tif 000106
を定義するとκは、
【０１４９】
【数１０５】
2007102477.tif 000107
として表現することができ、
【０１５０】
【数１０６】
2007102477.tif 000108
は、
【０１５１】
【数１０７】
2007102477.tif 000109
で表される。したがって、
【０１５２】
【数１０８】
2007102477.tif 000110
（式１２）
が得られる。また、（６）式中における、
【０１５３】
【数１０９】
2007102477.tif 000111
、
【０１５４】
【数１１０】
2007102477.tif 000112
の具体的な形は、
【０１５５】
【数１１１】
2007102477.tif 000113
、
【０１５６】
【数１１２】
2007102477.tif 000114
、
【０１５７】
【数１１３】
2007102477.tif 000115
の定義により決定される。本論文では、
【０１５８】
【数１１４】
2007102477.tif 000116
（式１３）
【０１５９】
【数１１５】
2007102477.tif 000117
（式１４）
で与えられる。また、
【０１６０】
【数１１６】
2007102477.tif 000118
であるので、最終的にΨjは式（１０）、（１１）、（１２）、（１３）、（１４）式を用いて算出することができる。
ここで、逐次予測誤差法におけるパラメータ推定時の初期推定値について述べる。
【０１６１】
【数１１７】
2007102477.tif 000119
での
【０１６２】
【数１１８】
2007102477.tif 000120
、
【０１６３】
【数１１９】
2007102477.tif 000121
、
【０１６４】
【数１２０】
2007102477.tif 000122
、
【０１６５】
【数１２１】
2007102477.tif 000123
、ｙj（の平均値と同じ値で良い）を与え、初期値ＰＬは対角要素に１〜１００程度の値を入れ、初期値
【０１６６】
【数１２２】
2007102477.tif 000124
においては、インパルス応答関数にはｈ0,0のみ１を入れそれ以外の要素については０を入れたものから
【０１６７】
【数１２３】
2007102477.tif 000125
を算出し、半因果的確率モデルに含まれる予測係数は経験的に
【０１６８】
【数１２４】
2007102477.tif 000126
を入れたものから
【０１６９】
【数１２５】
2007102477.tif 000127
を算出する。また、
【０１７０】
【数１２６】
2007102477.tif 000128
は
【０１７１】
【数１２７】
2007102477.tif 000129
に１０の−３〜−４乗の値を入れ計算したものを初期値とし、
【０１７２】
【数１２８】
2007102477.tif 000130
については
【０１７３】
【数１２９】
2007102477.tif 000131
に１０の−４〜−５乗の値を入れ計算したものを初期値とすれば良い。
以上により求まった未知パラメータ
【０１７４】
【数１３０】
2007102477.tif 000132
を用い、最適化問題を定式化することによって
【０１７５】
【数１３１】
2007102477.tif 000133
を推定する手法について述べる。
上述したアルゴリズムによって推定値として得られるパラメータは、
【０１７６】
【数１３２】
2007102477.tif 000134
、
【０１７７】
【数１３３】
2007102477.tif 000135
、
【０１７８】
【数１３４】
2007102477.tif 000136
、
【０１７９】
【数１３５】
2007102477.tif 000137
であり、これらから
【０１８０】
【数１３６】
2007102477.tif 000138
、
【０１８１】
【数１３７】
2007102477.tif 000139
、
【０１８２】
【数１３８】
2007102477.tif 000140
、
【０１８３】
【数１３９】
2007102477.tif 000141
を求めるために、以下の評価規範を最小化する。
【０１８４】
【数１４０】
2007102477.tif 000142
は、
【０１８５】
【数１４１】
2007102477.tif 000143
の関数であり、評価関数
【０１８６】
【数１４２】
2007102477.tif 000144
を、
【０１８７】
【数１４３】
2007102477.tif 000145
として定義し、パラメータ
【０１８８】
【数１４４】
2007102477.tif 000146
について最小化を行う。最小化のための初期値は、逐次予測誤差法で用いた初期値と同じで良い。
次に、得られた
【０１８９】
【数１４５】
2007102477.tif 000147
を用いて、
【０１９０】
【数１４６】
2007102477.tif 000148
を計算し、
【０１９１】
【数１４７】
2007102477.tif 000149
の関係より擬似逆行列を求めてシステム雑音の分散である
【０１９２】
【数１４８】
2007102477.tif 000150
を求める。
また、インパルス応答関数の推定は、
【０１９３】
【数１４９】
2007102477.tif 000151
なる関係式を用い、擬似逆行列により推定する。
最後に、観測雑音の分散は
【０１９４】
【数１５０】
2007102477.tif 000152
である
【０１９５】
【数１５１】
2007102477.tif 000153
についての算術平均として計算する。
次に、スペクトル因子分解について述べる。式で、右辺が得られている下で左辺を求めるために本手法を用いる。右辺の
【０１９６】
【数１５２】
2007102477.tif 000154
の係数が
【０１９７】
【数１５３】
2007102477.tif 000155
として得られているとする。このとき、（９）式の
【０１９８】
【数１５４】
2007102477.tif 000156
の係数の関係をベクトル・行列表現すれば、
【０１９９】
【数１５５】
2007102477.tif 000157
が成り立つ。ベクトル関数
【０２００】
【数１５６】
2007102477.tif 000158
を定義し、両辺を
【０２０１】
【数１５７】
2007102477.tif 000159
で微分すると、
【０２０２】
【数１５８】
2007102477.tif 000160
であるので、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法により
【０２０３】
【数１５９】
2007102477.tif 000161
を求める反復式は、ｎ回目の推定値を
【０２０４】
【数１６０】
2007102477.tif 000162
と表すと、初期値
【０２０５】
【数１６１】
2007102477.tif 000163
を与えることにより、
【０２０６】
【数１６２】
2007102477.tif 000164
として求めることができる。初期値は
【０２０７】
【数１６３】
2007102477.tif 000165
の第一要素を比較的大きな値（１００等）とし、他の要素を０で与えることで解を得ることができる。
【０２０８】
以下に、これまで述べたアルゴリズムをまとめる。
Ｓｔｅｐ１．劣化画像データ
【０２０９】
【数１６４】
2007102477.tif 000166
をＤＦＴ行列
【０２１０】
【数１６５】
2007102477.tif 000167
により変換し、
【０２１１】
【数１６６】
2007102477.tif 000168
を計算して、初期設定として
【０２１２】
【数１６７】
2007102477.tif 000169
での
【０２１３】
【数１６８】
2007102477.tif 000170
、
【０２１４】
【数１６９】
2007102477.tif 000171
、
【０２１５】
【数１７０】
2007102477.tif 000172
、
【０２１６】
【数１７１】
2007102477.tif 000173
を与え，初期値
【０２１７】
【数１７２】
2007102477.tif 000174
、
【０２１８】
【数１７３】
2007102477.tif 000175
、
【０２１９】
【数１７４】
2007102477.tif 000176
、
【０２２０】
【数１７５】
2007102477.tif 000177
を与える。さらに、
【０２２１】
【数１７６】
2007102477.tif 000178
を求め、スペクトル因子分解法により、
【０２２２】
【数１７７】
2007102477.tif 000179
、
【０２２３】
【数１７８】
2007102477.tif 000180
を求める。
ステップ２．
【０２２４】
【数１７９】
2007102477.tif 000181
とし、
【０２２５】
【数１８０】
2007102477.tif 000182
、
【０２２６】
【数１８１】
2007102477.tif 000183
を用いて一段予測誤差
【０２２７】
【数１８２】
2007102477.tif 000184
を
【０２２８】
【数１８３】
2007102477.tif 000185
により求める。
Ｓｔｅｐ３．逐次予測誤差法
【０２２９】
【数１８４】
2007102477.tif 000186
【０２３０】
【数１８５】
2007102477.tif 000187
【０２３１】
【数１８６】
2007102477.tif 000188
を実行する。
Ｓｔｅｐ４．推定値
【０２３２】
【数１８７】
2007102477.tif 000189
から
【０２３３】
【数１８８】
2007102477.tif 000190
を算出し、スペクトル因子分解法により、
【０２３４】
【数１８９】
2007102477.tif 000191
、
【０２３５】
【数１９０】
2007102477.tif 000192
を求める。
Ｓｔｅｐ５．新しく求まった
【０２３６】
【数１９１】
2007102477.tif 000193
、
【０２３７】
【数１９２】
2007102477.tif 000194
を用いて、再び一段予測誤差
【０２３８】
【数１９３】
2007102477.tif 000195
を
【０２３９】
【数１９４】
2007102477.tif 000196
により求める。
Ｓｔｅｐ６．予測誤差分散、システム雑音、観測雑音を
【０２４０】
【数１９５】
2007102477.tif 000197
【０２４１】
【数１９６】
2007102477.tif 000198
により求める。
【０２４２】
【数１９７】
2007102477.tif 000199
【０２４３】
【数１９８】
2007102477.tif 000200
を計算し
【０２４４】
【数１９９】
2007102477.tif 000201
を求め、
【０２４５】
【数２００】
2007102477.tif 000202
により
【０２４６】
【数２０１】
2007102477.tif 000203
を求めてStep２へ戻る。
Ｓｔｅｐ７．
【０２４７】
【数２０２】
2007102477.tif 000204
まで到達したときに収束していなければ、最終推定値
【０２４８】
【数２０３】
2007102477.tif 000205
を初期値
【０２４９】
【数２０４】
2007102477.tif 000206
に代入し、再びStep １から繰り返す。
【０２５０】
次に、状態空間表現からカルマン・固定区間スムーザを導出する。（３）式にはＬだけ遅れ項が存在するので、
【０２５１】
【数２０５】
2007102477.tif 000207
を
【０２５２】
【数２０６】
2007102477.tif 000208
のように定義する。ただし、
【０２５３】
【数２０７】
2007102477.tif 000209
を定義した。いま、
【０２５４】
【数２０８】
2007102477.tif 000210
と２Ｌの最大値を
【０２５５】
【数２０９】
2007102477.tif 000211
として定義し、状態変数
【０２５６】
【数２１０】
2007102477.tif 000212
をＲ次に拡大化した新たな状態変数
【０２５７】
【数２１１】
2007102477.tif 000213
を
【０２５８】
【数２１２】
2007102477.tif 000214
で定義する。このとき状態空間モデルを記述し直すと、
【０２５９】
【数２１３】
2007102477.tif 000215
に対して
【０２６０】
【数２１４】
2007102477.tif 000216
（式１５）
【０２６１】
【数２１５】
2007102477.tif 000217
（式１６）
を得ることができる。ここで、Ｇは、
【０２６２】
【数２１６】
2007102477.tif 000218
（式１７）
であり、
【０２６３】
【数２１７】
2007102477.tif 000219
であれば、
【０２６４】
【数２１８】
2007102477.tif 000220
（式１８）
【０２６５】
【数２１９】
2007102477.tif 000221
（式１９）
で定義される。もし、
【０２６６】
【数２２０】
2007102477.tif 000222
であれば、
【０２６７】
【数２２１】
2007102477.tif 000223
の第1行目が
【０２６８】
【数２２２】
2007102477.tif 000224
となる（２Ｌ＋１）×（２Ｌ＋１）サイズの行列として定義され、
【０２６９】
【数２２３】
2007102477.tif 000225
であれば
【０２７０】
【数２２４】
2007102477.tif 000226
が
【０２７１】
【数２２５】
2007102477.tif 000227
となる１×Ｑの行列として定義される。
１５式乃至１９式中の変数右肩にある行を示す添え字（ｉ）を省略し、条件付き期待値と誤差共分散行列を、
【０２７２】
【数２２６】
2007102477.tif 000228
【０２７３】
【数２２７】
2007102477.tif 000229
で定義すると、以下のカルマン・固定区間スムーザアルゴリズムを与えることができる。初期分布の平均と分散を
【０２７４】
【数２２８】
2007102477.tif 000230
【０２７５】
【数２２９】
2007102477.tif 000231
（通常は
【０２７６】
【数２３０】
2007102477.tif 000232
をデータの平均値、
【０２７７】
【数２３１】
2007102477.tif 000233
を定常分散値）として与え、
【０２７８】
【数２３２】
2007102477.tif 000234
に対してカルマンフィルタは
【０２７９】
【数２３３】
2007102477.tif 000235
（式２０）
【０２８０】
【数２３４】
2007102477.tif 000236
【０２８１】
【数２３５】
2007102477.tif 000237
【０２８２】
【数２３６】
2007102477.tif 000238
【０２８３】
【数２３７】
2007102477.tif 000239
の逐次的計算で実行される。また、カルマンフィルタで求まった
【０２８４】
【数２３８】
2007102477.tif 000240
、
【０２８５】
【数２３９】
2007102477.tif 000241
を初期値として、
【０２８６】
【数２４０】
2007102477.tif 000242
に対して固定区間スムーザは
【０２８７】
【数２４１】
2007102477.tif 000243
【０２８８】
【数２４２】
2007102477.tif 000244
【０２８９】
【数２４３】
2007102477.tif 000245
の逆逐次計算によって実行される。
【０２９０】
以下に，カルマン・固定区間スムーザのアルゴリズムをまとめる． 
Ｓｔｅｐ１．  劣化画像
【０２９１】
【数２４４】
2007102477.tif 000246
をＤＦＴし、
【０２９２】
【数２４５】
2007102477.tif 000247
を計算する。
Ｓｔｅｐ２．  （２０）式以下の カルマン・固定区間スムーザを実行する。
Ｓｔｅｐ３．  Ｓｔｅｐ２を適切な回数繰り返す。
Ｓｔｅｐ４．  状態推定値から
【０２９３】
【数２４６】
2007102477.tif 000248
を求める。
Ｓｔｅｐ５．  最終的に得られた
【０２９４】
【数２４７】
2007102477.tif 000249
を逆ＤＦＴすることにより、修復画像
【０２９５】
【数２４８】
2007102477.tif 000250
を得る。
【０２９６】
実施の形態２．
図８を参照し、本発明に最適な実施の形態２．について詳細に述べる。
【０２９７】
図８において、劣化画像は画像修復装置１０１における入力端子１１１から入力される。また劣化検出部５１では、角速度センサ出力を用いたカメラの変位量を下に、数学的画像モデルに含まれるインパルス応答関数を推定し、インパルス応答関数として出力する。画像修復装置１０１の入力端子１１２では、劣化検出部５１で推定されたインパルス応答関数を入力し、劣化画像とインパルス応答関数を境界値補正部１に入力する。また、原画像推定部から修復画像を境界値補正部１に入力し、境界値補正部１では、境界値補正劣化画像を計算し、原画像推定部３に出力する。原画像推定部３では画像修復アルゴリズムを用いて原画像を推定し、画像修復装置１０１の出力端子１１３から修復画像を出力する。また、境界値補正部１へも修復画像を出力する。画像修復装置１０１の入力端子１１２から入力したインパルス応答関数は、出力端子１１４からパラメータの最終推定値として出力される。
【０２９８】
上記のようにセンサ出力を用いることによって、さらに良好な修復画像の推定が可能である。
【図面の簡単な説明】
【０２９９】
【図１】本発明における撮像領域と手ぶれ等により、写り込む可能性のある領域を示した図である。
【図２】本発明における撮像領域と手ぶれ等により、写り込む可能性のある領域および上記写り込む可能性のある領域から写り込まれることで劣化を受ける領域を示す図である。
【図３】本発明における画像修復装置の境界値補正部において、０埋めされる領域を示す図である。
【図４】本発明における画像修復装置の境界値補正部において、フィルタリングされる領域を示す図である。
【図５】本発明における画像修復装置の境界値補正部において、劣化画像を埋める領域を示す図である。
【図６】本発明における画像修復装置の構成を示すブロック図である。
【図７】本発明における境界値補正部の構成を示すブロック図である。
【図８】本発明における画像修復のフロー図である。
【図９】本発明の実施の形態２における画像修復装置の構成を示すブロック図である。
【符号の説明】
【０３００】
１  境界値補正部
３  パラメータ推定部
４  原画像推定部
５  パラメータ制御部
【図１】
2007102477.tif 000251
【図２】
2007102477.tif 000252
【図３】
2007102477.tif 000253
【図４】
2007102477.tif 000254
【図５】
2007102477.tif 000255
【図６】
2007102477.tif 000256
【図７】
2007102477.tif 000257
【図８】
2007102477.tif 000258
【図９】
2007102477.tif 000259