(57)【要約】      （修正有）【課題】推定精度の向上に寄与できるように擬似的な荷電粒子の散乱点を生成することができる画像処理装置、画像処理システム及び画像処理方法を提供する。【解決手段】画像処理装置は、推定部２１と、第１生成部２３と、第２生成部２５と、を備える。推定部は、実荷電粒子の入射方向を示す実入射ベクトルと射出方向を示す実射出ベクトルとの組である実ベクトル対を用いて、実荷電粒子の散乱点である実散乱点を推定する。第１生成部は、疑似荷電粒子の散乱点である疑似散乱点を実散乱点の近傍領域に生成する。第２生成部は、疑似散乱点を通る疑似荷電粒子の入射方向を示す疑似入射ベクトルと射出方向を示す疑似射出ベクトルとの組である疑似ベクトル対を生成する。【選択図】図２
2016180664.tif 000002
【特許請求の範囲】
【請求項１】
実荷電粒子の入射方向を示す実入射ベクトルと射出方向を示す実射出ベクトルとの組である実ベクトル対を用いて、前記実荷電粒子の散乱点である実散乱点を推定する推定部と、
疑似荷電粒子の散乱点である疑似散乱点を前記実散乱点の近傍領域に生成する第１生成部と、
前記疑似散乱点を通る前記疑似荷電粒子の入射方向を示す疑似入射ベクトルと射出方向を示す疑似射出ベクトルとの組である疑似ベクトル対を生成する第２生成部と、
を備える画像処理装置。
【請求項２】
前記第１生成部は、前記実ベクトル対を用いて前記実荷電粒子の前記実散乱点における散乱角度を求め、前記散乱角度の絶対値が小さいほど密になるように前記近傍領域に前記疑似散乱点を複数生成する請求項１に記載の画像処理装置。
【請求項３】
前記第１生成部は、前記実散乱点に近いほど密になるように前記近傍領域に前記疑似散乱点を複数生成する請求項１に記載の画像処理装置。
【請求項４】
前記近傍領域は、１以上の実数値で構成される近接パラメータに基づいて、範囲が定まる請求項１に記載の画像処理装置。
【請求項５】
前記疑似入射ベクトルは、前記実入射ベクトルと同じ向きのベクトルであり、
前記疑似射出ベクトルは、前記実射出ベクトルと同じ向きのベクトルである請求項１に記載の画像処理装置。
【請求項６】
前記疑似入射ベクトルは、前記実荷電粒子の入射点から前記疑似散乱点へ向かうベクトルであり、
前記疑似射出ベクトルは、前記疑似散乱点から前記実荷電粒子の射出点へ向かうベクトルである請求項１に記載の画像処理装置。
【請求項７】
前記実散乱点は、前記実入射ベクトルと平行で前記実荷電粒子の入射点を通る直線と、前記実射出ベクトルと平行で前記実荷電粒子の射出点を通る直線と、に最も近い点である請求項１に記載の画像処理装置。
【請求項８】
前記推定部は、前記実荷電粒子の入射点及び射出点を更に用いて、前記実散乱点を推定する請求項１に記載の画像処理装置。
【請求項９】
前記推定部は、複数の実荷電粒子それぞれ毎に、前記実ベクトル対を用いて前記実散乱点を推定し、
前記第１生成部は、前記複数の実散乱点それぞれ毎に、前記実散乱点の近傍領域に前記疑似散乱点を生成し、
前記第２生成部は、前記複数の疑似散乱点それぞれ毎に、前記疑似ベクトル対を生成し、
前記複数の実荷電粒子の前記実ベクトル対、及び前記実散乱点、並びに複数の疑似散乱点の前記疑似ベクトル対、及び前記疑似散乱点を用いて、前記複数の実荷電粒子が透過した対象物の内部構造を示す画像を生成する画像生成部を更に備える請求項１に記載の画像処理装置。
【請求項１０】
前記画像生成部は、前記実荷電粒子の入射点及び射出点を更に用いて、前記画像を生成する請求項９に記載の画像処理装置。
【請求項１１】
前記画像生成部は、前記複数の疑似荷電粒子それぞれに与えられた前記実散乱点から前記疑似散乱点までの距離の増加に対して少なくとも増加しない重みを更に用いて、前記画像を生成する請求項９に記載の画像処理装置。
【請求項１２】
請求項９に記載の画像処理装置と、
前記複数の実荷電粒子それぞれの前記入射点、前記実入射ベクトル、前記射出点、及び前記実射出ベクトルを検出する複数の検出部と、
前記画像を表示する表示部と、
を備える画像処理システム。
【請求項１３】
実荷電粒子の入射方向を示す実入射ベクトルと射出方向を示す実射出ベクトルとの組である実ベクトル対を用いて、前記実荷電粒子の散乱点である実散乱点を推定する推定ステップと、
疑似荷電粒子の散乱点である疑似散乱点を前記実散乱点の近傍領域に生成する第１生成ステップと、
前記疑似散乱点を通る前記疑似荷電粒子の入射方向を示す疑似入射ベクトルと射出方向を示す疑似射出ベクトルとの組である疑似ベクトル対を生成する第２生成ステップと、
を含む画像処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、画像処理装置、画像処理システム及び画像処理方法に関する。
【背景技術】
【０００２】
荷電粒子の散乱を利用して可視化対象の内部を可視化する荷電粒子トモグラフィが知られている。このような荷電粒子トモグラフィでは、例えば、可視化の対象物を挟めるように荷電粒子を観測するための検出器を設置し、これらの検出器から観測された荷電粒子の入射位置、入射方向、射出位置、及び射出方向を得、観測された荷電粒子の散乱角度及び軌跡を推定し、推定した散乱角度及び軌跡に基づいて対象物の内部構造を推定し、可視化した画像を生成する。
【０００３】
ところで、荷電粒子トモグラフィでは、観測された荷電粒子の散乱角度の統計を用いて構造物の推定した内部構造を可視化するため、推定精度は、観測された荷電粒子の数に依存する。つまり、構造物の内部構造の推定精度を高めるためには、多くの荷電粒子を観測する必要があるが、検出器の位置や観測時間によっては、荷電粒子を十分に観測できない場合もある。
【０００４】
このため、観測された荷電粒子の入射方向を示す入射ベクトル及び射出方向を示す射出ベクトルの周囲に放射状に疑似的な入射ベクトル及び射出ベクトルをランダムに生成し、生成した疑似的な入射ベクトル及び射出ベクトルを用いて擬似的な荷電粒子の散乱点を生成し、荷電粒子の数を疑似的に増やす手法が提案されている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】国際公開第２０１２／１６７３８１号
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上述したような従来技術では、実在する荷電粒子の散乱点の位置を考慮せずに擬似的な荷電粒子の散乱点を生成するため、擬似的な荷電粒子の散乱点が実在する荷電粒子の散乱点から空間的に大きく離れて配置される可能性があり、このような場合、対象物の内部構造の推定精度を高めることは難しい。
【０００７】
本発明が解決しようとする課題は、推定精度の向上に寄与できるように擬似的な荷電粒子の散乱点を生成することができる画像処理装置、画像処理システム及び画像処理方法を提供することである。
【課題を解決するための手段】
【０００８】
実施形態の画像処理装置は、推定部と、第１生成部と、第２生成部と、を備える。推定部は、実荷電粒子の入射方向を示す実入射ベクトルと射出方向を示す実射出ベクトルとの組である実ベクトル対を用いて、前記実荷電粒子の散乱点である実散乱点を推定する。第１生成部は、疑似荷電粒子の散乱点である疑似散乱点を前記実散乱点の近傍領域に生成する。第２生成部は、前記疑似散乱点を通る前記疑似荷電粒子の入射方向を示す疑似入射ベクトルと射出方向を示す疑似射出ベクトルとの組である疑似ベクトル対を生成する。
【図面の簡単な説明】
【０００９】
【図１】本実施形態の画像処理システムの例を示す構成図。
【図２】本実施形態の画像処理部の例を示す構成図。
【図３】本実施形態の実散乱点の推定例の説明図。
【図４】本実施形態の疑似散乱点の生成例の説明図。
【図５】本実施形態の疑似ベクトル対の生成例の説明図。
【図６】本実施形態の疑似軌跡の例の説明図。
【図７】本実施形態の疑似軌跡の例の説明図。
【図８】本実施形態の処理例を示すフローチャート。
【図９】本実施形態の画像処理部のハードウェア構成例を示す図。
【発明を実施するための形態】
【００１０】
以下、添付図面を参照しながら、実施形態を詳細に説明する。なお、以下の実施形態では、実在する荷電粒子を実荷電粒子と称し、実在しない疑似的な荷電粒子を疑似荷電粒子と称する。
【００１１】
図１は、本実施形態の画像処理システム１の一例を示す構成図である。図１に示すように、画像処理システム１は、検出部１０Ａ、１０Ｂ（複数の検出部の一例）と、画像処理部２０（画像処理装置の一例）と、表示部３０と、を備える。
【００１２】
検出部１０Ａ、１０Ｂは、実荷電粒子を検出するものであり、例えば、公知の荷電粒子検出器により実現できる。また検出部１０Ａ、１０Ｂは、可視化対象の対象物５を挟めるように配置されている。図１に示す例では、検出部１０Ａ、１０Ｂは、対象物５を上下方向に挟むように配置されているが、これに限定されず、対象物５を挟めるように配置されていればその方向は問わない。また、対象物５は、構造物などの固定物を想定しているが、別の線源を用意するなど短時間で十分な量の線源を照射できれば、移動可能な移動体などであってもよい。対象物５が移動体の場合、検出部１０Ａ、１０Ｂは、例えば、対象物５の移動経路上のある地点において、移動中の対象物５を挟めるように配置されていればよい。
【００１３】
また本実施形態では、実荷電粒子が宇宙線ミュオン（以下、単に「ミュオン」と称する）である場合を例に取り説明するが、これに限定されず、宇宙線や他の線源から発生したミュオン以外の荷電粒子であってもよい。なお、図１に示す例では、実荷電粒子は、検出部１０Ａの上空から降り注いでおり、降り注がれた実荷電粒子は、主として、検出部１０Ａ、対象物５、検出部１０Ｂの順に透過することを想定している。
【００１４】
従って、図１に示す例では、検出部１０Ａは、主として、自身（検出部１０Ａ）を透過した複数の実荷電粒子それぞれの入射点である実入射点及び入射ベクトルである実入射ベクトルを検出し、検出部１０Ｂは、主として、自身（検出部１０Ｂ）を透過した複数の実荷電粒子それぞれの射出点である実射出点及び射出ベクトルである実射出ベクトルを検出する。
【００１５】
実入射点は、検出部１０Ａにおいて実荷電粒子が検出された位置である。実入射ベクトルは、実荷電粒子の実入射点における移動方向である入射方向を示す任意の長さのベクトルである。実射出点は、検出部１０Ｂにおいて実荷電粒子が検出された位置である。実射出ベクトルは、実荷電粒子の実射出点における移動方向である射出方向を示す任意の長さのベクトルである。なお、実荷電粒子は、対象物５内において散乱するため、実入射ベクトルと実射出ベクトルとが同じ方向を示すとは限らない。
【００１６】
但し、実荷電粒子は、必ずしも、検出部１０Ａ、対象物５、検出部１０Ｂの順に透過するとは限らず、検出部１０Ｂ、対象物５、検出部１０Ａの順に透過する場合も考えられる。この場合、検出部１０Ｂは、実荷電粒子の実入射点及び実入射ベクトルを検出し、検出部１０Ａは、実荷電粒子の実射出点及び実射出ベクトルを検出する。
【００１７】
つまり、検出部１０Ａ、１０Ｂは、複数の実荷電粒子それぞれの実入射点、実入射ベクトル、実射出点、及び実射出ベクトルを検出する。なお、実荷電粒子の速度は非常に高速であるため（ミュオンであれば光速の９９．５％の速度）、一定時間内（ほぼ同時）に検出された実荷電粒子を、同一の実荷電粒子とみなすことで、検出部１０Ａ、１０Ｂは、複数の実荷電粒子それぞれの実入射点、実入射ベクトル、実射出点、及び実射出ベクトルを検出できる。以下では、１つの実荷電粒子の実入射ベクトルと実射出ベクトルとの組を実ベクトル対と称し、１つの実荷電粒子の実入射点、実射出点、及び実ベクトル対をまとめて実荷電粒子データと称し、複数の実荷電粒子の実荷電粒子データをまとめて観測データと称する。
【００１８】
画像処理部２０は、検出部１０Ａ、１０Ｂから入力された観測データを用いて、対象物５の内部構造を推定し、可視化した画像を生成するものである。画像処理部２０の詳細については、後述する。
【００１９】
表示部３０は、画像処理部２０から入力された画像を表示するものであり、例えば、公知のディスプレイにより実現できる。
【００２０】
図２は、本実施形態の画像処理部２０の一例を示す構成図である。図２に示すように、画像処理部２０は、推定部２１と、第１生成部２３と、第２生成部２５と、画像生成部２７と、を備える。
【００２１】
推定部２１、第１生成部２３、第２生成部２５、及び画像生成部２７は、例えば、ＣＰＵ（Central  Processing  Unit）などの処理装置にプログラムを実行させること、即ち、ソフトウェアにより実現してもよいし、ＩＣ（Integrated  Circuit）などのハードウェアにより実現してもよいし、ソフトウェア及びハードウェアを併用して実現してもよい。
【００２２】
推定部２１には、検出部１０Ａ、１０Ｂから観測データが入力される。推定部２１は、検出部１０Ａ、１０Ｂにより検出された複数の実荷電粒子それぞれ毎に、入力された観測データに含まれる当該実荷電粒子の実荷電粒子データに含まれる実ベクトル対を用いて、当該実荷電粒子の散乱点である実散乱点を推定し、第１生成部２３及び画像生成部２７に出力する。詳細には、推定部２１は、実荷電粒子毎に、入力された観測データに含まれる当該実荷電粒子の実荷電粒子データを用いて、当該実荷電粒子の実散乱点を推定する。
【００２３】
実散乱点の推定には、例えば、L.  J.  Schultz,  "Cosmic  ray  muon  radiography."  Diss.  Portland  State  University,  2003.に開示されている最近点法（ＰｏＣＡ：Point  of  Closest  Approach）を用いることができる。ＰｏＣＡでは、実入射ベクトルと平行で実入射点を通る直線である実入射線と、実射出ベクトルと平行で実射出点を通る直線である実射出線と、に最も近い点である最近点を実散乱点とする。
【００２４】
例えば、実入射点１０１、実入射ベクトル１０２、実射出点１０３、及び実射出ベクトル１０４が、図３に示すような位置関係にある場合、推定部２１は、実入射ベクトル１０２と平行で実入射点１０１を通る実入射線１０５と、実射出ベクトル１０４と平行で実射出点１０３を通る実射出線１０６と、の最近点を実散乱点１０７に推定する。
【００２５】
但し、３次元空間では実入射線と実射出線とが必ずしも交わるとは限らないため、推定部２１は、実入射線及び実射出線上に実散乱点を推定できない場合がある。この場合、推定部２１は、実入射線と実射出線との交点が実散乱点となるように、実入射点及び実射出点の少なくともいずれかの位置を修正したり、実入射ベクトル及び実射出ベクトルの少なくともいずれかの方向を修正したりしてもよい。
【００２６】
第１生成部２３には、検出部１０Ａ、１０Ｂから観測データが入力されるとともに、推定部２１から複数の実荷電粒子それぞれの実散乱点が入力される。また、第１生成部２３には、画像処理部２０の外部から１以上の実数値で構成される近接パラメータが入力される。近接パラメータは、例えば、画像処理システム１を使用するユーザなどが図示せぬ入力装置を用いて入力できる。なお、近接パラメータの詳細については、後述する。
【００２７】
第１生成部２３は、推定部２１により入力された複数の実散乱点それぞれ毎に、疑似荷電粒子の散乱点である疑似散乱点を当該実散乱点の近傍領域に生成し、第２生成部２５に出力する。
【００２８】
ここで、第１生成部２３は、入力された近接パラメータに基づいて、近傍領域の範囲を定める。例えば、第１生成部２３は、近接パラメータを構成する１以上の実数値の値が小さいほど、近傍領域の範囲を小さくする。なお、近接パラメータを構成する１以上の実数値の値や数は、第１生成部２３の設計に依存する。
【００２９】
例えば、実荷電粒子データの検出誤差による実散乱点の位置の推定誤差が標準偏差σｐｏｃの正規分布で表されると考える場合、第１生成部２３は、実散乱点を中心とし、σｐｏｃに近接パラメータを乗じた距離の点を境界とする領域を近傍領域とすることができる。この場合、近接パラメータを構成する１以上の実数値は１つの変数（例えば、α）であり、近傍領域は、実散乱点を中心とする半径ασｐｏｃの円や球とすることができる。
【００３０】
また例えば、実荷電粒子データの検出誤差による実散乱点の位置の推定誤差の分布が空間的に非一様なばらつきを持つと考える場合、第１生成部２３は、実散乱点を中心、ばらつきの大きな２つ以上の方向を軸とし、軸方向の標準偏差に近傍パラメータを乗じた距離の点を境界とする領域を近傍領域とすることができる。この場合、近接パラメータを構成する１以上の実数値は、それぞれ、軸毎に定まる変数となる。従って、軸の数が２つの場合、変数は２つとなり、軸の数が３つの場合、変数は３つとなる。なお、近傍領域は、軸の数が２つの場合、楕円とすることができ、軸の数が３つの場合、楕円体とすることができる。
【００３１】
第１生成部２３は、このようにして定めた近傍領域内に１以上の疑似散乱点を生成（配置）する。例えば、第１生成部２３は、図４に示すように、実散乱点１０７を中心に楕円形状に定められた近傍領域１０８内に、疑似散乱点１１７、１２７、１３７、１４７を生成する。なお、第１生成部２３は、疑似散乱点を近傍領域内に、一様に生成してもよいし、ランダムに生成してもよい。
【００３２】
また、実荷電粒子は、実散乱点の散乱角度（但し散乱角度は、０を中心とするガウス分布に従うので、詳細には、「散乱角度の絶対値」、第１生成部２３の説明における散乱角度については、以下同様）が大きいほど対象物５内の複数の物体中で散乱する多重散乱を起こし、その影響を受けて実散乱点の位置の推定精度が低くなると考えて、第１生成部２３は、散乱角度が小さいほど近傍領域内に密になるように疑似散乱点を生成してもよい。
【００３３】
つまり、第１生成部２３は、推定部２１により入力された複数の実散乱点それぞれ毎に、検出部１０Ａ、１０Ｂから入力された観測データに含まれる当該実散乱点が推定された実荷電粒子の実荷電粒子データに含まれる実ベクトル対を用いて、当該実散乱点の散乱角度を求める。そして第１生成部２３は、求めた散乱角度が小さいほど密になるように当該実散乱点に対し定めた近傍領域に疑似散乱点を複数生成するようにしてもよい。
【００３４】
なお、散乱角度は、例えば、実ベクトル対を構成する実入射ベクトル及び実射出ベクトルを直交する２つの平面に射影し、各平面上での実入射ベクトル及び実射出ベクトルそれぞれの角度を求め、各平面での実入射ベクトルと実射出ベクトルとの角度の変位（差分）とすることができる。ここで、各平面の散乱角度を独立した別の散乱として扱うのであれば、１つの実荷電粒子を散乱角度以外が同一の２つの別の実荷電粒子として扱えばよく、各平面の散乱角度を何らかの方法で合成して１つの散乱角度として扱うのであれば、例えば、θｘ、θｙを各平面の散乱角度だとすると、θ２＝（θｘ２＋θｙ２）／２となるθを散乱角度とすればよい。但し、散乱角度の求め方は、これに限定されず、正負が必要なければ、単純に、実ベクトル対を構成する実入射ベクトルと実射出ベクトルとの成す角として求めてもよいし、他の方法で求めてもよい。
【００３５】
また、実散乱点に近いほど実荷電粒子の散乱と類似した散乱が起こりやすいと考えて、第１生成部２３は、実散乱点に近いほど密になるように近傍領域内に疑似散乱点を生成してもよい。
【００３６】
なお、第１生成部２３は、疑似荷電粒子に対し、実散乱点から当該疑似荷電粒子の疑似散乱点までの距離の増加に対して少なくとも増加しない重み、即ち、距離が大きくなると減少するか、距離が大きくなっても変化しない重みを与えるようにしてもよい。
【００３７】
例えば、第１生成部２３は、実散乱点から疑似散乱点までの距離を変数とする正規分布で疑似荷電粒子に重みを与えてもよい。また例えば、第１生成部２３は、実散乱点から疑似散乱点までの距離が閾値以下の場合は疑似荷電粒子に定数の重みを与え、実散乱点から疑似散乱点までの距離が閾値よりも大きい場合は距離が大きくなるほど定数を小さくした重みを疑似荷電粒子に与えてもよい。
【００３８】
疑似荷電粒子に重みを与える場合、当該重みは画像生成部２７で使用されるため、第１生成部２３は、当該重みを画像生成部２７に出力する。
【００３９】
第２生成部２５には、検出部１０Ａ、１０Ｂから観測データが入力されるとともに、第１生成部２３から複数の実散乱点それぞれに対し生成された１以上の疑似散乱点が入力される。第２生成部２５は、入力された疑似散乱点それぞれ毎に、当該疑似散乱点を通る疑似荷電粒子の入射方向を示す疑似入射ベクトルと射出方向を示す疑似射出ベクトルとの組である疑似ベクトル対を生成し、画像生成部２７に出力する。
【００４０】
なお、疑似荷電粒子の軌跡では、疑似入射点から疑似散乱点の方向は疑似入射ベクトルと一致し、疑似散乱点から疑似射出点の方向は疑似射出ベクトルと一致するので、第２生成部２５は、疑似ベクトル対ではなく、疑似入射点及び疑似射出点を出力するようにしても、画像生成部２７において同様の処理が可能である。
【００４１】
実散乱点に対し生成された疑似散乱点を通る疑似荷電粒子の軌跡である疑似軌跡が、当該実散乱点を通る実荷電粒子の軌跡である実軌跡と大きく異なると、疑似散乱点が実散乱点から大きく離れた場合と同様、対象物の内部構造の推定精度を高めることは難しい。なお本実施形態では、実軌跡は、実荷電粒子の実入射点及び実散乱点を結ぶ線分と、実散乱点及び実射出点を結ぶ線分と、を接続した軌跡とする。また、疑似軌跡は、疑似荷電粒子の入射点である疑似入射点及び疑似散乱点を結ぶ線分と、疑似散乱点及び疑似荷電粒子の射出点である疑似射出点を結ぶ線分と、を接続した軌跡とする。このため本実施形態では、第２生成部２５は、疑似軌跡が実軌跡の近傍を通るように、疑似入射点及び疑似射出点を生成して、疑似入射ベクトル及び疑似射出ベクトルを生成する。但し、第２生成部２５は、疑似軌跡が実軌跡の近傍を通るように、疑似入射ベクトル及び疑似射出ベクトルを生成して、疑似入射点及び疑似射出点を生成するようにしてもよい。
【００４２】
疑似入射点は、例えば、疑似散乱点が生成された実散乱点を通る実荷電粒子の実入射ベクトルと平行で当該疑似散乱点を通る直線と、当該実入射ベクトルを検出した検出部１０Ａ又は１０Ｂと、の交点とすることができる。同様に、疑似射出点は、例えば、疑似散乱点が生成された実散乱点を通る実荷電粒子の実射出ベクトルと平行で当該疑似散乱点を通る直線と、当該実射出ベクトルを検出した検出部１０Ａ又は１０Ｂと、の交点とすることができる。なお、疑似散乱点が生成された実散乱点を通る実荷電粒子の実入射ベクトルや実射出ベクトルは、検出部１０Ａ、１０Ｂから入力された観測データから特定できる。
【００４３】
そして、第２生成部２５は、疑似入射点を始点とし疑似散乱点を終点とするベクトルを疑似入射ベクトルとして生成し、疑似散乱点を始点とし疑似射出点を終点とするベクトルを疑似射出ベクトルとして生成する。
【００４４】
例えば、第２生成部２５は、図５に示すように、疑似散乱点１１７を通る疑似荷電粒子については、疑似入射点１１１を始点とし疑似散乱点１１７を終点とする疑似入射ベクトル１１２、及び疑似散乱点１１７を始点とし疑似射出点１１３を終点とする疑似射出ベクトル１１４の組である疑似ベクトル対を生成する。また第２生成部２５は、疑似散乱点１２７を通る疑似荷電粒子については、疑似入射点１２１を始点とし疑似散乱点１２７を終点とする疑似入射ベクトル１２２、及び疑似散乱点１２７を始点とし疑似射出点１２３を終点とする疑似射出ベクトル１２４の組である疑似ベクトル対を生成する。また第２生成部２５は、疑似散乱点１３７を通る疑似荷電粒子については、疑似入射点１３１を始点とし疑似散乱点１３７を終点とする疑似入射ベクトル１３２、及び疑似散乱点１３７を始点とし疑似射出点１３３を終点とする疑似射出ベクトル１３４の組である疑似ベクトル対を生成する。また第２生成部２５は、疑似散乱点１４７を通る疑似荷電粒子については、疑似入射点１４１を始点とし疑似散乱点１４７を終点とする疑似入射ベクトル１４２、及び疑似散乱点１４７を始点とし疑似射出点１４３を終点とする疑似射出ベクトル１４４の組である疑似ベクトル対を生成する。
【００４５】
この手法で疑似ベクトル対を生成した場合、図５に示すように、疑似入射ベクトル１１２、１２２、１３２、１４２は、いずれも実入射ベクトル１０２と同じ向きのベクトルとなり、疑似射出ベクトル１１４、１２４、１３４、１４４は、いずれも実射出ベクトル１０４と同じ向きのベクトルとなる。
【００４６】
また、この手法の場合、図６に示すように、実散乱点１０７に対し設定された疑似散乱点を通る疑似荷電粒子の疑似軌跡は、近傍領域１０８が図のように表される場合、近傍領域１０８を底面又は上面とし、高さが実入射点１０１と実散乱点１０７とを結ぶ線分の長さとなる円柱１５１と、近傍領域１０８を底面又は上面とし、高さが実散乱点１０７と実射出点１０３とを結ぶ線分の長さとなる円柱１５２と、を結合した空間内に生成される。つまり、疑似軌跡は、実軌跡の近傍を平行に通る軌跡となる。但し、図６は、３次元の状態を示した図であり、図４及び図５は、２次元の状態を示した図であるため、近傍領域や疑似散乱点の位置がこれらの図面間で厳密に対応しているわけではないことを付記しておく。
【００４７】
なお、第２生成部２５は、疑似入射点を、疑似散乱点が生成された実散乱点を通る実荷電粒子の実入射点と同一の点とし、疑似射出点を、疑似散乱点が生成された実散乱点を通る実荷電粒子の実射出点と同一の点として、疑似入射ベクトル及び疑似射出ベクトルの組である疑似ベクトル対を生成してもよい。
【００４８】
例えば、第２生成部２５は、図７に示すように、疑似散乱点１６７を通る疑似荷電粒子については、実入射点１０１を始点とし疑似散乱点１６７を終点とする疑似入射ベクトル１６２、及び疑似散乱点１６７を始点とし実射出点１０３を終点とする疑似射出ベクトル１６４の組である疑似ベクトル対を生成する。
【００４９】
この手法で疑似ベクトル対を生成した場合、疑似入射ベクトルは、いずれも実入射点から疑似散乱点へ向かうベクトルとなり、疑似射出ベクトルは、いずれも疑似散乱点から実射出点へ向かうベクトルとなる。
【００５０】
また、この手法の場合、図７に示すように、実散乱点１０７に対し設定された疑似散乱点を通る疑似荷電粒子の疑似軌跡は、近傍領域１０８が図のように表される場合、近傍領域１０８を底面とし、実入射点１０１を頂点とする円錐１７１と、近傍領域１０８を底面とし、実射出点１０３を頂点とする円錐１７２と、を結合した空間内に生成される。つまり、疑似軌跡は、実軌跡の近傍を放射状に通る軌跡となる。但し、図７は、３次元の状態を示した図であり、図４及び図５は、２次元の状態を示した図であるため、近傍領域や疑似散乱点の位置がこれらの図面間で厳密に対応しているわけではないことを付記しておく。
【００５１】
画像生成部２７には、検出部１０Ａ、１０Ｂから観測データが入力され、推定部２１から複数の実荷電粒子それぞれの実散乱点が入力され、第１生成部２３から複数の実散乱点それぞれに対し生成された１以上の疑似散乱点が入力され、第２生成部２５から実散乱点それぞれに対し生成された疑似ベクトル対が入力される。画像生成部２７は、複数の実荷電粒子の実ベクトル対、及び実散乱点、並びに複数の疑似散乱点の疑似ベクトル対、及び疑似散乱点を用いて、複数の実荷電粒子が透過した対象物５の内部構造を示す画像を生成し、表示部３０に出力する。
【００５２】
画像生成部２７は、例えば、特許第５６２３７４６号公報に開示された期待値最大化アルゴリズム、上述のL.  J.  Schultz,  "Cosmic  ray  muon  radiography."  Diss.  Portland  State  University,  2003.に開示されたＰｏＣＡを用いた画像再構成法、又は散乱角度の確率密度関数を尤度関数としたエネルギー最小化に基づく手法などを用いて、複数の実荷電粒子が透過した対象物５の内部構造を示す画像を再構成することで生成すればよい。なお、上述の例は、画像再構成手法の一例に過ぎず、荷電粒子の軌跡に基づく画像再構成法であればどのような手法を用いてもよい。
【００５３】
ここでは、実荷電粒子をミュオンとし、ＰｏＣＡを用いた画像再構成法を例に取り、対象物５の内部構造を示す画像の再構成手法について説明する。なお、以下のＰｏＣＡを用いた画像再構成法の説明におけるミュオンの散乱角度は、第１生成部２３の説明における散乱角度とは異なり、散乱角度の絶対値ではなく、散乱角度そのものとして扱う。
【００５４】
まず、ミュオンの散乱の特性について説明する。ミュオンの散乱角度θの中央９８％の分布は、数式（１）に示すように、平均値０の正規分布で近似することができる。
【００５５】
【数１】
2016180664.tif 000003
【００５６】
ここで、標準偏差σθは、物質特性の関数として近似できる。近似式は、様々な式が提案されているが，例えば、数式（２）を用いることができる。
【００５７】
【数２】
2016180664.tif 000004
【００５８】
ここで、ｐはミュオンの運動量を示す。Ｌは対象物５を透過するミュオンの経路長を示す。βｃはミュオンの速度を示す。Ｌｒａｄは、物質の放射長を示し、原子番号が大きくなるに従い減少する。ｅは変数であり、ここでは、ｅ＝１５Ｍｅｖとする。なお、ミュオンの経路長は、ミュオンの実際の軌跡としてもよいし、近似した軌跡の長さとしてもよい。例えば、後述のようにボクセル単位で考える場合、ミュオンの経路長は、ミュオンの軌跡のうちボクセルに含まれる軌跡の長さとしてもよいし、ボクセルの一辺の長さに近似してボクセルの一辺の長さとしてもよい。
【００５９】
ミュオンの経路長を、ミュオンの実際の軌跡とする場合、入射点から散乱点の方向が入射ベクトルと一致し、散乱点から射出点の方向が射出ベクトルと一致する条件下であれば、入射点、散乱点、及び射出点からでも、入射ベクトル、散乱点、及び射出ベクトルからでも求められる。一方、入射点から散乱点の方向が入射ベクトルと一致せず、散乱点から射出点の方向が射出ベクトルと一致しない条件下であれば、入射点から散乱点及び散乱点から射出点を軌跡とするならば、入射点、散乱点、及び射出点から求められ、散乱点から入射ベクトル及び射出ベクトル方向それぞれに伸ばした線を軌跡とするならば、入射ベクトル、散乱点、及び射出ベクトルから求められる。
【００６０】
ミュオンの標準的な運動量をｐ０とし、βｃを１に近似すると、物質の散乱密度λを数式（３）で定義できる。
【００６１】
【数３】
2016180664.tif 000005
【００６２】
ここで、散乱密度λは、放射長がＬｒａｄの物質を透過する際のミュオンの散乱のしやすさを示す指標になるため、散乱密度λの空間分布を推定して可視化することで、対象物５の内部構造を再構成できる。なお、数式（２）及び（３）から、ミュオンの散乱角度θの分布の分散σθ２は、数式（４）で表される。
【００６３】
【数４】
2016180664.tif 000006
【００６４】
ここで、ミュオンの運動量ｐが検出（観測）できない場合は、ｐ＝ｐ０とする。なお、数式（２）に別の式を採用した場合、数式（３）及び（４）は、数式（２）に採用された別の式に合わせて適宜修正すればよい。
【００６５】
ＰｏＣＡを用いた画像再構成法では、ミュオンの入射ベクトル及び射出ベクトルから求めた散乱角度θの分布の分散σθ２を数式（４）に適用することで、散乱密度λを推定し、対象物５の内部構造を示す画像を再構成する。本実施形態では、画像生成部２７は、実荷電粒子だけでなく疑似荷電粒子もミュオンとして扱い、実ベクトル対から求められる散乱角度θだけでなく疑似ベクトル対から求められる散乱角度θも用いて、散乱角度θの分布の分散σθ２を求める。以下では、実ベクトル対及び疑似ベクトル対を区別する必要がない場合は、単に、ベクトル対と称し、実散乱点及び疑似散乱点を区別する必要がない場合は、単に、散乱点と称し、実軌跡及び疑似軌跡を区別する必要がない場合は、単に、軌跡と称する。
【００６６】
ここで、検出部１０Ａ、１０Ｂで挟まれた空間を、３次元の格子点上の立方体の集合であるボクセル空間で表し、各立方体はボクセルと呼び、各ボクセルには、通し番号が与えられているとする。そして、ｉ番目のミュオンの軌跡上のｊ番目のボクセルにおける散乱角度をθｉｊ、ｊ番目のボクセルを通過したミュオン数をＮｊ、ｊ番目のボクセルを通過したミュオンの番号の集合をＳｊとすると、ｊ番目のボクセルのミュオンの散乱角度θｊの分布の分散σθｊ２は、数式（５）で求めることができる。なお、散乱角度θｉｊは、ミュオンのベクトル対及び散乱点から特定できる。また、ミュオン数Ｎｊは、ミュオンのベクトル対及び散乱点から特定してもよいし、ミュオンの軌跡から特定してもよい。
【００６７】
【数５】
2016180664.tif 000007
【００６８】
ここで、ＰｏＣＡでは、ミュオンが１点で散乱したと仮定するため、ｉ番目のミュオンの散乱角度をθｉとすれば、ｊ番目のボクセルがｉ番目のミュオンの散乱点を含む場合、θｉｊ＝θｉとなり、ｊ番目のボクセルがｉ番目のミュオンの散乱点を含まない場合、θｉｊ＝０となる。
【００６９】
なお、第１生成部２３において、ミュオン（疑似荷電粒子）に重みが与えられた場合、ｉ番目のミュオンの重みをｗｉとすれば、分散σθｊ２は、数式（６）に示す重み付き平均で求めることができる。
【００７０】
【数６】
2016180664.tif 000008
【００７１】
また、画像再構成法として、期待値最大化アルゴリズムやエネルギー最小化に基づく手法を用いる場合、尤度関数における各ミュオンの散乱の対数尤度の和で表される項をｗｉで重み付けすればよい。
【００７２】
また、検出部１０Ａ、１０Ｂにおいて、ミュオン（実荷電粒子）の運動量を検出でき、実荷電粒子データにミュオンの運動量が含まれる場合、数式（５）、（６）に代えて、それぞれ、数式（７）、（８）を用いることができる。
【００７３】
【数７】
2016180664.tif 000009
【００７４】
【数８】
2016180664.tif 000010
【００７５】
ここで、ｐｉは、ｉ番目のミュオンの運動量を示す。
【００７６】
図８は、本実施形態の画像処理部２０で行われる処理の手順の流れの一例を示すフローチャートである。
【００７７】
まず、推定部２１は、検出部１０Ａ、１０Ｂにより検出された複数の実荷電粒子それぞれ毎に、検出部１０Ａ、１０Ｂから入力された観測データに含まれる当該実荷電粒子の実荷電粒子データに含まれる実ベクトル対を用いて、当該実荷電粒子の実散乱点を推定する（ステップＳ１０１）。
【００７８】
続いて、第１生成部２３は、推定部２１により推定された実散乱点それぞれ毎に、当該実散乱点の近傍領域に１以上の疑似散乱点を生成する（ステップＳ１０３）。
【００７９】
続いて、第２生成部２５は、第１生成部２３により生成された疑似散乱点それぞれ毎に、検出部１０Ａ、１０Ｂから入力された観測データに含まれる、当該疑似散乱点が生成された実荷電粒子の実荷電粒子データを用いて、当該疑似散乱点を通る疑似荷電粒子の入射方向を示す疑似入射ベクトルと射出方向を示す疑似射出ベクトルとの組である疑似ベクトル対を生成する（ステップＳ１０５）。
【００８０】
続いて、画像生成部２７は、検出部１０Ａ、１０Ｂから入力された観測データに含まれる複数の実荷電粒子の実荷電粒子データに含まれる実ベクトル対、推定部２１により推定された当該複数の実荷電粒子の実散乱点、第１生成部２３により複数の実荷電粒子それぞれに対し生成された１以上の疑似散乱点、及び第２生成部２５により疑似散乱点それぞれに対し生成された疑似ベクトル対を用いて、当該複数の実荷電粒子が透過した対象物５の内部構造を示す画像を生成する（ステップＳ１０７）。
【００８１】
以上のように本実施形態によれば、実荷電粒子の実散乱点の近傍に擬似荷電粒子の擬似散乱点を生成するため、推定精度の向上に寄与できるように擬似散乱点を生成することができる。
【００８２】
特に本実施形態では、複数の実荷電粒子それぞれの実散乱点及び実ベクトル対に加え、各実荷電粒子に対し、上記のように生成された１以上の擬似散乱点、及び当該１以上の擬似散乱点それぞれに対し生成された擬似ベクトル対を用いて、当該複数の実荷電粒子が透過した対象物の内部構造を示す画像を生成する。このため本実施形態によれば、生成された画像におけるぼやけなどを抑えることができ、当該画像における対象物の内部構造の可視化精度を向上できる。この結果、例えば、ぼやけなどが原因で検出対象の物質を画像上で判別できないなどという事態の発生を低減させることができる。
【００８３】
また本実施形態では、荷電粒子にミュオンをはじめとする透過力の高い自然放射線を用いることで、人工の放射線源を用いずに非破壊で対象物の内部を可視化できる。
【００８４】
（ハードウェア構成）
図９は、本実施形態の画像処理部２０のハードウェア構成の一例を示す図である。図９に示すように、画像処理部２０は、ＣＰＵなどの制御装置９０１と、ＲＯＭ（Read  Only  Memory）やＲＡＭ（Random  Access  Memory）などの主記憶装置９０２と、ＨＤＤ（Hard  Disk  Drive）やＳＳＤ（Solid  State  Drive）などの補助記憶装置９０３と、ディスプレイなどの表示装置９０４と、キーボードやマウスなどの入力装置９０５と、通信インタフェースなどの通信装置９０６と、を備えており、通常のコンピュータを利用したハードウェア構成となっている。
【００８５】
画像処理部２０で実行されるプログラムは、インストール可能な形式又は実行可能な形式のファイルでＣＤ−ＲＯＭ、ＣＤ−Ｒ、メモリカード、ＤＶＤ（Digital  Versatile  Disk）、フレキシブルディスク（ＦＤ）等のコンピュータで読み取り可能な記憶媒体に記憶されて提供される。
【００８６】
また、画像処理部２０で実行されるプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するようにしてもよい。また、画像処理部２０で実行されるプログラムを、インターネット等のネットワーク経由で提供または配布するようにしてもよい。また、画像処理部２０で実行されるプログラムを、ＲＯＭ等に予め組み込んで提供するようにしてもよい。
【００８７】
画像処理部２０で実行されるプログラムは、上述した各部をコンピュータ上で実現させるためのモジュール構成となっている。実際のハードウェアとしては、ＣＰＵがＲＯＭやＨＤＤなどからプログラムをＲＡＭ上に読み出して実行することにより、上記各部がコンピュータ上で実現されるようになっている。また、画像処理部２０外部から入力される観測データなどの各種データをＨＤＤなどに記憶しておき、上記各部がＨＤＤなどから観測データなどをＲＡＭ上に読み出して上述した処理を実行する。
【００８８】
なお、上記各部の少なくとも一部をＩＣなどのハードウェアで実現する場合、画像処理部２０は、当該ＩＣなどを更に備えればよい。
【００８９】
なお、本発明は、上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化することができる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成することができる。例えば、実施形態に示される全構成要素からいくつかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせても良い。
【００９０】
例えば、上記実施形態のフローチャートにおける各ステップを、その性質に反しない限り、実行順序を変更し、複数同時に実施し、あるいは実施毎に異なった順序で実施してもよい。
【００９１】
以上のように、上記実施形態によれば、推定精度の向上に寄与できるように擬似的な荷電粒子の散乱点を生成することができる。
【符号の説明】
【００９２】
１  画像処理システム
１０Ａ、１０Ｂ  検出部
２０  画像処理部
２１  推定部
２３  第１生成部
２５  第２生成部
２７  画像生成部
３０  表示部
９０１  制御装置
９０２  主記憶装置
９０３  補助記憶装置
９０４  表示装置
９０５  入力装置
９０６  通信装置
【図１】
2016180664.tif 000011
【図２】
2016180664.tif 000012
【図３】
2016180664.tif 000013
【図４】
2016180664.tif 000014
【図５】
2016180664.tif 000015
【図６】
2016180664.tif 000016
【図７】
2016180664.tif 000017
【図８】
2016180664.tif 000018
【図９】
2016180664.tif 000019