JP3937660B2

JP3937660B2 - Image processing method and computer-readable storage medium

Info

Publication number: JP3937660B2
Application number: JP25418499A
Authority: JP
Inventors: 弘之新畠
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 1998-09-18
Filing date: 1999-09-08
Publication date: 2007-06-27
Anticipated expiration: 2019-09-08
Also published as: JP2000155838A

Description

【０００１】
【発明の属する技術分野】
本発明は、Ｘ線画像等の画像のダイナミックレンジ圧縮処理機能を有する画像処理方法及びコンピュータ読み取り可能な記憶媒体に関するものである。
【０００２】
【従来の技術】
例えば、Ｘ線胸部画像は、Ｘ線が透過しやすい肺野の画像及びＸ線が非常に透過しにくい縦隔部の画像より構成されるため、画素値の存在するレンジが非常に広い。このため、肺野及び縦隔部の両方を同時に観察することが可能なＸ線胸部画像を得ることは困難であるとされてきた。
【０００３】
そこで、この問題を回避する方法として、「ＳＰＩＥＶｏｌ．６２６ＭｅｄｉｃｉｎｅＸＩＶ／ＰＡＣＳＩＶ（１９８６）」に記載された方法（以下、「方法１」と言う）がある。この方法１は、処理後の画素値Ｓ_D、オリジナル画素値（入力画素値）Ｓ_org、オリジナル画像（入力画像）の低周波数画像の画素値Ｓ_US、定数Ａ，Ｂ，Ｃ（例えば、Ａ＝２，Ｂ＝０．７）を持って、
Ｓ_D＝Ａ［Ｓ_org−Ｓ_US］＋Ｂ［Ｓ_US］＋Ｃ・・・・（１）
なる式（１）で表されるものである。
【０００４】
上述のような方法１では、高周波成分（第１項）、低周波成分（第２項）の重み付けを変えることが可能で、例えば、定数Ａ，ＢをＡ＝２，Ｂ＝０．７とした場合、高周波成分を強調し、且つ全体のダイナミックレンジを圧縮する効果が得られる。これは、ある放射線医師達によれば、本処理なしの画像と比較して、診断に有効であるという評価が得られている。
【０００５】
また、第２５０９５０３特許公報等には、処理後の画素値Ｓ_D、オリジナル画素値（入力画素値）Ｓ_org、オリジナル画像（入力画像）のＹ方向のプロファイルの平均プロファイルＰｙとＸ方向のプロファイルの平均プロファイルＰｘを持って、
Ｓ_D＝Ｓ_org＋Ｆ［Ｇ（Ｐｘ，Ｐｙ）］・・・・（２）
なる式（２）で表される方法（以下、「方法２」と言う）が記載されている。
【０００６】
ここで、式（２）での関数ｆ（ｘ）が有する特性について説明すると、まず、「ｘ＞Ｄｔｈ」ではｆ（ｘ）が「０」となり、「０≦ｘ≦Ｄｔｈ」ではｆ（ｘ）が切片を「Ｅ」、傾きを「Ｅ／Ｄｔｈ」として単調減少するものとなり、次の式（３）で表される。
ｉｆ（ｘ≦Ｄｔｈ）ｆ（ｘ）＝Ｅ−（Ｅ／Ｄｔｈ）ｘ
ｅｌｓｅｆ（ｘ）＝０・・・・（３）
また、”Ｐｘ”及び”Ｐｙ”は、
Ｐｙ＝（ΣＰｙｉ）／ｎ・・・・（４）
Ｐｘ＝（ΣＰｘｉ）／ｎ・・・・（５）
なる式（４）及び（５）で表される。
【０００７】
式（４）及び（５）において、”ｉ”は”ｉ＝１〜ｎ”であり、”Ｐｙｉ”及び”Ｐｘｉ”はプロファイルであり、また、式（２）での”Ｇ（Ｐｘ，Ｐｙ）”は、例えば、
Ｇ（Ｐｘ，Ｐｙ）＝ｍａｘ（Ｐｘ，Ｐｙ）・・・・（６）
で示される。
【０００８】
上述のような方法２では、低周波画像の画素値で”Ｄｔｈ”以下の濃度レンジが圧縮される。
【０００９】
また、方法２（第２５０９５０３特許公報等に記載の方法）と同様の方法（以下、「方法３」と言う）が「日本放射線技術学会雑誌第４５巻第８号１９８９年８月１０３０頁阿部ほか」に記載されている。この方法３は、処理後の画素値Ｓ_D、オリジナル画素値（入力画素値）Ｓ_org、オリジナル画像（入力画像）をマスクサイズＭ×Ｍ画素で移動平均をとった時の平均画素値Ｓ_US、単調減少関数ｆ（ｘ）を持って、
Ｓ_D＝Ｓ_org＋ｆ（Ｓ_US）・・・・（７）
Ｓ_US＝ΣＳ_org／Ｍ² ・・・・（８）
なる式（７）及び（８）で表されるものである。
【００１０】
方法３は、式（２）によって表される方法２とは低周波画像の作成方法が異なり、方法２では１次元データで低周波画像を作成していたのに対し、２次元データで低周波画像を作成する方法である。このような方法３によっても、低周波画像の画素値で”Ｄｔｈ”以下の濃度レンジが圧縮される。
【００１１】
また、第２６６３１８９特許公報等には、処理後の画素値Ｓ_D、オリジナル画素値（入力画素値）Ｓ_org、オリジナル画像（入力画像）をマスクサイズＭ×Ｍ画素で移動平均をとった時の平均画素値Ｓ_US、単調増加関数ｆ１（ｘ）を持って、
Ｓ_D＝Ｓ_org＋ｆ１（Ｓ_US）・・・・（７）
Ｓ_US＝ΣＳ_org／Ｍ² ・・・・（８）
なる式（７）及び（８）で表される方法（以下、「方法４」と言う）が記載されている。
【００１２】
ここで、式（７）での関数ｆ１（ｘ）が有する特性について説明すると、まず、「ｘ＜Ｄｔｈ」ではｆ（ｘ）が「０」となり、「Ｄｔｈ≦ｘ」ではｆ（ｘ）が切片を「Ｅ」、傾きを「Ｅ／Ｄｔｈ」として単調減少するものとなり、次の式（１１）で表される。
ｆ１［ｘ］＝Ｅ−（Ｅ／Ｄｔｈ）Ｘ・・・・（１１）
【００１３】
上述のような方法４では、低周波画像の画素値Ｄｔｈ以上の画素値（濃度値）が圧縮される。このときのアルゴリズムについては、方法３（「日本放射線技術学会雑誌第４５巻第８号１９８９年８月１０３０頁阿部ほか」等に記載の方法）と同様である。
【００１４】
【発明が解決しようとする課題】
しかしながら、上述したような従来の方法では、取得画像（入力画像）を階調変換してフィルムやＣＲＴに対して出力した場合、視認確認できない濃度値もしくは輝度値の範囲に階調変換後の注目領域（被写体領域等）の濃度値もしくは輝度値がかかるのを防ぐことを目的としていた。
【００１５】
このため、従来の方法では、圧縮開始位置や圧縮率が固定であり、ダイナミックレンジ圧縮を施した取得画像を階調変換した場合、注目領域の濃度値もしくは輝度値が視認領域に入らないことがあった。これは診断等に支障をきたす問題につながる。
【００１６】
また、取得画像によっては、被写体部分の圧縮が必要以上に行なわれてしまい、本来被写体部分が保持している情報量が無意味に損出される場合があった。これは、診断能等が低下する問題につながる。
【００１７】
さらに、ダイナミックレンジを拡大するという発想が無かったので、フィルムやＣＲＴ上に画像を出力する場合に、注目領域のダイナミックレンジを有効に出力することができないという問題があった。すなわち、注目領域のダイナミックレンジを所定の領域に変換するという思想がなかった。
【００１８】
本発明は上記のような問題を解決するためになされたもので、高周波成分の振幅を保った状態で、注目領域のダイナミックレンジを所定の濃度領域に変換できるようにすることを目的とする。
【００１９】
【課題を解決するための手段】
本発明の画像処理方法は、原画像データに対し諧調変換処理及び周波数処理を施す画像処理方法であって、前記原画像データにおける被写体内の領域から、階調変換曲線を生成するための値を抽出する工程と、前記原画像データに基づく高周波成分を生成する工程と、前記高周波成分を、前記階調変換曲線の傾きに基づいて変更する工程と、前記変更した高周波波成分に基づき画像データを生成する工程とを備えることを特徴とする。
本発明のコンピュータ読み取り可能な記憶媒体は、原画像データに対し諧調変換処理及び周波数処理を施す画像処理方法をコンピュータに実行させるためのプログラムを記録したコンピュータ読み取り可能な記憶媒体であって、前記原画像データにおける被写体内の領域から、階調変換曲線を生成するための値を抽出する工程と、前記原画像データに基づく高周波成分を生成する工程と、前記高周波成分を、前記階調変換曲線の傾きに基づいて変更する工程と、前記変更した高周波波成分に基づき画像データを生成する工程とをコンピュータに実行させるためのプログラムを記憶したことを特徴とする。
【００４３】
【発明の実施の形態】
以下、本発明の実施の形態を図面と共に説明する。
【００４４】
（第１の実施の形態）
図１は、本発明の実施の形態による画像処理装置１００を示す。
この画像処理装置１００は、濃度値変換機能を有するＸ線画像の画像処理装置であり、前処理回路１０６、濃度値変換部１１３，ＣＰＵ１０８、メインメモリ１０９、操作パネル１１０、画像表示装置１１１を備えており、ＣＰＵバス１０７を介して互いにデータ授受されるようになされている。
【００４５】
また、画像処理装置１００は、前処理回路１０６に接続されたデータ収集回路１０５と、データ収集回路１０５に接続された２次元Ｘ線センサ１０４及びＸ線発生回路１０１とを備えており、これらの各回路はＣＰＵバス１０７にも接続されている。
【００４６】
図２は画像処理装置の処理の流れを示すフローチャートである。
図３（ａ）は入力画像を示し、図３（ｂ）は上記入力画像からＸ線の照射領域とす抜け領域（Ｘ線が素通りした領域）とを削除した図である。
図４は階調変換曲線上で視認できる濃度域と注目領域の濃度域との関係を示す図であり、横軸が入力画像の濃度値、縦軸が出力画像の濃度値である。
【００４７】
上記構成による画像処理装置１００において、メインメモリ１０９は、ＣＰＵ１０８の処理に必要な各種のデータなどが記憶されると共に、ＣＰＵ１０８の作業用としてのワークメモリを含む。ＣＰＵ１０８は、メインメモリ１０９を用いて、操作パネル１１０からの操作に従って装置全体の動作制御等を行う。これにより画像処理装置１００は、以下のように動作する。
【００４８】
まず、Ｘ線発生回路１０１は、被検査体１０３に対してＸ線ビーム１０２を放射する。このＸ線ビーム１０２は、被検査体１０３を減衰しながら透過して、２次元Ｘ線センサ１０４に到達し、２次元Ｘ線センサ１０４よりＸ線画像として出力される。ここでは、２次元Ｘ線センサ１０４から出力されるＸ線画像を、例えば図３（ａ）のような膝、胸椎画像とする。
【００４９】
データ収集回路１０５は、２次元Ｘ線センサ１０４から出力されたＸ線画像を電気信号に変換して前処理回路１０６に供給する。前処理回路１０６は、データ収集回路１０５からの信号（Ｘ線画像信号）に対して、オフセット補正処理やゲイン補正処理等の前処理を行う。前処理が行われたＸ線画像信号は入力画像として、ＣＰＵ１０８の制御により、ＣＰＵバス１０７を介してメインメモリ１０９、照射領域抽出回路１１２及び濃度値変換部１１３のす抜け削除回路１１３ａに転送される。
【００５０】
濃度値変換部１１３において、１１３ａはす抜け領域とす抜け領域と一定幅で接する体領域とを削除するす抜け削除回路、１１３ｂはす抜け削除回路１１３ａで削除されなかった領域から、濃度値の最大値と最小値を算出する最高値最小値抽出回路、１１３ｃは入力画像の平滑化画像を作成する平滑化画像作成回路、１１３ｄは入力画像と平滑化画像との差分から高周波成分を作成する高周波成分作成回路、１１３ｅは最大値最小値抽出回路１１３ｂで抽出された最大値と最小値に基づき平滑化画像の濃度値を変換する濃度値変換回路、１１３ｆは濃度値変換回路１１３ｅで変換された画像に高周波成分作成回路１１３ｄで作成した高周波を足し込む高周波成分足し込み回路である。
【００５１】
次に濃度値変換部１１３の動作について図２のフローチャートにより説明する。
ＣＰＵバス１０７を介して前処理回路１０６で処理された入力画像をＣＰＵ１０８の制御により受信した照射領域抽出回路１１２は、入力画像中の照射領域を抽出する（ステップＳ２０１）。また、同時に入力画像を受信したす抜け削除回路１１３ａは照射領域外及び照射領域内のす抜け領域とす抜け領域と一定間隔内で接する体領域とを例えば０画素で置き換える（ステップＳ２０２）。具体的には次の式（１２）で表されるような画像の変換を行う。
【００５２】
【数１】

