JP4054454B2

JP4054454B2 - Image processing apparatus, image processing method, and computer-readable storage medium

Info

Publication number: JP4054454B2
Application number: JP27228298A
Authority: JP
Inventors: 弘之新畠; 真野北
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 1998-09-25
Filing date: 1998-09-25
Publication date: 2008-02-27
Anticipated expiration: 2018-09-25
Also published as: JP2000101840A

Description

【０００１】
【発明の属する技術分野】
本発明は、入力画像のダイナミックレンジの拡大、圧縮を行う技術に関し、特に、Ｘ線撮影された画像の微細構造を保存したまま特定領域のダイナミックレンジを拡大、縮小する場合に用いて好適なものである。
【０００２】
【従来の技術】
例えば、Ｘ線四肢部画像は、Ｘ線が透過しやすい軟部組織の画像、及びＸ線が非常に透過しにくい骨部の画像より構成されるため、画素値の存在するレンジが非常に広い。このため、軟部組織及ぴ骨部の両方を同時に観察することが可能なＸ線四肢部画像を得ることは困難であるとされてきた。
【０００３】
そこで、この問題を回避する方法として、従来より特許第２６６３１８９号記載のダイナミックレンジ圧縮方法がある。この方法は、補償後（処理後）の画素値Ｓ_D、オリジナル画素値（入力画素値）Ｓ_ORG、オリジナル画像をマスクサイズＭ×Ｍ画素で移動平均をとった時の平均画素値Ｓ_US及び図１５に示されるような特性を有する関数ｆ（Ｘ）を用いて、
【０００４】
Ｓ_D＝Ｓ_ORG＋ｆ（Ｓ_US） ───（１）
Ｓ_US＝ΣＳ_ORG／Ｍ² ───（２）
なる式（１）、（２）で表わされるものである。
【０００５】
ここで、関数ｆ（Ｓ_US）が有する特性について説明すると、まず、図１５に示す特性は、「０≦Ｓ_US≦ＢＡＳＥ」ではｆ（Ｓ_US）が「０」となり、「Ｓ_US＞ＢＡＳＥ」ではｆ（Ｓ_US）が切片を「しきい値ＢＡＳＥ」、傾き「ＳＬＯＰＥ」として単調減少するものである。
【０００６】
従って、オリジナルの画素値Ｓ_ORGを濃度相当量として、上記（１）式を実行した際には、画像の平均濃度の高いところで濃度を引き下げる、という画像に対する効果が得られる。
【０００７】
【発明が解決しようとする課題】
しかしながら、上述の従来方法では、予めＢＡＳＥとＳＬＯＰＥが決まってるため、どのような濃度分布の入力画像に対しても一律にダイナミックレンジを圧縮する範囲と度合いが決まっており、ダイナミックレンジを圧縮したくない注目領域まで圧縮してしまうと共に、逆に圧縮したい領域を圧縮できないという問題があった。
【０００８】
また、従来方法では入力画像の濃度分布の幅が狭い場合に、逆にダイナミックレンジを拡大するという思想がなく、濃度範囲を有効活用できないという問題があった。
【０００９】
また、例えば人体の四肢の様に、主に２種の組織（骨と軟部組織）から構成される物体の画像をフィルム出力する場合、注目領域（骨）はＦ／Ｓ系のフイルム出力と同様な出力が望まれ、さらに、Ｆ／Ｓ系のフィルムでは見えない高濃度領域（軟部組織の上部濃度）も見える様にフィルム出力することが望まれている。さらに、どのような入力画像に対しても同様の階調（例えば骨の濃度１．０、軟部組織の濃度の上限は３．０として出力する）のフィルム出力が望まれている。
【００１０】
しかしながら、従来の方法では、階調変換のための特徴量を安定して抽出するという思想がなく、一律に一定濃度領域のダイナミックレンジ領域を圧縮するため、入力画像の濃度分布の違いにより、出力画像の階調が変動してしまい、入力画像に依存しない安定した濃度の階調変換が得られないという問題があった。
【００１１】
本発明は、上記の問題に鑑みて為されたもので、放射線画像中の微細構造を示す画像の階調を保存しつつ、該放射線画像中の任意の画像領域の階調を圧縮或は拡大し、放射線画像における所望の観察領域の画像の階調を適正な状態に保ちながら該放射線画像の階調を所望の階調にすることを目的とする。
【００１２】
【課題を解決するための手段】
本発明の画像処理装置は、被写体にＸ線を照射することにより得られた原画像の階調変換処理を行う画像処理装置において、前記被写体を透過せずにＸ線が素通りした領域であるす抜け領域を代表する画素の値以上の画素値を有する画素の前記原画像上での座標から前記原画像上で一定幅の領域に属さない前記原画像内の領域を、被写体領域として抽出する手段と、前記被写体領域内の画素値の出現頻度を示すヒストグラムを得る手段と、前記ヒストグラムに基づき特徴量を算出する算出手段と、前記特徴量の値を有する画素が階調変換処理後に目的とする値となるように定めた階調変換曲線で前記原画像を階調変換処理する手段と、を有することを特徴とする。
本発明の画像処理方法は、被写体にＸ線を照射することにより得られた原画像の階調変換処理を行う画像処理方法において、前記被写体を透過せずにＸ線が素通りした領域であるす抜け領域を代表する画素の値以上の画素値を有する画素の前記原画像上での座標から前記原画像上で一定幅の領域に属さない前記原画像内の領域を、被写体領域として抽出する工程と、前記被写体領域内の画素値の出現頻度を示すヒストグラムを得る工程と、前記ヒストグラムに基づき特徴量を算出する工程と、前記特徴量の値を有する画素が階調変換処理後に目的とする値となるように定めた階調変換曲線で前記原画像を階調変換処理する工程と、を有することを特徴とする。
本発明のコンピュータ読み取り可能な記憶媒体は、被写体にＸ線を照射することにより得られた原画像の階調変換処理を行うためにコンピュータを、前記被写体を透過せずにＸ線が素通りした領域であるす抜け領域を代表する画素の値以上の画素値を有する画素の前記原画像上での座標から前記原画像上で一定幅の領域に属さない前記原画像内の領域を、被写体領域として抽出する手段と、前記被写体領域内の画素値の出現頻度を示すヒストグラムを得る手段と、前記ヒストグラムに基づき特徴量を算出する手段と、前記特徴量の値を有する画素が階調変換処理後に目的とする値となるように定めた階調変換曲線で前記原画像を階調変換処理する手段として機能させることを特徴とする。
【００１７】
【発明の実施の形態】
以下、本発明の実施の形態を図面と共に説明する。
図１は、本発明の第１の実施の形態による画像処理装置の構成を示すブロック図である。
図１において、１０１は入力画像からす抜け領域（Ｘ線が素通りした領域）を削除するす抜け削除手段、１０２はす抜け削除手段１０１です抜け領域を削除した画像の領域の濃度ヒストグラムを作成するヒストグラム作成手段、１０３はヒストグラム作成手段１０２で作成された濃度ヒストグラムの凹凸部を抽出する凹凸部抽出手段、１０４は凹凸抽出手段１０３で得られた凹凸部の情報を用いてダイナミックレンジ変更のための特徴量と階調変換のための特徴量を計算する解析手段、１０６は入力画像を平滑化する平滑化手段、１０７は解析手段１０４で得られた特徴量と平滑化手段１０６で得られた平滑画像に基づき濃度を変更する変更手段である。
【００１８】
図２は第１の実施の形態による画像処理装置の処理の流れを示すフローチャートである。
図３は横軸に濃度値、縦軸に濃度値の出現頻度をとった凹部を持つヒストグラムで、横軸のｍａｘが最大濃度値、Ｔａ２が左側凸部の最大値、濃度Ｔａ２の頻度をｈ１、右側凸部の最大値をＴａ３、濃度Ｔａ３の頻度をｈ２、濃度Ｔａ２以下の平均濃度値をＴａ１、濃度Ｔａ３以下の平均濃度値をＴａ４とする。Ｔａ１、Ｔａ４が後述する階調変換に用いる濃度値となる。
図４は横軸に濃度値、縦軸に濃度値の出現頻度をとった凹部の無いヒストグラムで、横軸のｍａｘが最大濃度値、Ｔａ６が凸部の最大値、濃度Ｔａ６の頻度をｈ３とするものである。
【００１９】
図５は高濃度側のダイナミックレンジを圧縮する場合に用いる関数を示し、図６は高濃度側のダイナミックレンジを拡大する場合に用いる関数を示す。両図において、ＢＡＳＥはダイナミックレンジ変更開始の濃度値である。
図７は入力画像濃度をＸ線フイルムの濃度に変換する階調変換曲線を示し、横軸が入力画像濃度値、縦軸が出力画像濃度値を表し、ｄ１が濃度１．０、ｄ２が濃度１．５、ｄ３が濃度３．０に対応する。上記階調変換曲線は、例えば画像データを実際のフィルム出力に模して出力する場合等に用いる。
【００２０】
次に動作について説明する。
す抜け削除手段１０１は、入力画像からす抜け領域とす抜け領域と一定輻で接する領域とを削除する。具体的には、Ｘ線の照射領域内のす抜け領域及びす抜け領域と一定間隔内で接する体領域を例えば０画素で置き換える次のような画像の変換を行う。
【００２１】
【数１】

