JP6661491B2 - 画像処理装置および画像処理方法 - Google Patents

Description

本発明は、複数枚の画像からなる画像群（例えば焦点位置を変えて撮影した複数枚の画像あるいは動画像データ）中に含まれる固定パターンノイズを低減し、高画質化する技術に関する。

カメラや顕微鏡などの撮像光学系を通して、焦点位置を変えて被写体を撮影した複数枚の画像群（以降、焦点ぼけ画像群と呼ぶ）から、光線空間の復元を介して、任意視点画像や任意焦点ぼけ画像を生成する技術が知られている。例えば、特許文献１では、３次元のコンボリューションモデルに合致するように各画像の座標変換処理を行い、３次元周波数空間上でぼけを変える３次元フィルタリング処理を行うことで任意視点画像や任意焦点ぼけ画像を生成する手法を開示している。また特許文献１以外にも、数学的にほぼ等価な式を用いて、焦点ぼけ画像群の視線方向の積分画像に対するフィルタ処理により光線空間の復元を行い、任意視点画像や任意焦点ぼけ画像を生成する手法が提案されている。なお、本明細書では、焦点ぼけ画像群から光線空間の復元を介して任意視点画像や任意焦点ぼけ画像を生成する画像生成手法を総称して「フィルタ型方式」と呼ぶこととする。

カメラや顕微鏡などで撮影した焦点ぼけ画像群に対して、これらのフィルタ型方式を適用すれば、カメラや顕微鏡の光学系に変更を加えることなく、撮影後の後処理によって従来にない機能を提供できる。例えば、カメラの場合には撮影後にぼけ味を変更する新たな写真の楽しみ方を提供できる。顕微鏡の場合には視点を変えることで被写体の立体形状の直感的な把握を可能にできる。また、開口数の大きなレンズを用いた高倍観察時において、被写界深度を拡大することで焦点位置から外れた病変部の見落としを防ぐことができる。

特開２００７−１２８００９号公報 特開２０１４−０９０４０１号公報 特開２０１３−２０７７８８号公報

しかしながら上述した従来の技術においては、以下のような問題があった。焦点ぼけ画像群に対し、特許文献１に開示される手法に代表されるフィルタ型方式を適用すると、イメージセンサの固定パターンノイズ（以降、固定ノイズと呼ぶ）に起因して、生成される画像の画質が劣化するという問題が生じる場合がある。この問題は、焦点ぼけ画像群を構成する画像のあいだで像の拡縮（スケールの変化）が非常に小さい場合、つまり、固定ノイズが画像内の殆ど同じ位置に現れる画像群の場合に顕在化する。そのような焦点ぼけ画像群は、例えば両側テレセントリック光学系のように、焦点位置を変えても像とイメージセンサの相対的な位置・大きさが殆ど変化しない光学系で撮影した場合に得られる。

前記課題への対策として、特許文献２では撮像装置の制御によって固定ノイズを回避する手法を開示している。この手法は焦点位置を変えて被写体を撮影する際、焦点位置毎に像とイメージセンサの相対的な位置をずらして撮影することで、焦点ぼけ画像群の同一画素で固定ノイズが重畳されることを抑える。しかし、この手法は、すでに撮影済みの焦点ぼけ画像群から全焦点画像を求める場合の画質改善には適用できない。

ここまで焦点ぼけ画像群での問題を記載したが、同様の問題は、動画像データのように同じ光学系で撮影された画像群でも生じうる。特許文献３では動画像から固定ノイズを低減する手法が開示されている。しかし、この手法は比較的平坦な領域に対して平均値を減算した差分画像が固定ノイズであると推定する簡易的な処理であるため、固定ノイズの推定精度は低い。そのため、様々な被写体の像を含む画像データから高精度な固定ノイズ低減を実現することは原理的に困難である。また加法性の固定ノイズ以外のノイズ、例えば、イメージセンサのセルごとの受光部の感度ばらつきや増幅器のゲインばらつきなどに起因する乗法性の固定ノイズの低減を行うことはできない。なお、本明細書では、イメージセンサの画素と画像データの画素とを明確に区別するため、イメージセンサの画素を「セル」と表記する。

本発明は、このような問題点に鑑みてなされたものであり、同一のイメージセンサで取得された複数枚の画像に含まれる、イメージセンサに起因する固定ノイズを低減する新規な技術を提供することを目的とする。

本発明の第一態様は、コンピュータが、同一のイメージセンサで撮像することで取得された複数の入力画像のデータを取得する入力ステップと、コンピュータが、逐次近似法を用いた反復計算によって、前記複数の入力画像の画質を全体的に向上させる、各入力画像に対する輝度変更の最適解を求める最適化処理ステップと、を有し、前記最適化処理ステップでは、各入力画像の同じ位置の画素に対して共通の輝度変更を行うという制約の下で、各入力画像の画素ごとの輝度変更の最適解が求められることを特徴とする画像処理方法を提供する。

本発明の第二態様は、本発明に係る画像処理方法の各ステップをコンピュータに実行させることを特徴とするプログラムを提供する。

本発明の第三態様は、同一のイメージセンサで撮像することで取得された複数の入力画像のデータを取得する入力手段と、逐次近似法を用いた反復計算によって、前記複数の入力画像の画質を全体的に向上させる、各入力画像に対する輝度変更の最適解を求める最適化処理手段と、を有し、前記最適化処理手段は、各入力画像の同じ位置の画素に対して共通の輝度変更を行うという制約の下で、各入力画像の画素ごとの輝度変更の最適解を求めることを特徴とする画像処理装置を提供する。

本発明によれば、同一のイメージセンサで取得された複数枚の画像に含まれる、イメージセンサに起因する固定ノイズを低減することができる。

バーチャル・スライド・システムの構成を示す図。 本計測ユニット１０１の内部構成を示す図。 画像処理装置（ホストコンピュータ）１１０の内部構成を示す図。 固定ノイズ低減のアイデアを説明する図。 実施例１の固定ノイズ低減処理を示すフローチャート。 最適化処理ステップＳ５０２の内部処理を示すフローチャート。 解空間制約ステップＳ６０３の内部処理を示すフローチャート。 制約条件Ｃ１およびＣ２が表す意味を説明するための図。 制約条件Ｃ１適用ステップＳ７０２の内部処理を示すフローチャート。 制約条件Ｃ２適用ステップＳ７０３の内部処理を示すフローチャート。 実施例３の固定ノイズ低減処理を示すフローチャート。 実施例４の固定ノイズ低減処理を示すフローチャート。 固定ノイズの低減が可能な動画像データの一例。

本発明は、同じ光学系（撮像系）を用いて撮影された複数枚の画像からなる画像群（例えば、焦点ぼけ画像群、動画像データなど）に含まれる固定ノイズを、画像処理（後処理）によって低減する技術に関する。本発明の構成ないし方法は、デジタルカメラ、デジタルビデオカメラ、デジタル顕微鏡、バーチャル・スライド・システムを含む様々な撮像装置で得られた画像に対し適用可能である。中でも、固定ノイズの影響が顕著となる、両側テレセントリック光学系で撮影した焦点ぼけ画像群に対し、本発明は好ましく適用できる。
以下に詳しく述べるように本発明の画像処理によれば、撮像装置のイメージセンサに起因する固定ノイズを高精度に低減した画像を得ることができる。また、焦点ぼけ画像群における固定ノイズを低減することによって、焦点ぼけ画像群からフィルタ型方式で任意視点画像や任意焦点ぼけ画像を生成する際に問題となっていた画質の劣化を抑制することができる。さらに、本発明による画像処理系を撮像装置と組み合わせることで、撮像装置のイメージセンサに求められる性能要件を緩和でき、イメージセンサの設計の自由度向上や低コスト化が実現できるという利点もある。

＜実施例１＞
（システム全体構成）
実施例１では、一例として、図１のような構成のバーチャル・スライド・システムに本発明を適用した例について述べる。
バーチャル・スライド・システムは、検体の撮像データを取得する撮像装置（バーチャル・スライド・スキャナとも呼ぶ）１２０とそのデータ処理・制御を行う画像処理装置（ホストコンピュータとも呼ぶ）１１０およびその周辺機器から構成されている。

画像処理装置１１０には、ユーザからの入力を受け付けるキーボードやマウスなどの操作入力デバイス１１１、処理画像を表示するディスプレイ１１２が接続される。また画像処理装置１１０には、記憶装置１１３、他のコンピュータシステム１１４が接続されている。
多数の検体（プレパラート）の撮像をバッチ処理で行う場合、画像処理装置１１０の制御の下で撮像装置１２０が各検体を順に撮像し、画像処理装置１１０が各検体の画像データに対し必要な処理を施す。そして得られた各検体の画像データは、大容量のデータストレージである記憶装置１１３又は他のコンピュータシステム１１４へ伝送され、蓄積される。

撮像装置１２０での撮像（プレ計測および本計測）は、ユーザの入力を受けて画像処理装置１１０がコントローラ１０８へ指示を送り、次にコントローラ１０８が本計測ユニット１０１とプレ計測ユニット１０２を制御することで実現される。
本計測ユニット１０１は、プレパラート内の検体診断のための高精細な画像を取得する撮像ユニットである。プレ計測ユニット１０２は、本計測に先立って撮像を行う撮像ユニットで、本計測で精度の良い画像取得をするための撮像制御情報取得を目的とした画像取得を行う。

コントローラ１０８には変位計１０３が接続され、本計測ユニット１０１またはプレ計測ユニット１０２内のステージに設置されるプレパラートの位置や距離が測定できる構成となっている。変位計１０３は、本計測およびプレ計測を行うにあたり、プレパラート内の検体の厚みを計測するために用いられる。
またコントローラ１０８には、本計測ユニット１０１およびプレ計測ユニット１０２の撮像条件を制御するための開口絞り制御１０４、ステージ制御１０５、照明制御１０６、センサ制御１０７が接続されている。そして、それぞれはコントローラ１０８からの制御信号に従って、開口絞り、ステージ、照明、イメージセンサの動作を制御する構成となっている。

ステージはプレパラートを保持・移動・位置決めする機構である。ステージには、プレパラートを光軸に垂直な方向に移動するＸＹステージと、光軸に沿った方向に移動するＺステージがある。ＸＹステージは検体内の撮像エリアを光軸に垂直な方向（ｘ，ｙ方向）に移動するために用いられ、Ｚステージは検体内の奥行き方向（ｚ方向）に焦点位置を変えるために用いられる。図示しないが、撮像装置１２０には、複数のプレパラートをセットすることが可能なラックと、ラックからステージ上の撮像位置へとプレパラートを送り出す搬送機構とが設けられている。バッチ処理の場合は、コントローラ１０８が搬送機構を制御することで、ラックからプレパラートを１枚ずつ、プレ計測ユニット１０２のステージ、本計測ユニット１０１のステージの順に送り出す。

本計測ユニット１０１およびプレ計測ユニット１０２には、撮像した画像を用いてオートフォーカスを実現するＡＦユニット１０９が接続されている。ＡＦユニット１０９は、コントローラ１０８を介して、本計測ユニット１０１、プレ計測ユニット１０２のステージの位置を制御することで合焦位置を探し出すことが出来る。オートフォーカスの方式は画像を用いるパッシブ型であり、公知の位相差検出方式やコントラスト検出方式が用いられる。

（本計測ユニットの構成）
図２は、実施例１における本計測ユニット１０１の内部構成を示す図である。
光源２０１の光は、照明光学系２０２を通じ、光量ムラが無いように均一化され、ステージ２０３の上に設置されたプレパラート２０４を照射する。プレパラート２０４は観察対象となる組織の切片や塗抹した細胞をスライドグラス上に貼り付け、封入剤とともにカバーグラスの下に固定したものであり、検体（被写体）を観察可能な状態に準備したものである。

結像光学系２０５は、検体の像を拡大して撮像手段である撮像部２０７に導くものである。プレパラート２０４を通った光は、結像光学系２０５を介して、撮像部２０７上の撮像面で結像する。結像光学系２０５の中には開口絞り２０６が存在し、開口絞り２０６を調整することで被写界深度が制御できる。

撮像にあたっては、光源２０１を点灯させ、プレパラート２０４に光を照射する。そして、照明光学系２０２、プレパラート２０４、結像光学系２０５を通って撮像面に結像した像を撮像部２０７のイメージセンサで受光する。モノクロ（グレースケール）撮影時には光源２０１を白色で露光し、１回撮像を行う。カラー撮影時には、ＲＧＢの３つの光源２０１で順番に露光し、３回撮像を行うことで、カラー画像を取得する。

撮像面で結像した検体の像は、撮像部２０７で光電変換され、Ａ／Ｄ変換を実行した後、電気信号として画像処理装置１１０に送られる。撮像部２０７は、複数のイメージセンサから構成されることを想定しているが、単一のセンサで構成されていても良い。また、本実施例では、Ａ／Ｄ変換を実行した後、イメージセンサの固定ノイズ低減処理を含むノイズ低減や色変換処理、鮮鋭化処理に代表される現像処理を画像処理装置１１０内部で行うとする。しかし、現像処理は撮像部２０７に接続された専用の画像処理ユニット（不図示）で行い、その後、画像処理装置１１０にデータを送信することも可能であり、そのような形態での実施も本発明の範疇となる。

一回の撮影で検体全体の画像を取得できない場合には、ステージ２０３（ＸＹステージ）をｘ方向および／またはｙ方向に移動しながら複数回の分割撮影を行い、得られた複数の分割画像を合成（繋ぎ合わせ）して検体全体の画像を生成する。また、ステージ２０３（Ｚステージ）をｚ方向に移動しつつ複数回の撮影を行うことで、光軸方向（深さ方向）の焦点位置が異なる複数枚の画像を取得する。本明細書では、撮像光学系の光軸方向に焦点位置を異ならせて被写体を撮像することで得た複数枚の画像からなる画像群を「焦点ぼけ画像群」と呼ぶ。なお、焦点ぼけ画像群を「Ｚスタック画像」、焦点ぼけ画像群を構成する各々の画像を「レイヤー画像」と呼ぶこともある。

