(57)【特許請求の範囲】
【請求項１】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記領域設定部は、
前記第１画像の画素値に基づいて、前記第１領域を設定し、
前記第２画像の画素値に基づいて、前記第２領域を設定すること、
を特徴とする医用画像処理装置。
【請求項２】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記周波数分布発生部は、
前記第１領域周波数分布に対応するパワースペクトルを平滑化することにより第１平滑スペクトルを発生し、
前記第２領域周波数分布に対応するパワースペクトルを平滑化することにより第２平滑スペクトルを発生し、
前記補正部は、
前記第１平滑スペクトルを前記第２平滑スペクトルに変換する変換関数を決定し、
前記変換関数を用いて、前記第１周波数分布を補正すること、
を特徴とする医用画像処理装置。
【請求項３】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記周波数分布発生部は、
前記第１領域周波数分布に対応するパワースペクトルを平滑化することにより第１平滑スペクトルを発生し、
前記第２領域周波数分布に対応するパワースペクトルを平滑化することにより第２平滑スペクトルを発生し、
前記補正部は、
前記第２平滑スペクトルを前記第１平滑スペクトルに変換する変換関数を決定し、
前記変換関数を用いて、前記第２周波数分布を補正すること、
を特徴とする医用画像処理装置。
【請求項４】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記周波数分布発生部は、
前記第１領域周波数分布に対応するパワースペクトルを平滑化することにより第１平滑スペクトルを発生し、
前記第２領域周波数分布に対応するパワースペクトルを平滑化することにより第２平滑スペクトルを発生し、
前記補正部は、
前記第１平滑スペクトルを前記所定の周波数分布に対応するパワースペクトルに変換する第１変換関数を決定し、
前記第２平滑スペクトルを前記所定の周波数分布に対応するパワースペクトルに変換する第２変換関数を決定し、
前記第１変換関数を用いて前記第１周波数分布を補正し、
前記第２変換関数を用いて前記第２周波数分布を補正すること、
を特徴とする医用画像処理装置。
【請求項５】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記周波数分布発生部は、
前記第１領域周波数分布に対応する第１パワースペクトルの最大値と前記第２領域周波数分布に対応する第２パワースペクトルの最大値とを用いて、前記第１パワースペクトルと前記第２パワースペクトルとを正規化し、
前記補正部は、
前記第１パワースペクトル、前記第２パワースペクトル、または前記所定の周波数分布に対応するパワースペクトルを基準として、前記第１周波数分布と前記第２周波数分布とのうち少なくとも一方を補正すること、
を特徴とする医用画像処理装置。
【請求項６】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記第１画像において複数の第１エッジを検出し、前記第２画像において複数の第２エッジを検出するエッジ検出部を更に具備し、
前記領域設定部は、
前記第１領域における対向する２辺が前記第１エッジに直交し、かつ所定の閾値により規定される画素値範囲に対応する領域を、前記第１領域として設定し、
前記第２領域における前記２辺が前記第２エッジに直交し、かつ前記画素値範囲に対応する領域を、前記第２領域として設定すること、
を特徴とする医用画像処理装置。
【請求項７】
第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、
を具備し、
前記補正部は、
前記第１領域周波数分布および前記第２領域周波数分布に関するフィルター関数と補間関数と補間回数とに基づいて補正関数を決定し、決定した補正関数を用いて、前記第１周波数分布と前記第２周波数分布とのうち少なくとも一方を補正すること、
を特徴とする医用画像処理装置。
【請求項８】
前記領域設定部は、
前記第１画像における空間的に均一な領域または経時的に不変な領域を、前記第１領域として設定し、
前記第２画像における空間的に均一な領域、または経時的に不変な領域を、前記第２領域として設定すること、
を特徴とする請求項１乃至７のうちいずれか一項に記載の医用画像処理装置。
【請求項９】
第１画像の一部分に対応する第１領域を前記第１画像の画素値に基づいて設定し、
前記第１画像より過去に収集された第２画像の一部分に対応する第２領域を前記第２画像の画素値に基づいて設定し、
前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、
前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生し、
前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正すること、
を具備することを特徴とする医用画像処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、医用画像を補正する医用画像処理装置および医用画像処理方法に関する。
【背景技術】
【０００２】
脳磁気共鳴イメージング（Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  Ｉｍａｇｉｎｇ：ＭＲＩ）において、Ｔ１Ｗ、Ｔ２Ｗなどの脳の３次元のＶｏｌｕｍｅ  ｄａｔａに対して、ＳＰＭ（Ｓｔａｔｉｓｔｉｃａｌ  Ｐａｒａｍｅｔｒｉｃ  Ｍａｐｐｉｎｇ）などの統計解析ツールを用いたＶｏｘｅｌ−Ｂａｓｅｄ  Ｍｏｒｐｈｏｍｅｔｒｙ  （ＶＢＭ）と称する分野がある。ＶＢＭでは、位置合わせ後に個人脳の標準脳への解剖学的標準化が行われるが、一般に強いスムージングを施し微細な変化を消すようにしている。すなわち、ＶＢＭでは、対象画像を意図的にぼかすことで、標準脳への解剖学的標準化（位置合わせ）が実行される。
【０００３】
このＶＢＭという分野では、時間経過に伴う形態の微細な変化を追跡する必要がある認知症や精神神経疾患に応用されてきている。例として、アルツハイマー（Ａｌｚｈｅｉｍｅｒ：ＡＤ）など痴呆症の診断などでは、皮質の厚さをサブｍｍのオーダーで、また海馬などの脳容積を数％のオーダーで、ＭＲＩを用いて、経時的な微小な変化を、１〜数ヶ月単位で数年間〜生涯に亘って随時追跡する必要がある。
【０００４】
また、ＶＢＭと類似の解析手法は、構造データのみならず、拡散画像、機能的（ｆｕｎｃｔｉｏｎａｌ）ＭＲＩ（ｆＭＲＩ）、脳血流（Ａｒｔｅｒｉａｌ  Ｓｐｉｎ  Ｌａｂｅｌｉｎｇ：ＡＳＬ）などの機能データでも広く応用されてきている。モダリティ別ではＭＲＩ装置に限定されず単光子放出コンピュータ断層撮影（Ｓｉｎｇｌｅ  Ｐｈｏｔｏｎ  Ｅｍｉｓｓｉｏｎ  Ｃｏｍｐｕｔｅｄ  Ｔｏｍｏｇｒａｐｈｙ：ＳＰＥＣＴ）装置、陽電子放出コンピュータ断層撮影（Ｐｏｓｉｔｉｒｏｎ  Ｅｍｉｓｓｉｏｎ  Ｃｏｍｐｕｔｅｄ  Ｔｏｍｏｇｒａｐｈｙ：ＰＥＴ）装置、Ｘ線ＣＴ（Ｘ線コンピュータ断層撮影）装置などでの応用も多い。
【０００５】
胸部Ｘ線画像では、ある期間毎に撮像した画像を用いて、計時差分をおこなって、以前との変化分を把握しやすくする方法が実用化されている。この場合、骨などの情報を用いて、平面内での２次元的な撮像位置に加え、透過画像特有の奥行き方向の患者撮像体位の補正などが行われている。また、Ｘ線ＣＴ、ＳＰＥＣＴ、ＰＥＴ、ＭＲＩなどの断面像を取得可能な撮像装置においても、以前の断面に切りなおすｒｅ−ｓｌｉｃｉｎｇ程度のことは行われている。
【０００６】
しかしながら、撮像時刻が異なる複数の医用画像において、装置の違いや撮像および画像処理に関するフィルターやマトリクスサイズの違いがあると、患者由来の経時変化因子による正しい経時的変化の追跡が難しい問題がある。また、同一装置条件（撮像条件）であっても、撮像体位の変化を補正する場合、画像の補間を行う必要がある。たとえフィルターが等価であったとしても、例えば、補間回数の相違により、ボケなどの効果が現れてしまうことがあり、現在それらの影響は考慮されていない問題がある。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特許第４０２２５８７号公報
【特許文献２】特許第５５３２７３０号公報
【特許文献３】特開２０１４−０４２６８４号公報
【特許文献４】特開２０１０−０４６１０３号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
目的は、対象物の解剖学的変化を精度よく追跡可能にするために、医用画像を補正可能な医用画像処理装置と医用画像処理方法とを提供することにある。
【課題を解決するための手段】
【０００９】
本実施形態に係る医用画像処理装置は、第１画像の一部分に対応する第１領域と、前記第１画像より過去に収集された第２画像の一部分に対応する第２領域とを設定する領域設定部と、前記第１領域における複数の画素値に基づいて第１領域周波数分布を発生し、前記第２領域における複数の画素値に基づいて第２領域周波数分布を発生する周波数分布発生部と、前記第１領域周波数分布、前記第２領域周波数分布、または所定の周波数分布を基準として、前記第１画像に対応する第１周波数分布と前記第２画像に対応する第２周波数分布とのうち少なくとも一方を補正する補正部と、を具備し、前記領域設定部は、前記第１画像の画素値に基づいて、前記第１領域を設定し、前記第２画像の画素値に基づいて、前記第２領域を設定することを特徴とする。
【図面の簡単な説明】
【００１０】
【図１】図１は、第１の実施形態に係る医用画像処理装置の構成を示す構成図である。
【図２】図２は、第１の実施形態に係り、第１補間画像および第２補間画像において、第１領域または第２領域の設定の一例を示す図である。
【図３】図３は、第１の実施形態に係り、第１領域または第２領域において、位置（例えば、空気部分）に対する画素値の変化（ノイズプロファイル）の一例を示す図である。
【図４】図４は、第１の実施形態に係り、ノイズプロファイルの逆フーリエ変換を用いて発生された第１領域周波数分布に対応するパワースペクトル、第２領域周波数分布に対応するパワースペクトル、およびこれらのパワースペクトルをそれぞれ平滑化した第１平滑スペクトル、第２平滑スペクトルの一例を示す図である。
【図５】図５は、第１の実施形態に係り、第１平滑スペクトルを第２平滑スペクトルに変換する変換関数の一例を示す図である。
【図６】図６は、第１の実施形態に係り、変換関数を用いて補正された第１周波数分布に対応するパワースペクトルの一例を示す図である。
【図７】図７は、第１の実施形態に係り、周波数分布補正処理の手順の一例を示すフローチャートである。
【図８】図８は、第２の実施形態に係る医用画像処理装置の構成を示す構成図である。
【図９】図９は、第２の実施形態に係り、第１補間画像および第２補間画像において、第１領域または第２領域の設定の一例を示す図である。
【図１０】図１０は、第２の実施形態に係り、図９における線分ｅｅ’上の画素値の変化（実信号プロファイル）の一例を示す図である。
【図１１】図１１は、第２の実施形態に係り、図９における線分ｄｄ’または図１０における抽出部分（ｆｆ’）上の画素値の変化（経時不変プロファイル）の一例を示す図である。
【図１２】図１２は、第２の実施形態に係り、経時不変プロファイルの逆フーリエ変換の絶対値の正規化により発生された第１正規化周波数分布と第２正規化周波数分布との一例を示す図である。
【図１３】図１３は、第２の実施形態に係り、第１正規化周波数分布を第２正規化周波数分布に変換する変換関数の一例を示す図である。
【図１４】図１４は、第２の実施形態に係り、変換関数を用いて補正された第１周波数分布の絶対値の一例を示す図である。
【図１５】図１５は、第２の実施形態に係り、周波数分布補正処理の手順の一例を示すフローチャートである。
【発明を実施するための形態】
【００１１】
以下、本実施形態に係る医用画像処理装置について、図面を参照して説明する。なお、以下の説明において、略同一の機能および構成を有する構成要素については、同一符号を付し、重複説明は必要な場合に行う。
【００１２】
（第１の実施形態）  
図１は、第１の実施形態に係る医用画像処理装置１の構成を示すブロック図である。医用画像処理装置１は、補間処理部１１と、領域設定部１３と、周波数分布発生部１５と、補正部１７と、画像発生部１９と、表示部２１と、インターフェース部２３と、入力部２５と、記憶部２７と、制御部２９とを有する。図１におけるＡは、本医用画像処理装置１におけるインターフェース部２３を介して、外部の各種モダリティ、医用画像保管装置などと接続するネットワークを表している。
【００１３】
入力部２５を介した操作者の指示により、第１画像と第２画像とが選択されると、補間処理部１１は、第１画像と第２画像とに対して（補間関数を用いた）補間処理を実行する。補間処理とは、例えば、第１画像と第２画像とにおいて、ボクセル（またはピクセル）サイズの統一、空間位置補正、位置合わせなどである。第１画像は、第２画像に対して、比較対象となる医用画像（最新または現在の医用画像）である。第２画像は、第１画像に対して基準（Ｂａｓｅ）となる医用画像である。具体的には、第２画像は、通常は現在（第１画像）と比較され、第１画像よりも過去に取得された医用画像、または過去に収集されたボリュームデータを基準断面に対して再スライス（ｒｅ−ｓｌｉｃｅ）された医用画像（例えば、ＭＰＲ（Ｍｕｌｔｉ−Ｐｌａｎａｒ  Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ）画像）である。
【００１４】
第１画像および第２画像は、被検体に関する同一の解剖学的領域を有する医用画像である。以下、第１画像および第２画像は、脳画像として説明する。なお、第１画像および第２画像は、脳画像に限定されず、任意の解剖学的領域を有する画像であってもよい。第１画像と第２画像とは、例えば、癌に対応する領域を有する医用画像、骨粗鬆症に関する医用画像、癌の進行の評価（浸潤度）として用いられる医用画像、フォローアップ検査に関する医用画像などである。また、第１画像および第２画像は、同一モダリティにより発生された医用画像である。
【００１５】
具体的には、第２画像の選択の後に、補間処理部１１は、例えば、第２画像に基づいて、基準（Ｂａｓｅ）となる画像条件を決定する。画像条件とは、例えば、マトリクスサイズ（ｍａｔｒｉｘ  ｓｉｚｅ）、撮像視野（Ｆｉｅｌｄ  ｏｆ  ｖｉｅｗ：ＦＯＶ）に基づくボクセル（ｖｏｘｅｌ）サイズ、患者の撮像体位などである。なお、基準となる画像条件は、第１画像の付帯情報に基づいて決定されてもよい。また、基準となる画像条件は、第２画像の付帯情報と第１画像の付帯情報とに基づいて決定されてもよい。
【００１６】
例えば、第１画像および第２画像における空間位置の相違に対する補間処理は、第１画像および第２画像が２次元画像なら面内での位置の補正、第１画像および第２画像が３次元画像ならスライス断面の補正に相当する。また、ボクセルサイズおよび空間位置の相違に対する処理は「補間」であるため、第１画像および第２画像に対して同様な処理（同様な補間関数を用いた補間処理）を実行することが望ましい。
【００１７】
なお、基準画像（第2画像）に対して対象画像（第1画像）を位置合わせする場合の補間処理は、極力高精度の周波数劣化の小さな手法で行うことが重要であると同時にその補間処理時のパワースペクトルの劣化を加味して変換フィルター特性を決める必要がある。
【００１８】
領域設定部１３は、補間処理後の第１画像（以下、第１補間画像と呼ぶ）の一部分に対応する第１領域を設定する。領域設定部１３は、補間処理後の第２画像（以下、第２補間画像と呼ぶ）に対応する第２領域を設定する。第１領域は、第１補間画像において、空間的に均一な領域である。第２領域は、第２補間画像において、空間的に均一である。
【００１９】
すなわち、第１部域および第２領域は、例えば、ノイズ部分に関する領域であって、無信号部分（空気領域）および構造がない（非構造であって）空間的に均一な部分に関する領域である。第１領域および第２領域は、例えば、脳画像において、空気に対応する領域および脳室に対応する領域などである。以下、第１領域および第２領域は、非構造であって空間的に均一な部分（無信号部分、ノイズ部分）に関する領域、すなわち空気領域とする。なお、第１領域および第２領域が経時的に変化の小さい領域（経時的に不変な領域）である場合については、第２の実施形態で詳述する。
【００２０】
図２は、第１補間画像および第２補間画像において、第１領域または第２領域の設定の一例を示す図である。図２におけるａは、補間画像における空気の領域を示している。図２におけるｂは、補間画像における骨の領域を示している。図２におけるｃは、補間画像における頭皮の領域を示している。図２におけるＧＭは、補間画像における灰白質（Ｇｒａｙ  Ｍａｔｔｅｒ）の領域を示している。図２におけるＷＭは、補間画像における白質（Ｗｈｉｔｅ  Ｍａｔｔｅｒ）の領域を示している。図２におけるＣＳＦは、補間画像における脳脊髄液（ＣｅｒｅｂｒｏＳｐｉｎａｌ  Ｆｌｕｉｄ：ＣＳＦ）を有する領域を示している。図２において、左右対称であって画像中心部分のＣＳＦは、脳室に対応する。図２におけるＮは、例えば、第１領域および第２領域が設定される領域に対応する。
【００２１】
領域設定部１３は、第1補間画像の画素値に基づいて第１領域を設定し、第２補間画像の画素値に基づいて前記第２領域を設定する。領域設定部１３は、例えば、第１補間画像における複数の画素値に対する閾値により、第１領域を第１補間画像に設定する。領域設定部１３は、第２補間画像における複数の画素値に対する閾値により、第２領域を第２補間画像に設定する。上記閾値とは、例えば、空気を抽出可能な画素値である。この\閾値は、記憶部２７に記憶される。
【００２２】
以下、説明を簡単にするために、第１領域と第２領域との形状は、矩形であるものとする。なお、第１領域と第２領域との形状は、矩形に限定されず、円形、楕円形等の任意の形状であってもよい。また、第１領域および第２領域は、空気の領域またはＣＳＦの領域であってもよい。なお、第１領域または第２領域を空気の領域に設定できない場合、領域設定部１３は、例えば、ＣＳＦの領域に第１領域と第２領域とを設定してもよい。
【００２３】
第１領域の大きさと第２領域の大きさとは、同一の大きさが好ましい。なお、第１領域の大きさと第２領域の大きさとは、異なっていてもよい。第１領域および第２領域の形状および大きさは、予め設定されているものとする。なお、第１領域および第２領域の形状および大きさは、入力部２５を介した操作者の指示により、適宜変更可能である。なお、第１領域と第２領域とは、入力部２５を介した操作者の指示により、設定されてもよい。第１領域の大きさと第２領域の大きさとは、例えば、数百の画素を有する。
【００２４】
なお、第１補間画像および第２補間画像が振幅画像である場合、信号対雑音比（Ｓｉｇｎａｌ−Ｎｏｉｓｅ  Ｒａｔｉｏ：ＳＮＲ）は十分大きい（ＳＮＲ＞５）ため、非空気領域（信号部）の複数の信号値は、ガウシアン（Ｇａｕｓｓｉａｎ）分布を有する。一方、空気領域（無信号部）における複数の信号値（ｎｏｉｓｅ）は、例えば、負値がない（非負の）ライシアン（Ｒｉｃｉａｎ）分布となる。このため、第１領域（および第２領域）の設定後の処理において、信号部と無信号部とを別々に扱うか、第１補間画像（第２補間画像）から第１補間画像（第２補間画像）のＤＣ（Ｄｉｒｅｃｔ  Ｃｕｒｅｎｔ：直流）成分（平均画素値）を差分した差分画像において（ＤＣ成分除去後）、絶対値化（ｍａｇｎｉｔｕｄｅ化）することにより、第１領域（および第２領域）における無信号部分をライシアンノイズに変換する必要がある。
【００２５】
以下、第１補間画像に設定された第１領域に含まれる複数の画素にそれぞれ対応する複数の画素値を、Ｉｏｒｉｇ（ｘ）とする。また、第２補間画像に設定された第２領域に含まれる複数の画素にそれぞれ対応する複数の画素値を、Ｉｂａｓｅ（ｘ）とする。説明を簡単にするために、第１領域に含まれる複数の画素の位置と、第２領域に含まれる複数の画素の位置とは同じであるものとする。このとき、第１補間画像と第２補間画像とは、補間処理部１１により、予め位置合わせされる。また、第１領域と第２領域とのうち一方が設定されると、他方が自動的に設定されてもよい。Ｉｏｒｉｇ（ｘ）とＩｂａｓｅ（ｘ）とにおける引数ｘは、第１領域および第２領域に含まれる画素の位置を示す。また、引数ｘが含まれる範囲は同一となる。なお、一般的に、第１領域と第２領域とが異なる位置および異なる大きさで設定された場合、引数はそれぞれ異なる文字（例えば、Ｉｏｒｉｇ（ｘ）とＩｂａｓｅ（ｙ）など）により規定され、引数の範囲も異なる。
【００２６】
図３は、第１領域（または第２領域）において、位置（例えば、空気部分）に対する画素値の変化（ノイズプロファイル）の一例を示す図である。図３の縦軸は、画素値を示している。図３の縦軸におけるａｖｅは、第１領域（または第２領域）の平均画素値（ＤＣ成分）に相当する。
【００２７】
第１領域および第２領域が２次元の場合、図３に対応するノイズプロファイルは、３次元的なグラフ（ｘ軸及びｙ軸は画素の位置、ｚ軸はノイズの大きさ）となる。また、第１画像および第２画像がボリュームデータである場合、図３に対応するノイズプロファイルは、４次元的なグラフ（ｘ軸、ｙ軸、ｚ軸は画素の位置、ｗ軸はノイズの大きさ）となる。
【００２８】
周波数分布発生部１５は、第１領域に含まれる複数の画素値に基づいて、第１領域の周波数分布の大きさの平方を示す第１パワースペクトルを発生し、第２領域における複数の画素値に基づいて、第２領域の周波数分布の大きさの平方を示す第２パワースペクトルを発生する。具体的には、周波数分布発生部１５は、第１領域および第２領域各々に含まれる複数の画素値に対して、逆フーリエ変換を実行する。次いで、周波数分布発生部１５は、逆フーリエ変換の結果に対して、絶対値および平方をそれぞれ適用することにより、第１パワースペクトルと第２パワースペクトルとを発生する。
【００２９】
第１パワースペクトルをＰｏｒｉｇ（ｆ）、第２パワースペクトルをＰｂａｓｅ（ｆ）と表記すると、周波数分布発生部１５において実行される計算は、例えば、以下のように表すことができる。ここで、ｆは、周波数であり第１パワースペクトルと第２パワースペクトルとにおける引数である。  
Ｐｏｒｉｇ（ｆ）＝｜ＩＦＴ［Ｉｏｒｉｇ（ｘ）］｜２ 
Ｐｂａｓｅ（ｆ）＝｜ＩＦＴ［Ｉｂａｓｅ（ｘ）］｜２  ・・・（１）  
上式におけるＩＦＴは、逆フーリエ変換を示している。
【００３０】
すなわち、周波数分布発生部１５は、第１領域および第２領域各々における複数の画素値（画像データ）に対して、第１（補間）画像の次元に応じた２次元または３次元の逆フーリエ変換を実行する。逆フーリエ変換に対する絶対値演算及び平方演算により、第１領域及び第２領域にそれぞれ対応する第１パワースペクトルと第２パワースペクトルとを発生する。
【００３１】
一般に、画像化時のゲインに依存して、第１画像と第２画像とにおける画素値のスケーリング（ｓｃａｌｉｎｇ）は異なる。このため、周波数分布発生部１５は、第１パワースペクトルと第２パワースペクトルとを、基準となるパワースペクトルの最大値（通常はＤＣ（ｆ＝０）成分）で正規化（ｎｏｒｍａｌｉｚｅ）する。説明を簡単にするために、基準となるパワースペクトルは、第２パワースペクトルの最大値Ｐｂａｓｅ（０）であるものとする。なお、基準となるパワースペクトルが第１パワースペクトルである場合、および所定の周波数分布に対応するパワースペクトルである場合については、後ほど変形例で説明する。
【００３２】
具体的には、周波数分布発生部１５は、第１パワースペクトルＰｏｒｉｇ（ｆ）を第２パワースペクトルの最大値Ｐｂａｓｅ（０）で除することにより、第１パワースペクトルを正規化する。また、周波数分布発生部１５は、第２パワースペクトルＰｂａｓｅ（ｆ）を第２パワースペクトルの最大値Ｐｂａｓｅ（０）で除することにより、第２パワースペクトルを正規化する。
【００３３】
すなわち、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）は、Ｐｏｒｉｇ（ｆ）／Ｐｂａｓｅ（０）として計算される。正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）は、Ｐｂａｓｅ（ｆ）／Ｐｂａｓｅ（０）として計算される。周波数分布発生部１５は、正規化（ｎｏｒｍａｌｉｚｅ）係数αを計算する。正規化係数αは、例えば、以下に示す式のように、計算される。  
α＝Ｐｂａｓｅ（０）／Ｐｏｒｉｇ（０）・・・（２）
なお、第１領域および第２領域が脳室などのノイズ以外の信号成分を有する場合、第１パワースペクトルと第２パワースペクトルとにおいて、ＤＣ成分が発生する。このため、第１領域と第２領域とが同一部分（同位置）であれば、正規化係数αは、上式（２）と同一の式で計算される。また、第１領域および第２領域において、ｆ＝０における信号成分の有無にかかわらず正規化係数αを上式（２）で計算する場合、Ｐｂａｓｅ（０）とＰｏｒｉｇ（０）とは、ｆ＞０におけるＰｂａｓｅ（ｆ）とＰｏｒｉｇ（ｆ）とからそれぞれ外挿により生成する。
【００３４】
周波数分布発生部１５は、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）と正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）とに対して、強いスムージング（平滑化）を実行する。このスムージングにより、周波数分布発生部１５は、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）を平滑化した第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）と、正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）を平滑化した第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）とを発生する。なお、周波数分布発生部１５は、スムージングの代わりに、所定の関数の多項式、ガウシアン（Ｇａｕｓｓｉａｎ）関数、またはライシアン（Ｒｉｃｉａｎ）関数などを用いた曲線近似により、第１平滑スペクトルと第２平滑スペクトルとを発生してもよい。
【００３５】
また、周波数分布発生部１５は、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）と、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）との発生において、カットオフ周波数を設けてもよい。第１画像（第１補間画像）および第２画像（第２補間画像）がデジタル画像である場合、一般に高周波成分は、有限なサンプリング（ｓａｍｐｌｉｎｇ）周波数の制約などから、折り返し歪（エイリアジングによる歪み）を含みやすい。このため、周波数分布発生部１５は、ナイキスト（Ｎｙｑｕｉｓｔ）周波数ｆｎ以下にカットオフ周波数ｆｃを設定し（ｆｎ≧ｆｃ）、過剰な高周波（ｆ＞ｆｎ≧ｆｃ）をカットする。
【００３６】
周波数分布発生部１５は、通常、第１画像と第２画像とのうちサンプリング周波数の低い側にカットオフ周波数を合わせる。サンプリング周波数が低いのは過去画像であるｂａｓｅ画像（第２画像）である場合が多いため、周波数分布発生部１５は、第２画像のサンプリング周波数ｆｎ以下に、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）におけるカットオフ周波数を決定する。
【００３７】
図４は、例えば、第１領域および第２領域にそれぞれ対応する図３に示すノイズプロファイルの逆フーリエ変換を用いて発生され、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）、正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）、および第１、第２パワースペクトルをそれぞれ平滑化した第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）の一例を示す図である。図４に示すように、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）は、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）を平滑化したパワースペクトルである。第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）は、正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）を平滑化したパワースペクトルである。図４におけるｆｃは、カットオフ（ｃｕｔ  ｏｆｆ）周波数を示している。
【００３８】
周波数分布発生部１５は、第１画像に対応する第１周波数分布Ｓｏｒｉｇ（ｆ）を発生する。具体的には、周波数分布発生部１５は、第１補間画像または第１画像の画素値Ｉｏｒｉｇ（ｘ）に対して逆フーリエ変換することにより、第１周波数分布Ｓｏｒｉｇ（ｆ）を発生する。すなわち、上記記載によれば、周波数分布発生部１５は、例えば以下に示す式で、第１周波数分布Ｓｏｒｉｇ（ｆ）を計算する。
【００３９】
Ｓｏｒｉｇ（ｆ）＝ＩＦＴ［Ｉｏｒｉｇ（ｘ）］
周波数分布発生部１５は、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）と、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）と、第１周波数分布Ｓｏｒｉｇ（ｆ）とを、補正部１７に出力する。なお、周波数分布発生部１５は、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）と、正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）と、第１周波数分布Ｓｏｒｉｇ（ｆ）とを、補正部１７に出力してもよい。このとき、周波数分布発生部１５において、関数近似は実行されない。
【００４０】
補正部１７は、近似関数であるＰｎ＿ｅｓｔｏｒｉｇ（ｆ）とＰｎ＿ｅｓｔｂａｓｅ（ｆ）との比により、第１領域におけるノイズレベル（例えば、ノイズの標準偏差）を第２領域におけるノイズレベルに合わせるための関数ＰＮＲ（ｆ）を決定する。具体的には、補正部１７は、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）を、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）に変換する（合わせるまたは一致させる）変換関数ＰＮＲ（ｆ）を決定する。
【００４１】
すなわち、補正部１７は、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）を第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）で除することで、変換関数ＰＮＲ（ｆ）を決定する。変換関数ＰＮＲ（ｆ）は、例えば、ノイズのパワーの比（Ｐｏｗｅｒ  Ｎｏｉｓｅ  Ｒａｔｉｏ）に相当する。以下に示す式により、変換関数ＰＮＲ（ｆ）は決定される。なお、変換関数ＰＮＲ（ｆ）の高周波側は、カットオフ周波数ｆｃ以下に滑らかに低減される。  
ＰＮＲ（ｆ）＝Ｐｎ＿ｅｓｔｂａｓｅ（ｆ）／Ｐｎ＿ｅｓｔｏｒｉｇ（ｆ）・・・（３）  
上式（３）において、ゼロによる除算が生じないように工夫する必要がある。例えば、周波数分布発生部１５および図４では、カットオフ周波数ｆｃがナイキスト周波数ｆｎより大きいため、Ｐｎ＿ｅｓｔｏｒｉｇ（ｆ）は、非ゼロとなる。
【００４２】
なお、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）と、正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）とに対して関数近似が実行されない場合、補正部１７は、例えば、以下のようにして、変換関数ＰＮＲ（ｆ）を決定することができる。補正部１７は、Ｌ１、Ｌ２ノルム最適化などの最適化問題を解くためのアルゴリズム（最適化アルゴリズムなど）を適用した繰り返し演算（逐次近似などの演算）により、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）と、正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）とを比較しながら、変換関数ＰＮＲ（ｆ）を決定する。
【００４３】
図５は、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）を、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）に変換する変換関数の一例を示す図である。図５に示すｆｃは、カットオフ周波数である。図４に対応するＰＮＲ（ｆ）では、図５におけるＰＮＲ（ｆ）の最大値未満とＰＮＲ（ｆ）がαより大きい値との間において、ＰＮＲ（ｆ）＝１の直線が存在する。
【００４４】
補正部１７は、変換関数ＰＮＲ（ｆ）を用いて、第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。補正部１７は、補正した第１周波数分布Ｓｃｏｒ（ｆ）を画像発生部１９に出力する。具体的には、補正部１７は、第１周波数分布Ｓｏｒｉｇ（ｆ）に変換関数ＰＮＲ（ｆ）を乗ずることにより、第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。補正部１７は、例えば、以下に示す式に従って第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。
【００４５】
Ｓｃｏｒ（ｆ）＝Ｓｏｒｉｇ（ｆ）×ＰＮＲ（ｆ）・・・（４）
図６は、変換関数を用いて補正された第１周波数分布のパワースペクトルＰｎｏｒｉｇ＿ｃｏｒ（ｆ）の一例を示す図である。なお、図６におけるＰｎｏｒｉｇ＿ｃｏｒ（ｆ）は、第２パワースペクトルの最大値Ｐｂａｓｅ（０）で正規化されている。このため、図６に示すように、グラフの縦軸とＰｎｏｒｉｇ＿ｃｏｒ（ｆ）との交点は、１となる。
【００４６】
画像発生部１９は、補正された第１周波数分布Ｓｃｏｒ（ｆ）に対してフーリエ変換ＦＴを適用することで、第１画像（第１補間画像）を補正した画像Ｉｃｏｒ（ｘ）（以下、補正第１画像と呼ぶ）を発生する。具体的には、画像発生部１９は、以下に示す式に従って、補正第１画像を発生する。
【００４７】
Ｉｃｏｒ（ｘ）＝ＦＴ［Ｓｃｏｒ（ｆ）］・・・（５）  
上式（１）乃至（５）は説明を簡単にするため、１次元の式として示したが、補正対象の第１画像（または第２画像）は、２次元画像または３次元画像である。このため、上式（１）乃至（５）は画像の次元に対応する式となる。例えば、補正対象の次元が２次元である場合、本補正は、ｘ、ｙ方向ごとに実行される。また、補正対象の次元が３次元である場合、本補正は、ｘ、ｙ、ｚ方向ごとに実行されるため、高次元の方がより高精度な補正として実行される。
【００４８】
画像発生部１９は、発生した補正第１画像を、表示部２１に出力する。画像発生部１９は、補正第１画像を記憶部２７に出力してもよい。なお、画像発生部１９は、補正第１画像と第２画像とに基づいて、第２画像に対する補正第１画像の経時的変化を示す画像（以下、経時的変化画像と呼ぶ）を発生してもよい。画像発生部１９は、経時的変化画像を、表示部２１および記憶部２７に出力する。
【００４９】
表示部２１は、画像発生部１９により発生された各種医用画像、記憶部２７に記憶された医用画像を表示する。表示部２１は、入力操作に係る各種入力画面を表示する。
【００５０】
インターフェース部２３は、ネットワークＡ、図示していない外部記憶装置などに関するインターフェースである。インターフェース部２３は、図示していない医用画像撮影装置および医用画像保管装置に対して電子的通信回線、典型的にはローカルエリアネットワーク（Ｌｏｃａｌ  Ａｒｅａ  Ｎｅｔｗｏｒｋ：以下ＬＡＮと呼ぶ）などを介して接続される。医用画像撮影装置で発生されたボリュームデータおよび画像データ、または医用画像保管装置で保管されたボリュームデータおよび画像データは、インターフェース部２３を介して、記憶部２７に記憶される。また、本医用画像処理装置１で補正された医用画像のデータ、ボリュームデータなどは、ネットワークを介して図示していない医用画像保管装置などの他の装置に転送可能である。
【００５１】
入力部２５は、操作者などからの各種指示・命令・情報・選択・設定などを制御部２９に入力する。具体的には、入力部２５は、上記各種指示・命令・情報・選択・設定などを入力するためのトラックボール、スイッチボタン、マウス、マウス・ホイール、キーボード等の入力デバイスを有する。なお、入力デバイスは、表示部２１における表示画面を覆うタッチパネルでもよい。入力部２５は、入力デバイスを介して、補正対象の第１画像と、第１画像より過去に同種のモダリティにより収集され第２画像との選択操作を入力する。なお、入力部２５は、第１領域および第２領域の位置、形状および大きさを入力することも可能である。
【００５２】
記憶部２７は、インターフェース部２３を介して転送された第１画像および第２画像、周波数分布発生部１５により発生された各種周波数分布、パワースペクトル、補正部１７により決定された変換関数、画像発生部１９により発生された各種医用画像などを記憶する。なお、記憶部２７は、領域設定部１３で用いられる第１領域および第２領域の形状および大きさを記憶してもよい。記憶部２７は、制御部２９により用いられる各種制御プログラムを記憶する。また、記憶部２７は、後述する周波数分布補正処理に関する補正プログラム、例えば、補間処理、第１、第２領域の設定、周波数分布および各種パワースペクトルの発生、パワースペクトルの正規化および平滑化、変換関数の決定、周波数分布の補正、補正画像の発生等に関する少なくとも一つ以上のプログラムを記憶する。
【００５３】
制御部２９は、図示していない中央処理装置（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ：以下、ＣＰＵと呼ぶ）、メモリ等を有する。制御部２９は、本医用画像処理装置１の中枢として機能する。具体的には、制御部２９は、入力部２５を介して入力された入力操作に応じて、記憶部２７から各種プログラムなどを読み出して自身のメモリ上に展開する。制御部２９は、展開したプログラムを実行することにより、本医用画像処理装置１の各部を制御する。また、例えば、制御部２９は、記憶部２７から補正プログラムを読み出して実行することにより、補間処理部１１、領域設定部１３、周波数分布発生部１５、補正部１７、画像発生部１９などの各部を制御する。
【００５４】
なお、本医用画像処理装置１は、各種モダリティ、例えば、磁気共鳴イメージング（Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  Ｉｍａｇｉｎｇ：ＭＲＩ）装置、Ｘ線コンピュータ断層撮影（Ｃｏｍｐｕｔｅｄ  Ｔｏｍｏｇｒａｐｈｙ：ＣＴ）装置、単光子放出コンピュータ断層撮影（Ｓｉｎｇｌｅ  Ｐｈｏｔｏｎ  Ｅｍｉｓｓｉｏｎ  Ｃｏｍｐｕｔｅｄ  Ｔｏｍｏｇｒａｐｈｙ：ＳＰＥＣＴ）装置、陽電子放出コンピュータ断層撮影（Ｐｏｓｉｔｉｒｏｎ  Ｅｍｉｓｓｉｏｎ  Ｃｏｍｐｕｔｅｄ  Ｔｏｍｏｇｒａｐｈｙ：ＰＥＴ）装置、Ｘ線診断装置などのコンソール（操作卓等）に搭載されてもよい。また、本医用画像処理装置１は、ワークステーションとして実現することも可能である。また、医用画像処理装置１は、例えば、シンクライアント（Ｔｈｉｎ  Ｃｌｉｅｎｔ）およびクラウドシステムにおけるサーバとして実現されてもよい。
【００５５】
（周波数分布補正機能）  
周波数分布補正機能とは、補正対象の医用画像（第１画像）のノイズレベルと、基準となる医用画像（第２画像）のノイズレベルとを一致させるように、補正対象の医用画像の周波数分布を補正する機能である。以下、周波数分布補正機能に関する処理（以下、周波数分布補正処理と呼ぶ）について説明する。
【００５６】
図７は、周波数分布補正処理に関する手順の一例を示すフローチャートである。  
入力部２５を介した操作者の指示により、選択された第１画像と第２画像とが、記憶部２７から読み出される。なお、選択された第１画像および第２画像は、インターフェース部２３とネットワークとを介して、図示していない医用画像保管装置、各種モダリティ等から読み出されてもよい。
【００５７】
読み出された第１画像における複数の画素値に基づいて、空間的に均一な第１領域が設定される（ステップＳａ１）。読み出された第２画像における複数の画素値に基づいて、空間的に均一な第２領域が設定される（ステップＳａ２）。
【００５８】
第１領域における複数の画素値に基づいて、第１領域周波数分布に対応する第１パワースペクトルが発生される（ステップＳａ３）。第２領域における複数の画素値に基づいて、第２領域周波数分布に対応する第２パワースペクトルが発生される（ステップＳａ３）。第２パワースペクトルの最大値を用いて、第１、第２パワースペクトル各々が正規化される（ステップＳａ４）。
【００５９】
正規化された第１パワースペクトルを平滑化することにより、第１平滑スペクトルが発生される（ステップＳａ５）。正規化された第２パワースペクトルを平滑化することにより、第２平滑スペクトルが発生される（ステップＳａ５）。第１平滑スペクトルを第２平滑スペクトルに変換する変換関数が決定される（ステップＳａ６）。
【００６０】
第１画像（第１補間画像）の全域における複数の画素値に基づいて、第１周波数分布が発生される（ステップＳａ７）。変換関数を用いて、第１周波数分布が補正される（ステップＳａ８）。補正された第１周波数分布に基づいて、第１画像を補正した補正第１画像が発生される（ステップＳａ９）。
【００６１】
（第１の変形例）  
第１の実施形態との相違は、基準として第１パワースペクトルを用いることにある。  
周波数分布発生部１５は、第１パワースペクトルＰｏｒｉｇ（ｆ）を第１パワースペクトルの最大値Ｐｏｒｉｇ（０）で除することにより、第１パワースペクトルを正規化する。また、周波数分布発生部１５は、第２パワースペクトルＰｂａｓｅ（ｆ）を第１パワースペクトルの最大値Ｐｏｒｉｇ（０）で除することにより、第２パワースペクトルを正規化する。
【００６２】
すなわち、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）は、Ｐｏｒｉｇ（ｆ）／Ｐｏｒｉｇ（０）として計算される。正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）は、Ｐｂａｓｅ（ｆ）／Ｐｏｒｉｇ（０）として計算される。周波数分布発生部１５は、正規化（ｎｏｒｍａｌｉｚｅ）係数αを計算する。正規化係数αは、例えば、以下に示す式のように、計算される。
【００６３】
α＝Ｐｏｒｉｇ（０）／Ｐｂａｓｅ（０）
周波数分布発生部１５は、第２画像に対応する第２周波数分布Ｓｂａｓｅ（ｆ）を発生する。具体的には、周波数分布発生部１５は、第２補間画像または第２画像の画素値Ｉｂａｓｅ（ｘ）に対して逆フーリエ変換することにより、第２周波数分布Ｓｂａｓｅ（ｆ）を発生する。すなわち、上記記載によれば、周波数分布発生部１５は、例えば以下に示す式で、第２周波数分布Ｓｂａｓｅ（ｆ）を計算する。
【００６４】
Ｓｂａｓｅ（ｆ）＝ＩＦＴ［Ｉｂａｓｅ（ｘ）］
補正部１７は、近似関数であるＰｎ＿ｅｓｔｏｒｉｇ（ｆ）とＰｎ＿ｅｓｔｂａｓｅ（ｆ）との比により、第２領域におけるノイズレベルを第１領域におけるノイズレベルに合わせるための関数ＰＮＲ（ｆ）を決定する。具体的には、補正部１７は、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）を第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）に変換する（合わせるまたは一致させる）変換関数ＰＮＲ（ｆ）を決定する。
【００６５】
すなわち、補正部１７は、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）を第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）で除することで、変換関数ＰＮＲ（ｆ）を決定する。変換関数ＰＮＲ（ｆ）は、例えば、ノイズのパワーの比（Ｐｏｗｅｒ  Ｎｏｉｓｅ  Ｒａｔｉｏ）に相当する。以下に示す式により、変換関数ＰＮＲ（ｆ）は決定される。  
ＰＮＲ（ｆ）＝Ｐｎ＿ｅｓｔｏｒｉｇ（ｆ）／Ｐｎ＿ｅｓｔｂａｓｅ（ｆ）  
上式において、ゼロによる除算が生じないように、例えば、周波数分布発生部１５において、カットオフ周波数ｆｃを決定してもよい。なお、変換関数は、第１画像および第２画像の発生時における補間の周波数分布、補間回数、補間アルゴリズム、補間サイズ、適応されたフィルターの種類及び形状等の各種情報が既知であれば、補間の周波数分布、補間回数、適応されたフィルター等を用いて変換関数を決定してもよい。
【００６６】
補正部１７は、変換関数ＰＮＲ（ｆ）を用いて、第２周波数分布Ｓｂａｓｅ（ｆ）を補正する。補正部１７は、補正した第２周波数分布Ｓｃｏｒ２（ｆ）を画像発生部１９に出力する。具体的には、補正部１７は、第２周波数分布Ｓｂａｓｅ（ｆ）に変換関数ＰＮＲ（ｆ）を乗ずることにより、第２周波数分布Ｓｂａｓｅ（ｆ）を補正する。補正部１７は、例えば、以下に示す式に従って第２周波数分布Ｓｂａｓｅ（ｆ）を補正する。
【００６７】
Ｓｃｏｒ２（ｆ）＝Ｓｂａｓｅ（ｆ）×ＰＮＲ（ｆ）
画像発生部１９は、補正された第２周波数分布Ｓｃｏｒ２（ｆ）に対してフーリエ変換ＦＴを適用することで、第２画像（第２補間画像）を補正した画像Ｉｃｏｒ２（ｘ）（以下、補正第２画像と呼ぶ）を発生する。具体的には、画像発生部１９は、以下に示す式に従って、補正第２画像を発生する。
【００６８】
Ｉｃｏｒ２（ｘ）＝ＦＴ［Ｓｃｏｒ２（ｆ）］
画像発生部１９は、発生した補正第２画像を、表示部２１に出力する。画像発生部１９は、補正第２画像を記憶部２７に出力してもよい。なお、画像発生部１９は、補正第２画像と第１画像とに基づいて、第１画像に対する補正第２画像の経時的変化を示す画像を発生してもよい。
【００６９】
（第２の変形例）  
第１の実施形態および第１の変形例との相違は、基準として所定の周波数分布に対応するパワースペクトルを用いることにある。  
記憶部２７は、所定の周波数分布に対応するパワースペクトル（以下、所定のパワースペクトルと呼ぶ）Ｐｐｒｅｄｅｔ（ｆ）を記憶する。所定の周波数分布とは、例えば、ホワイトノイズに対応する周波数分布、または固定パターンノイズに対応する周波数分布などである。また、記憶部２７は、正規化された所定のパワースペクトル（以下、正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）を記憶する。記憶された所定のパワースペクトルＰｐｒｅｄｅｔ（ｆ）は、周波数分布発生部１５により、記憶部２７から周波数分布発生部１５に読み出される。記憶された正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）は、補正部１７により、記憶部２７から補正部１７に読み出される。
【００７０】
周波数分布発生部１５は、第１周波数分布Ｓｏｒｉｇ（ｆ）と、第２周波数分布Ｓｂａｓｅ（ｆ）を発生する。周波数分布発生部１５は、第１パワースペクトルＰｏｒｉｇ（ｆ）を所定のパワースペクトルの最大値Ｐｐｒｅｄｅｔ（０）で除することにより、第１パワースペクトルを正規化する。周波数分布発生部１５は、第２パワースペクトルＰｂａｓｅ（ｆ）を所定のパワースペクトルの最大値Ｐｐｒｅｄｅｔ（０）で除することにより、第２パワースペクトルを正規化する。正規化された第１パワースペクトルと、正規化された第２パワースペクトルとは、補正部１７に出力される。
【００７１】
具体的には、すなわち、正規化された第１パワースペクトルＰｎｏｒｉｇ（ｆ）は、Ｐｏｒｉｇ（ｆ）／Ｐｐｒｅｄｅｔ（０）として計算される。正規化された第２パワースペクトルＰｎｂａｓｅ（ｆ）は、Ｐｂａｓｅ（ｆ）／Ｐｐｒｅｄｅｔ（０）として計算される。周波数分布発生部１５は、第１パワースペクトルに対応する正規化（ｎｏｒｍａｌｉｚｅ）係数α１を計算する。周波数分布発生部１５は、第２パワースペクトルに対応する正規化（ｎｏｒｍａｌｉｚｅ）係数α２を計算する。正規化係数α１は、例えば、以下に示す式のように、計算される。
【００７２】
α１＝Ｐｐｒｅｄｅｔ（０）／Ｐｏｒｉｇ（０）
また、正規化係数α２は、例えば、以下に示す式のように、計算される。
【００７３】
α２＝Ｐｐｒｅｄｅｔ（０）／Ｐｂａｓｅ（０）
補正部１７は、正規化された第１パワースペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）と正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）との比により、第１領域におけるノイズレベルを、所定の周波数分布に関するノイズレベル（以下、所定のノイズレベルと呼ぶ）に合わせるための第１変換関数ＰＮＲ１（ｆ）を決定する。補正部１７は、正規化された第２パワースペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）と正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）との比により、第２領域におけるノイズレベルを、所定のノイズレベルに合わせるための第２変換関数ＰＮＲ２（ｆ）を決定する。
【００７４】
具体的には、補正部１７は、第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）を正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）に変換する（合わせるまたは一致させる）第１変換関数ＰＮＲ（ｆ）を決定する。すなわち、補正部１７は、正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）を第１平滑スペクトルＰｎ＿ｅｓｔｏｒｉｇ（ｆ）で除することで、第１変換関数ＰＮＲ１（ｆ）を決定する。補正部１７は、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）を正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）に変換する（合わせるまたは一致させる）第２変換関数ＰＮＲ（ｆ）を決定する。すなわち、補正部１７は、正規化パワースペクトルＰｎｐｒｅｄｅｔ（ｆ）を第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）で除することで、第２変換関数ＰＮＲ（ｆ）を決定する。
【００７５】
第１変換関数ＰＮＲ１（ｆ）および第２変換関数ＰＮＲ２（ｆ）は、例えば、ノイズのパワーの比（Ｐｏｗｅｒ  Ｎｏｉｓｅ  Ｒａｔｉｏ）に相当する。以下に示す２式により、第１変換関数ＰＮＲ１（ｆ）と第２変換関数ＰＮＲ２（ｆ）とは、それぞれ決定される。
【００７６】
ＰＮＲ１（ｆ）＝Ｐｎｐｒｅｄｅｔ（ｆ）／Ｐｎ＿ｅｓｔｏｒｉｇ（ｆ）
ＰＮＲ２（ｆ）＝Ｐｎｐｒｅｄｅｔ（ｆ）／Ｐｎ＿ｅｓｔｂａｓｅ（ｆ）
補正部１７は、第１変換関数ＰＮＲ１（ｆ）を用いて、第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。補正部１７は、第２変換関数ＰＮＲ２（ｆ）を用いて、第２周波数分布Ｓｂａｓｅ（ｆ）を補正する。補正部１７は、補正した第１周波数分布Ｓｃｏｒ１（ｆ）と、補正した第２周波数分布Ｓｃｏｒ２（ｆ）とを画像発生部１９に出力する。
【００７７】
具体的には、補正部１７は、第１周波数分布Ｓｏｒｉｇ（ｆ）に第１変換関数ＰＮＲ１（ｆ）を乗ずることにより、第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。補正部１７は、例えば、以下に示す式に従って第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。
【００７８】
Ｓｃｏｒ１（ｆ）＝Ｓｏｒｉｇ（ｆ）×ＰＮＲ１（ｆ）
また、補正部１７は、第２周波数分布Ｓｂａｓｅ（ｆ）に第２変換関数ＰＮＲ２（ｆ）を乗ずることにより、第２周波数分布Ｓｂａｓｅ（ｆ）を補正する。補正部１７は、例えば、以下に示す式に従って第２周波数分布Ｓｂａｓｅ（ｆ）を補正する。
【００７９】
Ｓｃｏｒ２（ｆ）＝Ｓｂａｓｅ（ｆ）×ＰＮＲ２（ｆ）
画像発生部１９は、補正された第１周波数分布Ｓｃｏｒ１（ｆ）に対してフーリエ変換ＦＴを適用することで、第１画像（第１補間画像）を補正した補正第１画像Ｉｃｏｒ１（ｘ）を発生する。具体的には、画像発生部１９は、以下に示す式に従って、補正第１画像を発生する。
【００８０】
Ｉｃｏｒ１（ｘ）＝ＦＴ［Ｓｃｏｒ１（ｆ）］
画像発生部１９は、補正された第２周波数分布Ｓｃｏｒ２（ｆ）に対してフーリエ変換ＦＴを適用することで、第２画像（第２補間画像）を補正した補正第２画像Ｉｃｏｒ２（ｘ）を発生する。具体的には、画像発生部１９は、以下に示す式に従って、補正第２画像を発生する。
【００８１】
Ｉｃｏｒ２（ｘ）＝ＦＴ［Ｓｃｏｒ２（ｆ）］
画像発生部１９は、補正第１画像と補正第２画像とを、表示部２１に出力する。画像発生部１９は、補正第１画像と補正第２画像とを記憶部２７に出力してもよい。なお、画像発生部１９は、補正第１画像と補正第２画像とに基づいて、補正第２画像に対する補正第１画像の経時的変化を示す画像を発生してもよい。  
以上に述べた構成によれば、以下の効果を得ることができる。  
本実施形態に係る医用画像処理装置１によれば、非構造的であって空間的に均一な領域である第１領域と第２領域とを、第１画像と第２画像とにそれぞれ設定し、設定された第１領域における第１パワースペクトルと第２パワースペクトルとに基づいて第２画像のノイズレベルに第１画像のノイズレベルを合わせた補正第１画像を発生することができる。特に、本実施形態に係る周波数分布補正処理は、画像再構成手法やフィルターが線形アルゴリズムの場合に有効である。
【００８２】
また、本実施形態に係る第１の変形例によれば、第１パワースペクトルと第２パワースペクトルとに基づいて第１画像のノイズレベルに第２画像のノイズレベルを合わせた補正第２画像を発生することができる。
【００８３】
また、本実施形態に係る第２の変形例によれば、第１パワースペクトルと第２パワースペクトルと所定の周波数分布に関するパワースペクトルに基づいて、所定のノイズレベルに第１画像のノイズレベルと第２画像のノイズレベルとをそれぞれ合わせた補正第１画像と補正第２画像とを発生することができる。
【００８４】
すなわち、本医用画像処理装置１によれば、第１画像および第２画像の収集装置（モダリティ）、および収集条件のパラメータなどの装置由来の画像の変化の影響を除外した補正画像を、画像以外の情報がなくても発生することができる。これにより、意図的に医用画像をぼかすことなく、被検体の経時的な解剖学的・機能的変化、すなわち、本質的な経時変化を表示することができる。
【００８５】
以上のことから、本医用画像処理装置１によれば、ＶＢＭなどの分野を含んで広く、経時的な画像変化を追跡する場合において、撮像装置由来の因子を除外し患者本来の因子による経時変化を追跡することができる。これにより、操作者は、被検体に関する医用画像において、本来の経時的変化を、より詳細に追跡することが可能となる。
【００８６】
（第２の実施形態）  
第１の実施形態との相違は、第１領域および第２領域が、経時的に変化の小さい領域または経時的に不変な領域（以下、経時的不変領域と呼ぶ）であって、ノイズ以外の信号成分を有することにある。  
図８は、本実施形態に係る医用画像処理装置１の構成を示す構成図である。図１との相違は、エッジ検出部１０を有することにある。
【００８７】
エッジ検出部１０は、第１画像（または第１補間画像）において、所定のエッジ検出処理により、複数の第１エッジを検出する。エッジ検出部１０は、第２画像（または第２補間画像）において、所定のエッジ検出処理により、複数の第２エッジを検出する。
【００８８】
記憶部２７は、所定の閾値を記憶する、所定の閾値とは、第１画像および第２画像において、経時的不変領域を規定するエッジに対応する画素値である。経時的不変領域は、二つのエッジで挟まれた略一定（平坦）の画素値を有する領域であって、例えば、第１画像および第２画像において、被検体とともに撮影された所定のファントムを示す領域、骨の形状を示す領域（例えば、脳画像において頭蓋骨に対応する領域）、主として水を示す領域（例えば、脳画像におけるＣＳＦなど）などである。以下、説明を簡単にするために、経時的不変領域は、ファントムに対応する領域またはＣＳＦに対応する領域であるものとする。
【００８９】
領域設定部１３は、第１領域における対向する２辺が第１エッジに直交し、かつ所定の閾値で挟まれた画素値範囲を、第１領域として、第１補間画像に設定する。領域設定部１３は、第２領域における対向する２辺が第２エッジに直交し、かつ所定の閾値で挟まれた画素値範囲を、第２領域として、第２補間画像に設定する。第１領域と第２領域とは、画素値範囲が同じ領域（例えば、ファントム、骨、ＣＳＦ、水など）が好ましい。なお、第１領域および第２領域としてファントムを用いる場合、安定した周波数分布補正を実行するため、第１画像および第２画像各々の収集時と同時に、またはできるだけ小さい時間差で、ファントム画像を収集した方が好ましい。なお、第１領域および第２領域は、経時的変化が小さい骨、水抽出などにより、設定されてもよい。また、第１領域および第２領域は、入力部２５を介した操作者の指示により設定されてもよい。
【００９０】
図９は、第１補間画像および第２補間画像において、第１領域または第２領域の設定の一例を示す図である。図９におけるｆはファントムを示している。図９における線分ｄｄ’は、ファントムｆに対して設定される第１領域または第２領域の一例を示している。図９における線分ｅｅ’は、頭蓋骨に対応する領域およびＣＳＦに対応する領域に対して設定される第１領域または第２領域の一例を示している。
【００９１】
図１０は、図９における線分ｅｅ’上において一次元に配列した複数の画素にそれぞれ対応する複数の画素値による画素値の変化（実信号プロファイル）の一例を示す図である。図１０の横軸は、画素の位置ｘを示している。図１０の縦軸は、画像Ｉにおける画素の位置ｘの画素値Ｉ（ｘ）を示している。図１０において、第１領域および第２領域は、例えば、線分ｆｆ’で挟まれた領域として設定される。
【００９２】
図１１は、図９における線分ｄｄ’または図１０における抽出部分（ｆｆ’）上の画素値の変化（経時不変プロファイル）の一例を示す図である。図１１の横軸は、画素の位置ｘを示している。図１１の縦軸は、画像Ｉにおける画素の位置ｘの画素値Ｉ（ｘ）を示している。図１１におけるＩｏｒｉｇ（ｘ）は、第１画像において線分ｄｄ’または線分ｆｆ’に位置する複数の画素にそれぞれ対応する複数の画素値のプロファイルを示している。図１１におけるＩｂａｓｅ（ｘ）は、第２画像において線分ｄｄ’または線分ｆｆ’に位置する複数の画素にそれぞれ対応する複数の画素値のプロファイルを示している。図９乃至図１１において、第１領域および第２領域は１次元プロファイルとして説明しているが、２次元または３次元プロファイルであってもよい。
【００９３】
周波数分布発生部１５は、第１領域に対応する第１領域周波数分布Ｓ１ｏｒｉｇ（ｆ）の絶対値｜Ｓ１ｏｒｉｇ（ｆ）｜を発生する。具体的には、周波数分布発生部１５は、第１領域における複数の画素値Ｉ１ｏｒｉｇ（ｘ）に対して逆フーリエ変換を実行して絶対値をとることにより、第１領域周波数分布の絶対値｜Ｓ１ｏｒｉｇ（ｆ）｜を発生する。すなわち、上記記載によれば、周波数分布発生部１５は、例えば以下に示す式で、第１領域周波数分布の絶対値｜Ｓ１ｏｒｉｇ（ｆ）｜を計算する。
【００９４】
｜Ｓ１ｏｒｉｇ（ｆ）｜＝｜ＩＦＴ［Ｉ１ｏｒｉｇ（ｘ）］｜
周波数分布発生部１５は、第２領域に対応する第２領域周波数分布Ｓ２ｂａｓｅ（ｆ）の絶対値｜Ｓ２ｂａｓｅ（ｆ）｜を発生する。具体的には、周波数分布発生部１５は、第２領域における複数の画素値Ｉ２ｂａｓｅ（ｘ）に対して逆フーリエ変換を実行して絶対値をとることにより、第２領域周波数分布の絶対値｜Ｓ２ｂａｓｅ（ｆ）｜を発生する。すなわち、上記記載によれば、周波数分布発生部１５は、例えば以下に示す式で、第２領域周波数分布の絶対値の｜Ｓ２ｂａｓｅ（ｆ）｜を計算する。
【００９５】
｜Ｓ２ｂａｓｅ（ｆ）｜＝｜ＩＦＴ［Ｉ２ｂａｓｅ（ｘ）］｜
周波数分布発生部１５は、第１領域周波数分布の絶対値を、第２領域周波数分布の絶対値の最大値｜Ｓ２ｂａｓｅ（０）｜で除することにより、第１領域周波数分布の絶対値を正規化する。また、周波数分布発生部１５は、第２領域周波数分布の絶対値を第２領域周波数分布の絶対値の最大値｜Ｓ２ｂａｓｅ（０）｜で除することにより、第２領域周波数分布の絶対値を正規化する。
【００９６】
すなわち、正規化された第１領域周波数分布の絶対値（以下、第１正規化周波数分布と呼ぶ）｜Ｓｎ１ｏｒｉｇ（ｆ）｜は、｜Ｓ１ｏｒｉｇ（ｆ）｜／｜Ｓ２ｂａｓｅ（０）｜として計算される。正規化された第２領域周波数分布の絶対値（以下、第２正規化周波数分布と呼ぶ）｜Ｓｎ２ｂａｓｅ（ｆ）｜は、｜Ｓ２ｂａｓｅ（ｆ）｜／｜Ｓ２ｂａｓｅ（０）｜として計算される。周波数分布発生部１５は、正規化（ｎｏｒｍａｌｉｚｅ）係数αを計算する。正規化係数αは、例えば、以下に示す式のように、計算される。  
α＝｜Ｓ２ｂａｓｅ（０）｜／｜Ｓ１ｏｒｉｇ（０）｜  
周波数分布発生部１５は、第１正規化周波数分布｜Ｓｎ１ｏｒｉｇ（ｆ）｜と第２正規化周波数分布｜Ｓｎ２ｂａｓｅ（ｆ）｜とを、補正部１７に出力する。
【００９７】
図１２は、経時不変プロファイルの逆フーリエ変換の絶対値の正規化により発生された第１正規化周波数分布｜Ｓｎ１ｏｒｉｇ（ｆ）｜と第２正規化周波数分布｜Ｓｎ２ｂａｓｅ（ｆ）｜との一例を示す図である。図１２における第１正規化周波数分布｜Ｓｎ１ｏｒｉｇ（ｆ）｜と第２正規化周波数分布｜Ｓｎ２ｂａｓｅ（ｆ）｜とにおいて、ノイズの影響は実施形態１に比べて比較的小さいため、実施形態１におけるスムージングや関数近似のコストは、実施形態１に比べて小さい。説明を簡単にするため、本実施形態ではスムージングまたは関数近似の処理について言及していないが、｜Ｓ１ｏｒｉｇ（ｆ）｜と｜Ｓ２ｂａｓｅ（ｆ）｜とに対して、スムージングまたは関数近似が適宜実行されてもよい。
【００９８】
補正部１７は、第１正規化周波数分布｜Ｓｎ１ｏｒｉｇ（ｆ）｜と第２正規化周波数分布｜Ｓｎ２ｂａｓｅ（ｆ）｜との比により、第１領域におけるノイズレベルを第２領域におけるノイズレベルに合わせるための関数Ｆｃｏｒ（ｆ）を決定する。具体的には、補正部１７は、第１正規化周波数分布｜Ｓｎ１ｏｒｉｇ（ｆ）｜を、第２平滑スペクトルＰｎ＿ｅｓｔｂａｓｅ（ｆ）に変換する（合わせるまたは一致させる）変換関数Ｆｃｏｒ（ｆ）を決定する。
【００９９】
すなわち、補正部１７は、第２正規化周波数分布｜Ｓｎ２ｂａｓｅ（ｆ）｜を第１正規化周波数分布｜Ｓｎ１ｏｒｉｇ（ｆ）｜で除することで、変換関数Ｆｃｏｒ（ｆ）を決定する。具体的には、以下に示す式により、変換関数Ｆｃｏｒ（ｆ）は決定される。なお、変換関数Ｆｃｏｒ（ｆ）の高周波側は、カットオフ周波数ｆｃ以下に滑らかに低減される。
【０１００】
Ｆｃｏｒ（ｆ）＝｜Ｓｎ２ｂａｓｅ（ｆ）｜／｜Ｓｎ１ｏｒｉｇ（ｆ）｜
上式において、ゼロによる除算が生じないように工夫する必要がある。例えば、周波数分布発生部１５および図１２乃至図１３では、カットオフ周波数ｆｃがナイキスト周波数ｆｎより大きいため、｜Ｓｎ１ｏｒｉｇ（ｆ）｜は、非ゼロとなる。
【０１０１】
図１３は、第１正規化周波数分布を第２正規化周波数分布に変換する変換関数Ｆｃｏｒ（ｆ）の一例を示す図である。図１３に示すｆｃは、カットオフ周波数である。
【０１０２】
補正部１７は、変換関数Ｆｃｏｒ（ｆ）を用いて、第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。補正部１７は、補正した第１周波数分布Ｓｃｏｒ（ｆ）を画像発生部１９に出力する。具体的には、補正部１７は、第１周波数分布Ｓｏｒｉｇ（ｆ）に変換関数Ｆｃｏｒ（ｆ）を乗ずることにより、第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。補正部１７は、例えば、以下に示す式に従って第１周波数分布Ｓｏｒｉｇ（ｆ）を補正する。
【０１０３】
Ｓｃｏｒ（ｆ）＝Ｓｏｒｉｇ（ｆ）×Ｆｃｏｒ（ｆ）
図１４は、変換関数Ｆｃｏｒ（ｆ）を用いて補正された第１周波数分布の絶対値｜Ｓｃｏｒ（ｆ）｜の一例を示す図である。図１４に示すｆｃは、カットオフ周波数である。
【０１０４】
（周波数分布補正機能）  
図１５は、周波数分布補正処理に関する手順の一例を示すフローチャートである。  
第１画像において、複数の第１エッジが検出される（ステップＳｂ１）。第２画像において、複数の第２エッジが検出される（ステップＳｂ２）。第１領域の対辺が第１エッジに直交し、経時的に変化の小さい領域が、第１画像に第１領域として設定される（ステップＳｂ３）。第２領域の対辺が第２エッジに直交し、経時的に変化の小さい領域が、第２画像に第２領域として設定される（ステップＳｂ３）。
【０１０５】
第１領域における複数の画素値に基づいて、第１領域周波数分布の絶対値が発生される（ステップＳｂ４）。第２領域における複数の画素値に基づいて、第２領域周波数分布の絶対値が発生される（ステップＳｂ４）。第２領域周波数分布の絶対値の最大値を用いて、第１、第２領域周波数分布の絶対値各々が正規化される（ステップＳｂ５）。
【０１０６】
第１正規化周波数分布を第２正規化周波数分布に変換する変換関数が決定される（ステップＳｂ７）。第１画像の全域における複数の画素値に基づいて、第１周波数分布が発生される（ステップＳｂ８）。変換関数を用いて、第１周波数分布が補正される（ステップＳｂ９）。補正された第１周波数分布に基づいて、第１画像を補正した補正第１画像が発生される（ステップＳｂ１０）。
【０１０７】
以上に述べた構成によれば、以下の効果を得ることができる。  
本実施形態に係る医用画像処理装置１によれば、非構造的であって経時的に不変な領域である第１領域と第２領域とを、第１画像と第２画像とにそれぞれ設定し、第１領域における第１領域周波数分布と第２領域における第２領域周波数分布とに基づいて第２画像のノイズレベルに第１画像のノイズレベルを合わせた補正第１画像を発生することができる。
【０１０８】
本実施形態による周波数分布補正処理は、ノイズの特性と画像構造部分の特性との類推が困難となる場合、具体的には、例えば構造の複雑な部分を除外して、画素値が平坦な部分だけを強く平滑化するような非線形のアルゴリズムを使用している場合において、第１の実施形態に係る周波数分布補正処理に比べて有効となる。
【０１０９】
なお、第１の実施形態に係る周波数分布補正処理と第２の実施形態に係る周波数分布補正処理とは、組み合わせて用いることができる。例えば、通常は、対象画像（第１画像）に対して過去に取得されたｂａｓｅ画像（第２画像）が従来の医用画像発生装置で生成されている場合、線形のアルゴリズムを用いて画像生成が実行されている可能性が大きい。また対象画像となりうる最近の装置による画像では、空気部分をマスクしている場合も多いため、第１の実施形態に係る周波数分布補正処理ではパワースペクトルの取得は困難となる。このような場合、第１の実施形態に係る周波数分布補正処理を用いて、ｂａｓｅ画像（第２画像）からパワースペクトルを取得し、第２の実施形態に係る周波数分布補正処理を用いて、対象画像（第１画像）からパワースペクトルを取得することで、第１周波数分布を補正することが可能である。
【０１１０】
以上のことから、第１、第２の実施形態に係る医用画像処理装置１によれば、経時的な変化を過去画像と比較する場合に、装置依存の変動を除外するため、すなわち経時変化における撮像装置に由来する因子を除外するために、周波数分布を補正することができる。これにより、本医用画像処理装置１によれば、操作者は、患者本来の因子による解剖学的な経時変化の追跡を、高精度で実行することができる。
【０１１１】
（第３の実施形態）  
第１の実施形態と第２実施形態との相違は、第１画像および第２画像のパワースペクトル特性に関係するフィルター関数や補間処理に関する補間関数と、補間回数とに関する各種情報が既知であるならば、この各種情報を用いて、補正部１７で用いられる変換関数Ｆｃｏｒ（ｆ）、第１変換関数ＰＮＲ１（ｆ）、または第２変換関数ＰＮＲ２（ｆ）などの補正関数を生成して、補正することにある。
【０１１２】
記憶部２７は、第１画像および第２画像のパワースペクトル特性に関係するフィルター関数や補間処理に関する補間関数と、補間回数とに関する各種情報を記憶する。
【０１１３】
補正部１７は、記憶部２７から上記各種情報を読み出して、上記各種情報に基づいて、変換関数Ｆｃｏｒ（ｆ）、第１変換関数ＰＮＲ１（ｆ）、または第２変換関数ＰＮＲ２（ｆ）などの補正関数を決定する。例えば、補正部１７は、第１領域周波数分布および第２領域周波数分布に関するフィルター関数と補間関数と補間回数等とに基づいて補正関数を決定する。次いで、補正部１７は、決定した補正関数を用いて、第１周波数分布と第２周波数分布とのうち少なくとも一方等を補正する。
【０１１４】
以上のことから、第３の実施形態に係る医用画像処理装置１によれば、第１、第２の実施形態に係る効果に加えて、以下に示す効果を更に有する。第３の実施形態に係る医用画像処理装置１は、第１、第２の実施形態に比べて、周波数分布発生部１５等における各種パワースペクトルの特性に係る各種計算（パワースペクトル特性測定）は、不要となり、より高精度に、周波数分布、パワースペクトル等を補正することができる。
【０１１５】
加えて、実施形態に係る各機能は、当該処理を実行するプログラムをワークステーション等のコンピュータにインストールし、これらをメモリ上で展開することによっても実現することができる。このとき、コンピュータに当該手法を実行させることのできるプログラムは、磁気ディスク（フロッピー（登録商標）ディスク、ハードディスクなど）、光ディスク（ＣＤ−ＲＯＭ、ＤＶＤなど）、半導体メモリなどの記憶媒体に格納して頒布することも可能である。
【０１１６】
以上、本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１１７】
１…医用画像処理装置、１０…エッジ検出部、１１…補間処理部、１３…領域設定部、１５…周波数分布発生部、１７…補正部、１９…画像発生部、２１…表示部、２３…インターフェース部、２５…入力部、２７…記憶部、２９…制御部。
【図１】
0006491011.tif 000002
【図２】
0006491011.tif 000003
【図３】
0006491011.tif 000004
【図４】
0006491011.tif 000005
【図５】
0006491011.tif 000006
【図６】
0006491011.tif 000007
【図７】
0006491011.tif 000008
【図８】
0006491011.tif 000009
【図９】
0006491011.tif 000010
【図１０】
0006491011.tif 000011
【図１１】
0006491011.tif 000012
【図１２】
0006491011.tif 000013
【図１３】
0006491011.tif 000014
【図１４】
0006491011.tif 000015
【図１５】
0006491011.tif 000016