【００２２】
ここで、ｆ（ｘ，ｙ）は画像データを示し、ｆ１（ｘ，ｙ）はす抜け領域及びす抜け領域と一定間隔内で接する体領域を削除した後の画像を示す。ｓｇｎ（ｘ，ｙ）は以下のように表される。Ｔｈ１は実験により定められる定数、ｄ１，ｄ２は体領域を削除する幅を決める定数である。
ｓｇｎ（ｘ，ｙ）＝０ｆ（ｘ，ｙ）≧Ｔｈ１のとき
ｓｇｎ（ｘ，ｙ）＝１その他 ───（４）
【００２３】
次にヒストグラム作成手段１０２は、す抜け削除手段１０１で０に置き換えられなかった濃度値のヒストグラムを作成する。そして凹凸部抽出手段１０３は、Ｆｈｉｓ（ｘ）をヒストグラム曲線、ｄを定数とし、
Ｆｈｉｓ（ｘ−ｄ）＞Ｆｈｉｓ（ｘ）かつＦｈｉｓ（ｘ＋ｄ）＞Ｆｈｉｓ（ｘ）
なる領域を凹部とし、その他の領域を凸部と判定する。また定数ｄは（５）式で示される濃度分布ｂ１で正規化され（６）式で表される。
ｂ１＝ｍａｘ−ｍｉｎ ───（５）
ｄ＝ｃ×ｂ１ ───（６）
【００２４】
次に、解析手段１０４の処理を図２の流れに従い説明する。
まず凹凸部抽出手段１０３で凹部が抽出されたかどうか判定し（ステップＳ２０１）、凹部が無い場合は、図４の凸部の最大値Ｔａ６を抽出する（ステップＳ２０２）。そして（７）式に従い階調変換のための濃度特徴量Ｔａ５を計算する。
【００２５】
【数２】

【００２６】
そして図７に示す階調変換曲線から所定の濃度値を抽出する。例えば、濃度値Ｔａ５が所望する縦軸の出力濃度値（例えば１．０）になる様に階調変換曲線を横軸に平行にシフトするか、入力濃度値を横軸に沿ってシフトする（例えば、これは実際のフイルム濃度上で注目領域Ｔａ５の濃度が示す領域の濃度値に対応する）。そして所望する出力濃度値の最高値（例えば３．０）に対応する入力画像ｄ３を抽出する（これは例えば一般日本人がフイルム上濃度で識別できる濃度限界である）。
【００２７】
次に、（８）、（９）、（１０）式に従いｒａｔｉｏ，ＳＬＯＰＥ１を計算する（ステップＳ２０４）。
ｒａｔｉｏ＝（ｄ３−Ｔａ６）／（ｍａｘ−Ｔａ６） ───（８）
もしｒａｔｉｏが１以下ならば
ＳＬＯＰＥ１＝１−ｒａｔｉｏ ───（９）
ｒａｔｉｏが１より大きければ
ＳＬＯＰＥ１＝ｒａｔｉｏ−１ ───（１０）
【００２８】
次に、ステップＳ２０１で凹部が見つかった場合は、図３の左凸部の最大値Ｔａ２とその頻度ｈ１を抽出する（ステップＳ２０５）。次に、右凸部の最大値Ｔａ３とその頻度ｈ２を抽出する（ステップＳ２０６）。そしてｈ１とｈ２の比率を計算し（ステップＳ２０７）、比率が定数ｃ１より大きい場合は、階調変換に用いる特徴量Ｔａ４を式（１１）に従い計算する（ステップＳ２０８）。
そして図７に示す階調変換曲線から所定の濃度値を抽出する。例えば、濃度値Ｔａ４が所望する縦軸の出力濃度値（例えば１．０）になる様に階調変換曲線を横軸に平行にシフトするか、入力濃度値を横軸に沿ってシフトする（例えば、これは実際のフィルム濃度上で注目領域Ｔａ４の濃度が示す領域の濃度値に対応する）。そして所望する出力濃度値の量高値（例えば３．０）に対応する入力画像ｄ３を抽出する（これは例えば一般日本人がフィルム上濃度で識別できる濃度限界である）。
【００２９】
【数３】

【００３０】
次に、（１２）、（１３）、（１４）式に従いｒａｔｉｏ，ＳＬＯＰＥ２を計算する（ステップＳ２０９）。
ｒａｔｉｏ＝（ｄ３−Ｔａ３）／（ｍａｘ−Ｔａ３） ───（１２）
もしｒａｔｉｏが１以下ならば
ＳＬＯＰＥ２＝１−ｒａｔｉｏ ───（１３）
ｒａｔｉｏが１より大きければ
ＳＬＯＰＥ２＝ｒａｔｉｏ−１ ──（１４）
【００３１】
ステップＳ２０７がＹｅｓの場合は、階調変換の特徴量に用いるＴａ１を（１５）式に従い計算する（ステップＳ２１０）。
そして図７に示す階調変換曲線から所定の濃度値を抽出する。例えば、濃度値Ｔａ１が所望する縦軸の出力濃度値（例えば１．０）になる様に階調変換曲線を横軸に平行にシフトするか、入力濃度値を横軸に沿ってシフトする（例えば、これは実際のフィルム濃度上で注目領域Ｔａ１の濃度が示す領域の濃度値に対応する）。そして所望する出力濃度値の最高値（例えば３．０）に対応する入力画像ｄ３を抽出する（これは例えば一般日本人がフィルム上濃度で識別できる濃度限界である）。
【００３２】
【数４】