（画像処理装置の構成）
図３は、本実施例における画像処理装置（ホストコンピュータ）１１０の内部構成を示す図である。
ＣＰＵ（プロセッサ）３０１は、ＲＡＭ３０２やＲＯＭ３０３に格納されているプログラムやデータを用いて画像処理装置全体の制御を行う。またＣＰＵ３０１は、各種演算処理やデータ処理、例えば固定ノイズ低減処理、現像・補正処理、合成処理、圧縮処理、任意視点・焦点ぼけ画像生成処理等を行う。

ＲＡＭ３０２は、記憶装置１１３からロードされたプログラムやデータ、他のコンピュータシステム１１４からネットワークＩ／Ｆ（インターフェース）３０４を介してダウンロードしたプログラムやデータを一時的に記憶するメモリである。ＲＡＭ３０２は、ＣＰＵ３０１が各種の処理を行うために必要とするワークエリアも備える。ＲＯＭ３０３は、コンピュータの機能プログラムや設定データなどを記憶するメモリである。ディスプレイ制御装置３０６は、画像や文字等をディスプレイ１１２に表示させるための制御処理を行う。ディスプレイ１１２は、ユーザに入力を求めるための画面表示を行うとともに、撮像装置１２０から取得しＣＰＵ３０１で処理した画像データを画像表示する。

操作入力デバイス１１１は、キーボードやマウスなどＣＰＵ３０１に各種の指示を入力することのできるデバイスにより構成される。ユーザは、撮像装置１２０の動作を制御する情報を操作入力デバイス１１１により入力する。３０８は、操作入力デバイス１１１を介して入力された各種の指示等をＣＰＵ３０１に通知するためのＩ／Ｏである。

記憶装置１１３は、ハードディスクなどの大容量情報記憶装置である。これは、ＯＳ（オペレーティングシステム）や後述する処理をＣＰＵ３０１に実行させるためのプログラム、またバッチ処理によりスキャンした画像データ、処理後の画像データなどを記憶する。
記憶装置１１３への情報の書き込みや記憶装置１１３からの情報の読み出しは、Ｉ／Ｏ３１０を介して行われる。制御Ｉ／Ｆ３１２は、撮像装置１２０を制御するためのコントローラ１０８と制御コマンド（信号）をやりとりするためのＩ／Ｆである。

コントローラ１０８は、本計測ユニット１０１およびプレ計測ユニット１０２を制御する機能を持つ。画像Ｉ／Ｆ（インターフェース）３１３には、上述以外のインターフェース、例えばＣＭＯＳイメージセンサやＣＣＤイメージセンサの出力データを取り込むための外部インターフェースが接続されている。なお、インターフェースとしてはＵＳＢ、ＩＥＥＥ１３９４などのシリアルインターフェースやカメラリンクなどのインターフェースを使うことが出来る。この画像Ｉ／Ｆ３１３を通じて本計測ユニット１０１やプレ計測ユニット１０２が接続される。

（任意視点・焦点ぼけ画像生成プログラム）
画像処理装置１１０には、特許文献１に開示される手法を一例とするフィルタ型方式に
よる画像生成処理をコンピュータに実行させるためのプログラム（任意視点・焦点ぼけ画像生成プログラムと呼ぶ）が実装されている。画像処理装置１１０は、撮像装置１２０や記憶装置１１３、他のコンピュータシステム１１４を介して取得した焦点ぼけ画像群から、任意視点画像や任意焦点ぼけ画像や全焦点画像などを生成することができる。

（フィルタ型方式で生成した全焦点画像で固定ノイズが目立つ理由の説明）
顕微鏡の撮像系（図１の本計測ユニット１０１）では一般的に両側テレセントリック光学系が用いられ、焦点ぼけ画像群を構成する各画像のあいだでは検体像の拡縮（スケールの変化）が非常に小さくなる。このような焦点ぼけ画像群に対し特許文献１のフィルタ型方式を適用する場合、３次元のコンボリューションモデルに合致するように各画像の座標変換処理を行う必要はない。そのため、全焦点画像ａ（ｘ，ｙ）は、焦点ぼけ画像群の光軸方向の積分ｂ（ｘ，ｙ）を、３次元ぼけの光軸方向の積分ｃ（ｘ，ｙ）でデコンボリューションすることで得られる。デコンボリューションは周波数空間での割算となるため、全焦点画像ａ（ｘ，ｙ）は以下の式により求まる。
ただし、Ｂ（ｕ，ｖ）＝Ｆ（ｂ（ｘ，ｙ））、Ｃ（ｕ，ｖ）＝Ｆ（ｃ（ｘ，ｙ））であり、Ｆはフーリエ変換、Ｆ^−１は逆フーリエ変換を表す。ｕ，ｖはそれぞれｘ，ｙに対応する空間周波数である。

式（１）は、Ｃ（ｕ，ｖ）^−１が、焦点ぼけ画像群の光軸方向の積分ｂ（ｘ，ｙ）に対する周波数フィルタとなり、全焦点画像を生成することを表している。３次元ぼけとして、ピント位置からの距離に従ってぼけの半径が大きくなるガウスぼけを考える。この場合、周波数フィルタＣ（ｕ，ｖ）^−１の値は、直流成分では焦点ぼけ画像群の枚数分の１となるが、周波数が高くになるにつれて上昇し、最大周波数では１に近づく特性を持っている。つまり、周波数フィルタＣ（ｕ，ｖ）^−１は低周波のノイズを抑制する作用をもつが、高周波のノイズを抑制する作用はもたない。

続いて、焦点ぼけ画像群の各画素に時間不変の固定ノイズとして標準偏差σ＝０．５の加法性のガウスノイズが存在する場合を考える。焦点ぼけ画像群の画像枚数が６４枚のとき、焦点ぼけ画像群の光軸方向の積分ｂ（ｘ，ｙ）では固定ノイズは重畳されて固定ノイズの標準偏差は３２となる。その後、周波数フィルタＣ（ｕ，ｖ）^−１によって、加法性のガウスノイズの低周波成分の振幅は２５５階調中の０．５に近い微小な値に戻る。しかし、高周波成分の振幅は依然、２５５階調中の３２に近い大きな値のままであり、全焦点画像中では大きなノイズとなって視認されることとなる。

以上の理由により、顕微鏡で得られた焦点ぼけ画像群からフィルタ型方式で全焦点画像を求める場合に、固定ノイズに起因する全焦点画像の品質の劣化が顕在化する。固定ノイズの影響は、全焦点画像の生成に用いる画像の枚数が増えるほど、大きくなる。

（任意視点・焦点ぼけ画像の高画質化）
高画質な任意視点画像、任意焦点ぼけ画像を生成するには、任意視点・焦点ぼけ画像生成プログラムに入力する焦点ぼけ画像群に対する前処理として、以降で述べる固定ノイズ低減処理を行い、焦点ぼけ画像群から予め固定ノイズを低減するとよい。なお、前述のように固定ノイズは微弱な場合もあるため、固定ノイズ低減後の焦点ぼけ画像群の画像の画素値は整数値に量子化せず、実数データとして任意視点・焦点ぼけ画像生成プログラムに入力することが望ましい。また、固定ノイズ低減処理は画像処理装置１１０で実現するため、コンピュータで実行可能な形で実装したプログラムを用いて行うとする。

（本実施例で用いる焦点ぼけ画像群に関して）
本実施例では焦点ぼけ画像群は、縦Ｎ画素、横Ｎ画素の、Ｍ枚の画像から構成されるとする。（本実施例では説明簡略化のため、縦および横の画素数が等しい例で説明をするが、縦と横の画素数が一致しない画像にも本発明の方法を適用可能である。）
撮像時にイメージセンサの固定ノイズが加わった焦点ぼけ画像群を
で表す。Ｒ^Ｎ×Ｎは、実数の要素からなるＮ行Ｎ列の行列を表す。添え字のｋは焦点ぼけ画像群を構成する画像の番号であり、ｋは焦点位置と対応する。焦点ぼけ画像群を構成する画像ｙ_ｋは、イメージセンサの出力信号（ビニング処理により複数のセルの値を加算した信号も含む）そのものでもよい。あるいは、画像ｙ_ｋは、イメージセンサの出力信号に対し、イメージセンサのセルとの対応関係を維持した階調変換や補間等の現像処理をした画像でもよい。以降、画像の画素の値を輝度と呼ぶ。

焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）に対して、固定ノイズの低減のための輝度変更をした焦点ぼけ画像群を
で表す。また、撮像系の光学的なぼけによる劣化のみを含む焦点ぼけ画像群を
で表す。
焦点ぼけ画像群を構成するｋ番目の画像ｙ_ｋおよびｘ_ｋの縦ｊ、横ｉの位置にある画素の輝度はそれぞれ、（ｙ_ｋ）_ｉ，ｊ、（ｘ_ｋ）_ｉ，ｊで表す。

本実施例では、焦点ぼけ画像群ｙ_１，…，ｙ_Ｍには、被写体の像の画像成分や時間変化には依存しない加法性の固定ノイズｎ∈Ｒ^Ｎ×Ｎが加わっているとする。すなわち、焦点ぼけ画像群ｙ_ｋは以下の式で表現できる。
固定ノイズｎの各画素（ｉ，ｊ）での値（ｎ）_ｉ，ｊはイメージセンサの様々な製造誤差に起因するが、本実施例では（ｎ）_ｉ，ｊは出現頻度分布が平均０、標準偏差σの正規分布に従うものとして扱う。ただし、固定ノイズｎの各要素の値の出現頻度分布が正規分布に従わない場合でも本実施例の方法では固定ノイズの低減が可能である。

（画像の滑らかさを表す指標について）
本実施例では、画像に含まれる固定ノイズが低減したかどうかを評価する指標として、画像の滑らかさを表す指標を用いる。画像の滑らかさは数値で定量的に表現出来ればよい。例えば、画像中の滑らかではない情報の総和、すなわち、エッジやテクスチャやノイズ等の強さを画像全体で積分した指標を用いてもよい。この場合、指標の値が小さい程、画像が滑らかであることを示す。画像の滑らかさを表す１つの指標として全変動（Ｔｏｔａｌ ＶａｒｉａｔｉｏｎまたはＴＶ）が知られている。全変動は画像の各画素の輝度の微分の絶対値を積分したものであり、以下の式で表される。
ここで、Ωは積分対象領域、ｕはｘ，ｙに関する関数、ｄｘ，ｄｙは微小区間である。

上記で定義される全変動は画像の変化が少ない平坦な領域ほど小さな値をとり、テクスチャやノイズが存在する画像の変化が多い領域では大きな値をとる。すなわち、全変動が小さいほど滑らかであることを示す。そのため、全変動は、画像からノイズを低減し、滑らかにする最適化問題において全変動最小化（Ｔｏｔａｌ Ｖａｒｉａｔｉｏｎ Ｍｉｎｉｍｉｚａｔｉｏｎ）の形で使われている。画像が離散的なデジタルデータの場合、輝度の微分は隣接画素間の輝度の差分で表されるため、全変動は以下のように隣接画素間の輝度の差分の絶対値の総和で定義できる。

（固定ノイズ低減のアイデア）
図４Ａ、図４Ｂは、本実施例における固定ノイズ低減のアイデアを説明する図である。４０１および４０４は、焦点位置を変えて撮影した４枚の画像からなる焦点ぼけ画像群の一例である。左から順に焦点位置ｚ＝０，１，２，３の画像である。各画像の網掛けは、各画像に同一の固定ノイズが加わっている様子を示している。焦点ぼけ画像群４０１および４０４は、少数の鮮鋭度の高い合焦位置の画像（ｚ＝０の画像）と比較的多数の滑らかな焦点がぼけた画像（ｚ＝１，２，３の３枚の画像）から構成されている。（図４Ａでは表現されていないが、ｚ＝０で合焦している被写体は合焦位置から離れるに従い、徐々にぼけていく。）

焦点ぼけ画像群の各画像の同じ位置の画素（対応画素と呼ぶ）は、イメージセンサの同じセルで取得されたデータに基づく輝度をもつため、共通の固定ノイズ（イメージセンサの同じセルに起因するノイズ）を含んでいる。したがって、各画像の対応画素に対して共通の輝度変更を行うという制約の下で、固定ノイズを低減することを考える。ただし、画像内のどの位置（画素）にどの程度の固定ノイズが加わっているかは不明のため、上記制約の下で、焦点ぼけ画像群の画質が全体的に向上するように、画像内の位置ごと（画素ごと）の輝度変更の最適解を求めるとよい。このとき、焦点ぼけ画像群の画質が全体的に向上したかどうかを評価する指標として、焦点ぼけ画像群全体の滑らかさを表す指標に注目する。

図４Ａは、焦点ぼけ画像群の各画像から固定ノイズでない輝度成分を減じる例である。４０２は、固定ノイズでない輝度成分、すなわち被写体の画像成分（例えばテクスチャ）を示す。４０３は焦点ぼけ画像群４０１を構成する各画像から画像成分４０２を減算して得られた焦点ぼけ画像群を示す。焦点ぼけ画像群４０３のｚ＝０の画像に注目すると、被写体の画像成分が無くなった分、画像の滑らかさは向上している。しかし、焦点ぼけによって滑らかだった焦点ぼけ画像群４０１のｚ＝１〜３の画像に関しては、焦点ぼけ画像群４０３のｚ＝１〜３に示すように、逆にテクスチャが発生し、画像の滑らかさは低下する。従って、焦点ぼけ画像群全体ではｚ＝０での少数枚での画像の滑らかさの向上よりも、ｚ＝１〜３に示す多数枚での画像の滑らかさの低下が支配的となる。よって、焦点ぼけ画
像群全体で評価した場合、滑らかさは低下する。

図４Ｂは、焦点ぼけ画像群の各画像から固定ノイズに相当する輝度成分を減じる例である。４０５は焦点ぼけ画像群４０４に含まれる固定ノイズを表す。４０６は焦点ぼけ画像群４０４を構成する全ての画像から画像成分４０５を減算して得られた焦点ぼけ画像群を示す。図４Ｂに示すように、ｚ＝０〜３の全ての画像において固定ノイズの打ち消しにより滑らかさが向上するため、焦点ぼけ画像群４０６全体の滑らかさが向上する。仮に、固定ノイズの強度が僅かであっても、焦点ぼけ画像群全体では滑らかさが向上する。

