(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-05-29
(45)【発行日】2024-06-06
(54)【発明の名称】画像処理装置及び画像処理方法。
(51)【国際特許分類】
A61B 8/14 20060101AFI20240530BHJP
【FI】
A61B8/14
(21)【出願番号】P 2020096755
(22)【出願日】2020-06-03
【審査請求日】2023-05-18
(73)【特許権者】
【識別番号】320011683
【氏名又は名称】富士フイルムヘルスケア株式会社
(74)【代理人】
【識別番号】110000350
【氏名又は名称】ポレール弁理士法人
(72)【発明者】
【氏名】吉川 秀樹
(72)【発明者】
【氏名】川畑 健一
【審査官】蔵田 真彦
(56)【参考文献】
【文献】米国特許出願公開第2019/0053780(US,A1)
【文献】特表2019-535448(JP,A)
【文献】特開2019-097795(JP,A)
(58)【調査した分野】(Int.Cl.,DB名)
A61B 8/00-8/15
(57)【特許請求の範囲】
【請求項1】
画像処理装置であって、
プロセッサとメモリとを含み、
前記メモリは、対象の第1の複数の画像の時系列を示す行列を保持し、
前記プロセッサは、
前記行列に対して、画像生成処理を実行して、前記対象の第1の複数の画像を生成し、
前記画像生成処理において、
前記行列を特異値分解して左行列、特異値行列、及び右行列を取得し、
前記右行列を
行方向にフーリエ変換することで周波数に変換し、
前記変換した周波数に基づいて、前記右行列
の各列に対応する周波数のパワー分布を算出し、
前記パワー分布に基づいて、
前記右行列の列ごとの周波数パワーを示すカーブを算出し、
所定の条件に基づいて、前記算出したカーブが示す周波数パワーが一定となる列番号を上側閾値に決定し、
前記特異値行列において前記上側閾値以上の行番号の行に位置する特異値を含む特異値群の値を0に変換し、
前記左行列、前記変換した特異値行列、及び前記右行列と、に基づいて、前記対象の第1の複数の画像を生成する、画像処理装置。
【請求項2】
請求項1に記載の画像処理装置であって、
前記プロセッサは、
前記画像生成処理において、
前記特異値行列の各特異値が示す情報量を算出し、
前記算出した情報量が所定値と最も近い特異値の前記特異値行列における行番号を下側閾値に決定し、
前記上側閾値
が前記下側閾値
に1を加えた値以下であると判定した場合には、
前記上側閾値と前記下側閾値
との少なくとも一方を再設定し、
前記特異値行列において前記上側閾値以上の行番号の行に含まれる特異値と、前記下側閾値以下の行番号の行に含まれる特異値と、を0に変換する、画像処理装置。
【請求項3】
請求項1又は2に記載の画像処理装置であって、
前記第1の複数の画像それぞれは、血管を含む組織が撮像された超音波画像であり、
前記左行列は、前記第1の複数の画像の空間に関する情報を示し、
前記右行列は、前記時系列における時間に関する情報を示す、画像処理装置。
【請求項4】
請求項2に記載の画像処理装置であって、
前記対象は、第1の対象と第2の対象とを含み、
前記メモリは、前記対象の第2の複数の画像の時系列を示す行列を保持し、
前記プロセッサは、
前記第2の複数の画像の時系列を示す行列に対して、前記画像生成処理を実行して、前記対象の第2の複数の画像を生成し、
前記生成した第1の複数の画像から前記対象の第1の画像を生成し、
前記生成した第2の複数の画像から前記対象の第2の画像を生成し、
所定の条件に基づいて、前記生成した第1の画像に含まれる前記第1の対象の高密度領域と、前記生成した第2の画像に含まれる前記第2の対象の高密度領域と、を特定し、
前記生成した第1の画像に含まれる前記第1の対象の高密度領域と、前記第2の画像に含まれる前記生成した第2の対象の高密度領域と、で前記生成した第1の画像と、前記生成した第2の画像と、を位置合わせして、前記対象の第3の画像を生成する、画像処理装置。
【請求項5】
請求項2に記載の画像処理装置であって、
前記対象は、第1の対象を含み、
前記プロセッサは、
前記左行列と、前記下側閾値以上の行番号の行に位置する特異値を0に変換した前記特異値行列と、及び前記右行列と、に基づいて、前記対象の第2の複数の画像を生成し、
前記生成した複数の画像に含まれる前記対象の移動解析をし、
前記移動解析の結果に基づいて、前記生成した第1の複数の画像それぞれに含まれる前記第1の対象を移動補正し、
前記移動補正した第1の対象を含む画像を生成する、画像処理装置。
【請求項6】
請求項2に記載の画像処理装置であって、
前記メモリは、前記第1の複数の画像を保持し、
前記プロセッサは、
前記カーブに基づいて、前記右行列の前記下側閾値以下の列番号の列それぞれが示す周波数の範囲における、平衡周波数及びピーク周波数の少なくとも一方を含む指標を算出し、
前記指標に基づいて、IIRフィルタを生成し、
前記第1の複数の画像に対して、流速解析を実行し、
前記流速解析の結果に前記IIRフィルタを適用して、前記第1の複数の画像を生成する、画像処理装置。
【請求項7】
請求項1に記載の画像処理装置であって、
出力装置を含み、
前記プロセッサは、前記パワー分布又は前記カーブの少なくとも一方を、前記出力装置に出力する、画像処理装置。
【請求項8】
画像処理装置による画像処理方法であって、
前記画像処理装置は、プロセッサとメモリとを含み、
前記メモリは、対象の第1の複数の画像の時系列を示す行列を保持し、
前記画像処理方法は、
前記プロセッサが、前記行列に対して、画像生成処理を実行して、前記対象の第1の複数の画像を生成し、
前記画像生成処理において、
前記プロセッサが、前記行列を特異値分解して左行列、特異値行列、及び右行列を取得し、
前記プロセッサが、前記右行列を
行方向にフーリエ変換することで周波数に変換し、
前記プロセッサが、前記変換した周波数に基づいて、前記右行列
の各列に対応する周波数のパワー分布を算出し、
前記プロセッサが、前記パワー分布に基づいて、
前記右行列の列ごとの周波数パワーを示すカーブを算出し、
前記プロセッサが、所定の条件に基づいて、前記算出したカーブが示す周波数パワーが一定となる列番号を上側閾値に決定し、
前記プロセッサが、前記特異値行列において前記上側閾値以上の行番号の行に位置する特異値を含む特異値群の値を0に変換し、
前記プロセッサが、前記左行列、前記変換した特異値行列、及び前記右行列と、に基づいて、前記対象の第1の複数の画像を生成する、画像処理方法。
【請求項9】
請求項8に記載の画像処理方法であって、
前記画像生成処理において、
前記プロセッサが、前記特異値行列の各特異値が示す情報量を算出し、
前記プロセッサが、前記算出した情報量が所定値と最も近い特異値の前記特異値行列における行番号を下側閾値に決定し、
前記上側閾値
が前記下側閾値
に1を加えた値以下であると判定した場合には、
前記上側閾値と前記下側閾値
との少なくとも一方を再設定し、
前記プロセッサが、前記特異値行列において前記上側閾値以上の行番号の行に含まれる特異値と、前記下側閾値以下の行番号の行に含まれる特異値と、を0に変換する、画像処理方法。
【請求項10】
請求項8又は9に記載の画像処理方法であって、
前記第1の複数の画像それぞれは、血管を含む組織が撮像された超音波画像であり、
前記左行列は、前記第1の複数の画像の空間に関する情報を示し、
前記右行列は、前記時系列における時間に関する情報を示す、画像処理方法。
【請求項11】
請求項9に記載の画像処理方法であって、
前記対象は、第1の対象と第2の対象とを含み、
前記メモリは、前記対象の第2の複数の画像の時系列を示す行列を保持し、
前記画像処理方法は、
前記プロセッサが、前記第2の複数の画像の時系列を示す行列に対して、前記画像生成処理を実行して、前記対象の第2の複数の画像を生成し、
前記プロセッサが、前記生成した第1の複数の画像から前記対象の第1の画像を生成し、
前記プロセッサが、前記生成した第2の複数の画像から前記対象の第2の画像を生成し、
前記プロセッサが、所定の条件に基づいて、前記生成した第1の画像に含まれる前記第1の対象の高密度領域と、前記生成した第2の画像に含まれる前記第2の対象の高密度領域と、を特定し、
前記プロセッサが、前記生成した第1の画像に含まれる前記第1の対象の高密度領域と、前記第2の画像に含まれる前記生成した第2の対象の高密度領域と、で前記生成した第1の画像と、前記生成した第2の画像と、を位置合わせして、前記対象の第3の画像を生成する、画像処理方法。
【請求項12】
請求項9に記載の画像処理方法であって、
前記対象は、第1の対象を含み、
前記画像処理方法は、
前記プロセッサが、前記左行列と、前記下側閾値以上の行番号の行に位置する特異値を0に変換した前記特異値行列と、及び前記右行列と、に基づいて、前記対象の第2の複数の画像を生成し、
前記プロセッサが、前記生成した複数の画像に含まれる前記対象の移動解析をし、
前記プロセッサが、前記移動解析の結果に基づいて、前記生成した第1の複数の画像それぞれに含まれる前記第1の対象を移動補正し、
前記プロセッサが、前記移動補正した第1の対象を含む画像を生成する、画像処理方法。
【請求項13】
請求項9に記載の画像処理方法であって、
前記メモリは、前記第1の複数の画像を保持し、
前記画像処理方法は、
前記プロセッサが、前記カーブに基づいて、前記右行列の前記下側閾値以下の列番号の列それぞれが示す周波数の範囲における、平衡周波数及びピーク周波数の少なくとも一方を含む指標を算出し、
前記プロセッサが、前記指標に基づいて、IIRフィルタを生成し、
前記プロセッサが、前記第1の複数の画像に対して、流速解析を実行し、
前記プロセッサが、前記流速解析の結果に前記IIRフィルタを適用して、前記第1の複数の画像を生成する、画像処理方法。
【請求項14】
請求項8に記載の画像処理方法であって、
前記画像処理装置は、出力装置を含み、
前記画像処理方法は、前記プロセッサが、前記パワー分布又は前記カーブの少なくとも一方を、前記出力装置に出力する、画像処理方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、画像処理装置及び画像処理方法に関する。
【背景技術】
【0002】
超音波撮像装置やMRI(Magnetic Resonance Imaging)装置、X線CT(Computed Tomography)装置に代表される医療用の画像表示装置は、目視できない生体内の情報を数値又は画像の形態で提示する装置として広く利用されている。特に、超音波撮像装置は、他の装置と比較して高い時間分解能を備えており、拍動下の心臓を滲みなく画像化できる性能を有する。
【0003】
非特許文献1に記載の技術は、微細な血管網を画像情報として表示するための技術であり、特異値分解(Singular Value Decomposition)に基づく血流映像化技術を開示している。非特許文献1に記載の技術は、複数枚分の画像のエコー信号から2次元行列(Casorati行列)を生成し、2次元行列に対して特異値分解を施して、空間に関する特異ベクトルと時間に関する特異ベクトルと、を取得する。また、非特許文献1に記載は、特異ベクトルを周波数方向に対するフィルタリングをする処理を開示している。
【0004】
また、特許文献1は、画素ごとの相関行列を平均化してから固有値分解を適用し、一定範囲の固有ベクトルを排除する方法を開示している。さらに、特許文献1に記載の技術は、処理範囲を小領域に分割し、各小領域にて相関行列の平均を取る処理を開示している。
【先行技術文献】
【特許文献】
【0005】
【非特許文献】
【0006】
【文献】C. Demene, T. Deffieux, M. Pernot, B. Osmanski, V. Biran, J. Gennisson, L. Sieu, A. Bergel, S. Franqui, J. Correas, I. Cohen, O. Baud, M. Tanter, IEEE Trans. Med. Img., vol. 34, no.11, p.2271, 2015年.
【発明の概要】
【発明が解決しようとする課題】
【0007】
非特許文献1には、雑音成分を抑制するための閾値を設定することについては記載されておらず、当該閾値の必要性についても記載されてない。特許文献1には、雑音成分について記載されていない。つまり、非特許文献1に記載の技術及び特許文献1に記載の技術は、雑音成分が除去された画像を生成することができないおそれがある。そこで、本発明の一態様は、雑音成分が高精度に除去された画像を生成することを目的とする。
【課題を解決するための手段】
【0008】
上記課題を解決するため、本発明の一態様は以下の構成を採用する。画像処理装置は、プロセッサとメモリとを含み、前記メモリは、対象の第1の複数の画像の時系列を示す行列を保持し、前記プロセッサは、前記行列に対して、画像生成処理を実行して、前記対象の第1の複数の画像を生成し、前記画像生成処理において、前記行列を特異値分解して左行列、特異値行列、及び右行列を取得し、前記右行列を各列が示す周波数に変換し、前記変換した周波数に基づいて、前記右行列に含まれる列群の列それぞれが示す周波数のパワー分布を算出し、前記パワー分布に基づいて、前記列群の列それぞれが示す周波数パワーを示すカーブを算出し、所定の条件に基づいて、前記算出したカーブが示す周波数パワーが一定となる列番号を上側閾値に決定し、前記特異値行列において前記上側閾値以上の行番号の行に位置する特異値を含む特異値群の値を0に変換し、前記左行列、前記変換した特異値行列、及び前記右行列と、に基づいて、前記対象の第1の複数の画像を生成する。
【発明の効果】
【0009】
本発明によれば、雑音成分が高精度に除去された画像を生成することができる。
【0010】
上記した以外の課題、構成及び効果は、以下の実施形態の説明により明らかにされる。
【図面の簡単な説明】
【0011】
【
図1】実施例1における画像処理装置のブロック図の一例を示すブロック図である。
【
図2】実施例1における超音波画像処理の一例を示すフローチャートである。
【
図3】実施例1における超音波画像処理に含まれる閾値設定処理の一例を示す説明図である。
【
図4】実施例2における治療効果判定処理の一例を示すフローチャートである。
【
図5】実施例2における治療効果判定処理の一例を示す説明図である。
【
図6】実施例2における位置合わせ処理の一例を示す説明図である。
【
図7】実施例3における移動補正を含む血管画像処理の一例を示すフローチャートである。
【
図8】実施例3における移動補正を含む血管画像処理の一例を示す説明図である。
【
図9】実施例4における周波数フィルタ最適化処理の一例を示すフローチャートである。
【
図10】実施例4における周波数フィルタ最適化処理の一例を示す説明図である。
【発明を実施するための形態】
【0012】
以下、本発明の実施形態を図面に基づいて詳細に説明する。本実施形態において、同一の構成には原則として同一の符号を付け、繰り返しの説明は省略する。なお、本実施形態は本発明を実現するための一例に過ぎず、本発明の技術的範囲を限定するものではないことに注意すべきである。
【0013】
本実施形態は、画像処理装置について説明する。微細血管網の映像情報は、腫瘍鑑別及び炎症検査への有効性が高いと期待される情報である。乳腺組織、腹部臓器、及び筋肉等、生体内の様々な組織が検査対象として想定される。検査対象が置かれている環境により体動の影響は異なり、血管径や流速も多様である。特に深部に位置する低速血流が検査対象である場合、エコー信号に含まれる血流成分は微弱であり、当該エコー信号に含まれる雑音成分との区別が難しくなる。
【0014】
従って、画像処理装置が特異値分解型フィルタにより血流情報を取得するためには、エコー信号から組織成分や体動成分などのクラッタ成分を排除する他、画像処理装置内に備える増幅器等に起因する雑音成分を排除する必要がある。そして、実用性の観点から、画像処理装置は、クラッタ成分および雑音成分を抑制するために抽出する特異ベクトルのランク範囲を指定するための閾値(後述するランク閾値)は、自動的に設定する。特に、検査者が探触子を手動で操作しながら検査対象を探索する、超音波検査に特有の扱いを想定すれば、撮像面の変化に応じてランク閾値が自動的に切り替わる機能の必要性は高い。
【0015】
本実施形態の画像処理装置は、エコー信号から、自動的に組織成分や体動成分などのクラッタ成分を排除するのみならず、雑音成分も排除することにより、血管情報を高精度に取得することができる。なお、本実施形態では、画像処理装置が、超音波画像に含まれる血管情報を取得する例を主に説明しているが、任意の画像(例えば、CT(Computed Tomography)画像やMRI(Magnetic Resonance Imaging)画像)の任意の情報(例えば組織情報やリンパ情報)を取得する場合にも本実施形態の画像処理装置を適用することができる。
【実施例1】
【0016】
図1は、画像処理装置のブロック図の一例を示すブロック図である。画像処理装置100は、例えば圧電素子によって実現される探触子200に接続され、探触子200によって得られた超音波画像に対する画像処理を行う。画像処理装置100は、例えば、CPU(Central Processing Unit)110、メモリ120、補助記憶装置130、通信装置140、入力装置150、及び出力装置160を有する計算機によって構成される。
【0017】
CPU110は、プロセッサを含み、メモリ120に格納されたプログラムを実行する。メモリ120は、不揮発性の記憶素子であるROM(Read Only Memory)及び揮発性の記憶素子であるRAM(Random Access Memory)を含む。ROMは、不変のプログラム(例えば、BIOS(Basic Input/Output System))などを格納する。RAMは、DRAM(Dynamic Random Access Memory)のような高速かつ揮発性の記憶素子であり、CPU110が実行するプログラム及びプログラムの実行時に使用されるデータを一時的に格納する。
【0018】
補助記憶装置130は、例えば、磁気記憶装置(HDD(Hard Disk Drive))、フラッシュメモリ(SSD(Solid State Drive))等の大容量かつ不揮発性の記憶装置であり、CPU110が実行するプログラム及びプログラムの実行時に使用されるデータを格納する。すなわち、プログラムは、補助記憶装置130から読み出されて、メモリ120にロードされて、CPU110によって実行される。
【0019】
なお、本実施形態において、画像処理装置100が使用する情報は、データ構造に依存せずどのようなデータ構造で表現されていてもよい。例えば、リスト、テーブル、データベース、又はキューから適切に選択したデータ構造体が、情報を格納することができる。
【0020】
入力装置150は、マウスやキーボード等の、オペレータからの入力を受けつけるため装置である。出力装置160は、ディスプレイ装置やプリンタ等の、プログラムの実行結果をオペレータが視認可能な形式で出力する装置である。なお、入力装置150や出力装置160が画像処理装置100に含まれずに、画像処理装置100が有する入力インターフェースに入力装置150が接続され、画像処理装置100が有する出力インターフェースに出力装置160が接続されていてもよい。
【0021】
通信装置140は、所定のプロトコルに従って、他の装置との通信を制御するネットワークインターフェース装置である。また、通信装置140は、例えば、USB(Univwsal Serial Bus)等のシリアルインターフェースを含む。なお、画像処理装置100は、通信装置140を含まなくてもよい。
【0022】
CPU110が実行するプログラムは、リムーバブルメディア(CD-ROM、フラッシュメモリなど)又はネットワークを介して画像処理装置100に提供され、非一時的記憶媒体である不揮発性の補助記憶装置130に格納される。このため、画像処理装置100は、リムーバブルメディアからデータを読み込むインターフェースを有するとよい。
【0023】
画像処理装置100は、物理的に一つの計算機上で、又は、論理的又は物理的に構成された複数の計算機上で構成される計算機システムであり、同一の計算機上で別個のスレッドで動作してもよく、複数の物理的計算機資源上に構築された仮想計算機上で動作してもよい。
【0024】
CPU110は、例えば、送受信制御部111及び信号処理部112を含む。送受信制御部111は、送信ビームフォーマ113及び受信ビームフォーマ114を含み、これらによる処理を制御する。送信ビームフォーマ113は、探触子200が検査対象に向けて送信する超音波信号を所定の時間遅延を素子毎に与えて生成する。受信ビームフォーマ114は、当該時間遅延の補正を含めて検査対象からの受信信号を生成し、当該受信信号からI/Qデータを生成する。
【0025】
信号処理部112は、行列変換部115、特異値解析部116、特異ベクトル解析部117、フィルタ解析部118、画像再構成部119、及び画質調整部1110を含む。行列変換部115は、送受信制御部111が取得した複数枚の超音波画像のI/Qデータを行列(例えば行方向が空間を、列方向が時間を示す行列)に変換する。特異値解析部116は、行列変換部115が生成した行列に対して特異値分解を施し、特異値分解によって得られた特異値を解析し、第1のランク閾値(下側閾値)を決定する。ランク閾値の詳細については後述する。
【0026】
特異ベクトル解析部117は、特異値解析部116による特異値分解によって得られた時間的な特異ベクトルを周波数解析し、第2のランク閾値(上側閾値)を決定する。フィルタ解析部118は、第1のランク閾値と第2のランク閾値とによって決定される範囲の特異ベクトルを抽出する。画像再構成部119は、抽出された特異ベクトルを、超音波画像の取得枚数分のI/Qデータに再変換する。画質調整部1110は、エッジ強調、背景除去、及びスムージング効果等の画像処理を行い、画像処理が行われた画像を出力装置160に表示する。
【0027】
例えば、CPU110は、メモリ120にロードされた送受信制御プログラムに従って動作することで、送受信制御部111として機能し、メモリ120にロードされた信号処理プログラムに従って動作することで、信号処理部112として機能する。CPU110に含まれる他の機能部(即ち、送受信制御部111及び信号処理部112それぞれに含まれる機能部)についても、プログラムと機能部の関係は同様である。
【0028】
なお、CPU110に含まれる機能部による機能の一部又は全部が、例えば、ASIC(Application Specific Integrated Circuit)やFPGA(Field-Programmable Gate Array)等のハードウェアによって実現されてもよい。例えば、送受信制御部111は、FPGAによって実現されることが望ましい。
【0029】
図2は、超音波画像処理の一例を示すフローチャートであり、
図3は、超音波画像処理に含まれる閾値設定処理(後述するS201~S207)の一例を示す説明図である。送受信制御部111は、探触子200から検査対象のNフレーム分の画像(必ずしも連続するフレームの画像である必要はないが複数の異なる時刻に撮影された画像)それぞれのI/Qデータを取得し、メモリ120及び/又は補助記憶装置130に格納する(S201)。各画像の空間サイズは(Z,X)である、つまり送受信制御部111は、検査対象のZ×X画素のNフレーム分の画像を取得する。
【0030】
行列変換部115は、空間サイズ(Z,X)のNフレーム分の画像を、Z行X列に変換する(S202)。具体的には、例えば、行列変換部115は、Nフレーム分のZ×X行列それぞれについて、全ての要素を縦に並べた列ベクトルを生成し、当該生成した列ベクトルをフレームの順序で列方向に並べた(例えば時間が古いフレームほど左側に位置するように並べた)ZX行N列のCasorati行列を生成する。つまり、このZX×N行列は、検査対象のNフレームの画像の時系列を示し、このZX×N行列を構成する各列ベクトルは、画像情報の1フレーム分のI/Qデータを示す。
【0031】
特異値解析部116は、このZX×N行列を特異値分解することで、VΛWを得る(S203)。ここで、Vは空間に関する情報を示すZX行ZX列の左特異ベクトル行列であり、Wは時間に関する情報を示すN行N列の右特異ベクトル行列である。
【0032】
ΛはZX行N列の特異値行列である。なお、特異値行列は、各特異値を対角成分とする対角行列が行方向(下側)及び/又は列方向(右側)に拡張されたZX行N列の行列であり、当該拡張された部分の成分は全て0である。特異値行列の対角成分をλ1,・・・,λq(q=min{ZX,N})としたとき、λ1≧λ2≧・・・≧λq-1≧λq≧0であるとする。また、特異値行列のランクは(0でない)特異値の個数と等しい。以下、説明の便宜のために、λk(1≦k≦q)をランクkの特異値とも呼び、行列の積であるVΛWを算出する際にλk(1≦k≦q)を含む成分に掛けられる行をランクkの行とも呼び、VΛWを算出する際にλk(1≦k≦q)を含む成分に掛ける列をランクkの列とも呼ぶ。
【0033】
特異ベクトル解析部117は、特異値行列Λの各対角成分を、特異値行列Λの対角成分の総和で除算した商(σi=λi/Σλi(Σの範囲はi=1からi=q))を算出する(S204)ことで、各ランクの成分が全てのI/Qデータに占める情報の割合を算出することができる。
【0034】
一般的に、生体からのエコー信号のうち、クラッタ成分と血流成分とには、約40dB(100対1)程度の強度差がある。従って、したがって、特異ベクトル解析部117は、情報量(前述した商)が0.01(予め定められた値)に最も近い値を示した特異値に対応するランクを、クラッタ成分と血流成分とを判別する第1のランク閾値kcr(単に第1閾値とも呼ぶ)として設定する(S205)。なお、送受信制御部111が取得する画像の種類や分析の目的等によって、前述した0.01という値を、画像処理装置100のユーザが変更可能であることが望ましい。
【0035】
なお、特異ベクトル解析部117は、
図3に示す、情報量とランクとの関係を示すグラフ、並びに当該グラフ上に示された血流閾値(上述した予め定められた値)及び第1のランク閾値を、出力装置160に出力してもよい。
【0036】
なお、特異ベクトル解析部117は、情報量(前述した商)が0.01(予め定められた値)に最も近い値を示した特異値に対応する複数のランクがあると判定した場合には、例えば、当該複数のランクのうち特異値行列において最も左上に位置する特異値に対応するランクを第1のランク閾値として設定してもよいし、最も右下に位置する特異値に対応するランクを第1のランク閾値として設定してもよいし、これらの中央に位置する特異値に対応するランクを第1のランク閾値として設定してもよい。
【0037】
続いて、特異ベクトル解析部117は、右特異ベクトル行列Wを行方向にフーリエ変換することで右特異ベクトル行列Wを周波数変換する(S206)。特異ベクトル解析部117は、フーリエ変換後の右特異ベクトル行列Wに対して周波数解析を行って、第2のランク閾値knr(単に第2閾値とも呼ぶ)を算出する(S207)。ステップS207の具体例について説明する。
【0038】
高周波域かつ特異値が低い高次側(即ち特異値行列においてより右下に位置する対角成分に対応するランク)のランク範囲(つまり当該ランク範囲の特異値に掛ける右特異ベクトル行列Wの列)が除去すべき雑音成分の支配域である。従って、例えば、特異ベクトル解析部117は、高周波域(例えば1kH(予め定められた値)以上の範囲、全範囲の周波数について算出してもよい)のランクごとの(即ち右特異ベクトル行列Wの各列に対応する)周波数のパワー分布を算出し、算出したパワー分布から、ランクごとの周波数パワーの平均値を計算することで、エネルギーカーブを算出する。
【0039】
エネルギーカーブが示す周波数パワーは、クラッタや血流の成分を含むランク範囲より高いランクにおいてほぼ一定、つまりクラッタや血流の成分を含むランク範囲より高いランクにおいては雑音成分が支配的である。つまりエネルギーカーブが示す周波数パワーが一定となる起点付近が、血流成分と雑音成分とを判別する第2のランク閾値として適当である。
【0040】
具体的には、特異ベクトル解析部117は、例えば、エネルギーカーブの微分係数が0となる始点のランクを第2のランク閾値に決定する。また、特異ベクトル解析部117は、例えば、エネルギーカーブの始点と終点を結ぶ直線が示す周波数パワーと、エネルギーカーブが示す周波数パワーと、の差分が最大となるランクを第2のランク閾値に決定してもよい。また、特異ベクトル解析部117は、例えば、エネルギーカーブの変曲点を通過する直線グラフ(例えばエネルギーカーブに対するフィッテイング処理を用いて算出される)と、最終ランク(即ち0でない特異値のうち最も小さい特異値に対応するランク)を通過する水平の直線グラフと、の交点のような、エネルギーカーブに対する解析結果から得られる起点を推定し、当該起点が示すランクによって第2のランク閾値を設定する。
【0041】
なお、特異ベクトル解析部117は、第2のランク閾値として複数の値の候補があると判定した場合には、例えば、当該複数のランクのうち特異値行列において最も左上に位置する特異値に対応するランクを第2のランク閾値として設定してもよいし、最も右下に位置する特異値に対応するランクを第2のランク閾値として設定してもよいし、これらの中央に位置する特異値に対応するランクを第2のランク閾値として設定してもよい。
【0042】
また、特異ベクトル解析部117は、第1のランク閾値+1≧第2のランク閾値、であると判定した場合、例えば、第1のランク閾値及び/又は第2のランク閾値を再設定するとよい。具体的には、例えば、特異ベクトル解析部117は、第1のランク閾値に2を加えた値以上のランクから、ステップS207と同様の方法で第2のランク閾値を再設定する。また、特異ベクトル解析部117は、単に第1のランク閾値又は第2のランク閾値の一方を固定して、第2のランク閾値-第1のランク閾値=2以上の所定値(例えば2)となるように、他方を再設定するようにしてもよい。
【0043】
なお、特異ベクトル解析部117は、
図3に示すパワー分布、高周波域の平均パワー推移(第2閾値を含む)、及び差分グラフ(第2閾値を含む)を出力装置160に出力してもよい。
【0044】
続いて、フィルタ解析部118は、特異値行列Λの、第1のランク閾値より小さい、又は第2のランク閾値よりも大きいランクに対応する特異値を全て0に変換した行列である血流抽出フィルタΛ(k)を生成する(S208)。フィルタ解析部118は、ZX×N行列であるS=V(Λ(k))WH(但し、WHはWの共役転置行列)を算出する(S209)。画像再構成部119は、ZX×N行列のSを、N個のZ×X行列へと分解する(ステップS202で実行された変換の逆変換をする)ことで、Z×X画素のNフレーム分の2次元画像を再構成する(S210)。
【0045】
これにより、フィルタ解析部118は、第1のランク閾値と第2のランク閾値との間で規定されるランク範囲によって血流成分を抽出し、元のI/Qデータと同じサイズのデータを生成することができる。
【0046】
画質調整部1110は、ステップS211~ステップS214において、血管画像としての視認性を向上する画像フィルタを適用する。なお、ステップS211~ステップS214に示す処理はあくまで画像フィルタ適用の例であり、ステップS211~ステップS214の少なくとも一部の処理が行われなくてもよいし、さらに追加の画質調整処理が行われてもよい。
【0047】
画質調整部1110は、再構成されたZ×X画素のNフレーム分の画像を時間軸方向に加算平均する(S211)ことで、画像の背景雑音を抑制することができる。さらに、画質調整部1110は、再構成されたZ×X画素のNフレーム分の画像の背景除去を行う(S212)ことで、血管領域の信号をより強調することができる。画像全体の平均値を各画素の輝度から一律に差分する方法は、背景除去の一例である。
【0048】
画質調整部1110は、再構成されたZ×X画素のNフレーム分の画像に対し、移動平均とメディアンフィルタを適用する(S213)。画質調整部1110は、さらに構造物の輪郭を強調する画像処理を適用する。例えば、画質調整部1110は、特定方向の輪郭を強調する機能を備えるモルフォロジフィルタを、入力画像に対して複数の角度について適用し、モルフォロジフィルタ適用後の画像に対して加算平均を行う等の多重解像度処理を行う(S214)。これにより、画質調整部1110は、画像内の血管部の繋がりを改善し、輪郭が明瞭化された画像を得ることができる。
【0049】
続いて、画質調整部1110は、画質調整後の画像を出力装置160に出力する(S215)。以上の超音波画像処理により、術者は探触子200を検査対象に向けて設置するだけで、明瞭な血管網の情報を出力装置160に表示される映像情報として取得することができる。
【0050】
なお、本実施例の画像処理装置100は、超音波画像から血流成分を得るために第1のランク閾値と第2のランク閾値とを算出したが、組織成分を主とするクラッタ成分のみを得たいのであれば第1のランク閾値のみを算出して第1のランク閾値以上のランクに対応する特異値を0に変換してもよいし、雑音のみを除去したいのであれば第2のランク閾値のみを算出して第2のランク閾値以上のランクに対応する特異値を0に変換してもよい。
【実施例2】
【0051】
本実施例では、治療効果判定処理について説明する。実施例1と同様の構成や処理については説明を一部省略する。
図4は、治療効果判定処理の一例を示すフローチャートであり、
図5は、治療効果判定処理の一例を示す説明図である。
図6は、後述するステップS407の位置合わせ処理の一例を示す説明図である。
【0052】
例えば癌の化学療法に用いられる薬剤は、抗血管新生薬などの血管を縮退させる作用を持つものが少なくない。癌に対する放射線治療において効果がある場合も、血管網の縮退が認められる。
図4乃至
図6に示す治療効果判定処理は、治療前後の血管画像に対する画像処理を行うことで、治療効果判定の定量化する。
【0053】
送受信制御部111は、ステップS201と同様の処理により、治療前の検査対象(血管を含む対象)の複数の画像と、治療後の当該検査対象の複数の画像と、を収集する(S401、S402)。行列変換部115は、ステップS202と同様の処理によって、治療前の検査対象の複数の画像に対して行列変換を施し、治療後の検査対象の複数の画像に対しても行列変換を施し、特異値解析部116は、行列変換部115が生成した行列それぞれを特異値分解する(S403)。
【0054】
特異ベクトル解析部117は、特異値分解された行列それぞれの特異値行列及び右特異ベクトル行列に対して、ステップS204~ステップS207と同様の処理によって、第1のランク閾値及び第2のランク閾値を決定する(S404)。フィルタ解析部118、画像再構成部119、及び画質調整部1110は、ステップS208~ステップS214の処理によって、治療前の検査対象及び治療後の検査対象の血管が明瞭化された画像を得る(S405)。これにより、
図5に示すように、治療前と治療後とそれぞれについて、血管503(第1の対象の一例)を含む腫瘍領域502(第2の対象の一例)の画像データ501が得られる。
【0055】
画質調整部1110は、治療前と治療後それぞれの画像を小領域に分割し(例えば5mm四方)、各小領域について輝値を加算し血管値として算出し、各小領域の面積に対する血管値を当該小領域の血管密度として算出することで、血管密度画像を生成する(S406)。画質調整部1110は、このようにして血管密度を算出することにより、僅かな治療効果が反映される血管密度画像を算出することができる。
【0056】
なお、画質調整部1110は、血管密度を正確に試算するために、一定値以上の輝度の画素を血管領域と定義し、血管領域に該当する画素の総数を血管値に決定することが望ましい。画質調整部1110は、このように血管密度を算出することにより、流量や流速の影響を含みにくく、血管の有無だけを判定しやすい血管密度画像を生成することができる。
【0057】
続いて、画質調整部1110は、治療前の血管密度画像と治療後の血管密度画像と、を血管密度が高い領域同士で位置合わせする(S407)。例えば、画質調整部1110は、治療前後の血管密度分布画像を重畳させる。具体的には、例えば、画質調整部1110は、最も簡便な方法として、血管密度が最大となる画素を位置整合の基準として治療前後の血管密度分布画像を重畳させる。
【0058】
また、画質調整部1110は、
図6に示すように、治療前後の血管密度分布画像それぞれに対して重心を定義し、当該重心を基準に位置整合を行なってもよい。一般的に重心は質点系における基準点からの距離と質量から算出されるが、ここでは、画質調整部1110は、質量を血管密度に置き換えて重心として定義する。
【0059】
具体的には、画質調整部1110は、治療前の血管密度分布画像の小領域Piそれぞれについて、ri=OPi(OPiはPiの位置ベクトル)としたとき、治療前の血管密度分布画像の重心ベクトルrG=Σmiri/Σmi(miはPiの血管密度、Σの範囲はいずれもi=1からi=小領域の個数)を算出して、重心PPreを決定する。画質調整部1110は、同様の方法で、治療後の血管密度分布画像の重心PPostを決定する。画質調整部1110は、治療前の血管密度画像と治療後の血管密度画像との重心の位置が一致するように、これら2つの血管密度画像の位置合わせをする。画質調整部1110は、このように位置合わせを行うことで、血管密度画像において局所的に血管密度が高い領域があっても、その影響を軽減した正確な位置合わせを行うことができる。
【0060】
続いて、画質調整部1110は、位置合わせをした治療前と治療後の血管密度画像を用いて、効果判定指標を算出し(S408)、位置合わせをした血管密度画像と効果判定指標とを出力装置160に出力する(S409)。
【0061】
具体的には、例えば、画質調整部1110は、血管密度が一定値(予め定められた値又は画像処理装置100のユーザが設定可能な値)以上の範囲を定義して、ユーザの関心領域を可視化することができる。例えば、画質調整部1110は、治療前後の関心領域の面積を算出し、関心領域の面積変化率を効果判定指標として算出する。
【0062】
また、画質調整部1110は、関心領域の血管密度変化率を効果判定指標として算出してもよい。また、画質調整部1110は、腫瘍領域を球体と仮定して、関心領域の体積変化率を効果判定指標として算出してもよい。なお、画質調整部1110は、関心領域ではなく腫瘍領域502全体の面積変化率、血管密度変化率、及び体積変化率を効果判定指標として算出してもよい。
【0063】
本実施例の画像処理装置100により、各種治療の効果判定の定量化が実現し、ひいては、画像処理装置100のユーザは、客観的情報に基づく治療最適化(治療方法や薬剤選定の最適化)を迅速に行うことができる。
【実施例3】
【0064】
本実施例では、移動補正を含む血管画像処理について説明する。実施例1と同様の構成や処理については説明を一部省略する。
図7は、移動補正を含む血管画像処理の一例を示すフローチャートであり、
図8は、移動補正を含む血管画像処理の一例を示す説明図である。
【0065】
送受信制御部111は、ステップS201と同様の処理により、検査対象(血管を含む対象)の複数の画像を収集する(S701)。行列変換部115は、ステップS202と同様の処理によって、検査対象の複数の画像に対して行列変換を施し、特異値解析部116は、行列変換部115が生成した行列を特異値分解する(S702)。特異ベクトル解析部117は、特異値分解された行列の特異値行列及び右特異ベクトル行列に対して、ステップS204~ステップS207と同様の処理によって、第1のランク閾値及び第2のランク閾値を決定する(S703)。
【0066】
フィルタ解析部118は、特異値行列Λの、第1のランク閾値以上のランクに対応する特異値を全て0に変換した行列である組織抽出フィルタΛ(j)を生成し、ZX×N行列であるS’=V(Λ(j))WH(但し、WHはWの共役転置行列)を算出する(S704)。組織抽出フィルタにより、血流成分を含まないクラッタ成分を抽出することができ、組織の移動を推定することができる。画像再構成部119は、ZX×N行列のS’を、N個のZ×X行列へと分解する(ステップS202で実行された変換の逆変換をする)ことで、Z×X画素のNフレーム分の2次元組織画像を再構成する(S705)。
【0067】
画質調整部1110は、再構成したNフレーム分の2次元組織画像に対して、フレーム間に生じるZ方向及びX方向の移動量を計測する等の移動解析を行う(S706)。具体的には、例えば、画質調整部1110は、一般的に知られる有効な手法であるトラッキングアルゴリズム(特に例えば、正規化相互相関演算)によって移動量を計測する。
【0068】
フィルタ解析部118は、ステップS208と同様の処理によって、~ステップS214の処理と同様の処理によって血管画像を明瞭化して抽出するための血流抽出フィルタを生成し、ステップS209と同様の処理によってZX×N行列であるSを算出する(S707)。画像再構成部119は、ステップ210と同様の処理によって、Z×X画素のNフレーム分の2次元血管画像を再構成する(S708)。
【0069】
画質調整部1110は、各フレームの2次元血管画像に対して、対応するフレームのステップS706において計測した2次元組織画像の移動量に応じた移動補正(例えば移動量をそのまま差し引く)を行う(S709)。画質調整部1110は、移動補正後のNフレームの2次元血管画像に対して加算平均を行う等(加算平均に加えてステップS212~ステップS214の処理を行ってもよい)の画像解析処理行う(S710)。
【0070】
なお、画質調整部1110は、ステップS709における加算平均処理に代えて又は加えて、移動補正後の2次元血管画像の画素ごとの輝度値の時間変化を算出し、最大の輝度値を選択して、最大輝度値保持画像を取得してもよいし、移動補正後の2次元血管画像から血管内の血流の流量の変化値を示す画像を生成してもよいし、移動補正後の2次元血管画像から血流の流速を示す画像を取得してもよい。
【0071】
画質調整部1110は、ステップS710で生成した画像を出力装置160に出力する(S711)。なお、画質調整部1110は、ステップS710で得られた解析結果をさらに出力装置160に出力してもよい。本実施例の画像処理装置100は、組織画像から組織の移動量を推定し、推定した移動量を血管画像に適用することで、フレーム間で大きな体動があった場合に平均加算等を行っても血管画像に滲みが発生しづらくなる上に、他の血流や血管の変化値画像についても正確な画像を生成することができる。なお、画質調整部1110は、血管画像から血管の移動量を推定し、推定した移動量を組織画像に適用してもよい。
【実施例4】
【0072】
本実施例の画像処理装置100は、特異値分解によって得られるクラッタ成分の特性を分析し、分析結果を周波数フィルタの設計に適応的に反映させる周波数フィルタ最適化処理を行う。実施例1と同様の構成や処理については説明を一部省略する。
図9は、周波数フィルタ最適化処理の一例を示すフローチャートであり、
図10は、周波数フィルタ最適化処理の一例を示す説明図である。
【0073】
送受信制御部111は、ステップS201と同様の処理により、検査対象(血管を含む対象)の複数の画像を収集する(S901)。行列変換部115は、ステップS202と同様の処理によって、検査対象の複数の画像に対して行列変換を施し、特異値解析部116は、行列変換部115が生成した行列を特異値分解する(S902)。
【0074】
特異ベクトル解析部117は、ステップS204及びステップS205と同様の処理によって、第1のランク閾値を決定する(S903)。特異ベクトル解析部117は、ステップS206及びステップS207で説明した処理と同様の処理によって、全範囲の周波数のパワー分布を算出し、算出したパワー分布からランクごとの周波数パワーの平均値を計算してエネルギーカーブを算出する(S904)。
【0075】
特異ベクトル解析部117は、
図10に示す、情報量とランクとの関係を示すグラフ、当該グラフ上に示された血流閾値(上述した予め定められた値)及び第1のランク閾値、平均勾配、並びにパワー分布を出力装置160に出力してもよい。
【0076】
特異ベクトル解析部117は、第1のランク閾値以下のランクに対応する周波数の範囲である組織周波数(組織成分を主とするクラッタ成分の周波数)解析を行う(S905)。具体的には、例えば、特異ベクトル解析部117は、組織周波数における周波数ごとの平均パワーを算出する。例えば、第1のランク閾値以下のランクがランク1からランク5、組織周波数範囲が±3000Hzである場合、特異ベクトル解析部117は、5つのデータ平均によって得られる0~3000Hzの平均周波数パワーを算出結果として得る。特異ベクトル解析部117は、この平均周波数パワーを微分することにより、平衡周波数b(平均周波数パワーが一定、つまり微分曲線が0に近接する(又は0になる)周波数)や、ピーク周波数a(平均周波数の最大勾配、つまり微分曲線の極値となる周波数)、を算出することができる。
【0077】
また、例えば、特異ベクトル解析部117は、最小のランクであるランク1の周波数である0Hzを起点に、各ランクの最大パワーを選択して得られる周波数ピークグラフを算出し、当該ピークグラフの勾配(角度c)を算出する。また、例えば、特異ベクトル解析部117は、組織周波数の範囲において、周波数毎の平均値を算出する。
【0078】
なお、特異ベクトル解析部117は、
図10に示すクラッタ範囲の平均周波数パワー(ピーク周波数a及び平衡周波数bを含む)及び組織領域の周波数ピークグラフ(角度cを含む)を出力装置160に出力してもよい。
【0079】
続いて、特異ベクトル解析部117は、ステップS906で算出した情報に基づく解析を行う(S906)。具体的には、例えば、特異ベクトル解析部117は、クラッタのパワー、内在する周波数レンジ、及び速度を推定する。ピーク周波数aや平衡周波数bは、クラッタ範囲に含まれる周波数範囲と強さを示す指標であり、角度cは、この勾配が大きいほど、クラッタの速度は高くなる傾向があるため、c=Δf/ΔR(ランク範囲)が速度関連因子であるため、特異ベクトル解析部117は、これらの指標や関連因子を用いてクラッタ解析を行うことができる。なお、特異ベクトル解析部117は、クラッタ解析結果を出力装置160に出力してもよい。
【0080】
フィルタ解析部118は、ステップS906のクラッタ解析によって得られたクラッタ成分の推定結果から、IIR(Infinite Impluse Response)型フィルタの種類やカットオフ周波数等のパラメータを決定することで、適応型IIRフィルタを設計する(S907)。
【0081】
フィルタ解析部118は、元のI/Qデータに自己相関演算を適用して流速解析の結果を取得し(S908)、画像再構成部119は、流速解析の結果に対して適応型IIRフィルタを適用することにより、クラッタ成分を抑制した血流画像(及び流速画像)を生成する(S909)。なお、画質調整部1110が、例えば、ステップS211~ステップS214等の画質調整処理を、ステップS908で得られた血流画像(及び流速画像)に対して行ってもよい。画質調整部1110は、クラッタ成分を抑制した血流画像(及び流速画像)を出力装置160に出力する(S911)。
【0082】
本実施例の画像処理装置100は、フレームの取り込みと同時にIIRフィルタも再設計できるため、検査者が探触子200を操作していても、常に適切な血流画像を表示することが可能であり、ひいては高精度な超音波診断が実現する。
【0083】
なお、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。例えば、上記した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることも可能であり、また、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。
【0084】
また、上記の各構成、機能、処理部、処理手段等は、それらの一部又は全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリや、ハードディスク、SSD(Solid State Drive)等の記録装置、または、ICカード、SDカード、DVD等の記録媒体に置くことができる。
【0085】
また、制御線や情報線は説明上必要と考えられるものを示しており、製品上必ずしも全ての制御線や情報線を示しているとは限らない。実際には殆ど全ての構成が相互に接続されていると考えてもよい。
【符号の説明】
【0086】
100 画像処理装置、110 CPU、111 送受信制御部、112 信号処理部、113 送信ビームフォーマ、114 受信ビームフォーマ、115 行列変換部、116 特異値解析部、117 特異ベクトル解析部、118 フィルタ解析部、119 画像再構成部、120 メモリ、130 補助記憶装置、140 通信装置、150 入力装置、160 出力装置、200 探触子、1110 画質調整部