JP6731919B2 - 光コヒーレンス断層撮影像におけるスペックル低減 - Google Patents

Description

本発明は、光コヒーレンス断層撮影(OCT: optical coherence tomography)像を処理してスペックルを低減した像を形成する方法とシステムとに関する。

光コヒーレンス断層撮影は、光信号の取得と処理の方法である。これは、マイクロメートルの分解能で、光散乱媒体(例えば、生体組織)内からの三次元像を撮影する。OCTを使用した眼像形成の原理は、光束(「A走査像」)が目に入ったとき、各光インターフェースから反射された光の時間遅延を測定することに基づく。構造物にわたる一連のA走査像は、前眼部または後眼部を通る平面の断面再現を可能にする。これは、B走査像と呼ばれる。

スペックルは、OCT像形成だけでなく、合成開口レーダー(SAR: synthetic aperture radar)、遠隔検出、超音波、ソナーなどを含む他の像形成法において問題となる。これらの分野では、スペックルノイズに対処するため、かなりの量の信号処理解析が実施されており、そのことが、エッジの鮮明さを保ちながらスペックルを減らすための強力なデジタルフィルタの開発につながっている。OCTにおけるスペックル低減は、像形成の設定を変更することにより、および、信号処理を使用して対処される。

既存のスペックル低減技術は、主に、1)Leeフィルタ[1]、Kuanフィルタ[2]、拡張LEEフィルタ[3]などのフィルタ、2)提案される、異方性拡散[4]、スペックル低減異方性拡散(SRAD: speckle reduction anisotropic diffusion)[5]、詳細保存異方性拡散[6]、および指向性スペックル低減異方性拡散[7]などの異方性拡散を使用したフレームに基づく。これらの方法は、最初、レーダー像のために提案されたが、OCT像に適用可能である。

Topconは、1,050nmの波長を使用した、後眼部像形成のための掃引光源OCTであるDRI OCT-1を開発した。これは、100,000 A走査像/秒の速い走査速度を示す。この1,050nmの波長を使用することで、DRI OCT-1は、850nmの範囲の波長を使用した従来受け入れられているOCTと比較してより深く突き抜け得る。脈絡膜または、さらには強膜などの眼組織を可視化することは、非常に短時間で実行され得る。深部像形成DRI OCT-1は、網膜だけでなく脈絡膜、さらには強膜の詳細な構造を表示する。DRI OCT-1は、散乱にあまり影響を与えずにより深い解剖学的構造を像形成する能力をもち、高感度で全体の断層写真を可視化し得る。

しかし、Topcon DRI OCT-1掃引光源OCTにおいて、像の単一のフレームまたはスライスは、通常、大きなスペックルノイズに起因して低品質である。したがって、従来、移動平均が計算される。例えば、256フレームの3Dモードでは、DRI OCT-1における3から5の隣接フレームから移動平均が計算される。移動平均は、スペックルノイズを低減させることができるものとして示されるが、少数の隣接フレームにわたる平均に起因して、スペックルの低減が制限される。しかし、平均するときにより多くのフレームが含まれると、特に、下方にある網膜構造が急速に変化する領域において、視認可能なアーチファクトとしてぼやけが生じる。

本発明は、複数のA走査像により構成される光コヒーレンス断層撮影(OCT)像を分析する新規かつ有用な方法とシステムとを提供することを目的とし、各A走査像は、既定の方向において像形成された、構造の一次元走査像である。

概略的には、本発明は、各A走査像を、そのA走査像の周辺において隣接したA走査像の組、すなわち像形成方向に直交するように広がる二次元領域においてそのA走査像の最も近い隣接部から再現することを提案する。

これを表現する他の方法は、所与のA走査像に対して、隣接した走査像の組が、集合的に像形成方向に延びた円柱である凸状外殻を画定し、A走査像は、凸状外殻内の線であり、凸状外殻の周辺部に囲まれることである。A走査像は、凸状外殻内の実質的な中心線であり得る。像形成方向に沿って見ると、隣接した走査像の組は、二次元の外形をもつ。凸状外殻は、例えば、長方形の外形または楕円形の外形をもち得る。

まず、隣接部のA走査像が像形成方向に位置合わせされ、次に、位置合わせされたA走査像の行列Xが形成され、行列補完は、低スペックルノイズ像を取得するために実行される。

行列低減は、第1の行列Xと第2の行列との間の差分を表すコスト関数を最小化する第2の行列を求めることにより実行され、ここで第2の行列は、複雑性制約に従う。