前述したように、焦点ぼけ画像群全体の滑らかさは、一例として、焦点ぼけ画像群全体の全変動の小ささで評価することができる。従って、焦点ぼけ画像群を構成する各画像の対応画素に対して共通の輝度変更を行う制約の下、焦点ぼけ画像群全体の滑らかさを最大化する最適化、すなわち、焦点ぼけ画像群全体の全変動最小化を行うとよい。これにより、焦点ぼけ画像群に含まれる被写体の像（テクスチャなど）を維持したまま、固定ノイズを低減することができる。そのため、本実施例では、上記の処理を最適化問題として定式化することを考える。
なお、「共通の輝度変更」とは、各画像の対応画素の輝度あるいは前記輝度を階調変換した値に対して同じ値を加算、減算、乗算、又は、除算することを意味する。加法性の固定ノイズの場合は、各画像の対応画素の輝度に対し同じ値（補正量）を加算又は減算すればよいし、後の実施例で述べる乗法性の固定ノイズの場合は、各画像の対応画素の輝度に対し同じ値（補正係数）を乗算又は除算すればよい。また、同じく後述する加法性・乗法性以外の固定ノイズの場合でも、各画像の対応画素の輝度を階調変換した値に対して同じ値（補正係数）を乗算したのち、減算することで補正できる。

（数式・記号に関して）
最適化問題として定式化する前に、本実施例で用いる数式・記号について説明する。
ベクトルｘがｘ＝（ｘ_１，ｘ_２，…，ｘ_Ｎ）で表されるとき、ベクトルｘのＬ１ノルムとＬ２ノルムは以下の式で定義される。

また本実施例では画像データ（行列）ｕ∈Ｒ^Ｎ×Ｎに対するＬ１ノルムとＬ２ノルムを以下のように定める。

また画像データｕ∈Ｒ^Ｎ×Ｎに対する全変動は、式（５）に示すＴＶノルムで表現する。
なお、式（５）の（∇ｕ）_ｉ，ｊは横（ｘ）方向，縦（ｙ）方向の偏微分からなるベクトルで、
であるとする。

上記の画像データに対するノルムは、画像データｕ∈Ｒ^Ｎ×Ｎの各行ベクトルを転置して列方向に並べるベクトル表現を行って求めたベクトル
に対するノルムと等しい。なお、ｖｅｃ［・］は行列を並べ替えてベクトル表現する操作を表すとする。

続いて、最適化問題の定式化で用いる式について説明する。式（６）は最適化問題の定式化の一例である。
上式のａｒｇｍｉｎは関数ｆ（ｘ）が最小となる変数ｘの値を求めることを意味する。最小化する対象の関数ｆ（ｘ）を目的関数と呼ぶ。「ｓ．ｔ．」は「ｓｕｂｊｅｃｔ ｔｏ」の略で、その後ろに解ｘが満たすべき制約（式（６）では０≦ｘ≦２５５）を記載する。ｓ．ｔ．の後ろに記載する条件を制約条件と呼ぶ。式（６）の左辺のｘ^（＊）は前述の制約条件を満たす範囲でｘの値を変えて求めた、目的関数ｆ（ｘ）が最小となるときのｘの値を表す。

本実施例では、上記の最適化問題を解く際に、計算の便宜のため、焦点ぼけ画像群の複数の画像を統合する操作を行って１枚の統合画像を生成し、この統合画像に対し最適化の処理を適用する。以降、前記の操作で生成した１枚の画像を展開画像と呼ぶ。また、最適化処理の後は、固定ノイズが低減された展開画像から元の複数枚の画像から構成される焦
点ぼけ画像群を求める操作も行う。そのため、焦点ぼけ画像群から展開画像を生成する操作と、展開画像から焦点ぼけ画像群を生成する操作に関する説明を行う。

操作Ｔは、焦点ぼけ画像群を構成する複数枚の画像を統合し、１枚の画像を生成する操作である。Ｍ枚の画像からなる焦点ぼけ画像群ｙ_１，…，ｙ_Ｍを入力とし、展開画像ｇを生成する操作Ｔを以下のように表す。

操作Ｔには様々な操作が存在しうるが、一例として、以下のように引数として設定される複数枚の画像を順に列方向に連結する操作が挙げられる。
（ただし行列の右上のＴは行列の転置を表す）
その場合、元の画像ｙ_１，…，ｙ_Ｍの要素と出力である展開画像ｇの要素のあいだには以下の関係が成立する。

続いて、操作Ｔと逆の操作である操作Ｅについて説明する。操作Ｅは展開画像ｇから焦点ぼけ画像群を構成するｋ番目の画像ｙ_ｋを求める操作であり、以下のように表す。

焦点ぼけ画像群を構成する各画像ｙ_１，…，ｙ_Ｍにおける対応画素（ｉ，ｊ）の輝度をｋ＝１からＭまで順に並べたＭ次元のベクトルを、画素（ｉ，ｊ）の画素ベクトルと定義する。例えば、焦点ぼけ画像群ｙ_１，…，ｙ_Ｍの画素（ｉ，ｊ）の画素ベクトルは
で表される。

操作Ｔを用いて焦点ぼけ画像群ｙ_１，…，ｙ_Ｍから生成した展開画像ｇから、画素ベクトルを生成する操作を操作Ｖで表すと、操作Ｖは以下のように表される。
なお、式（１１）の左辺ｇ_{［ｉ，ｊ］}は画素ベクトルを表し、以下となる。

操作Ｔの変形例として、焦点ぼけ画像群を構成する各画像を行方向に連結する方法や行、列それぞれの方向に決められた枚数分連結する方法がある。また各画像を連結する際の境界部分に緩衝領域として輝度が全て０の領域や輝度変化を少なくする領域を入れても良い。その場合には操作Ｔに対応して操作Ｅ，操作Ｖの操作も変化する。本実施例では操作
Ｔとして式（８）に示す連結方法を行う場合について説明する。

また、輝度変更後の画像は元画像と同じ輝度の範囲（例えば０〜２５５）に収まっている必要がある。輝度０〜２５５の範囲を表す閉凸集合を以下のとおりＣ_２５５で表す。

（最適化問題の定式化）
図４Ｂで説明した焦点ぼけ画像群から固定ノイズを低減する問題は、以下の最適化問題として定式化できる。
ただし、ｘ^（＊）はｋ＝１，…，Ｍの全てに対して求める。

次に式（１２）で定式化した最適化問題の目的関数、制約条件について詳しく説明する。

（目的関数）
式（１２）の目的関数は、焦点ぼけ画像群を構成するＭ枚の画像から構成した１枚の展開画像に対するＴＶノルムを表す。一方、全焦点画像における固定ノイズの影響を低減するには、図４Ｂに示すように、焦点ぼけ画像群を構成するＭ枚の画像それぞれのＴＶノルムを低減する必要がある。しかし、展開画像に対するＴＶノルムと、焦点ぼけ画像群を構成する各画像のＴＶノルムの和とは、画像を並べる際に発生する境界での僅かな差を除けば等しい。したがって、式（１２）の目的関数により展開画像に対するＴＶノルムを最小化することと、元の焦点ぼけ画像群を構成するＭ枚の画像のＴＶノルムを全体的に最小化することとは、実質的に等価である。

（制約条件）
式（１２）の３つの制約条件のうち、１番目の制約条件は、輝度変更後の焦点ぼけ画像群を構成する各画像の輝度が０〜２５５の範囲に収まることを表す。画像回復処理での一般的な制約条件である。

２番目の制約条件は、いずれの画像の番号ｋにおいても、輝度変更前の画像ｙ_ｋと輝度変更後の画像ｘ_ｋの差が等しいことを表す。図４Ｂで説明したように、イメージセンサに起因する加法性の固定ノイズの場合、画像ｙ_ｋに含まれるノイズ成分は焦点位置（ｋ）によらず等しい。したがって、２番目の制約条件を与えることで、全ての画像ｙ_１，…，ｙ_ｋに等しく重畳された加法性の固定ノイズ成分の除去が可能となる。

２番目の制約条件を画素ごとの制約で表すことを考える。以下のように方向ベクトルｄとして要素が全て１のＭ次元のベクトルを設定する。
次に、撮影時の焦点ぼけ画像群ｙ_１，…，ｙ_Ｍと輝度変更後の焦点ぼけ画像群ｘ_１，…，ｘ_Ｍからそれぞれの画素（ｉ，ｊ）の画素ベクトルを求め、その差分をｕ_{［ｉ，ｊ］}で表すと
となる。ｕ_{［ｉ，ｊ］}は元の画素ベクトルからの補正量を示す画素ベクトルとなるため、以降、画素補正ベクトルと呼ぶ。

画素補正ベクトルｕ_{［ｉ，ｊ］}の各要素は２番目の制約条件により常に等しい値となる。よって画素補正ベクトルは、
で表される（ただし、ｃ_ｉ，ｊは画素ごとに異なる実数）。すなわち、画素補正ベクトルｕ_{［ｉ，ｊ］}と方向ベクトルｄは平行である。また完全に固定ノイズを取り除ける場合、ｃ_ｉ，ｊは各画素に加わった固定ノイズ（ｎ）_ｉ，ｊに−１をかけた値と等しくなる。

式（１２）の３番目の制約条件は、元の焦点ぼけ画像群から生成した展開画像Ｔ［ｙ_１，…，ｙ_Ｎ］と輝度変更後の展開画像Ｔ［ｘ_１，…，ｘ_Ｎ］の差のＬ２ノルムがε以下であることを示す。すなわち、εは輝度変更量の上限値を定める定数であり、３番目の制約条件は、展開画像Ｔ［ｙ_１，…，ｙ_Ｎ］と展開画像Ｔ［ｘ_１，…，ｘ_Ｎ］の間の輝度変更量（総量）が所定値（ε）をこえないように抑える制約を表している。この条件は、焦点ぼけ画像群を構成する画像ごとに、元の画像と輝度変更後の画像の差のＬ２ノルムを所定値（ε_１）以下とする制約条件
よりも緩和された条件である。しかし、２番目の制約条件（画素ベクトルによる制約）があるため、特定の画像ｘ_ｋで大きく偏った輝度変更を許すことはない。

なお、固定ノイズｎ∈Ｒ^Ｎ×Ｎが、各輝度の出現頻度分布が平均０、標準偏差σの正規分布に従うガウスノイズの場合、前記固定ノイズｎのＬ２ノルムは
となる。従って、固定ノイズの標準偏差σが既知の場合、式（１２）の３番目の制約条件の場合では
のようにεを設定すれば過剰な輝度変更を抑えた適切なノイズ低減が実現できる。なお、式（１４）のε_１はε_１＝Ｎσである。

元の焦点ぼけ画像群ｙ_１，…，ｙ_Ｍから操作Ｔで生成した展開画像をｇ、輝度変更後の焦点ぼけ画像群ｘ_１，…，ｘ_Ｍから操作Ｔで生成した展開画像をｆとし、式（１２）を整理すると、以下の式（１５）、（１６）が得られる。

すなわち、展開画像ｇから固定ノイズを除いた展開画像ｆ^（＊）を求める最適化問題として定式化できる。ただし、閉凸集合Ｃ_{２５５，Ｍ}は、以下に示すように、輝度変更後の展開画像ｆの要素が［０，２５５］の範囲外には解を持たない制約とする。

式（１５）の２番目の制約条件にある画素ベクトルｆ_{［ｉ，ｊ］}、ｇ_{［ｉ，ｊ］}は、それぞれ以下の式
で求めることができる。式（１５）の２番目の制約条件は、任意の画素（ｉ，ｊ）の画素補正ベクトル（ｆ_{［ｉ，ｊ］}−ｇ_{［ｉ，ｊ］}）が方向ベクトルｄと平行である、という条件を表す。「‖」は平行を表す記号である。

最適化アルゴリズムを適用しやすくするため、式（１５）を式（１７）のように変形する。
式（１７）のＩ（ｆ）は閉凸集合Ｃ_{２５５，Ｍ}の指示関数（ｉｎｄｉｃａｔｏｒ ｆｕｎｃｔｉｏｎ）であり、ｆの要素が［０，２５５］の範囲内のときには０、［０，２５５］の範囲外のときには∞（無限大）となる。また、以降、式（１７）の１番目の制約条件を制約条件Ｃ１、２番目の制約条件を制約条件Ｃ２と呼ぶ。

詳細は後述するが、式（１５），（１７）で定式化した最適化問題は凸最適化問題となる。凸最適化問題とは目的関数が凸関数で制約条件が閉凸集合である最適化問題で、最適解は全体最適解となることが保証され、どのような初期値を与えても全体最適解に至る取扱いの容易な最適化問題となる。

最適化問題を解くアルゴリズムの１つに逐次近似法（反復法）がある。逐次近似法とは、初期値を元に反復計算式から近似解を求め、その後、反復終了条件を満たすまで反復計算式を用いて近似解を求める計算を繰り返すことで近似解の精度を高める方法である。すなわち、式（１５），（１７）で定式化した凸最適化問題に逐次近似法を適用することで反復計算式が導きだせ、反復計算を繰り返すことで、精度の良い近似解を求めることがで
きる。反復計算式は適用する最適化アルゴリズムに依存する。

このとき、最適化アルゴリズムのパラメータを適切に設定すれば、反復回数ｎの増加とともに近似解は目的関数の最小（または最大）値に収束する。解への収束判定基準としては様々な基準を用いることができる。例えば、反復計算前後の近似解の差のＬ２ノルム
又は、相対評価
などを収束判定基準として用いることができる。ただし、ｆ^（ｎ）は反復回数ｎ回における近似解を表す。

次に、凸最適化問題の基本的な用語である凸関数、凸集合、距離射影（凸射影）について説明する。