【００５３】
ここで、ｆ（ｘ，ｙ）は画像データを示し、ｆ１（ｘ，ｙ）はす抜け領域及びす抜け領域と一定間隔内で接する体領域を削除した後の画像を示す。Ｓｇｎ（ｘ，ｙ）は以下の式（４）のように表される。Ｔｈ１は実験により定められる定数、ｄ１，ｄ２は体領域を削除する幅を決める定数である。
ｓｇｎ（ｘ，ｙ）＝０ｆ（ｘ，ｙ）≧Ｔｈ１のとき
ｓｇｎ（ｘ，ｙ）＝１その他・・・・（１３）。
図３（ｂ）は、入力画像の照射領域外及びす抜け領域を０で置き換えた場合の画像である。
【００５４】
次に、最大最小値抽出回路１１３ｂでは、画像濃度値の最大値（図４の”ｄ２”）と最小値（図４の”ｄ３”）を算出する（ステップＳ２０３）。ここで、最大値、最小値を算出するのに、照射領域及びす抜け領域を除去した領域の平滑化画像から抽出してもよい。次に平滑化画像作成回路１１３ｃは次式に従い平滑化画像を作成する（ステップＳ２０４）。ｆｕｓ（ｘ，ｙ）は入力画像ｆ０（ｘ，ｙ）の平滑化（低周波）画像の画素値であり、例えば、式（１４）〜（１８）、又は式（１９）で示される。
【００５５】

【００５６】
”Ｄ（ｘ，ｙ）”は、円盤状フィルタであり、ｒ１を任意の定数とし、入力画像に応じて選択される。

【００５７】
ここで得られたｆｕｓ（ｘ，ｙ）のプロファイルはエッジ構造を保存しているものであり、従来ダイナミックレンジ圧縮の欠点であるオーバーシュート、アンダーシュートが起きないものである。
【００５８】
平滑化画像は、例えば式（１７）で示される平均濃度を用いてもよいし、次の式（１９）を用いてもよい。ここでｄは定数である。
【００５９】
【数２】

【００６０】
また、例えばＥｒｏｓｉｏｎ、Ｄａｉ１ａｔｉｏｎ、Ｏｐｅｎｉｎｇ、Ｃ１ｏｓｉｎｎｇ等のモルフォロジカルフィルタを用い平滑化画像を作成してもよい。
【００６１】
次に、高周波成分作成回路１１３ｃでは、入力画像ｆ０（ｘ，ｙ）と平滑化画像ｆｕｓ（ｘ，ｙ）とから次式に従い高周波画像ｆｈ（ｘ，ｙ）を作成する（ステップＳ２０５）。次に濃度値変換回路１１３ｅは、最大最小値値抽出回路１１３ｂで抽出した最大値（ｄ３）、最小値（ｄ２）、そして階調変換後の視認濃度値で決まる濃度値ｄ１，ｄ４に基づき次式により濃度値を変換し、濃度値変換した平滑化画像Ｓｕｓ０×（ｘ，ｙ）を作成する（ステップＳ２０６）。
【００６２】
図４において、ＳｉｇｈｔｍａｘとＳｉｇｈｔｍｉｎが視認できる最大最小濃度値であり、ｄ１，ｄ４に対応する。平滑化画像の濃度値がｍｉｎ（平滑化画像が取り得る最小値）からｄ１以下の場合、

【００６３】
ｄ２＜Ｓｕｓ（ｘ，ｙ）≦ｄ３の場合、

【００６４】
ｄ３＜Ｓｕｓ（ｘ，ｙ）≦ｍａｘの場合、