第1の例において、第2の行列は単一の行列Lであり、複雑性制約は、第2の行列Lが既定の値以下のランクをもつように制約されることである。

第2の例において、複雑性制約は、第2の行列を低ランク行列Lと疎行列Sとの和として表し、低ランク行列は、第1の所定の値未満のランクをもつように制約され、疎行列Sは、第2の所定の値未満のいくつかの非ゼロ元を含むように制約される。

いずれの例においても、低ランク行列Lは、スペックルノイズが低減された行列Xの1つのバージョンであり、XのA走査像に対応する低ランク行列Lの列は、A走査像の低スペックルノイズ再現である。

本方法の一実施形態である方法は、OCT像形成品質を改善するため、Topcon掃引光源OCT像形成システムなどのOCT像形成システムにより実施され得る、

本発明の特定の可能な一形態は、複数のA走査像により構成される光コヒーレンス断層撮影(OCT)像を改善するコンピュータ実装方法であって、各A走査像が、既定の方向において像形成された構造を表し、A走査像が、1つまたは複数のB走査像としてグループ化され、各A走査像に対して、
(a)A走査像の周辺における各隣接部を画定するステップと、
(b)各隣接部において既定の方向におけるA走査像を他のA走査像と位置合わせするステップと、
(c)位置合わせされたA走査像を使用して第1の行列Xを形成するステップであって、第1の行列の各列が位置合わせされたA走査像である、ステップと、
(d)第1の行列Xと第2の行列との間の差分を表すコスト関数を最小化する第2の行列を求めるステップであって、第2の行列が、複雑性制約に従うように制約される、ステップと、
を含む方法である。

改善されたバージョンのA走査像は、第2の行列の列から取得され得る。代替的に、第2の行列に対する複雑性制約が、第2の行列の複数の行列要素の各々に対する複雑性制約の形態である場合、改善されたバージョンのA走査像は行列要素の1つの列から取得され得る。いずれの例においても、フィルタ処理ステップなどのさらなる処理ステップが、改善されたバージョンのA走査像を取得するために実行され得る。

本発明は、本方法を実行するためのコンピュータシステムの観点からも表現され得る。コンピュータシステムのプロセッサにより実行されて、本方法をプロセッサが実行することをもたらし得るプログラム命令を含むコンピュータプログラム製品(例えば、非一時的な形態でプログラム命令を記憶する有形のデータ記憶装置)としても表現され得る。

本発明の一実施形態が、以下の図面を参照しながら例示のみを目的として以下で説明される。

本発明の一実施形態である方法を示すフロー図である。 三次元ボリュームのOCT像を示す図であり、OCT像は、N個のB走査像B₁、B₂、...、B_Nである。 (a)長方形の隣接部、(b)円形/楕円形の隣接部として、A走査像A_i,jの隣接部を画定する2つの方法を示す図である。 図1に示す方法のA走査像位置合わせステップを示す図である。 図1に示す方法の第1の例を使用した像再現の結果を示す図である。 図1に示す方法の第2の例を使用した像再現の結果を示す図である。 スペックル低減の3つの方法を比較したCNRの箱形図である。 図8(a)と図8(d)とは、生OCTデータを示す図であり、図8(b)と図8(e)とは、生データの移動平均を示す図であり、図8(c)と図8(f)とは、図1に示す方法により得られるデータを示す図である。

図1を参照して、本発明による方法が説明される。本方法の入力は、従来のOCT法により生成され得る種類のOCTボリュームである。OCTボリュームは、図2に示されるようにN個のB走査像(本分野において、「フレーム」または「スライス」とも呼ばれる)B_iにより作られ、i=1、2、...、Nである。各B走査像B_iは、M個のA走査像、A_i,jからなり、j=1、2、...、Mであり、各A走査像が、図1において垂直である方向に延びた一次元光束を使用して被写体を像形成することにより形成された一次元走査像である。言い換えると、OCTボリュームは、A走査像線A_i,jの組からなり、i=1、2、...、Nであり、j=1、2、...、Mであり、Nは、ボリューム内のB走査像の数であり、Mは、各B走査像内のA走査像の数である。