（凸関数）
任意のｘ，ｙ∈Ｒ^Ｎとλ∈［０，１］に対して
が常に成り立つとき、ｆ（ｘ）はＲ^Ｎ上で凸関数であるという。言わば、ｆ（ｘ）上の任意の２点を結んだ線分が常にｆ（ｘ）の上側にあるような関数が凸関数である。なお、式（１７）の目的関数は凸関数であり、また、目的関数を構成する２つの式のいずれも凸関数である。

（凸集合）
ＣをＲ^Ｎの部分集合としたとき、全てのｘ，ｙ∈Ｃと全てのλ∈［０，１］に対して
となるとき、Ｃを凸集合という。言わば、ｘとｙとを結ぶ線分上の各点が必ずＣに属すことを表す。また凸集合Ｃが閉集合であるとき、Ｃは閉凸集合であるという。

式（１７）の制約条件Ｃ１を満たす集合Ｃ_１は、以下の理由から閉凸集合となる。すなわち、画素（ｉ，ｊ）の画素ベクトルがとりうる集合をＣ_{１（ｉ，ｊ）}で表すと、制約条件Ｃ１は方向ベクトルｄ上の直線で表せる。任意の２点を結ぶ線分上の各点は必ず線分に属すため、Ｃ_{１（ｉ，ｊ）}は閉凸集合となり、その直積集合であるＣ_１も閉凸集合となる。
また、制約条件Ｃ２を満たす集合Ｃ_２も閉凸集合である。

（距離射影）
凸最適化問題では、制約条件を満たす解を求めるために距離射影と呼ばれる写像を適用する。距離射影とは、空でない閉凸集合Ｃ⊂Ｒ^Ｎと任意のｘ∈Ｒ^Ｎに対して
を満たす唯一の点Ｐ_Ｃ（ｘ）∈Ｃが存在し、ｘにＰ_Ｃ（ｘ）を対応させる写像を表す。すなわち、距離射影とは、Ｃを閉凸集合としたとき、ｘ∈Ｒ^Ｎを、閉凸集合Ｃ上にありかつ最短距離のｚに写す写像である。

（固定ノイズ低減処理）
逐次近似法を用いた近似解の求め方の詳細を説明する前に、本実施例における焦点ぼけ画像群から固定ノイズを低減する処理である、固定ノイズ低減処理の全体概要について説明する。

図５は本実施例における固定ノイズ低減処理の処理フローを示す。固定ノイズ低減処理は、画像処理装置１１０で実行可能なプログラムとして実装されており、画像処理装置１１０は、撮像装置１２０や記憶装置１１３、他のコンピュータシステム１１４から焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）を取得することが可能である。

まず画像入力ステップＳ５０１では、画像処理装置１１０が、焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）を入力画像として取得し、前述の操作Ｔを用いて展開画像ｇを生成する。続いて、最適化処理ステップＳ５０２では、画像処理装置１１０が、式（１７）で定式化された最適化問題に最適化アルゴリズムを適用して求めた反復計算式を用いて、固定ノイズを低減した展開画像ｆ^（＊）を求める。詳細は後述する。ステップＳ５０２の計算はＣＰＵ３０１あるいはＧＰＵを用いた演算ユニットや専用の計算ハードウェア（不図示）などを用いて行い、結果はＲＡＭ３０２に出力する。
続いて、画像出力ステップＳ５０３では、画像処理装置１１０は、式（１６）に示すように操作Ｅを用いて、固定ノイズを低減した展開画像ｆ^（＊）から焦点ぼけ画像群ｘ_１ ^（＊），…，ｘ_ｋ ^（＊）を生成する。このようにして固定ノイズを低減した焦点ぼけ画像群ｘ_１ ^（＊），…，ｘ_ｋ ^（＊）は、出力画像として、ＲＡＭ３０２や記憶装置１１３、あるいは他のコンピュータシステム１１４に出力される。

（最適化アルゴリズム適用時の処理フロー）
図６は最適化処理ステップＳ５０２の内部処理を示すフローチャートである。
まず、初期値設定ステップＳ６０１では、画像処理装置１１０が、以降の計算で用いる初期値やパラメータを設定する。具体的には、ステップＳ６０２で用いる初期値や最適化アルゴリズムのパラメータ、ステップＳ６０３で用いる方向ベクトルｄやεの値、ステップＳ６０４で用いる反復回数の上限値ｎ_ｍａｘや解への収束判定基準の閾値ｅ_ｔｈが設定される。また、反復回数ｎのカウンタの初期値としてｎ＝０が設定される。

反復計算ステップＳ６０２では、画像処理装置１１０が、最適化アルゴリズムを適用して求めた反復計算式を用い、近似解ｆ^（ｎ）から近似解ｆ^{（ｎ＋１）}を求める。ただし、ｎ＝０の場合、近似解ｆ^（０）には初期値設定ステップＳ６０１で与えられた初期値が用いられる。

解空間制約ステップＳ６０３では、画像処理装置１１０が、ステップＳ６０２で求めた近似解ｆ^{（ｎ＋１）}が制約条件を満たすように、制約条件毎に距離射影（凸射影）と呼ばれる処理を行い、近似解ｆ^{（ｎ＋１）}を変更する。簡単に言うと、近似解が制約条件で定める解空間を外れる場合に、近似解を前記の解空間に引き戻す処理を行う。詳細は後述する。

終了判定ステップＳ６０４では、画像処理装置１１０が、反復処理を終了するか否かを判定する。具体的には、画像処理装置１１０は反復回数ｎがステップＳ６０１で定めた上限値ｎ_ｍａｘに達したかを判定する。次に、画像処理装置１１０は、解への収束判定基準（例えば前述の反復計算前後の近似解の差のＬ２ノルム）を計算し、ステップＳ６０１で定めた閾値ｅ_ｔｈを下回ったかを判定する。前記２つのいずれの判定条件も満たさない場合には反復回数更新ステップＳ６０５に進み、現在の反復回数ｎに１を加え、再び反復計算ステップＳ６０２に戻る。いずれか一方の判定条件に合致すれば、近似解出力ステップＳ６０６に進む。近似解出力ステップＳ６０６では、画像処理装置１１０は、直近のステップＳ６０３で求めた近似解ｆ^{（ｎ＋１）}を最適解、すなわち展開画像ｆ^（＊）に設定し、最適化処理ステップＳ５０２を終了する。

（制約条件Ｃ１について）
図８Ａを用いて、制約条件Ｃ１が表す意味について説明する。図８Ａは制約条件Ｃ１、すなわち
による解空間の制約を示す図である。これは、任意の画素位置の画素補正ベクトルｕ_{［ｉ，ｊ］}＝ｆ_{［ｉ，ｊ］}−ｇ_{［ｉ，ｊ］}が方向ベクトルｄと平行である、という制約を表す。

例として３枚の画像から構成される焦点ぼけ画像群を想定する。図８Ａは、元の焦点ぼけ画像群ｙ_１，ｙ_２，ｙ_３から生成した展開画像ｇから求めた画素位置（ｉ，ｊ）での画素ベクトルｇ_{［ｉ，ｊ］}を原点として、展開画像ｆから求めた画素ベクトルｆ_{［ｉ，ｊ］}が存在しうる３次元の解空間を表す。ｆ_{［ｉ，ｊ］}＝（（ｘ_１）_ｉ，ｊ，（ｘ_２）_ｉ，ｊ，（ｘ_３）_ｉ，ｊ）である。尚、焦点ぼけ画像群を構成する画像の枚数がＭ枚のときはＭ次元の解空間となる。破線は方向ベクトルｄ＝（１，１，１）を表す。

制約条件Ｃ１が存在しない場合、画素補正ベクトルｕ_{［ｉ，ｊ］}＝ｆ_{［ｉ，ｊ］}−ｇ_{［ｉ，ｊ］}は様々な方向を取りうる。しかし制約条件Ｃ１が存在する場合、画素補正ベクトルｕ_{［ｉ，ｊ］}は方向ベクトルｄに平行、すなわち方向ベクトルｄの実数倍となる必要がある。この制約条件Ｃ１を満足するように、画素ベクトルｆ_{［ｉ，ｊ］}を変更し、ｆ′_{［ｉ，ｊ］}とする。

変更前の画素ベクトルｆ_{［ｉ，ｊ］}と変更後の画素ベクトルｆ′_{［ｉ，ｊ］}の距離を最短とする画素補正ベクトルｕ_{［ｉ，ｊ］}は、以下に示すベクトルの正射影（直交射影）の式を用いて計算できる。
ただし、＜・，・＞はベクトルの内積を表す。

制約条件Ｃ１を満たす画素ベクトルｆ′_{［ｉ，ｊ］}は、式（１９）のとおり画素補正ベクトルｕ′_{［ｉ，ｊ］}と画素ベクトルｇ_{［ｉ，ｊ］}の足し算で生成できる。
ただし、式（１８），（１９）の計算は全ての画素（ｉ，ｊ）に対して行う必要がある。

（制約条件Ｃ２について）
続いて、図８Ｂを用いて制約条件Ｃ２が表す意味について説明する。図８Ｂは制約条件Ｃ２、すなわち
による解空間の制約を示す図である。この制約は、入力画像（展開画像ｇ）に対する輝度変更の総量が上限値ε以下であるという制約である。

図８Ｂでのｆ′_ｖは、制約条件Ｃ１による解空間の制約を満たす近似解ｆ′のベクトル表現ｆ′_ｖ＝ｖｅｃ［ｆ′］を表し、ｇ_ｖは展開画像ｇのベクトル表現ｇ_ｖ＝ｖｅｃ［ｇ］を表す。ｕ′_ｖはベクトルｆ′_ｖとｇ_ｖの差を表す。

制約条件Ｃ２を満たすｆ′_ｖの範囲は、ｇ_ｖを原点とする半径εの多次元（Ｎ^２Ｍ次元）の超球の内部となる。半径εの値は、前述したように撮影に用いたイメージセンサの固定ノイズの標準偏差σから計算する。固定ノイズの標準偏差σはイメージセンサの特性分析から事前に求めることができる。本実施例では加法性の固定ノイズを想定するため、
とする。

制約条件Ｃ２を満たす閉凸集合Ｃ_２への距離射影は以下の式で計算できる。
すなわち、ｕ′_ｖが半径εの超球の内部に存在しない場合にのみ、ｆ′_ｖは、ｆ′_ｖとｇ_ｖを結ぶ線分と、半径εの超球との交点ｆ″_ｖに射影される。

図８Ｂに示すように、閉凸集合Ｃ_２への距離射影では、変更前後のｕ′_ｖとｕ″_ｖではベクトルの方向は変化していない。つまり、制約条件Ｃ１で求めた各画素の画素補正ベクトルの方向は変化せず、方向ベクトルｄと平行のままであり、制約条件Ｃ１、Ｃ２を同時に満たす。

（解空間制約ステップＳ６０３）
続いて、図７を用いて解空間制約ステップＳ６０３の内部処理フローについて説明する。図６のフローチャートの順序では反復計算ステップＳ６０２が先になるが、解空間制約ステップＳ６０３は本実施例の特徴的な構成であるため、先に説明する。

まず差分展開画像計算ステップＳ７０１では、画像処理装置１１０が、近似解ｆ^{（ｎ＋１）}と展開画像ｇの差分展開画像ｕ＝ｆ^{（ｎ＋１）}−ｇを計算する。次に、制約条件Ｃ１適用ステップＳ７０２では、画像処理装置１１０が、ステップＳ７０１で計算した差分展開画像ｕに対し、制約条件Ｃ１を満たすように閉凸集合Ｃ_１への距離射影を行い、差分展開画像ｕ′を生成する。詳細は後述する。
次に、制約条件Ｃ２適用ステップＳ７０３では、画像処理装置１１０が、ステップＳ７
０２で求めた差分展開画像ｕ′に対し、制約条件Ｃ２を満たすように閉凸集合Ｃ_２への距離射影を行い、差分展開画像ｕ″を生成する。詳細は後述する。
最後に、近似解更新ステップＳ７０４では、画像処理装置１１０が、
の計算を行い、ステップＳ７０３で求めた差分展開画像ｕ″と展開画像ｇから制約条件Ｃ１，Ｃ２を満たすｆ″を求め、その値を近似解ｆ^{（ｎ＋１）}に設定する。

なお、図７では、差分展開画像ｕを最初に計算し、差分展開画像ｕに対し２つの距離射影を適用した後、変更後の差分展開画像ｕ″から変更後の近似解ｆ^{（ｎ＋１）}を計算したが、処理手順はこれに限られない。例えば、差分展開画像ｕを計算することなく、近似解ｆ^{（ｎ＋１）}に対して、閉凸集合Ｃ_１への距離射影を求めた後、閉凸集合Ｃ_２への距離射影を求めても、同じ結果を得ることができる。

（制約条件Ｃ１適用ステップＳ７０２）
次に、制約条件Ｃ１適用ステップＳ７０２の詳細について説明する。図９は制約条件Ｃ１適用ステップＳ７０２の内部処理を示すフローチャートである。画素ループＳ９０１、Ｓ９０６はループ端記号を表し、全ての画素位置（ｉ，ｊ）に対する処理が終了するまでＳ９０２〜Ｓ９０５の処理を繰り返すことを示す。

方向ベクトル取得ステップＳ９０２では、画像処理装置１１０が、式（１８）の計算に必要な方向ベクトルｄを取得する。既に述べたように本実施例では方向ベクトルｄは全ての要素が１である要素数Ｍ個のベクトルである。

次に、画素補正ベクトル取得ステップＳ９０３では、画像処理装置１１０が、
の操作を行い、画素位置（ｉ，ｊ）における画素補正ベクトルを求める。
次に、正射影計算ステップＳ９０４では、画像処理装置１１０が、式（１８）を用いて、方向ベクトルｄに対する画素補正ベクトルｕ_{［ｉ，ｊ］}の正射影ベクトルｕ′_{［ｉ，ｊ］}を計算する。ｕ′_{［ｉ，ｊ］}は方向ベクトルの実数倍となる。

次に、差分展開画像更新Ｓ９０５では、操作Ｖと逆の操作により、計算した画素補正ベクトルｕ′_{［ｉ，ｊ］}の値を差分展開画像ｕ′に書き込んでいく。差分展開画像ｕ′と画素補正ベクトルｕ′_{［ｉ，ｊ］}の間には
という対応関係があるため、予めＲＡＭ３０２に確保した差分展開画像ｕ′のデータ領域に、対応する画素補正ベクトルｕ′_{［ｉ，ｊ］}のデータを書き込む操作で実現できる。

