JP2013102951A - 撮像装置および画像処理方法 - Google Patents

撮像装置および画像処理方法 Download PDF

Info

Publication number
JP2013102951A
JP2013102951A JP2011248717A JP2011248717A JP2013102951A JP 2013102951 A JP2013102951 A JP 2013102951A JP 2011248717 A JP2011248717 A JP 2011248717A JP 2011248717 A JP2011248717 A JP 2011248717A JP 2013102951 A JP2013102951 A JP 2013102951A
Authority
JP
Japan
Prior art keywords
phase
data
differential phase
image
phase data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2011248717A
Other languages
English (en)
Other versions
JP5868132B2 (ja
Inventor
Kentaro Nagai
健太郎 長井
Mitsuo Takeda
光夫 武田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Canon Inc
University of Electro Communications NUC
Original Assignee
Canon Inc
University of Electro Communications NUC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Canon Inc, University of Electro Communications NUC filed Critical Canon Inc
Priority to JP2011248717A priority Critical patent/JP5868132B2/ja
Priority to US14/347,270 priority patent/US9019479B2/en
Priority to EP12798431.8A priority patent/EP2779902A1/en
Priority to PCT/JP2012/079051 priority patent/WO2013073453A1/en
Publication of JP2013102951A publication Critical patent/JP2013102951A/ja
Application granted granted Critical
Publication of JP5868132B2 publication Critical patent/JP5868132B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/041Phase-contrast imaging, e.g. using grating interferometers
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/48Diagnostic techniques
    • A61B6/484Diagnostic techniques involving phase contrast X-ray imaging
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01JMEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J9/00Measuring optical phase difference; Determining degree of coherence; Measuring optical wavelength
    • G01J9/02Measuring optical phase difference; Determining degree of coherence; Measuring optical wavelength by interferometric methods
    • G01J9/0215Measuring optical phase difference; Determining degree of coherence; Measuring optical wavelength by interferometric methods by shearing interferometric methods
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/30Accessories, mechanical or electrical features
    • G01N2223/345Accessories, mechanical or electrical features mathematical transformations on beams or signals, e.g. Fourier
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/40Imaging
    • G01N2223/401Imaging image processing

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Pathology (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Animal Behavior & Ethology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Surgery (AREA)
  • Optics & Photonics (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)
  • General Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Immunology (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Biochemistry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

【課題】位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成可能な技術を提供する。
【解決手段】干渉パターンのデータからx方向とy方向の微分位相データを抽出し、抽出されたx方向の微分位相データをx方向に関して微分しx方向の二階微分位相データを算出し、y方向の微分位相データをy方向に関して微分しy方向の二階微分位相データを算出する。そして、x方向とy方向の二階微分位相データを関数として含む二階微分方程式を解くことによって、被検体の位相像のデータを算出する。
【選択図】図6

Description

本発明はシアリング干渉計により得られる干渉パターンから位相像や微分位相像のデータを生成するための技術に関する。
シアリング干渉計とは、同一光源から照射された光線などのコヒーレントな光を分割し、一方に被検体による波面の歪曲を生成し、もう一方をわずかにずらすことで干渉縞を形成し、この干渉縞のずれから被検体による光の位相変化を可視化する手段である。このようなシアリング干渉計は、例えば位相差顕微鏡などの分野で被検体の視認性を向上させる目的などに利用されている。
シアリング干渉計の一つとして、トールボット干渉法を利用したものが知られている。トールボット干渉法では、光源からの光を被検体に照射し、被検体を透過することで位相が変化した光を回折格子に導く。回折格子により回折された光は、回折格子からトールボット距離と呼ばれる所定の距離だけ離れた位置に干渉パターンを形成する。X線を用いた場合、この干渉パターンの周期は非常に小さいため、干渉パターンを直接検出することが難しいことがある。そこで、干渉パターンが形成される位置に吸収格子を配置し、吸収格子によって干渉パターンの一部を遮ることでモアレを形成し、このモアレ像を検出器で検出する方法が提案されている。このモアレ像を解析することにより、被検体の位相の情報を得ることができる。このように光の位相を利用したイメージングは一般に位相イメージング(X線を用いる場合はX線位相イメージング)と呼ばれている。
ただし、干渉パターンによって直接解析できるのは、あくまで被検体を透過することによって変化した光の位相の空間的な微分量(微分位相像、シア像とも呼ばれる)である。空間的な微分はエッジを強調する効果があるため、微分位相像は光学顕微鏡などの分野では被検体像を強調するために用いられる。また、この微分位相像から元の位相を計算し、位相像を得ることも可能である。
非特許文献1には、二次元トールボット干渉法を用いたX線位相イメージング手法が開示されている。具体的には、二次元の回折格子により形成された干渉パターンからX方向の微分位相像とY方向の微分位相像を取得し、それらに対しフーリエ積分を行うことによって、位相像を計算している。
C. Kottler, C. David, F. Pfeiffer and O. Bunk, "A two-directional approach for grating based differential phase contrast imaging using hard x-rays," Opt. Express 15 (2007) 1175-1181.
シアリング干渉計を利用した位相イメージングは、画像診断や各種の非破壊検査への応用が可能である。中でもX線位相イメージングは、CT(コンピュータ・トモグラフィ)などX線吸収量を画像化する従来の手法に比べて被ばく量を大幅に低減できる、筋肉や内蔵などの軟組織も画像化できるといった利点から、新たな画像診断手法として期待されている。
この種の位相イメージングにおいては、ノイズやアーチファクトの低減による画質向上が基本的な課題として存在する。特に、光源としてX線を用いる場合には、量子ノイズが不可避であり、しかも被ばく量低減のために線量を下げるほどノイズが増える傾向にあるため、ノイズ対策は重要となる。また、X線ではなく可視光を用いる装置であっても、検査対象によっては(例えば眼球)低い光量しか用いることができない場合があり、SN比の悪化が問題となる。
また、位相イメージングにおいては「位相折りたたみ(位相ラッピング)」という問題があることが知られている。位相折りたたみは、干渉パターンから得られる微分位相の値が2πの範囲に折りたたまれて(ラップされて)しまう現象である。例えば−πからπの間で折りたたまれる場合、実際には3/2πという位相が正しい値であったとしても、それは−1/2πとして測定される。非特許文献1においてはこのような位相折りたたみの問題は無視されている。そのため、非特許文献1の手法で生成した位相像には、本来何もない、すなわち位相が一定になる空間に、位相のゆがみやアーチファクトが出現する場合がある。
位相折りたたみへの対策として、「位相接続(位相のアンラッピング)」と呼ばれる処理が知られている。位相接続とは、位相が折りたたまれた箇所を検出し、位相を元に戻す操作である。位相接続に関しては様々な手法が提案されている。最も基本的な操作は隣接ピクセルとの位相の差を計算する手法である。その差の絶対値が一定以上(例えば3/2π以上)であるならば、そこはもともとスムーズにつながっていた位相が折りたたまれたものと疑うことができ、位相を元に戻す操作をする。具体的には位相飛びが生じているピクセル以降のピクセル全ての位相に全て2πか、−2πを足すことで位相隣接ピクセルとの位相差を3/2π以下にする操作を行う。従来の位相接続は、積分とは独立した操作として、微分位相像全体に対して適用する必要があるため、位相像を得るまでに要する処理時間が長くなるという問題があった。
本発明は上記実情に鑑みてなされたものであり、位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成可能な技術を提供することを目的とする。
本発明のさらなる目的は、位相イメージングにおいて、位相像あるいは微分位相像におけるノイズの低減、並びに、位相折りたたみによる位相のゆがみやアーチファクトの抑制を、簡易な処理で実現するための技術を提供することにある。
本発明の第1態様は、シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する撮像装置であって、被検体を透過した電磁波により形成された干渉パターンのデータから、第1の方向に関する位相の変化を表す第1の微分位相データと、前記第1の方向に交差する第2の方向に関する位相の変化を表す第2の微分位相データとを抽出する微分位相データ抽出手段と、抽出された前記第1の微分位相データを前記第1の方向に関して微分することによって第1の二階微分位相データを算出すると共に、抽出された前記第2の微分位相データを前記第2の方向に関して微分することによって第2の二階微分位相データを算出する二階微分位相データ算出手段と、前記第1および第2の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する位相データ算出手段と、を有する撮像装置である。
本発明の第2態様は、シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する画像処理方法であって、情報処理装置が、被検体を透過した電磁波により形成された干渉パターンのデータから、第1の方向に関する位相の変化を表す第1の
微分位相データと、前記第1の方向に交差する第2の方向に関する位相の変化を表す第2の微分位相データとを抽出する工程と、情報処理装置が、抽出された前記第1の微分位相データを前記第1の方向に関して微分することによって第1の二階微分位相データを算出すると共に、抽出された前記第2の微分位相データを前記第2の方向に関して微分することによって第2の二階微分位相データを算出する工程と、情報処理装置が、前記第1および第2の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する工程と、を有する画像処理方法である。
本発明の第3態様は、上記画像処理方法の各工程を情報処理装置に実行させるプログラムである。
本発明によれば、位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成することができる。また、位相イメージングにおいて、位相像あるいは微分位相像におけるノイズの低減、並びに、位相折りたたみによる位相のゆがみやアーチファクトの抑制を、簡易な処理で実現することも可能となる。
撮像装置の構成を示す図である。 回折格子、干渉パターン、吸収格子の構成例を示す図である。 情報処理装置の機能構成を示す図である。 干渉パターンから微分位相像を抽出する手法を示す図である。 従来例における微分位相像から位相像を取得する処理を示す図である。 実施例における微分位相像から位相像を取得する処理を示す図である。 シミュレーションで用いた被検体と微分位相像を示す図である。 従来例の手法で得られた位相像を示す図である。 実施例の手法で得られた位相像を示す図である。 実施例の手法で得られた位相像から生成した微分位相像を示す図である。
以下、本発明の好ましい実施の形態について添付の図面に基づいて説明する。なお、各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。
本実施形態では微分位相像を得る手段としてX線トールボット干渉計を使用した例について説明する。ただし、本発明が適用可能な微分位相像の取得手法はこれに限られるわけではなく、被検体(被検物)の微分位相像を得ることができれば、如何なる形態のシアリング干渉計を用いてもよい。また測定に用いる光はX線に限らず、如何なる波長の電磁波を用いてもよい。なお本明細書においてX線とはエネルギーが2以上200keV以下の光を指す。
(撮像装置の構成)
図1は本実施形態における撮像装置の構成を示した図である。図1に示した撮像装置1は、概略、干渉計10と情報処理装置11とから構成される。干渉計10は、2次元X線トールボット干渉計であり、X線を発生させるX線源(光源)110、X線を回折する回折格子310、X線の一部を遮る吸収格子410、X線を検出する検出器510を備えている。情報処理装置11は、演算部(コンピュータ本体)610と画像表示装置710を備えている。情報処理装置11は干渉計10にケーブルで接続されており、干渉計10の各部の制御を行う。また干渉計10の検出器510で得られた干渉パターンのデータは情報処理装置11に伝送され、演算部610において後述する各種の演算処理に供される。演算部610による演算結果に基づいた画像は画像表示装置710に出力される。
測定対象である被検体210は、X線源110と回折格子310の間に配置される。X線源110から照射され被検体210を透過したX線は回折格子310により回折され、トールボット距離と呼ばれる所定の距離Lの面上に明部と暗部が配列方向に並んだ干渉パターン810を形成する。本明細書では、X線(光)の強度が大きい所を明部、小さい所を暗部とする。本実施形態に用いられる回折格子310は位相型の回折格子である。回折格子として振幅型の回折格子を用いることもできるが、位相型の回折格子の方がX線量(光量)の損失が少ないので有利である。
図2(a)は位相型の回折格子310(以下単に位相格子310と呼ぶ)の平面上の構成の一例を表した図である。311は位相の基準となる領域、312は基準の領域311に対してπほどX線の位相が変化する領域である。図2(b)は位相格子310によってトールボット距離に生じる干渉パターン810の明部と暗部を表している。これらの組み合わせには何種類かあるが、本発明では公知のものを含めて如何なる構成のものを用いてもよいため、ここでは詳しい説明を割愛する。
この干渉パターン810から微分位相像を抽出する方法には大きく分けて二種類考えられる。一つは、吸収格子410を少しずつずらしながら複数回撮像し、その結果から微分位相像を計算する縞走査法である。もう一つは、自己像としての干渉パターン810の周期と少しずらした吸収格子410を利用して撮像する方法で、フーリエ変換法と呼ばれる。以下ではフーリエ変換法に関して説明するが、どちらの手法でも微分位相像を取得できる。なお吸収格子は、遮蔽格子または振幅格子とも呼ばれることもある。
フーリエ変換法では図2(c)で示すような吸収格子410を用いる。吸収格子410には、干渉パターン810の明部に対応する透過部411と、暗部に対応する遮蔽部412とが設けられている。ただし、透過部411の周期(ピッチ)は干渉パターン810の明部の周期とは異なる。このため、吸収格子410を透過した光像を実際に撮像すると、干渉パターン810と吸収格子410の周期の違いによるモアレが観測される。
通常、位相格子310による干渉パターン810の周期は数μmから十数μm程度であるが、吸収格子410によりモアレを作ることにより、その周期を数十μm以上に拡大することができる。これによって数十μm程度の分解能の検出器510で干渉パターンを撮像することができる。しかしながら、例えば、検出器510の分解能が十分高い場合や、位相格子310と検出器510の距離が大きい場合など、検出器510で干渉パターン810そのものを解像可能な場合には、吸収格子410を設けなくてもよい。
(情報処理装置の機能)
図3は、情報処理装置11が有する機能を示すブロック図である。情報処理装置11は、概略、干渉計10の制御(例えばX線の照射、検出器510からのデータ取り込みなど)を行う制御部30と、干渉計10から取り込まれたデータに基づいて位相像のデータを生成するデータ処理部31とを有する。データ処理部31は、微分位相データ抽出部32、二階微分位相データ算出部33、及び、位相データ算出部34から構成されている。微分位相データ抽出部32は、干渉パターン(モアレ像)のデータを数値的に解析して微分位相データを抽出する機能であり、二階微分位相データ算出部33は、微分位相データを微分することによって二階微分位相データを算出する機能である。位相データ算出部34は、二階微分位相データを関数として含む二階微分方程式を解くことによって、被検体の位相像のデータを算出する機能である。これらの機能の詳細については後述する。
情報処理装置11の演算部610は、CPU(中央演算処理装置)、主記憶装置(RAMなど)、補助記憶装置(HDD、SSDなど)、各種I/Fを有するコンピュータで構
成することができる。図3に示した機能は、補助記憶装置に格納されたプログラムが主記憶装置にロードされ、CPUにより実行されることで実現されるものである。もちろん、この構成はあくまで一例であり、本発明の範囲を限定するものではない。
(撮像装置の動作)
以下、本実施形態の撮像装置1において実行される画像処理方法、並びに従来例との差異を説明する。ただしここでは、干渉パターンのデータとして、干渉計10により実際に被検物を計測して得たデータではなく、コンピュータシミュレーションによって生成したデータを用いる。干渉パターンのデータから位相像のデータを生成する一連の処理については、情報処理装置11のデータ処理部31によって行った。
シミュレーションでは、X線として17.5keVの単色光を仮定し、位相格子として図2(a)に示す構造の位相格子310を仮定した。位相格子310には、X線を十分透過させる基盤部311に対して柱状の構造物312がピッチp1で並んでいると仮定する。構造物312を透過したX線は、基盤部311を透過したX線に対して、位相がπ変化するようにした。このような構造の位相格子310は例えば構造物312がシリコンで作成された場合、22.6μm程度の奥行きにすることによって実現される。
さらに図2(a)においてピッチp1を8μmと仮定すると図1における位相格子310と吸収格子410の距離Lを約11.3cmに設定することにより、吸収格子410の直前で干渉パターン810が形成される。この干渉パターン810の周期p2は4μmである。吸収格子410によりモアレを形成することで、干渉パターン810を拡大し、検出器510でその強度を検出する。簡単のため、検出器510に含まれる画素素子の大きさは4μmとして計算した。
トールボットのような干渉計の場合、入射光には高いコヒーレンス性が必要とされる。そのような高いコヒーレンス性を実験室や医療装置のような大きさがせいぜい数m程度までの施設で用いる場合、微小な点光源を用いてコヒーレンス性を出す手法が用いられる。このような場合、入射光は平行波ではなく球面波になる。しかし、この違いは本発明の本質を変えるものではないため、本発明は平行光でも非平行光(球面波など)でも適用可能である。以下の実施例では平行光を仮定してシミュレーションを行った。
被検体としては二種類用意した。一つ目は図7(a)に示すような直径400μmのカルシウム球を想定した。もう一つは図7(b)に示すような直径200μmのカルシウム球が4つ重なっている被検体である。計算の範囲は512μm×512μmの範囲となる。
微分位相像の抽出方法としては前述の通りフーリエ変換法を使用した。そのために、吸収格子410の透過部411の周期p3(図2(c)参照)は4.5μmとした。これによりほぼ4画素で一周期のモアレが実現することとなる。
図4(a)から図4(d)はフーリエ変換法により干渉パターンから微分位相像を抽出する方法、すなわち微分位相データ抽出部32の処理内容を概念的に示した図である。図2(b)のような干渉パターン810を拡大した結果、図4(a)のようなモアレ像が取得できたとする。フーリエ変換法ではこのモアレ像をフーリエ変換する。これによってフーリエ空間内では9つのスペクトルに分離される。この9つのスペクトルのうち1次スペクトルにあたるスペクトル911、912、913、914のいずれかの周辺を窓関数で切り取って原点に戻し(図4(c))、これを逆フーリエ変換することで位相微分像を得ることができる。図4(d)の左側は、スペクトル914から生成された微分位相像の例であり、x方向(第1の方向)に関する位相の変化を表している。図4(d)の右側は、
スペクトル911から生成された微分位相像の例であり、y方向(第2の方向)に関する位相の変化を表している。
図7(c)は、図7(a)の被検体から得られたx方向、y方向それぞれの微分位相像を示す。また図7(d)は、図7(b)の被検体から得られたx方向、y方向それぞれの微分位相像を示す。球体のエッジ部分は最も急激に位相の微分が変化する個所となるため、この部分で位相飛びが発生しやすい。また、フーリエ変換法では、前述の窓関数の取り方によっては微分位相像にアーチファクト(本来は存在しない像)が発生する場合がある。図7(c)、および図7(d)の微分位相像でも若干のアーチファクトが見られる。
二次元の場合はx方向とy方向の微分位相像がそれぞれ独立に取得できる。この方向の違いは、図4(b)のフーリエ空間においてどのスペクトルを切り取るかによって制御できる。例えばスペクトル911又は912を切り取ればy方向の微分位相像が得られ、スペクトル913又は914を切り取ればx方向の微分位相像が得られる。もちろんx、yの直交する2方向ではなく、互いに交差する任意の2方向についての微分位相像を得ることも可能である。以下では取得されたx方向微分位相像をφx(x,y)、y方向微分位相像をφy(x,y)とおく。(x,y)は座標である。これらのx方向微分位相像とy方向微分位相像から、積分された位相像を求める。
(1)従来の位相像取得方法
従来例では、x方向微分位相像とy方向微分位相像を別々に積分するか、フーリエ変換などの手段を用いて積分する。以下ではより誤差の少ない位相が求められるフーリエ変換による積分法(フーリエ積分法)を従来例として述べる。
図5はフーリエ積分法による位相像取得方法の手順を示すフローチャートである。まず前述のように、x方向微分位相像φx(x,y)とy方向微分位相像φy(x,y)を取得する(ステップS50)。g(x,y)を次の式で定義する(ステップS51)。
Figure 2013102951
ここで、iは虚数単位を示す。
g(x,y)をフーリエ変換することでG(kx,ky)を得て(ステップS52)、次の(式3)によって位相P(x,y)を取得する(ステップS53、S54)。
Figure 2013102951
ここで、F[・・・]は括弧の中をフーリエ変換することを示す演算子であり、F−1[・・・]は逆フーリエ変換を示す演算子である。kxとkyはそれぞれx方向、y方向の波数を表す。(式3)はすなわち(式1)で定義されたgをフーリエ変換し積分演算子1/(2πi×(kx+ky))を作用させて逆フーリエ変換をかけることと一致する。こうして得たP(x,y)は被検体の位相像を表す(ステップS55)。
図8(a)は、図7(c)で示した微分位相像より、図5に示した従来例による方法で取得した位相像の例である。カルシウム球の位相が球形に回復されているのがわかる。しかし、位相の折りたたみにより被検体の外周部で位相が不連続になっており、取得した位相像にゆがみが生じている。特に何もない真空部(一定)であるはずの被検体の外周部に
、本来ありえない位相のゆがみが生じていることがわかる。
図8(b)は、図7(d)で示した微分位相像より、図5に示した従来例による方法で取得した位相像の図である。4つのカルシウム球の位相が復元されていることがわかる。しかし、先の例と同様に本来は真空であるはずの領域に位相のゆがみができている。また、アーチファクトの影響で本来は球形である位相がジャギーを持っているように見える。
(2)実施例の位相像取得方法
図6は、本発明の実施例における位相像取得方法の手順を示すフローチャートである。従来例では、微分位相像を直接積分することで位相像を得ていたのに対し、本実施例では、微分位相像をx,yそれぞれの方向について微分し二階微分位相像を生成した後、二階微分位相像からなる二階微分方程式を解いて位相像を得る点に特徴がある。
まず、微分位相データ抽出部32が、干渉パターンのデータ(モアレ像のデータ)から、x方向微分位相像φx(x,y)とy方向微分位相像φy(x,y)を取得する(ステップS60)。微分位相データ抽出部32が生成したx方向微分位相像のデータ(第1の微分位相データ)およびy方向微分位相像のデータ(第2の微分位相データ)は、二階微分位相データ算出部33に引き渡される。
二階微分位相データ算出部33は、x方向微分位相像をx方向に関して微分することによってx方向二階微分位相像を算出すると共に、y方向微分位相像をy方向に関して微分することによってy方向二階微分位相像を算出する(ステップS61)。x方向二階微分位相像ψx(x,y)、y方向二階微分位相像ψy(x,y)は、下記式で表される。
Figure 2013102951
x方向二階微分位相像のデータ(第1の二階微分位相データ)およびy方向二階微分位相像のデータ(第2の二階微分位相データ)は、位相データ算出部34に引き渡される。
位相データ算出部34は、既知関数としてx方向、y方向それぞれの二階微分位相像を含む二階微分方程式を解くことで、未知関数であるx方向、y方向それぞれの位相像を求める(ステップS62〜S66)。そのアルゴリズムについて以下詳しく述べる。
まず位相データ算出部34は、x方向とy方向の二次微分位相像を用いて、ρ(x,y)を次のように定義する(ステップS62)。
Figure 2013102951
ここで、W[・・・]は、[・・・]の値のうち、所定の範囲を外れた値をその所定の範囲内に収める処理を行う演算子である。微分位相像φx、φyの値は基本的に連続的に変化するため、二階微分位相像ψx、ψyの値は0を中心とする狭い範囲(せいぜい−π〜πの範囲)の内に収まる。しかし、位相の折りたたみが生じると、微分位相像φx、φyの値が−πからπ、又はπから−πへと不連続に変化するため、二階微分位相像ψx、ψyの値が2π又は−2πとなる。すなわち、位相の折りたたみが生じた箇所は、二階微分位相像ψx、ψyにおいてスパイク状の異常値となって現れる。W[・・・]の操作は
、この異常値を除去(修正)する操作、すなわち、位相の折りたたみが生じた状態を元の位相の状態に戻す操作(位相接続)に相当する。以下、W[・・・]をアンラッピング演算子と呼ぶ。
アンラッピング演算子にはさまざまなデザインが考えられるが、本実施例では、具体的かつもっとも簡単な演算子として、次のルールを持つ演算子W[・・・」を用いる。
(a)[・・・]の値がπ以上のものに関して、その値から2πを引く
(b)[・・・]の値が−π以下のものに関して、その値に2πを足す
以上のルールを適用することで元関数φx(x,y)やφy(x,y)に位相折りたたみが生じていても、位相接続がなされたのと同じ効果を持つ。
実際の位相接続では、上記のみならずノイズによる特異点の処理などを含めた様々な手法が存在する。W[・・・]もそれらの処理方法を応用した様々な演算アルゴリズムが考えられるが、本発明はそれらの手法に本質的には依存せず様々な方法を適用可能である。
位相データ算出部34は、(式6)のρ(x,y)から次のポワソン方程式を解くことによって位相P(x,y)を計算する。
Figure 2013102951
このポワソン方程式を解く具体的アルゴリズムも様々な手法が考えられるがその一つとしては、フーリエ変換による二階の積分法が考えられる。すなわち、位相データ算出部34は、ρ(x,y)をフーリエ変換することでH(kx,ky)を得て(ステップS63)、次の(式9)によって位相P(x,y)を算出する(ステップS64、S65)。
Figure 2013102951
(式9)では(式6)で定義されたρ(x,y)に対してフーリエ変換を求め、1/(4π(kx+ky))の二階の積分演算子を作用させた後に逆フーリエ変換して、取得した値を位相P(x,y)とする。以上の処理により、二階微分位相像の二階積分である位相像を得ることができる。
図9(a)は、図7(c)で示した微分位相像より、図6に示した本実施例の手法によって取得した位相像の例である。従来例の図8(a)と比べると、位相折りたたみによる歪みが改善されていることがわかる。これは従来例では位相が不連続になっていた箇所が、ある程度修復された影響によるものであり、本実施例の手法の有効性を示している。
図9(b)は、図7(d)で示した微分位相像より、図6に示した本実施例の手法によって取得した位相像の例である。先の例と同様に、従来例の図8(b)と比べると位相折りたたみによる歪みが改善されており、本実施例の手法の有効性を示している。また、従来例の図8(b)では目立っているアーチファクトについても、図9(b)では低減していることがわかる。
また本実施例の手法はノイズ耐性に関しても優れている。表1は、従来例の位相像取得方法と本実施例の位相像取得方法とで、同じノイズに対するノイズリダクションの機能を
比べた結果である。入力画像として平均0、分散が0.5のランダムパターンの画像2枚を微分位相像とみなして位相を計算し、その位相像にどれだけのノイズが存在するかを比較した表である。この実験では、純粋にノイズリダクションを評価するため、被検体を含まない画像を用いている。表1の数値は分散であり、値が小さいほどノイズが少ない(ノイズが低減されている)ことを示す。10回の計算を行なった結果、すべての回において本実施例の手法でノイズが減少していた。
Figure 2013102951
このように本実施例の手法は従来の手法よりもノイズ低減の効果もあることが分かる。この実験では被検体を含まない画像を用いたので、このノイズ低減の効果は、前述したアンラッピング演算子Wによる位相接続による画像の改善とは、本質的に別の効果である。すなわち位相接続を行わなくとも(つまり(式6)において演算子Wを省略しても)、本実施例の手法が従来の手法に比べて有利な効果を有していることがわかる。
図10(a)、図10(b)はノイズ低減効果の応用例を示している。図10(a)は、図7(a)の被検体の像に対し分散0.25のノイズを重畳した場合のx方向及びy方向の微分位相像である。図10(b)は、図10(a)の微分位相像に対し本実施例の位相像取得方法を適用した後、位相像をx方向、y方向それぞれに関し微分することで生成した微分位相像である。図10(a)の元の微分位相像においてみられたアーチファクトやノイズが、図10(b)では改善されており、高品位な微分位相像が得られている。すなわち本実施例の手法は、微分位相像の画質を向上するという目的にも利用できる。画像診断や検査の内容によっては、位相像をみるよりも、エッジが強調された微分位相像のほうが有用な場合もある。よって、高画質の微分位相像を生成できる機能も実用上非常に有利である。構成としては、例えば、図3の位相データ算出部34の後段に、位相像のデータを微分する微分位相データ算出部を設ければよい。
以上述べたように、本実施例の位相像取得方法によれば、位相イメージングにおいて、ノイズが低減された高品位の位相像および微分位相像を生成することが可能である。そして、位相接続処理を施したことにより、位相折りたたみに起因する位相のゆがみやアーチファクトを改善することもできる。しかも本実施例では、二階微分位相像に対して位相接続処理を施す構成としたので、従来(微分位相像に対して位相接続処理を施す)に比べて、位相接続処理が格段に簡単となり、また高速に処理できる。なぜなら、二階微分位相像では位相飛びがスパイク状の異常値となって現れるため、本実施例のように簡易なアルゴリズムで修正できると共に、近傍ピクセルの値を調べる必要がないため積分演算の中に演算子Wとして組み込むことも容易だからである。
なお、上記実施形態および実施例は本発明の一具体例を示したものにすぎず、本発明の
範囲を限定する趣旨のものではない。例えば、上記実施例では、直交する2方向について演算を行ったが、第1の方向及び第2の方向については互いに交差する限りにおいて任意に設定することができる。また上記実施例ではトールボット干渉計を例示したが、本発明は他のシアリング干渉計にも好適に適用することができる。
1:撮像装置、10:干渉計、11:情報処理装置、32:微分位相データ抽出部、33:二階微分位相データ算出部、34:位相データ算出部、210:被検体、810:干渉パターン

Claims (7)

  1. シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する撮像装置であって、
    被検体を透過した電磁波により形成された干渉パターンのデータから、第1の方向に関する位相の変化を表す第1の微分位相データと、前記第1の方向に交差する第2の方向に関する位相の変化を表す第2の微分位相データとを抽出する微分位相データ抽出手段と、
    抽出された前記第1の微分位相データを前記第1の方向に関して微分することによって第1の二階微分位相データを算出すると共に、抽出された前記第2の微分位相データを前記第2の方向に関して微分することによって第2の二階微分位相データを算出する二階微分位相データ算出手段と、
    前記第1および第2の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する位相データ算出手段と、
    を有することを特徴とする撮像装置。
  2. 前記位相データ算出手段は、位相像のデータを算出する際に、前記第1および第2の二階微分位相データに対して位相接続処理を施すことを特徴とする請求項1に記載の撮像装置。
  3. 前記位相接続処理は、前記第1および第2の二階微分位相データの値のうち所定の範囲を外れた値を、前記所定の範囲内に収める処理であることを特徴とする請求項2に記載の撮像装置。
  4. 前記二階微分方程式は、既知関数である前記第1の二階微分位相データと前記第2の二階微分位相データとの和が、未知関数である前記被検体の位相像のデータを前記第1の方向に二階微分したものと前記第2の方向に二階微分したものとの和に等しいとする方程式であることを特徴とする請求項1〜3のうちいずれか1項に記載の撮像装置。
  5. 前記位相データ算出手段によって算出された前記被検体の位相像のデータを微分することによって、前記被検体の微分位相像のデータを生成する微分位相データ算出手段をさらに有することを特徴とする請求項1〜4のうちいずれか1項に記載の撮像装置。
  6. シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する画像処理方法であって、
    情報処理装置が、被検体を透過した電磁波により形成された干渉パターンのデータから、第1の方向に関する位相の変化を表す第1の微分位相データと、前記第1の方向に交差する第2の方向に関する位相の変化を表す第2の微分位相データとを抽出する工程と、
    情報処理装置が、抽出された前記第1の微分位相データを前記第1の方向に関して微分することによって第1の二階微分位相データを算出すると共に、抽出された前記第2の微分位相データを前記第2の方向に関して微分することによって第2の二階微分位相データを算出する工程と、
    情報処理装置が、前記第1および第2の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する工程と、
    を有することを特徴とする画像処理方法。
  7. 請求項6に記載の画像処理方法の各工程を情報処理装置に実行させることを特徴とするプログラム。
JP2011248717A 2011-11-14 2011-11-14 撮像装置および画像処理方法 Expired - Fee Related JP5868132B2 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2011248717A JP5868132B2 (ja) 2011-11-14 2011-11-14 撮像装置および画像処理方法
US14/347,270 US9019479B2 (en) 2011-11-14 2012-11-02 Imaging apparatus and image processing method
EP12798431.8A EP2779902A1 (en) 2011-11-14 2012-11-02 Imaging apparatus and image processing method
PCT/JP2012/079051 WO2013073453A1 (en) 2011-11-14 2012-11-02 Imaging apparatus and image processing method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2011248717A JP5868132B2 (ja) 2011-11-14 2011-11-14 撮像装置および画像処理方法

Publications (2)

Publication Number Publication Date
JP2013102951A true JP2013102951A (ja) 2013-05-30
JP5868132B2 JP5868132B2 (ja) 2016-02-24

Family

ID=47324324

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2011248717A Expired - Fee Related JP5868132B2 (ja) 2011-11-14 2011-11-14 撮像装置および画像処理方法

Country Status (4)

Country Link
US (1) US9019479B2 (ja)
EP (1) EP2779902A1 (ja)
JP (1) JP5868132B2 (ja)
WO (1) WO2013073453A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016087473A (ja) * 2014-11-06 2016-05-23 キヤノン株式会社 2次元のフリンジ模様の復調において軸外周波数を低減するための非線形処理の方法、記憶媒体、および撮影システム
JP2016214615A (ja) * 2015-05-21 2016-12-22 コニカミノルタ株式会社 位相画像処理方法、位相画像処理装置、画像処理装置及び画像処理プログラム
JP2018048996A (ja) * 2016-09-14 2018-03-29 リコーエレメックス株式会社 検査システム
US10422744B2 (en) * 2016-10-04 2019-09-24 Industrial Technology Research Institute Interferometer and imaging method therefor

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103604508B (zh) * 2013-12-02 2016-03-02 青岛大学 一种自适应消除倾斜误差的波前重建方法
DE102014224638A1 (de) * 2014-12-02 2016-06-02 Olympus Soft Imaging Solutions Gmbh Digitales Bilderfassungssystem und Verfahren zur Fehlerkorrektur in einem solchen System
TWI639809B (zh) * 2016-10-04 2018-11-01 財團法人工業技術研究院 干涉儀及其成像方法
EP3327413B1 (en) * 2016-11-24 2022-04-27 IMEC vzw A method, an apparatus and a system for holographic wavefront sensing
CN112683848B (zh) * 2020-12-21 2022-09-02 中国科学院上海光学精密机械研究所 光学相干层析成像系统色散补偿方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0784202A2 (en) * 1996-01-10 1997-07-16 Hitachi, Ltd. Phase-contrast X-ray CT apparatus
US20040069949A1 (en) * 2002-10-04 2004-04-15 Fuji Photo Film Co., Ltd. Method, apparatus and program for restoring phase information
WO2010050483A1 (ja) * 2008-10-29 2010-05-06 キヤノン株式会社 X線撮像装置およびx線撮像方法
WO2010146503A1 (en) * 2009-06-16 2010-12-23 Koninklijke Philips Electronics N. V. Correction method for differential phase contrast imaging
JP2011136156A (ja) * 2009-12-04 2011-07-14 Canon Inc X線撮像装置およびx線撮像方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5481356A (en) * 1994-04-25 1996-01-02 Northwestern University Apparatus and method for nondestructive testing using additive-subtractive phase-modulated interferometry
US5424743A (en) * 1994-06-01 1995-06-13 U.S. Department Of Energy 2-D weighted least-squares phase unwrapping
WO2006058187A2 (en) * 2004-11-23 2006-06-01 Robert Eric Betzig Optical lattice microscopy
EP1696201A1 (de) * 2005-02-23 2006-08-30 Leica Geosystems AG Phasenrauschkompensation für interferometrische Absolutdistanzmesser
US7538891B1 (en) * 2005-09-30 2009-05-26 California Institute Of Technology Surface characterization based on lateral shearing of diffracted wave fronts to measure in-plane and out-of-plane displacement gradient fields
US8184298B2 (en) * 2008-05-21 2012-05-22 The Board Of Trustees Of The University Of Illinois Spatial light interference microscopy and fourier transform light scattering for cell and tissue characterization

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0784202A2 (en) * 1996-01-10 1997-07-16 Hitachi, Ltd. Phase-contrast X-ray CT apparatus
US20040069949A1 (en) * 2002-10-04 2004-04-15 Fuji Photo Film Co., Ltd. Method, apparatus and program for restoring phase information
WO2010050483A1 (ja) * 2008-10-29 2010-05-06 キヤノン株式会社 X線撮像装置およびx線撮像方法
WO2010146503A1 (en) * 2009-06-16 2010-12-23 Koninklijke Philips Electronics N. V. Correction method for differential phase contrast imaging
JP2011136156A (ja) * 2009-12-04 2011-07-14 Canon Inc X線撮像装置およびx線撮像方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
JPN6015049278; C. Kottler, C. David, F. Pfeiffer, and O. Bunk: 'A two-directional approach for grating based differential phase contrast imaging using hard x-rays' OPTICS EXPRESS Vol. 15, Issue 3, 20070205, 1175-1181 *
JPN6015049280; A. Momose, W. Yashiro, Y. Takeda, Y. Suzuki and T. Hattori: 'Phase Tomography by X-ray Talbot Interferometry for Biological Imaging' Japanese Journal of Applied Physics Vol. 45, No. 6A, 2006, 5254-5262 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016087473A (ja) * 2014-11-06 2016-05-23 キヤノン株式会社 2次元のフリンジ模様の復調において軸外周波数を低減するための非線形処理の方法、記憶媒体、および撮影システム
JP2016214615A (ja) * 2015-05-21 2016-12-22 コニカミノルタ株式会社 位相画像処理方法、位相画像処理装置、画像処理装置及び画像処理プログラム
JP2018048996A (ja) * 2016-09-14 2018-03-29 リコーエレメックス株式会社 検査システム
US10422744B2 (en) * 2016-10-04 2019-09-24 Industrial Technology Research Institute Interferometer and imaging method therefor

Also Published As

Publication number Publication date
JP5868132B2 (ja) 2016-02-24
WO2013073453A1 (en) 2013-05-23
US9019479B2 (en) 2015-04-28
US20140241632A1 (en) 2014-08-28
EP2779902A1 (en) 2014-09-24

Similar Documents

Publication Publication Date Title
JP5868132B2 (ja) 撮像装置および画像処理方法
JP6105586B2 (ja) エネルギー高感度検出の微分位相コントラストイメージング
JP5777360B2 (ja) X線撮像装置
JP2017506925A (ja) 微分位相コントラスト撮像からの位相回復
JP6214819B1 (ja) 微分位相コントラストx線撮像における暗視野信号の最適なエネルギ加重
EP2745265B1 (en) Frequency dependent combination of x-ray images of different modalities
JP2013050441A (ja) 波面測定装置、波面測定方法、及びプログラム並びにx線撮像装置
JP6670398B2 (ja) 暗視野又は位相コントラストx線撮像における特徴抑制
JP2013536723A (ja) 微分位相差イメージングにおける正規化位相回復
US20150362444A1 (en) Phase information acquisition apparatus and imaging system
Sadi et al. Removal of ring artifacts in computed tomographic imaging using iterative center weighted median filter
JP2015190776A (ja) 画像処理装置および撮像システム
JP2014140632A (ja) 演算装置、画像取得方法、プログラム、及びx線撮像システム
US20140114615A1 (en) Imaging apparatus and program and method for analyzing interference pattern
JP5220383B2 (ja) 放射線撮像装置
Li et al. Spatial resolution characterization of differential phase contrast CT systems via modulation transfer function (MTF) measurements
JP6116222B2 (ja) 演算装置、プログラム及び撮像システム
US9330456B2 (en) Systems and methods for regularized Fourier analysis in x-ray phase contrast imaging
JP2016106721A (ja) 画像処理装置および画像処理方法
JP5527481B2 (ja) X線診断装置およびx線診断用プログラム
Harasse et al. X-ray phase laminography with Talbot interferometer
CN115398471A (zh) 基于滑动窗口相位恢复的针对暗场成像的偏差校正
Gelikonov et al. Numerical method for axial motion correction in optical coherence tomography
JP2016061608A (ja) 画像処理方法、画像処理装置、撮像システム
Ling et al. Patch based grid artifact suppressing in digital mammography

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20141113

A711 Notification of change in applicant

Free format text: JAPANESE INTERMEDIATE CODE: A711

Effective date: 20141113

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20141113

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20151208

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160105

R151 Written notification of patent or utility model registration

Ref document number: 5868132

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151

LAPS Cancellation because of no payment of annual fees