図1に示す方法の元にある基本概念は、その元となるきれいな像l_i,jとノイズg_i,jとの和として、すなわち次式のようにA_i,jを定式化することである。
A_i,j=l_i,j+g_i,j
目的は、A_i,jからl_i,jを再現することである。l_i,jとn_i,jとに関するいくつかの仮定がなければこの問題を解くことはできないが、本実施形態は、小さな隣接部内でl_i,jが非常に似ているという仮定により動機付けされる。この仮定に基づいて、隣接したA走査像を、類似モデルを使用してモデル化する。
A_i+p,j+q=l_i+p,j+q+g_i+p,j+q
式中、(p,q)∈Ωであり、Ωは「隣接部」(または、「隣接領域」)を表し、A走査像に隣接してA走査像を囲むA走査像の組である。したがって、A走査像の組は、像形成方向に見たとき二次元配列を形成する。

隣接領域を定義する多くの方法がある。第1の可能な定義は、図3(a)に示されるように、長方形領域{(p,q)|-n≦p≦n,-m≦q≦m}である。別の可能な定義は、図3(b)に示されるように楕円形領域

である。n=mである場合、隣接部は、それぞれ、正方形と円形とである。以下で報告される実験において、隣接部はn=1とm=2とを使用して、図3(a)に示す長方形の隣接部として定義された。

この概念を実装する1つの方法は、Xを、列としてΩからのすべてのA走査像を積み重ねることにより形成された行列と定義することであり、すなわち、X=[A_i-n,j-m,...,A_i+n,j+m]、L=[l_i-n,j-m,...,l_i+n,j+m]、G=[g_i-n,j-m,...,g_i+n,j+m]である。L内の各列は、非常に似ていると仮定され得るので、Lのランクは、低いと予想され、したがって、本方法は、所与のXからLを解くために低ランク行列補完を使用し得る。

残念ながら、上記の技術を生OCTボリュームに直接適用することは、隣接したA走査像における内容の変化に起因して、A走査像の何らかのぼやけをもたらし得る。加えて、避けられない眼球運動も何らかの位置ずれをもたらす。この理由により、図1は、ステップ3において行列の形成によるA走査像の再現が実行される前に、A走査像を位置合わせするステップ2を含む。

さらに、元のOCT A走査像が大きなスペックルノイズによる影響を受け、このことが、ステップ2において、これらのA走査像を位置合わせすることを困難な作業とする。したがって、前処理ステップ1は、ステップ2が実行される前にノイズを低減するために使用される。

多くのアルゴリズム[1]〜[7]が、ノイズ低減ステップ1に使用され得る。本発明者らは、SRADアルゴリズム[5]を使用する。SRADは、各走査像B_iに適用される。具体的には、強度像I₀(x,y)が与えられると、PDEに従って、次式のように出力像I(x,y;t)が生成される。

式中、c(q)は、

または、

および

により計算された拡散係数を表す。後述の実験結果において、繰り返しが50回実行された。

A走査像位置合わせステップ2において、(p,q)∈Ωとしたとき、隣接したA走査像A_i+p,j+qは、間の差分が最小化されるように、垂直な平行移動Δ_p,q(すなわち、上または下への移動)によりA_i,jに位置合わせされ、Ωは隣接したA走査像の位置を表す。発明者らの経験によると、A走査像のみから差分を計算することは、局所スペックルノイズにより支配されるという結果をもたらし得る。したがって、本発明者らは、A_i,jとA_i+p,j+qとの間の差分を、それらに隣接したA走査像も使用して計算する。

ステップ2において、pとqとの各々に対して、A_i+p,j+qとA_i,jとを中心として配置されたA走査像のグループ(「窓」)、すなわち、G_i+p,j+q=[A_i+p,j+q-w,A_i+p,j+q-w+1,...,A_i+p,j+q+w]およびG_i,j=[A_i,j-w,A_i,j-w+1,...,A_i,j+w]を定義する。式中、wは、窓寸法を規定し、位置合わせの計算をするために使用される。平行移動されたG_i+p,j+qとG_i,jとの間の差分が最小化されるように、G_i+p,j+qにおける垂直な平行移動Δ_p,qを得る。多くのブロックマッチングアルゴリズムが、位置合わせを計算するために使用され得る。実施形態における発明者らの数値解析において、その効率と実装の容易さとを理由としてダイヤモンド探索法[8]を使用して、垂直な平行移動を得る。

ここで、本発明者らは2つの方法を提示し、この2つの方法によりステップ3、すなわち2つの例のステップ3が実行され得る。両方の方法における第1のサブステップは、行列Xの列として、位置合わせされたA走査像を組み合わせることである。すなわち、X=[A_i-n,j-m,...,A_i+n,j+m]であり、式中、角括弧の内容は、位置合わせされたA走査像であると理解される。