画素ループＳ９０６では全ての画素に対する処理が終了したかを判定し、全ての画素に対する処理が終了した場合には処理を終了する。全ての画素に対する処理が終了していない場合には再びＳ９０２に戻り、処理を繰り返す。

（制約条件Ｃ２適用ステップＳ７０３）
次に、制約条件Ｃ２適用ステップＳ７０３の詳細について説明する。図１０は制約条件Ｃ２適用ステップＳ７０３の内部処理を示すフローチャートである。

まずＬ２ノルム計算ステップＳ１００１では、画僧処理装置１１０が、差分展開画像ｕ′のベクトル表現ｕ′_ｖ＝ｖｅｃ［ｕ′］のＬ２ノルムを計算する。
次に、距離射影ステップＳ１００２では、画像処理装置１１０が、式（２１）を用いて、ｕ″_ｖ＝Ｂａｌｌ（ｕ′_ｖ，ε）を求める。すなわち、画像処理装置１１０は、ステップＳ１００１で求めたＬ２ノルムが
を満たすか否かを判定し、満たす場合には
とし、満たさない場合には
とする。
次に、差分展開画像更新ステップＳ１００３では、画像処理装置１１０が、ベクトル表現したｕ″_ｖを再び並べ替え、その結果を差分展開画像ｕ″とする。

（凸最適化問題への最適化アルゴリズムの適用）
図６の初期値設定ステップＳ６０１および反復計算ステップＳ６０２の処理を説明するため、式（１７）で表される凸最適化問題に最適化アルゴリズムを適用して解く方法について説明する。
式（１７）で表される凸最適化問題は様々な逐次近似法を用いて求めることができるが、本実施例では一例としてＰｒｉｍａｌ−ｄｕａｌ Ｓｐｌｉｔｔｉｎｇ Ｍｅｔｈｏｄ（以降、ＰＤＳ法と記述する）を用いた解法について説明する。

ＰＤＳ法は目的関数にＬ１ノルムのような微分できない関数を含む場合に有効な逐次近似法の１つである。ＰＤＳ法では以下の式（２５）で定式化される凸最適化問題を解くための反復計算式を求めることができる。
ただし、関数Ａ，Ｂは凸関数で、Ｌは線形演算子（Ｌｉｎｅａｒ Ｏｐｅｒａｔｏｒ）である。ＰＤＳ法によれば、式（２５）で定式化された凸最適化問題の反復計算式は以下となる。
式（２６ａ）のσや式（２６ｂ）のτは初期値設定ステップＳ６０１で定めるパラメータである。

式（２６ａ），（２６ｂ）中のＰｒｏｘで始まりσまたはτのパラメータおよび関数名の添え字がついた関数は、添え字の関数の近接写像（Ｐｒｏｘｉｍｉｔｙ Ｏｐｅｒａｔｏｒ）と呼ばれる関数であり、以下で定義される。
関数Ｂの近接写像も同様に定義されるため、その説明を省略する。

式（２６ａ）の関数Ａ^＊は関数ＡのＦｅｎｃｈｅｌ−Ｒｏｃｋａｆｅｌｌａｒ共役関数であり、以下の式（２８）で定義される。共役関数Ａ^＊の近接写像は、Ｍｏｒｅａｕ分解を用いることで、式（２９）のように関数Ａの近接写像を用いて計算できる。
また式（２６ｂ）のＬ^＊はＬのａｄｊｏｉｎｔ ｏｐｅｒａｔｏｒであり、以下の関係を満たす。

式（１７）と式（２５）を対応づけると、関数ＡはＬ１ノルム、関数Ｂは指示関数Ｉ、Ｌは∇（ナブラ）となる。またＬが∇のとき、Ｌ^＊は−ｄｉｖ（ｄｉｖ：発散）となることが知られている。
離散的な∇、離散的な発散（ｄｉｖ）は、以下の演算で求めることができる。
ｕ∈Ｒ^Ｎ×Ｎとなるデータに対する∇ｕ∈Ｒ^{Ｎ×Ｎ×２}は、
である（画素ごとにｉ，ｊ方向の微分値を持つベクトルが得られる）。

また、ｐ＝（ｐ^１，ｐ^２）∈Ｒ^{Ｎ×Ｎ×２}となるデータに対するｄｉｖ（ｐ）∈Ｒ^Ｎ×Ｎは、
である。

従って、ＰＤＳ法では以下に示す反復計算によって近似解を求めることができる。
ただし、近接写像のパラメータを分かりやすくするため、式（３１ａ）〜（３１ｃ）では近接写像を
で表現している。

続いて、式（３１ａ）〜（３１ｃ）の計算方法について述べる。

（Ｌ１ノルムの共役関数の近接写像の計算について）
Ｌ１ノルムの共役関数の近接写像は
で求まり、内部で以下のＬ１ノルムの近接写像を呼び出すことで計算できる。
ｚ_ｉがベクトルの場合、｜ｚ_ｉ｜は、ベクトルｚ_ｉのＬ２ノルムを求めた後に絶対値を求める計算に相当するため、｜ｚ_ｉ｜は、ベクトルｚ_ｉのＬ２ノルムと等しくなる。

つまり、上記のＬ１ノルムの近接写像のzにｘ／σを、λに１／σを代入することで、
式（３２）中の
の値が計算できる。

（指示関数Ｉの近接写像の計算について）
指示関数Ｉの近接写像は、τとは無関係な距離射影となり、以下の式で求まる。
すなわち、ｘの要素が０未満の要素は０、２５５を超える要素は２５５に変更し、それ以外では変更はしないことを表す。

（ＰＤＳ法を用いた場合の最適化処理フロー）
続いて、ＰＤＳ法を用いた場合の最適化処理フローについて説明する。既に図６を用いて逐次近似法を用いた最適化処理（図５のステップＳ５０２）について説明したため、差分となる個所のみ説明する。

まず、初期値設定ステップＳ６０１では、画像処理装置１１０が、ＰＤＳ法のパラメータであるσおよびτを設定する。また画像処理装置１１０は、反復計算を行うための初期値の一例として、ｚ^（０）＝０、ｔ^（０）＝ｇ、ｆ^（０）＝ｇを設定する。σやτはＰＤＳ法での近似解の収束要件に従った値、例えば、σ＝１０、τ＝０．９／（８σ）を設定する。収束要件は
となる。

続いて、反復計算ステップＳ６０２では、画像処理装置１１０が、式（３１ａ）および式（３１ｂ）を用いてｎ＝１以降のｆ^（ｎ）を求める。なお、ｎ＝１の場合にはステップＳ６０１で設定された初期値ｆ^（０）、ｔ^（０）、ｚ^（０）を用いる。次に、解空間制約ステップＳ６０３では、図７を用いて説明したように、画像処理装置１１０が、式（３１ｂ）で計算したｆ^{（ｎ＋１）}の値を制約条件Ｃ１およびＣ２を満たす値に更新する。なお、ＰＤＳ法を用いる場合には、近似解更新ステップＳ７０４においてｆ^{（ｎ＋１）}を更新した後、追加の処理として、更新したｆ^{（ｎ＋１）}を式（３１ｃ）に代入し、次の反復計算に必要なｔ^{（ｎ＋１）}を求める処理を行う。ステップＳ６０４以降の処理に関しては、既に説明した内容と同一であるため、説明は省略する。
以上で本実施例の固定ノイズ低減処理について述べた。

なお、式（１５）をＰＤＳ法で解く例として、式（１７）に示すように関数Ｂに閉凸集合Ｃ_{２５５，Ｍ}の指示関数Ｉ（ｆ）を当てはめて解く場合について説明したが、それ以外にも多くの変形例があることは言うまでもない。例えば、以下の式（１７′）に示すように、関数Ｂに制約条件Ｃ１と制約条件Ｃ２を満たす集合である閉凸集合Ｃ_１∩Ｃ_２の指示関数Ｉ′（ｆ）を割り当てて解くこともできる。

このとき、指示関数Ｉ′（ｆ）の近接写像は、閉凸集合Ｃ_１∩Ｃ_２への距離射影で求まる。そして、閉凸集合Ｃ_１∩Ｃ_２への距離射影は、閉凸集合Ｃ_１への距離射影を行った後に、閉凸集合Ｃ_２への距離射影を行うことで実装できる。それぞれの距離射影の計算方法は既に述べているため、説明は省略する。なお、式（１７′）の制約条件ｆ∈Ｃ_{２５５，Ｍ}を省き、展開画像ｆ^（＊）を求めた後に輝度を［０，２５５］の範囲に収める処理で代用する場合には制約条件なしの最適化問題とすることができる。

また、式（１５）をＰＤＳ法で解く場合の変形例として、ＰＤＳ法の関数Ａの一部に指示関数を含むように割り当てることで、制約条件ｆ∈Ｃ_{２５５，Ｍ}と制約条件ｆ∈Ｃ_１∩Ｃ_２を目的関数に取り込み、制約条件なしの最適化問題に変更することもできる。例えば、式（２５）において、Ｌ：ｆ→（∇ｆ，ｆ）、Ａ：（ｚ_１，ｚ_２）→||ｚ_１||_１＋Ｉ′（ｚ_２）、Ｂ：ｆ→Ｉ（ｆ）となるようにそれぞれ線形演算子Ｌ、関数Ａ，関数Ｂを当て
はめる。関数Ａでは関数Ａを構成する関数の変数が異なる（変数分離されている）ため、関数Ａの近接写像は関数Ａを構成する関数の近接写像に分解できる。それぞれの近接写像は既知であるため、反復計算式は求められる。なお、指示関数と関数Ａ，関数Ｂの対応のさせ方として、指示関数Ｉを関数Ａの一部に、指示関数Ｉ′を関数Ｂに対応づけても良く、この場合でも同様に反復計算式が求められる。

（その他の目的関数について）
これまで焦点ぼけ画像群全体の滑らかさを評価する指標として全変動を用いる場合について説明してきた。しかし、既に述べたように、全変動以外にも様々な指標を用いることができる。
例えば、Ｔｏｔａｌ Ｇｅｎｅｒａｌｉｚｅｄ Ｖａｒｉａｔｉｏｎ（以降、ＴＧＶと略す）という指標がある。ＴＧＶは画像の１次微分（あるいは差分）だけでなく、より高次の微分（あるいは差分に対する複数回の差分）に対する変動量も加味した指標となっている。そのため、ＴＧＶを目的関数の指標として用いると、全変動を用いる場合に比べ、画像のグラデーション領域に発生する階段状のアーティファクトを抑制することができる。

例えば、２次のＴＧＶは以下の式で表現される。
式（３５）の中には明示的には２次微分（Ｇ∇）は含まれないが、式（３５）は１次微分（∇）と２次微分（Ｇ∇）の２つの関数のｉｎｆｉｍａｌ ｃｏｎｖｏｌｕｔｉｏｎとみなすことができる。よって、式（３５）は１次微分のＬ１ノルムと２次微分のＬ１ノルムの和の最小値に近い値を与える。なお、α_１とα_２は１次微分と２次微分のバランスを制御するパラメータである。

式（１７）において全変動の代わりに式（３５）のＴＧＶを用いれば、式（１７）は以下のように変形できる。
式（３６）の目的関数と式（２５）のＬ，Ａ，Ｂを以下のように対応づけることで、式（３６）はＰＤＳ法で解くことができる。
（→は写像を表す）

式（３６）の反復計算式の求め方については本発明の本質的な部分ではなく、また画像
のランダムノイズ低減にＴＧＶを適用した公知文献にも示されているため、詳細は省略する。なお、ＴＧＶ以外にも高次の微分（差分）を用いる画像の滑らかさに関する指標は複数知られている。式（１７）の全変動の代わりにそのような指標を用いる方法も本発明の範疇とする。

（微分可能な指標）
画像全体の滑らかさを表す指標として、一例として以下に示す微分可能な関数を用いてもよい。
（ただし、βは微小な正の実数とする）
式（１５）や式（１７）において全変動の代わりに式（３７）で定めるＪ_β（ｆ）を用いる場合、Ｊ_β（ｆ）は微分可能な関数のため、反復計算式で微分を用いる最急降下法やニュートン法等の逐次近似法が適用できる。

式（１５）や式（１７）で定式化した凸最適化問題において、指示関数やラグランジュの未定乗数法などにより制約条件を目的関数に移すことができる。また公知の技術によって目的関数であるＴＶやＴＧＶを制約条件に移すことも可能である。また、微分できない目的関数を含む凸最適化問題を解く場合には、ＰＤＳ法以外にも様々なアルゴリズムを適用可能である。例えば、ｆｏｒｗａｒｄ−ｂａｃｋｗａｒｄ ｓｐｌｉｔｔｉｎｇ ｍｅｔｈｏｄ、Ｄｏｕｇｌａｓ−Ｒａｃｈｆｏｒｄ ｓｐｌｉｔｔｉｎｇ ｍｅｔｈｏｄ、交互方向乗数法（ＡＤＭＭ）などがある。

本実施例の特徴は、焦点ぼけ画像群を構成する各画像の対応画素に対して同じ補正量を減算するという制約の下、焦点ぼけ画像群全体の滑らかさを最大化するように、画素位置ごとの輝度変更量（ノイズ補正量）を最適化する点にある。なお、目的関数において、画像の滑らかさと逆の関係にある指標（例えば画像内での輝度の変動量など）を用い、「焦点ぼけ画像群全体の滑らかさを最大化」する代わりに「焦点ぼけ画像群全体での変動量の総和を最小化」してもよい。その概念を実現するどのような目的関数、制約条件による最適化問題の定式化、および前記最適化問題を解くどのような最適化アルゴリズムを用いても本発明の範疇である。

以上述べたとおり、本実施例の方法を用いれば、焦点ぼけ画像群に含まれるイメージセンサに起因する加法性の固定ノイズを高精度に低減することができる。このような固定ノイズが低減された焦点ぼけ画像群を利用すれば、焦点ぼけ画像群からフィルタ型方式で任意視点画像や任意焦点ぼけ画像を生成する際に問題となっていた画質の劣化を抑制することができる。

