画像処理装置、画像処理方法、及び画像処理プログラム
【課題】画像に対する反復フィルタ処理の過程において、反復フィルタ処理を適切なタイミングで停止することができる画像処理装置、画像処理方法、及び画像処理プログラムを提供することである。
【解決手段】本実施形態は、画像処理装置及びその方法に関するものであり、該画像処理装置は、入力画像に対して反復フィルタ処理を実施するフィルタ装置と、前記反復フィルタ処理過程において各フィルタリング毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定する反復フィルタ処理停止装置とを具備する。
【解決手段】本実施形態は、画像処理装置及びその方法に関するものであり、該画像処理装置は、入力画像に対して反復フィルタ処理を実施するフィルタ装置と、前記反復フィルタ処理過程において各フィルタリング毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定する反復フィルタ処理停止装置とを具備する。
【発明の詳細な説明】
【技術分野】
【０００１】
本実施形態は画像処理分野に関し、具体的には、画像（例えば、医学画像）に対してノイズを低減するための画像処理装置及びその方法に関する。
【背景技術】
【０００２】
画像のノイズ低減は、画像処理における重要なアイテムの一つである。例えば、医学画像処理において、装置の信号検出性能や画像生成アルゴリズムが適合していない、又は、造影剤量が少ない等の諸原因によって、ＣＴ（Computed Tomography）装置やＭＲＩ（Magnetic Resonance Imaging）装置・超音波（Ultrasound）診断装置等の医療診断装置により取得された医学画像にノイズが含まれてしまうことがある。これら画像におけるノイズは、医者の診断に多大な影響を及ぼし、後に続く画像処理（例えば、画像分割、計測、ボリュームレンダリング（Volume Rending）等）の正確性にも影響を及ぼしてしまう。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】Ｐ．Ｐｅｒｏｎａ，Ｊ．Ｍａｌｉｋ，“Scale-Space and Edge Detection Using Anisotropic Diffusion”，Proceedings of IEEE，Computer Society Workshop on Computer Vision刊行，１９８７年１１月，第１６−２２頁
【発明の概要】
【発明が解決しようとする課題】
【０００４】
本発明が解決しようとする課題は、画像に対する反復フィルタ処理の過程において、反復フィルタ処理を適切なタイミングで停止することができる画像処理装置、画像処理方法、及び画像処理プログラムを提供することである。
【課題を解決するための手段】
【０００５】
実施形態に係る画像処理装置は、フィルタ装置と、反復フィルタ処理停止装置とを備える。フィルタ装置は、入力画像に対して反復フィルタ処理を実施する。反復フィルタ処理停止装置は、前記反復フィルタ処理の過程における各反復毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する。
【図面の簡単な説明】
【０００６】
【図１】図１は、一実施形態に係る画像処理装置の構成を示すブロック図である。
【図２】図２は、上記実施形態に係る画像処理方法を示すフロー図である。
【図３】図３は、本実施形態の一つの具体的実施例に係る画像処理装置の構成を示すブロック図である。
【図４】図４は、上記具体的実施例に係る画像処理方法を示すフロー図である。
【図５】図５は、反復フィルタ処理を停止する時間を確定するための一つの具体的例を示すフロー図である。
【図６】図６は、それぞれ異なる類型の画像に対して反復フィルタ処理過程におけるフィルタリング画像とノイズ画像の関係係数の変化曲線を示す図である。
【図７】図７は、本実施形態の実施例に基づく応用場面の一つの例を示す図である。
【図８】図８は、本実施形態の他の一具体的実施例に基づく画像処理装置の構成を示すブロック図である。
【図９】図９は、本実施形態の他の一具体的実施例に基づく画像処理装置の構成を示すブロック図である。
【図１０】図１０は、図８に示す実施例の画像処理方法を示すフロー図である。
【図１１】図１１は、図９に示す実施例の画像処理方法を示すフロー図である。
【図１２】図１２は、それぞれ反復フィルタ処理過程におけるフィルタリング画像とノイズ画像の分散の変化曲線を示す図である。
【図１３】図１３は、本実施形態の他の一具体的実施例に基づく画像処理装置の構成を示すブロック図である。
【図１４】図１４は、上記具体的実施例に基づく画像処理方法を示すフロー図である。
【図１５】図１５は、上記具体的実施例の一つの具体的応用例を示す図である。
【図１６】図１６は、本実施形態の他の一実施例に基づく画像処理装置の構成を示すブロック図である。
【図１７】図１７は、上記実施例に基づく画像処理方法を示すフロー図である。
【図１８】図１８は、画像における各画像単位の勾配値分布の例を示す図である。
【図１９】図１９は、本実施形態を実現するためのコンピュータの構成を示すブロック図である。
【発明を実施するための形態】
【０００７】
以下において、本実施形態の基本的理解のために、その一例の概要について説明する。ただし、この概要は、本実施形態のキーワードや重要な部分を確定するものではなく、本実施形態の範囲を限定するものでもない。その目的は、単に、本実施形態の概念を簡潔に説明するためのものであり、その後に論述する更に詳細な説明の序章にすぎない。
【０００８】
本実施形態に係る画像処理装置は、入力画像に対して反復フィルタ（iterative filter）処理を実施するフィルタ装置と、前記反復フィルタ処理の過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する反復フィルタ停止装置と、を備えたことを特徴とする。
【０００９】
本実施形態に係る画像処理方法は、入力画像に対して反復フィルタ処理を実施するステップと、前記反復フィルタ処理過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定するステップを含んだことを特徴とする。
【００１０】
また、本実施形態は、更に上述方法を実現するためのコンピュータプログラムを提供する。
【００１１】
更に、本実施形態は、更に少なくともコンピュータ読み取り可能な媒体形式のコンピュータプログラム製品を提供し、その上、上述方法を実現するためのコンピュータプログラムコードを提供する。
【００１２】
以下、本実施形態の構成や目的・特徴・優位性についての理解を助けるため、図面を参照しながら説明する。図面中の各構成は単に本実施形態の原理を示すためのものである。図面において、同じ又は類似する技術特徴や各構成については、同じ又は類似する附番や表記を用いて説明する。
【００１３】
以下、本実施形態について図面を参照しながら説明する。本実施形態の一つの図面又は一つの実施例において説明される構成や特徴は、一つ又は複数の他の図面又は実施例における構成や特徴と組み合わせることが可能である。なお、説明をわかりやすくするために、図面と説明において本実施形態と関係の無い内容や、当業者にとって既知の構成や処理に関する表示や説明については省略する。
【００１４】
画像処理を実施している時、通常は画像に対して反復フィルタ処理を実施して、できるだけ画像中のノイズを除去することができる。もし反復の回数が少なすぎると、ノイズ低減の効果は十分期待できるものにはならず、かと言って、もし反復の回数が多すぎると、画像における有用な情報まで平滑化されてしまい、画像が不鮮明になって、その後の処理に影響を及ぼす。そこで、本実施形態では、画像の反復フィルタ処理の過程において、反復フィルタ処理の適切なタイミングでの停止を実現することができる画像処理装置及びその方法を提供する。
【００１５】
ここで、本実施形態に基づく装置及び方法は、各類型の画像に対するノイズ低減処理に適用することが可能である。例えば、本実施形態に基づく装置及び方法は、医学画像に対するノイズ低減処理に適用することが可能である。ここでの医学画像は、医療画像診断装置で取得した被験者のデータに基づいて生成された画像である。ここでの医療画像診断装置は、Ｘ線画像診断装置、超音波画像診断装置、ＣＴ装置、ＭＲＩ装置、ＰＥＴ（Positron Emission Tomography）装置等に限定されない。図７は、医療画像診断装置を利用して医学画像の取得を示す図である。図７に示すとおり、各種画像生成装置（例えば、上記医療画像診断装置）に予め定められた或いはユーザが定義したスキャンプロトコルを持たせる。ＣＴ装置を例にして説明すると、例えば、心臓同期スキャン（cardiac synchronization）、腹部４相強調スキャン（abdominal 4 phases enhance）、ダイナミックスキャン等のプロトコルを持たせることが可能であり、図７に示すプロトコル１、プロトコル２、・・・・・、プロトコルＮ（Ｎ≧１）のように、いくつかのプロトコルを利用して、複数の画像シーケンス（例えば、図７に示すシーケンス１、シーケンス２、・・・・・、シーケンスＮ等）を得ることができる。それゆえ、処理を必要とする画像は、医療画像診断装置により対応するプロトコルに基づいて患者データを取得して生成された医学画像であって良い。
【００１６】
図１は、本実施形態の一つの実施例に基づく画像処理装置の構成を示すブロック図であり、図２は、この実施例に基づく画像処理方法を示すフロー図である。
【００１７】
図１に示すとおり、画像処理装置１００は、フィルタ装置１０１と反復フィルタ処理停止装置１０３を具備することが可能である。画像処理装置１００は、図２に示す方法を用いて画像フィルタ処理を実施することができる。以下、図２に示す方法フローを参考にしながら、画像処理装置１００における各構成の機能について説明する。
【００１８】
フィルタ装置１０１は、該画像処理装置に入力された入力画像に対して反復フィルタ処理を実施する（図２のステップ２０２）。実際の応用場面に基づいて、フィルタ装置１０１（即ち、ステップ２０２）は、いずれかの適切なフィルタリング方法、例えば、非線形拡散（Non-Linear Diffusion：ＮＬＤ）フィルタリング方法、線形フィルタリング方法、異方性拡散（Anisotropic Diffusion：ＡＴＤ）フィルタリング方法等を採用することができる。なお、ここではこれらに限定しない。
【００１９】
フィルタ装置１０１は、各フィルタ処理後、得られたフィルタリング画像を反復フィルタ処理停止装置１０３へ出力し、反復フィルタ処理停止装置１０３は、各フィルタ処理後に得られたフィルタリング画像に基づいて反復フィルタ処理を停止するか否かを判断する。具体的には、反復フィルタ処理停止装置１０３は、反復フィルタ処理過程におけるフィルタリング画像の反復回数に伴う増大変化の変化率に基づいて反復フィルタ処理を停止するか否かを判断することができる（図２のステップ２０６）。例えば、もし反復フィルタ処理停止装置１０３が、フィルタリング画像の変化速度が予め定められた閾値に等しいか或いはそれより小さいと判断した場合、これは、この時のフィルタ処理操作の画像に対する影響は小さくなり、つまり軽減されたノイズの程度は小さくなったと言える。そうなれば、フィルタ装置１０１に対して反復フィルタ処理を停止するように指示をすることが可能である。そうでなければ、フィルタ装置１０１に対して継続して次の反復を実施するように指示をする。
【００２０】
ここでのフィルタリング画像の“変化”は、直前のフィルタ処理で得られるフィルタリング画像とその前のフィルタ処理（“その前”とは、直前の反復の前に実施された、ある一回の反復フィルタリング処理（例えば、その一回前の反復等）を意味する。なおこれに限定されない）で得られたフィルタリング画像との間の差を意味する。画像に対する反復フィルタ処理過程において、反復回数の増大に伴って、各フィルタ処理毎に得られるフィルタリング画像とその前のフィルタ処理にて得られたフィルタリング画像と比較した変化（即ち、差）が反復回数に伴って小さくなる。つまり、反復回数の増大に伴って、フィルタリング画像の変化速度が緩慢になる。図１と図２に示す装置と方法はこの原理に基づいており、フィルタリング画像の変化速度に基づいて反復フィルタ処理を停止する時間を確定する。このような方法と装置を利用して、反復フィルタ処理を停止する適切な時間を有効に確定することができ、画像における大部分のノイズを除去できると共に、画像における有用な情報が不鮮明になってしまわないことも保証することがきる。さらに、人の手により反復フィルタ処理を停止する時間を設定するといったような煩雑さも回避することができ、画像処理の更なる自動化を実現することができる。
【００２１】
フィルタリング画像の反復過程における変化は、その異なる特徴や特性にまで反映することができ、例えば、各反復毎に得られるフィルタリング画像とノイズ画像との相関性の反復過程における変化、フィルタリング画像又はノイズ画像のある画像特徴（例えば、分散又は標準差）の値の反復過程における変化等にも反映できる。なお、これらに限らず、例えば、上記の各種変化の数学的組み合わせでも良い。以下、図３〜図１５を参照しながら、反復フィルタ処理過程におけるフィルタリング画像の変化速度に基づいて反復フィルタ処理の停止をするか否かの判断をする画像処理装置及びその方法について、更に具体的な一実施例について説明する。
【００２２】
図３は、一つの具体的実施例に基づく画像処理装置の構成を示すブロック図であり、図４は、その具体的実施例に基づく画像処理方法を示すフロー図である。
【００２３】
上述のとおり、反復フィルタ処理過程において、各フィルタ処理で得られるフィルタリング画像は、その前のフィルタ処理で得られるフィルタ画像と比較して変化が小さくなる、すなわち変化速度が緩やかになる。反復フィルタ処理の開始段階において、各フィルタ処理で軽減されるノイズは比較的多く、反復回数の増大に伴って、画像におけるノイズはどんどん少なくなり、画像自身もだんだん平滑になっていく。フィルタリング画像の変化速度はフィルタリング画像の勾配（gradient）に正比例し、画像自身がだんだん平滑になるのに伴って、画像の勾配の変化が小さくなるので、フィルタリング画像の変化速度は緩やかになる。それに対応して、フィルタリング画像におけるノイズ量の変化に伴って、それとフィルタ処理にて除去されるノイズ画像との相関性も変化が発生する。つまり、前記相関性の変化もフィルタリング画像の変化を反映し、相関性の変化速度もフィルタリング画像の変化速度を反映する。図６の（Ａ）、（Ｂ）及び（Ｃ）は、それぞれこの相関性の反復フィルタ処理過程における変化曲線を示している。図６の（Ａ）、（Ｂ）及び（Ｃ）における曲線は、それぞれ３つの異なる画像（図に示す例では、３つの画像はそれぞれ超音波診断装置（ＵＬ）、Ｘ線ＣＴ装置（ＣＴ）及び磁気共鳴診断装置（ＭＲＩ）で異なる人体部位を撮影して得られた医学画像である。）に対する反復フィルタ処理過程における各フィルタ処理にて得られるフィルタリング画像と相対するノイズ画像の相関係数の時間（反復回数）に伴う変化動向を示している。図６の（Ａ）、（Ｂ）及び（Ｃ）に示すとおり、反復フィルタ処理の開始段階において、フィルタリング画像とノイズ画像の相関係数の値は減少傾向を示し、この傾向はノイズが除去されればされるほど大きくなる。一定時間経過後、相関係数の曲線は底となり（図中の“ｘ”で表記された底点Ｔ１、Ｔ２）、そして底点の後、フィルタリング画像とノイズの相関係数は再び緩やかな増大傾向を示し、この傾向は底点の後、反復回数の増大に伴い、画像中の有用な情報もだんだん平滑化されて失われ、つまり、有用な情報も除去されてしまう、もしくはノイズ画像が入り込み、ひいては、フィルタリング画像とフィルタ処理後のノイズ画像との間の相関性を再び増大させることになる。既存の方法では、上記底点を探し出して、そこを反復フィルタ処理を停止する時間として確定している。この既存の方法では、フィルタリング画像とノイズ画像の相関係数の最小となる時間を反復フィルタ処理を停止する時間点として確定している（例えば、Ｐａｖｅｌ Ｍｒａｚｅｋの論文“Selection of Optimal Stopping Time for Nonlinear Diffusion Filtering”（“Scale-Space and Morphology in Computer Vision”刊行，第三回国際会議，２００１年を参照）（以下、この論文文献を、関連文献１と称す。）。図６の（Ａ）に示す状況においては、底点Ｔ１はただ一点存在しており、且つその底点の位置は正確な停止時間に対応しているので、このような状況においては、上述の既存方法を利用して正確な停止時間を確定することができる。ただし、上記の既存方法は、図６の（Ｂ）又は（Ｃ）に示す状況においては適用できない。なぜなら、例えば、図６の（Ｂ）に示す曲線において、底点Ｔ２が出現している位置は正確な停止時間より明らかに遅れており、且つ底点Ｔ２の後も相関係数曲線はフラットな状態を保持（即ち、底点Ｔ２の後、相関係数値はほとんど変化がない状態）しており、したがって、このような状況においてもし上述の既存方法を利用しようとすれば、確定した時間点を底点Ｔ２としたとしても、それは正確な停止時間ではない。また、例えば、図６の（Ｃ）に示す曲線においては、単調に減少し続け、底点が無いので、このような状況においてたとえ上述の既存方法を利用したとしても、相関係数の最小の時間点が探し出すことができず、正確な停止時間を確定することができない。
【００２４】
本出願の発明者は、図６の（Ａ）、（Ｂ）及び（Ｃ）の３つの曲線に示すとおり、底点の有無に関係なく、この３つの曲線は、反復フィルタ処理の開始段階において、全て単調減少を示し、減少変化率は時間の経過に伴って緩やかになっていくことを見出した。図３と図４に示す装置と方法では、フィルタリング画像とノイズ画像の間の相関性が反復フィルタ処理過程において減少するその変化率（相関性が最小となる時間点ではない）を利用して、反復フィルタ処理を停止する時間を確定する。
【００２５】
図３に示すとおり、画像処理装置３００は、フィルタ装置３０１と反復フィルタ処理停止装置３０３を備え、反復フィルタ処理停止装置３０３は、ノイズ画像取得部３０３−１、相関性取得部３０３−２及び停止判断部３０３−３を備える。画像処理装置３００は、図４に示す方法フローを用いて画像フィルタ処理を実施する。以下、図４に示す方法フロー（流れ）を参考にしながら、画像処理装置３００における各装置及び各部の機能について説明する。
【００２６】
フィルタ装置３０１と図１に示すフィルタ装置１０１は類似するので、いかなる適切な方法を用いて該画像処理装置に入力された入力画像に対して反復フィルタ処理を実施するかについては（図４のステップ４０２）、ここでは説明を省略する。
【００２７】
フィルタ装置３０１は、各フィルタ処理後、得られるフィルタリング画像をノイズ画像取得部３０３−１へ出力する。ノイズ画像取得部３０３−１は、各フィルタ処理後、得られるフィルタリング画像と入力画像に基づいてノイズ画像を推定する（図４のステップ４０６−１）。
【００２８】
一つの例として、各フィルタ処理後、得られるフィルタリング画像と入力画像との間の差分画像を計算し、それをノイズ画像とすることができる。入力画像をＩ0、直前のフィルタ処理後のフィルタリング画像をＩstとした場合、対応するノイズ画像Ｉntは、以下の式で推定できる。
【００２９】
【数１】
【００３０】
なお、各フィルタ処理後、いずれかの適切なその他の方法を用いてノイズ画像を推定することができる（例えば、式（１）に基づいて計算して得られた差分画像の絶対値、即ち、Ｉnt＝｜Ｉ0−Ｉst｜をノイズ画像とすることができる。もちろん他の方法によりノイズ画像を指定しても良いが、ここでは、いちいち列挙して詳述はしない）。
【００３１】
ノイズ画像取得部３０３−１は、得られたノイズ画像を相関性取得部３０３−２へ出力し、相関性取得部３０３−２は、該ノイズ画像と直前の反復にて得られたフィルタリング画像の相関性を推定し（図４のステップ４０６−２）、推定結果を停止判断部３０３−３へ出力する。
【００３２】
一つの例として、両者画像間の相関性を表徴するため、ノイズ画像とフィルタリング画像の間の相関係数を計算する。例えば、両者画像間の相関係数は、以下の式を用いて計算することができる。
【００３３】
【数２】
【００３４】
ここで、ρtは、直前の反復で得られたフィルタリング画像Ｉstと対応するノイズ画像Ｉntとの間の相関係数を表し、Ｅ（・）は、平均値関数（mean value）を表す。
【００３５】
なお、相関係数以外の他のパラメータを用いて両者画像間の相関性を表して、上記の式（２）に対していずれかの形式の数学的変化を求めるようにしても良く、それゆえ、本実施形態は上述の例に限定されることは無い。
【００３６】
停止判断部３０３−３は、各フィルタ処理で得られたフィルタリング画像とノイズ画像との間の相関性が前記反復フィルタ処理過程において減少するその変化率に基づいて反復フィルタ処理を停止するか否かを確定する（図４のステップ４０６−３）。もし、停止判断部３０３−３が、相関性の減少変化率が既に予め定めた閾値に等しいか又はそれより小さいと判断した場合、フィルタ装置３０１に対して反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置３０１に対して次の反復の実施を指示する。
【００３７】
図５は、相関性の減少変化率を表す特徴値に基づいて、反復の停止時間を確定する一例を示している。図５に示すとおり、ステップ４０６−３１において、停止判断部３０３−３は、現在のフィルタ処理で得られたフィルタリング画像とノイズ画像との間の相関性の値及びその前のフィルタ処理後に推定された相関性の値に基づいて、相関性が反復回数に伴って増大又は減少する変化率を表すのに用いる特徴値を計算する。その後、ステップ４０６−３２において、停止判断部３０３−３は、該特徴値と予め定められた閾値との間で予め定めた条件を満足するか否かを判断し、反復フィルタ処理を停止するか、そうでなければ、次の反復を継続して実施するようにする。
【００３８】
一つの具体的例として、前記特徴値は、フィルタリング画像とノイズ画像の相関性の反復回数に伴う変化曲線（図６の（Ａ）、（Ｂ）又は（Ｃ）に示す相関係数曲線）の傾斜角度である。例えば、前記傾斜角度θは、以下の式を用いて計算することができる。
【００３９】
【数３】
【００４０】
ここで、ρtは、時間ｔにおいて得られたフィルタリング画像Ｉstと対応するノイズ画像Ｉntの間の相関係数を表し、ρt-Δtは、時間ｔ−Δｔにおいて得られたフィルタリング画像Ｉs(t-Δt)と対応するノイズ画像Ｉn(t-Δt)の間の相関係数を表し、Δｔは、時間間隔を表す。
【００４１】
その他の一具体的例として、前記傾斜角度は弧度（ラジアン）を用いても表現できる。
【００４２】
【数４】
【００４３】
停止判断部３０３−３が、計算して得られた傾斜角度θrad又はθが、予め定められた閾値に等しい或いはそれより小さい場合、フィルタ装置３０１に対して、反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置３０１に対して、次の反復の実施を指示する。
【００４４】
その他の具体例として、前記特徴値は、フィルタリング画像とノイズ画像の関連性の反復回数に伴う変化曲線（図６の（Ａ）、（Ｂ）又は（Ｃ）に示す相関係数曲線）の傾斜率であっても良い。例えば、前記傾斜率ｋslopeは、以下の式を用いて計算することができる。
【００４５】
【数５】
【００４６】
ここで、ρtは、時間ｔにおいて得られるフィルタリング画像Ｉstと対応するノイズ画像Ｉntとの間の相関係数を表し、ρt-Δtは、時間ｔ−Δｔにおいて得られるフィルタリング画像Ｉs(t-Δt)と対応するノイズ画像Ｉn(t-Δt)の間の相関係数を表し、Δｔは、時間間隔を表す。
【００４７】
停止判断部３０３−３が、計算して得られた傾斜率ｋslopeが、予め定められた閾値と等しい又はそれより小さい場合、フィルタ装置３０１に対して反復フィルタ処理を停止するように指示する。そうでなければ、フィルタ装置３０１に対して次の反復を実施するように指示する。
【００４８】
ここで、フィルタリング画像とノイズ画像との間の相関性の前記反復フィルタ処理過程における減少変化率を表す特徴値としては、上述の相関係数曲線の傾斜角度や傾斜率に限られない。前記特徴値として、相関性の減少変化率を表す他の特徴を用いても構わない。例えば、直前のフィルタ処理後に推定した相関係数とその前のフィルタ処理後に推定した相関係数ρt-Δtとの間の差（又は、該差の逆数）を前記特徴値として用いてもかまわない。上述の差が予め定められた閾値と等しい又はそれより小さい場合（又は、上述の差の逆数が予め定められた閾値に等しい又はそれより大きい場合）、反復フィルタ処理を停止する、そうでなければ、継続して次の反復を実施する。この他、上述の各種特徴値の計算方式に対して、いかなる形式の数学的変化を求めるようにしても良い。
【００４９】
その他、ここで、上述の実施例／例示において前記の予め定められた閾値は、実際の応用場面に基づいて予め確定されたものであり、異なる類型（又は特性）に対する画像により分析して得られた実験値や経験値であっても良い。異なる応用場面や採用する異なる特徴値パラメータに従って、該閾値は異なる値を採用することができる。例えば、前記入力画像が医療画像診断装置を介して得られたデータに基づいて生成された医学画像である場合、実験又は経験に基づいて、各種スキャンプロトコル又は画像生成装置の類型に従って予め対応する閾値を確定することができる。又は、複数のスキャンプロトコル又は画像生成装置に従う複数の閾値（例えば、プロファイル等の形式でもって）をメモリ部（該メモリ部は、画像処理装置に内在するメモリ装置であっても良いし、画像処理装置外部のメモリ装置であって該画像処理装置を経由して送るようにしても良い）へ保存する。このように、反復フィルタ処理過程において、入力画像のスキャンプロトコル又は画像生成装置の類型に基づいて対応する閾値を取得して（例えば、メモリ部から取り込む）、相関性の減少変化率を判断する根拠とすることができるので、更に画像処理の自動化の程度を向上させることができる。
【００５０】
図３〜５に示す画像処理装置又はその方法を参照する中において、フィルタリング画像とノイズ画像の相関性の減少変化率を利用することにより、反復フィルタ処理を停止する時間を有効に確定することができ、画像における大部分のノイズを除去できると共に、画像における有用な情報を平滑化により失われることがないようにすることも可能である。
【００５１】
図８及び図９は、それぞれ他の２つの具体的実施例に基づく画像処理装置の構造を示すブロック図であり、図１０及び図１１は、それぞれこの２つの具体的実施例に基づく画像処理方法を示すフロー図である。図３〜図５に示す実施例では、フィルタリング画像とノイズ画像の相関性の減少する変化率を利用して反復フィルタ処理過程におけるフィルタリング画像の変化速度を表す。図８〜図１１に示す実施例では、フィルタリング画像又はノイズ画像の予め定められた画像特徴の反復フィルタ処理過程における変化速度を利用してフィルタリング画像の変化速度を表す。
【００５２】
ここでの予め定められた画像特徴は、反復フィルタ処理過程におけるフィルタリング画像の変化のいかなる画像特徴値をも反映することができ、例えば、フィルタリング画像又はノイズ画像の分散又は標準差、フィルタリング画像又はノイズ画像の画素の平均勾配、フィルタリング画像又はノイズ画像のグレイスケール（gray scale）のダイナミックレンジ（即ち、最大グレイスケールと最小グレイスケールとの差）等であっても良い。ここでは、いちいち列挙して詳述はしない。
【００５３】
図１２の（Ａ）及び（Ｂ）は、それぞれ、反復フィルタ処理過程において、フィルタリング画像とノイズ画像の分散を例にした変化曲線を示している。図１２の（Ａ）に示すとおり、反復回数の増大に伴って、フィルタリング画像の分散は単調に減少し、且つ減少変化率は緩やかになっていく。図１２の（Ｂ）に示すとおり、反復回数の増大に伴って、ノイズ画像の分散は単調に増大し、且つ増大変化率は緩やかになっていく。上記の方法を利用してノイズ画像を推定することができる。説明は重複するので省略する。
【００５４】
以下、図８及び図１０を参照しながら、フィルタリング画像の予め定められた画像特徴の変化速度を利用してフィルタ処理停止時間を確定する実施例について説明する。
【００５５】
図８に示すとおり、画像処理装置８００は、フィルタ装置８０１と反復フィルタ処理停止装置８０３を備え、更に、反復フィルタ処理停止装置８０３は、画像特徴計算部８０３−４と停止判断部８０３−３を備える。画像処理装置８００は、図１０に示す画像処理方法を用いて画像フィルタ処理を実施することができる。
【００５６】
フィルタ装置８０１は、該画像処理装置へ入力された入力画像に対して反復フィルタ処理（図１０のステップ１００２）を実施するためのものであって、上記におけるフィルタ装置１０１及び３０１と類似するものであり、ここでは、詳細な説明は省略する。
【００５７】
フィルタ装置８０１は、各フィルタ処理後、得られたフィルタリング処理を画像特徴計算部８０３−４へ出力し、画像特徴計算部８０３−４を介して、各フィルタ処理で得られたフィルタリング画像の予め定められた画像特徴値を計算する（図１０に示すステップ１０００６−１）。例えば、フィルタリング画像の分散又は標準差を計算し、計算結果を停止判断部８０３−３へ出力する。停止判断部８０３−３は、反復フィルタ処理過程における前記画像特徴（例えば、分散又は標準差）値の反復回数の増加に伴う変化率（例えば、分散又は標準差の減少変化率）に基づいて反復フィルタ処理を停止するか否かを確定する。もし、停止判断部８０３−３が、前記画像特徴値の変化速度が予め定められた閾値と等しい又はそれよりも小さいと判断した場合、フィルタ装置８０１に対して反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置８０１に対して、次の反復の実施を継続する指示をする。
【００５８】
具体的例として、停止判断部８０３−３は、画像特徴計算部８０３−４の計算結果に基づいて、予め定められた画像特徴値の反復フィルタ処理過程における変化速度を表す特徴値を計算する。例えば、該特徴値は、予め定められた画像特徴値が反復回数の増加に伴って変化する曲線の傾斜角度又は傾斜率（例えば、図１２の（Ａ）に示す曲線の傾斜角度又は傾斜率であって、上述の公式（３）−（５）を用いて算出することができる。詳細については省略）、又は、該特徴値は、直前のフィルタ処理で得られたフィルタリング画像の予め定められた画像特徴値とその前のフィルタ処理で得られたフィルタリング画像の予め定められた画像特徴値との間の差（又は、該差の逆数）であって良い。ここでは、いちいち列挙して詳述はしない。その後、停止判断部８０３−３は、計算して得られた特徴値と予め定めた閾値との間で予め定めた関係を満足するか否かを判断し（例えば、傾斜角度、傾斜率又は上述差を特徴値として用いる状況下で、特徴値が予め定められた閾値に等しい又はそれより小さいか否かを判断する、又は、上述差の逆数を特徴値として用いる状況下で、特徴値が予め定められた閾値に等しい又はそれより大きいか否かを判断する）、フィルタ装置８０１に対して反復フィルタ処理の停止を指示する、そうでなければ、フィルタ装置８０１に対して次の反復を継続して実施するように指示をする。上述実施例／例示のように、ここでの閾値は、実際の応用場面に基づいて予め確定されるものである。異なる応用場面又はそこで用いられる異なる特徴パラメータに従って該閾値は異なる値をとる（例えば、図１２の（Ａ）に示す画像特徴値の変化曲線により予め閾値を確定し、これらの曲線は、異なる類型（又は特性）の画像に対して分析することにより取得できる。ここでは、更なる詳細な説明は省略する。）。なお、ここでは、それは、具体的数値に限定されない。
【００５９】
以下、図９と図１１を参照しながら、ノイズ画像の予め定めた画像特徴値の変化速度を利用してフィルタ処理の停止時間を確定する実施例について説明する。
【００６０】
図９に示すとおり、画像処理装置９００は、フィルタ装置９０１と反復フィルタ処理停止装置９０３を備える。反復フィルタ処理停止装置９０３は、図８に示す実施例と同様に、画像特徴計算部９０３−４と停止判断部９０３−３を備え、図８に示す実施例と異なるのは、反復フィルタ処理停止装置９０３は、更にノイズ画像取得部９０３−１を備えることにある。画像処理装置９００は、図１１に示すような画像処理方法を用いて画像フィルタ処理を実施することができる。
【００６１】
フィルタ装置９０１は、上記のフィルタ装置１０１や３０１と同様に、画像処理装置に入力された入力画像に対して反復フィルタ処理を実施するものである（図１１のステップ１１０２）。なお、ここでは重複する部分は説明を省略する。
【００６２】
フィルタ装置９０１は、各フィルタ処理後、得られたフィルタリング画像をノイズ画像取得部９０３−１へ出力する。ノイズ画像取得部９０３−１は、各フィルタ処理後に、得られたフィルタリング画像と入力画像に基づいてノイズ画像を推定する（図１１のステップ１１０６−１）。ノイズ画像取得部９０３−１は、上述の図３や図４に示す方法を用いてノイズ画像を推定する。なお、ここでは重複する部分は説明を省略する。
【００６３】
ノイズ画像取得部９０３−１は、得られたノイズ画像を画像特徴計算部９０３−４へ出力し、画像特徴計算部９０３−４を介してノイズ画像の予め定められた画像特徴値を計算する（図１１に示すステップ１１０６−２）、例えば、ノイズ画像の分散又は標準差を計算し、計算結果を停止判断部９０３―３へ出力する。
【００６４】
停止判断部９０３−３は、反復フィルタ処理過程においてノイズ画像の画像特徴値（例えば、分散や標準差）の反復回数に伴う増大変化の変化率（例えば、分散又は標準差の増大の変化率）に基づいて反復フィルタ処理を停止するか否かを確定する。もし、停止判断部９０３−３が、ノイズ画像の前記画像特徴値の変化速度が予め定めたレベルに達したと判断した場合、フィルタ装置９０１に対して反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置９０１に対して次の反復の実施を継続する指示をする。
【００６５】
具体的例として、停止判断部９０３−３は、画像特徴計算部９０３−４の計算結果に基づいて、予め定められた画像特徴値の反復フィルタ処理過程における変化速度を表す特徴値を計算する。例えば、該特徴値は、予め定められた画像特徴値の反復回数に伴い増大変化する曲線の傾斜角又は傾斜率を反映するもの（例えば、図１２の（Ｂ）に示す曲線の傾斜角度又は傾斜率であり、上述の式（３）−（５）を用いて計算できる。詳細は重複するので省略）、又は、該特徴値は、直前のフィルタ処理で得られたノイズ画像の予め定められた画像特徴値と、その前のフィルタ処理で得られたノイズ画像の予め定められた画像特徴値との間の差（又は、該差の逆数）である。ここでは、いちいち列挙して詳述はしない。その後、停止判断部９０３−３は、計算して得られた特徴値と予め定められた閾値との間で予め定められた関係を満足するか否かを判断し（例えば、傾斜角、傾斜率又は上述差を特徴値とする場合、特徴値が予め定められた閾値と等しい又はそれより小さいかどうかを判断する。上述差の逆数を特徴値とする場合、特徴値が予め定められた閾値と等しいか又はそれより大きいかどうかを判断する）、フィルタ装置９０１に対して反復フィルタ処理の停止を指示する、又は、フィルタ装置９０１に対して次の反復の実施を継続するよう指示する。上述実施例／例示のように、ここでは、前記の閾値は実際の応用場面に基づいて予め決められたものである。異なる応用場面とそこで用いられる異なる特徴パラメータに従って、該閾値は異なる値（例えば、図１２の（Ｂ）に示す画像特徴値の変化曲線により予め閾値を確定することができる。この曲線は、異なる類型（又は特性）の画像に対して分析を実施することにより得ることができる。）を用いることができる。なお、ここではその具体的数値に限定されない。
【００６６】
図８〜図１１に示す画像処理装置や方法において、反復フィルタ処理過程においてフィルタリング画像又はノイズ画像の予め定められた画像特徴値を利用することにより、反復フィルタ処理を停止する時間を有効に確定することができ、画像中の大部分のノイズを除去できると共に、画像中の有用な情報が平滑化により失うことを防ぐことができる。
【００６７】
具体的例として、図８〜図１１に示す画像処理装置又は方法を参考にして、以下の式により画像の分散を計算し、前記予め定められた画像特徴値とすることができる。
【００６８】
【数６】
【００６９】
ここで、Ｉは画像を表し（例えば、上述のフィルタリング画像やノイズ画像）、Ｄは該画像の分散を表す。
【００７０】
図１３は、他の一具体的実施例に基づく画像処理装置の構造を示すブロック図であり、１４は、該具体的実施例に基づく画像処理方法を示すフロー図である。上記の実施例はそれぞれノイズ画像とフィルタリング画像との間の相関性、又は、ノイズ画像（又はフィルタリング画像）の予め定められた画像特徴値の反復フィルタ処理過程における変化速度を利用し、そして図１３と図１４に示す実施例は両者を組み合わせることを考える。例えば、図６の（Ａ）、（Ｂ）、（Ｃ）に示すとおり、ノイズ画像とフィルタリング画像の間の相関性が反復フィルタ処理の開始段階は単調減少し、そのあと緩やかに増大（例えば、図６の（Ａ）及び（Ｂ）に示すとおり）、又は、緩やかな減少が続き（例えば、図６の（Ｃ）のとおり）、図１２の（Ｂ）に示すとおり、ノイズ画像の分散は反復フィルタ処理過程において単調に増大する。もし、両者が適切な方法により組み合わさる場合（両者が組み合わさった場合の具体的実施例としては、例えば、以下の文中で図１３−１５を参照しながら説明するようなものを挙げることができる。）、生成される組み合わせ特徴値の反復回数に伴う変化曲線も反復フィルタ処理過程において画像の変化を反映する。例えば、図６の（Ｂ）に示すような相関性曲線の底点が時間方向下流へ延びていく状況下では、ノイズ画像の分散は単調増大するので、組み合わせ特徴値の曲線は底点を時間方向上流へ適切な位置まで遡るようシフトする。また、例えば、図６の（Ｃ）に示すような相関性が単調減少する状況下では、ノイズ画像の分散は単調増大するので、組み合わせ特徴値の曲線には底点が出現すると共に該底点は正確な停止時間に対応する。図１５の（Ａ）は、ノイズ画像とフィルタリング画像の相関係数を示したものであり、ノイズ画像の分散及び両者の組み合わせ特徴値の反復回数に伴う変化を示す図である。図１５の（Ａ）において、縦軸は相関係数を表し、分散及び組み合わせ特徴の値を取り、横軸は反復時間を表す。図１５の（Ａ）に示すとおり、組み合わせ特徴値の曲線のｔ＝０．４の時間点にて底点（極小値）が出現し、前記組み合わせ特徴値の反復回数に伴って増大して出現する底点（即ち、極小値）を用いて、反復フィルタ処理を停止する時間を確定することができる。
【００７１】
説明をわかりやすくするために、以下では、ノイズ画像の分散を前記予め定められた画像特徴とし、ノイズ画像とフィルタリング画像の相関係数が前記相関性を表すものとして説明する。なお、その他の例としては、ノイズ画像のその他の画像特徴（例えば、標準差など）、又は、フィルタリング画像の分散や標準差の逆数など、又は、フィルタリング画像とノイズ画像の間の相関性のその他のパラメータでも良い。ここでは、いちいち列挙して詳述はしない。
【００７２】
図１３に示すとおり、画像処理装置１３００は、フィルタ装置１３０１と反復フィルタ処理停止装置１３０３を備え、更に、反復フィルタ処理停止装置１３０３は、ノイズ画像取得部１３０３−１、相関性取得部１３０３−２、画像特徴計算部１３０３−４、組み合わせ部１３０３−５及び停止判断部１３０３−３を備える。画像処理装置１３００は、図１４に示すような画像処理方法を用いて画像フィルタ処理を実施する。
【００７３】
上述のフィルタ装置１０１や３０１のように、フィルタ装置１３０１は、該画像処理装置に入力された入力画像に対して反復フィルタ処理を実施するものである（図１４のステップ１４０２）。ここでは、重複する部分については説明を省略する。
【００７４】
フィルタ装置１３０１は、各フィルタ処理後、得られたフィルタリング画像をノイズ画像取得部１３０３−１に出力する。ノイズ画像取得部１３０３−１は、各フィルタ処理後に得られたフィルタリング画像と入力画像とに基づいてノイズ画像を推定し（図１４のステップ１４０６−１）、ノイズ画像取得部１３０３−１は、上述の図３と図４に示す方法を参考にしてノイズ画像を推定する。ここでは、重複する部分については説明を省略する。
【００７５】
ノイズ画像取得部１３０３−１は、得られたノイズ画像を、それぞれ相関性取得部１３０３−２と画像特徴計算部１３０３−４へ出力する。相関性取得部１３０３−２は、該ノイズ画像と直前のフィルタ処理で得られたフィルタリング画像との相関性を推定し（図１４のステップ１４０６−２）、推定結果を組み合わせ部１３０３−５へ出力する。例えば、相関性取得部１３０３−２は、２つの画像間の相関性を表徴するため、ノイズ画像とフィルタリング画像との間の相関係数を計算する（例えば、上述の方法を用いて計算できる。重複する説明は省略。）。
【００７６】
画像特徴計算部１３０３−４は、ノイズ画像の予め定められた画像特徴値を計算し（例えば、図１４のステップ１４０６−４）、例えば、ノイズ画像の分散を計算し（例えば、上記式（６）を用いて）、計算結果を組み合わせ部１３０３−５へ出力する。
【００７７】
組み合わせ部１３０３−５は、相関性取得部１３０３−２により取得した相関性値と画像特徴計算部１３０３−４により計算して得られた分散とを組み合わせて、組み合わせ特徴値を生成する（例えば、図１４のステップ１４０６−５）。例えば、Ｃ1を単調増大するノイズ画像の分散曲線とし、Ｃ2をノイズ画像とフィルタリング画像の相関曲線とすると、組み合わせ特徴値は次のように表せる。
【００７８】
【数７】
【００７９】
ここで、Ｃcombineは組み合わせ特徴値を表し、ｆ（・）は曲線Ｃ1に対する数学的変化を表し、ｇ（・）は曲線Ｃ2に対する数学的変化を表す。
【００８０】
一例として、組み合わせ部１３０３−５は、相関性取得部１３０３−２により取得した相関性値と画像特徴計算部１３０３−４により計算して得られた分散値を計算して和を求め、該和を組み合わせ特徴値とすることができる。ただし、状況によっては、フィルタリング画像とノイズ画像の相関係数の値を取る範囲は、しばしばノイズ画像の分散の値を取る範囲と異なる。それゆえ、相関性値と分散値の和曲線は、上述の底点（極小値）が存在しない可能性がある。したがって、２つの曲線をマッチングさせるように調整する必要がある。以下、相関性取得部１３０３−２により得られた相関性値と画像特徴計算部１３０３−４により計算して得られた分散値とを組み合わせることについて、その一例を示す。
【００８１】
一例として、フィルタ処理の前に、画像処理装置（例えば、フィルタ装置１３０１）は、先ずオリジナルの入力画像に対してノーマライズ処理を行う。言い換えると、図１４のステップ１４０２の前に、入力画像のノーマライズ処理を行うステップを含む。このように、各反復の後で、ノイズ画像取得部１３０３−１は、ノーマライズされた入力画像を用いて得られた前記相関性値（相関性取得部１３０３−２により得られる）と前記分散値（画像特徴計算部１３０３−４により得られる）の和を求め、該和を組み合わせ特徴値とする。
【００８２】
その他の一例として、画像特徴計算部１３０３−４は、各反復の後、ノイズ画像取得部１３０３−１にて推定されたノイズ画像に対してノーマライズを実施し（該例においては、オリジナル入力画像に対してノーマライズを実施する必要はない）、ノーマライズされたノイズ画像の分散を計算し、ノーマライズされたノイズ画像の分散を組み合わせ部１３０３−５へ出力する。ノイズ画像のダイナミック範囲（例えば、［Ｎmin，Ｎmax］を用いて画像のダイナミック範囲を表すことができる）を［０，１］へノーマライズする。異なる類型の入力画像に対しては、各反復で生成されるノイズ画像のダイナミック範囲は異なり、ノーマライズが一つの同じダイナミック範囲に基づいて実施されない状態を引き起こす。それゆえ、一つの標準ダイナミック範囲を指定し、各反復時に全てその標準ダイナミック範囲でもって［０，１］へのノーマライズを実施する。例えば、ある種類の入力画像に対して、相応する標準ダイナミック範囲を［Ｓmin，Ｓmax］と設定し（ここで、異なる入力画像に対しては、異なる標準ダイナミック範囲を用いることができる）、ノーマライズ処理は、以下の式（８）により実施することができる。
【００８３】
【数８】
ここで、Ｎはノイズ画像の画素値を表し、Ｎ’はノーマライズ処理後の画素値を表す。
【００８４】
ここで、実際の応用場面では、実際の必要性に応じて上述の標準ダイナミックレンジを確定するものであって、本実施形態を上記実施例に限るべきではない。例えば、医学画像の場合、標準ダイナミックレンジは、入力画像を取得した医療画像診断装置の異なる類型に基づいて確定するものである。その他の例では、入力画像のスキャンプロトコル等に基づいて確定しても良い。ここでは、いちいち列挙して詳述はしない。
【００８５】
該例において、相関性取得部１３０３−２は、ノイズ画像取得部１３０３−１により推定されたノイズ画像（ノーマライズされたノイズ画像ではない）に基づいて上記相関性を計算する。
【００８６】
その他の一例として、組み合わせ部１３０３−５は、計算して得られたノイズ画像の分散に対して変換処理を実施し、その後、変換後の分散と相関性取得部１３０３−２にて得られた相関性値の和を求め、該和を組み合わせ特徴値とする。例えば、組み合わせ部１３０３−５は、画像特徴計算部１３０３−４により計算して得られたノイズ画像の分散に対して変換処理を実施し、その後、変換処理後の分散値に相関性取得部１３０３−２により得られた相関性値を加え、得られた和を上記組み合わせ特徴値とする。例えば、組み合わせ部１３０３−５は、以下の式を介して、画像特徴計算部１３０３−４により計算して得られたノイズ画像の分散に対して変換処理を実施する。
【００８７】
【数９】
【００８８】
ここで、Ｖiは画像特徴計算部１３０３−４により計算して得られたノイズ画像の分散を表し、Ｖi’は組み合わせ部１３０３−５により分散値Ｖiに対して変換を実施して得られた変換値を表し、ｉ（ｉ＝１，２，・・・）は反復時間を表し、ｏｆｆｓｅｔは平行移動係数を表し、ｏｆｆｓｅｔ＝√Ｖ1のＶ1は一回目の反復後に画像特徴計算部１３０３−４にて計算して得られたノイズ画像の分散を表し、ａiは反復回数の増大に伴って増大する係数であり、該係数は１より大きく、例えば、ａi＝１．０＋ｃ×ｉ（ｃはゼロより大きい予め定められた定数）である。式（９）において、分散値Ｖiの平方根を求めることにより、変換後の分散値Ｖi’は大きくなり、平行移動係数ｏｆｆｓｅｔにより，変換処理後の分散曲線を平行移動し、ゼロから開始するようにする。係数ａiにより，変換処理後の分散曲線の増大変化率は増大し、反復回数の増大に伴って、該増大変化率も増大する。図１５の（Ｂ）は、上記変換後のノイズ画像の分散、ノイズ画像のフィルタリング画像との相関係数及び両者の和が反復回数に伴って変化する様子を示した図である。図１５の（Ｂ）に示すとおり、相関係数曲線と比べ、得られる和曲線の底点は時間方向上流へ遡るようシフトされる。何回もの検証により、もし、図１５の（Ｂ）に示す相関係数曲線の底点を反復フィルタ処理を停止する時間とすると、反復フィルタ処理を経た出力画像における有用な情報は曖昧となってしまう。もし、和曲線の底点を用いれば、反復フィルタ処理を経た出力画像においても、オリジナルの入力画像における有用な情報は保たれると共に、大部分のノイズは除去される。
【００８９】
その他、上記の式（９）におけるパラメータｃは実際の応用場面に基づいて予め確定した値であり、例えば、異なる類型（又は特性）の画像に対して分析を実施して得られる実験値又は経験値であって良い。上記の式（９）により、ｃ値が小さければ小さいほど、変換後の分散曲線はなだらかになり、反対に、ｃ値が大きければ大きいほど、変換後の分散曲線の増大変化率は速くなることを知ることができる。ただし、ｃ値は小さすぎてはいけない。小さすぎてしまうと、組み合わせ部１３０３−５により得られる組み合わせ特徴値の反復回数に伴う変化曲線は単調減少を引き起こす可能性があり、その場合、極小値が出現しないことになってしまう。かといって、ｃ値は大きすぎてもいけない。大きすぎると、組み合わせ部１３０３−５により得られる組み合わせ特徴値の反復回数に伴う変化曲線の増大変化率は速過ぎることになってしまう可能性があり、そうなると、曲線の底点を適切な反復フィルタ処理停止時間点よりも早く出現させてしまうことになってしまう。それゆえ、異なる類型（又は特性）の画像に対して分析を実施する時には、異なるｃの値によって導かれる異なる変換分散曲線に基づいて最適なｃ値を確定し、変換分散曲線と相関性曲線の和曲線に底点を出現させるようにし、且つ底点の位置を反復停止時間点が適切になるようにちょうど良く対応するようにする。
【００９０】
実際の適用においては、各類型（特性）毎の画像に対して予め設定された対応するｃ値を選択する。各類型（例えば、各種画像生成装置又はスキャンプロトコル）毎の画像は、同じｃ値を対応させることができ、異なる類型の画像は異なるｃ値を対応させることができる。具体的例として、複数類型の画像に対して、複数の異なるｃ0値を予め決め、これらの値（例えば、プロファイル形式）をメモリ部に記憶させる（該メモリ部は、画像処理装置に内在するメモリ装置であっても良いし、画像処理装置外部のメモリ装置であって該画像処理装置を経由して送るようにしても良い）。このように、実際の画像処理過程においては、人が事前になにかを行うといった必要はなく、複数類型の画像に対して処理を実施し、画像処理の自動化レベルを大きく向上させることができる。
【００９１】
その他の具体的な一例として、組み合わせ部１３０３−５は、以下の式により、画像特徴計算部１３０３−４にて計算して得られた分散に対して変換処理を実施する。
【００９２】
【数１０】
【００９３】
ここで、Ｖiは画像特徴計算部１３０３−４にて計算して得られたノイズ画像の分散を表し、Ｖi’は、組み合わせ部１３０３−５にて分散値Ｖiに対して変換処理を実施した後に得られる変換値を表し、αは、予め定められた変換パラメータを表す。パラメータαは、実際の応用場面に基づいて予め決められた値であり、例えば、異なる類型（又は特性）の画像に対して分析を実施して得られる実験値又は経験値であって、異なる類型（又は特性）の画像に対して分析を実施する時には、異なるα値により導かれる異なる変換分散曲線に基づいて、最適なα値を確定し、変換分散曲線と相関性曲線の和曲線に底点を出現させるようにし、且つ底点の位置を反復停止時間点が適切になるようにちょうど良く対応するようにする。もちろん、これに限定されるものではない。実際の適用においては、各類型（特性）毎の画像に対して予め設定された対応するα値を選択する。各類型毎の画像（例えば、各種画像生成装置で生成された画像又は各種スキャンプロトコルに基づいて生成された画像）は、同じα値を対応させることができ、異なる類型の画像は異なるα値を対応させることができる。
【００９４】
組み合わせ部１３０３−５は、得られた組み合わせ特徴値を停止判断部１３０３−３へ出力する。停止判断部１３０３−３は、前記組み合わせ特徴値が反復フィルタ処理過程において極小値（底点）に達した時に、フィルタ装置１３０１に対して反復フィルタ処理を停止する指示をする（図１４のステップ１４０６−３）ものである。具体的には、停止判断部１３０３−３は、組み合わせ特徴値が極小値に達したか否かを判断し、フィルタ装置１３０１に対して反復フィルタ処理の停止を指示するか、又は、フィルタ装置１３０１に対して次の反復を継続して実施するよう指示をするものである。停止判断部１３０３−３は、適した方法を用いて組み合わせ特徴値が底点（極小値）に達したか否かを判断することができ、例えば、直前のフィルタ処理で得られた組み合わせ特徴値とその前のフィルタ処理で得られた組み合わせ特徴値を比較し、もし前者が後者に等しいか又は大きい場合、組み合わせ特徴値は増大する傾向を示し、それゆえ、反復フィルタ処理を停止することができる。当然ながら、停止判断部１３０３−３は、その他の方法を用いて組み合わせ特徴値が底点に達したか否かを判断しても良い。ここでは、いちいち列挙して詳述はしない。
【００９５】
図１３〜図１５に示す画像処理装置又は方法において、フィルタリング画像とノイズ画像との相関係数及びノイズ画像の予め定められた画像特徴値の組み合わせ特徴値を利用することにより、反復フィルタ処理を停止する時間を有効に確定することができ、画像における大部分のノイズを除去できると共に、画像における有用な情報を平滑化で失うこともない。
【００９６】
本実施形態の一つの実施例として、非線形拡散フィルタ法を用いて反復フィルタ処理を実施する画像処理装置又は方法を提供することができる。
【００９７】
Ｐ．ＰｅｒｏｎａとＪ．Ｍａｌｉｋの論文“Scale-Space and Edge Detection Using Anisotropic Diffusion”（Proceedings of IEEE，Computer Society Workshop on Computer Vision刊行，１９８７年１１月，第１６−２２頁）（以下、関連文献２と称す。）には、非線形拡散フィルタ法（ここで、該関連文献２では、その開示している方法を異方性拡散（ＡＴＤ）と称しているが、該方法の本質は非線形拡散であるということであって、正真正銘の異方性拡散であるということではない。にもかかわらず、本分野の論文やその他の文献によっては、依然として関連文献２の方法を異方性拡散と呼ぶものもある）が開示され、その中に以下の等式（Ｐ-Ｍ等式と呼ぶ）が開示されている。
【００９８】
【数１１】
【００９９】
上記の式では、Ｉは画像を表し、ｇ拡散係数を表す。拡散係数gは画像勾配｜∇Ｉ｜の関数であって、以下のように表すことができる。
【０１００】
【数１２】
【０１０１】
式（１１）におけるＫは勾配閾値と称す。勾配閾値Ｋは画像における境界の感度を表すパラメータである。
【０１０２】
図１６、図１７に示す実施例は、画像の勾配分布を利用して非線形拡散フィルタにおいて用いる勾配閾値Ｋを確定する画像処理装置及びその方法を提供するものである。図１６は該実施例に基づく画像処理装置の構成を示すブロック図であり、図１７は該実施例に基づく画像処理方法を示すフロー図である。
【０１０３】
図１６に示すとおり、画像処理装置１６００は、フィルタ装置１６０１、勾配閾値確定装置１６０５、及び反復フィルタ処理停止装置１６０３を備える。画像処理装置１６００は、図１７に示す画像処理方法を用いて画像フィルタ処理を実施することができる。
【０１０４】
フィルタ装置１６０１は、非線形拡散フィルタ法に基づき、入力画像に対して反復フィルタ処理を実施し（図１７のステップ１７０２）、且つ、各反復を実施後、得られるフィルタリング画像を反復フィルタ処理停止装置１６０３へ出力する。フィルタ装置１６０１は、いかなる非線形拡散フィルタ計算方法を用いてフィルタ処理を実施しても良く、ここでは、これに限定されるものではない。
【０１０５】
反復フィルタ処理停止装置１６０３は、フィルタリング画像の反復フィルタ処理過程における反復回数に伴って増大変化する変化率に基づいて、反復フィルタ処理を停止するか否かを確定する（図１７のステップ１７０６）。具体的には、反復フィルタ処理停止装置１６０３は、上述実施例／例示における反復フィルタ処理停止装置１０３、３０３、８０３、９０３又は１３０３と同様な機能と構造（例えば、上述の図２〜図１５に示す方法を参考にして反復フィルタ処理停止時間を確定）を備える。ここでは、重複する説明は省略する。
【０１０６】
勾配閾値確定装置１６０５は、フィルタ装置１６０１にて、各フィルタ処理前に、入力画像又は一回前のフィルタ処理で得られたフィルタリング画像における各画像単位の勾配値（図１７のステップ１７１０）を計算するためのものであり、画像中のあらゆる画像単位の勾配分布に基づいて勾配閾値Ｋを確定する（図１７のステップ１７１２）。ここで、“画像単位の勾配値”とは、その画像単位の隣接領域に対する階調（Gray scale）変化値をいう。
【０１０７】
ここでの画像単位とは、画像の画素であり、例えば、２次元画像の画素（pixel）又は３次元画像の画素（voxel）等である。図１８は、画像の各画像単位の勾配値分布を示す一例である。図１８では、横軸は画像単位の勾配値を表し、縦軸は、勾配値に対応する画像単位の個数を表す。一例として、勾配閾値確定装置１６０５は、入力画像に対応する比例値を取得し、該比例値に従って、勾配分布に基づいて勾配閾値Kを計算し、画像における勾配値が勾配閾値Ｋよりも小さい画像単位の個数とあらゆる画像単位の個数との比が前記比例値に近づくようにする（当然ながら、該比例に対してその他の数学的変化をさせることは可能で、例えば、その反映画像の勾配値が勾配閾値Ｋよりも小さい画像単位の個数と勾配値が勾配閾値Ｋよりも小さくない画像単位の個数の比であっても良い）。
【０１０８】
一具体的例として、画像単位の勾配値を、小さいほうから大きいほう、又は大きいほうから小さいほうへ順番に並べ、その並べたものの中で前記比例値に対応する勾配値を勾配閾値とし、勾配値が勾配閾値Ｋよりも小さい画像単位の個数とすべての画像単位の個数との比が前記比例値へ近づくようにする。他の具体的例として、画像の画像単位の勾配のヒストグラム（Histogram）を計算し、ヒストグラムに基づいて勾配閾値を決め、勾配値が勾配閾値Ｋよりも小さい画像単位の個数とすべての画像単位の個数との比が前記比例値へ近づくようにする。
【０１０９】
ここでの比例値とは、入力画像の類型又は特性（例えば、医学画像の画像生成装置又はスキャンプロトコル等）に基づいて予め決められた値であり、Ｐ（０＜Ｐ＜１）を使って表すことができる。各種類型の画像は同じ比例値に対応可能で、異なる画像類型は異なる比例値に対応可能である。ここで、実際の応用場面では、各種画像類型に対して、複数の該類型の画像を分析して予め相応の比例値を確定することにより、後続の非線形拡散フィルタ処理において使用できる。ここでは詳細な説明は省略する。例えば、予め決めた複数の画像類型に対応する複数の比例値（例えば、プロファイル形式）をメモリ部（該メモリ部は、画像処理装置に内在するメモリ装置であっても良いし、画像処理装置外部のメモリ装置であって該画像処理装置を経由して送るようにしても良い）へ保存する。実際の非線形拡散フィルタ処理において、勾配閾値確定装置１６０５は、入力画像の類型又は特性（例えば、画像生成装置の種類又はスキャンプロトコル等）に対応づけて保存している相応の比例値を取得する。こうすることで、人が手動で該パラメータを設定することなく、更に画像処理の自動化レベルを向上させることができる。異なる応用場面では該比例値の大小は異なるので、ここでも、該比例値の数字は具体的に限定されない。
【０１１０】
勾配閾値確定装置１６０５は、いかなる適切な方法を用いて画像の画像単位の勾配値を計算し、画像における画像単位の勾配分布に基づいて勾配閾値を確定しても良く、一例として、例えば、画像をｎ次元画像（ｎ（n≧１）とすれば、以下の式に基づいて画像単位の勾配値を計算できる。
【０１１１】
【数１３】
【０１１２】
ここで、Ｄ1，Ｄ2，・・・，Ｄnは、画像の次元数（ｎ次元）を表し、｜Ｉgrad（Ｄ1，Ｄ2，・・・，Ｄn）｜は、画像上の位置（Ｄ1，Ｄ2，・・・，Ｄn）における画素の勾配値を表し、さらに、ＩＤk（１≦ｋ≦ｎ）は、その画素の第Ｄk次元上の導関数を表す。
【０１１３】
各反復において、勾配閾値確定装置１６０５は、計算して得られる勾配閾値をフィルタ装置１６０１へ出力し、フィルタ装置１６０１は、勾配閾値確定装置１６０５にて決めた勾配閾値に基づいて非線形拡散フィルタ処理を実施する（図１７のステップ１７０２）。
【０１１４】
本実施形態の実施例に基づく画像処理方法と装置は、各種画像のフィルタ処理に応用することができる。これらの画像は１次元又は多次元（例えば、２次元又は３次元等）であって、例えば、本実施形態の実施例に基づく画像処理方法と装置は、医学画像（例えば、医療画像診断装置により取得されたデータに基づいて生成された画像）のノイズ軽減処理に応用可能である。一例として、上述方法の各ステップ及び上述装置の各構成及び／又は部分は、医療画像診断装置（例えば、Ｘ線診断装置、超音波診断装置、ＣＴ装置、磁気共鳴診断装置、ＰＥＴ装置等）のソフトウエア、ファームウエア、ハードウエア又はそれらの組み合わせとして、また、該医療画像診断装置の中の一部分として実施することができる。一例として、既存の医療画像診断装置において、本実施形態の上述方法及び／又は装置に基づいて実施することができ、しかも、既存の医療画像診断装置の各構成部分に対してある程度改良をするだけで即実施可能となる。その他の例として、上述方法の各ステップ及び上述装置の各構成及び／又は部分は、前記医療画像診断装置とは独立した装置として実施してもよい。上述装置の各構成や部分は、ソフトウエア、ファームウエア、ハードウエア又はそれらの組み合わせた方式により構成することができるが、使用する具体的手段や方式は当業者の周知のものであり、ここでは特に述べることはしない。
【０１１５】
一例として、上述方法の各ステップ及び上述装置の各構成及び／又は部分はソフトウエア、ファームウエア、ハードウエア又はそれらの組み合わせとして実施しても良い。ソフトウエア又はファームウエアを介して実現した場合、上述方法のソフトウエアプログラムを実施するため、メモリ媒体から又はネットワークを介して専用のハードウエア構造のコンピュータ（例えば、図１９に示す汎用コンピュータ１９００）へダウンロードして構成することができ、該コンピュータに各種プログラムがダウンロードされた状態で、各種機能等を実施することができる。
【０１１６】
図１９において、演算処理部（即ち、ＣＰＵ）１９０１は、読み出し専用メモリ（ＲＯＭ）１９０２の中に記憶されているプログラム、又は、記憶部１９０８から読み書き兼用メモリ（ＲＡＭ）１９０３へ書き込まれたプログラムに基づいて、各種処理を実施する。ＲＡＭ１９０３では、必要に応じて、ＣＰＵ１９０１が各種処理等を実施するときに必要なデータも記憶しておく。ＣＰＵ１９０１、ＲＯＭ１９０２及びＲＡＭ１９０３は、綜合ライン１９０４を経由してそれぞれ接続されている。入力／出力インターフェース１９０５も、綜合ライン１９０４につながっている。
【０１１７】
下記の各部は、入力／出力インターフェース１９０５に接続されている：入力部１９０６（キーボード、マウス等を含む）、出力部１９０７（モニタ、例えば、ブラウン管（ＣＲＴ）、液晶モニタ（ＬＣＤ）等や、スピーカ等を含む）、記憶部１９０８（キーボードを含む）、通信部１９０９（ネットワークインターフェースカード、例えば、ＬＡＮカード、モデム等）。通信部１９０９は、ネットワーク（例えば、インターネット）を介して通信処理を実施する。必要に応じて、駆動器１９１０も入力／出力インターフェース１９０５に接続可能である。取り外し可能な媒体１９１１は、例えば、磁気ディスク、光ディスク、ＭＯ、半導体メモリ等であって、必要に応じて駆動器１９１０に装着され、必要に応じてコンピュータプログラムを読み出して、記憶部１９０８へダウンロードされる。
【０１１８】
ソフトウエアを介して上述システム処理を実施する場合、ネットワーク（例えば、インターネット又は記憶媒体（例えば、取り外し可能な媒体１９１１））からプログラムをダウンロードしても良い。
【０１１９】
当業者においては、このような記憶媒体は図１９に示すようなプログラムを記憶した記憶媒体は、装置とは離れたところからユーザにプログラムを提供する取り外し可能な媒体１９１１に限らない。取り外し可能な媒体１９１１の例としては、磁気ディスク（フロッピー（登録商標）ディスク、光ディスク（ＣＤ−ＲＯＭやＤＶＤを含む）、磁気光ディスク（ＭｉｎｉＤｉｓｃ（ＭＤ、登録商標）を含む）らを含む。また、記憶媒体はＲＯＭ１９０２であっても良く、記憶部１９０８に含まれるハードディスク等、その中にプログラムが記憶され、それらを含む装置からユーザへプログラムが送られる形態でも良い。
【０１２０】
本実施形態では、更に、メモリとして、機器読み取り可能なコマンドコードを記憶しているプログラム製品でも応用でき、前記コマンドコードが機器を介して読み取られると、本実施形態の実施例における方法が実施される。
【０１２１】
上述機器読み取り可能なコマンドコードを記憶しているプログラム製品を受け入れるための記憶媒体も本実施形態に適用できる。その記憶媒体は、ハードディスク、光ディスク、磁気光ディスク、メモリカード、メモリスティックには限定されない。
【０１２２】
上記の本実施形態の具体的実施例においては、一つの実施方式に示す特徴について、同様の方法を一つ又は複数のほかの実施方法方の中で適用したり、その他の実施方法と組み合わせたり、又はその他の実施方法における特徴に替えるといったことも可能である。
【０１２３】
さらに、“包含する／含む”といった用語を使用したときは、特徴・構成・ステップ又は構造の存在を指し示す。ただし、その他の特徴・構成・ステップ又は構造の存在や付加の排除を意味するものではない。
【０１２４】
上記実施例においては、数字構成の図番記号を用いて各ステップや構成を表記している。ただし、これらの図番記号は単なる説明や画図の都合への考慮によるものであって、その順序やいかなるほかの限定を表すものではない、と当業者は理解すべきである。
【０１２５】
このほか、本実施形態の方法は、詳細な説明の欄において説明された時間順序に沿って実施されるものに限らず、その他の時間順序に沿って、同時に、又は独立して実施されても良い。それゆえ、本願の詳細な説明において説明された方法の実施順序は、本実施形態の技術範囲に対する構成を制限するものではない。
【０１２６】
上記では、既に、本実施形態の具体的実施例の説明をもって、本実施形態の説明を行っているものの、上述のすべての実施例はすべて単なる例示に過ぎず、限定するものではない。当業者は、特許請求の主旨や範囲において、本実施形態の各種手直し・改良又は同等物の設計を行うことが可能である。これらの手直し・改良又は同等物は、本実施形態の保護範囲内に含まれるものである。
【符号の説明】
【０１２７】
１００ 画像処理装置
１０１ フィルタ装置
１０３ 反復フィルタ処理停止装置
【技術分野】
【０００１】
本実施形態は画像処理分野に関し、具体的には、画像（例えば、医学画像）に対してノイズを低減するための画像処理装置及びその方法に関する。
【背景技術】
【０００２】
画像のノイズ低減は、画像処理における重要なアイテムの一つである。例えば、医学画像処理において、装置の信号検出性能や画像生成アルゴリズムが適合していない、又は、造影剤量が少ない等の諸原因によって、ＣＴ（Computed Tomography）装置やＭＲＩ（Magnetic Resonance Imaging）装置・超音波（Ultrasound）診断装置等の医療診断装置により取得された医学画像にノイズが含まれてしまうことがある。これら画像におけるノイズは、医者の診断に多大な影響を及ぼし、後に続く画像処理（例えば、画像分割、計測、ボリュームレンダリング（Volume Rending）等）の正確性にも影響を及ぼしてしまう。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】Ｐ．Ｐｅｒｏｎａ，Ｊ．Ｍａｌｉｋ，“Scale-Space and Edge Detection Using Anisotropic Diffusion”，Proceedings of IEEE，Computer Society Workshop on Computer Vision刊行，１９８７年１１月，第１６−２２頁
【発明の概要】
【発明が解決しようとする課題】
【０００４】
本発明が解決しようとする課題は、画像に対する反復フィルタ処理の過程において、反復フィルタ処理を適切なタイミングで停止することができる画像処理装置、画像処理方法、及び画像処理プログラムを提供することである。
【課題を解決するための手段】
【０００５】
実施形態に係る画像処理装置は、フィルタ装置と、反復フィルタ処理停止装置とを備える。フィルタ装置は、入力画像に対して反復フィルタ処理を実施する。反復フィルタ処理停止装置は、前記反復フィルタ処理の過程における各反復毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する。
【図面の簡単な説明】
【０００６】
【図１】図１は、一実施形態に係る画像処理装置の構成を示すブロック図である。
【図２】図２は、上記実施形態に係る画像処理方法を示すフロー図である。
【図３】図３は、本実施形態の一つの具体的実施例に係る画像処理装置の構成を示すブロック図である。
【図４】図４は、上記具体的実施例に係る画像処理方法を示すフロー図である。
【図５】図５は、反復フィルタ処理を停止する時間を確定するための一つの具体的例を示すフロー図である。
【図６】図６は、それぞれ異なる類型の画像に対して反復フィルタ処理過程におけるフィルタリング画像とノイズ画像の関係係数の変化曲線を示す図である。
【図７】図７は、本実施形態の実施例に基づく応用場面の一つの例を示す図である。
【図８】図８は、本実施形態の他の一具体的実施例に基づく画像処理装置の構成を示すブロック図である。
【図９】図９は、本実施形態の他の一具体的実施例に基づく画像処理装置の構成を示すブロック図である。
【図１０】図１０は、図８に示す実施例の画像処理方法を示すフロー図である。
【図１１】図１１は、図９に示す実施例の画像処理方法を示すフロー図である。
【図１２】図１２は、それぞれ反復フィルタ処理過程におけるフィルタリング画像とノイズ画像の分散の変化曲線を示す図である。
【図１３】図１３は、本実施形態の他の一具体的実施例に基づく画像処理装置の構成を示すブロック図である。
【図１４】図１４は、上記具体的実施例に基づく画像処理方法を示すフロー図である。
【図１５】図１５は、上記具体的実施例の一つの具体的応用例を示す図である。
【図１６】図１６は、本実施形態の他の一実施例に基づく画像処理装置の構成を示すブロック図である。
【図１７】図１７は、上記実施例に基づく画像処理方法を示すフロー図である。
【図１８】図１８は、画像における各画像単位の勾配値分布の例を示す図である。
【図１９】図１９は、本実施形態を実現するためのコンピュータの構成を示すブロック図である。
【発明を実施するための形態】
【０００７】
以下において、本実施形態の基本的理解のために、その一例の概要について説明する。ただし、この概要は、本実施形態のキーワードや重要な部分を確定するものではなく、本実施形態の範囲を限定するものでもない。その目的は、単に、本実施形態の概念を簡潔に説明するためのものであり、その後に論述する更に詳細な説明の序章にすぎない。
【０００８】
本実施形態に係る画像処理装置は、入力画像に対して反復フィルタ（iterative filter）処理を実施するフィルタ装置と、前記反復フィルタ処理の過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する反復フィルタ停止装置と、を備えたことを特徴とする。
【０００９】
本実施形態に係る画像処理方法は、入力画像に対して反復フィルタ処理を実施するステップと、前記反復フィルタ処理過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定するステップを含んだことを特徴とする。
【００１０】
また、本実施形態は、更に上述方法を実現するためのコンピュータプログラムを提供する。
【００１１】
更に、本実施形態は、更に少なくともコンピュータ読み取り可能な媒体形式のコンピュータプログラム製品を提供し、その上、上述方法を実現するためのコンピュータプログラムコードを提供する。
【００１２】
以下、本実施形態の構成や目的・特徴・優位性についての理解を助けるため、図面を参照しながら説明する。図面中の各構成は単に本実施形態の原理を示すためのものである。図面において、同じ又は類似する技術特徴や各構成については、同じ又は類似する附番や表記を用いて説明する。
【００１３】
以下、本実施形態について図面を参照しながら説明する。本実施形態の一つの図面又は一つの実施例において説明される構成や特徴は、一つ又は複数の他の図面又は実施例における構成や特徴と組み合わせることが可能である。なお、説明をわかりやすくするために、図面と説明において本実施形態と関係の無い内容や、当業者にとって既知の構成や処理に関する表示や説明については省略する。
【００１４】
画像処理を実施している時、通常は画像に対して反復フィルタ処理を実施して、できるだけ画像中のノイズを除去することができる。もし反復の回数が少なすぎると、ノイズ低減の効果は十分期待できるものにはならず、かと言って、もし反復の回数が多すぎると、画像における有用な情報まで平滑化されてしまい、画像が不鮮明になって、その後の処理に影響を及ぼす。そこで、本実施形態では、画像の反復フィルタ処理の過程において、反復フィルタ処理の適切なタイミングでの停止を実現することができる画像処理装置及びその方法を提供する。
【００１５】
ここで、本実施形態に基づく装置及び方法は、各類型の画像に対するノイズ低減処理に適用することが可能である。例えば、本実施形態に基づく装置及び方法は、医学画像に対するノイズ低減処理に適用することが可能である。ここでの医学画像は、医療画像診断装置で取得した被験者のデータに基づいて生成された画像である。ここでの医療画像診断装置は、Ｘ線画像診断装置、超音波画像診断装置、ＣＴ装置、ＭＲＩ装置、ＰＥＴ（Positron Emission Tomography）装置等に限定されない。図７は、医療画像診断装置を利用して医学画像の取得を示す図である。図７に示すとおり、各種画像生成装置（例えば、上記医療画像診断装置）に予め定められた或いはユーザが定義したスキャンプロトコルを持たせる。ＣＴ装置を例にして説明すると、例えば、心臓同期スキャン（cardiac synchronization）、腹部４相強調スキャン（abdominal 4 phases enhance）、ダイナミックスキャン等のプロトコルを持たせることが可能であり、図７に示すプロトコル１、プロトコル２、・・・・・、プロトコルＮ（Ｎ≧１）のように、いくつかのプロトコルを利用して、複数の画像シーケンス（例えば、図７に示すシーケンス１、シーケンス２、・・・・・、シーケンスＮ等）を得ることができる。それゆえ、処理を必要とする画像は、医療画像診断装置により対応するプロトコルに基づいて患者データを取得して生成された医学画像であって良い。
【００１６】
図１は、本実施形態の一つの実施例に基づく画像処理装置の構成を示すブロック図であり、図２は、この実施例に基づく画像処理方法を示すフロー図である。
【００１７】
図１に示すとおり、画像処理装置１００は、フィルタ装置１０１と反復フィルタ処理停止装置１０３を具備することが可能である。画像処理装置１００は、図２に示す方法を用いて画像フィルタ処理を実施することができる。以下、図２に示す方法フローを参考にしながら、画像処理装置１００における各構成の機能について説明する。
【００１８】
フィルタ装置１０１は、該画像処理装置に入力された入力画像に対して反復フィルタ処理を実施する（図２のステップ２０２）。実際の応用場面に基づいて、フィルタ装置１０１（即ち、ステップ２０２）は、いずれかの適切なフィルタリング方法、例えば、非線形拡散（Non-Linear Diffusion：ＮＬＤ）フィルタリング方法、線形フィルタリング方法、異方性拡散（Anisotropic Diffusion：ＡＴＤ）フィルタリング方法等を採用することができる。なお、ここではこれらに限定しない。
【００１９】
フィルタ装置１０１は、各フィルタ処理後、得られたフィルタリング画像を反復フィルタ処理停止装置１０３へ出力し、反復フィルタ処理停止装置１０３は、各フィルタ処理後に得られたフィルタリング画像に基づいて反復フィルタ処理を停止するか否かを判断する。具体的には、反復フィルタ処理停止装置１０３は、反復フィルタ処理過程におけるフィルタリング画像の反復回数に伴う増大変化の変化率に基づいて反復フィルタ処理を停止するか否かを判断することができる（図２のステップ２０６）。例えば、もし反復フィルタ処理停止装置１０３が、フィルタリング画像の変化速度が予め定められた閾値に等しいか或いはそれより小さいと判断した場合、これは、この時のフィルタ処理操作の画像に対する影響は小さくなり、つまり軽減されたノイズの程度は小さくなったと言える。そうなれば、フィルタ装置１０１に対して反復フィルタ処理を停止するように指示をすることが可能である。そうでなければ、フィルタ装置１０１に対して継続して次の反復を実施するように指示をする。
【００２０】
ここでのフィルタリング画像の“変化”は、直前のフィルタ処理で得られるフィルタリング画像とその前のフィルタ処理（“その前”とは、直前の反復の前に実施された、ある一回の反復フィルタリング処理（例えば、その一回前の反復等）を意味する。なおこれに限定されない）で得られたフィルタリング画像との間の差を意味する。画像に対する反復フィルタ処理過程において、反復回数の増大に伴って、各フィルタ処理毎に得られるフィルタリング画像とその前のフィルタ処理にて得られたフィルタリング画像と比較した変化（即ち、差）が反復回数に伴って小さくなる。つまり、反復回数の増大に伴って、フィルタリング画像の変化速度が緩慢になる。図１と図２に示す装置と方法はこの原理に基づいており、フィルタリング画像の変化速度に基づいて反復フィルタ処理を停止する時間を確定する。このような方法と装置を利用して、反復フィルタ処理を停止する適切な時間を有効に確定することができ、画像における大部分のノイズを除去できると共に、画像における有用な情報が不鮮明になってしまわないことも保証することがきる。さらに、人の手により反復フィルタ処理を停止する時間を設定するといったような煩雑さも回避することができ、画像処理の更なる自動化を実現することができる。
【００２１】
フィルタリング画像の反復過程における変化は、その異なる特徴や特性にまで反映することができ、例えば、各反復毎に得られるフィルタリング画像とノイズ画像との相関性の反復過程における変化、フィルタリング画像又はノイズ画像のある画像特徴（例えば、分散又は標準差）の値の反復過程における変化等にも反映できる。なお、これらに限らず、例えば、上記の各種変化の数学的組み合わせでも良い。以下、図３〜図１５を参照しながら、反復フィルタ処理過程におけるフィルタリング画像の変化速度に基づいて反復フィルタ処理の停止をするか否かの判断をする画像処理装置及びその方法について、更に具体的な一実施例について説明する。
【００２２】
図３は、一つの具体的実施例に基づく画像処理装置の構成を示すブロック図であり、図４は、その具体的実施例に基づく画像処理方法を示すフロー図である。
【００２３】
上述のとおり、反復フィルタ処理過程において、各フィルタ処理で得られるフィルタリング画像は、その前のフィルタ処理で得られるフィルタ画像と比較して変化が小さくなる、すなわち変化速度が緩やかになる。反復フィルタ処理の開始段階において、各フィルタ処理で軽減されるノイズは比較的多く、反復回数の増大に伴って、画像におけるノイズはどんどん少なくなり、画像自身もだんだん平滑になっていく。フィルタリング画像の変化速度はフィルタリング画像の勾配（gradient）に正比例し、画像自身がだんだん平滑になるのに伴って、画像の勾配の変化が小さくなるので、フィルタリング画像の変化速度は緩やかになる。それに対応して、フィルタリング画像におけるノイズ量の変化に伴って、それとフィルタ処理にて除去されるノイズ画像との相関性も変化が発生する。つまり、前記相関性の変化もフィルタリング画像の変化を反映し、相関性の変化速度もフィルタリング画像の変化速度を反映する。図６の（Ａ）、（Ｂ）及び（Ｃ）は、それぞれこの相関性の反復フィルタ処理過程における変化曲線を示している。図６の（Ａ）、（Ｂ）及び（Ｃ）における曲線は、それぞれ３つの異なる画像（図に示す例では、３つの画像はそれぞれ超音波診断装置（ＵＬ）、Ｘ線ＣＴ装置（ＣＴ）及び磁気共鳴診断装置（ＭＲＩ）で異なる人体部位を撮影して得られた医学画像である。）に対する反復フィルタ処理過程における各フィルタ処理にて得られるフィルタリング画像と相対するノイズ画像の相関係数の時間（反復回数）に伴う変化動向を示している。図６の（Ａ）、（Ｂ）及び（Ｃ）に示すとおり、反復フィルタ処理の開始段階において、フィルタリング画像とノイズ画像の相関係数の値は減少傾向を示し、この傾向はノイズが除去されればされるほど大きくなる。一定時間経過後、相関係数の曲線は底となり（図中の“ｘ”で表記された底点Ｔ１、Ｔ２）、そして底点の後、フィルタリング画像とノイズの相関係数は再び緩やかな増大傾向を示し、この傾向は底点の後、反復回数の増大に伴い、画像中の有用な情報もだんだん平滑化されて失われ、つまり、有用な情報も除去されてしまう、もしくはノイズ画像が入り込み、ひいては、フィルタリング画像とフィルタ処理後のノイズ画像との間の相関性を再び増大させることになる。既存の方法では、上記底点を探し出して、そこを反復フィルタ処理を停止する時間として確定している。この既存の方法では、フィルタリング画像とノイズ画像の相関係数の最小となる時間を反復フィルタ処理を停止する時間点として確定している（例えば、Ｐａｖｅｌ Ｍｒａｚｅｋの論文“Selection of Optimal Stopping Time for Nonlinear Diffusion Filtering”（“Scale-Space and Morphology in Computer Vision”刊行，第三回国際会議，２００１年を参照）（以下、この論文文献を、関連文献１と称す。）。図６の（Ａ）に示す状況においては、底点Ｔ１はただ一点存在しており、且つその底点の位置は正確な停止時間に対応しているので、このような状況においては、上述の既存方法を利用して正確な停止時間を確定することができる。ただし、上記の既存方法は、図６の（Ｂ）又は（Ｃ）に示す状況においては適用できない。なぜなら、例えば、図６の（Ｂ）に示す曲線において、底点Ｔ２が出現している位置は正確な停止時間より明らかに遅れており、且つ底点Ｔ２の後も相関係数曲線はフラットな状態を保持（即ち、底点Ｔ２の後、相関係数値はほとんど変化がない状態）しており、したがって、このような状況においてもし上述の既存方法を利用しようとすれば、確定した時間点を底点Ｔ２としたとしても、それは正確な停止時間ではない。また、例えば、図６の（Ｃ）に示す曲線においては、単調に減少し続け、底点が無いので、このような状況においてたとえ上述の既存方法を利用したとしても、相関係数の最小の時間点が探し出すことができず、正確な停止時間を確定することができない。
【００２４】
本出願の発明者は、図６の（Ａ）、（Ｂ）及び（Ｃ）の３つの曲線に示すとおり、底点の有無に関係なく、この３つの曲線は、反復フィルタ処理の開始段階において、全て単調減少を示し、減少変化率は時間の経過に伴って緩やかになっていくことを見出した。図３と図４に示す装置と方法では、フィルタリング画像とノイズ画像の間の相関性が反復フィルタ処理過程において減少するその変化率（相関性が最小となる時間点ではない）を利用して、反復フィルタ処理を停止する時間を確定する。
【００２５】
図３に示すとおり、画像処理装置３００は、フィルタ装置３０１と反復フィルタ処理停止装置３０３を備え、反復フィルタ処理停止装置３０３は、ノイズ画像取得部３０３−１、相関性取得部３０３−２及び停止判断部３０３−３を備える。画像処理装置３００は、図４に示す方法フローを用いて画像フィルタ処理を実施する。以下、図４に示す方法フロー（流れ）を参考にしながら、画像処理装置３００における各装置及び各部の機能について説明する。
【００２６】
フィルタ装置３０１と図１に示すフィルタ装置１０１は類似するので、いかなる適切な方法を用いて該画像処理装置に入力された入力画像に対して反復フィルタ処理を実施するかについては（図４のステップ４０２）、ここでは説明を省略する。
【００２７】
フィルタ装置３０１は、各フィルタ処理後、得られるフィルタリング画像をノイズ画像取得部３０３−１へ出力する。ノイズ画像取得部３０３−１は、各フィルタ処理後、得られるフィルタリング画像と入力画像に基づいてノイズ画像を推定する（図４のステップ４０６−１）。
【００２８】
一つの例として、各フィルタ処理後、得られるフィルタリング画像と入力画像との間の差分画像を計算し、それをノイズ画像とすることができる。入力画像をＩ0、直前のフィルタ処理後のフィルタリング画像をＩstとした場合、対応するノイズ画像Ｉntは、以下の式で推定できる。
【００２９】
【数１】
【００３０】
なお、各フィルタ処理後、いずれかの適切なその他の方法を用いてノイズ画像を推定することができる（例えば、式（１）に基づいて計算して得られた差分画像の絶対値、即ち、Ｉnt＝｜Ｉ0−Ｉst｜をノイズ画像とすることができる。もちろん他の方法によりノイズ画像を指定しても良いが、ここでは、いちいち列挙して詳述はしない）。
【００３１】
ノイズ画像取得部３０３−１は、得られたノイズ画像を相関性取得部３０３−２へ出力し、相関性取得部３０３−２は、該ノイズ画像と直前の反復にて得られたフィルタリング画像の相関性を推定し（図４のステップ４０６−２）、推定結果を停止判断部３０３−３へ出力する。
【００３２】
一つの例として、両者画像間の相関性を表徴するため、ノイズ画像とフィルタリング画像の間の相関係数を計算する。例えば、両者画像間の相関係数は、以下の式を用いて計算することができる。
【００３３】
【数２】
【００３４】
ここで、ρtは、直前の反復で得られたフィルタリング画像Ｉstと対応するノイズ画像Ｉntとの間の相関係数を表し、Ｅ（・）は、平均値関数（mean value）を表す。
【００３５】
なお、相関係数以外の他のパラメータを用いて両者画像間の相関性を表して、上記の式（２）に対していずれかの形式の数学的変化を求めるようにしても良く、それゆえ、本実施形態は上述の例に限定されることは無い。
【００３６】
停止判断部３０３−３は、各フィルタ処理で得られたフィルタリング画像とノイズ画像との間の相関性が前記反復フィルタ処理過程において減少するその変化率に基づいて反復フィルタ処理を停止するか否かを確定する（図４のステップ４０６−３）。もし、停止判断部３０３−３が、相関性の減少変化率が既に予め定めた閾値に等しいか又はそれより小さいと判断した場合、フィルタ装置３０１に対して反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置３０１に対して次の反復の実施を指示する。
【００３７】
図５は、相関性の減少変化率を表す特徴値に基づいて、反復の停止時間を確定する一例を示している。図５に示すとおり、ステップ４０６−３１において、停止判断部３０３−３は、現在のフィルタ処理で得られたフィルタリング画像とノイズ画像との間の相関性の値及びその前のフィルタ処理後に推定された相関性の値に基づいて、相関性が反復回数に伴って増大又は減少する変化率を表すのに用いる特徴値を計算する。その後、ステップ４０６−３２において、停止判断部３０３−３は、該特徴値と予め定められた閾値との間で予め定めた条件を満足するか否かを判断し、反復フィルタ処理を停止するか、そうでなければ、次の反復を継続して実施するようにする。
【００３８】
一つの具体的例として、前記特徴値は、フィルタリング画像とノイズ画像の相関性の反復回数に伴う変化曲線（図６の（Ａ）、（Ｂ）又は（Ｃ）に示す相関係数曲線）の傾斜角度である。例えば、前記傾斜角度θは、以下の式を用いて計算することができる。
【００３９】
【数３】
【００４０】
ここで、ρtは、時間ｔにおいて得られたフィルタリング画像Ｉstと対応するノイズ画像Ｉntの間の相関係数を表し、ρt-Δtは、時間ｔ−Δｔにおいて得られたフィルタリング画像Ｉs(t-Δt)と対応するノイズ画像Ｉn(t-Δt)の間の相関係数を表し、Δｔは、時間間隔を表す。
【００４１】
その他の一具体的例として、前記傾斜角度は弧度（ラジアン）を用いても表現できる。
【００４２】
【数４】
【００４３】
停止判断部３０３−３が、計算して得られた傾斜角度θrad又はθが、予め定められた閾値に等しい或いはそれより小さい場合、フィルタ装置３０１に対して、反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置３０１に対して、次の反復の実施を指示する。
【００４４】
その他の具体例として、前記特徴値は、フィルタリング画像とノイズ画像の関連性の反復回数に伴う変化曲線（図６の（Ａ）、（Ｂ）又は（Ｃ）に示す相関係数曲線）の傾斜率であっても良い。例えば、前記傾斜率ｋslopeは、以下の式を用いて計算することができる。
【００４５】
【数５】
【００４６】
ここで、ρtは、時間ｔにおいて得られるフィルタリング画像Ｉstと対応するノイズ画像Ｉntとの間の相関係数を表し、ρt-Δtは、時間ｔ−Δｔにおいて得られるフィルタリング画像Ｉs(t-Δt)と対応するノイズ画像Ｉn(t-Δt)の間の相関係数を表し、Δｔは、時間間隔を表す。
【００４７】
停止判断部３０３−３が、計算して得られた傾斜率ｋslopeが、予め定められた閾値と等しい又はそれより小さい場合、フィルタ装置３０１に対して反復フィルタ処理を停止するように指示する。そうでなければ、フィルタ装置３０１に対して次の反復を実施するように指示する。
【００４８】
ここで、フィルタリング画像とノイズ画像との間の相関性の前記反復フィルタ処理過程における減少変化率を表す特徴値としては、上述の相関係数曲線の傾斜角度や傾斜率に限られない。前記特徴値として、相関性の減少変化率を表す他の特徴を用いても構わない。例えば、直前のフィルタ処理後に推定した相関係数とその前のフィルタ処理後に推定した相関係数ρt-Δtとの間の差（又は、該差の逆数）を前記特徴値として用いてもかまわない。上述の差が予め定められた閾値と等しい又はそれより小さい場合（又は、上述の差の逆数が予め定められた閾値に等しい又はそれより大きい場合）、反復フィルタ処理を停止する、そうでなければ、継続して次の反復を実施する。この他、上述の各種特徴値の計算方式に対して、いかなる形式の数学的変化を求めるようにしても良い。
【００４９】
その他、ここで、上述の実施例／例示において前記の予め定められた閾値は、実際の応用場面に基づいて予め確定されたものであり、異なる類型（又は特性）に対する画像により分析して得られた実験値や経験値であっても良い。異なる応用場面や採用する異なる特徴値パラメータに従って、該閾値は異なる値を採用することができる。例えば、前記入力画像が医療画像診断装置を介して得られたデータに基づいて生成された医学画像である場合、実験又は経験に基づいて、各種スキャンプロトコル又は画像生成装置の類型に従って予め対応する閾値を確定することができる。又は、複数のスキャンプロトコル又は画像生成装置に従う複数の閾値（例えば、プロファイル等の形式でもって）をメモリ部（該メモリ部は、画像処理装置に内在するメモリ装置であっても良いし、画像処理装置外部のメモリ装置であって該画像処理装置を経由して送るようにしても良い）へ保存する。このように、反復フィルタ処理過程において、入力画像のスキャンプロトコル又は画像生成装置の類型に基づいて対応する閾値を取得して（例えば、メモリ部から取り込む）、相関性の減少変化率を判断する根拠とすることができるので、更に画像処理の自動化の程度を向上させることができる。
【００５０】
図３〜５に示す画像処理装置又はその方法を参照する中において、フィルタリング画像とノイズ画像の相関性の減少変化率を利用することにより、反復フィルタ処理を停止する時間を有効に確定することができ、画像における大部分のノイズを除去できると共に、画像における有用な情報を平滑化により失われることがないようにすることも可能である。
【００５１】
図８及び図９は、それぞれ他の２つの具体的実施例に基づく画像処理装置の構造を示すブロック図であり、図１０及び図１１は、それぞれこの２つの具体的実施例に基づく画像処理方法を示すフロー図である。図３〜図５に示す実施例では、フィルタリング画像とノイズ画像の相関性の減少する変化率を利用して反復フィルタ処理過程におけるフィルタリング画像の変化速度を表す。図８〜図１１に示す実施例では、フィルタリング画像又はノイズ画像の予め定められた画像特徴の反復フィルタ処理過程における変化速度を利用してフィルタリング画像の変化速度を表す。
【００５２】
ここでの予め定められた画像特徴は、反復フィルタ処理過程におけるフィルタリング画像の変化のいかなる画像特徴値をも反映することができ、例えば、フィルタリング画像又はノイズ画像の分散又は標準差、フィルタリング画像又はノイズ画像の画素の平均勾配、フィルタリング画像又はノイズ画像のグレイスケール（gray scale）のダイナミックレンジ（即ち、最大グレイスケールと最小グレイスケールとの差）等であっても良い。ここでは、いちいち列挙して詳述はしない。
【００５３】
図１２の（Ａ）及び（Ｂ）は、それぞれ、反復フィルタ処理過程において、フィルタリング画像とノイズ画像の分散を例にした変化曲線を示している。図１２の（Ａ）に示すとおり、反復回数の増大に伴って、フィルタリング画像の分散は単調に減少し、且つ減少変化率は緩やかになっていく。図１２の（Ｂ）に示すとおり、反復回数の増大に伴って、ノイズ画像の分散は単調に増大し、且つ増大変化率は緩やかになっていく。上記の方法を利用してノイズ画像を推定することができる。説明は重複するので省略する。
【００５４】
以下、図８及び図１０を参照しながら、フィルタリング画像の予め定められた画像特徴の変化速度を利用してフィルタ処理停止時間を確定する実施例について説明する。
【００５５】
図８に示すとおり、画像処理装置８００は、フィルタ装置８０１と反復フィルタ処理停止装置８０３を備え、更に、反復フィルタ処理停止装置８０３は、画像特徴計算部８０３−４と停止判断部８０３−３を備える。画像処理装置８００は、図１０に示す画像処理方法を用いて画像フィルタ処理を実施することができる。
【００５６】
フィルタ装置８０１は、該画像処理装置へ入力された入力画像に対して反復フィルタ処理（図１０のステップ１００２）を実施するためのものであって、上記におけるフィルタ装置１０１及び３０１と類似するものであり、ここでは、詳細な説明は省略する。
【００５７】
フィルタ装置８０１は、各フィルタ処理後、得られたフィルタリング処理を画像特徴計算部８０３−４へ出力し、画像特徴計算部８０３−４を介して、各フィルタ処理で得られたフィルタリング画像の予め定められた画像特徴値を計算する（図１０に示すステップ１０００６−１）。例えば、フィルタリング画像の分散又は標準差を計算し、計算結果を停止判断部８０３−３へ出力する。停止判断部８０３−３は、反復フィルタ処理過程における前記画像特徴（例えば、分散又は標準差）値の反復回数の増加に伴う変化率（例えば、分散又は標準差の減少変化率）に基づいて反復フィルタ処理を停止するか否かを確定する。もし、停止判断部８０３−３が、前記画像特徴値の変化速度が予め定められた閾値と等しい又はそれよりも小さいと判断した場合、フィルタ装置８０１に対して反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置８０１に対して、次の反復の実施を継続する指示をする。
【００５８】
具体的例として、停止判断部８０３−３は、画像特徴計算部８０３−４の計算結果に基づいて、予め定められた画像特徴値の反復フィルタ処理過程における変化速度を表す特徴値を計算する。例えば、該特徴値は、予め定められた画像特徴値が反復回数の増加に伴って変化する曲線の傾斜角度又は傾斜率（例えば、図１２の（Ａ）に示す曲線の傾斜角度又は傾斜率であって、上述の公式（３）−（５）を用いて算出することができる。詳細については省略）、又は、該特徴値は、直前のフィルタ処理で得られたフィルタリング画像の予め定められた画像特徴値とその前のフィルタ処理で得られたフィルタリング画像の予め定められた画像特徴値との間の差（又は、該差の逆数）であって良い。ここでは、いちいち列挙して詳述はしない。その後、停止判断部８０３−３は、計算して得られた特徴値と予め定めた閾値との間で予め定めた関係を満足するか否かを判断し（例えば、傾斜角度、傾斜率又は上述差を特徴値として用いる状況下で、特徴値が予め定められた閾値に等しい又はそれより小さいか否かを判断する、又は、上述差の逆数を特徴値として用いる状況下で、特徴値が予め定められた閾値に等しい又はそれより大きいか否かを判断する）、フィルタ装置８０１に対して反復フィルタ処理の停止を指示する、そうでなければ、フィルタ装置８０１に対して次の反復を継続して実施するように指示をする。上述実施例／例示のように、ここでの閾値は、実際の応用場面に基づいて予め確定されるものである。異なる応用場面又はそこで用いられる異なる特徴パラメータに従って該閾値は異なる値をとる（例えば、図１２の（Ａ）に示す画像特徴値の変化曲線により予め閾値を確定し、これらの曲線は、異なる類型（又は特性）の画像に対して分析することにより取得できる。ここでは、更なる詳細な説明は省略する。）。なお、ここでは、それは、具体的数値に限定されない。
【００５９】
以下、図９と図１１を参照しながら、ノイズ画像の予め定めた画像特徴値の変化速度を利用してフィルタ処理の停止時間を確定する実施例について説明する。
【００６０】
図９に示すとおり、画像処理装置９００は、フィルタ装置９０１と反復フィルタ処理停止装置９０３を備える。反復フィルタ処理停止装置９０３は、図８に示す実施例と同様に、画像特徴計算部９０３−４と停止判断部９０３−３を備え、図８に示す実施例と異なるのは、反復フィルタ処理停止装置９０３は、更にノイズ画像取得部９０３−１を備えることにある。画像処理装置９００は、図１１に示すような画像処理方法を用いて画像フィルタ処理を実施することができる。
【００６１】
フィルタ装置９０１は、上記のフィルタ装置１０１や３０１と同様に、画像処理装置に入力された入力画像に対して反復フィルタ処理を実施するものである（図１１のステップ１１０２）。なお、ここでは重複する部分は説明を省略する。
【００６２】
フィルタ装置９０１は、各フィルタ処理後、得られたフィルタリング画像をノイズ画像取得部９０３−１へ出力する。ノイズ画像取得部９０３−１は、各フィルタ処理後に、得られたフィルタリング画像と入力画像に基づいてノイズ画像を推定する（図１１のステップ１１０６−１）。ノイズ画像取得部９０３−１は、上述の図３や図４に示す方法を用いてノイズ画像を推定する。なお、ここでは重複する部分は説明を省略する。
【００６３】
ノイズ画像取得部９０３−１は、得られたノイズ画像を画像特徴計算部９０３−４へ出力し、画像特徴計算部９０３−４を介してノイズ画像の予め定められた画像特徴値を計算する（図１１に示すステップ１１０６−２）、例えば、ノイズ画像の分散又は標準差を計算し、計算結果を停止判断部９０３―３へ出力する。
【００６４】
停止判断部９０３−３は、反復フィルタ処理過程においてノイズ画像の画像特徴値（例えば、分散や標準差）の反復回数に伴う増大変化の変化率（例えば、分散又は標準差の増大の変化率）に基づいて反復フィルタ処理を停止するか否かを確定する。もし、停止判断部９０３−３が、ノイズ画像の前記画像特徴値の変化速度が予め定めたレベルに達したと判断した場合、フィルタ装置９０１に対して反復フィルタ処理の停止を指示する。そうでなければ、フィルタ装置９０１に対して次の反復の実施を継続する指示をする。
【００６５】
具体的例として、停止判断部９０３−３は、画像特徴計算部９０３−４の計算結果に基づいて、予め定められた画像特徴値の反復フィルタ処理過程における変化速度を表す特徴値を計算する。例えば、該特徴値は、予め定められた画像特徴値の反復回数に伴い増大変化する曲線の傾斜角又は傾斜率を反映するもの（例えば、図１２の（Ｂ）に示す曲線の傾斜角度又は傾斜率であり、上述の式（３）−（５）を用いて計算できる。詳細は重複するので省略）、又は、該特徴値は、直前のフィルタ処理で得られたノイズ画像の予め定められた画像特徴値と、その前のフィルタ処理で得られたノイズ画像の予め定められた画像特徴値との間の差（又は、該差の逆数）である。ここでは、いちいち列挙して詳述はしない。その後、停止判断部９０３−３は、計算して得られた特徴値と予め定められた閾値との間で予め定められた関係を満足するか否かを判断し（例えば、傾斜角、傾斜率又は上述差を特徴値とする場合、特徴値が予め定められた閾値と等しい又はそれより小さいかどうかを判断する。上述差の逆数を特徴値とする場合、特徴値が予め定められた閾値と等しいか又はそれより大きいかどうかを判断する）、フィルタ装置９０１に対して反復フィルタ処理の停止を指示する、又は、フィルタ装置９０１に対して次の反復の実施を継続するよう指示する。上述実施例／例示のように、ここでは、前記の閾値は実際の応用場面に基づいて予め決められたものである。異なる応用場面とそこで用いられる異なる特徴パラメータに従って、該閾値は異なる値（例えば、図１２の（Ｂ）に示す画像特徴値の変化曲線により予め閾値を確定することができる。この曲線は、異なる類型（又は特性）の画像に対して分析を実施することにより得ることができる。）を用いることができる。なお、ここではその具体的数値に限定されない。
【００６６】
図８〜図１１に示す画像処理装置や方法において、反復フィルタ処理過程においてフィルタリング画像又はノイズ画像の予め定められた画像特徴値を利用することにより、反復フィルタ処理を停止する時間を有効に確定することができ、画像中の大部分のノイズを除去できると共に、画像中の有用な情報が平滑化により失うことを防ぐことができる。
【００６７】
具体的例として、図８〜図１１に示す画像処理装置又は方法を参考にして、以下の式により画像の分散を計算し、前記予め定められた画像特徴値とすることができる。
【００６８】
【数６】
【００６９】
ここで、Ｉは画像を表し（例えば、上述のフィルタリング画像やノイズ画像）、Ｄは該画像の分散を表す。
【００７０】
図１３は、他の一具体的実施例に基づく画像処理装置の構造を示すブロック図であり、１４は、該具体的実施例に基づく画像処理方法を示すフロー図である。上記の実施例はそれぞれノイズ画像とフィルタリング画像との間の相関性、又は、ノイズ画像（又はフィルタリング画像）の予め定められた画像特徴値の反復フィルタ処理過程における変化速度を利用し、そして図１３と図１４に示す実施例は両者を組み合わせることを考える。例えば、図６の（Ａ）、（Ｂ）、（Ｃ）に示すとおり、ノイズ画像とフィルタリング画像の間の相関性が反復フィルタ処理の開始段階は単調減少し、そのあと緩やかに増大（例えば、図６の（Ａ）及び（Ｂ）に示すとおり）、又は、緩やかな減少が続き（例えば、図６の（Ｃ）のとおり）、図１２の（Ｂ）に示すとおり、ノイズ画像の分散は反復フィルタ処理過程において単調に増大する。もし、両者が適切な方法により組み合わさる場合（両者が組み合わさった場合の具体的実施例としては、例えば、以下の文中で図１３−１５を参照しながら説明するようなものを挙げることができる。）、生成される組み合わせ特徴値の反復回数に伴う変化曲線も反復フィルタ処理過程において画像の変化を反映する。例えば、図６の（Ｂ）に示すような相関性曲線の底点が時間方向下流へ延びていく状況下では、ノイズ画像の分散は単調増大するので、組み合わせ特徴値の曲線は底点を時間方向上流へ適切な位置まで遡るようシフトする。また、例えば、図６の（Ｃ）に示すような相関性が単調減少する状況下では、ノイズ画像の分散は単調増大するので、組み合わせ特徴値の曲線には底点が出現すると共に該底点は正確な停止時間に対応する。図１５の（Ａ）は、ノイズ画像とフィルタリング画像の相関係数を示したものであり、ノイズ画像の分散及び両者の組み合わせ特徴値の反復回数に伴う変化を示す図である。図１５の（Ａ）において、縦軸は相関係数を表し、分散及び組み合わせ特徴の値を取り、横軸は反復時間を表す。図１５の（Ａ）に示すとおり、組み合わせ特徴値の曲線のｔ＝０．４の時間点にて底点（極小値）が出現し、前記組み合わせ特徴値の反復回数に伴って増大して出現する底点（即ち、極小値）を用いて、反復フィルタ処理を停止する時間を確定することができる。
【００７１】
説明をわかりやすくするために、以下では、ノイズ画像の分散を前記予め定められた画像特徴とし、ノイズ画像とフィルタリング画像の相関係数が前記相関性を表すものとして説明する。なお、その他の例としては、ノイズ画像のその他の画像特徴（例えば、標準差など）、又は、フィルタリング画像の分散や標準差の逆数など、又は、フィルタリング画像とノイズ画像の間の相関性のその他のパラメータでも良い。ここでは、いちいち列挙して詳述はしない。
【００７２】
図１３に示すとおり、画像処理装置１３００は、フィルタ装置１３０１と反復フィルタ処理停止装置１３０３を備え、更に、反復フィルタ処理停止装置１３０３は、ノイズ画像取得部１３０３−１、相関性取得部１３０３−２、画像特徴計算部１３０３−４、組み合わせ部１３０３−５及び停止判断部１３０３−３を備える。画像処理装置１３００は、図１４に示すような画像処理方法を用いて画像フィルタ処理を実施する。
【００７３】
上述のフィルタ装置１０１や３０１のように、フィルタ装置１３０１は、該画像処理装置に入力された入力画像に対して反復フィルタ処理を実施するものである（図１４のステップ１４０２）。ここでは、重複する部分については説明を省略する。
【００７４】
フィルタ装置１３０１は、各フィルタ処理後、得られたフィルタリング画像をノイズ画像取得部１３０３−１に出力する。ノイズ画像取得部１３０３−１は、各フィルタ処理後に得られたフィルタリング画像と入力画像とに基づいてノイズ画像を推定し（図１４のステップ１４０６−１）、ノイズ画像取得部１３０３−１は、上述の図３と図４に示す方法を参考にしてノイズ画像を推定する。ここでは、重複する部分については説明を省略する。
【００７５】
ノイズ画像取得部１３０３−１は、得られたノイズ画像を、それぞれ相関性取得部１３０３−２と画像特徴計算部１３０３−４へ出力する。相関性取得部１３０３−２は、該ノイズ画像と直前のフィルタ処理で得られたフィルタリング画像との相関性を推定し（図１４のステップ１４０６−２）、推定結果を組み合わせ部１３０３−５へ出力する。例えば、相関性取得部１３０３−２は、２つの画像間の相関性を表徴するため、ノイズ画像とフィルタリング画像との間の相関係数を計算する（例えば、上述の方法を用いて計算できる。重複する説明は省略。）。
【００７６】
画像特徴計算部１３０３−４は、ノイズ画像の予め定められた画像特徴値を計算し（例えば、図１４のステップ１４０６−４）、例えば、ノイズ画像の分散を計算し（例えば、上記式（６）を用いて）、計算結果を組み合わせ部１３０３−５へ出力する。
【００７７】
組み合わせ部１３０３−５は、相関性取得部１３０３−２により取得した相関性値と画像特徴計算部１３０３−４により計算して得られた分散とを組み合わせて、組み合わせ特徴値を生成する（例えば、図１４のステップ１４０６−５）。例えば、Ｃ1を単調増大するノイズ画像の分散曲線とし、Ｃ2をノイズ画像とフィルタリング画像の相関曲線とすると、組み合わせ特徴値は次のように表せる。
【００７８】
【数７】
【００７９】
ここで、Ｃcombineは組み合わせ特徴値を表し、ｆ（・）は曲線Ｃ1に対する数学的変化を表し、ｇ（・）は曲線Ｃ2に対する数学的変化を表す。
【００８０】
一例として、組み合わせ部１３０３−５は、相関性取得部１３０３−２により取得した相関性値と画像特徴計算部１３０３−４により計算して得られた分散値を計算して和を求め、該和を組み合わせ特徴値とすることができる。ただし、状況によっては、フィルタリング画像とノイズ画像の相関係数の値を取る範囲は、しばしばノイズ画像の分散の値を取る範囲と異なる。それゆえ、相関性値と分散値の和曲線は、上述の底点（極小値）が存在しない可能性がある。したがって、２つの曲線をマッチングさせるように調整する必要がある。以下、相関性取得部１３０３−２により得られた相関性値と画像特徴計算部１３０３−４により計算して得られた分散値とを組み合わせることについて、その一例を示す。
【００８１】
一例として、フィルタ処理の前に、画像処理装置（例えば、フィルタ装置１３０１）は、先ずオリジナルの入力画像に対してノーマライズ処理を行う。言い換えると、図１４のステップ１４０２の前に、入力画像のノーマライズ処理を行うステップを含む。このように、各反復の後で、ノイズ画像取得部１３０３−１は、ノーマライズされた入力画像を用いて得られた前記相関性値（相関性取得部１３０３−２により得られる）と前記分散値（画像特徴計算部１３０３−４により得られる）の和を求め、該和を組み合わせ特徴値とする。
【００８２】
その他の一例として、画像特徴計算部１３０３−４は、各反復の後、ノイズ画像取得部１３０３−１にて推定されたノイズ画像に対してノーマライズを実施し（該例においては、オリジナル入力画像に対してノーマライズを実施する必要はない）、ノーマライズされたノイズ画像の分散を計算し、ノーマライズされたノイズ画像の分散を組み合わせ部１３０３−５へ出力する。ノイズ画像のダイナミック範囲（例えば、［Ｎmin，Ｎmax］を用いて画像のダイナミック範囲を表すことができる）を［０，１］へノーマライズする。異なる類型の入力画像に対しては、各反復で生成されるノイズ画像のダイナミック範囲は異なり、ノーマライズが一つの同じダイナミック範囲に基づいて実施されない状態を引き起こす。それゆえ、一つの標準ダイナミック範囲を指定し、各反復時に全てその標準ダイナミック範囲でもって［０，１］へのノーマライズを実施する。例えば、ある種類の入力画像に対して、相応する標準ダイナミック範囲を［Ｓmin，Ｓmax］と設定し（ここで、異なる入力画像に対しては、異なる標準ダイナミック範囲を用いることができる）、ノーマライズ処理は、以下の式（８）により実施することができる。
【００８３】
【数８】
ここで、Ｎはノイズ画像の画素値を表し、Ｎ’はノーマライズ処理後の画素値を表す。
【００８４】
ここで、実際の応用場面では、実際の必要性に応じて上述の標準ダイナミックレンジを確定するものであって、本実施形態を上記実施例に限るべきではない。例えば、医学画像の場合、標準ダイナミックレンジは、入力画像を取得した医療画像診断装置の異なる類型に基づいて確定するものである。その他の例では、入力画像のスキャンプロトコル等に基づいて確定しても良い。ここでは、いちいち列挙して詳述はしない。
【００８５】
該例において、相関性取得部１３０３−２は、ノイズ画像取得部１３０３−１により推定されたノイズ画像（ノーマライズされたノイズ画像ではない）に基づいて上記相関性を計算する。
【００８６】
その他の一例として、組み合わせ部１３０３−５は、計算して得られたノイズ画像の分散に対して変換処理を実施し、その後、変換後の分散と相関性取得部１３０３−２にて得られた相関性値の和を求め、該和を組み合わせ特徴値とする。例えば、組み合わせ部１３０３−５は、画像特徴計算部１３０３−４により計算して得られたノイズ画像の分散に対して変換処理を実施し、その後、変換処理後の分散値に相関性取得部１３０３−２により得られた相関性値を加え、得られた和を上記組み合わせ特徴値とする。例えば、組み合わせ部１３０３−５は、以下の式を介して、画像特徴計算部１３０３−４により計算して得られたノイズ画像の分散に対して変換処理を実施する。
【００８７】
【数９】
【００８８】
ここで、Ｖiは画像特徴計算部１３０３−４により計算して得られたノイズ画像の分散を表し、Ｖi’は組み合わせ部１３０３−５により分散値Ｖiに対して変換を実施して得られた変換値を表し、ｉ（ｉ＝１，２，・・・）は反復時間を表し、ｏｆｆｓｅｔは平行移動係数を表し、ｏｆｆｓｅｔ＝√Ｖ1のＶ1は一回目の反復後に画像特徴計算部１３０３−４にて計算して得られたノイズ画像の分散を表し、ａiは反復回数の増大に伴って増大する係数であり、該係数は１より大きく、例えば、ａi＝１．０＋ｃ×ｉ（ｃはゼロより大きい予め定められた定数）である。式（９）において、分散値Ｖiの平方根を求めることにより、変換後の分散値Ｖi’は大きくなり、平行移動係数ｏｆｆｓｅｔにより，変換処理後の分散曲線を平行移動し、ゼロから開始するようにする。係数ａiにより，変換処理後の分散曲線の増大変化率は増大し、反復回数の増大に伴って、該増大変化率も増大する。図１５の（Ｂ）は、上記変換後のノイズ画像の分散、ノイズ画像のフィルタリング画像との相関係数及び両者の和が反復回数に伴って変化する様子を示した図である。図１５の（Ｂ）に示すとおり、相関係数曲線と比べ、得られる和曲線の底点は時間方向上流へ遡るようシフトされる。何回もの検証により、もし、図１５の（Ｂ）に示す相関係数曲線の底点を反復フィルタ処理を停止する時間とすると、反復フィルタ処理を経た出力画像における有用な情報は曖昧となってしまう。もし、和曲線の底点を用いれば、反復フィルタ処理を経た出力画像においても、オリジナルの入力画像における有用な情報は保たれると共に、大部分のノイズは除去される。
【００８９】
その他、上記の式（９）におけるパラメータｃは実際の応用場面に基づいて予め確定した値であり、例えば、異なる類型（又は特性）の画像に対して分析を実施して得られる実験値又は経験値であって良い。上記の式（９）により、ｃ値が小さければ小さいほど、変換後の分散曲線はなだらかになり、反対に、ｃ値が大きければ大きいほど、変換後の分散曲線の増大変化率は速くなることを知ることができる。ただし、ｃ値は小さすぎてはいけない。小さすぎてしまうと、組み合わせ部１３０３−５により得られる組み合わせ特徴値の反復回数に伴う変化曲線は単調減少を引き起こす可能性があり、その場合、極小値が出現しないことになってしまう。かといって、ｃ値は大きすぎてもいけない。大きすぎると、組み合わせ部１３０３−５により得られる組み合わせ特徴値の反復回数に伴う変化曲線の増大変化率は速過ぎることになってしまう可能性があり、そうなると、曲線の底点を適切な反復フィルタ処理停止時間点よりも早く出現させてしまうことになってしまう。それゆえ、異なる類型（又は特性）の画像に対して分析を実施する時には、異なるｃの値によって導かれる異なる変換分散曲線に基づいて最適なｃ値を確定し、変換分散曲線と相関性曲線の和曲線に底点を出現させるようにし、且つ底点の位置を反復停止時間点が適切になるようにちょうど良く対応するようにする。
【００９０】
実際の適用においては、各類型（特性）毎の画像に対して予め設定された対応するｃ値を選択する。各類型（例えば、各種画像生成装置又はスキャンプロトコル）毎の画像は、同じｃ値を対応させることができ、異なる類型の画像は異なるｃ値を対応させることができる。具体的例として、複数類型の画像に対して、複数の異なるｃ0値を予め決め、これらの値（例えば、プロファイル形式）をメモリ部に記憶させる（該メモリ部は、画像処理装置に内在するメモリ装置であっても良いし、画像処理装置外部のメモリ装置であって該画像処理装置を経由して送るようにしても良い）。このように、実際の画像処理過程においては、人が事前になにかを行うといった必要はなく、複数類型の画像に対して処理を実施し、画像処理の自動化レベルを大きく向上させることができる。
【００９１】
その他の具体的な一例として、組み合わせ部１３０３−５は、以下の式により、画像特徴計算部１３０３−４にて計算して得られた分散に対して変換処理を実施する。
【００９２】
【数１０】
【００９３】
ここで、Ｖiは画像特徴計算部１３０３−４にて計算して得られたノイズ画像の分散を表し、Ｖi’は、組み合わせ部１３０３−５にて分散値Ｖiに対して変換処理を実施した後に得られる変換値を表し、αは、予め定められた変換パラメータを表す。パラメータαは、実際の応用場面に基づいて予め決められた値であり、例えば、異なる類型（又は特性）の画像に対して分析を実施して得られる実験値又は経験値であって、異なる類型（又は特性）の画像に対して分析を実施する時には、異なるα値により導かれる異なる変換分散曲線に基づいて、最適なα値を確定し、変換分散曲線と相関性曲線の和曲線に底点を出現させるようにし、且つ底点の位置を反復停止時間点が適切になるようにちょうど良く対応するようにする。もちろん、これに限定されるものではない。実際の適用においては、各類型（特性）毎の画像に対して予め設定された対応するα値を選択する。各類型毎の画像（例えば、各種画像生成装置で生成された画像又は各種スキャンプロトコルに基づいて生成された画像）は、同じα値を対応させることができ、異なる類型の画像は異なるα値を対応させることができる。
【００９４】
組み合わせ部１３０３−５は、得られた組み合わせ特徴値を停止判断部１３０３−３へ出力する。停止判断部１３０３−３は、前記組み合わせ特徴値が反復フィルタ処理過程において極小値（底点）に達した時に、フィルタ装置１３０１に対して反復フィルタ処理を停止する指示をする（図１４のステップ１４０６−３）ものである。具体的には、停止判断部１３０３−３は、組み合わせ特徴値が極小値に達したか否かを判断し、フィルタ装置１３０１に対して反復フィルタ処理の停止を指示するか、又は、フィルタ装置１３０１に対して次の反復を継続して実施するよう指示をするものである。停止判断部１３０３−３は、適した方法を用いて組み合わせ特徴値が底点（極小値）に達したか否かを判断することができ、例えば、直前のフィルタ処理で得られた組み合わせ特徴値とその前のフィルタ処理で得られた組み合わせ特徴値を比較し、もし前者が後者に等しいか又は大きい場合、組み合わせ特徴値は増大する傾向を示し、それゆえ、反復フィルタ処理を停止することができる。当然ながら、停止判断部１３０３−３は、その他の方法を用いて組み合わせ特徴値が底点に達したか否かを判断しても良い。ここでは、いちいち列挙して詳述はしない。
【００９５】
図１３〜図１５に示す画像処理装置又は方法において、フィルタリング画像とノイズ画像との相関係数及びノイズ画像の予め定められた画像特徴値の組み合わせ特徴値を利用することにより、反復フィルタ処理を停止する時間を有効に確定することができ、画像における大部分のノイズを除去できると共に、画像における有用な情報を平滑化で失うこともない。
【００９６】
本実施形態の一つの実施例として、非線形拡散フィルタ法を用いて反復フィルタ処理を実施する画像処理装置又は方法を提供することができる。
【００９７】
Ｐ．ＰｅｒｏｎａとＪ．Ｍａｌｉｋの論文“Scale-Space and Edge Detection Using Anisotropic Diffusion”（Proceedings of IEEE，Computer Society Workshop on Computer Vision刊行，１９８７年１１月，第１６−２２頁）（以下、関連文献２と称す。）には、非線形拡散フィルタ法（ここで、該関連文献２では、その開示している方法を異方性拡散（ＡＴＤ）と称しているが、該方法の本質は非線形拡散であるということであって、正真正銘の異方性拡散であるということではない。にもかかわらず、本分野の論文やその他の文献によっては、依然として関連文献２の方法を異方性拡散と呼ぶものもある）が開示され、その中に以下の等式（Ｐ-Ｍ等式と呼ぶ）が開示されている。
【００９８】
【数１１】
【００９９】
上記の式では、Ｉは画像を表し、ｇ拡散係数を表す。拡散係数gは画像勾配｜∇Ｉ｜の関数であって、以下のように表すことができる。
【０１００】
【数１２】
【０１０１】
式（１１）におけるＫは勾配閾値と称す。勾配閾値Ｋは画像における境界の感度を表すパラメータである。
【０１０２】
図１６、図１７に示す実施例は、画像の勾配分布を利用して非線形拡散フィルタにおいて用いる勾配閾値Ｋを確定する画像処理装置及びその方法を提供するものである。図１６は該実施例に基づく画像処理装置の構成を示すブロック図であり、図１７は該実施例に基づく画像処理方法を示すフロー図である。
【０１０３】
図１６に示すとおり、画像処理装置１６００は、フィルタ装置１６０１、勾配閾値確定装置１６０５、及び反復フィルタ処理停止装置１６０３を備える。画像処理装置１６００は、図１７に示す画像処理方法を用いて画像フィルタ処理を実施することができる。
【０１０４】
フィルタ装置１６０１は、非線形拡散フィルタ法に基づき、入力画像に対して反復フィルタ処理を実施し（図１７のステップ１７０２）、且つ、各反復を実施後、得られるフィルタリング画像を反復フィルタ処理停止装置１６０３へ出力する。フィルタ装置１６０１は、いかなる非線形拡散フィルタ計算方法を用いてフィルタ処理を実施しても良く、ここでは、これに限定されるものではない。
【０１０５】
反復フィルタ処理停止装置１６０３は、フィルタリング画像の反復フィルタ処理過程における反復回数に伴って増大変化する変化率に基づいて、反復フィルタ処理を停止するか否かを確定する（図１７のステップ１７０６）。具体的には、反復フィルタ処理停止装置１６０３は、上述実施例／例示における反復フィルタ処理停止装置１０３、３０３、８０３、９０３又は１３０３と同様な機能と構造（例えば、上述の図２〜図１５に示す方法を参考にして反復フィルタ処理停止時間を確定）を備える。ここでは、重複する説明は省略する。
【０１０６】
勾配閾値確定装置１６０５は、フィルタ装置１６０１にて、各フィルタ処理前に、入力画像又は一回前のフィルタ処理で得られたフィルタリング画像における各画像単位の勾配値（図１７のステップ１７１０）を計算するためのものであり、画像中のあらゆる画像単位の勾配分布に基づいて勾配閾値Ｋを確定する（図１７のステップ１７１２）。ここで、“画像単位の勾配値”とは、その画像単位の隣接領域に対する階調（Gray scale）変化値をいう。
【０１０７】
ここでの画像単位とは、画像の画素であり、例えば、２次元画像の画素（pixel）又は３次元画像の画素（voxel）等である。図１８は、画像の各画像単位の勾配値分布を示す一例である。図１８では、横軸は画像単位の勾配値を表し、縦軸は、勾配値に対応する画像単位の個数を表す。一例として、勾配閾値確定装置１６０５は、入力画像に対応する比例値を取得し、該比例値に従って、勾配分布に基づいて勾配閾値Kを計算し、画像における勾配値が勾配閾値Ｋよりも小さい画像単位の個数とあらゆる画像単位の個数との比が前記比例値に近づくようにする（当然ながら、該比例に対してその他の数学的変化をさせることは可能で、例えば、その反映画像の勾配値が勾配閾値Ｋよりも小さい画像単位の個数と勾配値が勾配閾値Ｋよりも小さくない画像単位の個数の比であっても良い）。
【０１０８】
一具体的例として、画像単位の勾配値を、小さいほうから大きいほう、又は大きいほうから小さいほうへ順番に並べ、その並べたものの中で前記比例値に対応する勾配値を勾配閾値とし、勾配値が勾配閾値Ｋよりも小さい画像単位の個数とすべての画像単位の個数との比が前記比例値へ近づくようにする。他の具体的例として、画像の画像単位の勾配のヒストグラム（Histogram）を計算し、ヒストグラムに基づいて勾配閾値を決め、勾配値が勾配閾値Ｋよりも小さい画像単位の個数とすべての画像単位の個数との比が前記比例値へ近づくようにする。
【０１０９】
ここでの比例値とは、入力画像の類型又は特性（例えば、医学画像の画像生成装置又はスキャンプロトコル等）に基づいて予め決められた値であり、Ｐ（０＜Ｐ＜１）を使って表すことができる。各種類型の画像は同じ比例値に対応可能で、異なる画像類型は異なる比例値に対応可能である。ここで、実際の応用場面では、各種画像類型に対して、複数の該類型の画像を分析して予め相応の比例値を確定することにより、後続の非線形拡散フィルタ処理において使用できる。ここでは詳細な説明は省略する。例えば、予め決めた複数の画像類型に対応する複数の比例値（例えば、プロファイル形式）をメモリ部（該メモリ部は、画像処理装置に内在するメモリ装置であっても良いし、画像処理装置外部のメモリ装置であって該画像処理装置を経由して送るようにしても良い）へ保存する。実際の非線形拡散フィルタ処理において、勾配閾値確定装置１６０５は、入力画像の類型又は特性（例えば、画像生成装置の種類又はスキャンプロトコル等）に対応づけて保存している相応の比例値を取得する。こうすることで、人が手動で該パラメータを設定することなく、更に画像処理の自動化レベルを向上させることができる。異なる応用場面では該比例値の大小は異なるので、ここでも、該比例値の数字は具体的に限定されない。
【０１１０】
勾配閾値確定装置１６０５は、いかなる適切な方法を用いて画像の画像単位の勾配値を計算し、画像における画像単位の勾配分布に基づいて勾配閾値を確定しても良く、一例として、例えば、画像をｎ次元画像（ｎ（n≧１）とすれば、以下の式に基づいて画像単位の勾配値を計算できる。
【０１１１】
【数１３】
【０１１２】
ここで、Ｄ1，Ｄ2，・・・，Ｄnは、画像の次元数（ｎ次元）を表し、｜Ｉgrad（Ｄ1，Ｄ2，・・・，Ｄn）｜は、画像上の位置（Ｄ1，Ｄ2，・・・，Ｄn）における画素の勾配値を表し、さらに、ＩＤk（１≦ｋ≦ｎ）は、その画素の第Ｄk次元上の導関数を表す。
【０１１３】
各反復において、勾配閾値確定装置１６０５は、計算して得られる勾配閾値をフィルタ装置１６０１へ出力し、フィルタ装置１６０１は、勾配閾値確定装置１６０５にて決めた勾配閾値に基づいて非線形拡散フィルタ処理を実施する（図１７のステップ１７０２）。
【０１１４】
本実施形態の実施例に基づく画像処理方法と装置は、各種画像のフィルタ処理に応用することができる。これらの画像は１次元又は多次元（例えば、２次元又は３次元等）であって、例えば、本実施形態の実施例に基づく画像処理方法と装置は、医学画像（例えば、医療画像診断装置により取得されたデータに基づいて生成された画像）のノイズ軽減処理に応用可能である。一例として、上述方法の各ステップ及び上述装置の各構成及び／又は部分は、医療画像診断装置（例えば、Ｘ線診断装置、超音波診断装置、ＣＴ装置、磁気共鳴診断装置、ＰＥＴ装置等）のソフトウエア、ファームウエア、ハードウエア又はそれらの組み合わせとして、また、該医療画像診断装置の中の一部分として実施することができる。一例として、既存の医療画像診断装置において、本実施形態の上述方法及び／又は装置に基づいて実施することができ、しかも、既存の医療画像診断装置の各構成部分に対してある程度改良をするだけで即実施可能となる。その他の例として、上述方法の各ステップ及び上述装置の各構成及び／又は部分は、前記医療画像診断装置とは独立した装置として実施してもよい。上述装置の各構成や部分は、ソフトウエア、ファームウエア、ハードウエア又はそれらの組み合わせた方式により構成することができるが、使用する具体的手段や方式は当業者の周知のものであり、ここでは特に述べることはしない。
【０１１５】
一例として、上述方法の各ステップ及び上述装置の各構成及び／又は部分はソフトウエア、ファームウエア、ハードウエア又はそれらの組み合わせとして実施しても良い。ソフトウエア又はファームウエアを介して実現した場合、上述方法のソフトウエアプログラムを実施するため、メモリ媒体から又はネットワークを介して専用のハードウエア構造のコンピュータ（例えば、図１９に示す汎用コンピュータ１９００）へダウンロードして構成することができ、該コンピュータに各種プログラムがダウンロードされた状態で、各種機能等を実施することができる。
【０１１６】
図１９において、演算処理部（即ち、ＣＰＵ）１９０１は、読み出し専用メモリ（ＲＯＭ）１９０２の中に記憶されているプログラム、又は、記憶部１９０８から読み書き兼用メモリ（ＲＡＭ）１９０３へ書き込まれたプログラムに基づいて、各種処理を実施する。ＲＡＭ１９０３では、必要に応じて、ＣＰＵ１９０１が各種処理等を実施するときに必要なデータも記憶しておく。ＣＰＵ１９０１、ＲＯＭ１９０２及びＲＡＭ１９０３は、綜合ライン１９０４を経由してそれぞれ接続されている。入力／出力インターフェース１９０５も、綜合ライン１９０４につながっている。
【０１１７】
下記の各部は、入力／出力インターフェース１９０５に接続されている：入力部１９０６（キーボード、マウス等を含む）、出力部１９０７（モニタ、例えば、ブラウン管（ＣＲＴ）、液晶モニタ（ＬＣＤ）等や、スピーカ等を含む）、記憶部１９０８（キーボードを含む）、通信部１９０９（ネットワークインターフェースカード、例えば、ＬＡＮカード、モデム等）。通信部１９０９は、ネットワーク（例えば、インターネット）を介して通信処理を実施する。必要に応じて、駆動器１９１０も入力／出力インターフェース１９０５に接続可能である。取り外し可能な媒体１９１１は、例えば、磁気ディスク、光ディスク、ＭＯ、半導体メモリ等であって、必要に応じて駆動器１９１０に装着され、必要に応じてコンピュータプログラムを読み出して、記憶部１９０８へダウンロードされる。
【０１１８】
ソフトウエアを介して上述システム処理を実施する場合、ネットワーク（例えば、インターネット又は記憶媒体（例えば、取り外し可能な媒体１９１１））からプログラムをダウンロードしても良い。
【０１１９】
当業者においては、このような記憶媒体は図１９に示すようなプログラムを記憶した記憶媒体は、装置とは離れたところからユーザにプログラムを提供する取り外し可能な媒体１９１１に限らない。取り外し可能な媒体１９１１の例としては、磁気ディスク（フロッピー（登録商標）ディスク、光ディスク（ＣＤ−ＲＯＭやＤＶＤを含む）、磁気光ディスク（ＭｉｎｉＤｉｓｃ（ＭＤ、登録商標）を含む）らを含む。また、記憶媒体はＲＯＭ１９０２であっても良く、記憶部１９０８に含まれるハードディスク等、その中にプログラムが記憶され、それらを含む装置からユーザへプログラムが送られる形態でも良い。
【０１２０】
本実施形態では、更に、メモリとして、機器読み取り可能なコマンドコードを記憶しているプログラム製品でも応用でき、前記コマンドコードが機器を介して読み取られると、本実施形態の実施例における方法が実施される。
【０１２１】
上述機器読み取り可能なコマンドコードを記憶しているプログラム製品を受け入れるための記憶媒体も本実施形態に適用できる。その記憶媒体は、ハードディスク、光ディスク、磁気光ディスク、メモリカード、メモリスティックには限定されない。
【０１２２】
上記の本実施形態の具体的実施例においては、一つの実施方式に示す特徴について、同様の方法を一つ又は複数のほかの実施方法方の中で適用したり、その他の実施方法と組み合わせたり、又はその他の実施方法における特徴に替えるといったことも可能である。
【０１２３】
さらに、“包含する／含む”といった用語を使用したときは、特徴・構成・ステップ又は構造の存在を指し示す。ただし、その他の特徴・構成・ステップ又は構造の存在や付加の排除を意味するものではない。
【０１２４】
上記実施例においては、数字構成の図番記号を用いて各ステップや構成を表記している。ただし、これらの図番記号は単なる説明や画図の都合への考慮によるものであって、その順序やいかなるほかの限定を表すものではない、と当業者は理解すべきである。
【０１２５】
このほか、本実施形態の方法は、詳細な説明の欄において説明された時間順序に沿って実施されるものに限らず、その他の時間順序に沿って、同時に、又は独立して実施されても良い。それゆえ、本願の詳細な説明において説明された方法の実施順序は、本実施形態の技術範囲に対する構成を制限するものではない。
【０１２６】
上記では、既に、本実施形態の具体的実施例の説明をもって、本実施形態の説明を行っているものの、上述のすべての実施例はすべて単なる例示に過ぎず、限定するものではない。当業者は、特許請求の主旨や範囲において、本実施形態の各種手直し・改良又は同等物の設計を行うことが可能である。これらの手直し・改良又は同等物は、本実施形態の保護範囲内に含まれるものである。
【符号の説明】
【０１２７】
１００ 画像処理装置
１０１ フィルタ装置
１０３ 反復フィルタ処理停止装置
【特許請求の範囲】
【請求項１】
入力画像に対して反復フィルタ処理を実施するフィルタ装置と、
前記反復フィルタ処理の過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する反復フィルタ処理停止装置と、
を備えたことを特徴とする画像処理装置。
【請求項２】
前記反復フィルタ処理停止装置は、
前記フィルタ装置が各フィルタ処理を実施する毎に、前記フィルタ処理によって得られたフィルタリング画像及び前記入力画像に基づいて、ノイズ画像を推定するノイズ画像取得部と、
前記ノイズ画像と前記フィルタリング画像との相関性を推定する相関性取得部と、
前記反復フィルタ処理の過程において前記相関性が減少する変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する停止判断部と、
を備えたことを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記停止判断部は、前記反復フィルタ処理の過程において前記相関性が減少する変化率を示す特徴値を計算し、前記特徴値と所定の閾値との間で所定の関係を満足するか否かを判定し、前記特徴値と前記所定の閾値との間で前記所定の関係を満足していた場合に、前記反復フィルタ処理を停止することを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記相関性取得部は、前記相関性を示すために、前記ノイズ画像と前記フィルタリング画像との間の相関係数を計算することを特徴とする請求項３に記載の画像処理装置。
【請求項５】
前記特徴値は、前記反復フィルタ処理の過程における反復回数の増加に伴う前記相関係数の減少を表す曲線の傾斜度を示す値であることを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記所定の閾値は、前記入力画像に基づく特性に基づいて予め決められた値であることを特徴とする請求項３〜５のいずれか一つに記載の画像処理装置。
【請求項７】
前記反復フィルタ処理停止装置は、
各フィルタ処理毎に、フィルタリング画像の所定の画像特徴値を計算する画像特徴計算部と、
前記所定の画像特徴の反復回数に伴う増大変化の変化率に基づいて、前記反復フィルタ処理を停止するか否かを確定する停止判断部と、
を具備することを特徴とする請求項１に記載の画像処理装置。
【請求項８】
前記反復フィルタ処理停止装置は、
前記フィルタ装置が各フィルタ処理を実施する毎に、前記フィルタ処理によって得られたフィルタリング画像と前記入力画像に基づいて、ノイズ画像を推定するノイズ画像取得部をさらに備え、
前記画像特徴計算部は、前記フィルタリング画像の前記所定の画像特徴値を示すために、前記ノイズ画像の所定の画像特徴値を計算することを特徴とする請求項７に記載の画像処理装置。
【請求項９】
前記所定の画像特徴は、画像の分散又は標準差であることを特徴とする請求項７又は８に記載の画像処理装置。
【請求項１０】
前記停止判断部は、前記反復フィルタ処理過程における前記所定の画像特徴の変化の変化率を表す特徴値を計算し、前記特徴値と所定の閾値との間で所定の関係を満足するか否かを判断し、前記特徴値と前記所定の閾値との間で前記所定の関係を満足していた場合に、前記反復フィルタ処理を停止することを特徴とする請求項７又は８に記載の画像処理装置。
【請求項１１】
前記特徴値は、前記反復フィルタ処理過程における反復回数の増加に伴う前記所定の画像特徴の変化を表す曲線の傾斜度を示す値であることを特徴とする請求項１０に記載の画像処理装置。
【請求項１２】
前記反復フィルタ処理停止装置は、
前記ノイズ画像の所定の画像特徴値を計算する画像特徴計算部と、
前記相関性と前記ノイズ画像の前記所定の画像特徴値とを組み合わせて、組み合わせ特徴値を得る組み合わせ部とさらに備え、
前記停止判断部は、前記反復フィルタ処理の過程において前記組み合わせ特徴値が極小値に達した時に、前記反復フィルタ処理を停止することを特徴とする請求項２に記載の画像処理装置。
【請求項１３】
前記所定の画像特徴は、画像の分散又は標準差であることを特徴とする請求項１２に記載の画像処理装置。
【請求項１４】
前記画像特徴計算部は、変換後の所定の画像特徴値を得るために、計算して得られた前記ノイズ画像の所定の画像特徴値に対して変換を実施し、
前記組み合わせ部は、前記相関性と変換後の所定の画像特徴値との和を計算することを特徴とする請求項１２又は１３に記載の画像処理装置。
【請求項１５】
前記相関性取得部は、前記相関性を示すために、前記ノイズ画像と前記フィルタリング画像との間の相関係数を計算することを特徴とする請求項１２又は１３に記載の画像処理装置。
【請求項１６】
前記フィルタ装置は、非線形拡散方法に基づいてフィルタ処理を実施することを特徴とする請求項１〜５のいずれか一つに記載の画像処理装置。
【請求項１７】
フィルタ装置が各フィルタ処理を実施する前に、前記入力画像又は前回のフィルタ処理で得られたフィルタリング画像における各一画像単位の勾配値を計算し、前記画像における全ての画像単位の勾配分布に基づいて勾配閾値を確定する勾配閾値確定装置をさらに備え、
前記フィルタ装置は、前記勾配閾値確定装置が確定した勾配閾値を利用して、非線形拡散方式に基づいてフィルタ処理を実施することを特徴とする請求項１６に記載の画像処理装置。
【請求項１８】
前記勾配閾値確定装置は、前記入力画像に応じた比例値を取得し、勾配度分布及び前記比例値に基づいて前記勾配閾値を計算し、
画像における前記勾配閾値より小さい勾配値を有する画像単位の個数と、画像における全ての画像単位の個数との比が、前記比例値と略同一であることを特徴とする請求項１７に記載の画像処理装置。
【請求項１９】
前記比例値は、前記入力画像の特性に基づいて予め決められた値であることを特徴とする請求項１８に記載の画像処理装置。
【請求項２０】
前記入力画像は、医療画診断装置により取得されたデータに基づいて生成された医学画像であることを特徴とする請求項１〜５のいずれか一つに記載の画像処理装置。
【請求項２１】
入力画像に対して反復フィルタ処理を実施するステップと、
前記反復フィルタ処理過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定するステップと、
を含んだことを特徴とする画像処理方法。
【請求項２２】
入力画像に対して反復フィルタ処理を実施する手順と、
前記反復フィルタ処理過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定する手順と、
をコンピュータに実行させることを特徴とする画像処理プログラム。
【請求項１】
入力画像に対して反復フィルタ処理を実施するフィルタ装置と、
前記反復フィルタ処理の過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する反復フィルタ処理停止装置と、
を備えたことを特徴とする画像処理装置。
【請求項２】
前記反復フィルタ処理停止装置は、
前記フィルタ装置が各フィルタ処理を実施する毎に、前記フィルタ処理によって得られたフィルタリング画像及び前記入力画像に基づいて、ノイズ画像を推定するノイズ画像取得部と、
前記ノイズ画像と前記フィルタリング画像との相関性を推定する相関性取得部と、
前記反復フィルタ処理の過程において前記相関性が減少する変化速度に基づいて、前記反復フィルタ処理を停止するか否かを確定する停止判断部と、
を備えたことを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記停止判断部は、前記反復フィルタ処理の過程において前記相関性が減少する変化率を示す特徴値を計算し、前記特徴値と所定の閾値との間で所定の関係を満足するか否かを判定し、前記特徴値と前記所定の閾値との間で前記所定の関係を満足していた場合に、前記反復フィルタ処理を停止することを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記相関性取得部は、前記相関性を示すために、前記ノイズ画像と前記フィルタリング画像との間の相関係数を計算することを特徴とする請求項３に記載の画像処理装置。
【請求項５】
前記特徴値は、前記反復フィルタ処理の過程における反復回数の増加に伴う前記相関係数の減少を表す曲線の傾斜度を示す値であることを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記所定の閾値は、前記入力画像に基づく特性に基づいて予め決められた値であることを特徴とする請求項３〜５のいずれか一つに記載の画像処理装置。
【請求項７】
前記反復フィルタ処理停止装置は、
各フィルタ処理毎に、フィルタリング画像の所定の画像特徴値を計算する画像特徴計算部と、
前記所定の画像特徴の反復回数に伴う増大変化の変化率に基づいて、前記反復フィルタ処理を停止するか否かを確定する停止判断部と、
を具備することを特徴とする請求項１に記載の画像処理装置。
【請求項８】
前記反復フィルタ処理停止装置は、
前記フィルタ装置が各フィルタ処理を実施する毎に、前記フィルタ処理によって得られたフィルタリング画像と前記入力画像に基づいて、ノイズ画像を推定するノイズ画像取得部をさらに備え、
前記画像特徴計算部は、前記フィルタリング画像の前記所定の画像特徴値を示すために、前記ノイズ画像の所定の画像特徴値を計算することを特徴とする請求項７に記載の画像処理装置。
【請求項９】
前記所定の画像特徴は、画像の分散又は標準差であることを特徴とする請求項７又は８に記載の画像処理装置。
【請求項１０】
前記停止判断部は、前記反復フィルタ処理過程における前記所定の画像特徴の変化の変化率を表す特徴値を計算し、前記特徴値と所定の閾値との間で所定の関係を満足するか否かを判断し、前記特徴値と前記所定の閾値との間で前記所定の関係を満足していた場合に、前記反復フィルタ処理を停止することを特徴とする請求項７又は８に記載の画像処理装置。
【請求項１１】
前記特徴値は、前記反復フィルタ処理過程における反復回数の増加に伴う前記所定の画像特徴の変化を表す曲線の傾斜度を示す値であることを特徴とする請求項１０に記載の画像処理装置。
【請求項１２】
前記反復フィルタ処理停止装置は、
前記ノイズ画像の所定の画像特徴値を計算する画像特徴計算部と、
前記相関性と前記ノイズ画像の前記所定の画像特徴値とを組み合わせて、組み合わせ特徴値を得る組み合わせ部とさらに備え、
前記停止判断部は、前記反復フィルタ処理の過程において前記組み合わせ特徴値が極小値に達した時に、前記反復フィルタ処理を停止することを特徴とする請求項２に記載の画像処理装置。
【請求項１３】
前記所定の画像特徴は、画像の分散又は標準差であることを特徴とする請求項１２に記載の画像処理装置。
【請求項１４】
前記画像特徴計算部は、変換後の所定の画像特徴値を得るために、計算して得られた前記ノイズ画像の所定の画像特徴値に対して変換を実施し、
前記組み合わせ部は、前記相関性と変換後の所定の画像特徴値との和を計算することを特徴とする請求項１２又は１３に記載の画像処理装置。
【請求項１５】
前記相関性取得部は、前記相関性を示すために、前記ノイズ画像と前記フィルタリング画像との間の相関係数を計算することを特徴とする請求項１２又は１３に記載の画像処理装置。
【請求項１６】
前記フィルタ装置は、非線形拡散方法に基づいてフィルタ処理を実施することを特徴とする請求項１〜５のいずれか一つに記載の画像処理装置。
【請求項１７】
フィルタ装置が各フィルタ処理を実施する前に、前記入力画像又は前回のフィルタ処理で得られたフィルタリング画像における各一画像単位の勾配値を計算し、前記画像における全ての画像単位の勾配分布に基づいて勾配閾値を確定する勾配閾値確定装置をさらに備え、
前記フィルタ装置は、前記勾配閾値確定装置が確定した勾配閾値を利用して、非線形拡散方式に基づいてフィルタ処理を実施することを特徴とする請求項１６に記載の画像処理装置。
【請求項１８】
前記勾配閾値確定装置は、前記入力画像に応じた比例値を取得し、勾配度分布及び前記比例値に基づいて前記勾配閾値を計算し、
画像における前記勾配閾値より小さい勾配値を有する画像単位の個数と、画像における全ての画像単位の個数との比が、前記比例値と略同一であることを特徴とする請求項１７に記載の画像処理装置。
【請求項１９】
前記比例値は、前記入力画像の特性に基づいて予め決められた値であることを特徴とする請求項１８に記載の画像処理装置。
【請求項２０】
前記入力画像は、医療画診断装置により取得されたデータに基づいて生成された医学画像であることを特徴とする請求項１〜５のいずれか一つに記載の画像処理装置。
【請求項２１】
入力画像に対して反復フィルタ処理を実施するステップと、
前記反復フィルタ処理過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定するステップと、
を含んだことを特徴とする画像処理方法。
【請求項２２】
入力画像に対して反復フィルタ処理を実施する手順と、
前記反復フィルタ処理過程において各フィルタ処理毎に得られるフィルタリング画像の反復回数に伴う変化速度に基づいて前記反復フィルタ処理を停止するか否かを確定する手順と、
をコンピュータに実行させることを特徴とする画像処理プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【公開番号】特開２０１２−１７４２７５（Ｐ２０１２−１７４２７５Ａ）
【公開日】平成２４年９月１０日（２０１２．９．１０）
【国際特許分類】
【出願番号】特願２０１２−３４２２５（Ｐ２０１２−３４２２５）
【出願日】平成２４年２月２０日（２０１２．２．２０）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【出願人】（５９４１６４５４２）東芝メディカルシステムズ株式会社 (4,066)
【Ｆターム（参考）】
【公開日】平成２４年９月１０日（２０１２．９．１０）
【国際特許分類】
【出願日】平成２４年２月２０日（２０１２．２．２０）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【出願人】（５９４１６４５４２）東芝メディカルシステムズ株式会社 (4,066)
【Ｆターム（参考）】
[ Back to top ]