方法1:まず、本実施形態は、低ランク部分とノイズ部分との和として、Xを分解する。すなわち、
X=L+G (1)
式中、Lは、rank(L)≦rの低ランク部分であり、Gはノイズである。以下で報告される実験において、本発明者らは、r=2に設定する。

上記の分解は、分解誤差を最小化することにより解かれる。

あらゆる行列Xの表記||C||_Fは、行列のフロベニウスノルム(すなわち、行列の列を積み重ねることにより形成されたベクトルのユークリッドノルム)を表す。

上記の最適化問題は、バイラテラルランダムプロジェクション[9]アルゴリズムにより解かれる。アルゴリズムは、次のように要約される。

Lを解くためのアルゴリズム1

ランダム行列A₁を生成する。
A₁∈R^Dxrであり、式中、DはA走査像における画素数である。

を計算する(本発明者らは、q=2を使用する)。

標準的な方法を使用して、QRの分解Y₁=Q₁R₁、Y₂=Q₂R₂を計算する。式中、Q₁とQ₂とは直交行列であり、R₁とR₂とは上三角行列である。

方法2:第2の方法において、本実施形態は、低ランク部分と疎部分とノイズ部分との和として、Xを分解する。すなわち、
X=L+S+G (3)
式中、Lは、rank(L)≦rの低ランク部分であり、Sは、card(S)≦kの疎部分であり、Gはノイズである。以下で報告される実験において、rは2に等しく、kは行列Sのすべての元の30%に設定された。

上記の分解は、分解誤差を最小化することにより解かれる。

上記の最適化問題は、収束するまで、次の2つの副問題を代わりに解くことにより解かれる。

アルゴリズムは、次のように要約される。

Lを解くためのアルゴリズム2

t:=0;L₀=X;S₀=0に初期化する。

While文開始

が満たされる間、以下を実行する。

t:=t+1

ランダム行列A₁∈R^Dxrを生成する。式中、DはA走査像における画素数である。

を計算する。

QRの分解Y₁=Q₁R₁、Y₂=Q₂R₂を計算する。

If文開始

である場合、

とし、第1のステップに飛ぶ。If文終了。

とする。式中、Ωは、|X-L_t|の最初のk個の最大要素の非ゼロの部分集合であり、

は行列の、Ωへの射影である。

while文終了

本明細書において「列」という用語は、行列の行も含む意味で使用されること、言い換えると、原理的には、Xは、Xの行として、位置合わせされたA走査像を積み重ねることにより形成され得る一方で、方法1と方法2とにおいて式が解かれる前にXの転置が適用されることに注意されたい。

第1の方法は、S=0とした第2の方法の特別な場合であることが理解され得る。主な相違点は、第2のモデルがノイズから、内容の変化に起因した差分を分離することである。

図1は、さらなる処理ステップを示さないが、任意選択的に、再現されたA走査像ボリュームに様々なフィルタが適用され得る。適切なフィルタは、拡張Lee、適応ウィーナーフィルタなどを含む。例えば、適応ウィーナーフィルタが使用され得る。加えて、ノイズ除去を補助するため、閾値が適用され得る。数学的に、本発明者らは、以下の工程を適用する。

式中、T_lowとT_highとは、2つの適応的に決定された閾値である。

結果
まず、方法1(すなわち、ステップ3の第1の例)に関する実施形態の使用の結果を示す。便宜上、本発明者らは、再現されたA走査像から形成されたB走査像について検証する。図5は、以下を含む像復元の例を示す。

(i)元データのB走査像(「生の像」と表示される)。

(ii)復元された像(「きれいな部分」と表示される)。各列は、方法1のそれぞれの実行により取得された。すなわち、生の像における所与の列に対して、対応する行列Xが形成される。方法1により、各Xから対応する行列LとGとが形成される。生の像の各列に対して、「きれいな部分」の像の対応する列は、対応する行列Xにおいて生の像の列がとる位置と同じ位置を対応する行列Lにおいてとる、対応する行列Lの列である。

(iii)「ノイズ部分」の像。生の像の各列に対して、「ノイズ部分」の像の対応する列は、対応する行列Xにおいて生の像の列がとる位置と同じ位置を対応する行列Gにおいてとる、対応する行列Gの列である。

図6は、ステップ2の第2の例(すなわち、方法2)を使用してステップ3が実行されたときの結果を示し、XをX=L+S+Gとモデル化する。

(i)元データのB走査像(「生の像」と表示される)。