【００６５】
そして、濃度値変換した画像Ｓｕｓｕｏ（ｘ，ｙ）に高周波画像ｆｈ（ｘ，ｙ）を足し込み、最終画像ｆｅ（ｘ，ｙ）を得る（ステップＳ２０７）。
ｆｅ（ｘ，ｙ）＝Ｓｕｓｕ０（ｘ，ｙ）十ｆｈ（ｘ，ｙ）・・・（２３）
最終的には、得られた画像ｆｅ（ｘ，ｙ）を階調変換回路１１４で階調変換してフィルム出力又は画像表示装置１１１に表示する。
【００６６】
尚、本実施の形態では、す抜け領域がある画像に対する場合について説明したが、す抜け領域がない場合は、す抜け削除回路１１３ａを経由せず、照射領域内の画像から、最大最小値抽出回路１１３ｂにより最大値と最小値を抽出すればよい。
【００６７】
本実施の形態によれば、高周波成分を保持したまま、注目領域の画像を一定幅に圧縮、拡大することが可能であり、画像情報を有効に使える効果がある。さらに、上記一定幅をフィルムなどの視認領域とすれば、フィルム上の視認領域の幅に注目領域を拡大縮小することができる効果が有る。また、最大値、最小値を平滑化画像から抽出する場合は、濃度値変換部１１３での変換をより精度よく行える効果が有る。
【００６８】
（第２の実施の形態）
図５は、本発明の実施の形態による画像処理装置を適用したＸ線撮影装置６００の構成を示す。
【００６９】
Ｘ線撮影装置６００は、画像処理機能を有するものであり、前処理回路６０６、ＣＰＵ６０８、メインメモリ６０９、操作パネル６１０、画像表示器６１１、及び画像処理回路６１２を備えており、これらの各回路はＣＰＵバス６０７を介して互いにデータ授受されるようになされている。
【００７０】
また、Ｘ線撮影装置６００は、前処理回路６０６に接続されたデータ収集回路６０５と、データ収集回路６０５に接続された２次元Ｘ線センサ６０４及びＸ線発生回路６０１とを備えており、これらの各回路もＣＰＵバス６０７に接続されている。
【００７１】
画像処理回路６１２は、階調変換のための特徴量を算出する特徴抽出回路６１２ａと、後述する第１の階調変換回路６１２ｄの階調変換曲線からダイナミックレンジを変更する範囲と変更量を算出する制御回路６１２ｂと、制御回路６１２ｂで算出されたダイナミックレンジを変更する範囲と変更量から原画像のダイナミックレンジを変更するＤＲＣ回路６１２ｃと、ＤＲＣ回路６１２ｃにてダイナミックレンジが変更された原画像の階調変換を行なう第１の階調変換回路６１２ｄとを備えている。
【００７２】
ＤＲＣ回路６１２ｃは、原画像の階調変換を行なう第２の階調変換回路６１３と、第２の階調変換回路６１３にて階調変換された画像に対して原画像の高周波成分を足し込む高周波成分調整回路とを備えている。
【００７３】
図６は、例えば、頸椎側面で撮影を行なった場合に得られる画像（Ｘ線撮影画像）において、特徴量を抽出する領域ａを示したものである。
【００７４】
図７は、第１の階調変換回路６１２ｄの階調変換曲線を示したものであり、この図７において、横軸は画素値を示し、縦軸は濃度値を示す。”Ｓｌ”は、濃度上で視認識限界の下限濃度値に対応する画素値を示し、”Ｓｈ”は、濃度上で視認識限界の上限濃度値に対応する画素値を示す。”Ｓｍｉｎ”及び”Ｓｍａｘ”は、注目領域の最小画素値及び最大画素値を示す。”Ｓａ”は、特徴抽出回路６１２ａにて算出された特徴量に対応する画素値を示し、”ｄｌ”及び”ｄｈ”はそれぞれ、ダイナミックレンジの変更しない範囲の下限濃度値及び上限濃度値に対応する画素値を示す。
【００７５】
図８は、第２の階調変換回路６１３の階調変換曲線を示したものであり、この図８において、横軸は入力画素値を示し、縦軸は出力画素値を示す。”ｄｌ”、”ｄｈ”、及び”Ｓａ”は、図７での画素値と同じ画素値を示す。
図８に示す階調変換曲線ａは、直線成分のみから構成され、同図に示す階調変換曲線ｂは、階調変換曲線ａの微分値が連続となるようななだらかな曲線から構成されている。
【００７６】
ここで、上述のようなＸ線撮影装置６００において、メインメモリ６０９は、ＣＰＵ６０８による本装置全体の動作制御のための処理に必要な処理プログラムや各種のデータ等が記憶されるものであると共に、ＣＰＵ６０８の作業用メモリとしてのワークメモリを含むものである。
したがって、ＣＰＵ６０８は、メインメモリ６０９及びメインメモリ６０９に記憶された情報を用いて、操作パネル６１０からの操作に従った本装置全体の動作制御等を行なう。これにより、Ｘ線撮影装置６００は次のように動作する。
【００７７】
先ず、Ｘ線発生回路６０１が被写体６０３に対してＸ線ビーム６０２を放射する。Ｘ線発生回路６０１から放射されたＸ線ビーム６０２は、被写体６０３を減衰しながら透過して、２次元Ｘ線センサ６０４に到達し、２次元Ｘ線センサ６０４によりＸ線画像として出力される。この２次元Ｘ線センサ６０４から出力されるＸ線画像は、例えば、人体部画像である。
【００７８】
データ収集回路６０５は、２次元Ｘ線センサ６０４から出力されたＸ線画像を電気信号に変換して前処理回路６０６に供給する。前処理回路６０６は、データ収集回路６０５からの信号（Ｘ線画像信号）に対して、オフセット補正処理やゲイン補正処理等の前処理を施す。この前処理回路６０６にて前処理が施されたＸ線画像信号は入力画像として、ＣＰＵ６０８の制御により、ＣＰＵバス６０７を介して、メインメモリ６０９及び画像処理回路６１２にそれぞれ転送される。
【００７９】
画像処理回路６１２は、例えば、図９に示すフローチャートに従って次のように動作する。
【００８０】
先ず、ＣＰＵバス６０７を介して前処理回路６０６から供給された入力画像（前処理が施された原画像）ｆ０（ｘ，ｙ）を、ＣＰＵ６０８の制御により受信した特徴抽出回路６１２ａは、階調変換のための特徴量Ｓ_aを算出する（ステップＳ７０１）。
ここでの特徴量の抽出方法は、入力画像ｆ０（ｘ，ｙ）での被写体の部位毎に異なり、本出願人から複数の方法が提案されている。例えば、入力画像ｆ０（ｘ，ｙ）が、頸椎の撮影により得られた画像である場合には、特願平１０−２７２２８４号に記載された方法を用いる。この方法は、特徴量Ｓａとして、図６の”ａ”に示した領域内の画素値の平均値を算出する方法である。
【００８１】
次に、第１の階調変換回路６１２ｄでの階調変換曲線を、特徴抽出回路６１２ａにて得られた特徴量Ｓａを用いて規定する（ステップＳ７０２）。
例えば、図７に示すように、特徴量Ｓａの画素値が濃度値１．０に変更されるような階調変換曲線を規定する。
【００８２】
次に、制御回路６１２ｂは、ステップＳ７０２にて規定された階調変換曲線から入力画像ｆ０（ｘ，ｙ）のダイナミックレンジを変更する範囲と変更量を算出する（ステップＳ７０３）。この算出方法としての一例を次に述べる。
【００８３】
例えば、通常の日本人は、画像での視認識できるの濃度値が、濃度値０．２から３．０の範囲となっている。
そこで、先ず、ステップＳ７０２にて規定された階調変換曲線に基づき、濃度値０．２、３．０に階調変換される画素値Ｓｌ、Ｓｈを算出する。
次に、ダイナミックレンジを変更しない幅を、
ｄｌ＝Ｓａ−ｈｌ・・・・（２４）
ｄｈ＝Ｓａ＋ｈｈ・・・・（２５）
なる式（２４）及び（２５）に従って算出する。
式（２４）及び（２５）において、”ｈｌ”及び”ｈｈ”は定数であり、”ｄｌ”及び”ｄｈ”はダイナミックレンジを変更する起点となる画素である。
【００８４】
そして、”ｄｌ”以下の画素値、”ｄｌ”以上の画素値の変更量Ｒｌ，Ｒｈを、
Ｒｌ＝（Ｓｌ−ｄｌ）／（Ｓｍｉｎ−ｄｌ）・・・・（２６）
Ｒｈ＝（Ｓｈ−ｄｈ）／（Ｓｍａｘ−ｄｈ）・・・・（２７）
なる式（２６）及び（２７）に従って算出する。
式（２６）及び（２７）において、”Ｓｍｉｎ”及び”Ｓｍａｘ”は、被写体領域（注目領域）の画素値のうちの最小値及び最大値を示す。すなわち、ｄｌ以下の画素値のダイナミックレンジをＲｌ倍すると、画素値Ｓｍｉｎは画素値Ｓｌに変更される。同様に、画素値Ｓｍａｘは画素値Ｓｈに変更される。換言すれば、起点をｄｌ，ｄｈとし、変更量をＲｌ，Ｒｈとすると、階調変換後の被写体領域の濃度値は、視野領域（ここでは例えば、濃度値０．２〜３．０）の幅と一致する。
【００８５】
上述のように、変更量Ｒｌ，Ｒｈを求めることにより、ダイナミックレンジを圧縮しても、撮影画像の注目領域の全てを出力画像上に再現することができる。したがって、特徴抽出の不安定性を吸収することができ、診断能力の向上を図ることができる。
【００８６】
また、変更量Ｒｌ，Ｒｈに対する拘束条件として、
Ｒｌ＜Ｃｌ・・・・（２８）
Ｒｈ＜Ｃｈ・・・・（２９）
なる式（２８）及び（２９）に示すような条件を設けるようにしてもよい。
式（２８）及び（２９）において、”Ｃｌ”及び”Ｃｈ”はそれぞれ定数を示す。これらの定数Ｃｌ，Ｃｈ、及び上述した定数ｈｌ，ｈｈは、例えば、撮影対象となる被写体の部位毎に実験的に決まる定数である。
このように、拘束条件を設けることによって、圧縮のしすぎを防ぐことができるため、圧縮のしすぎによる診断能力の低下を確実に防止することができる。
【００８７】
上述のようにして、ダイナミックレンジを変更する起点となる画素ｄｌ，ｄｈ、及び変更量Ｒｌ，Ｒｈが求まると、これらのｄｌ，ｄｈ，Ｒｌ，Ｒｈに基づき、図８に示したような第２の階調変換回路６１３の階調変換曲線を規定する。
図８に示す階調変換曲線ａでの画素値ｄｌ以下の傾きは”Ｒｌ”であり、画素値ｄｌから画素値ｄｈの範囲の傾きは”１”であり、画素値ｄｈ以上の傾きは”Ｒ_h”である。一方、同じく図８に示す階調変換曲線ｂは、特願昭１１−７６８８２号等に記載されている技術によって、階調変換曲線ａの微分値が連続となるようになだらかにした曲線である。第２の階調変換回路６１３の階調変換曲線として、階調変換曲線ａを用いる場合には、微分不連続点で偽輪郭が発生することがあるが、階調変換曲線ｂを用いる場合には、当該偽輪郭が発生することはない。ここでは一例として、階調変換曲線ｂを第２の階調変換回路６１３で用いるものとする。
【００８８】
次に、ＤＲＣ回路６１２ｃは、以降の処理で、高周波成分を保持するか否かを選択する（ステップＳ７０４）。
高周波成分を保持せずに処理実行した場合に得られる結果画像は、低周波成分が重要である軟部組織等を診断する場合に用いて有効である。一方、高周波成分を保持して処理実行した場合に得られる結果画像は、高周波成分が重要である骨部、肺野等を診断する場合に用いて有効である。すなわち、処理対象の画像が軟部組織等の画像であれば、高周波成分を保持しないと選択し、処理対象の画像が骨部、肺野等の画像であれば、高周波成分を保持すると選択することによって、軟部組織や、骨部、肺野等のそれぞれの診断に有効な結果画像が得られる。
このステップＳ７０４での選択は、撮影部位（軟部組織、骨部、肺野等）に基づき自動選択する構成としてもよいし、ユーザがマニュアル選択する構成としてもよい。
【００８９】
ステップＳ７０４での選択の結果、高周波成分を保持しない場合、第２の階調変換回路６１３により、入力画像（原画像）に対して階調変換曲線ｂ（図８参照）を用いた階調変換を施す（ステップＳ７０８）。
そして、第１の階調変換回路６１２ｄにより、第２の階調変換回路６１３での階調変換後の入力画像に対して階調変換を施す（ステップＳ７０９）。
【００９０】
ここで、入力画像（原画像）の座標を”ｘ，ｙ”、原画像の画素値を”ｆ０（ｘ，ｙ）”、第２の階調変換回路６１３の階調変換曲線ｂを”Ｆ（ｘ）”、階調変換曲線ｂを用いた階調変換後の画像の画素値を”ｆ１（ｘ，ｙ）”、第１の階調変換回路６１２ｄの階調変換曲線を”Ｆ１（ｘ）”とすると、結果的に得られる画像（最終画像）ｆ２（ｘ，ｙ）は、
ｆ２（ｘ，ｙ）＝Ｆ１（Ｆ（ｆ０（ｘ，ｙ））・・・・（３０）
なる式（３０）で表される。
この場合、”Ｆ１（Ｆ（ｘ））を、１つの階調変換曲線と考えてよい。
【００９１】
一方、ステップＳ７０４での選択の結果、高周波成分を保持する場合、先ず、高周波成分調整回路６１４により、原画像ｆ０（ｘ，ｙ）の高周波成分ｆｈ（ｘ，ｙ）を、
ｆｈ（ｘ，ｙ）＝ｆ０（ｘ，ｙ）−ｆｕｓ（ｘ，ｙ）・・・・（３１）
なる式（３１）により算出する（ステップＳ７０５）。
式（３１）において、”ｆｕｓ（ｘ，ｙ）”は、平滑化画像であり、マスクサイズの大きさを示す定数ｄを持って、
【００９２】
【数３】