＜実施例２＞
実施例１では、焦点ぼけ画像群に対し、被写体の像の画像成分には依存しない加法性の固定ノイズが加わっている場合の固定ノイズ低減方法について説明した。しかし、固定ノイズは加法性に限らず、イメージセンサのセルごとの感度ばらつきや増幅器のゲインばらつきなどの元の信号データ、すなわち、被写体の像の画像成分に依存する固定ノイズも存在する。そのような固定ノイズは乗法性の固定ノイズと呼ばれる。乗法性の固定ノイズが加わった焦点ぼけ画像群ｙ_ｋ∈Ｒ^Ｎ×Ｎ（ｋ＝１，…，Ｍ）は以下の式で表現できる。

式（３８）のｎ_ｋ∈Ｒ^Ｎ×Ｎ（ｋ＝１，…，Ｍ）は乗法性の固定ノイズを表す。同一のイメージセンサの同一のセルで撮影した場合でも、被写体の像が変化すれば乗法性の固定ノイズの値（ｎ_ｋ）_ｉ，ｊは異なる。
式（３９）のｖ∈Ｒ^Ｎ×Ｎは乗法性の固定ノイズのノイズ源を表す。ｖの各画素の値（ｖ）_ｉ，ｊはイメージセンサの様々な製造誤差に起因するが、以降、（ｖ）_ｉ，ｊは出現頻度分布が平均０、標準偏差σ_ｍの正規分布に従うものとして扱う。ただし、（ｖ）_ｉ，ｊが正規分布に従わない場合でも本実施例の方法では固定ノイズの低減が可能である。

焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）に存在する乗法性の固定ノイズは、式（３８），（３９）から以下の式（４０）に変形できる。

実施例１では、図８Ａに示す固定ノイズを低減するための方向ベクトルｄとして、式（１３）に示すように要素が全て１のＭ次元のベクトルを設定した。しかし、実施例２では式（３８）で与えられる焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）から固定ノイズを低減する場合、式（４０）に示すように固定ノイズの強さは（ｙ_ｋ）_ｉ，ｊに比例する。したがって、方向ベクトルｄとして、焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）（あるいは展開画像ｇ）の画素ベクトル
を設定すればよいことが分かる。

また、実施例１では、図８Ｂに示す展開画像ｇ（すなわち焦点ぼけ画像群全体）での固定ノイズの総量を表すεとして
を用いた。このεは、制約条件Ｃ２で用いる輝度変更量の上限値である。
実施例２では乗法性の固定ノイズは式（３８）で表される。そのため、εとして
を用いる。ｙ_ａは、焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）の全画素の平均値である。
式（４２）は乗法性の固定ノイズｎ_ｋの各要素を１次元に並べたベクトルのＬ２ノルムの近似値の一例である。

以上より、乗法性ノイズの低減を実現するには、実施例１における初期値設定ステップＳ６０１でのεの計算式、および制約条件Ｃ１適用ステップＳ７０２の方向ベクトル取得ステップＳ９０２での方向ベクトルｄの取得方法を修正するだけで良いことが分かる。つまり、制約条件Ｃ１について、乗法性ノイズの場合は、各画像の対応画素の輝度値に同じ値（補正係数）を乗算するという制約条件を設定する。また、制約条件Ｃ２で用いる輝度変更量の上限値εについては、入力画像である焦点ぼけ画像群の輝度に応じて設定する（本実施例では式（４３）で示す平均輝度ｙ_ａに比例した値とする）。

実施例２の最適化処理ステップＳ５０２の内部処理について説明する。
初期値設定ステップＳ６０１では、画像処理装置１１０が、以降の計算で用いる初期値やパラメータを設定する。具体的には、反復計算ステップＳ６０２で用いる初期値や最適化アルゴリズムのパラメータ、終了判定ステップＳ６０４で用いる反復回数の上限値ｎ_ｍａｘや解への収束判定基準の閾値ｅ_ｔｈが設定される。また、反復回数ｎのカウンタの初期値としてｎ＝０が設定される。そして、画像処理装置１１０は、式（４２），式（４３）を用いてεの値を計算する。なお、方向ベクトルｄは固定ではないため、ステップＳ６０１では設定しない。
次の反復計算ステップＳ６０２は実施例１と同一の処理となるため、説明は省略する。

次の解空間制約ステップＳ６０３については、制約条件Ｃ１適用ステップＳ７０２の処理のみが実施例１と異なる。以下、実施例２の制約条件Ｃ１適用ステップＳ７０２について説明を行う。

図９に示すように、制約条件Ｃ１適用ステップＳ７０２では、画像処理装置１１０が、展開画像ｆ^{（ｎ＋１）}に対し、各画素位置（ｉ，ｊ）についてＳ９０２〜Ｓ９０５の処理を繰り返す。まず、方向ベクトル取得ステップＳ９０２では、画像処理装置１１０は、画素ループで処理対象となっている画素位置の画素ベクトルを取得し、それを方向ベクトルとする。すなわち、画像処理装置１１０は、
の計算を行い、画素位置ごとに異なる方向ベクトルを求める。

画素補正ベクトル取得ステップＳ９０３では、実施例１と同じく、画像処理装置１１０は、ｕ_{［ｉ，ｊ］}＝Ｖ［ｕ，（ｉ，ｊ）］の計算により、画素補正ベクトルｕ_{［ｉ，ｊ］}を求める。次に正射影計算ステップＳ９０４では、実施例１と同じく、画像処理装置１１０は、式（１８）を用いて、方向ベクトルｄと平行な画素補正ベクトルｕ′_{［ｉ，ｊ］}を求める。次に差分展開画像更新ステップＳ９０５では、実施例１と同様、画像処理装置１１０は、計算した画素補正ベクトルｕ′_{［ｉ，ｊ］}の値を差分展開画像ｕ′に書き込んでいく。なお、画素ループＳ９０１〜Ｓ９０６の処理は各々の画素で独立に計算できるため、並列処理をしても良い。

以降の処理（図６のステップＳ６０４〜Ｓ６０６）は実施例１と同様のため、説明は省略する。近似解として求められた展開画像ｆ^（＊）は後段の画像出力ステップＳ５０３に入力される。ステップＳ５０３では、画像処理装置１１０は、式（１６）に示すように操作Ｅを用いて、固定ノイズを低減した展開画像ｆ^（＊）から焦点ぼけ画像群ｘ_１ ^（＊），…，ｘ_ｋ ^（＊）を生成する。

本実施例の特徴は、焦点ぼけ画像群を構成する各画像の対応画素に対して同じ補正係数を乗算するという制約の下、焦点ぼけ画像群全体の滑らかさを最大化するように、画素位
置ごとの輝度変更係数を最適化する点にある。なお、目的関数において、画像の滑らかさと逆の関係にある指標（例えば画像内での輝度の変動量など）を用い、「焦点ぼけ画像群全体の滑らかさを最大化」する代わりに「焦点ぼけ画像群全体での変動量の総和を最小化」してもよい。
以上述べたとおり、本実施例の方法を用いれば、焦点ぼけ画像群に含まれるイメージセンサに起因する乗法性の固定ノイズを高精度に低減することができる。このような固定ノイズが低減された焦点ぼけ画像群を利用すれば、焦点ぼけ画像群からフィルタ型方式で任意視点画像や任意焦点ぼけ画像を生成する際に問題となっていた画質の劣化を抑制することができる。

＜実施例３＞
実施例１，２では撮影時の焦点ぼけ画像群全体を最適化することで固定ノイズを低減する方法について説明した。しかし、焦点ぼけ画像群を構成する画像の枚数Ｍが増加するほど、焦点ぼけ画像群全体の最適化に必要な計算負荷（演算量やメモリ）が増加する。

従って、本実施例では、焦点ぼけ画像群を構成する画像から一部の画像を選択し、最適化することで固定ノイズを推定し、その後焦点ぼけ画像群全体から固定ノイズを低減する方法について述べる。本実施例の方法を用いれば、焦点ぼけ画像群を構成する画像の枚数Ｍが多数である場合でも計算負荷を抑え、高速に処理ができる。

（加法性の固定ノイズの低減）
図１１は本実施例における固定ノイズ低減処理を示すフローチャートである。
はじめに、加法性の固定ノイズが加わった焦点ぼけ画像群全体から固定ノイズを低減する場合について説明する。なお、本実施例で述べる固定ノイズ低減処理においては、あらかじめ、乗法性の固定ノイズ低減か加法性の固定ノイズ低減を行うかをユーザに設定させ、その設定に従って、ステップＳ１１０１〜Ｓ１１０５の処理を行うとよい。

まず、画像入力ステップＳ１１０１では、画像処理装置１１０が、固定ノイズが加わった焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）を取得する。次に、画像選択処理ステップＳ１１０２では、画像処理装置１１０が、焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）のうちの一部のＭ_ｓ枚（Ｍ＞Ｍ_ｓ＞１）の画像を選択し、新たな焦点ぼけ画像群ｂ_ｋ∈Ｒ^Ｎ×Ｎ（ｋ＝１，…，Ｍ_ｓ）を生成する。本実施例では、焦点ぼけ画像群ｂ_ｋ（ｋ＝１，…，Ｍ_ｓ）を入力画像として用いる。画像処理装置１１０は、操作Ｔを用いて焦点ぼけ画像群ｂ_ｋ（ｋ＝１，…，Ｍ_ｓ）から展開画像ｇを生成する。詳細は後述する。

最適化処理ステップＳ１１０３では、画像処理装置１１０が、展開画像ｇから加法性の固定ノイズを低減した展開画像ｆ^（＊）を求める。最適化処理ステップＳ１１０３の内部処理は実施例１と同じとなるため、説明は省略する。

推定固定ノイズ生成処理Ｓ１１０４では、画像処理装置１１０が、展開画像ｆ^（＊）に対して操作Ｅを適用し、固定ノイズを低減した焦点ぼけ画像群ａ_ｋ∈Ｒ^Ｎ×Ｎ（ｋ＝１，…，Ｍ_ｓ）を生成する。次に、画像処理装置１１０は、以下に示す式から推定固定ノイズｎ′∈Ｒ^Ｎ×Ｎを生成する。

続いて、推定固定ノイズ補正処理Ｓ１１０５では、画像処理装置１１０が、式（４５）
を用いて焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）から推定固定ノイズｎ′を取り除く。これにより、固定ノイズを低減した焦点ぼけ画像群ｘ_ｋ ^（＊）が生成される。

（サンプリングによる選択）
続いて、画像選択処理ステップＳ１１０１の処理について述べる。
焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）から、ノイズ推定に用いる焦点ぼけ画像群ｂ_ｋ（ｋ＝１，…，Ｍ_ｓ）を選択する方法には様々な方法がある。例えば最も簡単な例としては、予め決めた間隔でサンプリングする方法である。例えば、サンプリング間隔をｓｐとした場合、以下の式に示す関係
により焦点ぼけ画像群ｂｋ（ｋ＝１，…，Ｍ_ｓ）を生成することができる。ただし、Ｍ_ｓはＭ／ｓｐ以下の最大の整数である。例えば、Ｍ＝６４、ｓｐ＝８とすると、Ｍ_ｓ＝８となり、焦点ぼけ画像群ｂｋ（ｋ＝１，…，Ｍ_ｓ）は
となる。

なお、サンプリング間隔は均等でなくてもよく、不均等なサンプリングも可能である。本実施例で説明する固定ノイズ低減方法においては、図４Ｂで説明したように、焦点ぼけを伴うサンプリング位置を選択することが望ましい。

（サンプリングと加重平均）
またランダムノイズの影響を抑えるために、サンプリング間隔ｓｐで指定される位置の前後ｍｐ枚の画像に対して単純な平均や加重平均を求めた画像を生成し、その画像を固定ノイズ推定に用いても良い。その場合、焦点ぼけ画像群ｂｋは、例えば、以下の式で求めることができる。
ただし、ｗ（ｌ）は加重平均の重みを表す係数であり、
とする。なお、式（４６）の（（ｋ−１）・ｓｐ＋ｌ）ｍｏｄＭ＋１は、１≦ｋ・ｓｐ＋ｌ≦Ｍの範囲を外れた場合に巡回した位置を選択することを表す。

（合焦位置近傍の除外）
図４Ｂで述べた想定とは異なり、焦点ぼけ画像群のぼけ量が小さく、焦点ぼけ画像群の端でも被写体の像が十分にぼけていない場合がある。焦点ぼけ画像群の各画像において被
写体の像が同じような輝度で現れていると、最適化処理において被写体の像成分と固定ノイズ成分との区別が困難となり、結果的にステップＳ１１０４での推定固定ノイズの推定精度が低下する場合がある。

そのような場合、焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）を構成する画像のうち、被写体に合焦した位置の近傍の画像を除外することで、推定固定ノイズの推定精度を向上させることができる。そのため、ステップＳ１１０２では、焦点ぼけ画像群ｙ_ｋ（ｋ＝１，…，Ｍ）のうちから鮮鋭度が相対的に高い画像を除外する（又は鮮鋭度が相対的に低い画像を選択する）ことにより、焦点ぼけ画像群ｂ_ｋ（ｋ＝１，…，Ｍ_ｓ）を生成すると良い。

例えば、合焦位置近傍の画像の値が低くなるような以下の関数
を定め、Ｓｅｌ（ｋ）の値が所定の閾値ｓ_ｔｈ以上となる複数の画像ｙ_ｋを選択し、焦点ぼけ画像群ｂ_ｋ（ｋ＝１，…，Ｍ_ｓ）を生成する。なお、式（４７）において、ｓｔｄ（ｙ_ｋ）は画像ｙ_ｋの標準偏差を求める関数、αはＳｅｌ（１）〜Ｓｅｌ（Ｍ）の合計を１とするための定数、β（＞０）は関数パラメータである。