(ii)復元された像(「きれいな部分」と表示される)。各列は、方法1のそれぞれの実行により取得された。すなわち、生の像における所与の列に対して、対応する行列Xが形成される。方法2により、各Xから対応する行列L、S、およびGが形成される。生の像の各列に対して、「きれいな部分」の像の対応する列は、対応する行列Xにおいて生の像の列がとる位置と同じ位置を対応する行列Lにおいてとる、対応する行列Lの列である。

(iii)「疎部分」の像。生の像の各列に対して、「疎部分」の像の対応する列は、対応する行列Xにおいて生の像の列がとる位置と同じ位置を対応する行列Sにおいてとる、対応する行列Sの列である。

(iv)「ノイズ部分」の像。生の像の各列に対して、「ノイズ部分」の像の対応する列は、対応する行列Xにおいて生の像の列がとる位置と同じ位置を対応する行列Gにおいてとる、対応する行列Gの列である。

両方の場合において、像Lは元の像Xよりはるかに鮮明である。視覚的に、ステップ3の2つの変形例による2つの結果は、視認可能な大きな違いはなく非常に似ている。

像の品質を評価するため、コントラスト-ノイズ比(CNR: contrast to noise ratio)の尺度[6]が使用される。これは、対象領域と背景領域との組の間のコントラストを測定する。

式中、μ_rとσ_rとは、すべての対象領域の平均と分散とであり、μ_bとσ_bとは、背景領域の平均と分散とである。様々な網膜層からの10個の対象領域と、各B走査像からの1つの背景領域とが、CNRを計算するために使用された。

試験のため、256個のB走査像である992×512の像を含む2つの典型的なOCTボリュームが使用される。Table 1(表1)は、(i)生ボリュームデータと、(ii)移動アルゴリズムを使用した従来の方法と、(iii)方法2を使用した実施形態との間の比較を示す。見てわかるとおり、本実施形態は、著しく大きなCNRをとる。図7は、箱形図を使用して、様々なフレームにおけるCNRの変形例を示す。

Table 1 生データと、従来の移動平均スペックル低減アルゴリズムと、本発明の一実施形態との間CNRの比較

図8は、2つの例示的な生データを使用した、提案される方法と基準となる方法との間の視覚的な違いを示す(図8(a)および図8(d))。図8(b)と図8(e)とに、移動平均法により生成された像が示され、図8(c)と図8(d)とに、方法2を使用した本発明の一実施形態により生成された像が示される。

生フレームとそれらの移動平均と比較すると、本実施形態は、構造検査において、はるかに良好なボリューム品質を提供する。網膜層および篩板などの構造に関するより多くの詳細事項が、より快適に検査され得る。

さらに、バイラテラルランダムプロジェクションを使用した、低ランク復元によりステップ3を実行する場合は、ロバストPCA[8]などの従来のアルゴリズムより約20倍速い。

計算工程は、プロセッサとプロセッサによる実装のためのプログラム命令を記憶するデータ記憶装置とを含み、実質的に自動的に(すなわち、任意選択的に方法の開始を除いて、人が関与せず)プロセッサが上記の方法を実施することをもたらすコンピュータシステムを使用して実装され得る。コンピュータシステムは、方法の結果を表示するための画面を含む表示装置、および/または、方法の結果を印刷された形態で生成するための印刷装置、および/または、方法の結果を電子的に伝達するためのインターフェースをさらに備え得る。コンピュータシステムは、インターフェースから生の像を受信する汎用コンピュータであり得る。コンピュータシステムは、代替的に、まず像を撮影した後、実施形態を使用して像を処理するOCT像形成装置の一部であり得る。

結果として得られる像は、患者が病気を患っていることの診断を行うため、および/または、被術者に実行される手術工程などの医療処置を選択するための道具の1つとして医療関係者により使用され得る。