【００９３】
なる式（３２）により算出される。
【００９４】
次に、第２の階調変換回路６１３により、原画像ｆ０（ｘ，ｙ）に対して、
ｆ１（ｘ，ｙ）＝Ｆ（ｆ０（ｘ，ｙ））・・・・（３３）
なる式（３３）で示される階調変換を施し、当該階調変換後の画像ｆ１（ｘ，ｙ）を取得する（ステップＳ７０６）。
【００９５】
次に、高周波成分調整回路６１４により、第２の階調変換回路６１３の階調変換曲線ｂの微係数から係数ｃ（ｘ）を、
【００９６】
【数４】

【００９７】
なる式（３４）により算出する。この式（３４）に示されるように、係数ｃ（ｘ）は、階調変換曲線ｂの傾きから”１”を引いた値である。
そして、係数ｃ（ｘ）に基づいて、ステップＳ７０６にて取得された第２の階調変換回路６１３での階調変換後の画像ｆ１（ｘ，ｙ）に対して、ステップＳ７０６にて取得された原画像ｆ０（ｘ，ｙ）の高周波成分ｆｈ（ｘ，ｙ）を足し込み、その処理後画像ｆ３（ｘ，ｙ）を得る。すなわち、処理後画像ｆ３（ｘ，ｙ）は、

なる式（３５）で表される（ステップＳ７０７）。
【００９８】
そして、第１の階調変換回路６１２ｄにより、ステップＳ７０７にて得られた処理後画像ｆ３（ｘ，ｙ）に対して、
ｆ２（ｘ，ｙ）＝Ｆ１（ｆ３（ｘ，ｙ））・・・・（３６）
なる式（３６）に従った階調変換を行い、この結果得られた画像ｆ２（ｘ，ｙ）を最終画像として得る（ステップＳ７０９）。
【００９９】
尚、ステップＳ７０５の式（３１）における平滑化画像ｆｕｓ（ｘ，ｙ）を、例えば、次の式（３７）〜（４１）で示されるモルフォジ演算を用いて算出するようにしてもよい。

式（３７）〜（４０）において、”Ｄ（ｘ１，ｙ１）”は、円盤状フィルタであり、入力画像に応じて選択される任意の定数ｒ１を持って、