Ｓｅｌ（ｋ）の値が閾値ｓ_ｔｈ以上となる複数の画像ｙ_ｋに対して、前述のサンプリング間隔を定めて選択する手法や加重平均を求めて選択する手法を適用することで、固定ノイズ推定に用いる焦点ぼけ画像群の枚数をさらに減らしても良い。

（乗法性の固定ノイズの低減）
次に、乗法性の固定ノイズが加わった焦点ぼけ画像群全体から固定ノイズを低減する場合について説明する。前述のとおり、あらかじめ、乗法性の固定ノイズ低減か加法性の固定ノイズ低減を行うかをユーザに設定させ、その設定に従って、ステップＳ１１０１〜Ｓ１１０５の処理を行うとよい。

ステップＳ１１０１およびＳ１１０２の処理は加法性ノイズを低減する場合と同一であるため省略する。最適化処理ステップＳ１１０３では、画像処理装置１１０は、実施例２で述べた乗法性の固定ノイズの低減を行う。推定固定ノイズ生成処理ステップＳ１１０４では、画像処理装置１１０は、以下の式に示す計算を用い、固定ノイズのノイズ源の推定値ｖ′∈Ｒ^Ｎ×Ｎを求める。
ただし、ｎ_ｓ ^{（ａｖｅ）}は推定する乗法性の固定ノイズの平均値、ａ^{（ａｖｅ）}は固定ノイズ低減後の焦点ぼけ画像群ａ_ｋ（ｋ＝１，…，Ｍ_ｓ）の平均画像で、それぞれ以下の式から求める。

あるいは式（４８）の代わりに、以下の式を用いて固定ノイズのノイズ源の推定値ｖ′を求めてもよい。以下の式は、輝度変更後の画像と輝度変更前の入力画像との比に基づいて乗法性の固定ノイズのノイズ源を推定することを意味する。

推定固定ノイズ補正処理ステップＳ１１０５では、画像処理装置１１０は、以下の式に示す計算を用いて乗法性の固定ノイズを低減する。
なお、以下の式を用いて乗法性の固定ノイズを低減しても良い。

本実施例で述べた方法を用いれば、焦点ぼけ画像群の枚数が増加した場合に計算負荷を抑えた固定ノイズの低減処理ができる。そのため、計算資源の限られた画像処理装置や専用演算装置でも固定ノイズ低減処理を行うことができる。また焦点ぼけ画像群での焦点ぼけが十分でない場合には一部の画像を除外することで固定ノイズの推定精度を高めることもできる。

また本実施例では、被写体の異なる複数の焦点ぼけ画像群から複数の推定固定ノイズを生成し、前記の複数の推定固定ノイズからさらに精度を高めた推定固定ノイズを生成することも可能である。精度の高い推定固定ノイズを用いて推定固定ノイズ補正処理（図１１のステップＳ１１０５）を行えば、より高画質な焦点ぼけ画像群が生成できる。

＜実施例４＞
実施例２では固定ノイズ低減処理における乗法性の固定ノイズを低減するための最適化処理の実現方法について説明した。本実施例では実施例２と異なる方法で乗法性の固定ノイズを低減する方法について説明する。なお、本実施例では表記の簡略化のため、乗法性の固定ノイズ（ｖ）_ｉ，ｊの標準偏差σ_ｍをσと記す。

固定ノイズが乗法性のノイズの場合、画像に対する適切な階調変換を施すことで、乗法
性のノイズを加法性のノイズに近づけることができる。階調変換としては、例えば、対数変換、Ａｎｓｃｏｍｂ変換などが適用できる。本実施例では一例として対数変換を用いる処理を説明する。
式（３８）および式（３９）から乗法性の固定ノイズが加わった焦点ぼけ画像群ｙ_１，…，ｙ_Ｍは以下の式で表現できる。
式（５１）の両辺の対数をとると、以下の式が導ける。
すなわち、焦点ぼけ画像群の対数をとることで、乗法性の固定ノイズは加法性のノイズとして扱うことができる。

図１２は本実施例における固定ノイズ低減処理を示すフローチャートである。図１２を用いて本実施例の処理を説明する。
まず、画像入力ステップＳ１２０１では、画像処理装置１１０が、焦点ぼけ画像群ｙ_１，…，ｙ_Ｍを入力し、操作Ｔを用いて展開画像ｇを生成する。輝度変換ステップＳ１２０２では、画像処理装置１１０が、展開画像ｇの各画素の輝度に対数変換を行い、展開画像ｇ_ｃを求める。

次に、最適化処理ステップＳ１２０３では、画像処理装置１１０は、実施例１の最適化処理ステップＳ５０２と同様の処理を行い、展開画像ｇ_ｃから加法性の固定ノイズを低減した展開画像ｆ_ｃ ^（＊）を求める。ただし、本実施例では、初期値設定ステップＳ６０１において、εの値として輝度変換後の固定ノイズのＬ２ノルムに相当する値を設定する（例えば、ε＝（Ｎ（Ｍ）^１／２）σ）。ｌｏｇ（１＋ｘ）はｘが０に近いとき、ｘに近似できるためである。

次に、輝度逆変換ステップＳ１２０４では、画像処理装置１１０は、展開画像ｆ_ｃ ^（＊）の各画素の輝度に対数変換の逆変換である指数関数を適用し、展開画像ｆ^（＊）を求める。最後に画像出力ステップＳ１２０５では、画像処理装置１１０は、式（１６）に示すように操作Ｅを用いて、展開画像ｆ^（＊）から固定ノイズを低減した焦点ぼけ画像群ｘ_１ ^（＊），…，ｘ_Ｍ ^（＊）を生成する。
以上の処理により、乗法性の固定ノイズを低減することができる。

続いて、乗法性の固定ノイズを精度良く低減する乗法性の固定ノイズ低減方法について説明する。

（対数変換によるノイズの平均、分散の変化）
画素毎の乗法性の固定ノイズ（ｖ）_ｉ，ｊをｗで表すと、式（５２）の右辺の第二項はｌｏｇ（１＋ｗ）と表せる。ｗが０に十分近い、即ちｗの標準偏差σが０に十分近い場合、ｌｏｇ（１＋ｗ）はマクローリン展開による１次近似式（ｌｏｇ（１＋ｗ）＝ｗ）で十分近似できるため、輝度を対数変換した後の固定ノイズは平均０、分散σ^２のままで変化しないと見なせる。そのため、前述のように加法性の固定ノイズ低減手法がそのまま適用できる。

しかし、ｗが０から離れるにつれて、画像の輝度を対数変換した後に得られる加法性の
固定ノイズ、すなわち、ｌｏｇ（１＋ｗ）の平均および分散は上記の近似から離れ、平均０、分散σ^２のノイズとして見なせなくなる。その結果、固定ノイズ低減性能は低下する。

階調変換を適用した後のノイズの平均・分散は、元のノイズの確率密度分布および階調変換の式によって決まる。イメージセンサの乗法性の固定ノイズ（ｗ）は画素毎の画素感度ばらつき（ＰＲＮＵ）に起因し、その確率密度分布は平均０、分散σ^２の正規分布Ｎ（０，σ^２）で表せることが知られている。しかし、乗法性の固定ノイズ（ｗ）の確率密度分布が正規分布であるとき、入力画像の輝度を対数変換した後の固定ノイズに相当する成分、即ち、ｌｏｇ（１＋ｗ）の平均や分散を解析的に求めることは難しい。そのため、マクローリン展開を用いて平均や分散の近似値を求めることが有効である。

ｌｏｇ（１＋ｗ）のマクローリン展開によるＮ次近似の式は
（ただし、｜ｗ｜＜１）
で表され、Ｎの数を大きくすることでより精度の高い近似が得られる。

表１に、ｗの確率密度分布が正規分布Ｎ（０，σ^２）に従う場合の、マクローリン展開による１次〜４次の近似式を使って求めたｌｏｇ（１＋ｗ）の平均μ_ｔ、分散σ_ｔ ^２の近似値（μ_ｔａ、σ_ｔａ ^２）を示す。
（表１では、確率変数ＷがＮ（０，σ^２）に従うとき、確率変数Ｗのｎ次モーメント（Ｗ^ｎの期待値）Ｅ［Ｗ^ｎ］は、ｎが奇数のとき０、ｎが偶数のときσ^ｎ（ｎ−１）・（ｎ−３）・…・３・１となる式を用いて、各次数での平均と分散の近似式を求めた。）

当然、次数Ｎを大きくするほど精度は高まるが、実用上は標準偏差σ＝０．２程度で十分近似できる次数Ｎを用いれば問題はない。（乗法性の固定ノイズが正規分布で表現できるのは１＋ｗのほとんど全てが正に収まる、即ちσ≦０．２程度に収まる場合である。）そのため、予め、計算や実験等により十分な精度が得られる近似次数Ｎ_ｍを求めるとよい。

例えば、次数ＮをＮ＋１に更新したときの平均の近似値（μ_ｔａ）、および分散の近似値（σ_ｔａ ^２）のそれぞれにおける差が所定の閾値以下になる次数ＮをＮ_ｍとする。ある
いは、後述する乗法性の固定ノイズの低減処理フローを用いて多数の画像を用いたシミュレーションを行い、実験的にＮ_ｍを定めてもよい。具体的には次数ＮとＮ＋１での平均・分散の近似値を用いて固定ノイズを低減したときに画質の差が知覚できないレベルに収まる次数ＮをＮ_ｍと定める。

そして、Ｎ_ｍ次の平均μ_ｔａと分散σ_ｔａ ^２を求める標準偏差σに関する多項式（例えば表１のような式）を画像処理装置１１０中のプログラムとして実装する。そうすれば、入力画像に含まれる乗法性の固定ノイズの標準偏差σ（あるいは分散σ^２）から、対数変換後の固定ノイズ（ｌｏｇ（１＋ｗ））の平均と分散の近似値を求めることができる。

（平均シフト分の補正）
ここで注意すべき点として、表１に示すように、２次以上の近似では対数変換後の固定ノイズの平均は０ではなくなるという現象がある。そのため、従来の加法性のノイズを低減する最適化処理を用いると、閉凸集合Ｃ_２（対数変換して得られた展開画像ｇ_ｃを中心とするＬ２ノルム球）内には固定ノイズを低減した解は見つからず、適切にノイズが低減できない場合が発生する。

これに対し、本実施例で述べた加法性のノイズ低減方式を適用して固定ノイズを低減するには、対数変換後に平均μ_ｔのずれを補正することで変換後の平均を０に補正する階調変換を用いる方法が考えられる。具体的には、乗法性の固定ノイズが加わった焦点ぼけ画像群の各画素の輝度に対して、
で表される階調変換ｆを適用することで実現できる。μ_ｔは乗法性の固定ノイズの標準偏差σに基づく値であるため、式（５４）は対数変換した輝度から、輝度に依存しないバイアス成分を減算する処理を表す。

乗法性の固定ノイズが加わった焦点ぼけ画像群に前記階調変換ｆを適用すると以下の式（５５）が得られる。

式（５５）において、右辺の第一項を除く部分、即ち、ｌｏｇ（１＋ｗ）−μ_ｔを新たな固定ノイズと見なせば、階調変換ｆ適用後の固定ノイズは平均０、分散σ_ｔ ^２を持つ加法性の固定ノイズに変換されることが分かる。従って、その後の加法性の固定ノイズ低減処理（最適化処理Ｓ１２０３）において、閉凸集合Ｃ_２のＬ２ノルム球の半径εをε＝（Ｎ（Ｍ）^１／２）σ_ｔとすれば、展開画像ｇ_ｃを中心とするＬ２ノルム球に解が見つかり、適切に固定ノイズが低減できる。なお、平均μ_ｔ、標準偏差σ_ｔはそれぞれ次数Ｎ_ｍでのμ_ｔａ、σ_ｔａを用いる。

一方、最適化処理Ｓ１２０３の後の輝度逆変換Ｓ１２０４では、各画素の輝度に対数変換の逆変換である指数関数を適用すればよく、バイアス成分μ_ｔに関する処理は不要である。なぜなら、最適化処理Ｓ１２０３で既に固定ノイズは低減されており、逆変換による固定ノイズの平均、分散の変化を考慮する必要はないためである。
上記の処理により、乗法性の固定ノイズｗが十分０に近くない場合でも、実施例１で述べた加法性の固定ノイズ低減処理を適用して精度良く固定ノイズが低減できる。

（性能を改善した乗法性固定ノイズの低減処理フロー）
続いて、図１２を用いて、性能を改善した乗法性の固定ノイズの低減処理を説明する。
まず、初期設定ステップ（不図示）において、十分な精度で近似できる次数Ｎ_ｍでの対数変換後の固定ノイズの平均μ_ｔと分散σ_ｔ ^２の各々のσに関する多項式から、各々の近似値μ_ｔａ、およびσ_ｔａ ^２を求める。前記の近似値は、前述のように画像処理装置１１０中のプログラムを用いて、入力画像の乗法性固定ノイズの標準偏差σから容易に計算できる。

続く、画像入力ステップＳ１２０１では、画像処理装置１１０は焦点ぼけ画像群ｙ_１，…，ｙ_Ｍを入力し、操作Ｔを用いて展開画像ｇを生成する。

次に、輝度変換ステップＳ１２０２では、画像処理装置１１０は展開画像ｇの各画素の輝度に対し、式（５４）で示す階調変換ｆを適用し、展開画像ｇ_ｃを生成する。
次に、最適化処理ステップＳ１２０３では、画像処理装置１１０は展開画像ｇ_ｃから加法性の固定ノイズを低減した展開画像ｆ_ｃ ^（＊）を求める。このとき、最適化処理ステップ内の初期値設定ステップＳ６０１において、εの値として階調変換ｆ適用後の固定ノイズのＬ２ノルムに相当する値を設定する。具体的にはε＝（Ｎ（Ｍ）^１／２）σ_ｔａとなる。

次に、輝度逆変換ステップＳ１２０４では、画像処理装置１１０は、展開画像ｆ_ｃ ^（＊）の各画素の輝度に対数変換の逆変換である指数関数を適用し、展開画像ｆ^（＊）を求める。
最後に、画像出力ステップＳ１２０５では、画像処理装置１１０は、式（１６）に示すように操作Ｅを用いて、展開画像ｆ^（＊）から固定ノイズを低減した焦点ぼけ画像群ｘ_１ ^（＊），…，ｘ_Ｍ ^（＊）を生成する。

