JP5932184B1 - 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム - Google Patents

超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム Download PDF

Info

Publication number
JP5932184B1
JP5932184B1 JP2016505643A JP2016505643A JP5932184B1 JP 5932184 B1 JP5932184 B1 JP 5932184B1 JP 2016505643 A JP2016505643 A JP 2016505643A JP 2016505643 A JP2016505643 A JP 2016505643A JP 5932184 B1 JP5932184 B1 JP 5932184B1
Authority
JP
Japan
Prior art keywords
frequency
unit
ultrasonic
candidate
frequency band
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2016505643A
Other languages
English (en)
Other versions
JPWO2016103848A1 (ja
Inventor
知宏 中辻
知宏 中辻
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Olympus Corp
Original Assignee
Olympus Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Olympus Corp filed Critical Olympus Corp
Priority claimed from PCT/JP2015/078547 external-priority patent/WO2016103848A1/ja
Application granted granted Critical
Publication of JP5932184B1 publication Critical patent/JP5932184B1/ja
Publication of JPWO2016103848A1 publication Critical patent/JPWO2016103848A1/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

本発明にかかる超音波観測装置は、超音波信号の周波数を解析することによって超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部と、周波数解析部が算出した各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連したパラメータを用いて設定する周波数帯域設定部と、周波数帯域設定部が設定した周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出部と、特徴量算出部が算出した特徴量に関する情報を表示する特徴量画像データを生成する特徴量画像データ生成部と、を備える。

Description

本発明は、超音波を用いて観測対象の組織を観測する超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラムに関する。
従来、超音波を用いて観測対象の組織を観察する超音波観測装置において、組織性状に応じた特性を有する超音波信号の周波数スペクトルの特徴量を算出し、この特徴量に基づいて、組織性状を鑑別するための特徴量画像を生成する技術が知られている(例えば、特許文献1を参照)。この技術では、受信した超音波信号の周波数を解析することによって周波数スペクトルを取得した後、所定の周波数帯域における周波数スペクトルの近似式を算出し、この近似式から特徴量を抽出している。
国際公開第2012/011414号
上述した特許文献1に記載の技術では、周波数スペクトルの特徴量を算出する際に設定する周波数帯域が、超音波の受信深度や観測対象の組織によらず一定である。このため、超音波の受信深度や組織に応じて適切な周波数帯域を設定しているとは限らず、観測対象の組織性状を鑑別する際の精度に影響を及ぼすおそれがあった。
本発明は、上記に鑑みてなされたものであって、観測対象の組織性状を精度よく鑑別することを可能にする超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラムを提供することを目的とする。
上述した課題を解決し、目的を達成するために、本発明に係る超音波観測装置は、観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置であって、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部と、前記周波数解析部が算出した各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定部と、前記周波数帯域設定部が設定した周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出部と、前記特徴量算出部が算出した特徴量に関する情報を表示する特徴量画像データを生成する特徴量画像データ生成部と、を備えたことを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記周波数帯域設定部は、前記周波数スペクトルの包絡線を検出する包絡線検出部と、前記複数の候補帯域において前記包絡線を近似する近似直線に関連するパラメータをそれぞれ抽出するパラメータ抽出部と、前記パラメータ抽出部が前記候補帯域ごとに抽出した複数のパラメータのうち所定の関係にある2つの前記候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出する変化量算出部と、前記変化量算出部が算出した前記変化量に基づいて前記周波数帯域を決定する決定部と、を有することを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記パラメータ抽出部は、前記候補帯域の周波数の最小値を固定したまま最大値を順次大きくして前記複数の候補帯域を設定し、1つの前記候補帯域を設定するごとに前記包絡線を直線で近似してパラメータを抽出した後、前記変化量算出部へ前記パラメータを出力し、前記変化量算出部は、前記パラメータ抽出部が少なくとも2回のパラメータ抽出処理を実行済みである場合には、最新のパラメータとその1つ前のパラメータ抽出処理で抽出されたパラメータとの差を前記変化量として算出し、前記決定部は、前記変化量算出部が少なくとも2回の変化量算出処理を実行済みである場合には、最新の前記変化量とその1つ前の変化量算出処理で算出された前記変化量とを比較し、該比較の結果、前記変化量の符号が0または一方の符号から他方の符号に変化した場合、前記1つ前のパラメータ抽出処理で前記パラメータを抽出した際の前記候補帯域を前記周波数帯域に決定し、前記パラメータ抽出部、前記変化量算出部および前記決定部は、前記決定部が前記周波数帯域を決定するまで処理を繰り返し実行することを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記パラメータ抽出部は、前記包絡線を回帰分析することによって前記直線を算出し、該算出した前記直線の傾きを前記パラメータとすることを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記パラメータ抽出部は、前記候補帯域の周波数の最小値を固定したまま最大値を順次大きくして前記複数の候補帯域を設定し、各候補帯域において前記包絡線を直線で近似してパラメータを抽出し、前記変化量算出部は、帯域幅の値が隣接する2つの前記候補帯域においてそれぞれ抽出された2つの前記パラメータの差を前記変化量として算出し、前記決定部は、前記変化量の絶対値を2つの前記候補帯域のうち帯域幅が小さい方の候補帯域の最大周波数の関数として近似曲線を算出し、該近似曲線において前記絶対値が最小となるときの前記最大周波数をもとに前記周波数帯域を決定することを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記パラメータ抽出部は、前記候補帯域の周波数の最小値を固定したまま最大値を順次大きくして前記複数の候補帯域を設定し、各候補帯域において前記包絡線を直線で近似してパラメータを抽出し、前記変化量算出部は、帯域幅の値が隣接する2つの前記候補帯域においてそれぞれ抽出された2つの前記パラメータの差を前記変化量として算出し、前記決定部は、前記変化量を2つの前記候補帯域のうち帯域幅が小さい方の候補帯域の最大周波数の関数として近似曲線を算出し、該近似曲線において値がゼロとなるときの前記最大周波数をもとに前記周波数帯域を決定することを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記パラメータ抽出部は、前記候補帯域における前記包絡線上のサンプリング点を複数の投影直線にそれぞれ投影したときに各投影直線における度数分布を算出し、該度数分布の統計的なばらつきが最も小さい投影直線に基づいて前記直線の傾きを算出し、該算出した前記直線の傾きを前記パラメータとすることを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記パラメータ抽出部は、前記候補帯域における前記包絡線上のサンプリング点を複数の投影直線にそれぞれ投影したときに各投影直線における度数分布を算出し、該度数分布の最頻値および統計的なばらつきを用いて定義される指標であって前記最頻値が大きいほど大きい値をとるとともに前記統計的なばらつきが小さいほど大きい値をとる指標に基づいて前記直線の傾きを算出する際に用いる投影直線を選択し、該選択した投影直線に基づいて前記直線の傾きを算出し、該算出した前記直線の傾きを前記パラメータとすることを特徴とする。
本発明に係る超音波観測装置は、上記発明において、前記特徴量算出部は、前記周波数帯域における前記包絡線を回帰分析によって近似した後、超音波が伝搬する際に発生する減衰の寄与を削減する減衰補正を行うことによって前記周波数スペクトルの特徴量を算出することを特徴とする。
本発明に係る超音波観測装置の作動方法は、観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置の作動方法であって、周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、周波数帯域設定部が、各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定ステップと、特徴量算出部が、前記周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出ステップと、特徴量画像データ生成部が、前記特徴量に基づいて特徴量画像データを生成する特徴量画像データ生成ステップと、を有することを特徴とする。
本発明に係る超音波観測装置の作動プログラムは、観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置に、周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、周波数帯域設定部が、各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定ステップと、特徴量算出部が、前記周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出ステップと、特徴量画像データ生成部が、前記特徴量に基づいて特徴量画像データを生成する特徴量画像データ生成ステップと、を実行させることを特徴とする。
本発明によれば、周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連したパラメータを用いて設定し、設定した周波数帯域をもとに各周波数スペクトルの特徴量を算出するため、周波数スペクトルを近似する際、解析に最適な周波数帯域を適応的に決定することができる。したがって、観測対象の組織性状を精度よく鑑別することが可能となる。
図1は、本発明の実施の形態1に係る超音波観測装置を備えた超音波診断システムの機能構成を示すブロック図である。 図2は、本発明の実施の形態1に係る超音波観測装置の信号増幅部が行う増幅処理における受信深度と増幅率との関係を示す図である。 図3は、本発明の実施の形態1に係る超音波観測装置の増幅補正部が行う増幅補正処理における受信深度と増幅率との関係を示す図である。 図4は、超音波信号の1つの音線におけるデータ配列を模式的に示す図である。 図5は、本発明の実施の形態1に係る超音波観測装置の包絡線検出部が検出した周波数スペクトルの包絡線の例を示す図である。 図6は、本発明の実施の形態1に係る超音波観測装置のパラメータ抽出部が実行するパラメータ抽出処理の概要を模式的に示す図である。 図7は、本発明の実施の形態1に係る超音波観測装置の決定部が実行する周波数帯域決定処理の概要を示す図である。 図8は、本発明の実施の形態1に係る超音波観測装置の近似部が実行する近似処理の概要を示す図である。 図9は、本発明の実施の形態1に係る超音波観測装置が実行する処理の概要を示すフローチャートである。 図10は、本発明の実施の形態1に係る超音波観測装置の周波数解析部が実行する処理の概要を示すフローチャートである。 図11は、本発明の実施の形態1に係る超音波観測装置の周波数帯域設定部が実行する周波数帯域設定処理の概要を示すフローチャートである。 図12は、本発明の実施の形態1の効果を説明する図である。 図13は、従来技術の概要を説明する図である。 図14は、本発明の実施の形態1の変形例1に係る超音波観測装置の周波数帯域設定部が実行する周波数帯域設定処理の概要を示す図である。 図15は、本発明の実施の形態1の変形例2に係る超音波観測装置の周波数帯域設定部が実行する周波数帯域設定処理の概要を示す図である。 図16は、本発明の実施の形態2に係る超音波観測装置のパラメータ抽出部が実行するパラメータ抽出処理の概要を示すフローチャートである。 図17Aは、本発明の実施の形態2に係る超音波観測装置のパラメータ抽出部が実行するパラメータ抽出処理の概要(第1例)を模式的に示す図である。 図17Bは、本発明の実施の形態2に係る超音波観測装置のパラメータ抽出部が実行するパラメータ抽出処理の概要(第2例)を模式的に示す図である。 図18は、本発明の実施の形態2の超音波観測装置のパラメータ抽出部が投影直線と直交する直線の傾きを算出する処理を模式的に示す図である。
以下、添付図面を参照して、本発明を実施するための形態(以下、「実施の形態」という)を説明する。
(実施の形態1)
図1は、本発明の実施の形態1に係る超音波観測装置を備えた超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム1は、観測対象である被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波内視鏡2と、超音波内視鏡2が取得した超音波信号に基づいて超音波画像を生成する超音波観測装置3と、超音波観測装置3が生成した超音波画像を表示する表示装置4と、を備える。
超音波内視鏡2は、その先端部に、超音波観測装置3から受信した電気的なパルス信号を超音波パルス(音響パルス)に変換して被検体へ照射するとともに、被検体で反射された超音波エコーを電圧変化で表現する電気的なエコー信号に変換して出力する超音波振動子21を有する。超音波振動子21は、コンベックス振動子、リニア振動子およびラジアル振動子のいずれでも構わない。超音波内視鏡2は、超音波振動子21をメカ的に走査させるものであってもよいし、超音波振動子21として複数の素子をアレイ状に設け、送受信にかかわる素子を電子的に切り替えたり、各素子の送受信に遅延をかけたりすることで、電子的に走査させるものであってもよい。
超音波内視鏡2は、通常は撮像光学系および撮像素子を有しており、被検体の消化管(食道、胃、十二指腸、大腸)または呼吸器(気管、気管支)へ挿入され、消化管、呼吸器やその周囲臓器(膵臓、胆嚢、胆管、胆道、リンパ節、縦隔臓器、血管等)を撮像することが可能である。また、超音波内視鏡2は、撮像時に被検体へ照射する照明光を導くライトガイドを有する。このライトガイドは、先端部が超音波内視鏡2の被検体への挿入部の先端まで達している一方、基端部が照明光を発生する光源装置に接続されている。
超音波観測装置3は、超音波内視鏡2と電気的に接続され、所定の波形および送信タイミングに基づいて高電圧パルスからなる送信信号(パルス信号)を超音波振動子21へ送信するとともに、超音波振動子21から電気的な受信信号であるエコー信号を受信してデジタルの高周波(RF:Radio Frequency)信号のデータ(以下、RFデータという)を生成、出力する送受信部31と、送受信部31から受信したRFデータをもとにデジタルのBモード用受信データを生成する信号処理部32と、送受信部31から受信したRFデータに対して所定の演算を施す演算部33と、各種画像データを生成する画像処理部34と、キーボード、マウス、タッチパネル等のユーザインタフェースを用いて実現され、各種情報の入力を受け付ける入力部35と、超音波診断システム1全体を制御する制御部36と、超音波観測装置3の動作に必要な各種情報を記憶する記憶部37と、を備える。
送受信部31は、エコー信号を増幅する信号増幅部311を有する。信号増幅部311は、受信深度が大きいエコー信号ほど高い増幅率で増幅するSTC(Sensitivity Time Control)補正を行う。図2は、信号増幅部311が行う増幅処理における受信深度と増幅率との関係を示す図である。図2に示す受信深度zは、超音波の受信開始時点からの経過時間に基づいて算出される量である。図2に示すように、増幅率β(dB)は、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴ってβからβth(>β0)へ線型に増加する。また、増幅率β(dB)は、受信深度zが閾値zth以上である場合、一定値βthをとる。閾値zthの値は、観測対象から受信する超音波信号がほとんど減衰してしまい、ノイズが支配的になるような値である。より一般に、増幅率βは、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴って単調増加すればよい。図2に示す関係は、予め記憶部37に記憶されている。
送受信部31は、信号増幅部311によって増幅されたエコー信号に対してフィルタリング等の処理を施した後、A/D変換することによって時間ドメインのRFデータを生成し、信号処理部32および演算部33へ出力する。超音波内視鏡2が複数の素子をアレイ状に設けた超音波振動子21を電子的に走査させる構成を有する場合、送受信部31は、複数の素子に対応したビーム合成用の多チャンネル回路を有する。
送受信部31が送信するパルス信号の周波数帯域は、超音波振動子21におけるパルス信号の超音波パルスへの電気音響変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。また、信号増幅部311におけるエコー信号の各種処理周波数帯域は、超音波振動子21による超音波エコーのエコー信号への音響電気変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。これらにより、後述する周波数スペクトルの近似処理を実行する際、精度のよい近似を行うことが可能となる。
送受信部31は、制御部36が出力する各種制御信号を超音波内視鏡2に対して送信するとともに、超音波内視鏡2から識別用のIDを含む各種情報を受信して制御部36へ送信する機能も有する。
信号処理部32は、RFデータに対してバンドパスフィルタ、包絡線検波、対数変換など公知の処理を施し、デジタルのBモード用受信データを生成する。対数変換では、RFデータを基準電圧で除した量の常用対数をとってデシベル値で表現する。信号処理部32は、生成したBモード用受信データを、画像処理部34へ出力する。信号処理部32は、CPU(Central Processing Unit)や各種演算回路等を用いて実現される。
演算部33は、送受信部31が出力したRFデータに対して受信深度によらず増幅率を一定とするよう増幅補正を行う増幅補正部331と、増幅補正を行ったRFデータに高速フーリエ変換(FFT:Fast Fourier Transform)を施して周波数解析を行うことにより超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部332と、周波数解析部332が算出した各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定部333と、周波数帯域設定部333が設定した周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出部334と、を有する。演算部33は、CPUや各種演算回路等を用いて実現される。
図3は、増幅補正部331が行う増幅補正処理における受信深度と増幅率との関係を示す図である。図3に示すように、増幅補正部331が行う増幅処理における増幅率β(dB)は、受信深度zがゼロのとき最大値βth−β0をとり、受信深度zがゼロから閾値zthに達するまで線型に減少し、受信深度zが閾値zth以上のときゼロである。なお、図3に示す関係は、予め記憶部37に記憶されている。増幅補正部331が図3に示す関係に基づいてデジタルRF信号を増幅補正することにより、信号増幅部311におけるSTC補正の影響を相殺し、一定の増幅率βthの信号を出力することができる。なお、増幅補正部331が行う受信深度zと増幅率βの関係は、信号増幅部311における受信深度と増幅率の関係に応じて異なることは勿論である。
このような増幅補正を行う理由を説明する。STC補正は、アナログ信号波形の振幅を全周波数帯域にわたって均一に、かつ、深度に対しては単調増加する増幅率で増幅させることで、アナログ信号波形の振幅から減衰の影響を排除する補正処理である。このため、エコー信号の振幅を輝度に変換して表示するBモード画像を生成する場合、かつ、一様な組織を走査した場合には、STC補正を行うことによって深度によらず輝度値が一定になる。すなわち、Bモード画像の輝度値から減衰の影響を排除する効果を得ることができる。
一方、本実施の形態1のように超音波の周波数スペクトルを算出して解析した結果を利用する場合、STC補正でも超音波の伝播に伴う減衰の影響を正確に排除できるとは限らない、という問題がある。なぜなら、一般に減衰量は周波数によって異なるが(後述する式(1)を参照)、STC補正の増幅率は距離だけに応じて変化し、周波数依存性がないためである。
上述した問題を解決するには、Bモード画像を生成する際にSTC補正を施した受信信号を出力する一方、周波数スペクトルに基づいた画像を生成する際に、Bモード画像を生成するための送信とは異なる新たな送信を行い、STC補正を施していない受信信号を出力することが考えられる。ところがこの場合には、受信信号に基づいて生成される画像データのフレームレートが低下してしまうという問題がある。
そこで、本実施の形態1では、生成される画像データのフレームレートを維持しつつ、Bモード画像用にSTC補正を施した信号に対してSTC補正の影響を排除するために、増幅補正部331によって増幅率の補正を行う。
周波数解析部332は、増幅補正部331が増幅補正した各音線のRFデータ(ラインデータ)を所定の時間間隔でサンプリングし、サンプルデータを生成する。周波数解析部332は、サンプルデータ群にFFT処理を施すことにより、RFデータ上の複数の箇所(データ位置)における周波数スペクトルを算出する。
図4は、超音波信号の1つの音線におけるデータ配列を模式的に示す図である。同図に示す音線SRkにおいて、白または黒の長方形は、1つのサンプル点におけるデータを意味している。また、音線SRkにおいて、右側に位置するデータほど、超音波振動子21から音線SRkに沿って計った場合の深い箇所からのサンプルデータである(図4の矢印を参照)。音線SRkは、送受信部31が行うA/D変換におけるサンプリング周波数(例えば50MHz)に対応した時間間隔で離散化されている。図4では、番号kの音線SRkの8番目のデータ位置を受信深度zの方向の初期値Z(k) 0として設定した場合を示しているが、初期値の位置は任意に設定することができる。周波数解析部332による算出結果は複素数で得られ、記憶部37に格納される。
図4に示すデータ群Fj(j=1,2,・・・,K)は、FFT処理の対象となるサンプルデータ群である。一般に、FFT処理を行うためには、サンプルデータ群が2のべき乗のデータ数を有している必要がある。この意味で、サンプルデータ群Fj(j=1,2,・・・,K−1)はデータ数が16(=24)で正常なデータ群である一方、サンプルデータ群FKは、データ数が12であるため異常なデータ群である。異常なデータ群に対してFFT処理を行う際には、不足分だけゼロデータを挿入することにより、正常なサンプルデータ群を生成する処理を行う。この点については、周波数解析部332の処理を説明する際に詳述する(図10を参照)。
周波数解析部332が算出する「周波数スペクトル」とは、サンプルデータ群にFFT処理を施すことによって得られた「ある受信深度zにおける強度の周波数分布」を意味する。また、ここでいう「強度」とは、例えばエコー信号の電圧、エコー信号の電力、超音波エコーの音圧、超音波エコーの音響エネルギー等のパラメータ、これらパラメータの振幅や時間積分値やその組み合わせのいずれかを指す。
一般に、周波数スペクトルは、観測対象が生体組織である場合、超音波が走査された生体組織の性状によって異なる傾向を示す。これは、周波数スペクトルが、超音波を散乱する散乱体の大きさ、数密度、音響インピーダンス等と相関を有しているためである。ここでいう「生体組織の性状」とは、例えば悪性腫瘍(癌)、良性腫瘍、内分泌腫瘍、粘液性腫瘍、正常組織、嚢胞、脈管などのことである。
周波数帯域設定部333は、周波数スペクトルの包絡線を検出する包絡線検出部333aと、包絡線検出部333aが検出した包絡線の最大値に対応する周波数を含む複数の候補帯域において包絡線を近似する近似曲線に関連するパラメータを候補帯域ごとに抽出するパラメータ抽出部333bと、パラメータ抽出部333bが候補帯域ごとに抽出した複数のパラメータのうち所定の関係にある2つの候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出する変化量算出部333cと、変化量算出部333cが算出したパラメータの変化量に基づいて周波数帯域を決定する決定部333dと、を有する。
包絡線検出部333aは、公知の包絡線検出処理を用いて周波数スペクトルの包絡線を検出する。図5は、包絡線検出部333aが検出した周波数スペクトルの包絡線の例を示す図である。同図に示す包絡線Cは、強度I0を基準強度Ic(定数)で除した量の常用対数(デシベル表現)I=10log10(I0/Ic)の周波数fに対する分布を示している。なお、本実施の形態1において、曲線および直線は、離散的な点の集合からなる。
パラメータ抽出部333bは、包絡線検出部333aが検出した包絡線を候補帯域で回帰分析することによって包絡線を直線近似し、近似した直線の傾きをパラメータとして抽出する。図6は、パラメータ抽出部333bが実行するパラメータ抽出処理の概要を模式的に示す図である。パラメータ抽出部333bは、包絡線Cに対して候補帯域Un={f|fs<f<fn}を近似する直線を回帰分析によって算出し、得られた近似直線(回帰直線)Lnの傾きをパラメータとする。候補帯域Unの開始周波数fsは、包絡線Cの最大値をとるピーク周波数fpよりも所定値だけ小さい周波数(定数)として設定される。また、候補帯域Unの最大周波数fn(n=0,1,2,・・・)はfn>fmaxであり、かつfn=fn-1+Δf(Δf>0)である。ここで、Δfは所定の定数であり、開始周波数fsよりも十分小さい値である。パラメータ抽出部333bは、決定部333dが周波数帯域を決定するまでnの値を増加させながら直線の傾きを算出する。なお、開始周波数fsおよび候補帯域の変更に伴う最大周波数の増加量Δfは任意に設定することが可能であり、ユーザが入力部35を介して設定入力するようにしてもよい。
変化量算出部333cは、候補帯域Un+1で算出したパラメータPn+1から候補帯域Unで算出したパラメータPnを減算した値Pn+1−Pnをパラメータの変化量ΔPnとして算出する。変化量算出部333cは、決定部333dが周波数帯域を決定するまでパラメータの変化量ΔPnを算出する。
決定部333dは、変化量算出部333cが算出したパラメータの変化量の符号に基づいて周波数帯域を決定する。図7は、決定部333dが実行する周波数帯域決定処理の概要を示す図である。図7において、横軸はパラメータの変化量を算出する2つの候補帯域(Un,Un+1)(n=0,1,・・・)を示しており、縦軸は2つの候補帯域(Un,Un+1)におけるパラメータの変化量ΔPn=Pn+1−Pnを示している。図7に示す場合、2つの候補帯域(Un,Un+1)を示すnの値がn=0,1,・・・,N−1である場合にはΔPn≦0である一方、nの値がn=Nである場合にはΔPn(=ΔPN)>0となっている。決定部333dは、パラメータの変化量ΔPが「0以下」から「正」に変化したときの2つの候補帯域(UN,UN+1)のうち帯域幅が小さい方の候補帯域UNを周波数帯域とする。以下、周波数帯域となった候補帯域UNの最大周波数fNを終了周波数feという。
特徴量算出部334は、回帰分析により周波数スペクトルの包絡線を直線で近似することによって周波数スペクトルの補正前特徴量を算出する近似部334aと、近似部334aが算出した補正前特徴量に対して複数の減衰率候補値の各々に基づいた減衰補正を行うことによって補正特徴量を算出する減衰補正部334bと、を有する。
近似部334aは、所定周波数帯域における周波数スペクトルの包絡線に対して回帰分析を行って一次式で近似することにより、この近似した一次式を特徴付ける補正前特徴量を算出する。例えば、図5に示す包絡線Cの場合、近似部334aは、周波数帯域設定部333が設定した周波数帯域(候補帯域)で回帰分析を行うことによって包絡線Cの近似直線を得る。図8は、近似部334aが実行する近似処理の概要を示す図である。図8に示す直線LU0は、候補帯域Uで包絡線Cを一次式I=a0f+b0で近似することによって得られた近似直線である。近似部334aは、直線LU0に対応する補正前特徴量として、傾きa0、切片b0、および候補帯域Uの中心周波数fM=(fs+fe)/2における強度Iの値であるミッドバンドフィット(Mid-band fit)c0=a0M+b0を算出する。なお、近似部334aは、回帰分析によって二次以上の多項式で周波数スペクトルを近似するようにしてもよい。
3つの補正前特徴量のうち、傾きa0は、超音波の散乱体の大きさと相関を有し、一般に散乱体が大きいほど傾きが小さな値を有すると考えられる。また、切片b0は、散乱体の大きさ、音響インピーダンスの差、散乱体の数密度(濃度)等と相関を有している。具体的には、切片b0は、散乱体が大きいほど大きな値を有し、音響インピーダンスの差が大きいほど大きな値を有し、散乱体の数密度が大きいほど大きな値を有すると考えられる。ミッドバンドフィットc0は、傾きa0と切片b0から導出される間接的なパラメータであり、有効な周波数帯域内の中心におけるスペクトルの強度を与える。このため、ミッドバンドフィットc0は、散乱体の大きさ、音響インピーダンスの差、散乱体の数密度に加えて、Bモード画像の輝度とある程度の相関を有していると考えられる。
減衰補正部334bは、単位長さおよび単位周波数あたりの超音波の減衰量を与える減衰率を用いて減衰補正を行う。一般に、超音波の減衰量A(f,z)は、超音波が受信深度0と受信深度zとの間を往復する間に生じる減衰であり、往復する前後の強度変化(デシベル表現での差)として定義される。減衰量A(f,z)は、一様な組織内では周波数に比例することが経験的に知られており、以下の式(1)で表現される。
A(f,z)=2αzf ・・・(1)
ここで、比例定数αは減衰率と呼ばれる量であり、単位長さおよび単位周波数あたりの超音波の減衰量を与える。また、zは超音波の受信深度であり、fは周波数である。減衰率αの具体的な値は、観測対象が生体である場合、生体の部位に応じて定まる。減衰率αの単位は、例えばdB/cm/MHzである。
減衰補正部334bは、近似部334aが抽出した補正前特徴量(傾きa0、切片b0、ミッドバンドフィットc0)に対し、以下に示す式(2)〜(4)にしたがって減衰補正を行うことにより、補正特徴量a、b、cを算出する。
a=a0+2αz ・・・(2)
b=b0 ・・・(3)
c=c0+A(fM,z)=c0+2αzfM(=afM+b) ・・・(4)
式(2)、(4)からも明らかなように、減衰補正部334bは、超音波の受信深度zが大きいほど、補正量が大きい補正を行う。また、式(3)によれば、切片に関する補正は恒等変換である。これは、切片が周波数0(Hz)に対応する周波数成分であって減衰の影響を受けないためである。
図8に示す直線LUは、減衰補正部334bが補正した補正特徴量a、b、cをパラメータとして有する直線を示す図である。直線LUの式は、
I=af+b=(a0+2αz)f+b0 ・・・(5)
で表される。この式(5)からも明らかなように、直線LUは、減衰補正前の直線LU0と比較して、傾きが大きく(a>a0)、かつ切片が同じ(b=b0)である。
画像処理部34は、エコー信号の振幅を輝度に変換して表示する超音波画像であるBモード画像データを生成するBモード画像データ生成部341と、特徴量算出部334が算出した特徴量に関する情報を表示する特徴量画像データを生成する特徴量画像データ生成部342と、を有する。
Bモード画像データ生成部341は、信号処理部32から受信したBモード用受信データに対してゲイン処理、コントラスト処理等の公知の技術を用いた信号処理を行うとともに、表示装置4における画像の表示レンジに応じて定まるデータステップ幅に応じたデータの間引き等を行うことによってBモード画像データを生成する。Bモード画像は、色空間としてRGB表色系を採用した場合の変数であるR(赤)、G(緑)、B(青)の値を一致させたグレースケール画像である。
Bモード画像データ生成部341は、Bモード用受信データに走査範囲を空間的に正しく表現できるよう並べ直す座標変換を施した後、Bモード用受信データ間の補間処理を施すことによってBモード用受信データ間の空隙を埋め、Bモード画像データを生成する。Bモード画像データ生成部341は、生成したBモード画像データを特徴量画像データ生成部342へ出力する。
特徴量画像データ生成部342は、特徴量算出部334が算出した特徴量に関連する視覚情報をBモード画像データにおける画像の各画素に対して重畳することによって特徴量画像データを生成する。特徴量画像データ生成部342は、例えば図4に示す1つのサンプルデータ群Fj(j=1、2、・・・、K)のデータ量に対応する画素領域に対し、その振幅データ群Fjから算出される周波数スペクトルの特徴量に対応する視覚情報を割り当てる。特徴量画像データ生成部342は、例えば上述した傾き、切片、ミッドバンドフィットのいずれか一つに視覚情報としての色相を対応付けることによって特徴量画像データを生成する。なお、特徴量画像データ生成部342が、傾き、切片、ミッドバンドフィットから選択される2つの特徴量の一方に色相を対応付けるとともに、他方に明暗を対応付けることによって特徴量画像データを生成するようにしてもよい。特徴量に関連する視覚情報としては、色相や明暗(明度)のほか、例えば彩度、輝度値、R(赤)、G(緑)、B(青)などの所定の表色系を構成する色空間の変数を挙げることができる。
制御部36は、演算および制御機能を有するCPU(Central Processing Unit)や各種演算回路等を用いて実現される。制御部36は、記憶部37が記憶、格納する情報を記憶部37から読み出し、超音波観測装置3の作動方法に関連した各種演算処理を実行することによって超音波観測装置3を統括して制御する。なお、制御部36を信号処理部32および演算部33と共通のCPU等を用いて構成することも可能である。
記憶部37は、周波数解析部332が算出した周波数スペクトルの情報を受信深度および受信方向とともに記憶するスペクトル情報記憶部371と、パラメータ抽出部333bが候補帯域ごとに抽出したパラメータの情報を記憶するパラメータ情報記憶部372と、周波数帯域設定部333が設定した周波数帯域の情報を記憶する周波数帯域情報記憶部373と、を有する。
記憶部37は、上記以外にも、例えば増幅処理に必要な情報(図2に示す増幅率と受信深度との関係)、増幅補正処理に必要な情報(図3に示す増幅率と受信深度との関係)、減衰補正処理に必要な情報(式(1)参照)、周波数解析処理に必要な窓関数(Hamming、Hanning、Blackman等)の情報等を記憶する。
また、記憶部37は、超音波観測装置3の作動方法を実行するための作動プログラムを含む各種プログラムを記憶する。作動プログラムは、ハードディスク、フラッシュメモリ、CD−ROM、DVD−ROM、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体に記録して広く流通させることも可能である。なお、上述した各種プログラムは、通信ネットワークを介してダウンロードすることによって取得することも可能である。ここでいう通信ネットワークは、例えば既存の公衆回線網、LAN(Local Area Network)、WAN(Wide Area Network)などによって実現されるものであり、有線、無線を問わない。
以上の構成を有する記憶部37は、各種プログラム等が予めインストールされたROM(Read Only Memory)、および各処理の演算パラメータやデータ等を記憶するRAM(Random Access Memory)等を用いて実現される。
図9は、以上の構成を有する超音波観測装置3が実行する処理の概要を示すフローチャートである。具体的には、超音波観測装置3が超音波内視鏡2からエコー信号を受信する以降の処理の概要を示すフローチャートである。以下、図9を参照して、超音波観測装置3が行う処理を説明する。まず、超音波観測装置3は、超音波内視鏡2から超音波振動子21による観測対象の測定結果としてのエコー信号を受信する(ステップS1)。
超音波振動子21からエコー信号を受信した信号増幅部311は、そのエコー信号の増幅を行う(ステップS2)。ここで、信号増幅部311は、例えば図2に示す増幅率と受信深度との関係に基づいてエコー信号の増幅(STC補正)を行う。
続いて、Bモード画像データ生成部341は、信号増幅部311が増幅したエコー信号を用いてBモード画像データを生成して、表示装置4へ出力する(ステップS3)。Bモード画像データを受信した表示装置4は、そのBモード画像データに対応するBモード画像を表示する。
増幅補正部331は、送受信部31から出力されたRFデータに対して受信深度によらず増幅率が一定となるように増幅補正を行う(ステップS4)。ここで、増幅補正部331は、例えば図3に示す増幅率と受信深度との関係に基づいて増幅補正を行う。
この後、周波数解析部332は、増幅補正後の各音線のRFデータに対してFFTによる周波数解析を行うことによって全てのサンプルデータ群に対する周波数スペクトルを算出し、スペクトル情報記憶部371へ格納する(ステップS5)。図10は、ステップS5において周波数解析部332が実行する処理の概要を示すフローチャートである。以下、図10に示すフローチャートを参照して、周波数解析処理を詳細に説明する。
まず、周波数解析部332は、解析対象の音線を識別するカウンタkをk0とする(ステップS11)。
続いて、周波数解析部332は、FFT演算用に生成した一連のデータ群(サンプルデータ群)を代表するデータ位置(受信深度に相当)Z(k)の初期値Z(k) 0を設定する(ステップS12)。例えば、図4では、上述したように、音線SRkの8番目のデータ位置を初期値Z(k) 0として設定した場合を示している。
その後、周波数解析部332は、サンプルデータ群を取得し(ステップS13)、取得したサンプルデータ群に対し、記憶部37が記憶する窓関数を作用させる(ステップS14)。このようにサンプルデータ群に対して窓関数を作用させることにより、サンプルデータ群が境界で不連続になることを回避し、アーチファクトが発生するのを防止することができる。
続いて、周波数解析部332は、データ位置Z(k)のサンプルデータ群が正常なデータ群であるか否かを判定する(ステップS15)。図4を参照した際に説明したように、サンプルデータ群は、2のべき乗のデータ数を有している必要がある。以下、正常なサンプルデータ群のデータ数を2n(nは正の整数)とする。本実施の形態では、データ位置Z(k)が、できるだけZ(k)が属するサンプルデータ群の中心になるよう設定される。具体的には、サンプルデータ群のデータ数は2nであるので、Z(k)はそのサンプルデータ群の中心に近い2n/2(=2n-1)番目の位置に設定される。この場合、サンプルデータ群が正常であるとは、データ位置Z(k)より浅部側に2n-1−1(=Nとする)個のデータがあり、データ位置Z(k)より深部側に2n-1(=Mとする)個のデータがあることを意味する。図4に示す場合、サンプルデータ群Fj(j=1、2、・・・、K−1)は正常である。なお、図4ではn=4(N=7,M=8)の場合を例示している。
ステップS15における判定の結果、データ位置Z(k)のサンプルデータ群が正常である場合(ステップS15:Yes)、周波数解析部332は、後述するステップS17へ移行する。
ステップS15における判定の結果、データ位置Z(k)のサンプルデータ群が正常でない場合(ステップS15:No)、周波数解析部332は、不足分だけゼロデータを挿入することによって正常なサンプルデータ群を生成する(ステップS16)。ステップS15において正常でないと判定されたサンプルデータ群(例えば図4のサンプルデータ群FK)は、ゼロデータを追加する前に窓関数が作用されている。このため、サンプルデータ群にゼロデータを挿入してもデータの不連続は生じない。ステップS16の後、周波数解析部332は、後述するステップS17へ移行する。
ステップS17において、周波数解析部332は、サンプルデータ群を用いてFFT演算を行うことにより、振幅の周波数分布である周波数スペクトルを得る(ステップS17)。
続いて、周波数解析部332は、データ位置Z(k)をステップ幅Dで変化させる(ステップS18)。ステップ幅Dは、記憶部37が予め記憶しているものとする。図4では、D=15の場合を例示している。ステップ幅Dは、Bモード画像データ生成部341がBモード画像データを生成する際に利用するデータステップ幅と一致させることが望ましいが、周波数解析部332における演算量を削減したい場合には、ステップ幅Dとしてデータステップ幅より大きい値を設定してもよい。
その後、周波数解析部332は、データ位置Z(k)が音線SRkにおける最大値Z(k) maxより大きいか否かを判定する(ステップS19)。データ位置Z(k)が最大値Z(k) maxより大きい場合(ステップS19:Yes)、周波数解析部332はカウンタkを1増加させる(ステップS20)。これは、処理をとなりの音線へ移すことを意味する。一方、データ位置Z(k)が最大値Z(k) max以下である場合(ステップS19:No)、周波数解析部332はステップS13へ戻る。
ステップS20の後、周波数解析部332は、カウンタkが最大値kmaxより大きいか否かを判定する(ステップS21)。カウンタkが最大値kmaxより大きい場合(ステップS21:Yes)、周波数解析部332は一連の周波数解析処理を終了する。一方、カウンタkが最大値kmax以下である場合(ステップS21:No)、周波数解析部332はステップS12に戻る。この最大値kmaxは、術者等のユーザが入力部35を通じて任意に指示入力した値、もしくは、記憶部37にあらかじめ設定された値とする。
このようにして、周波数解析部332は、解析対象領域内の(kmax−k0+1)本の音線の各々について複数回のFFT演算を行う。FFT演算の結果として得られる周波数スペクトルは、受信深度、受信方向および識別情報とともにスペクトル情報記憶部371に格納される。ここで識別情報は、後述する周波数帯域設定処理の際に周波数スペクトルを識別するための情報であり、例えば各周波数スペクトルに対して順次付される識別番号i=0,1,2,・・・,imaxを含む情報である。本実施の形態1では、この識別番号iが、後述する周波数帯域設定処理におけるカウンタの機能を有している。
なお、以上の説明では、周波数解析部332が超音波信号を受信したすべての領域に対して周波数解析処理を行うものとしたが、入力部35が特定の深度幅および音線幅で区切られる部分領域の設定入力を受け付け可能な構成とし、設定された部分領域内においてのみ周波数解析処理を行うようにすることも可能である。
以上説明したステップS5の周波数解析処理に続いて、周波数帯域設定部333は、周波数スペクトルごとに特徴量算出用の周波数帯域を設定する(ステップS6)。図11は、このステップS6における周波数帯域設定処理の概要を示すフローチャートである。以下、図11を参照して、周波数帯域設定処理を詳細に説明する。
まず、周波数帯域設定部333は、周波数スペクトルを識別するためのカウンタiを0に設定する(ステップS31)。
続いて、包絡線検出部333aは、カウンタiに対応する周波数スペクトルをスペクトル情報記憶部371から取得し、取得した周波数スペクトルの包絡線を検出する(ステップS32)。この結果、例えば図5に示す包絡線Cが検出される。
この後、パラメータ抽出部333bは、開始周波数fsを設定する(ステップS33)。具体的には、パラメータ抽出部333bは、処理対象の包絡線の強度のピーク値に対応するピーク周波数fpを抽出し、ピーク周波数fpよりも所定値だけ小さい周波数を開始周波数fsとする。なお、このステップS33は、あらかじめユーザが入力部35を介して設定入力した開始周波数fsの値を記憶部37から読み出して設定する処理であってもよい。
パラメータ抽出部333bは、繰り返しのカウンタnを0と設定する(ステップS34)。続いて、パラメータ抽出部333bは、候補帯域Un={f|fs<f<fn}における包絡線の近似直線を回帰分析によって算出し、その近似直線の傾きをパラメータPnとして抽出し、パラメータ情報記憶部372に書き込んで記憶させる(ステップS35)。
その後、パラメータ抽出部333bは、カウンタnが0であるか否かを判定する(ステップS36)。カウンタnが0である場合(ステップS36:Yes)、パラメータ抽出部333bは、カウンタnを1増加させ(ステップS37)、ステップS35の処理へ戻る。ここで、カウンタnが1増加すると、候補帯域の最大周波数の値は、上述した通りΔfだけ増加する。カウンタnが0でない場合(ステップS36:No)、パラメータ抽出部333bは、ステップS38の処理へ移行する。
ステップS38において、変化量算出部333cは、パラメータの変化量として、2つの候補帯域(Un-1,Un)における近似直線の傾きの差を算出する(ステップS38)。具体的には、変化量算出部333cは、パラメータ抽出部333bが算出した候補帯域Unにおける近似直線の傾き(パラメータPn)から、その候補帯域Un-1における近似直線の傾き(パラメータPn-1)を減算してパラメータの変化量ΔPn=Pn+1−Pnを算出し、算出結果をパラメータ情報記憶部372に書き込んで記憶させる。
決定部333dは、変化量算出部333cが算出したパラメータの変化量の符号が「0以下」から「正」に変化したか否かを判定する(ステップS39)。具体的には、決定部333dは、変化量算出部333cが算出した最新のパラメータの変化量ΔPnの符号が「正」である場合、その一つ前に変化量算出部333cが算出したパラメータの変化量ΔPn-1の符号が「0以下」であるか否かを確認することにより、パラメータの変化量の符号が「0以下」から「正」に変化したか否かを判定する。判定の結果、パラメータの変化量の符号が「0以下」から「正」に変化した場合(ステップS39:Yes)、決定部333dは、候補帯域Un-1を周波数帯域Uとする(ステップS40)。例えば、図7では、候補帯域UN(N=n−1)が周波数帯域Uとして設定される場合を示している。一方、判定の結果、パラメータの変化量の符号が「0以下」から「正」に変化していない場合(ステップS39:No)、周波数帯域設定部333は、ステップS37の処理へ移行する。
この後、周波数帯域設定部333は、カウンタiが所定値imaxより大きいか否かを判定する(ステップS41)。判定の結果、i>imaxである場合(ステップS41:Yes)、周波数帯域設定部333は、周波数帯域設定処理を終了する。一方、判定の結果、i≦imaxである場合(ステップS41:No)、周波数帯域設定部333は、カウンタiを1増加させ(ステップS42)、ステップS32の処理に戻る。
以上説明したステップS6の周波数帯域設定処理の後、近似部334aは、周波数帯域U={f|fs<f≦fe}における近似直線をもとに補正前特徴量を算出する(ステップS7)。補正前特徴量としては、例えば上述した近似直線の傾き、切片、およびミッドバンドフィットの少なくともいずれか一つを挙げることができる。
続いて、減衰補正部334bは、補正前特徴量に対して減衰補正を行う(ステップS8)。これにより、超音波の伝播における減衰の影響を排除した特徴量が算出されることとなる。図8に示す場合には、この減衰補正処理の結果、直線LUが得られる。
この後、特徴量画像データ生成部342は、ステップS8で算出した特徴量とステップS3で生成したBモード画像データを用いて特徴量画像データを生成する(ステップS9)。特徴量画像データ生成部342は、生成した特徴量画像データを表示装置4へ送信する。特徴量画像データを受信した表示装置4は、受信した特徴量画像データに対応する特徴量画像を表示する。
ステップS9の後、超音波観測装置3は一連の処理を終了する。なお、超音波観測装置3は、ステップS1〜S9の処理を周期的に繰り返し実行する。
図12は、本実施の形態1の効果を説明する図である。図12において、包絡線C’は、包絡線Cと同じ受信深度で異なる組織性状を有する組織の周波数スペクトルから検出されたものである。具体的には、包絡線Cの組織の粒径が包絡線C’の組織の粒径より大きい場合を示している。直線LU0’は、包絡線C’に対して設定される周波数帯域U’={f|fs<f≦fe’}に対して近似部334aが算出した近似直線である。本実施の形態1では、特徴量算出用の周波数帯域を適応的に決定するため、包絡線Cで設定される周波数帯域Uと包絡線C’で設定される周波数帯域U’は異なる。その結果、2つの近似直線LU0、LU0’の各周波数帯域U、U’における傾きは明らかに異なっている。
これに対して、従来技術では、同じ受信深度であれば同じ周波数帯域を用いる。図13は、比較のために図12と同じ包絡線C、C’に対して従来技術を適用した場合に得られる近似直線を例示した図である。図13に示す場合、包絡線C、C’に対して同じ周波数帯域U(fix)={f|fs<f≦fe(fix)}で近似直線を求めている。その結果、包絡線Cの近似直線LU0(fix)の傾きと、包絡線C’の近似直線LU0(fix)’の傾きに大きな差が生じず、2つの組織を明確に峻別することが不可能である。
このように、本発明の実施の形態1によれば、周波数スペクトルを近似する際、解析に最適な周波数帯域を適応的に決定するため、周波数帯域を少なくとも受信深度ごとに固定している従来技術と比較して、同一の受信深度であっても組織性状が異なる場合を明確に峻別することができる。したがって、本実施の形態1によれば、観測対象の組織性状を精度よく鑑別することが可能となる。
また、本実施の形態1によれば、複数の候補帯域において包絡線を近似する近似直線に関連するパラメータとして近似直線の傾きをそれぞれ抽出し、候補帯域ごとに抽出した複数のパラメータのうち候補帯域の帯域幅の値が隣接する2つの候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出し、この変化量の符号の変化に基づいて周波数帯域を決定することにより、組織性状に応じて周波数帯域をきめ細かく設定することができる。
なお、本実施の形態1において、パラメータ抽出部333bが、包絡線検出部333aが検出した包絡線を候補帯域で回帰分析することによって包絡線を直線近似し、近似した直線の切片をパラメータとして抽出するようにしてもよい。
(実施の形態1の変形例1)
図14は、本実施の形態1の変形例1における周波数帯域設定処理の概要を示す図である。本変形例1において、パラメータ抽出部333bは、複数の候補帯域Un(n=1,2,・・・,nmax)に対してパラメータPnを抽出する。続いて、変化量算出部333cは、複数の候補帯域Unに対して実施の形態1と同様に2つの候補領域(Un,Un+1)におけるパラメータの変化量ΔPnを算出する。決定部333dは、パラメータの変化量ΔPnの絶対値|ΔPn|を2つの候補帯域のうち帯域幅が小さい方の候補帯域Unの最大周波数fnの関数として近似曲線を求める。この近似曲線は、例えば回帰分析を行うことによって得られる。図14に示す曲線K1は、本変形例1における近似曲線の例である。決定部333dは、曲線K1において、絶対値|ΔPn|の最小値|ΔPnminにおける周波数fminを終了周波数feとすることによって周波数帯域を決定する。
本変形例1によれば、実施の形態1と同様の効果を得ることができる。また、本変形例1によれば、一段ときめ細かく周波数帯域を設定することができるため、組織性状の鑑別精度をさらに向上させることができる。
(実施の形態1の変形例2)
図15は、本実施の形態1の変形例2における周波数帯域設定処理の概要を示す図である。本変形例2において、パラメータ抽出部333bは、複数の候補帯域Un(n=1,2,・・・,nmax)に対してパラメータPnを抽出する。続いて、変化量算出部333cは、実施の形態1と同様に2つの候補領域(Un,Un+1)におけるパラメータの変化量ΔPnを算出する。その後、決定部333dは、パラメータの変化量ΔPnを2つの候補帯域のうち帯域幅が小さい方の候補帯域Unの最大周波数fnの関数として近似曲線を求める。この近似曲線も、例えば回帰分析を行うことによって得られる。図15に示す曲線K2は本変形例2における近似曲線の例である。決定部333dは、曲線K2において、ΔPn=0となる周波数fminを終了周波数feとすることによって周波数帯域を設定する。
本変形例2によれば、実施の形態1と同様の効果を得ることができる。また、本変形例2によれば、変形例1と同様、一段ときめ細かく周波数帯域を設定することができるため、組織性状の鑑別精度をさらに向上させることができる。
(実施の形態2)
本発明の実施の形態2は、パラメータ抽出部333bが行うパラメータの算出方法が実施の形態1と異なる。具体的には、本実施の形態2において、パラメータ抽出部333bは、包絡線と同一平面内に設定される複数の投影直線へサンプリング点を投影した場合の度数分布を求め、求めた度数分布における統計的なばらつきが最小である投影直線を決定し、この投影直線と直交する直線の傾きを包絡線のパラメータとして抽出する。なお、本実施の形態2に係る超音波観測装置の構成は、実施の形態1で説明した超音波観測装置3の構成と同様である。
図16は、本実施の形態2においてパラメータ抽出部333bが行う包絡線のパラメータ抽出処理(図11のステップS35)の概要を示すフローチャートである。図17Aおよび図17Bは、パラメータ抽出部333bが行うパラメータ抽出処理の概要を模式的に示す図である。以下、図16、17Aおよび17Bを参照して包絡線のパラメータ抽出処理を説明する。
まず、パラメータ抽出部333bは、投影直線を識別するためのカウンタμを0と設定する(ステップS51)。
パラメータ抽出部333bは、角度パラメータθμの投影直線を作成する(ステップS52)。この角度パラメータθμは、包絡線を与える平面(横軸:周波数f、縦軸:強度I)における投影直線の傾きに関連した角度である。換言すれば、投影直線の傾きはtanθμで与えられる。図17Aには、カウンタμ=μ1(角度パラメータθμ=θμ1)である場合の投影直線FLμ1が記載されている。また、図17Bには、カウンタμ=μ2(角度パラメータθμ=θμ2)である場合の投影直線FLμ2が記載されている。ここで、μ1≠μ2であり、カウンタμの最大値をμmaxとすると、0≦μ1,μ2≦μmaxを満たす。
続いて、パラメータ抽出部333bは、候補帯域Unの包絡線Cのサンプリング点を投影直線に投影して投影直線上における度数分布を取得する(ステップS53)。図17Aに示す曲線Hμ1は、パラメータ抽出部333bが投影直線FLμ1に対して取得した度数分布に対応するヒストグラムである。また、図17Bに示す曲線Hμ2は、パラメータ抽出部333bが投影直線FLμ2に対して取得した度数分布に対応するヒストグラムである。
この後、パラメータ抽出部333bは、ステップS53で取得した度数分布の分散を算出し、パラメータ情報記憶部372へ書き込んで記憶させる(ステップS54)。
この後、パラメータ抽出部333bは、カウンタμが最大値μmaxより大きいか否かを判定する(ステップS55)。判定の結果、μ>μmaxである場合(ステップS55:Yes)、パラメータ抽出部333bはステップS56の処理へ移行する。一方、判定の結果、μ≦μmaxである場合(ステップS55:No)、パラメータ抽出部333bは、カウンタμを1増加させ(ステップS57)、ステップS52の処理に戻る。ここでパラメータ抽出部333bがカウンタμを1増加させたときの角度パラメータの増加量Δθはμによらず一定である。この増加量Δθは任意に設定することが可能であり、ユーザが入力部35を介して設定入力するようにしてもよい。
ステップS56において、パラメータ抽出部333bは、パラメータ情報記憶部372を参照して、サンプリング点を投影したときの度数分布の統計的ばらつきを与える分散が最小である投影直線を抽出する(ステップS56)。続いて、パラメータ抽出部333bは、ステップS56で抽出した投影直線と直交する直線の傾きを算出し、包絡線CのパラメータPnとしてパラメータ情報記憶部372へ書き込んで記憶させる(ステップS58)。
図18は、このステップS58の処理を模式的に示す図である。図18では、度数分布の分散が最小である投影直線が投影直線FLμ2である場合を例示している。この場合、パラメータ抽出部333bは、投影直線FLμ2と直交する直線Lμ2の傾きをパラメータPnとして抽出する。直線Lμ2の傾きは、角度パラメータθμ2を用いてtan(θμ2+90°)と表される。
以上説明した本発明の実施の形態2によれば、実施の形態1と同様、周波数スペクトルを近似する際、解析に最適な周波数帯域を適応的に決定するため、周波数帯域を少なくとも受信深度ごとに固定している従来技術と比較して、同一の受信深度であっても組織性状が異なる場合を明確に峻別することができる。したがって、本実施の形態2によれば、観測対象の組織性状を精度よく鑑別することが可能となる。
また、本実施の形態2においても、複数の候補帯域において包絡線を近似する近似直線のパラメータをそれぞれ抽出し、候補帯域ごとに抽出した複数のパラメータのうち候補帯域の帯域幅の値が隣接する2つの候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出し、この変化量の符号の変化に基づいて周波数帯域を決定することにより、組織性状に応じて周波数帯域をきめ細かく設定することができる。
なお、本実施の形態2において、候補帯域内のサンプリング点を投影直線に投影したときの度数分布の統計的なばらつきを与える量として、例えば標準偏差、母集団における特徴量の最大値と最小値の差、特徴量の分布の半値幅のいずれかを適用することも可能である。また、統計的なばらつきを与える量として分散の逆数を適用する場合も考えられるが、この場合には、その値が最大となる投影直線を用いてパラメータが決定されることはいうまでもない。
また、本実施の形態2において、パラメータ抽出部333bが、包絡線と同一平面内に設定される複数の投影直線へサンプリング点を投影した場合の度数分布を求め、求めた度数分布における統計的なばらつきが最小である投影直線を決定し、この投影直線と直交するとともにその投影直線上で最も度数が高い点を通過する直線の切片を包絡線のパラメータとして抽出するようにしてもよい。図18に示す直線Lμ2は、この2つの条件を満たす直線である。
また、本実施の形態2において、度数分布の統計的なばらつきに加えて該度数分布に対応するヒストグラムの最頻値を用いて定義される指標であって最頻値が大きいほど大きい値をとるとともに統計的なばらつきが小さいほど大きい値をとる指標に基づいてパラメータを算出する投影直線を設定するようにしてもよい。具体的には、統計的なばらつきを分散とする場合、(最頻値)/(分散)を指標として投影直線ごとに算出し、この指標の値が最大となる投影直線を用いて包絡線のパラメータを算出するようにしてもよい。
(その他の実施の形態)
ここまで、本発明を実施するための形態を説明してきたが、本発明は上述した実施の形態1、2によってのみ限定されるべきものではない。例えば、周波数スペクトルに対して減衰補正を先に行った後、周波数帯域設定処理および特徴量算出処理を行うようにしてもよい。この場合には、特徴量算出処理において減衰補正処理を行わなくてよい。
また、超音波内視鏡以外の超音波プローブに対しても適用することが可能である。超音波プローブとして、例えば光学系のない細径の超音波ミニチュアプローブを適用してもよい。超音波ミニチュアプローブは、通常、胆道、胆管、膵管、気管、気管支、尿道、尿管へ挿入され、その周囲臓器(膵臓、肺、前立腺、膀胱、リンパ節等)を観察する際に用いられる。また、超音波プローブとして、被検体の体表から超音波を照射する体外式超音波プローブを適用してもよい。体外式超音波プローブは、通常、腹部臓器(肝臓、胆嚢、膀胱)、乳房(特に乳腺)、甲状腺を観察する際に用いられる。
このように、本発明は、請求の範囲に記載した技術的思想を逸脱しない範囲内において、様々な実施の形態を含みうるものである。
以上のように、本発明にかかる超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラムは、観測対象の組織性状を精度よく鑑別するのに有用である。
1 超音波診断システム
2 超音波内視鏡
3 超音波観測装置
4 表示装置
21 超音波振動子
31 送受信部
32 信号処理部
33 演算部
34 画像処理部
35 入力部
36 制御部
37 記憶部
311 信号増幅部
331 増幅補正部
332 周波数解析部
333 周波数帯域設定部
333a 包絡線検出部
333b パラメータ抽出部
333c 変化量算出部
333d 決定部
334 特徴量算出部
334a 近似部
334b 減衰補正部
341 Bモード画像データ生成部
342 特徴量画像データ生成部
371 スペクトル情報記憶部
372 パラメータ情報記憶部
373 周波数帯域情報記憶部

Claims (10)

  1. 観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置であって、
    前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部と、
    前記周波数解析部が算出した各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定部と、
    前記周波数帯域設定部が設定した周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出部と、
    前記特徴量算出部が算出した特徴量に関する情報を表示する特徴量画像データを生成する特徴量画像データ生成部と、
    を備え
    前記周波数帯域設定部は、
    前記周波数スペクトルの包絡線を検出する包絡線検出部と、
    前記複数の候補帯域において前記包絡線を近似する近似直線に関連するパラメータをそれぞれ抽出するパラメータ抽出部と、
    前記パラメータ抽出部が前記候補帯域ごとに抽出した複数のパラメータのうち所定の関係にある2つの前記候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出する変化量算出部と、
    前記変化量算出部が算出した前記変化量に基づいて前記周波数帯域を決定する決定部と、
    を有することを特徴とする超音波観測装置。
  2. 前記パラメータ抽出部は、前記候補帯域の周波数の最小値を固定したまま最大値を順次大きくして前記複数の候補帯域を設定し、1つの前記候補帯域を設定するごとに前記包絡線を直線で近似してパラメータを抽出した後、前記変化量算出部へ前記パラメータを出力し、
    前記変化量算出部は、前記パラメータ抽出部が少なくとも2回のパラメータ抽出処理を実行済みである場合には、最新のパラメータとその1つ前のパラメータ抽出処理で抽出されたパラメータとの差を前記変化量として算出し、
    前記決定部は、前記変化量算出部が少なくとも2回の変化量算出処理を実行済みである場合には、最新の前記変化量とその1つ前の変化量算出処理で算出された前記変化量とを比較し、該比較の結果、前記変化量の符号が0または一方の符号から他方の符号に変化した場合、前記1つ前のパラメータ抽出処理で前記パラメータを抽出した際の前記候補帯域を前記周波数帯域に決定し、
    前記パラメータ抽出部、前記変化量算出部および前記決定部は、前記決定部が前記周波数帯域を決定するまで処理を繰り返し実行することを特徴とする請求項に記載の超音波観測装置。
  3. 前記パラメータ抽出部は、
    前記包絡線を回帰分析することによって前記直線を算出し、該算出した前記直線の傾きを前記パラメータとすることを特徴とする請求項に記載の超音波観測装置。
  4. 前記パラメータ抽出部は、前記候補帯域の周波数の最小値を固定したまま最大値を順次大きくして前記複数の候補帯域を設定し、各候補帯域において前記包絡線を直線で近似してパラメータを抽出し、
    前記変化量算出部は、帯域幅の値が隣接する2つの前記候補帯域においてそれぞれ抽出された2つの前記パラメータの差を前記変化量として算出し、
    前記決定部は、前記変化量の絶対値を2つの前記候補帯域のうち帯域幅が小さい方の候補帯域の最大周波数の関数として近似曲線を算出し、該近似曲線において前記絶対値が最小となるときの前記最大周波数をもとに前記周波数帯域を決定することを特徴とする請求項に記載の超音波観測装置。
  5. 前記パラメータ抽出部は、前記候補帯域の周波数の最小値を固定したまま最大値を順次大きくして前記複数の候補帯域を設定し、各候補帯域において前記包絡線を直線で近似してパラメータを抽出し、
    前記変化量算出部は、帯域幅の値が隣接する2つの前記候補帯域においてそれぞれ抽出された2つの前記パラメータの差を前記変化量として算出し、
    前記決定部は、前記変化量を2つの前記候補帯域のうち帯域幅が小さい方の候補帯域の最大周波数の関数として近似曲線を算出し、該近似曲線において値がゼロとなるときの前記最大周波数をもとに前記周波数帯域を決定することを特徴とする請求項に記載の超音波観測装置。
  6. 前記パラメータ抽出部は、
    前記候補帯域における前記包絡線上のサンプリング点を複数の投影直線にそれぞれ投影したときに各投影直線における度数分布を算出し、該度数分布の統計的なばらつきが最も小さい投影直線に基づいて前記直線の傾きを算出し、該算出した前記直線の傾きを前記パラメータとすることを特徴とする請求項に記載の超音波観測装置。
  7. 前記パラメータ抽出部は、
    前記候補帯域における前記包絡線上のサンプリング点を複数の投影直線にそれぞれ投影したときに各投影直線における度数分布を算出し、該度数分布の最頻値および統計的なばらつきを用いて定義される指標であって前記最頻値が大きいほど大きい値をとるとともに前記統計的なばらつきが小さいほど大きい値をとる指標に基づいて前記直線の傾きを算出する際に用いる投影直線を選択し、該選択した投影直線に基づいて前記直線の傾きを算出し、該算出した前記直線の傾きを前記パラメータとすることを特徴とする請求項に記載の超音波観測装置。
  8. 前記特徴量算出部は、
    前記周波数帯域における前記包絡線を回帰分析によって近似した後、超音波が伝搬する際に発生する減衰の寄与を削減する減衰補正を行うことによって前記周波数スペクトルの特徴量を算出することを特徴とする請求項1に記載の超音波観測装置。
  9. 観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置の作動方法であって、
    周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、
    周波数帯域設定部が、各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定ステップと、
    特徴量算出部が、前記周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出ステップと、
    特徴量画像データ生成部が、前記特徴量に基づいて特徴量画像データを生成する特徴量画像データ生成ステップと、
    含み、
    前記周波数帯域設定ステップは、
    包絡線検出部が、前記周波数スペクトルの包絡線を検出する包絡線検出ステップと、
    パラメータ抽出部が、前記複数の候補帯域において前記包絡線を近似する近似直線に関連するパラメータをそれぞれ抽出するパラメータ抽出ステップと、
    変化量算出部が、前記候補帯域ごとに抽出された複数のパラメータのうち所定の関係にある2つの前記候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出する変化量算出ステップと、
    決定部が、前記変化量に基づいて前記周波数帯域を決定する決定ステップと、
    を有することを特徴とする超音波観測装置の作動方法。
  10. 観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置に、
    周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、
    周波数帯域設定部が、各周波数スペクトルの特徴量の算出に使用する周波数帯域を、互いに異なる帯域幅を有する複数の候補帯域の各々における周波数スペクトルに関連するパラメータを用いて設定する周波数帯域設定ステップと、
    特徴量算出部が、前記周波数帯域をもとに各周波数スペクトルの特徴量を算出する特徴量算出ステップと、
    特徴量画像データ生成部が、前記特徴量に基づいて特徴量画像データを生成する特徴量画像データ生成ステップと、
    を実行させ
    前記周波数帯域設定ステップは、
    包絡線検出部が、前記周波数スペクトルの包絡線を検出する包絡線検出ステップと、
    パラメータ抽出部が、前記複数の候補帯域において前記包絡線を近似する近似直線に関連するパラメータをそれぞれ抽出するパラメータ抽出ステップと、
    変化量算出部が、前記候補帯域ごとに抽出された複数のパラメータのうち所定の関係にある2つの前記候補帯域でそれぞれ抽出した2つのパラメータの変化量を算出する変化量算出ステップと、
    決定部が、前記変化量に基づいて前記周波数帯域を決定する決定ステップと、
    を有することを特徴とする超音波観測装置の作動プログラム。
JP2016505643A 2014-12-22 2015-10-07 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム Active JP5932184B1 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2014259473 2014-12-22
JP2014259473 2014-12-22
PCT/JP2015/078547 WO2016103848A1 (ja) 2014-12-22 2015-10-07 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Publications (2)

Publication Number Publication Date
JP5932184B1 true JP5932184B1 (ja) 2016-06-08
JPWO2016103848A1 JPWO2016103848A1 (ja) 2017-04-27

Family

ID=56102945

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2016505643A Active JP5932184B1 (ja) 2014-12-22 2015-10-07 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Country Status (1)

Country Link
JP (1) JP5932184B1 (ja)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013166059A (ja) * 2010-11-11 2013-08-29 Olympus Medical Systems Corp 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013166059A (ja) * 2010-11-11 2013-08-29 Olympus Medical Systems Corp 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Also Published As

Publication number Publication date
JPWO2016103848A1 (ja) 2017-04-27

Similar Documents

Publication Publication Date Title
JP5974210B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016151951A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5948527B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2018142937A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2017110756A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6289772B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016103848A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6892320B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2016202567A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2020044044A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
US20170150943A1 (en) Ultrasound observation apparatus, method for operating ultrasound observation apparatus, and computer-readable recording medium
JP6157790B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6138402B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5932184B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5981072B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2017113145A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6253572B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5927367B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6010274B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5953457B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016103849A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JPWO2015198713A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2017217313A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Legal Events

Date Code Title Description
TRDD Decision of grant or rejection written
A975 Report on accelerated examination

Free format text: JAPANESE INTERMEDIATE CODE: A971005

Effective date: 20160407

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20160412

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160427

R151 Written notification of patent or utility model registration

Ref document number: 5932184

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250