(57)【要約】【課題】本発明は時間的に連続したディジタル医用画像の間の経時変化を検知するために強調した画像を用いて診断の精度を高めるためのコンピュータ化された方法及び装置を提供することを目的とする。【解決手段】１組の画像はディジタル化（ステップ１０、２０）され、２枚の画像における対応する位置が等しくなるように一方の画像を非線形に歪ませることを含む画像位置合わせが行なわれる（ステップ４０）。この後、歪みを受けた画像と他方の画像がサブトラクション処理される（ステップ５０）。後から撮影した画像にだけ存在する僅かに不透明になった部分が、サブトラクション画像において検知される。【選択図】      図１
2005176402.tif 000002
【特許請求の範囲】
【請求項１】
時間的に連続するデジタル化された画像についてのコンピュータ分析のためのシステムであって、
分析対象についての第１および第２のデジタル化された画像を得る手段と、
前記第１および第２画像内に存在する対応する構造的特徴に関連させ、前記第１および第２画像の対応する部分の位置を調整する手段と、
ここで、前記位置を調整する手段は、前記第１および第２画像の対応する画素の位置に関して、前記第１および第２画像の対応する構造的特徴間の位置の差を示す複数のシフト値を生成するために相互相関マッピングを実行する手段と、
前記相互相関マッピングに導入された誤差を補償するために前記シフト値を訂正するため所定の面を前記シフト値に数学的にフィッティングする手段と、
そして訂正されたシフト値に基づき前記第１および第２画像のうちの一方の構造的特徴をシフトする手段とを含み、
前記第１および第２画像間でサブトラクション処理を実行する手段とを含み、
それにより、前記第１および第２画像を得た間の時間経過の間に、前記分析対象について生じた構造的変化を示す前記第１および第２画像間の差が決定されるコンピュータ分析のためのシステム。
【請求項２】
ここで、前記相互相関マッピングを実行する手段は、
前記第１および第２画像の一方のテンプレートＲＯＩと他方の画像の探索領域ＲＯＩを選択する手段と、ここで前記探索領域ＲＯＩは前記テンプレートＲＯＩより大きく、
前記一方の画像の前記テンプレートＲＯＩと、前記他方の画像の前記探索領域ＲＯＩ内の異なる位置の複数のサブ領域との間の相互相関値を決定する手段と、
前記一方の領域の前記テンプレートＲＯＩに対し、前記他方の画像の前記探索領域ＲＯＩのどのサブ領域が最大の相互相関値となるかを決定する手段と、そして
前記一方の画像のテンプレート領域の位置と、前記テンプレートＲＯＩに関し最大の相互相関値となる前記探索領域ＲＯＩの前記サブ領域の他方の画像の位置との間の位置の差に基づき、前記テンプレートＲＯＩに関する前記シフト値を決定する手段と、
を含む請求項１に記載のシステム。
【請求項３】
ここで、前記所定の面をフィッティングする前記手段は、
０〜１の範囲の重み付け係数を使用する前記シフト値の前記相互相関マッピングに対し、２次元ｎ次の多項式関数を使用する手段を含む請求項２に記載のシステム。
【請求項４】
ここで、前記所定の面をフィッティングする前記手段は、
最大相互相関値に基づく前記重み付け係数を選択する手段を含む請求項３に記載のシステム。
【請求項５】
ここで、２次元ｎ次の多項式関数を使用する前記手段は、
前記相互相関マッピングにより得られた第１および第２画像の第１の相関部分に対する重み付け係数として１を、および第１および第２画像の前記第１の相関部分より低い相関性を有する第２の相関部分に対する重み付け係数として０を選択する手段を含む請求項３に記載のシステム。
【請求項６】
さらに、前記重み付け係数を生成するためのルックアップテーブル手段を含む請求項３に記載のシステム。
【請求項７】
さらに、前記サブトラクション処理により得られたサブトラクション画像の出力画像表示を表示する手段を含む請求項１に記載のシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、例えば胸部Ｘ線画像のようなディジタル医用画像のコンピュータ解析の分野に関し、２枚の画像の間の違いがディジタル画像処理技術を用いて検出することができるように１枚を他方に先だって得た一人の患者についての１組の画像の間に生じる微妙な経時変化を強調するために用いられるコンピュータ解析に関する。
【背景技術】
【０００２】
臨床医療では、同じ患者の過去に撮影した胸部Ｘ線画像との比較は重要であると考えられてきた。ほとんどの後方−前方(POSTERO-ANTERIOR:Ｐ−Ａ) 胸部Ｘ線画像は過去に撮影したＸ線画像と並べて観察される。比較しながらの読影は、読影者が異常を発見し、その程度を決めることの助けとなる。また、比較しながらの読影は、既知の病巣の時間変化を検知し、処置の効果を評価するためにも重要である。しかし、Greene, "Missed Lung Nodules…", Radiology, vol.182, pp.1, 8-9, 1992と、Austin, "Missed Bronchogenic Carcinoma …", Radiology,vol.182 , p.1, 115-122に記述されているように、連続した胸部Ｘ線画像を観察しているときでさえ、読影者は重要な経時変化を見落としてしまうことがある。このような誤りが生じるのは、部分的には２枚のＸ線画像を交互に見て比較することが難しいことや、２枚のＸ線画像の間で濃度、コントラスト、または患者の姿勢が違っていること等に原因がある。数多くの異常がある患者の経時変化は見落とされることがある。というのは、いくつかの異常は変化を示さない他の異常でカムフラージュされてしまうからである。したがって、読影者に注意を促すとともに病理学上重要な変化を検知する率を高めるために、時間的に連続した胸部Ｘ線画像における経時変化に関連した情報を得ることは重要なことである。経時変化とは、ここでは前回の検査の後で今回の検査の前に生じた病理学的変化と定義している。
【０００３】
一般に、患者の姿勢、Ｘ線の照射の状態やその他の露光条件等に関してＸ線画像を再現することは困難である。また、患者の呼吸や心拍は、通常２枚の画像の間では異なった状態にあり、そのため、肺、横隔膜、心臓の大きさや形が違ってしまう。そのため、１枚の画像の画素が他方の画像で全く同一の解剖学的な構造に対応する画素から差し引かれるような画像の位置合わせ技術が必要である。さらに、３次元の対象物を２次元に投影していることになるＸ線画像は、単純な方法で正確に位置合わせをすることが困難である。図１１は時間的に連続した複数組のＰ−Ａ胸部画像での位置ずれの様々な原因を示している。これらは４０組以上の画像を調べて得られたものである。ほとんどの組で、これらの原因の組み合わせによって位置ずれが起こっている。
【発明の開示】
【発明が解決しようとする課題】
【０００４】
本発明は上述した事情に対処すべくなされたもので、その目的は、時間的に連続したディジタル医用画像の間の経時変化を検知するために強調した画像を用いて診断の精度を高めるためのコンピュータ化された方法及び装置を提供することである。
本発明の関連する目的は、１枚の画像を非線形に歪ませた後に、歪ませた画像と他方の画像の差分をとることによって時間的に連続した１組の胸部画像の位置合わせ処理を行うためのコンピュータ化された方法及び装置を提供することである。
【０００５】
さらなる目的は、患者の時間的に連続した１組の画像のうちの１枚を非線形に歪ませることが２枚の画像での多くの小さな関心領域(region of interest:ＲＯＩ) に適用されるローカルマッチング技術に基づいている画像の位置合わせとサブトラクションを行なうためのコンピュータ化された方法及び装置を提供することである。
【０００６】
本発明のさらなる目的は、時間的に連続した１組の胸部画像をマッチングさせ、２枚の画像間のサブトラクションを行ない、２枚の胸部画像の画像解析によって得られたデータに基づいて画像上の多くの小さなＲＯＩを選択することを含む画像位置合わせ技術を用いて肺疾患を強調し検知するためのコンピュータ化された方法及び装置を提供することである。
【０００７】
本発明の別の目的は、時間的に連続した２枚の画像上の多くの小さなＲＯＩにローカルマッチング技術を適用した結果を示すマッチングシフト値を決定することに基づいて、２枚の画像の間の位置合わせを行うために非線形歪み技術を用いた時間的に連続した１組の医用画像の間の経時変化を検知するためのコンピュータ化された方法及び装置を提供することである。
【０００８】
本発明の別の目的は、画像の中で異常の疑いのある領域を強調し検知するためにディジタル化された１組の画像のサブトラクションを行ない、ディジタル画像のサブトラクション処理を可能とするために２枚の画像のうちの一方に適用される非線形歪み技術に基づいて２枚の画像が互いにマッチングされ、その結果、非線形歪み技術は２枚の画像の間の対応する場所の間のシフト値のマッピングを用いて行なわれ、シフト値のマッピングはシフト値にフィットさせた２次元の多項式関数に基づいている方法及び装置を提供することである。
【０００９】
本発明の別の目的は、重み係数が各シフト値の多項式フィッティングへの寄与の程度を表している多項式関数を複数のシフト値にフィッティングさせる重み付けフィッティング技術を用いる画像の位置合わせとサブトラクションを行なうための方法及び装置を提供することである。
【００１０】
本発明の別の目的は、フィッティング技術はローカルマッチング処理に関する情報、および／または画像解析結果に基づいて決定された重み係数を有する重み付け技術であり、２枚の画像で対応する位置の間の位置ズレの量を示す複数のシフト値に多項式関数をフィットさせる重み付けフィッティング技術を用い、時間的に連続した１組の医用画像において異常を強調し検知することを可能にするコンピュータ化された方法を実行するための方法及び装置を提供することである。
【課題を解決するための手段】
【００１１】
上記のまたは他の目的を達成するために、本発明は、既知のＸ線技術を用いることによって撮影され、さらにディジタル化され、ディジタルコンピュータによって処理される１組の胸部画像が生成されるコンピュータ化された方法を含んでいる。本技術は、サブトラクション処理が行えるように複数の画像の中で対応する位置どうしを位置合わせする画像位置合わせ技術を含んでいる。この２枚の画像のサブトラクションを行なうことによって、間質性疾患や腫瘍等の異常を示していると考えられる２枚の画像の間の微妙な違いが強調される。２枚の画像の位置合わせを行うために、非線形歪みが２枚の画像の一方に適用されて、２枚の画像の対応する場所の間の必要な位置合せが得られる。この非線形歪み法は、胸部画像の解剖学的構造に関する検知された情報に基づいて選択された多くの小さな関心領域（ＲＯＩ）のローカルマッチング技術に基づいている。このローカルマッチング技術は、２枚の画像の対応する場所の間のシフト値をマッピングするものである。ここで、シフト値のマッピングはローカルマッチング量、またはＲＯＩに対して行われる画像データ解析の結果に基づく重み係数を用いるカーブフィッティング技術によってなされる。
【発明の効果】
【００１２】
本発明によるコンピュータ化された方法及び装置によれば、強調した画像を用いて時間的に連続したディジタル医用画像の間の経時変化を検知し、診断の精度を高めることができる。
【発明を実施するための最良の形態】
【００１３】
以下、図面を参照して本発明によるコンピュータ化された経時変化検出方法及び装置の実施例を説明する。図１は本発明に基づく基本的な概念がフローチャートとして図示している。この図を見てわかるように、本発明の方法は第１、第２の１組の画像のディジタル化（ステップ１０，２０）と、画像位置合わせ（ステップ４０）と、画像のサブトラクションを行なうこと（ステップ５０）とを含んでいる。ディジタルＸ線画像の代わりに、蓄積型の蛍光管を用いたシステムやＣＣＤを用いたシステムを使って直接得られるディジタル胸部画像を使用することができる。括弧付きで書かれているステップ３０と６０は行うことが望ましいが、本発明の方法を実行するのに不可欠なものではない。差分画像はステップ７０で出力される。
【００１４】
図２は時間的に連続した一組の胸部画像の間の経時変化を強調するときの全体像をより詳しく図示したものである。この方法は１組の胸部Ｘ線画像をディジタル化するステップ（ステップ１０，２０）と、多くの小さな関心領域（ＲＯＩ）を選択するステップ（ステップ４１，４２）と、選択されたＲＯＩを用いて画像の位置合わせをするステップからなる。画像の位置合わせは、２枚の画像の対応するＲＯＩの各組に対して適用されるローカルマッチング（ステップ４３）と、対応するＲＯＩの各組の間のローカルマッチングの量を示すシフト値のマッピング（ステップ４４）からなっている。カーブフィッティング（ステップ４５）がこのマッピングされたシフト値に適用され、このカーブフィッティング解析の結果に基づいて２枚の画像のうちの一方に非線形の歪みが与えられる（ステップ４６）。最後に対応するＲＯＩの各組に対して、歪まされた画像と歪んでいない画像との間のサブトラクションが行なわれる（ステップ５０）。このサブトラクション画像が出力され（ステップ７０）、読影者はこのサブトラクション処理によって得られた結果と時間的に連続した画像とを横に並べて比較することができる。ステップ３０，６０は実施した方が好ましいが本発明にとって本質的なものではない。
【００１５】
図３はふたつのディジタル化された２枚の画像内でＲＯＩを選択する過程（ステップ４１，４２）を示している。この過程は胸部画像の解剖学的構造の解析（ステップ１１、例えば胸郭の端の位置等）と、各ＲＯＩの中心間距離の決定（ステップ１４）と、その中心の位置の決定（ステップ１３）からなっている。ＲＯＩの大きさと形もまた選択される（ステップ１５）。このようにしてＲＯＩの場所が決定される（ステップ４７）。ステップ１２も必要に応じて行うステップであり、本発明に必要なものではない。
【００１６】
本発明に従ったローカルマッピング技術を用いたシフトマッピング技術は図４に示してある。この図からわかるように、シフトマップが決定され（ステップ８０）、胸部画像のテンプレートＲＯＩと探索領域ＲＯＩの間の類似性の解析（ステップ８１）、解剖学的構造の解析（ステップ８２）、画素値分布の解析（ステップ８３）の後で重み付けカーブフィッティング技術で用いられる重み係数が決定される（ステップ８４）。重み付けフィッティングが行われた後で（ステップ８５）シフトマッピング計算では直接計算できない値を推定するために補間（内挿／外挿）技術が適用される（ステップ８６）。次に、フィットされたシフトマップが作られる（ステップ８７）。ステップ８１〜８３は必要に応じて行えばよい。
【００１７】
ローカルマッチング技術は２枚の画像、すなわちテンプレートＲＯＩ画像と探索領域ＲＯＩ画像におけるＲＯＩの選択を含んでいる。このテンプレートＲＯＩはもうひとつの画像中の探索領域ＲＯＩに対して最もよく一致するように探索領域ＲＯＩの中を移動される。図５に示されているように、この探索領域ＲＯＩとテンプレートＲＯＩはコンピュータによって自動的に選択される。図５では、各画像におけるＲＯＩの中心が数多くの点で示されている。
［画像位置合わせ技術］
本発明で用いられている画像位置合わせ技術の基本的な概念は、サブトラクション処理が実行できるように２枚の画像の間の位置合わせを改善するために、２枚の画像のうちのひとつを非線形に歪ませる方法を用いている。非線形歪み処理は２枚の画像で対応している場所にある多くの小さな関心領域（ＲＯＩ）のローカルマッチングに基づいて行なわれる。ローカルマッチング処理を行うためには、ひとつの画像の中で小さなテンプレートＲＯＩが選択される。もうひとつの画像では探索領域が選択される。そこでは、テンプレートサブ画像に一致するサブ領域が探索領域の中で探索される。言い替えると、第１の画像のテンプレートサブ画像に最も類似するサブ領域はもうひとつの画像の中の探索領域の中に必ず存在していなければならないことになる。この探索領域を画像全体とすることは可能であるが、解剖学的構造から見て対応すると思われる位置にあるテンプレートＲＯＩよりやや大きい領域を探索領域とすることが好ましい。これは精度を高め、また計算時間を短くするためである。
【００１８】
ＲＯＩは胸郭を含む肺領域で主に選択することが好ましい。何故ならば、肺のこの領域に最も重要な診断情報があるからである。胸郭の端や横隔膜、心臓領域、縦隔の場所等の解剖学的構造の特徴を検知することができる既知の画像解析技術によって肺領域は特定することが可能である。こういった既知の技術は、例えば米国特許第４，８５１，９５４号に述べられている。
【００１９】
テンプレートと探索領域の中心位置も上述した画像解析技術に基づいて決定することができる。探索領域の中心は、対応するテンプレートの中心と解剖学的に見てほぼ同じ場所に位置しているべきである。テンプレートＲＯＩの中心と探索領域ＲＯＩの中心はそれぞれの画像に独立に上述した画像解析技術を用いて選択することができる。別のやり方として、同じ方法で一方の画像でテンプレートＲＯＩを選択し、胸郭の端や肺の中心線や肺の上端のような所定の解剖学的構造の相対的な位置の情報に基づいてテンプレートの中心点を全体的にシフトさせたり、及び／または回転させたりして他方の画像で探索領域の中心の場所を決めることもできる。この方法は肺領域のある断片的な部分に独立に適用するとさらに正確になる。例えば、図５に示されているように中心線によって分けられた右肺領域と左肺領域に独立に適用した場合である。図５はテンプレートの中心が胸郭領域の中に、肺の中心線と上辺の交点の相対的な位置に基づいてテンプレートの中心をシフトさせることによって得られた対応する探索領域の中心として選択されているところを図示している。
【００２０】
テンプレートＲＯＩの大きさは好ましくは直径５〜５０（ｍｍ）であり、さらに好ましくは直径１５〜４０（ｍｍ）である。ここで選択される大きさは使用されるローカルマッチングアルゴリズムの特性にもよる。図６、図７に示されているように、テンプレートＲＯＩのサイズが小さすぎると、マッチングの正確さが低下してしまう。一方、サイズが大きすぎると処理時間が長くなってしまう。また、サイズが大きすぎると、ＲＯＩ全体に対するマッチングはＲＯＩの中心部分でのマッチングと同じではなくなってしまう。言い替えると、ひとつの大きなＲＯＩがテンプレートＲＯＩと探索領域ＲＯＩのそれぞれに用いられたとすると、広い領域に対して解析が行われるためにＲＯＩの中の構造が変わってしまう。したがって、そのような２枚の領域でマッチングを行なおうとすると、２つのＲＯＩのある部分はよく重なるが、他の部分がマッチングしないことがおきる可能性がある。探索領域ＲＯＩのサイズは対応するテンプレートＲＯＩより大きくすべきであり、好ましくは直径１０〜４０（ｍｍ）である。図６、図７では検索領域サイズはＲＯＩサイズ＋２２．４とされている。探索領域ＲＯＩのサイズもまた検知されるべき経時変化の大きさに依存している。探索領域ＲＯＩのサイズが小さすぎると、マッチングの正確さが低下する。一方、サイズが大きすぎると、対応する領域はマッチするが、心臓の大きさの変化や横隔膜の変化等の病理学的な経時変化が減少させられてしまう。２枚の領域の中心間の距離は好ましくは５〜３０（ｍｍ）である。しかし、この距離が小さすぎると、処理時間が相対的に長くなる。一方、この距離が大きすぎると、歪み処理の正確さが低下してしまう。テンプレートと探索領域ＲＯＩの形は図８に示されている多くのバリエーションのうちのいずれのものでも良い。
【００２１】
ＲＯＩは胸郭全体を覆うように位置している必要がないことに気をつけるべきである。例えば、縦隔領域の非常に低濃度の部分では肺領域に較べてＲＯＩは少なくて良い。何故ならば、その部分は診断では重要ではないからである。縦隔の領域では中心間の距離を長くしたり、縦隔領域の中心部分の低濃度領域ではＲＯＩを選択しないようにすることができる。縦隔は胸部画像の中心近くに見られる白い領域であり、これは脊椎及び関連する解剖学的構造のためである。したがってこの領域の濃度は低く、そのためこの領域からは多くの情報は得られず、ここでは少ない数のＲＯＩしか選択されない。
【００２２】
本発明のローカルマッチング技術は、対応するテンプレートＲＯＩに「最も良く」マッチする探索領域ＲＯＩの中のサブ領域を決定するために用いられる。「最も良く」マッチするということは、決定されたサブ領域の中の画像データとテンプレートの中の画像データの間の類似性が探索領域の中の他の可能なサブ領域よりも高いという意味である。２枚のサブ画像の間の類似性を評価し、「最も良く」マッチすることを決定できる、現在知られているいずれのアルゴリズムもローカルマッチングのために用いることができる。以下の式と図９は正規化相互相関法を示したものである。
【００２３】
【数１】
2005176402.tif 000003
ここで、ＡとＢが等しい場合は、Ｃ＝１である。
上式は正規化相互相関値を最大にするサブ領域を見つけるために使用される。Ｈａｌｌ著「コンピュータ画像処理と認識(Computer Image Processing and Recognition) 」アカデミックプレス(Academic Press)、１９７９年、として知られている正規化しない一般的な相互相関法も使用することができるが、２枚のサブ画像間の平均画素値の違いの影響を低減させるためには正規化する方が好ましい。
【００２４】
本発明に従ってローカルマッチング処理を行うことができる第２の方法は最小二乗差分法を用いることである。この方法によれば、次式に基づいて２枚の画像の対応する位置の間の差の二乗和を最小にするサブ領域が決定される。
【００２５】
【数２】
2005176402.tif 000004
この式は正規化相互相関法を用いて得られる結果と同じ結果を与える。
本発明に従った２枚の画像位置の間でローカルマッチングを行う第３の方法は最小差分法である。この技術によれば、以下の式を用いて２枚の画像で対応する各場所（画素）の間の差の絶対値の総和を最小にするサブ領域が決定される。
【００２６】
【数３】
2005176402.tif 000005
本発明に従った第４のローカルマッチング法はフーリエ変換位相相関技術を用いるものであり、２枚のサブ画像のフーリエ変換の位相項の間の相関を最大にするサブ画像が以下の式にしたがって決定される。
【００２７】
【数４】
2005176402.tif 000006
正規化相互相関法、最小二乗法差分法、及びフーリエ変換位相相関法ではマッチングの正確さに関して最も良い結果が得られる。上記のローカルマッチング技術を用いて良い結果を得ることができたパラメータの例は次のようなものである。画素の大きさは０．７ｍｍ、テンプレートＲＯＩのサイズは２２．４ｍｍ×２２．４ｍｍ（３２×３２画素）、探索領域のサイズは４４．８ｍｍ×４４．８ｍｍ（６４×６４画素）、中心間距離は７．０ｍｍ（１０画素）、正規化相互相関法を使用している。
【００２８】
ローカルマッチングアルゴリズムを用いた計算は画素毎に行う必要がないということにも注意すべきである。つまり、処理時間を節約するためにテンプレートＲＯＩと探索領域ＲＯＩの両方で画素を飛ばすことができるということである。上記のローカルマッチング技術を用いて良い結果を得ることができたパラメータの例は次のようなものである。画素の大きさは０．７ｍｍ、テンプレートＲＯＩは２２．４ｍｍ×２２．４ｍｍ（３２×３２画素）、探索領域の大きさは４４．８ｍｍ×４４．８ｍｍ（６４×６４画素）、中心間距離は７．０ｍｍ（１０画素）、正規化相互相関法を使用している。計算は２画素毎に１画素を飛ばし各テンプレートについて１６×１６画素で行なわれた。計算の結果得られる位置合わせは画素を飛ばさない場合に得られる結果とほぼ同じであり、処理時間は画素を飛ばさない場合の１／４に減少した。
【００２９】
上述のローカルマッチング技術のうちのひとつをテンプレートの各中心に対して２あるいはそれ以上のステップで行うことができる。例えば、第１のステップとして、比較的大きなテンプレートＲＯＩと比較的大きな探索領域ＲＯＩを選択することができる。次に、「最も良く」マッチするサブ領域の中心のおおよその位置を見つけるために画素を飛ばして粗いローカルマッチング法を行なうことができる。次に、おおよその「最も良く」マッチする場所近くで小さなテンプレートＲＯＩと探索領域ＲＯＩが、ローカルマッチングの第２のステップとして厳密なローカルマッチングを行うために選択される。適切なパラメータが選択されたならば、この方法は正確さを高め、また処理時間を節約する。
【００３０】
上記の方法に使用されたパラメータの例は以下の通りである。画素の大きさは０．７ｍｍ、テンプレートＲＯＩは第１のステップでは４４．８ｍｍ×４４．８ｍｍ（６４×６４画素）で、第２のステップでは１１．２ｍｍ×１１．２ｍｍ（１６×１６画素）、探索領域の大きさは第１のステップでは６７．２ｍｍ×６７．２ｍｍ（９６×９６画素）で、第２のステップでは１９．６ｍｍ×１９．６ｍｍ（２８×２８画素）、中心間距離は７．０ｍｍ（１０画素）、正規化相互相関法を使用している。第１のステップでは４画素毎に１画素を計算するように画素を飛ばして各テンプレートについて各１６×１６画素に対して計算が行われた。この結果得られた位置合わせは上述の結果とほぼ同じであり、画素を飛ばさない場合に較べて計算時間は１／３に減少した。
【００３１】
テンプレートＲＯＩと探索領域ＲＯＩにローカルマッチング技術を適用することにより、「最も良く」マッチするサブ領域の中心の場所のカーテシアン座標を見つけることができる。これは探索領域で（ｘ′，ｙ′）で表される。図９に示されているように、（ｘ′，ｙ′）とテンプレート（ｘ，ｙ）の中心の場所の間のｘ方向とｙ方向の間の距離はそれぞれシフト値Δｘ，Δｙと呼ばれる。多くのテンプレートと探索領域の組から得られた胸部画像の重要な領域上のこれらのシフト値の分布は、画像のサブトラクション処理に先だって座標（ｘ，ｙ）の適切な幾何学的歪みを得るために必要となる重要な情報である。
【００３２】
図１０はシフト値ΔｘとΔｙの分布を示すシフトマップの例を示しており、これは正規化相互相関法を用いたローカルマッチングによって得られたものである。ここで、画素の大きさは０．７ｍｍ、テンプレートの形状は２２．４ｍｍ×２２．４ｍｍ（３２×３２画素）の正方形となるように選択され、探索領域のサイズは４４．８ｍｍ×４４．８ｍｍ（６４×６４画素）となるように選択され、中心間距離は７．０ｍｍ（１０画素）となるように選択された。
【００３３】
（ｘ，ｙ）座標を適正に歪ませるためには、上述のシフトマップを補間またはフィットさせ、胸部画像の各（ｘ，ｙ）座標に対応するΔｘとΔｙの値を得ることが必要である。これはΔｘまたはΔｙの値を補間することによって行うことができる。最隣接補間や線形補間、スプライン補間等のいずれの既存の補間法もこの目的のために用いることができる。しかし、補間の結果は、いくつかの「悪い」マッチングのテンプレートを含む全ての結果に依存しており、また本当に経時変化した部分にも依存しており、この部分ではいずれにしろマッチングはできない。言い替えると、腫瘍やその他の間質性疾患の進行等の本当の経時変化があったとき、胸部画像に本当の変化が現れることになる。この場合、２枚の画像の特徴が違うため２枚の画像は完全にはマッチングしない。ローカルマッチングの理由は「マッチする」領域を可能な限り特定することであり、その後お互いをサブトラクションすることである。一方または両方の画像の低濃度領域、付加的なアーチファクト、あるいは疾患の進行による本当の変化が存在することにより、マッチングがうまくとれない領域がでてくる。例えば、図１０では、ほとんどの領域では比較的滑らかであるにも関わらず、ΔｘとΔｙの分布でシフト値が異常に大きいか小さい特異点がいくつか見られる。これらの特異点は一般に経時変化を起こして組織、位置が変わってしまった解剖学的構造や器官、または心臓領域内の著しく濃度が低い領域に対応している。これらの特異点では、これらのテンプレートに対する「良い」マッチングを見いだすことができないので、シフト値は信頼できるものではなくなる。したがって、補間法が用いられたとすると、これらのテンプレート近くの補間によって得られたシフト値は正しくないシフト値になりがちである。
【００３４】
この影響を取り除くために、シフト値にカーブフィッティング技術が適用される。どのようなタイプのカーブフィッティング法もここで用いることができるが、米国特許第４，８５１，９８４号で知られるように、２次元のｎ次多項式関数（ｎ≧２）によるフィッティングが滑らかさの点から好ましい。２次元のｎ次多項式関数の一般的な形は次のようなものである。
【００３５】
【数５】
2005176402.tif 000007
以下の式がシフト値へのカーブフィッティングを行うために用いることができる。
【００３６】
【数６】
2005176402.tif 000008
ここで、ｋ＝１，２，…，ｎ＋１
ｍ＝１，２，…，ｋ
正確さを求める時には、３次以上の高い次数の多項式が好ましい。何故ならば、図１９、図２０に示されている原因によって生じるＸ線画像の位置ずれの本質的な性質から考えられるように、ｘ，ｙ座標を適正に歪ませることはかなり複雑な変換であると考えられるからである。
【００３７】
２枚の画像の位置合わせを向上させるためには、２次元のｎ次多項式を用いた重み付けフィッティングが好ましい。そのような重み付けフィッティングによって、ひとつの画像領域の中の各点について０〜１の範囲の重み係数が得られる。この重み係数は上の式（７）と（８）での係数ａi （ｉ＝１，２，…）を決定する際の関与の程度を表している。良い（信頼できる）マッチングを示しているシフト値には大きな重み係数を使用し、悪い（信頼できない）マッチングを示しているシフト値には小さい重み係数を用いるのが好ましい。この重み係数は、解剖学的構造の画像解析によるのと同様にローカルマッチング技術を用いて２枚の画像の間の類似性の解析、及び以下に述べる画素値の分布の解析に基づいて決定することができる。
【００３８】
ローカルマッチング法は、「最も良い」マッチングを示すサブ領域を決定するためにテンプレートＲＯＩと探索領域ＲＯＩの中で選択されたサブ領域との間の類似性を評価する。いずれのローカルマッチング法も相互相関法における正規化相互相関値のような類似性の尺度を持っている。この場合、前述の式（１）に示されているように、正規化相互相関値が大きくなればなるほど、テンプレートに対するマッチングは良くなる。したがって、重み係数は、各テンプレートに対して「最も良く」マッチングしているサブ領域から得られる類似性の尺度に基づいて決定することができる。このとき、大きな類似性を有する良くマッチングしているシフト値は大きな重み係数を有するようにする。例えば、相互相関法がローカルマッチングに使用されたとすると、重み係数と正規化相互相関値の間の関係を定義づけるルックアップテーブルを作成することができる。図１１に示されているように、シフトマップの中の全ての重み係数は、したがってこのルックアップテーブルを使って求めることができる。ＰとＱの可能な値の例としては、Ｐは典型的には約０．６５の値であり、Ｑは約０．９の値である。
【００３９】
肺領域は他の領域に較べてより重要な診断情報を含んでいるので、肺領域のテンプレートから得られたシフト値のフィッティングには大きな重み係数が与えられる。肺領域の場所は、米国特許第４，８５１，９８４号、第５，０７２，３８４号で公開されている胸郭の端や横隔膜、心臓領域、縦隔等の特定の解剖学的構造を検知するための画像解析技術によって特定することができる。例えば、重み係数は横隔膜の下や縦隔の中央では０にされ、他の領域では１にされる。また、この方法が胸部画像における特定の領域での経時変化を調べるために用いられる場合は、それらの特定の領域には大きな重み係数が与えられる。例えば、縦隔、心臓、横隔膜周辺といった様々な場所の濃度の分布はしばしば非常に低く、そのためこれらの領域からはほとんど情報は得られない。したがって、これらの領域からの寄与は無視され、重み係数も０にされる。ヒストグラム解析からは、これらの領域についての最大と最小の画素値がこれらの場所を示す。低い画素値のおよそ１０％は無視される。これらの領域は濃度が低いために、読影者に対してほとんどあるいは全く診断情報を与えないからである。また、重み係数は解剖学的構造の位置にも依存している。胸郭の端の外側ではほとんど診断情報は無く、この領域の外側でも相互相関、すなわちマッチングは無視される。
【００４０】
テンプレートから得られるシフト値のフィッティングのための大きな重み係数は診断上重要な領域から画素値の分布を解析することによって得られる。画像のディジタル化で高濃度の領域は大きな画素値へ変換されるならば、一般的には肺領域は比較的大きな画素値を持ち、一方、縦隔の領域は小さい画素値を持つ。それゆえ、縦隔領域の中央に対応する極めて低い平均画素値を持つテンプレートから得られたシフト値に対する重み係数は小さくするのが望ましい。胸部画像の中の異なった場所に対する重み係数を決定するために、重み係数と画素値の間の関係を定義するルックアップテーブルを作ることができる。例えば、重み係数は最小値と｛最小値＋０．１×（最大値−最小値）｝という画素値の間は０にされ、それより大きい画素値に対しては１とする。重み係数はこれら３つのタイプの解析の組み合わせに基づいて決定することができる。
【００４１】
図１２は図１０に示されたシフトマップから生成されたフィッティングされたシフト値ΔｘとΔｙの分布を示しているシフトマップである。重み付けフィッティングは、Ｐ＝０．６５，Ｑ＝０．９５とした図１１に示されたルックアップテーブルを使って決定された重み係数を用いた２次元の１０次多項式によってなされた。胸郭の外側の点のシフト値は線形補間（外挿）法によって決定された。
【００４２】
画素（ｘ，ｙ）の場所の関数としてシフト値ΔｘとΔｙが得られた後、その（ｘ，ｙ）座標の非線形変換を決定することができる。図１３、図１４に示されているように、得られた（Δｘ，Δｙ）の組はひとつの画像における座標（ｘ，ｙ）がもうひとつの画像における座標（ｘ＋Δｘ，ｙ＋Δｙ）に対応していることを示している。
【００４３】
図１３は本発明に基づく非線形歪み処理を示す。左側の点線の格子はひとつの画像のカーテシアン座標（ｘ，ｙ）を示し、この点線の格子に重なって示されている実線の格子はもうひとつの画像の対応する点（ｘ＋Δｘ，ｙ＋Δｙ）を示している。左側の実線の格子は適正な位置合わせを行うために右側に実線で示した格子に変形させられる。下に示した式はこの非線形な変換の一般的な形態を示している。
【００４４】
【数７】
2005176402.tif 000009
図１４は図１２のシフト値ΔｘとΔｙの組によって作られた（ｘ＋Δｘ，ｙ＋Δｙ）の格子パターンの例を示している。
上の式の係数ａiiとｂijは決定されていなくても良いということに注意が必要である。何故ならば、サブトラクション画像は、計算されたΔｘとΔｙの値を使ってひとつの画像の点（ｘ，ｙ）の画素値をもうひとつの画像の対応する点（ｘ＋Δｘ，ｙ＋Δｙ）の画素値から減算し、サブトラクション画像における点（ｘ，ｙ）の画素値としてサブトラクション画素値を得ることにより生成されるからである。フィッティングされたシフト値ΔｘとΔｙは実数であるので、点（ｘ＋Δｘ，ｙ＋Δｙ）はディジタル化された画像で画素の間に位置させることができる。そのような場合、その画素値は隣接する画素を用いる補間画素値で代用できる。最隣接補間や線形補間、スプライン補間等のいずれの既存の補間法もこの目的のために用いることができる。
【００４５】
本発明に基づいた全体のシステムは図１５、図１６に示されている。これらの図に見られるように、ディジタル画像入力装置１００は同じディジタル画像をメモリ１１０、１２０に供給する。これらのメモリの各々はＲＯＩ位置選択装置１３０によって受け取られる第１の出力を有している。第１画像メモリ１１０は歪み画像メモリ２００によって受け取られる第２の出力を有している。第２画像メモリ１２０はサブトラクション演算装置１２０によって受け取られる第２の出力を有している。ＲＯＩ位置選択装置１３０はＲＯＩの位置を選択した結果をＲＯＩ画像メモリ１４０とシフトマップ発生装置１６０に出力する。シフトマップ発生装置１６０はシフトマップ計算の結果をシフトマップメモリ１７０へ出力する。メモリ１７０には２枚の画像の間の値のマップが保存される。カーブフィッティング演算装置１８０はシフトマップメモリ１７０の出力を受け取り、上述のカーブフィッティングを行う。この結果はフィットシフトマップメモリ１９０に保存され、次に歪み画像メモリ２００に出力される。歪み画像メモリ２００に保存された歪まされた画像はサブトラクション処理を行うために、サブトラクション演算装置２１０に出力される。この計算の結果はサブトラクション画像メモリ２２０に保存され、続いて読影者が見られるように表示装置２３０に表示される。
【００４６】
図１７、図１８はカーブフィッティング技術を実行する際に用いられた重み係数を使った本発明に基づいたシステムの全体像を示している。この図に示されているように、画素値分布解析装置１３５は第１、第２画像メモリ１１０，１２０とＲＯＩ位置選択装置１３０の間に直列に接続されている。さらに付加してある類似度メモリ１５５はローカルマッチング演算装置１５０で実行されたローカルマッチング計算の結果を保存するために設けられている。これらの結果は、次に、選択されたＲＯＩに適用される適当な重み係数を発生するために重み係数発生装置１７５へ出力される。図示されているように、これらの重み係数は重み付けカーブフィッティング演算装置１８５へ出力され、そこで上述したある次数の多項式へのフィッティングが行われる。図１７、図１８の残りの部分は図１５、図１６に示されているものと同じであるのでここでは説明を省略する。
【００４７】
本発明の技術による効果は図２１〜図３４に示されている。図２１と図２２に示されているように、図２２の画像には空気の浸潤は存在していないが、２年後の画像では浸潤の存在がはっきりわかった。本発明の技術を使うことによって、２枚の画像は位置合わせされ、ディジタル的にサブトラクションされ、その結果患者の右肺の浸潤が強調された。
【００４８】
図２４〜図２７は歪まされた画像と歪まされていない画像、この２枚の画像のサブトラクション画像（右下）を図示している。図２８と図２９はある患者の１組の画像を図示したものであり、図２９は図２８より２年前に撮影された古い画像である。これらの画像の差画像（サブトラクション画像）は図３０に示されている。この図では、サブトラクション技術によって、２枚の画像の間で肥大した左上側の肺節のみが強調されて示されている。図３１は心臓肥大を示し、２年前に撮影された同じ画像を図３２に示している。図３３に示されているこれらのサブトラクション画像は心臓肥大の程度と肥大に伴う心臓の欠陥に由来するプロフュージョン(profusion) （画像の下に示されている）を明確に示している。図３４は１組の画像のうちのひとつを全体的にシフトさせて得られたサブトラクション画像を示している。この図では解剖学的に特徴のある部分の位置合わせがうまく行われないことによる様々なアーチファクトの存在が認められる。
【００４９】
図２４〜図２７は上述した画像位置合わせ処理法を用いて原画像のうちの１枚が歪まされている２枚の時間的に連続した原画像を用いた本発明で用いられているディジタル画像のサブトラクション処理の例を示している。処理の結果得られたサブトラクション画像も図示されている。図２４〜図２７の例では、最隣接補間法が使用された。
［前処理技術］
位置合わせ、及び／またはサブトラクション処理を改善するために、ＲＯＩの選択に先だって、またはローカルマッチングを行うのに先だって前処理を実行することができる。このような前処理技術の例は米国特許庁に継続中の米国特許出願第０７／７８６，００８号に開示されているような非線形濃度補正、マトリクスサイズリダクション、コントラスト強調、及び／またはエッジぼかし等である。前処理として非線形濃度補正を行うのが望ましい。何故ならば、この方法ではサブトラクション画像内で経時変化が起きなかった領域では均一な濃度が得られ、経時変化が生じた領域を強調するに有用だからである。
［後処理技術］
より効果的に経時変化を可視化するために、または検知された経時変化について定量的な情報を得るために、２枚の画像のサブトラクション処理後に後処理を行うことができる。この例は、ウィンドウ処理、エッジ強調またはぼかし、及び／または閾値処理等である。また、例えば米国特許第４，９０７，１５６号、第４，８３９，８０７号、第５，０７２，３８４号及び米国特許出願第０７／８４３，７１５号で述べられているような胸部画像の異常を検知するためのコンピュータ支援診断（ＣＡＤ）技術もサブトラクション画像の後処理として適用することができる。
［サブトラクション画像の出力］
サブトラクション画像は最終的なデータ出力として読影者が見ることができる。原画像をサブトラクション画像と比較するためにサブトラクション画像は原画像と並べて表示されるのが好ましい。サブトラクション画像はビデオ表示装置のようなソフトコピー、またはレーザフィルムプリンタやサーマルプリンタといったプリントアウトのようなハードコピーのいずれかの形で表示することができる。画像は、例えばウィンドウ処理のような適当なコントラスト強調処理を行った後表示するのが好ましい。フィルタリング、拡大といったその他のタイプの画像処理も画質を向上させるために表示に先だって実行することができる。
【００５０】
サブトラクション画像で経時変化を発見することが可能になる。サブトラクション画像で読影者はそれらの変化の起こっている可能性のある場所についての情報を得ることができる。実際に起こっている可能性のある経時変化を検知するために、既存のＣＡＤ技術で使用されている閾値処理、フィルタリング、ヒストグラム解析、その他の画像解析技術が種々の画像解析処理技術が用いられる。経時変化が起こっている可能性のある検知された場所は、サブトラクション画像、または原画像に重ねて表示される矢印のような記号を使って独立した図形または数値データとして示される。経時変化が起こっている可能性のある検知された部分についてのサイズ、形状、経時変化についてのその他の重要な特徴といった場所以外の特徴もまた出力される。起こっている可能性のある経時変化についての情報を、サブトラクション画像上で得る代わりに原画像の組で得ることも可能である。
【００５１】
自動化された画像位置合わせ技術を用いたサブトラクション技術では、相互相関値（平均値）と処理時間はどちらもテンプレートのサイズが大きくなると増加することが実験によって示された。２４×２４よりも小さいテンプレートは望ましい結果を示さなかった。これらのサイズではシフト値にかなり大きいエラーが含まれるからである。４８×４８や６４×６４といった大きなサイズのテンプレートでは一般的に結果として得られる画像で良い位置合わせが見られた。しかし、ある領域では微細な構造のミスマッチが生じた。何故ならば、大きなテンプレートにおける最良のマッチングは常に中心における最良のマッチングを示すことに対応する訳では無いからである。３２×３２のサイズのテンプレートで充分な結果が得られる。非線形濃度補正技術の効果も露光量を１．４１倍と０．８４倍にした原画像を用いて実験的に調べられた。その結果、次のようなことがわかった。
【００５２】
（ａ）露光量が０．５倍以下の画像についてのみ、濃度補正は胸郭の端の検出能を高める。
（ｂ）濃度補正はローカルマッチングに影響しない。
（ｃ）濃度補正は背景を均一にするので、経時変化が起こっている可能性のある場所を強調するのに効果的である。
【００５３】
自動位置合わせプログラムは前処理としての濃度補正に適用できるように変更された。２０組のサブトラクション画像（１８組にはinterbval 変化があり、２組には経時変化が無かった）が非線形濃度補正を含んだ自動化画像位置合わせ技術を用いて得られた。
【００５４】
その結果得られた画像は８×１０インチのフィルムにレーザフィルムプリンタでプリントされた。平均的な処理時間はＶＡＸステーション３５００で約３５分であった。得られた結果はほとんどの画像が非常に良好な位置合わせができたことを示していた。２枚のサンプルでは大きな位置合わせの誤りが生じた。そのうちのひとつは激しい前後（Ａ−Ｐ）の傾きが原因であり、もう一つは探索領域ＲＯＩの選択に誤りがあったためである。以前に調べた手動による位置合わせと較べて、サブトラクション画像は様々な経時変化をより強調することがわかった。この技術は陰性を調べるときにも有用であろう。変化がない画像のサブトラクション画像は、経時変化を示す画像に較べて、極めて均一な画像濃度を示す。ローカルマッチング処理で得られた相互相関値を定量的に解析するために、この解析を行うためのプログラムが作成された。これらの値の累積ヒストグラムが計算され、小さな相互相関値を持つテンプレートの場所が原画像上で特定される。このプログラムを使って先の２０組が解析され、その結果、低い相互相関しか示さないある領域の周りにテンプレートが位置していることが分かった。これらの領域は以下のようなところである。
【００５５】
（ａ）経時変化がある大きな病変。
（ｂ）経時変化がある心臓または横隔膜の境界。
（ｃ）経時変化があるかまたはそれがない大きな浸潤。
（ｄ）電極等身体に取り付けられた医療用装置。
（ｅ）角度が変わった鎖骨。
（ｆ）非常に濃度が低い領域。
【００５６】
（ｇ）激しいＡ−Ｐの傾きを持つ画像組の肺領域全体。
本発明は上述した実施例に限定されず、種々変形して実施可能である。上述の教示内容を照らして、本発明に対して様々な修飾及び変形が可能であることは自明なことである。
以上説明したように本発明によれば、時間的に連続したディジタル医用画像の間の経時変化を検知するために強調した画像を用い診断の精度を高めることができるコンピュータ化された方法及び装置が提供される。
【図面の簡単な説明】
【００５７】
【図１】本発明による画像位置合わせ及びサブトラクション処理の全体を示す図。
【図２】本発明による２枚の胸部画像のディジタル化からサブトラクション画像を出力するまでの基本的な流れを示す図。
【図３】ＲＯＩの選択処理を実行する方法を示す図。
【図４】マッピングされたシフト値に適用されるカーブフィッティングを実行する各ステップを示す図。
【図５】探索領域とテンプレートを含む第１、第２のディジタル化画像を示すＸ線写真。
【図６】テンプレートの大きさに対する相互相関値の関係を示す図。
【図７】テンプレートの大きさに対するローカルマッピングに要する処理時間の関係を示す図。
【図８】ローカルマッチング技術に使用される様々な形のテンプレートと探索領域ＲＯＩを示す図。
【図９】相互相関法を用いたローカルマッチング技術の一般的な原理を示す図。
【図１０】２枚の画像の対応する位置間のシフトマッピングの結果を示す図。
【図１１】シフトマッピング結果に対する重み付けフィッティングの一般的な原理を示す図。
【図１２】図１０に示されているシフトマップから２次元のｎ次多項式を用いて生成されたフィッティングされたシフト値ΔｘとΔｙの分布のシフトマップを示す図。
【図１３】本発明に従ったｘ，ｙ座標の非線形歪みの概念を示す図。
【図１４】図１３に示されているような非線形歪みを与えられた画像に対応する格子状パターンのコンピュータ出力結果を示す図。
【図１５】図１６とともに本発明の技術を用いるためのシステムの全体像を示す図。
【図１６】図１５とともに本発明の技術を用いるためのシステムの全体像を示す図。
【図１７】図１８とともに本発明のシステムにおけるカーブフィッティング演算装置とともに使用される重み係数を説明する図。
【図１８】図１７とともに本発明のシステムにおけるカーブフィッティング演算装置とともに使用される重み係数を説明する図。
【図１９】１組の時間的に連続した胸部画像の間で位置合わせの誤りを起こす様々な原因を示す図。
【図２０】１組の時間的に連続した胸部画像の間で位置合わせの誤りを起こすさらに別の様々な原因を示す図。
【図２１】患者の右上の肺葉に見られた空気空間浸潤を持つ胸部Ｘ線画像を示すＸ線写真。
【図２２】図２１に対応し３ヶ月前に撮影された空気空間浸潤が見られない胸部Ｘ線画像を示すＸ線写真。
【図２３】図２１の画像と図２２の画像から得られたサブトラクション画像を示すＸ線写真。
【図２４】患者の胸部Ｘ線画像を示すＸ線写真。
【図２５】図２４と同じ患者の肺の第２の画像を示すＸ線写真。
【図２６】図２４の画像を歪ませた画像を示すＸ線写真。
【図２７】図２６に示す歪まされた画像と図２５に示す第２の画像から得られたサブトラクション画像を示すＸ線写真。
【図２８】複数の結節を持つ胸部画像を示すＸ線写真であり、肺の中にある左側の不透明な領域として見えているように左上の肺葉にある結節の１つは空洞を発達させている。
【図２９】図２８の２年前に撮影されたＸ線写真であり、複数の結節は図２８に見られるよりも目立っていない。
【図３０】図２８の画像と図２９の画像から得られたサブトラクション画像を表すＸ線写真。
【図３１】心臓肥大の存在を示すＸ線写真。
【図３２】図３１に見られる心臓肥大がある患者の２年前に撮影された肺の胸部Ｘ線画像を示すＸ線写真。
【図３３】本発明の技術にしたがって図３１の画像と図３２の画像から得られたサブトラクション画像を示すＸ線写真。
【図３４】本発明の技術を用いずに画像全体をシフトさせて位置合わせをすることによって得られたサブトラクション画像に現れる様々なアーチファクトによる影響を示すＸ線写真であり、ここではその影響は骨と血管のミスマッチにより非常に大きなものとなっている。
【符号の説明】
【００５８】
１００…ディジタル画像入力装置、１１０，１２０…ディジタル画像メモリ、１３０…ＲＯＩ位置選択装置、１４０…ＲＯＩ画像メモリ、１６０…シフトマップ発生装置、１７０…シフトマップメモリ、１８０…カーブフィッティング演算装置、１９０…フィットシフトマップメモリ、２００…歪み画像メモリ、２１０…サブトラクション演算装置、２２０…サブトラクション画像メモリ、２３０…表示装置。
【図１】
2005176402.tif 000010
【図２】
2005176402.tif 000011
【図３】
2005176402.tif 000012
【図４】
2005176402.tif 000013
【図５】
2005176402.tif 000014
【図６】
2005176402.tif 000015
【図７】
2005176402.tif 000016
【図８】
2005176402.tif 000017
【図９】
2005176402.tif 000018
【図１０】
2005176402.tif 000019
【図１１】
2005176402.tif 000020
【図１２】
2005176402.tif 000021
【図１３】
2005176402.tif 000022
【図１４】
2005176402.tif 000023
【図１５】
2005176402.tif 000024
【図１６】
2005176402.tif 000025
【図１７】
2005176402.tif 000026
【図１８】
2005176402.tif 000027
【図１９】
2005176402.tif 000028
【図２０】
2005176402.tif 000029
【図２１】
2005176402.tif 000030
【図２２】
2005176402.tif 000031
【図２３】
2005176402.tif 000032
【図２４】
2005176402.tif 000033
【図２５】
2005176402.tif 000034
【図２６】
2005176402.tif 000035
【図２７】
2005176402.tif 000036
【図２８】
2005176402.tif 000037
【図２９】
2005176402.tif 000038
【図３０】
2005176402.tif 000039
【図３１】
2005176402.tif 000040
【図３２】
2005176402.tif 000041
【図３３】
2005176402.tif 000042
【図３４】
2005176402.tif 000043