なる式（４１）で表される。
【０１００】
式（３７）〜（４１）により得られる平滑化画像ｆｕｓ（ｘ，ｙ）のプロファイルは、エッジ構造を保存しているものであり、従来ダイナミックレンジ圧縮の欠点であるオーバーシュート、アンダーシュートが起きないものである。
【０１０１】
本実施の形態によれば、入力画像の任意の階調領域の濃度分布幅を圧縮、伸張することができ、且つ階調変換後の高周波成分の振幅を階調変換前の画像の高周波成分の振幅を保持できる効果がある。
また、階調変換曲線で決まる画素値に基づき、画像のダイナミックレンジを変更するので、階調変換後の画像の濃度値を想定して、ダイナミックレンジを変更できるので、階調変換後の濃度値の範囲を一意に調整できる効果がある。
【０１０２】
また、階調変換曲線から視認限界で決まる画素値を求め、その画素値の値に基づきダイナミックレンジを変更する範囲、量を決定できるので、注目領域の画素値の濃度値の範囲を視認範囲にすることができる効果がある。
また、注目領域が視認領域にあることから、注目領域全体を観察することが可能であり、診断能が上がる効果がある。
また、注目領域を視認領域の幅に一致することが可能であり、注目領域の情報量を濃度値として最大限に表現できる。
また、画像のダイナミックレンジを変更しない範囲を有するので、診断等で重要となる画素範囲をダイナミックレンジを変更しない領域とし、従来通りの濃度値で表現し、さらに、従来では視認できなかった領域も濃度値として観察できるため、診断能が上がる効果がある。
また、ダイナミックレンジの変更手段において高周波成分を保持するため、高周波成分の情報量を落とすことがなく、さらに従来では濃度値上で視認できなかった領域の濃度値も観察できるため、診断能が上がる効果がある。
【０１０３】
（他の実施の形態）
上述の実施の形態は、ハード的に実現してもよく、また、ソフト的に構成してもよい。例えば、ＣＰＵやメインメモリ等のメモリ等からなるコンピュータシステムを用いて、図２又は図９に示したフローチャートによる処理を実行するためのプログラムを実現することにより、上述の実施の形態における構成を実現してもよい。
【０１０４】
上記プログラムを記憶する記憶媒体としては、ＲＯＭ、ＲＡＭ等の半導体メモリ、光ディスク、光磁気ディスク、磁気記憶媒体等を用いてよく、これらをＣＤ−ＲＯＭ、ＦＤ、磁気カード、磁気テープ、不揮発性メモリカード等を用いることができる。
【０１０５】
また、コンピュータ上で稼働しているＯＳ等が処理の一部又は全部を行う場合、あるいは記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された拡張機能ボードやコンピュータに接続された拡張機能ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づいて、上記拡張機能ボードや拡張機能ユニットに備わるＣＰＵ等が処理の一部又は全部を行う場合にも、上記実施の形態と同等の機能を実現できると共に、同等の効果を得ることができ、本発明に含まれることは言うまでもない。
【０１０６】
【発明の効果】
以上説明したように、本発明によれば、画像の微細な部分である高周波成分を保持したまま、注目領域の画像を一定幅に圧縮、拡大することが可能であり、これにより、画像情報を有効に使うことができる効果がある。また、上記一定幅をフィルムなどの視認領域とすれば、フィルム上の視認領域の幅に注目領域を拡大縮小することができる効果がある。
【０１０７】
また、画像にす抜け領域がある場合に、す抜け領域を除去することにより、注目領域をさらに精度よく抽出することが可能であり、濃度値変換の精度が上がる効果がある。
【０１０８】
さらに、平滑化画像から最大値、最小値を抽出することにより、濃度値変換をより精度よく行うことができる効果がある。
【図面の簡単な説明】
【図１】第１の実施の形態に係る画像処理装置の構成を示すブロック図である。
【図２】第１の実施の形態での処理手順を示すフローチャートである。
【図３】入力画像と照射領域とす抜け領域とを除去した画像を示す構成図である。
【図４】第１の実施の形態での動作を説明するための特性図である。
【図５】第２の実施の形態に係る画像処理装置を適用したＸ線撮影装置の構成を示すブロック図である。
【図６】特徴量を抽出する領域を説明するための図である。
【図７】上記画像処理装置の第１の階調変換曲線の階調変換曲線を説明するための図である。
【図８】上記画像処理装置の第２の階調変換曲線の階調変換曲線を説明するための図である。
【図９】第２の実施の形態での処理手順を示すフローチャートである。
【符号の説明】
１０５データ収集回路
１０８ＣＰＵ
１０９メインメモリ
１１０操作パネル
１１２照射領域抽出回路
１１３濃度値変換回路
１１３ａす抜け削除回路
１１３ｂ最大最小値抽出回路
１１３ｃ平滑化画像作成回路
１１３ｄ高周波成分作成回路
１１３ｅ濃度値変換回路
１１２ｆ高周波成分足し込み回路[0001]
BACKGROUND OF THE INVENTION
The present invention relates to an image processing method having a dynamic range compression processing function for an image such as an X-ray image and a computer-readable storage medium.
[0002]
[Prior art]
For example, since an X-ray chest image is composed of an image of a lung field through which X-rays are easily transmitted and an image of a mediastinal portion through which X-rays are very difficult to transmit, the range in which pixel values exist is very wide. For this reason, it has been considered difficult to obtain an X-ray chest image capable of simultaneously observing both the lung field and the mediastinum.
[0003]
Therefore, as a method for avoiding this problem, there is a method described in “SPIE Vol. 626 Medicine XIV / PACS IV (1986)” (hereinafter referred to as “Method 1”). This method 1 uses the processed pixel value S _D Original pixel value (input pixel value) S _org , Pixel value S of low frequency image of original image (input image) _US , With constants A, B, C (eg, A = 2, B = 0.7)
S _D = A [S _org -S _US ] + B [S _US ] + C (1)
It is represented by the following formula (1).
[0004]
In the method 1 as described above, the weighting of the high frequency component (first term) and the low frequency component (second term) can be changed. For example, the constants A and B are set to A = 2 and B = 0.7. In this case, the effect of emphasizing the high frequency component and compressing the entire dynamic range can be obtained. According to some radiologists, evaluation that it is effective for diagnosis is obtained compared with an image without this processing.
[0005]
In addition, in the 2509503 patent publication, etc., the pixel value S after processing is described. _D Original pixel value (input pixel value) S _org , Having an average profile Py of the Y-direction profile and an average profile Px of the X-direction profile of the original image (input image),
S _D = S _org + F [G (Px, Py)] (2)
A method represented by the following formula (2) (hereinafter referred to as “method 2”) is described.
[0006]
Here, the characteristics of the function f (x) in the expression (2) will be described. First, f (x) is “0” when “x> Dth”, and f (x) when “0 ≦ x ≦ Dth”. ) Monotonously decreases with the intercept being “E” and the slope being “E / Dth”, and is expressed by the following equation (3).
if (x ≦ Dth) f (x) = E− (E / Dth) x
else f (x) = 0 (3)
“Px” and “Py” are
Py = (ΣPyi) / n (4)
Px = (ΣPxi) / n (5)
It is represented by the following formulas (4) and (5).
[0007]
In equations (4) and (5), “i” is “i = 1 to n”, “Pyi” and “Pxi” are profiles, and “G (Px, Py) in equation (2) ) ”
G (Px, Py) = max (Px, Py) (6)
Indicated by
[0008]
In the method 2 as described above, the density range equal to or lower than “Dth” in the pixel value of the low frequency image is compressed.
[0009]
In addition, the same method (hereinafter referred to as “Method 3”) as Method 2 (the method described in Patent Publication No. 2509503 Patent Gazette) is described in “The Journal of Japanese Society of Radiological Technology, Vol. 45, No. 8, Aug. 1989, page 1030, Abe et al. "It is described in. This method 3 uses the processed pixel value S _D Original pixel value (input pixel value) S _org , Average pixel value S when moving average of original image (input image) with mask size M × M pixels _US Have a monotonically decreasing function f (x),
S _D = S _org + F (S _US (7)
S _US = ΣS _org / M ² .... (8)
It is represented by the following formulas (7) and (8).
[0010]
Method 3 is different from Method 2 represented by Expression (2) in that a low-frequency image is created. In Method 2, a low-frequency image is created using one-dimensional data, whereas a low-frequency image is created using two-dimensional data. This is a method of creating an image. According to the method 3 as described above, the density range equal to or lower than “Dth” is compressed as the pixel value of the low-frequency image.
[0011]
In addition, in the 2663189 patent publication and the like, the pixel value S after processing is described. _D Original pixel value (input pixel value) S _org , Average pixel value S when moving average of original image (input image) with mask size M × M pixels _US Have a monotonically increasing function f1 (x),
S _D = S _org + F1 (S _US (7)
S _US = ΣS _org / M ² .... (8)
A method represented by the following formulas (7) and (8) (hereinafter referred to as “method 4”) is described.
[0012]
Here, the characteristics of the function f1 (x) in Expression (7) will be described. First, f (x) is “0” when “x <Dth”, and f (x) is “Dth ≦ x”. The intercept is “E” and the slope is “E / Dth”, which monotonously decreases, and is expressed by the following equation (11).
f1 [x] = E− (E / Dth) X (11)
[0013]
In the method 4 as described above, a pixel value (density value) equal to or higher than the pixel value Dth of the low-frequency image is compressed. The algorithm at this time is the same as method 3 (the method described in “Journal of Japanese Society of Radiological Technology, Vol. 45, No. 8, August 1989, page 1030, Abe et al.”).
[0014]
[Problems to be solved by the invention]
However, in the conventional method as described above, when the acquired image (input image) is subjected to gradation conversion and output to a film or CRT, attention after gradation conversion is applied to a range of density values or luminance values that cannot be visually confirmed. The object is to prevent the density value or luminance value of an area (subject area, etc.) from being applied.
[0015]
For this reason, in the conventional method, the compression start position and compression rate are fixed, and when the acquired image subjected to dynamic range compression is subjected to gradation conversion, the density value or luminance value of the attention area may not enter the viewing area. there were. This leads to problems that impede diagnosis and the like.
[0016]
Depending on the acquired image, the subject portion may be compressed more than necessary, and the amount of information originally held by the subject portion may be lost meaninglessly. This leads to a problem that the diagnostic ability and the like are lowered.
[0017]
Furthermore, since there was no idea of expanding the dynamic range, there was a problem that the dynamic range of the region of interest could not be output effectively when outputting an image on a film or CRT. That is, there was no idea of converting the dynamic range of the attention area into a predetermined area.
[0018]
The present invention has been made to solve the above-described problems, and an object thereof is to convert the dynamic range of a region of interest into a predetermined density region while maintaining the amplitude of a high-frequency component.
[0019]
[Means for Solving the Problems]
An image processing method of the present invention is an image processing method for performing gradation conversion processing and frequency processing on original image data, and a value for generating a gradation conversion curve from a region within a subject in the original image data. Extracting, generating a high-frequency component based on the original image data, changing the high-frequency component based on the slope of the gradation conversion curve, and image data based on the changed high-frequency wave component And a generating step.
The computer-readable storage medium of the present invention is a computer-readable storage medium recording a program for causing a computer to execute an image processing method for performing gradation conversion processing and frequency processing on original image data. Extracting a value for generating a gradation conversion curve from a region within the subject in the image data, generating a high-frequency component based on the original image data, and converting the high-frequency component into the gradation conversion curve. A program for causing a computer to execute a step of changing based on an inclination and a step of generating image data based on the changed high-frequency wave component is stored.
[0043]
DETAILED DESCRIPTION OF THE INVENTION
Embodiments of the present invention will be described below with reference to the drawings.
[0044]
(First embodiment)
FIG. 1 shows an image processing apparatus 100 according to an embodiment of the present invention.
This image processing apparatus 100 is an X-ray image processing apparatus having a density value conversion function, and includes a preprocessing circuit 106, a density value conversion unit 113, a CPU 108, a main memory 109, an operation panel 110, and an image display device 111. The data is exchanged via the CPU bus 107.
[0045]
The image processing apparatus 100 also includes a data acquisition circuit 105 connected to the preprocessing circuit 106, a two-dimensional X-ray sensor 104 and an X-ray generation circuit 101 connected to the data acquisition circuit 105, and these Each circuit is also connected to the CPU bus 107.
[0046]
FIG. 2 is a flowchart showing a processing flow of the image processing apparatus.
FIG. 3A shows an input image, and FIG. 3B is a diagram in which an X-ray irradiation region and a through region (region through which X-rays pass) are deleted from the input image.
FIG. 4 is a diagram showing the relationship between the density range that can be visually recognized on the tone conversion curve and the density range of the attention area, where the horizontal axis represents the density value of the input image and the vertical axis represents the density value of the output image.
[0047]
In the image processing apparatus 100 configured as described above, the main memory 109 stores various data necessary for the processing of the CPU 108 and includes a work memory for work of the CPU 108. The CPU 108 uses the main memory 109 to perform operation control of the entire apparatus in accordance with an operation from the operation panel 110. As a result, the image processing apparatus 100 operates as follows.
[0048]
First, the X-ray generation circuit 101 emits an X-ray beam 102 to the object 103 to be inspected. The X-ray beam 102 passes through the object 103 while being attenuated, reaches the two-dimensional X-ray sensor 104, and is output as an X-ray image from the two-dimensional X-ray sensor 104. Here, the X-ray image output from the two-dimensional X-ray sensor 104 is, for example, a knee or thoracic vertebra image as shown in FIG.
[0049]
The data acquisition circuit 105 converts the X-ray image output from the two-dimensional X-ray sensor 104 into an electrical signal and supplies it to the preprocessing circuit 106. The preprocessing circuit 106 performs preprocessing such as offset correction processing and gain correction processing on the signal (X-ray image signal) from the data acquisition circuit 105. The preprocessed X-ray image signal is transferred as an input image to the main memory 109, the irradiation area extraction circuit 112, and the omission deletion circuit 113a of the density value conversion unit 113 via the CPU bus 107 under the control of the CPU 108. The
[0050]
In the density

value conversion unit

113, 113 a is a omission area and a omission area that deletes the omission area and a body area that is in contact with a certain width, and 113 b is an area of density value from the area that is not deleted by the omission area deletion circuit 113 a. A maximum / minimum value extraction circuit for calculating the maximum value and the minimum value, 113c is a smoothed image creation circuit for creating a smoothed image of the input image, and 113d is a high frequency for creating a high frequency component from the difference between the input image and the smoothed image. A component creation circuit 113e is a density value conversion circuit that converts the density value of the smoothed image based on the maximum value and minimum value extracted by the maximum value / minimum value extraction circuit 113b, and 113f is an image converted by the density value conversion circuit 113e. This is a high frequency component adding circuit that adds the high frequency generated by the high frequency component generating circuit 113d.
[0051]
Next, the operation of the density value converter 113 will be described with reference to the flowchart of FIG.
The irradiation area extraction circuit 112 that has received the input image processed by the preprocessing circuit 106 via the CPU bus 107 under the control of the CPU 108 extracts the irradiation area in the input image (step S201). At the same time, the omission deletion circuit 113a that has received the input image replaces the omission area and the omission area within the irradiation area and the body area that contacts the omission area within a predetermined interval with, for example, 0 pixels (step S202). Specifically, image conversion as represented by the following equation (12) is performed.
[0052]
[Expression 1]

[0053]
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 the following equation (4). 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 Others (13).
FIG. 3B is an image in the case where the outside of the irradiation area and the omission area of the input image are replaced with 0.
[0054]
Next, the maximum / minimum value extraction circuit 113b calculates the maximum value (“d2” in FIG. 4) and the minimum value (“d3” in FIG. 4) of the image density value (step S203). Here, in order to calculate the maximum value and the minimum value, they may be extracted from the smoothed image of the region from which the irradiation region and the void region are removed. Next, the smoothed image creating circuit 113c creates a smoothed image according to the following equation (step S204). fus (x, y) is a pixel value of a smoothed (low frequency) image of the input image f0 (x, y), and is represented by, for example, Expressions (14) to (18) or Expression (19).
[0055]

[0056]
“D (x, y)” is a disk-shaped filter, and is selected according to the input image, with r1 being an arbitrary constant.

[0057]
The profile of fus (x, y) obtained here preserves the edge structure and does not cause overshoot and undershoot, which are disadvantages of conventional dynamic range compression.
[0058]
For the smoothed image, for example, the average density represented by the equation (17) may be used, or the following equation (19) may be used. Here, d is a constant.
[0059]
[Expression 2]

[0060]
Further, for example, a smoothed image may be created using a morphological filter such as Erosion, Dai1ation, Opening, or C1osing.
[0061]
Next, the high frequency component creation circuit 113c creates a high frequency image fh (x, y) from the input image f0 (x, y) and the smoothed image fus (x, y) according to the following equation (step S205). Next, the density value conversion circuit 113e performs the next processing based on the density values d1 and d4 determined by the maximum value (d3), the minimum value (d2) extracted by the maximum / minimum value extraction circuit 113b, and the visual density value after gradation conversion. The density value is converted by the equation, and a smoothed image Sus0 × (x, y) having the density value converted is created (step S206).
[0062]
In FIG. 4, Lightmax and Lightmin are the maximum and minimum density values that can be visually recognized, and correspond to d1 and d4. When the density value of the smoothed image is less than d1 from min (minimum value that the smoothed image can take),

[0063]
If d2 <Sus (x, y) ≦ d3,

[0064]
If d3 <Sus (x, y) ≦ max,

[0065]
Then, the high-frequency image fh (x, y) is added to the density-converted image Susuo (x, y) to obtain the final image fe (x, y) (step S207).
fe (x, y) = Susu0 (x, y) + fh (x, y) (23)
Finally, the obtained image fe (x, y) is subjected to gradation conversion by the gradation conversion circuit 114 and displayed on the film output or the image display device 111.
[0066]
In this embodiment, the case of an image having a void area has been described. However, when there is no void area, the maximum and minimum values are extracted from the image in the irradiation area without passing through the void deletion circuit 113a. The maximum value and the minimum value may be extracted by the circuit 113b.
[0067]
According to the present embodiment, it is possible to compress and enlarge an image of a region of interest to a certain width while retaining a high frequency component, and there is an effect that image information can be used effectively. Furthermore, if the fixed width is set as a visual recognition area such as a film, there is an effect that the attention area can be enlarged or reduced to the width of the visual recognition area on the film. In addition, when the maximum value and the minimum value are extracted from the smoothed image, there is an effect that the conversion by the density value conversion unit 113 can be performed with higher accuracy.
[0068]
(Second Embodiment)
FIG. 5 shows a configuration of an X-ray imaging apparatus 600 to which the image processing apparatus according to the embodiment of the present invention is applied.
[0069]
The X-ray imaging apparatus 600 has an image processing function, and includes a preprocessing circuit 606, a CPU 608, a main memory 609, an operation panel 610, an image display 611, and an image processing circuit 612. Are exchanged with each other via a CPU bus 607.
[0070]
The X-ray imaging apparatus 600 includes a data acquisition circuit 605 connected to the preprocessing circuit 606, a two-dimensional X-ray sensor 604 and an X-ray generation circuit 601 connected to the data acquisition circuit 605. These circuits are also connected to the CPU bus 607.
[0071]
The image processing circuit 612 calculates a range and change amount for changing the dynamic range from a feature extraction circuit 612a that calculates a feature amount for tone conversion and a tone conversion curve of a first tone conversion circuit 612d described later. The control circuit 612b, the DRC circuit 612c that changes the dynamic range of the original image from the range and amount of change of the dynamic range calculated by the control circuit 612b, and the original image whose dynamic range has been changed by the DRC circuit 612c. A first gradation conversion circuit 612d that performs gradation conversion.
[0072]
The DRC circuit 612c adds a high-frequency component of the original image to the second gradation conversion circuit 613 that performs gradation conversion of the original image, and the image that has undergone gradation conversion by the second gradation conversion circuit 613. And a high-frequency component adjustment circuit.
[0073]
FIG. 6 shows, for example, a region a from which a feature amount is extracted in an image (X-ray image) obtained when imaging is performed on the side surface of the cervical spine.
[0074]
FIG. 7 shows a gradation conversion curve of the first gradation conversion circuit 612d. In FIG. 7, the horizontal axis indicates the pixel value and the vertical axis indicates the density value. “Sl” indicates a pixel value corresponding to the lower limit density value of the visual recognition limit on the density, and “Sh” indicates a pixel value corresponding to the upper limit density value of the visual recognition limit on the density. “Smin” and “Smax” indicate the minimum pixel value and the maximum pixel value of the region of interest. “Sa” indicates a pixel value corresponding to the feature amount calculated by the feature extraction circuit 612a, and “dl” and “dh” respectively correspond to the lower limit density value and the upper limit density value of the range in which the dynamic range is not changed. The pixel value to be shown.
[0075]
FIG. 8 shows a gradation conversion curve of the second gradation conversion circuit 613. In FIG. 8, the horizontal axis indicates the input pixel value, and the vertical axis indicates the output pixel value. “Dl”, “dh”, and “Sa” indicate the same pixel values as those in FIG.
The gradation conversion curve a shown in FIG. 8 is composed only of linear components, and the gradation conversion curve b shown in FIG. 8 is composed of a gentle curve in which the differential value of the gradation conversion curve a is continuous. Yes.
[0076]
Here, in the X-ray imaging apparatus 600 as described above, the main memory 609 stores processing programs and various data necessary for processing for operation control of the entire apparatus by the CPU 608, and the like. A work memory as a work memory of the CPU 608 is included.
Therefore, the CPU 608 controls the operation of the entire apparatus according to the operation from the operation panel 610 using the main memory 609 and the information stored in the main memory 609. Thereby, the X-ray imaging apparatus 600 operates as follows.
[0077]
First, the X-ray generation circuit 601 emits an X-ray beam 602 to the subject 603. The X-ray beam 602 emitted from the X-ray generation circuit 601 passes through the subject 603 while being attenuated, reaches the two-dimensional X-ray sensor 604, and is output as an X-ray image by the two-dimensional X-ray sensor 604. The X-ray image output from the two-dimensional X-ray sensor 604 is, for example, a human body part image.
[0078]
The data acquisition circuit 605 converts the X-ray image output from the two-dimensional X-ray sensor 604 into an electrical signal and supplies it to the preprocessing circuit 606. The preprocessing circuit 606 performs preprocessing such as offset correction processing and gain correction processing on the signal (X-ray image signal) from the data acquisition circuit 605. The X-ray image signal preprocessed by the preprocessing circuit 606 is transferred as an input image to the main memory 609 and the image processing circuit 612 via the CPU bus 607 under the control of the CPU 608.
[0079]
The image processing circuit 612 operates as follows, for example, according to the flowchart shown in FIG.
[0080]
First, the feature extraction circuit 612a that receives the input image (original image subjected to preprocessing) f0 (x, y) supplied from the preprocessing circuit 606 via the CPU bus 607 is controlled by the CPU 608. Feature S for conversion _a Is calculated (step S701).
The feature amount extraction method here is different for each part of the subject in the input image f0 (x, y), and a plurality of methods have been proposed by the present applicant. For example, when the input image f0 (x, y) is an image obtained by photographing the cervical spine, the method described in Japanese Patent Application No. 10-272284 is used. This method is a method of calculating an average value of pixel values in the region indicated by “a” in FIG. 6 as the feature amount Sa.
[0081]
Next, a gradation conversion curve in the first gradation conversion circuit 612d is defined using the feature amount Sa obtained in the feature extraction circuit 612a (step S702).
For example, as shown in FIG. 7, a gradation conversion curve is defined such that the pixel value of the feature value Sa is changed to a density value of 1.0.
[0082]
Next, the control circuit 612b calculates a range and a change amount for changing the dynamic range of the input image f0 (x, y) from the gradation conversion curve defined in step S702 (step S703). An example of this calculation method will be described next.
[0083]
For example, a normal Japanese has a density value that can be visually recognized in an image in a range of a density value of 0.2 to 3.0.
Therefore, first, pixel values Sl and Sh that are to be converted to density values of 0.2 and 3.0 are calculated based on the gradation conversion curve defined in step S702.
Next, the width that does not change the dynamic range,
dl = Sa-hl (24)
dh = Sa + hh (25)
It calculates according to the formulas (24) and (25).
In Expressions (24) and (25), “hl” and “hh” are constants, and “dl” and “dh” are pixels serving as starting points for changing the dynamic range.
[0084]
Then, the amount of change Rl, Rh of the pixel value below “dl” and the pixel value above “dl”
Rl = (Sl−dl) / (Smin−dl) (26)
Rh = (Sh-dh) / (Smax-dh) (27)
It calculates according to the following formulas (26) and (27).
In Expressions (26) and (27), “Smin” and “Smax” indicate the minimum value and the maximum value among the pixel values of the subject area (target area). That is, when the dynamic range of pixel values less than or equal to dl is multiplied by R1, the pixel value Smin is changed to the pixel value S1. Similarly, the pixel value Smax is changed to the pixel value Sh. In other words, assuming that the starting points are dl and dh and the change amounts are Rl and Rh, the density value of the subject area after gradation conversion is that of the visual field area (for example, density value 0.2 to 3.0 here). Matches the width.
[0085]
As described above, by obtaining the change amounts Rl and Rh, even if the dynamic range is compressed, all of the attention area of the captured image can be reproduced on the output image. Therefore, the instability of feature extraction can be absorbed, and the diagnostic ability can be improved.
[0086]
In addition, as a constraint condition for the change amounts Rl and Rh,
Rl <Cl (28)
Rh <Ch (29)
Conditions as shown in the following equations (28) and (29) may be provided.
In the expressions (28) and (29), “Cl” and “Ch” each represent a constant. These constants Cl and Ch and the above-described constants hl and hh are constants determined experimentally for each part of the subject to be imaged, for example.
In this way, by providing the constraint condition, it is possible to prevent excessive compression, and thus it is possible to reliably prevent a deterioration in diagnostic ability due to excessive compression.
[0087]
As described above, when the pixels dl and dh as starting points for changing the dynamic range and the change amounts Rl and Rh are obtained, the second values as shown in FIG. 8 are obtained based on these dl, dh, Rl, and Rh. The gradation conversion curve of the gradation conversion circuit 613 is defined.
In the gradation conversion curve a shown in FIG. 8, the slope below the pixel value dl is “Rl”, the slope from the pixel value dl to the pixel value dh is “1”, and the slope above the pixel value dh is “ R _h On the other hand, the gradation conversion curve b shown in FIG. 8 is also smoothed so that the differential value of the gradation conversion curve a becomes continuous by the technique described in Japanese Patent Application No. 11-762882. When the gradation conversion curve a is used as the gradation conversion curve of the second gradation conversion circuit 613, a false contour may occur at the differential discontinuity point. In the case where b is used, the false contour does not occur, and the gradation conversion curve b is used in the second gradation conversion circuit 613 as an example here.
[0088]
Next, the DRC circuit 612c selects whether or not to retain the high frequency component in the subsequent processing (step S704).
The result image obtained when the processing is executed without retaining the high frequency component is effective for diagnosing soft tissue or the like where the low frequency component is important. On the other hand, a result image obtained when processing is executed while holding a high-frequency component is effective when diagnosing a bone part, a lung field, or the like where the high-frequency component is important. That is, if the image to be processed is an image of a soft tissue or the like, select that the high frequency component is not retained, and if the image to be processed is an image of a bone or lung field, select to retain the high frequency component. As a result, a result image effective for diagnosis of soft tissue, bone, lung field, and the like can be obtained.
The selection in step S704 may be configured to be automatically selected based on the imaging region (soft tissue, bone, lung field, etc.), or may be configured to be manually selected by the user.
[0089]
If the high-frequency component is not retained as a result of the selection in step S704, the second gradation conversion circuit 613 performs gradation conversion using the gradation conversion curve b (see FIG. 8) for the input image (original image). (Step S708).
Then, the first gradation conversion circuit 612d performs gradation conversion on the input image after the gradation conversion by the second gradation conversion circuit 613 (step S709).
[0090]
Here, the coordinates of the input image (original image) are “x, y”, the pixel value of the original image is “f0 (x, y)”, and the gradation conversion curve b of the second gradation conversion circuit 613 is “F”. (X) ”, the pixel value of the image after gradation conversion using the gradation conversion curve b is“ f1 (x, y) ”, and the gradation conversion curve of the first gradation conversion circuit 612d is“ F1 (x ) ", The resulting image (final image) f2 (x, y) is
f2 (x, y) = F1 (F (f0 (x, y)) (30)
It is represented by the following formula (30).
In this case, “F1 (F (x)) may be considered as one gradation conversion curve.
[0091]
On the other hand, when the high frequency component is held as a result of the selection in step S704, first, the high frequency component adjustment circuit 614 converts the high frequency component fh (x, y) of the original image f0 (x, y) to
fh (x, y) = f0 (x, y) −fus (x, y) (31)
This is calculated by the following equation (31) (step S705).
In Expression (31), “fus (x, y)” is a smoothed image, and has a constant d indicating the size of the mask size.
[0092]
[Equation 3]

[0093]
It is calculated by the following equation (32).
[0094]
Next, the second gradation conversion circuit 613 applies the original image f0 (x, y) to
f1 (x, y) = F (f0 (x, y)) (33)
Gradation conversion represented by the following equation (33) is performed, and the image f1 (x, y) after the gradation conversion is acquired (step S706).
[0095]
Next, the high frequency component adjustment circuit 614 calculates the coefficient c (x) from the derivative of the gradation conversion curve b of the second gradation conversion circuit 613.
[0096]
[Expression 4]

[0097]
It calculates with the type | formula (34) which becomes. As shown in the equation (34), the coefficient c (x) is a value obtained by subtracting “1” from the gradient of the gradation conversion curve b.
Then, based on the coefficient c (x), the image f1 (x, y) after the gradation conversion by the second gradation conversion circuit 613 acquired in step S706 is acquired in step S706. The high-frequency component fh (x, y) of the original image f0 (x, y) is added to obtain a processed image f3 (x, y). That is, the processed image f3 (x, y) is

This is expressed by the following equation (35) (step S707).
[0098]
Then, with respect to the processed image f3 (x, y) obtained in step S707 by the first gradation conversion circuit 612d,
f2 (x, y) = F1 (f3 (x, y)) (36)
Gradation conversion according to the following equation (36) is performed, and an image f2 (x, y) obtained as a result is obtained as a final image (step S709).
[0099]
Note that the smoothed image fus (x, y) in the equation (31) in step S705 may be calculated using, for example, the morphological operation represented by the following equations (37) to (41).

In Expressions (37) to (40), “D (x1, y1)” is a disk-shaped filter, and has an arbitrary constant r1 selected according to the input image,

It is represented by the following formula (41).
[0100]
The profile of the smoothed image fus (x, y) obtained by the equations (37) to (41) preserves the edge structure, and overshoot and undershoot, which are disadvantages of conventional dynamic range compression, occur. There is nothing.
[0101]
According to the present embodiment, the density distribution width of an arbitrary gradation region of the input image can be compressed and expanded, and the amplitude of the high-frequency component after gradation conversion is set to the high-frequency component of the image before gradation conversion. There is an effect that the amplitude can be maintained.
Also, since the dynamic range of the image is changed based on the pixel value determined by the gradation conversion curve, the dynamic range can be changed assuming the density value of the image after gradation conversion. The range can be adjusted uniquely.
[0102]
In addition, since the pixel value determined by the visibility limit is obtained from the gradation conversion curve, and the range and amount for changing the dynamic range can be determined based on the pixel value, the density value range of the pixel value in the region of interest can be set as the viewing range. There is an effect that can be done.
In addition, since the attention area is in the visual recognition area, it is possible to observe the entire attention area, which has the effect of improving the diagnostic ability.
Further, the attention area can be matched with the width of the visual recognition area, and the information amount of the attention area can be expressed to the maximum as the density value.
Also, since it has a range that does not change the dynamic range of the image, the pixel range that is important for diagnosis etc. is set as a region where the dynamic range is not changed, and it is expressed by the conventional density value, and there are also regions that could not be visually recognized in the past Since it can be observed as a concentration value, it has the effect of improving diagnostic ability.
In addition, since the high frequency component is retained in the dynamic range changing means, the amount of information of the high frequency component is not reduced, and the density value of the region that could not be visually recognized on the density value can be observed conventionally. effective.
[0103]
(Other embodiments)
The above-described embodiment may be realized in hardware or may be configured in software. For example, the configuration in the above-described embodiment is realized by realizing a program for executing the processing according to the flowchart shown in FIG. 2 or FIG. 9 using a computer system including a memory such as a CPU and a main memory. May be.
[0104]
As a storage medium for storing the program, a semiconductor memory such as a ROM or a RAM, an optical disk, a magneto-optical disk, a magnetic storage medium, or the like may be used. These may be a CD-ROM, an FD, a magnetic card, a magnetic tape, or a nonvolatile memory. A card etc. can be used.
[0105]
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 the CPU or the like provided in the extension function board or extension 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 the above embodiment. It is needless to say that these functions can be realized, and equivalent effects can be obtained, which are included in the present invention.
[0106]
【The invention's effect】
As described above, according to the present invention, it is possible to compress and expand the image of the region of interest to a certain width while retaining the high-frequency component that is a fine part of the image. There is an effect that can be used effectively. Further, if the constant width is set as a visual recognition region such as a film, the attention region can be enlarged or reduced to the width of the visual recognition region on the film.
[0107]
In addition, when there is a void region in the image, it is possible to extract the region of interest more accurately by removing the void region, which has the effect of increasing the accuracy of density value conversion.
[0108]
Furthermore, there is an effect that density value conversion can be performed more accurately by extracting the maximum value and the minimum value from the smoothed image.
[Brief description of the drawings]
FIG. 1 is a block diagram illustrating a configuration of an image processing apparatus according to a first embodiment.
FIG. 2 is a flowchart showing a processing procedure in the first embodiment.
FIG. 3 is a configuration diagram illustrating an image from which an input image, an irradiation region, and a void region are removed.
FIG. 4 is a characteristic diagram for explaining an operation in the first embodiment;
FIG. 5 is a block diagram showing a configuration of an X-ray imaging apparatus to which an image processing apparatus according to a second embodiment is applied.
FIG. 6 is a diagram for explaining a region from which a feature amount is extracted.
FIG. 7 is a diagram for explaining a gradation conversion curve of a first gradation conversion curve of the image processing apparatus.
FIG. 8 is a diagram for explaining a gradation conversion curve of a second gradation conversion curve of the image processing apparatus.
FIG. 9 is a flowchart showing a processing procedure in the second embodiment;
[Explanation of symbols]
105 Data collection circuit
108 CPU
109 Main memory
110 Operation panel
112 Irradiation area extraction circuit
113 Density value conversion circuit
113a Through hole deletion circuit
113b Maximum and minimum value extraction circuit
113c Smoothed image creation circuit
113d High frequency component creation circuit
113e Density value conversion circuit
112f High frequency component addition circuit

Claims

An image processing method for performing gradation conversion processing and frequency processing on original image data,
Extracting a value for generating a gradation conversion curve from a region in the subject in the original image data;
Generating a high frequency component based on the original image data;
Changing the high-frequency component based on the slope of the gradation conversion curve;
And a step of generating image data based on the changed high-frequency wave component.

2. The image according to claim 1, wherein the value for generating the gradation conversion curve is at least one of a maximum value, a minimum value, and a feature amount extracted from a region in the subject. Processing method.

The image processing method according to claim 1, further comprising a step of performing gradation conversion processing on the original image data based on the gradation conversion curve.

The region within the subject is extracted based on a hollow region that is a region directly irradiated with X-rays and a region that contacts the hollow region with a certain width. 2. An image processing method according to item 1.

A computer-readable storage medium storing a program for causing a computer to execute an image processing method for performing gradation conversion processing and frequency processing on original image data,
Extracting a value for generating a gradation conversion curve from a region in the subject in the original image data;
Generating a high frequency component based on the original image data;
Changing the high-frequency component based on the slope of the gradation conversion curve;
A computer-readable storage medium storing a program for causing a computer to execute the step of generating image data based on the changed high-frequency wave component.