【００３３】
次に、（１６）、（１７）、（１８）式に従いｒａｔｉｏ，ＳＬＯＰＥ３を計算する（ステップＳ２１１）。
ｒａｔｉｏ＝（ｄ３−Ｔａ２）／（ｍａｘ−Ｔａ２） ───（１６）
もしｒａｔｉｏが１以下ならば
ＳＬＯＰＥ３＝１−ｒａｔｉｏ ───（１７）
ｒａｔｉｏが１より大きければ
ＳＬＯＰＥ３＝ｒａｔｉｏ−１ ──（１８）
【００３４】
次に、平滑手段１０６はモルフォロジカルフィルタによる演算を（１９）、（２０）、（２１）、（２２）（２３）式に従い平滑画像ｆｕｓ（ｘ，ｙ）画像を計算する。ここでｆｏ（ｘ，ｙ）を入力画像とする。
【００３５】
ｆ１（ｘ，ｙ）＝ｍｉｎ｛ｆｏ（ｘ＋ｘ１，ｙ＋ｙ１）−Ｄ（ｘ１，ｙ１）｜ｘ１×ｘ１＋ｙ１×ｙ１≦ｒ１×ｒ１｝ ─────（１９）
ｆ２（ｘ，ｙ）＝ｍａｘ｛ｆ１（ｘ＋ｘ１，ｙ＋ｙ１）＋Ｄ（ｘ１，ｙ１）｜ｘ１×ｘ１＋ｙ１×ｙ１≦ｒ１×ｒ１｝ ─────（２０）
ｆ３（ｘ，ｙ）＝ｍａｘ｛ｆ２（ｘ＋ｘ１，ｙ＋ｙ１）＋Ｄ（ｘ１，ｙ１）｜ｘ１×ｘ１＋ｙ１×ｙ１≦ｒ１Ｘｒ１｝ ─────（２１）
ｆｕｓ（ｘ，ｙ）＝ｍｉｎ｛ｆ３（ｘ＋ｘ１，ｙ＋ｙ１）−Ｄ（ｘ１，ｙ１）｜ｘ１×ｘ１＋ｙ１×ｙ１≦ｒ１×ｒ１｝ ─────（２２）
【００３６】
Ｄ（ｘ，ｙ）を半球状フィルタ、ｒ１を任意の定数とし、入力画像に応じて選択される。
ｘ²＋ｙ²≦ｒ1²の場合、Ｄ（ｘ，ｙ）＝ｓｑｒｔ（（ｒ１−ｘ）²＋（ｒ１−ｙ）²）
その他の場合、Ｄ（ｘ，ｙ）＝−∞ ───（２３）
【００３７】
次に、ｒａｔｉｏが１以下の場合には、図５で示される単調減少関数ｆ１（）を用い、式（２４）に従い濃度変更を行う。ここでのＢＡＳＥにはＳＬＯＰＥ１を用いる場合はＴａ６、ＳＬＯＰＥ２を用いる場合はＴａ３、ＳＬＯＰＥ３を用いる場合はＴａ２を用いる。
ｆｅ（ｘ，ｙ）＝ｆｏ（ｘ，ｙ）＋ｆ３（ｆｕｓ（ｘ，ｙ））；ｆｕｓ（ｘ，ｙ）＞ＢＡＳＥ＝ｆｏ（ｘ，ｙ）；ｆｕｓ（ｘ，ｙ）≦ＢＡＳＥ ───（２４）
【００３８】
また、ｒａｔｉｏが１より大きい場合には、図６で示される単調増加関数ｆ２（）を用い、式（２５）に従い濃度変更を行う。
ｆｅ（ｘ，ｙ）＝ｆｏ（ｘ，ｙ）＋ｆ４（ｆｕｓ（ｘ，ｙ））；ｆｕｓ（ｘ・ｙ）＞ＢＡＳＥ＝ｆｏ（ｘ，ｙ）；ｆｕｓ（ｘ，ｙ）≦ＢＡＳＥ ───（２５）
【００３９】
そして、ＳＬＯＰＥ１を用いる場合はＴａ５、ＳＬＯＰＥ２を用いる場合はＴａ４、ＳＬＯＰＥ３を用いる場合はＴａ１を濃度値ｄ１に一致させ、図７で示される階調変換曲線で階調変換を行う（ステップＳ２０８）。従って、どのような入力データが入力されても階調変換のための濃度特徴量はフィルム上で１．０となり、ｍａｘは３．０となる。そのため、注目領域（ＢＡＳＥ以下の入力濃度値）は実際のＸ線フィルムの見方に近く、それ以外の濃度領域も上部の濃度限界が３．０となり、肉眼で識別できるように階調変換できる。
【００４０】
以上のように第１の実施の形態によれば、入力画像の最大濃度値が所定濃度値より小さい場合には、細部構造を不変のまま、所定濃度値までダイナミックレンジを拡大することができ、濃度領域を有効に活用することができる。また、入力画像の最大濃度値が所定濃度値より小さい場合には、細部構造を不変のまま、所定濃度値までダイナミックレンジを圧縮することができ、濃度領域を有効に活用することができる。
【００４１】
さらに、注目領域を実際のフィルム画像のように階調変換することができ、それ以外の濃度領域も所定濃度範囲に階調変換することが可能である。また、入力画像の濃度分布が主に２種の組織（例えば骨と軟部組織）より構成される場合に、濃度分布の形状により、ダイナミックレンジを変更するＢＡＳＥ位置を変更するので、低濃度側の分布（例えば骨）のダイナミックレンジの変更量は少なく、高濃度側（軟部組織）のダイナミックレンジを主に変更できる効果がある。
【００４２】
図８は本発明の第２の実施の形態による画像処理装置の構成を示すブロック図である。
図８において、３０１は入力画像の所定領域から濃度特徴量を抽出する第１の特徴量抽出手段、３０２は入力画像からす抜け領域を削除するす抜け削除手段、３０３はす抜け削除手段３０２です抜け領域を削除した領域から濃度特徴量を抽出する第２の特徴量抽出手段、３０４は第１の特徴量抽出手段３０１で抽出された濃度特徴量と第２の特徴量抽出手段３０３で抽出された濃度特徴量からダイナミックレンジ圧縮に用いるパラメータとしての濃度基準値と比率を計算し決定するパラメータ決定手段としての第３の特徴量抽出手段、３０５は入力画像を平滑化する平滑化手段、３０６は第３の特徴量抽出手段３０４で得られた上記パラメータと平滑化手段３０５で得られた平滑化画像に基づいて濃度を変更する変更手段である。
【００４３】
図９はす抜け領域を削除した肺側面の画像で、Ａが第１の特徴量抽出手段３０１が特徴量を抽出する所定領域、Ｂが肺領域と体領域の境界、Ｃがす抜け領域と接する部分を一定の割合で削除した体領域を示す。
【００４４】
図１０は高濃度側のダイナミックレンジを変更する場合に用いる関数を示し、ｄ１が第１の特徴量抽出手段３０１が抽出した特徴、実線が高濃度側のダイナミックレンジを圧縮する場合に用いる関数、点線が高濃度側のダイナミックレンジを拡大する場合に用いる関数を示す。
【００４５】
図１１は低濃度側のダイナミックレンジを変更する場合に用いる関数を示し、ｄ１が第１の特徴量抽出手段が抽出した特徴、実線が低濃度側のダイナミックレンジを圧縮する場合に用いる関数、点線が低濃度側のダイナミックレンジを拡
大する場合に用いる関数を示す。
【００４６】
図１２は入力画像濃度をｘ線フィルムの濃度に変換する階調変換曲線を示し、横軸が入力画像濃度値、縦軸が出力画像濃度値を表し、ｄ１が濃度１．５，ｄ２が濃度３．０，ｄ３が濃度Ｏ．２５に対応する。上記階調曲線は、例えば画像データを実際のフィルム出力に模して出力する場合等で用いる。
【００４７】
次に動作について説明する。
第１の特徴量抽出手段３０１は、図９のＡの領域の様に予め決められた所定領域から濃度特徴量を抽出し、この特徴量の値をｄ１とする。そして、このｄ１を階調変換及びダイナミックレンジ圧縮のための特徴量の計算に用いる。本実施の形態では、例えば所定領域Ａ内の濃度平均値を特徴量とするものである。また、例えば所定領域内の濃度中間値、最高値を用いてもよい。
【００４８】
次に、す抜け削除手段１０２はす抜け領域とす抜け領域と一定幅で接する領域を削除する。具体的には、Ｘ線の照射領域内のす抜け領域及びす抜け領域と一定間隔内で接する体領域を例えば０画素で置き換える前記（３）式による画像の変換を行う。
【００４９】
次に、第２の特徴量抽出手段３０３では、（２６）、（２７）式で示されるす抜け除去後の画像中で０でない領域の最上部座標Ｙ１と最下部Ｙｈを抽出し、（２８）式で示されるＹｌとＹｈの中間座標Ｙｍを計算する。
ＹＩ＝ｍｉｎ｛ｙ｜ｆ１（ｘ，ｙ）＞０｝ ───（２６）
Ｙｈ＝ｍａｘ｛ｙ｜ｆ１（ｘ，ｙ）＞０｝ ───（２７）
Ｙｍ＝（Ｙ１＋Ｙｈ）／２ ───（２８）
【００５０】
そして、（２９）式で示されるＹｌからＹｍ（頭部側半分）の間のｆ１（ｘ，ｙ）の最大値ｍａｘ１と（３０）式で示されるＹｍからＹｈ（腹側半分）の間の最大値ｍａｘ２を抽出する。ここでｍａｘ１は例えば腕前面上部の最大濃度値、ｍａｘ２は肺内の心臓下部の最大濃度値を示す。
【００５１】
ｍａｘ１＝ｍａｘ｛ｆ１（ｘ，ｙ）｜Ｙ１≦ｙ≦Ｙｍｌ───（２９）
ｍａｘ２＝ｍａｘ｛ｆ１（ｘ，ｙ）｜Ｙｍ≦ｙ≦Ｙｈ｝───（３０）
さらに、（３１）式で示されるｆ１（ｘ，ｙ）の最小値をＭｉｎとする。
ｍｉｎ＝ｍｉｎ｛ｆ１（ｘ，ｙ）｜ｆ１（ｘ，ｙ）≠０｝──（３１）
【００５２】
そして図１２に示す階調変換曲線から所定の濃度値を抽出する。例えば、濃度値ｄ１が所望の縦軸の出力濃度値（例えば１．５）になる様に階調変換曲線を横軸に平行にシフトするか、入力濃度値を横軸に沿ってシフトする。そして所望の出力濃度値の最高値（例えば３．０）に対応する入力画像値ｄ２を抽出する（これは一般日本人がフイルム上濃度で識別できる濃度限界である）。そして、所望の出力濃度値の最低値（例えば０．２５）に対応する入力画像値ｄ３を抽出する。
【００５３】
次に、（３２）、（３３）、（３４）式に従い高濃度のダイナミックレンジの圧縮、拡大に用いるｒａｔｉｏ、ＳＬＯＰＥを計算する。ここでのｍａｘはｍａｘ１またはｍａｘ２のいずれかを示す。
ｒａｔｉｏ＝（ｄ２−ｄ１）／（ｍａｘ−ｄ１） ───（３２）
もしｒａｔｉｏが１以下ならば
ＳＬＯＰＥ＝１−ｒａｔｉｏ ───（３３）
ｒａｔｉｏが１より大きければ
ＳＬＯＰＥ＝ｒａｔｉｏ−１ ───（３４）
【００５４】
次に、（３５）、（３６）、（３７）式に従い低濃度のダイナミックレンジの圧縮、拡大に用いるｒａｔｉｏ１，ＳＬＯＰＥ１を計算する。
ｒａｔｉｏ１＝（ｄ１−ｄ３）／（ｄ１−ｍｉｎ） ───（３５）
【００５５】
もしｒａｔｉｏ１が１以下ならば
ＳＬＯＰＥｌ＝１−ｒａｔｉｏ ───（３６）
ｒａｔｉｏが１より大きければ
ＳＬＯＰＥｌ＝ｒａｔｉｏ−１ ───（３７）
【００５６】
そして、平滑化手段３０５はモルフォロジカルフィルタによる演算を前記（１９）〜（２２）式と同様の演算に従い、ｆ１（ｘ，ｙ）〜ｆ３（ｘ，ｙ）を順次求めながら平滑画像ｆｕｓ（ｘ，ｙ）画像を計算する。ここでｆｏ（ｘ，ｙ）を入力画像とする。また、Ｄ（ｘ，ｙ）のフィルタも前記（２３）式と同様である。
【００５７】
そして、変更手段１０６は高濃度側と低濃度側のダイナミックレンジ圧縮、拡大を行う。
ｒａｔｉｏが１以下の場合には、図１０の実線で示される単調減少関数、ｒａｔｉｏが１以上の場合には、単調増加関数となるｆ１（）を用い、式（３８）に従い濃度変更を行う。ここでのＢＡＳＥにはｄ１を用いる。変換後の画像をｆｅ（ｘ，ｙ）とする。
【００５８】
ｆｅ（ｘ，ｙ）＝ｆｏ（ｘ，ｙ）＋ｆ１（ｆｕｓ（ｘ，ｙ））；ｆｕｓ（ｘ，ｙ）＞ＢＡＳＥ＝ｆｏ（ｘ，ｙ）；ｆｕｓ（ｘ，ｙ）≦ＢＡＳＥ ───（３８）
【００５９】
さらに、低濃度側のダイナミックレンジ圧縮、拡大を行う。ｒａｔｉｏ１が１以下の場合に、図１１の実線で示される単調減少関数、ｒａｔｉｏが１以上の場合には単調増加関数となるｆ３（）を用い、式（３９）に従い濃度変更を行う。ここでのＢＡＳＥにはｄ１を用いる。変換後の画像をｆｅ（ｘ，ｙ）とする。
【００６０】
ｆｅ（ｘ，ｙ）＝ｆｏ（ｘ，ｙ）＋ｆ１（ｆｕｓ（ｘ，ｙ））；ｆｕｓ（ｘ，ｙ）≦ＢＡＳＥ＝ｆｏ（ｘ，ｙ）；ｆｕｓ（ｘ，ｙ）＞ＢＡＳＥ ───（３９）
【００６１】
そして、入力画像のｄ１を濃度値１．５に変換される位置に一致させ、図１２の階調変換曲線で階調変換を行う。従って、どのような入力データが入力されても階調変換のための濃度特徴量はフィルム上で１．５となり、ｍａｘは３．０、ｍｉｎは０．２５となる。そのため、注目領域（ダイナミックレンジの変換をしない領域）は実際のＸ線フィルムの見方に近く、それ以外の濃度領域も上部の濃度限界が３．０、下部の濃度値が０．２５となり、肉眼で識別できるように階調変換できる。
【００６２】
以上のように第２の実施の形態によれば、入力画像の最大濃度値が所定濃度値より小さい場合には、細部構造を不変のまま、所定濃度値までダイナミックレンジを拡大することができ、濃度領域を有効に活用することができる。さらに、入力画像の最大濃度値が所定濃度値より大きい場合には、細部構造を不変のまま、所定濃度値までダイナミックレンジを圧縮することができ、濃度領域を有効に活用することができる。さらに注目領域は実際のフィルム画像のように階調変換することができ、それ以外の濃度領域も所定濃度範囲に階調変換することが可能である。
【００６３】
次に、第３の実施の形態を説明する。
図１３は高濃度側のダイナミックレンジを圧縮する場合に用いる関数を示し、ｄ１が第１の特徴量抽出手段３０１が抽出した特徴量、実線が高濃度側のダイナミックレンジを圧縮する場合に用いる関数、点線が高濃度側のダイナミックレンジを拡大する場合に用いる関数を示す。ＳＬＯＰＥの開始点をｍａｘとｄ１の間の任意の位置ｄ４とする。
【００６４】
図１４は低濃度側のダイナミックレンジを圧縮する場合に用いる関数を示し、ｄ１が第１の特徴量抽出手段が抽出した特徴量、実線が低濃度側のダイナミックレンジを圧縮する場合に用いる関数、点線が低濃度側のダイナミックレンジを拡大する場合に用いる関数を示す。ＳＬＯＰＥの開始点をｍｉｎとｄ１の間の任意の位置ｄ５とする。
【００６５】
次に動作について説明する。
（４０）、（４１）、（４２）式に従い高濃度のダイナミックレンジの圧縮、拡大に用いるｒａｔｉｏ，ＳＬＯＰＥを計算する。ここでのｍａｘはｍａｘ１またはｍａｘ２のいずれかを示す。
【００６６】
ｒａｔｉｏ＝（ｄ２−ｄ４）／（ｍａｘ−ｄ４） ───（４０）
もしｒａｔｉｏが１以下ならば
ＳＬＯＰＥ＝１−ｒａｔｉｏ ───（４１）
ｒａｔｉｏが１より大きければ
ＳＬＯＰＥ＝ｒａｔｉｏ−１ ───（４２）
【００６７】
次に、（４３）、（４４）、（４５）式に従い低濃度のダイナミックレンジの圧縮、拡大に用いるｒａｔｉｏ１，ＳＬＯＰＥ１を計算する。
ｒａｔｉｏ１＝（ｄ５−ｄ３）／（ｄ５−ｍｉｎ） ───（４３）
【００６８】
もしｒａｔｉｏ１が１以下ならば
ＳＬＯＰＥ１＝１−ｒａｔｉｏ１ ───（４４）
ｒａｔｉｏ１が１より大きければ
ＳＬＯＰＥ１＝ｒａｔｉｏ１−１ ───（４５）
以降は第２の実施の形態と同様の処理を行う。
【００６９】
以上のように第３の実施の形態によればＳＬＯＰＥの開始点を任意に移動できるため、ダイナミックレンジを変更する範囲を調整できる効果がある。
【００７０】
次に、本発明による記憶媒体について説明する。
図１、図８の各機能ブロックによるシステムは，ハード的に構成してもよく、また、ＣＰＵやメモリ等からなるコンピュータシステムに構成してもよい。コンピュータシステムに構成する場合、上記メモリは本発明による記憶媒体を構成する。この記憶媒体には、前述した動作を制御するための図２のフローチャートを含む処理手順を実行するためのプログラムが記憶される。
【００７１】
この記憶媒体としては、ＲＯＭ、ＲＡＭ等の半導体メモリ、光ディスク、光磁気ディスク、磁気媒体等を用いてよく、これらをＣＤ−ＲＯＭ、ＦＤ、磁気カード、磁気テープ、不揮発性メモリカード等に構成して用いてよい。
【００７２】
従って、この記憶媒体を図１、図８に示したシステム以外の他のシステムあるいは装置で用い、そのシステムあるいはコンピュータがこの記憶媒体に格納されたプログラムコードを読み出し、実行することによっても、前述した各実施の形態と同等の機能を実現できると共に、同等の効果を得ることができ、本発明の目的を達成することができる。
【００７３】
また、コンピュータ上で稼働しているＯＳ等が処理の一部又は全部を行う場合、あるいは記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された拡張機能ボードやコンピュータに接続された拡張機能ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づいて、上記拡張機能ボードや拡張機能ユニットに備わるＣＰＵ等が処理の一部又は全部を行う場合にも、各実施の形態と同等の機能を実現できると共に、同等の効果を得ることができ、本発明の目的を達成することができる。
【００７４】
【発明の効果】
本発明によれば、原画像の階調変換処理に必要とする特徴量を算出するのに必要な被写体領域を算出できる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態による画像処理装置を示すブロック図である。
【図２】第１の実施の形態による処理の流れを示すフローチャートである。
【図３】濃度ヒストグラムを示す特性図である。
【図４】濃度ヒストグラムを示す特性図である。
【図５】関数ｆ１（）の形態を示す特性図である。
【図６】関数ｆ２（）の形態を示す特性図である。
【図７】階調変換曲線を示す特性図である。
【図８】本発明の第２の実施の形態による画像処理装置を示すブロック図である。
【図９】す抜け領域を削除した肺側面画像を示す構成図である。
【図１０】関数ｆ１の形態を示す特性図である。
【図１１】関数ｆ３の形態を示す特性図である。
【図１２】階調変換曲線を示す特性図である。
【図１３】関数ｆ２の形態を示す特性図である。
【図１４】関数ｆ４の形態を示す特性図である。
【図１５】従来の関数の形態を示す特性図である。
【符号の説明】
１０１す抜け削除手段
１０２ヒストグラム作成手段
１０３凹凸部抽出手段
１０４解析手段
１０６平滑化手段
１０７変更手段
３０１第１の特徴量抽出手段
３０２す抜け削除手段
３０３第２の特徴量抽出手段
３０４第３の特徴量抽出手段
３０５平滑化手段
３０６変更手段[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a technique for expanding and compressing the dynamic range of an input image, and particularly suitable for expanding and reducing the dynamic range of a specific area while preserving the fine structure of an X-ray image. It is.
[0002]
[Prior art]
For example, an X-ray limb image is composed of an image of a soft tissue that easily transmits X-rays and an image of a bone that is very difficult to transmit X-rays, and thus has a very wide range of pixel values. For this reason, it has been considered difficult to obtain an X-ray limb image capable of simultaneously observing both soft tissue and bone.
[0003]
Therefore, as a method of avoiding this problem, there is a conventional dynamic range compression method described in Japanese Patent No. 2663189. This method uses the pixel value S after compensation (after processing). _D Original pixel value (input pixel value) S _ORG The average pixel value S when the moving average of the original image is taken with a mask size of M × M pixels _US And a function f (X) having characteristics as shown in FIG.
[0004]
S _D = S _ORG + F (S _US ) (1)
S _US = ΣS _ORG / M ² ─── (2)
It is represented by the following formulas (1) and (2).
[0005]
Here, the function f (S _US ) Will be described. First, the characteristic shown in FIG. _US ≦ BASE ”, f (S _US ) Becomes "0" and "S _US > BASE ”is f (S _US ) Monotonously decreases with the intercept as “threshold BASE” and slope “SLOPE”.
[0006]
Therefore, the original pixel value S _ORG When the above equation (1) is executed with the density equivalent amount, the effect of reducing the density at a high average density of the image is obtained.
[0007]
[Problems to be solved by the invention]
However, in the above-described conventional method, since BASE and SLOPE are determined in advance, the range and degree of compression of the dynamic range are determined uniformly for an input image having any density distribution, and it is desired to compress the dynamic range. There is a problem in that the region to be compressed is not compressed and the region to be compressed cannot be compressed.
[0008]
Further, the conventional method has a problem that when the density distribution width of the input image is narrow, there is no idea of expanding the dynamic range, and the density range cannot be effectively used.
[0009]
In addition, for example, when an image of an object composed mainly of two types of tissues (bone and soft tissue) such as the extremities of a human body is output as a film, the attention area (bone) is the same as the F / S system film output. In addition, it is desired to output the film so that a high density region (upper density of the soft tissue) that cannot be seen with the F / S film is also visible. Further, it is desired to output a film having the same gradation (for example, the bone density is 1.0 and the upper limit of the soft tissue density is 3.0) for any input image.
[0010]
However, in the conventional method, there is no idea of stably extracting feature quantities for tone conversion, and the dynamic range area of a constant density area is uniformly compressed. The gradation of the image fluctuates, and there is a problem that gradation conversion with a stable density independent of the input image cannot be obtained.
[0011]
The present invention has been made in view of the above problems, and compresses or enlarges the gradation of an arbitrary image area in the radiation image while preserving the gradation of the image showing the fine structure in the radiation image. Then, it is an object to make the gradation of the radiation image a desired gradation while maintaining the gradation of the image of the desired observation region in the radiation image in an appropriate state.
[0012]
[Means for Solving the Problems]
The image processing apparatus of the present invention is an area in which X-rays pass through without passing through the subject in an image processing device that performs gradation conversion processing of an original image obtained by irradiating the subject with X-rays. Pixel values that represent missing areas A region in the original image that does not belong to a region of a certain width on the original image from the coordinates on the original image of pixels having the above pixel values is defined as a subject region. Means for extracting; Within the subject area Means for obtaining a histogram indicating the appearance frequency of pixel values; and calculating means for calculating a feature amount based on the histogram; The original image is represented by a gradation conversion curve determined so that a pixel having the feature value has a target value after gradation conversion processing. And means for performing gradation conversion processing.
The image processing method of the present invention is an area in which X-rays pass through the subject without passing through the subject in an image processing method for performing gradation conversion processing of an original image obtained by irradiating the subject with X-rays. Pixel values that represent missing areas A region in the original image that does not belong to a region of a certain width on the original image from the coordinates on the original image of pixels having the above pixel values is defined as a subject region. Extracting, and Within the subject area Obtaining a histogram indicating the appearance frequency of pixel values; calculating a feature amount based on the histogram; The original image is represented by a gradation conversion curve determined so that a pixel having the feature value has a target value after gradation conversion processing. And a step of performing gradation conversion processing.
The computer-readable storage medium of the present invention is a region through which X-rays pass through a computer without passing through the subject to perform gradation conversion processing of an original image obtained by irradiating the subject with X-rays. The value of the pixel representing the missing area A region in the original image that does not belong to a region of a certain width on the original image from the coordinates on the original image of pixels having the above pixel values is defined as a subject region. Means for extracting; Within the subject area Means for obtaining a histogram indicating the appearance frequency of pixel values; means for calculating a feature amount based on the histogram; The original image is represented by a gradation conversion curve determined so that a pixel having the feature value has a target value after gradation conversion processing. It is characterized by functioning as means for gradation conversion processing.
[0017]
DETAILED DESCRIPTION OF THE INVENTION
Embodiments of the present invention will be described below with reference to the drawings.
FIG. 1 is a block diagram showing the configuration of the image processing apparatus according to the first embodiment of the present invention.
In FIG. 1, 101 is a void deletion means for deleting a void area (area through which X-rays pass) from an input image, and 102 is a void deletion means 101. A density histogram of an image area from which the void area is deleted is created. Histogram creation means 103 is a concavo-convex part extraction means for extracting the concavo-convex part of the density histogram created by the histogram creation means 102, and 104 is for changing the dynamic range using the information of the concavo-convex part obtained by the concavo-convex extraction means 103 Analysis means for calculating the feature quantity and the feature quantity for gradation conversion, 106 a smoothing means for smoothing the input image, 107 a feature quantity obtained by the analysis means 104 and smoothing obtained by the smoothing means 106 This is a changing means for changing the density based on the image.
[0018]
FIG. 2 is a flowchart showing a processing flow of the image processing apparatus according to the first embodiment.
FIG. 3 is a histogram having a concave portion with the horizontal axis representing density values and the vertical axis representing density value appearance frequency, where the horizontal axis max is the maximum density value, Ta2 is the maximum value of the left convex portion, and the frequency of the density Ta2 is h1. The maximum value of the right convex portion is Ta3, the frequency of concentration Ta3 is h2, the average concentration value of concentration Ta2 or less is Ta1, and the average concentration value of concentration Ta3 or less is Ta4. Ta1 and Ta4 are density values used for gradation conversion described later.
FIG. 4 is a histogram without a concave portion in which the horizontal axis indicates density values, and the vertical axis indicates the appearance frequency of density values. The maximum on the horizontal axis is the maximum density value, Ta6 is the maximum value of the convex portions, and the frequency of density Ta6 is h3. To do.
[0019]
FIG. 5 shows a function used when compressing the dynamic range on the high density side, and FIG. 6 shows a function used when expanding the dynamic range on the high density side. In both figures, BASE is a density value at the start of dynamic range change.
FIG. 7 shows a gradation conversion curve for converting the input image density into the X-ray film density. The horizontal axis represents the input image density value, the vertical axis represents the output image density value, d1 is density 1.0, and d2 is density. 1.5 and d3 correspond to a density of 3.0. The gradation conversion curve is used, for example, when image data is output in the form of actual film output.
[0020]
Next, the operation will be described.
The void deletion means 101 deletes a void region, a void region, and a region that is in contact with a certain distance from the input image. More specifically, the following image conversion is performed in which, for example, a pixel region in an X-ray irradiation region and a body region in contact with the void region within a certain interval are replaced with, for example, 0 pixels.
[0021]
[Expression 1]

[0022]
Here, f (x, y) indicates image data, and f1 (x, y) indicates an image after deleting a void region and a body region that contacts the void region within a certain interval. sgn (x, y) is expressed as follows. Th1 is a constant determined by experiments, and d1 and d2 are constants that determine the width for deleting the body region.
When sgn (x, y) = 0 f (x, y) ≧ Th1
sgn (x, y) = 1 Other ─── (4)
[0023]
Next, the histogram creation unit 102 creates a histogram of density values that have not been replaced with 0 by the skip deletion unit 101. And the uneven | corrugated | grooved part extraction means 103 makes Fhis (x) a histogram curve, d is a constant,
Fhis (x−d)> Fhis (x) and Fhis (x + d)> Fhis (x)
A region to be formed is defined as a concave portion, and the other regions are determined as convex portions. The constant d is normalized by the concentration distribution b1 shown by the equation (5) and expressed by the equation (6).
b1 = max-min (5)
d = c × b1 (6)
[0024]
Next, the processing of the analysis unit 104 will be described according to the flow of FIG.
First, it is determined whether or not the concave portion is extracted by the concave and convex portion extraction means 103 (step S201). If there is no concave portion, the maximum value Ta6 of the convex portion in FIG. 4 is extracted (step S202). Then, the density feature amount Ta5 for gradation conversion is calculated according to the equation (7).
[0025]
[Expression 2]

[0026]
Then, a predetermined density value is extracted from the gradation conversion curve shown in FIG. For example, the gradation conversion curve is shifted parallel to the horizontal axis or the input density value is shifted along the horizontal axis so that the density value Ta5 becomes the desired output density value on the vertical axis (for example, 1.0) ( For example, this corresponds to the density value of the area indicated by the density of the attention area Ta5 on the actual film density). Then, an input image d3 corresponding to the maximum desired output density value (for example, 3.0) is extracted (this is a density limit that can be identified by, for example, the density on the film for general Japanese).
[0027]
Next, ratio and SLOPE1 are calculated according to the equations (8), (9), and (10) (step S204).
ratio = (d3-Ta6) / (max-Ta6) (8)
If ratio is less than 1,
SLOPE1 = 1-ratio ─── (9)
If the ratio is greater than 1,
SLOPE1 = ratio-1 (10)
[0028]
Next, if a concave portion is found in step S201, the maximum value Ta2 of the left convex portion in FIG. 3 and its frequency h1 are extracted (step S205). Next, the maximum value Ta3 and the frequency h2 of the right convex part are extracted (step S206). Then, the ratio between h1 and h2 is calculated (step S207). If the ratio is larger than the constant c1, the feature amount Ta4 used for gradation conversion is calculated according to the equation (11) (step S208).
Then, a predetermined density value is extracted from the gradation conversion curve shown in FIG. For example, the gradation conversion curve is shifted parallel to the horizontal axis or the input density value is shifted along the horizontal axis so that the density value Ta4 becomes the desired output density value on the vertical axis (for example, 1.0) ( For example, this corresponds to the density value of the area indicated by the density of the attention area Ta4 on the actual film density). Then, an input image d3 corresponding to a desired high value (for example, 3.0) of the output density value is extracted (this is a density limit that can be identified by, for example, a density on a general Japanese film).
[0029]
[Equation 3]

[0030]
Next, ratio and SLOPE2 are calculated according to the equations (12), (13), and (14) (step S209).
ratio = (d3-Ta3) / (max-Ta3) (12)
If ratio is less than 1,
SLOPE2 = 1-ratio (13)
If the ratio is greater than 1,
SLOPE2 = ratio-1 (14)
[0031]
When step S207 is Yes, Ta1 used for the feature value of gradation conversion is calculated according to the equation (15) (step S210).
Then, a predetermined density value is extracted from the gradation conversion curve shown in FIG. For example, the gradation conversion curve is shifted in parallel with the horizontal axis or the input density value is shifted along the horizontal axis so that the density value Ta1 becomes the desired output density value on the vertical axis (for example, 1.0) ( For example, this corresponds to the density value of the area indicated by the density of the attention area Ta1 on the actual film density). Then, an input image d3 corresponding to the desired maximum value of the output density value (for example, 3.0) is extracted (this is a density limit that can be identified by, for example, the density on the film for general Japanese).
[0032]
[Expression 4]

[0033]
Next, ratio and SLOPE3 are calculated according to the equations (16), (17), and (18) (step S211).
ratio = (d3-Ta2) / (max-Ta2) (16)
If ratio is less than 1,
SLOPE3 = 1-ratio ─── (17)
If the ratio is greater than 1,
SLOPE3 = ratio-1 (18)
[0034]
Next, the smoothing means 106 calculates a smoothed image fus (x, y) image according to the equations (19), (20), (21), (22), and (23) for the calculation by the morphological filter. Here, let fo (x, y) be an input image.
[0035]
f1 (x, y) = min {fo (x + x1, y + y1) −D (x1, y1) | x1 × x1 + y1 × y1 ≦ r1 × r1} (19)
f2 (x, y) = max {f1 (x + x1, y + y1) + D (x1, y1) | x1 × x1 + y1 × y1 ≦ r1 × r1} (20)
f3 (x, y) = max {f2 (x + x1, y + y1) + D (x1, y1) | x1 × x1 + y1 × y1 ≦ r1Xr1} (21)
fus (x, y) = min {f3 (x + x1, y + y1) −D (x1, y1) | x1 × x1 + y1 × y1 ≦ r1 × r1} (22)
[0036]
D (x, y) is a hemispherical filter, r1 is an arbitrary constant, and is selected according to the input image.
x ² + Y ² ≤r1 ² Then D (x, y) = sqrt ((r1-x) ² + (R1-y) ² )
In other cases, D (x, y) = − ∞ (23)
[0037]
Next, when the ratio is 1 or less, the density is changed according to the equation (24) using the monotone decreasing function f1 () shown in FIG. In this BASE, Ta6 is used when SLOPE1 is used, Ta3 is used when SLOPE2 is used, and Ta2 is used when SLOPE3 is used.
fe (x, y) = fo (x, y) + f3 (fus (x, y)); fus (x, y)> BASE = fo (x, y); fus (x, y) ≦ BASE (24)
[0038]
When the ratio is larger than 1, the density is changed according to the equation (25) using the monotonically increasing function f2 () shown in FIG.
fe (x, y) = fo (x, y) + f4 (fus (x, y)); fus (x · y)> BASE = fo (x, y); fus (x, y) ≦ BASE (25)
[0039]
Then, Ta5 is used when SLOPE1 is used, Ta4 is used when SLOPE2 is used, and Ta1 is used when SLOPE3 is used, and gradation conversion is performed using the gradation conversion curve shown in FIG. 7 (step S208). Therefore, regardless of what input data is input, the density feature amount for gradation conversion is 1.0 on the film, and max is 3.0. Therefore, the attention area (input density value below BASE) is close to the actual view of the X-ray film, and other density areas have an upper density limit of 3.0, so that gradation conversion can be performed so that they can be identified with the naked eye.
[0040]
As described above, according to the first embodiment, when the maximum density value of the input image is smaller than the predetermined density value, the dynamic range can be expanded to the predetermined density value without changing the detailed structure. The concentration region can be used effectively. In addition, when the maximum density value of the input image is smaller than the predetermined density value, the dynamic range can be compressed to the predetermined density value without changing the detailed structure, and the density area can be effectively used.
[0041]
Furthermore, it is possible to perform gradation conversion of the attention area like an actual film image, and it is also possible to perform gradation conversion of other density areas into a predetermined density range. Also, when the density distribution of the input image is mainly composed of two types of tissues (for example, bone and soft tissue), the BASE position for changing the dynamic range is changed depending on the shape of the density distribution. The amount of change in the dynamic range of the distribution (eg, bone) is small, and there is an effect that the dynamic range on the high concentration side (soft tissue) can be mainly changed.
[0042]
FIG. 8 is a block diagram showing a configuration of an image processing apparatus according to the second embodiment of the present invention.
In FIG. 8, reference numeral 301 denotes a first feature amount extraction unit that extracts a density feature amount from a predetermined region of the input image, 302 denotes a omission deletion unit that deletes a omission region from the input image, and 303 denotes an omission deletion unit 302. A second feature amount extraction unit 304 extracts a density feature amount from the region from which the missing region is deleted, and 304 is extracted by the density feature amount extracted by the first feature amount extraction unit 301 and the second feature amount extraction unit 303. A third feature quantity extraction means as a parameter determination means for calculating and determining a density reference value and a ratio as parameters used for dynamic range compression from the obtained density feature quantity; 305, a smoothing means for smoothing the input image; and 306, This is a changing means for changing the density based on the parameter obtained by the third feature quantity extracting means 304 and the smoothed image obtained by the smoothing means 305.
[0043]
FIG. 9 is an image of the lung side surface from which the omission region has been deleted. A body region in which the touching part is deleted at a certain rate is shown.
[0044]
FIG. 10 shows a function used when changing the dynamic range on the high density side, d1 is a feature extracted by the first feature quantity extraction unit 301, and a solid line is a function used when the dynamic range on the high density side is compressed. A dotted line indicates a function used when the dynamic range on the high density side is expanded.
[0045]
FIG. 11 shows a function used when changing the dynamic range on the low density side, d1 is a feature extracted by the first feature quantity extraction means, and a solid line is a function used when compressing the dynamic range on the low density side, a dotted line Expands the dynamic range on the low concentration side
Indicates the function to be used when increasing.
[0046]
FIG. 12 shows a gradation conversion curve for converting input image density to x-ray film density, the horizontal axis represents the input image density value, the vertical axis represents the output image density value, d1 is density 1.5, and d2 is density. 3.0, d3 is the concentration O.D. 25. The gradation curve is used, for example, when outputting image data imitating actual film output.
[0047]
Next, the operation will be described.
The first feature quantity extraction unit 301 extracts a density feature quantity from a predetermined area such as the area A in FIG. 9, and sets the value of this feature quantity to d1. This d1 is used to calculate the feature amount for gradation conversion and dynamic range compression. In the present embodiment, for example, an average density value in the predetermined area A is used as a feature amount. Further, for example, an intermediate density value or a maximum value in a predetermined area may be used.
[0048]
Next, the void deletion means 102 deletes the void region and the region that contacts the void region with a certain width. Specifically, the image conversion is performed according to the above equation (3) in which the void region in the X-ray irradiation region and the body region in contact with the void region within a certain interval are replaced with, for example, 0 pixels.
[0049]
Next, the second feature quantity extraction unit 303 extracts the uppermost coordinates Y1 and the lowermost part Yh of the non-zero region in the image after the omission removal shown by the equations (26) and (27), and (28 ) The intermediate coordinates Ym between Yl and Yh shown by the equation are calculated.
YI = min {y | f1 (x, y)> 0} (26)
Yh = max {y | f1 (x, y)> 0} (27)
Ym = (Y1 + Yh) / 2 (28)
[0050]
Then, the maximum value max1 of f1 (x, y) between Yl and Ym (head side half) represented by the equation (29) and between Ym and Yh (ventral half) represented by the equation (30) The maximum value max2 is extracted. Here, for example, max1 indicates the maximum density value at the upper front of the arm, and max2 indicates the maximum density value at the lower heart in the lung.
[0051]
max1 = max {f1 (x, y) | Y1 ≦ y ≦ Yml (29)
max2 = max {f1 (x, y) | Ym ≦ y ≦ Yh} (30)
Furthermore, let the minimum value of f1 (x, y) shown by the equation (31) be Min.
min = min {f1 (x, y) | f1 (x, y) ≠ 0} (31)
[0052]
Then, a predetermined density value is extracted from the gradation conversion curve shown in FIG. For example, the gradation conversion curve is shifted parallel to the horizontal axis or the input density value is shifted along the horizontal axis so that the density value d1 becomes the desired output density value (for example, 1.5) on the vertical axis. Then, the input image value d2 corresponding to the maximum value (for example, 3.0) of the desired output density value is extracted (this is a density limit that can be identified by the density on the film for general Japanese). Then, an input image value d3 corresponding to the lowest desired output density value (for example, 0.25) is extracted.
[0053]
Next, the ratio and SLOPE used for compression and expansion of the high density dynamic range are calculated according to the equations (32), (33), and (34). Here, max indicates either max1 or max2.
ratio = (d2-d1) / (max-d1) (32)
If ratio is less than 1,
SLOPE = 1-ratio ─── (33)
If the ratio is greater than 1,
SLOPE = ratio-1 (34)
[0054]
Next, ratio1 and SLOPE1 used for compression and expansion of the low density dynamic range are calculated according to equations (35), (36), and (37).
ratio1 = (d1-d3) / (d1-min) (35)
[0055]
If ratio1 is 1 or less
SLOPEl = 1-ratio ─── (36)
If the ratio is greater than 1,
SLOPEl = ratio-1 (37)
[0056]
Then, the smoothing means 305 performs the calculation by the morphological filter in accordance with the same calculation as the expressions (19) to (22), and sequentially obtains f1 (x, y) to f3 (x, y) while smoothing the image fus (x , Y) Calculate the image. Here, let fo (x, y) be an input image. Further, the filter of D (x, y) is the same as the equation (23).
[0057]
The changing means 106 performs dynamic range compression and expansion on the high density side and the low density side.
When the ratio is 1 or less, the density is changed according to the equation (38) using the monotonically decreasing function indicated by the solid line in FIG. 10, and when the ratio is 1 or more, f1 () that is the monotonically increasing function is used. Here, d1 is used for BASE. Let the image after conversion be fe (x, y).
[0058]
fe (x, y) = fo (x, y) + f1 (fus (x, y)); fus (x, y)> BASE = fo (x, y); fus (x, y) ≦ BASE (38)
[0059]
Furthermore, the dynamic range is compressed and expanded on the low density side. The density is changed according to the equation (39) using the monotonic decreasing function indicated by the solid line in FIG. 11 when the ratio 1 is 1 or less, and f3 () which is the monotonically increasing function when the ratio is 1 or more. Here, d1 is used for BASE. Let the image after conversion be fe (x, y).
[0060]
fe (x, y) = fo (x, y) + f1 (fus (x, y)); fus (x, y) ≦ BASE = fo (x, y); fus (x, y)> BASE (39)
[0061]
Then, d1 of the input image is made to coincide with the position where the density value is converted to 1.5, and gradation conversion is performed using the gradation conversion curve of FIG. Therefore, regardless of the input data, the density feature amount for gradation conversion is 1.5 on the film, max is 3.0, and min is 0.25. Therefore, the attention area (area where the dynamic range is not converted) is close to the actual view of the X-ray film, and other density areas have an upper density limit of 3.0 and a lower density value of 0.25. The tone can be converted so that it can be identified.
[0062]
As described above, according to the second embodiment, when the maximum density value of the input image is smaller than the predetermined density value, the dynamic range can be expanded to the predetermined density value without changing the detailed structure. The concentration region can be used effectively. Further, when the maximum density value of the input image is larger than the predetermined density value, the dynamic range can be compressed to the predetermined density value without changing the detailed structure, and the density area can be effectively used. Furthermore, the attention area can be converted into gradation as in an actual film image, and other density areas can be converted into gradation within a predetermined density range.
[0063]
Next, a third embodiment will be described.
FIG. 13 shows a function used when compressing the dynamic range on the high density side, d1 is a feature quantity extracted by the first feature quantity extraction means 301, and a solid line is a function used when compressing the dynamic range on the high density side. The dotted line indicates a function used when the dynamic range on the high density side is expanded. The starting point of SLOPE is an arbitrary position d4 between max and d1.
[0064]
FIG. 14 shows a function used when compressing the dynamic range on the low density side, d1 is a feature quantity extracted by the first feature quantity extraction means, and a solid line is a function used when compressing the dynamic range on the low density side, A dotted line indicates a function used when the dynamic range on the low density side is expanded. The starting point of SLOPE is an arbitrary position d5 between min and d1.
[0065]
Next, the operation will be described.
The ratio and SLOPE used for compression and expansion of the high density dynamic range are calculated according to the equations (40), (41), and (42). Here, max indicates either max1 or max2.
[0066]
ratio = (d2−d4) / (max−d4) (40)
If ratio is less than 1,
SLOPE = 1-ratio ─── (41)
If the ratio is greater than 1,
SLOPE = ratio-1 (42)
[0067]
Next, ratio1 and SLOPE1 used for compression and expansion of the low density dynamic range are calculated according to equations (43), (44), and (45).
ratio1 = (d5-d3) / (d5-min) (43)
[0068]
If ratio1 is 1 or less
SLOPE1 = 1-ratio1 (44)
If ratio1 is greater than 1,
SLOPE1 = ratio1-1 ─── (45)
Thereafter, the same processing as in the second embodiment is performed.
[0069]
As described above, according to the third embodiment, since the start point of SLOPE can be arbitrarily moved, there is an effect that the range in which the dynamic range is changed can be adjusted.
[0070]
Next, the storage medium according to the present invention will be described.
1 and 8 may be configured in hardware, or may be configured as a computer system including a CPU, a memory, and the like. When configured in a computer system, the memory constitutes a storage medium according to the present invention. This storage medium stores a program for executing a processing procedure including the flowchart of FIG. 2 for controlling the above-described operation.
[0071]
As this storage medium, a semiconductor memory such as ROM or RAM, an optical disk, a magneto-optical disk, a magnetic medium or the like may be used, and these are configured as a CD-ROM, FD, magnetic card, magnetic tape, nonvolatile memory card, or the like. May be used.
[0072]
Therefore, this storage medium is used in other systems or apparatuses other than the systems shown in FIGS. 1 and 8, and the system or computer reads out and executes the program code stored in the storage medium, as described above. The functions equivalent to those of the embodiments can be realized, the same effects can be obtained, and the object of the present invention can be achieved.
[0073]
Further, when an OS or the like running on the computer performs part or all of the processing, or an extended function board in which a program code read from a storage medium is inserted into the computer or an extended function connected to the computer Even when a CPU or the like provided in the extended function board or the extended function unit performs part or all of the processing based on the instruction of the program code after being written in the memory provided in the unit, it is equivalent to each embodiment. In addition to realizing the above functions, the same effect can be obtained and the object of the present invention can be achieved.
[0074]
【The invention's effect】
According to the present invention, it is possible to calculate a subject area necessary for calculating a feature amount necessary for gradation conversion processing of an original image.
[Brief description of the drawings]
FIG. 1 is a block diagram showing an image processing apparatus according to a first embodiment of the present invention.
FIG. 2 is a flowchart showing a flow of processing according to the first embodiment.
FIG. 3 is a characteristic diagram showing a density histogram.
FIG. 4 is a characteristic diagram showing a density histogram.
FIG. 5 is a characteristic diagram showing a form of a function f1 ().
FIG. 6 is a characteristic diagram showing the form of a function f2 ().
FIG. 7 is a characteristic diagram showing a gradation conversion curve.
FIG. 8 is a block diagram showing an image processing apparatus according to a second embodiment of the present invention.
FIG. 9 is a configuration diagram showing a lung side image from which a void region is deleted.
FIG. 10 is a characteristic diagram showing a form of a function f1.
FIG. 11 is a characteristic diagram showing a form of a function f3.
FIG. 12 is a characteristic diagram showing a gradation conversion curve.
FIG. 13 is a characteristic diagram showing a form of a function f2.
FIG. 14 is a characteristic diagram showing a form of a function f4.
FIG. 15 is a characteristic diagram showing the form of a conventional function.
[Explanation of symbols]
101 Clearing means
102 Histogram creation means
103 Concavity and convexity extraction means
104 Analysis means
106 Smoothing means
107 Change means
301 1st feature-value extraction means
302 Clearing means
303 second feature amount extraction means
304 3rd feature-value extraction means
305 Smoothing means
306 Change means

Claims

In an image processing apparatus that performs gradation conversion processing of an original image obtained by irradiating a subject with X-rays,
From a coordinate on the original image of a pixel having a pixel value equal to or greater than a value of a pixel representative of a void region that is a region through which X-rays pass without passing through the subject to a region having a certain width on the original image Means for extracting an area in the original image that does not belong as a subject area ;
Means for obtaining a histogram indicating the frequency of appearance of pixel values in the subject area ;
Calculating means for calculating a feature amount based on the histogram;
Means for performing gradation conversion processing on the original image with a gradation conversion curve determined so that the pixel having the feature value has a target value after gradation conversion processing;
An image processing apparatus comprising:

A convex region of the histogram is extracted by comparing the appearance frequencies at three points located at equal intervals on the horizontal axis of the histogram, and a pixel value indicating the maximum appearance frequency in the convex region is first The image processing apparatus according to claim 1, wherein a pixel average value equal to or less than the first pixel value in the image in the subject region is used as the feature amount.

In the histogram, the appearance frequency of the pixel value indicating a value equal to or lower than the first pixel value and the value of the pixel value are summed, and the product sum is calculated based on the appearance frequency of a pixel value indicating a value equal to or lower than the first pixel value. The image processing apparatus according to claim 2, wherein the feature amount is calculated by dividing.

Smoothing means for obtaining a smooth image of the original image;
Means for changing the dynamic range of the original image by adding, to the original image, an image obtained by gradation-converting the smoothed image so that the maximum value of the pixels in the subject region becomes a target value; ,
The image processing apparatus according to claim 1, further comprising:

In an image processing method for performing gradation conversion processing of an original image obtained by irradiating a subject with X-rays,
From a coordinate on the original image of a pixel having a pixel value equal to or greater than a value of a pixel representative of a void region that is a region through which X-rays pass without passing through the subject to a region having a certain width on the original image Extracting a region in the original image that does not belong as a subject region ;
Obtaining a histogram indicating the frequency of appearance of pixel values in the subject area ;
Calculating a feature amount based on the histogram;
A step of performing gradation conversion processing on the original image with a gradation conversion curve determined so that a pixel having the feature value has a target value after gradation conversion processing;
An image processing method comprising:

In order to perform gradation conversion processing of the original image obtained by irradiating the subject with X-rays,
From a coordinate on the original image of a pixel having a pixel value equal to or greater than a value of a pixel representative of a void region that is a region through which X-rays pass without passing through the subject to a region having a certain width on the original image Means for extracting an area in the original image that does not belong as a subject area ;
Means for obtaining a histogram indicating the frequency of appearance of pixel values in the subject area ;
Means for calculating a feature amount based on the histogram;
Computer-readable recording of a program for causing the original image to function as means for performing gradation conversion processing with a gradation conversion curve determined so that the pixel having the feature value has a target value after gradation conversion processing Possible storage medium.