なお、これまで乗法性の固定ノイズｗの確率密度分布が正規分布で表現できる場合について説明したが、それ以外の分布で表現される場合にも適用可能である。例えば、乗法性のノイズｗの確率密度分布がガンマ分布で表現できる場合、対数変換後のノイズの平均や分散の解析的な解が求まることが知られている。その場合、初期設定ステップ（不図示）では解析的に平均μ_ｔと分散σ_ｔ ^２を求める式を用いて解析的な解を計算すればよい。

なお、入力画像データの輝度が０となる場合には対数変換は適用できないが、画像の入出力処理に追加の例外処理を導入することで対処できる。例えば、画像入力ステップＳ１２０１において画像処理装置１１０内で輝度が０の画素に対してはその画素位置情報を記憶するとともに、輝度を１に変更する。（輝度０から１への変更は微小であるため、画質への影響はほとんどない。）また、画像出力ステップＳ１２０５では前記の画素位置情報を参照し、該当する画素の輝度を０に変更する。（輝度が０の場合、乗法性のノイズも０となるため、０への変更は合理的である。）

なお、本実施例の手法を用いれば、加法性や乗法性以外の固定ノイズが存在しても低減できる。例えば、固定ノイズの分散σ^２が平均輝度となるポアソンノイズの場合には、分散安定化変換（ｖａｒｉａｎｃｅ−ｓｔａｂｉｌｉｚｉｎｇ ｔｒａｎｓｆｏｒｍａｔｉｏｎ）の一種であるＡｎｓｃｏｍｂｅ ｔｒａｎｓｆｏｒｍを適用する。これによりポアソンノイズを加法性のノイズとして扱える。輝度変換ステップＳ１２０２では、画像処理装置１１０は、展開画像ｇの各画素の輝度に以下の変換式を適用し、展開画像ｇ_ｃを求める。

また、輝度逆変換ステップＳ１２０４では、画像処理装置１１０は、展開画像ｆ_ｃ ^（＊）の各画素の輝度に以下の逆変換式を適用し、展開画像ｆ^（＊）を求める。
この場合、最適化処理の中の初期値設定ステップＳ６０１ではεの値として輝度変換後の固定ノイズのＬ２ノルムに相当する値を設定する（例えば、ε＝Ｎ（Ｍ）^１／２）。Ａｎｓｃｏｍｂｅ ｔｒａｎｓｆｏｒｍでは輝度変換後のポアソンノイズの標準偏差は１となるためである。

＜実施例５＞
これまでの実施例では、焦点ぼけ画像群から固定ノイズを低減する方法について説明した。しかし、本発明の適用対象は必ずしも焦点ぼけ画像群に限定されず、動画像データにも適用可能である。
複数枚からなる画像群の同一の画素位置で被写体の像の画像成分（例えばテクスチャ）が大きく変化するならば、被写体の像の画像成分を打ち消す輝度変更は必ずしも画像群全体の滑らかさを向上させることにはならないため、固定ノイズの低減が可能である。

図１３Ａ、図１３Ｂに本実施例の固定ノイズ低減処理を用いて固定ノイズの低減が可能な動画像データの一例を示す。
図１３Ａはカメラを固定して撮影した泳ぐ魚の動画を構成する一部の画像群を示している。背景には海や水槽の水が映っている。画像の左にあるｔ＝０，ｔ＝１，ｔ＝２は撮影時刻の違いを表しており、魚が水中を移動している。図１３Ｂは時間とともにカメラを動かして、静止物体である花を撮影した動画を構成する一部の画像群を示している。背景は空である。ｔ＝０，ｔ＝１，ｔ＝２は撮影時刻の違いを表し、画像ごとに花の位置が変化している。図示しないがいずれの動画像データも３枚の画像（フレーム）だけでなく、時間変化する多数枚（Ｍ枚）の画像（フレーム）から構成されている。領域毎に見れば、図４Ｂで示した焦点ぼけ画像群の場合と同様に、少数の鮮鋭度の高い画像と多数の滑らかな画像から構成されることになる。また仮に、複数枚からなる画像群の中に鮮鋭度の高い画像が複数存在しても全く異なる被写体の像となる。上述のような画像群に対しても実施例１〜４で述べた固定ノイズ低減処理は適用できる。

以下、実施例１で説明した固定ノイズ低減処理を適用する場合について説明する。
画像処理装置１１０は、撮像装置１２０や記憶装置１１３、他のコンピュータシステム１１４から動画像データｙ_ｔ∈Ｒ^Ｎ×Ｎ（ｔ＝１，…，Ｍ）を取得することが可能である。画像処理装置１１０は、ＲＡＭ３０２やＲＯＭ３０３に格納されている固定ノイズ低減処理プログラムを用いて固定ノイズ低減処理を実行する。

以下、図５のフローチャートを用いて動画像データｙ_１，…，ｙ_Ｍから固定ノイズを低減する処理を説明する。
画像入力ステップＳ５０１では、画像処理装置１１０が、イメージセンサの固定ノイズが加わった動画像データｙ_１，…，ｙ_Ｍを入力画像として取得し、操作Ｔを用いて展開画像ｇを求める。次に、最適化処理ステップＳ５０２では、画像処理装置１１０が、展開画像ｇから加法性の固定ノイズを低減した展開画像ｆ^（＊）を求める。最後に、画像出力ステップＳ５０３では、画像処理装置１１０は、操作Ｅを用いて、展開画像ｆ^（＊）から固定ノイズを低減した動画像データｘ_１ ^（＊），…，ｘ_Ｍ ^（＊）を生成する。

動画像データでは時間方向の画像の枚数（フレーム数）Ｍが非常に大きくなるため、所定の時間毎又は所定のフレーム数毎に動画像データを区切り、それぞれを入力画像として、独立に固定ノイズ低減処理を行っても良い。
また実施例３で説明した図１１に示す固定ノイズ低減処理を適用しても良い。その場合、一部の画像を選択して推定固定ノイズを生成し、全体の動画像データから推定固定ノイズを低減するため、計算負荷を低減することができる。また実施例２で説明した固定ノイズ低減処理を適用することで、乗法性の固定ノイズを低減することもできる。
以上述べた本実施例で述べる固定ノイズ低減処理を用いれば、動画像データに含まれる固定ノイズを低減することができる。

＜他の実施例＞
上述した実施例１〜５は本発明の一具体例を示したものにすぎず、本発明の範囲をそれらの実施例の構成に限定する趣旨のものではない。各実施例の構成を適宜変更してもよいし、技術的な矛盾のない限り、実施例同士を組み合わせてもよい。

また、実施例１〜５で述べた固定ノイズ低減処理において、焦点ぼけ画像群や動画像データを構成する画像を複数のブロックに分割し、各ブロックで固定ノイズ低減処理を並列に行うことができる。これによって処理が高速化できる。
また実施例１〜５では焦点ぼけ画像群や動画像データが１チャネルからなる場合についてのみ説明したが、複数の色チャネルからなるカラー画像（例えばＲＧＢ）の場合にも適用可能である。例えば、色チャネルごとに固定ノイズ低減処理を適用することが可能である。また実施例２では、標準偏差が輝度に比例する乗法性の固定ノイズを低減する場合についてのみ説明したが、方向ベクトルｄとして輝度ではなく、輝度を階調変換した値を用いることで、様々な種類の固定ノイズの低減に対応できる。例えば、方向ベクトルｄの要素として輝度ではなく、輝度の平方根を用いれば、固定ノイズの分散が輝度に比例するような固定ノイズの低減にも対応できる。

また、実施例１や実施例２の初期値設定ステップＳ６０１で用いる初期値に、実施例３で計算した、固定ノイズを低減した焦点ぼけ画像群および前記から計算したデータを設定しても良い。その場合には、計算負荷を抑えつつ、かつ精度の良い固定ノイズの低減が実現できる。また、実施例３に限らず、他の固定ノイズ低減手法で固定ノイズを低減した焦点ぼけ画像群および前記から計算したデータを、初期値設定ステップＳ６０１の初期値に設定しても良い。

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

１１０：画像処理装置
２０７：撮像部

Claims (23)

1. コンピュータが、同一のイメージセンサで撮像することで取得された複数の入力画像のデータを取得する入力ステップと、
   コンピュータが、逐次近似法を用いた反復計算によって、前記複数の入力画像の画質を全体的に向上させる、各入力画像に対する輝度変更の最適解を求める最適化処理ステップと、
   を有し、
   前記最適化処理ステップでは、各入力画像の同じ位置の画素に対して共通の輝度変更を行うという制約の下で、各入力画像の画素ごとの輝度変更の最適解が求められる
   ことを特徴とする画像処理方法。
2. 前記複数の入力画像のデータは、撮像光学系の光軸方向に焦点位置を異ならせて被写体を撮像することで取得された複数の画像のデータである
   ことを特徴とする請求項１に記載の画像処理方法。
3. 前記撮像光学系は、両側テレセントリック光学系である
   ことを特徴とする請求項２に記載の画像処理方法。
4. 前記撮像光学系は、顕微鏡を構成する光学系である
   ことを特徴とする請求項２又は３に記載の画像処理方法。
5. 前記複数の入力画像のデータは、動画像を構成する複数のフレームのデータである
   ことを特徴とする請求項１に記載の画像処理方法。
6. 前記制約は、各入力画像の同じ位置の画素の輝度に対して同じ値を加算又は減算するという制約である
   ことを特徴とする請求項１〜５のうちいずれか１項に記載の画像処理方法。
7. 前記制約は、各入力画像の同じ位置の画素の輝度に対して同じ値を乗算又は除算するという制約である
   ことを特徴とする請求項１〜５のうちいずれか１項に記載の画像処理方法。
8. コンピュータが、各入力画像の各画素の輝度を階調変換する輝度変換ステップをさらに有し、
   前記最適化処理ステップでは、前記輝度変換ステップで階調変換された入力画像を用いて輝度変更の最適解が求められる
   ことを特徴とする請求項１〜５のうちいずれか１項に記載の画像処理方法。
9. 前記輝度変換ステップは、各入力画像の各画素の輝度を対数変換するステップを有することを特徴とする請求項８に記載の画像処理方法。
10. 前記輝度変換ステップは、対数変換された各入力画像の各画素の輝度を、前記イメージセンサの乗法性の固定ノイズの標準偏差に基づく値で補正するステップをさらに有する
    ことを特徴とする請求項９に記載の画像処理方法。
11. 前記最適化処理ステップでは、各入力画像の同じ位置の画素に対して共通の輝度変更を行うという制約に加え、前記複数の入力画像に対する輝度変更の総量が所定値をこえないという制約の下で、各入力画像の画素ごとの輝度変更の最適解が求められる
    ことを特徴とする請求項１〜６のうちいずれか１項に記載の画像処理方法。
12. 前記最適化処理ステップでは、前記所定値が、前記複数の入力画像の輝度に応じて設定される
    ことを特徴とする請求項１１に記載の画像処理方法。
13. 前記最適化処理ステップにおける逐次近似法では、前記複数の入力画像の画質が全体的に向上したかどうかを評価する指標を含む目的関数を用いる
    ことを特徴とする請求項１〜１２のうちいずれか１項に記載の画像処理方法。
14. 前記指標は、輝度変更後の画像の滑らかさを表す指標である
    ことを特徴とする請求項１３に記載の画像処理方法。
15. 前記指標は、輝度変更後の画像における各画素の輝度の微分の絶対値の積分を含む指標である
    ことを特徴とする請求項１３に記載の画像処理方法。
16. 前記指標は、輝度変更後の画像における各画素の輝度の１次よりも高次の微分の絶対値の積分を含む指標である
    ことを特徴とする請求項１３に記載の画像処理方法。
17. コンピュータが、前記最適化処理ステップで求められた最適解としての輝度変更を施した輝度変更後の画像と、輝度変更前の前記入力画像との差または比に基づいて、前記入力画像の各画素の輝度に含まれている前記イメージセンサに起因する固定ノイズを推定するステップをさらに有する
    ことを特徴とする請求項１〜１６のうちいずれか１項に記載の画像処理方法。
18. コンピュータが、前記イメージセンサで撮像された画像から、前記推定された固定ノイズを差し引くことにより、前記イメージセンサに起因する固定ノイズが低減された画像を生成するステップをさらに有する
    ことを特徴とする請求項１７に記載の画像処理方法。
19. 前記入力ステップでは、前記イメージセンサで撮像することで取得されたＭ枚の画像のうちからＭ_ｓ枚（Ｍ＞Ｍ_ｓ＞１）の画像を選択し、前記選択されたＭ_ｓ枚の画像を前記複数の入力画像とする
    ことを特徴とする請求項１〜１８のうちいずれか１項に記載の画像処理方法。
20. 前記入力ステップでは、Ｍ枚の画像のうちから鮮鋭度が相対的に高い画像を除外し、又は、Ｍ枚の画像のうちから鮮鋭度が相対的に低い画像を選択することにより、Ｍ_ｓ枚の画像を選択する
    ことを特徴とする請求項１９に記載の画像処理方法。
21. 前記最適化処理ステップでは、逐次近似法で用いる初期値に、前記複数の入力画像から前記イメージセンサに起因する固定ノイズを予め低減した複数の画像を設定する
    ことを特徴とする請求項１〜２０のうちいずれか１項に記載の画像処理方法。
22. 請求項１〜２１のうちいずれか１項に記載の画像処理方法の各ステップをコンピュータに実行させるためのプログラム。
23. 同一のイメージセンサで撮像することで取得された複数の入力画像のデータを取得する入力手段と、
    逐次近似法を用いた反復計算によって、前記複数の入力画像の画質を全体的に向上させる、各入力画像に対する輝度変更の最適解を求める最適化処理手段と、
    を有し、
    前記最適化処理手段は、各入力画像の同じ位置の画素に対して共通の輝度変更を行うという制約の下で、各入力画像の画素ごとの輝度変更の最適解を求める
    ことを特徴とする画像処理装置。