(参考文献)
以下の参考文献の開示内容は、参照により本明細書に組み込まれる。
[1] J. S. Lee、「Speckle analysis and smoothing of synthetic aperture radar images」、Computer Graphics Image Processing、17巻、24〜32頁、1981
[2] D. T. Kuanm、A. A. Sawchuk、T. C. Strand、P. Chavel、「 Adaptive noise smoothing filter for images with signal dependent noise」、IEEE Trans. Pat. Anal. Mach. Intell.、vol. 7、165〜177頁、1985
[3] A. Lopes、R. Touzi、E. Nesby、「Adaptive speckle filter and scene heterogeneity」、IEEE Trans. Geosci. Remote Sens.、28巻、992〜1000頁、1990
[4] P. Perona、J. Malik、「Scale-space and edge detection using anisotropic diffusion」、IEEE Trans. Patt. Ana. Mach. Intell.、12巻、629〜 639頁、1990
[5] Y. Yu、S. T. Acton、「Speckle reducing anisotropic diffusion」、IEEE Trans. Image Processing、11巻、1260〜1270頁、2002
[6] S. Aja-Fernandez、C. Alberola-Lopez、「On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering」、IEEE Trans. Image Processing」、15巻、9号、2006
[7] K. Krissian、C. F. Westin、R. Kikinis、K. G. Vosburgh、「Oriented speckle reducing anisotropic diffusion」、IEEE Trans. Image Processing、16(5)巻、1412〜1424頁、2007
[8] S. Zhu、K. Ma、「A new diamond search algorithm for fast block-matching motion estimation」、IEEE Trans. Image Processing、9巻、287〜290頁、2000
[9] T. Zhou、D. Tao、「GoDec: Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case」、International Conference on Machine Learning、33〜40頁、2011

Claims (11)

1. 複数のA走査像により構成される光コヒーレンス断層撮影(OCT)像を改善するためのコンピュータで実行される方法であって、各々のA走査像が、既定の方向において像形成された構造を表し、前記A走査像は、1つまたは複数のB走査像としてグループ化され、前記方法は、各々のA走査像に対して、
   (a)前記各々のA走査像の周辺におけるA走査像の組からなる各隣接部を画定するステップであって、前記A走査像の組は、前記既定の方向に見たときの二次元配列を形成する、ステップと、
   (b)前記各隣接部における前記A走査像の組を前記各隣接部の前記各々のA走査像に位置合わせするステップと、
   (c)前記位置合わせされたA走査像を使用して第1の行列Xを形成するステップであって、前記第1の行列Xの各列が前記位置合わせされたA走査像である、ステップと、
   (d)前記第1の行列Xと第2の行列との間の差分を表すコスト関数を最小化する前記第2の行列を求めるステップであって、前記第2の行列が複雑性制約に従うように制約される、ステップとを含む、
   方法。
2. 前記複雑性制約は、前記第2の行列が既定の値以下のランクをもつように制約されることである、
   請求項1に記載の方法。
3. 前記複雑性制約は、低ランク行列Lと疎行列Sとの和として前記第2の行列を表すものであり、前記複雑性制約は、前記低ランク行列Lが第1の所定の値未満のランクをもつように制約されることであり、前記疎行列Sが、第2の所定の値未満のいくつかの非ゼロ元を含むように制約される、
   請求項1に記載の方法。
4. 前記複数のA走査像におけるノイズ低減の事前ステップを含む、
   請求項1から請求項3のいずれか一項に記載の方法。
5. 前記ノイズ低減が、スペックル低減異方性拡散アルゴリズムである、
   請求項4に記載の方法。
6. ステップ(b)が、
   前記A走査像を含む1つまたは複数のA走査像からなる第1のA走査像グループを定義するステップと、
   各隣接したA走査像に対して、
   (i)第2の各A走査像グループを定義するステップであって、前記第2の各A走査像グループが、前記各隣接部の前記各隣接したA走査像を含む1つまたは複数のA走査像からなる、ステップと、
   (ii)前記各第2のA走査像グループの前記1つまたは複数のA走査像が平行移動値分だけ前記既定の方向に平行移動される場合、前記第1のA走査像グループと前記第2のA走査像グループとの間の差分を最小化する各平行移動値を求めるステップと、
   (iii)前記隣接したA走査像を前記各平行移動値分だけ平行移動させるステップと、
   により実行される、
   請求項1から請求項5のいずれか一項に記載の方法。
7. ステップ(b)が、ブロックマッチングアルゴリズムにより実行される、
   請求項6に記載の方法。
8. 前記ブロックマッチングアルゴリズムが、ダイヤモンド探索アルゴリズムである、
   請求項7に記載の方法。
9. 前記二次元配列が、長方形領域または楕円形領域である、
   請求項7に記載の方法。
10. プロセッサと、前記プロセッサにより実行されたとき前記プロセッサに請求項1から請求項9のいずれか一項に記載の方法を実行させるように動作するプログラム命令を記憶するように構成されたメモリ装置とを備える、
    コンピュータシステム。
11. プロセッサにより実行されたとき前記プロセッサに請求項1から請求項9のいずれか一項に記載の方法を実行させるためのプログラム命令を有する、
    コンピュータプログラム。