JP4899170B2 - 画像処理システム、画像処理方法、プログラムおよび記録媒体 - Google Patents

Description

本発明は、画像に現れている多数の粒子像を検出し得る画像処理技術に関し、特に、電子顕微鏡像のような低Ｓ／Ｎ比の画像に現れている粒子像を検出し得る画像処理技術に関する。

タンパク質の複合体、核酸、脂質あるいは多糖類などの生体高分子の立体構造を解析する従来技術として、生体高分子を結晶化させた試料を作成し、この試料をＸ線結晶構造解析法により解析し評価するという技術が使用されている。ところが、かかる技術は試料の結晶化を必要とするため、膜タンパク質などの結晶化が困難な試料の解析にこの技術を適用することができなかった。なお、核磁気共鳴（ＮＭＲ：Nuclear Magnetic Resonance）法は、試料の結晶化を必要としないものの、大きな分子量を有する生体高分子を解析することが一般に難しい。

そこで、試料の結晶化を必要としない方法として、極低温透過型電子顕微鏡（Cryo-ＴＥＭ：Cryo-Transmission Electron Microscope）を用いてタンパク質の構造解析を行う単粒子構造解析法が提案されている。この単粒子構造解析法は、極低温透過型電子顕微鏡中で極低温に冷却された試料から得た高解像度のＴＥＭ画像を使用する。試料には、多数の粒子状の生体高分子（以下「粒子」と呼ぶ。）が含まれており、これら生体高分子はランダムな方向を向いている。単粒子構造解析法は、このＴＥＭ画像から個々の粒子を表す局所画像（以下「粒子画像」と呼ぶ。）を切り出し、これら粒子画像に基づいて粒子の３次元構造を推定する方法である。この方法によれば、同一方向を向いた粒子の粒子画像を加算平均することにより、画像のノイズ成分を低減させることができる。

単粒子構造解析法に関する先行技術文献としては、たとえば、非特許文献１（佐藤主税他：「単粒子解析法による結晶を用いない蛋白質構造解析：電圧依存性Ｎａ^＋チャンネルの構造を例として」，電子顕微鏡，37(1), 40-44, 2002）が挙げられる。

以下、図１を参照しつつ、単粒子解析法の一例を説明する。図１は、非特許文献１に記載されている単粒子解析法の概要を説明するための図である。生体高分子の電子顕微鏡像は、電子線照射による試料の損傷を防ぐために分解能が制限され、ノイズ成分を多く含む画像となる。タンパク質のように電子線損傷を受けやすい生体高分子を撮像した場合、電子顕微鏡像のＳ／Ｎ比が低いので、電子顕微鏡の持つ本来の分解能を十分に発揮させることが難しい。それ故、１枚の電子顕微鏡像から十分な情報が得られないだけでなく、生体高分子の立体構造の再構成も容易ではないという問題がある。単粒子解析法は、かかる問題を克服する方法として注目されている。

単粒子解析法では、先ず、試料の電子顕微鏡像（ELECTRON MICROGRAPH）に現れている多数の粒子を特定する。粒子の特定は、肉眼により行われる。次いで、電子顕微鏡像からこれら粒子それぞれの局所画像（粒子画像）を切り出す。次に、これら局所画像について位置合わせと回転合わせが実行される。その後、粒子画像は、粒子が向いている方向別に複数のグループに分類される。同じグループ内の粒子画像を平均化することで高Ｓ／Ｎ比の２次元平均化像（投影像）を得、当該２次元平均化像を用いて粒子の３次元立体構造を推定する。凍結法による無染色電子顕微鏡像は分子内密度を反映した投影像を与えると考えられる。よって、高Ｓ／Ｎ比の２次元平均化像から１次元に投影した像（シノグラム）を生成し、これらシノグラムを相互に比較することで各２次元平均化像に現れている粒子の３次元上の方位角（オイラー角）や共有回転軸を推定することができる（図１左下）。そして、２次元平均化像、方位角および共有回転軸に基づいて粒子の立体構造が再構成される。立体構造の再構成方法は、たとえば、非特許文献２（Joachim Frank, "Three-Dimensional Electron Microscopy of Macromolecular Assemblies", Oxford University Press, 2006, pp. 193-276）に開示されている。

さらに、こうして得られた立体構造から再投影像を生成し、この再投影像を利用して２次元平均化像の位置合わせをした後に、さらに平均化することで新たな２次元平均化像を生成する。そして、新たな２次元平均化像に基づいて粒子の立体構造が再構成される。このような再投影像を利用した位置合わせと立体構造の再構成というプロセスを繰り返し実行することで高い分解能で立体構造を推定することが可能となる。

単粒子解析法に関する先行技術文献としては、非特許文献１や非特許文献２のほかに特許文献１（特開２００８−１３４２１２号公報）が挙げられる。
佐藤主税，上野豊，小椋俊彦，藤吉好則：「単粒子解析法による結晶を用いない蛋白質構造解析：電圧依存性Ｎａ＋チャンネルの構造を例として」、電子顕微鏡、第３７巻、第１号、第４０〜４４ページ（２００２年）。 Joachim Frank, "Three-Dimensional Electron Microscopy of Macromolecular Assemblies", Oxford University Press, pp. 193-276 (2006). van Heel M., "Multivariate statistical classification of noisy images (randomly oriented biological macromolecules)", Ultramicroscopy, vol.13, pp.165-183 (1984). van Heel M., "Classification of very large electron microscopical image data sets", Optik, vol. 82, pp.114-126 (1989). 特開２００８−１３４２１２号公報

前述の如く、単粒子解析法では、電子顕微鏡像から粒子画像を切り出すために、電子顕微鏡像に現れている多数の粒子を特定する必要があるが、各粒子の特定は肉眼により行われていた。しかしながら、電子顕微鏡像のＳ／Ｎ比が低く、そのコントラスト（電子顕微鏡像の最も明るい部分と最も暗い部分との輝度差の大きさ）が低い場合には、電子顕微鏡像の中から、ノイズに埋もれた数千〜数万もの粒子を肉眼で識別することは非常に難しく、また、その識別には多大な時間と熟練を要するという問題がある。

肉眼で粒子を識別する代わりに、パターンマッチング法により、参照画像（テンプレート）を用いて電子顕微鏡像の中から粒子を探し出す画像処理を利用すれば、短時間で粒子を特定することができる。しかしながら、単なるパターンマッチング法では、粒子の検出精度が非常に低いという問題がある。

上記に鑑みて本発明の目的は、生体高分子のような粒子像を低Ｓ／Ｎ比の画像の中から高精度かつ短時間で検出し得る画像処理システム、画像処理方法、プログラムおよび記録媒体を提供することである。

本発明によれば、複数の対象粒子像を含む全体画像から複数の部分画像を抽出する部分画像抽出部と、前記対象粒子像を模した参照粒子像の画像データが格納されている画像記憶部と、前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するピーク検出部と、前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するベクトル算出部と、前記全体画像の領域に対して前記位置設定ベクトルにより指定された位置に前記ピークの位置を設定して２次元ピーク分布を生成するピーク分布生成部と、を備えた画像処理システムが提供される。この画像処理システムでは、前記部分画像抽出部は、前記複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように前記全体画像から前記複数の部分画像を抽出する。

また、本発明によれば、（ａ）複数の対象粒子像を含む全体画像から複数の部分画像を抽出するステップと、（ｂ）前記対象粒子像を模した参照粒子像の画像データを画像記憶部から読み出し、前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するステップと、（ｃ）前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するステップと、（ｄ）前記全体画像の領域に対して前記位置設定ベクトルを用いて指定された位置に前記ピークの位置を設定して２次元ピーク分布を生成するステップと、を備えた画像処理方法が提供される。この画像処理方法では、複数の部分画像は、当該複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように全体画像から抽出される。

また、本発明によれば、コンピュータ読み取り可能な記録媒体から読み出されてコンピュータに画像処理を実行させるプログラムが提供される。この画像処理は、複数の対象粒子像を含む全体画像から複数の部分画像を抽出する部分画像抽出処理と、前記対象粒子像を模した参照粒子像の画像データを画像記憶部から読み出し、前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するピーク検出処理と、前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するベクトル算出処理と、前記全体画像の領域に対して前記位置設定ベクトルを用いて指定された位置に前記ピークの位置を設定して２次元ピーク分布を生成するピーク分布生成処理と、を含み、前記部分画像抽出処理では、前記複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように前記全体画像から前記複数の部分画像が抽出される。

さらに、本発明によれば、前記プログラムを記録した記録媒体が提供される。

本発明による画像処理システム、画像処理方法およびプログラムは、いずれも、複数の対象粒子像を含む全体画像から抽出された複数の部分画像を利用する。これら部分画像は、所定方向に隣り合う部分画像同士が重複するように抽出されるので、隣り合う部分画像が同一の対象粒子像を含み得るものである。本発明による画像処理システム、画像処理方法およびプログラムでは、各部分画像と参照粒子像との間の類似度分布のピークが検出されるとともに、各部分画像内の基準位置からピークの位置に至る位置設定ベクトルが算出される。そして、全体画像の領域に対して位置設定ベクトルにより指定された位置にピークの位置を設定することで２次元ピーク分布が生成される。位置設定ベクトルは、対象粒子像の存在確率が比較的高い領域を指すベクトルであるので、２次元ピーク分布では、対象粒子像の位置にピークが集中しやすい。それ故、低Ｓ／Ｎ比の画像からでも、２次元ピーク分布に基づいて高精度かつ短時間で対象粒子像を検出することが可能となる。

以下、本発明の実施の形態について、図面を用いて説明する。なお、すべての図面において、同様な構成要素には同様の符号を付し、適宜説明を省略する。

（画像処理システムの構成）
図２は、本発明に係る一実施形態の画像処理システムである粒子構造解析装置１の概略構成を示す機能ブロック図である。この粒子構造解析装置１は、粒子像検出部２と立体構造解析部３とを有する。図３は、粒子像検出部２の概略構成を示す機能ブロック図である。

図２に示されるように、粒子像検出部２には、電子顕微鏡像である全体画像が入力されている。この電子顕微鏡像は、タンパク質の複合体、核酸、脂質あるいは多糖類などの生体高分子を含む試料を透過型電子顕微鏡で撮像したＴＥＭ（Transmission Electron Microscope）画像である。粒子像検出部２は、ＴＥＭ画像に現れている観測粒子像（以下「対象粒子像」と呼ぶ。）の位置を高精度かつ短時間で特定し、その結果を立体構造解析部３に与える。立体構造解析部３は、粒子像検出部２で特定された対象粒子像の立体構造を推定する機能を有する。立体構造解析部３の機能は、たとえば、非特許文献１や特許文献１に開示されている技術を利用したものであればよく、特に制限されるものではない。

なお、粒子像検出部２と立体構造解析部３は、単一の装置内に組み込まれていてもよいし、あるいは、小規模または大規模なコンピュータネットワーク（たとえば、ＬＡＮ（Local Area Network）やＷＡＮ（Wide Area Network））を介して相互接続されていてもよい。

図３に示されるように、粒子像検出部２は、画像入力部１０、部分画像抽出部１１、ピーク検出部１２、ベクトル算出部１３、ピーク分布生成部１４、粒子像探索部１５、データ記憶部１６、粒子画像抽出部１７、画像分類部１８、平均化処理部１９、更新処理部２０および画像メモリ２１を有する。コントローラ２２は、機能ブロック１０〜２１の動作を制御して全体の処理を司る機能を有している。データ記憶部１６は、画像記憶部１６１、参照画像記憶部１６２、ピーク分布画像記憶部１６３および派生画像記憶部１６４を含む。

図３の粒子像検出部２の機能ブロック１０〜１５，１７〜２０，２２の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック１０〜１５，１７〜２０，２２の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。

データ記憶部１６は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。各種記憶部１６１〜１６４は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、粒子像検出部２の動作時に割り当てられる適当な記憶領域上に構成されてもよい。

図３を参照すると、画像入力部１０には、多数の対象粒子像を含む電子顕微鏡像が全体画像として入力される。画像入力部１０は、全体画像を表す転送信号を受信する受信器であってもよいし、あるいは、記録媒体から全体画像のデータを読み出す装置であってもよい。画像入力部１０は、入力された全体画像のデータを画像記憶部１６１に格納する。

部分画像抽出部１１は、コントローラ２２からの指示に応じて、複数の対象粒子像を含む全体画像を画像記憶部１６１から読み出し、当該読み出された全体画像から部分画像Ｐ_１〜Ｐ_Ｍ（Ｍは正整数）を抽出する。後述するように、部分画像Ｐ_１〜Ｐ_Ｍのうち所定方向に隣り合う部分画像Ｐ_ｉ，Ｐ_ｉ＋１が互いに重複するように部分画像Ｐ_１〜Ｐ_Ｍが抽出される。これにより、隣り合う部分画像Ｐ_ｉ，Ｐ_ｉ＋１が同一の対象粒子像の全部または一部を含むことが可能となる。

参照画像記憶部１６２は、全体画像に含まれ得る対象粒子像を模した参照粒子像を含む参照画像Ｒ_１〜Ｒ_Ｎ（Ｎは正整数）を格納している。本実施形態では、これら参照画像Ｒ_１〜Ｒ_Ｎは、それぞれ、粒子が向いている方向が異なる同一種の粒子像を表す画像である。たとえば、参照画像Ｒ_１は対象粒子の正面を模した像、参照画像Ｒ_２は対象粒子の底面を模した像、参照画像Ｒ_３は粒子の側面を模した像とすることができる。

ピーク検出部１２は、各部分画像Ｐ_ｊについて、当該各部分画像Ｐ_ｊと参照粒子像との間の類似度の２次元分布（以下「類似度分布」と呼ぶ。）を算出し、類似度分布に現れるピークを検出する。

より具体的には、先ず、ピーク検出部１２は、参照画像記憶部１６２から参照画像Ｒ_１〜Ｒ_Ｎを読み出す。上述した通り、参照画像Ｒ_１〜Ｒ_Ｎは、それぞれ、対象粒子像を模した参照粒子像を含む。次に、ピーク検出部１２は、各部分画像Ｐ_ｊに対して各参照画像Ｒ_ｋを移動させつつ、各部分画像Ｐ_ｊと各参照画像Ｒ_ｋとの間の類似度を順次算出することにより類似度分布ＤＳ_ｊ，ｋ（ｘ，ｙ，θ）を得る。ここで、ｘ、ｙは、各部分画像Ｐ_ｊ対する各参照画像Ｒ_ｋのｘ方向とｙ方向への移動量（シフト量）である。また、θは、各部分画像Ｐ_ｊ対する各参照画像Ｒ_ｋの点（ｘ、ｙ）周りの回転移動量である。そして、ピーク検出部１２は、この類似度分布ＤＳ_ｊ，ｋ（ｘ，ｙ，θ）に現れるピークを検出する。

たとえば、ピーク検出部１２は、各部分画像Ｐ_ｊの中から着目画素領域を順次選択するとともに、当該着目画素領域ごとに、当該着目画素領域を含む画像領域と参照画像Ｒ_ｋとの間の類似度を算出することにより類似度分布ＤＳ_ｊ，ｋ（ｘ，ｙ，θ）を得る。この場合、ｘ、ｙは、各部分画像Ｐ_ｊに対する各参照画像Ｒ_ｋの基準位置から着目画素領域へのシフト量であり、θは、着目画素領域の周りの回転移動量である。着目画素領域は、１画素単位で順次選択されてもよいし、あるいは、４画素×４画素といった複数画素単位で順次選択されてもよい。類似度としては、公知の相互相関係数が好適であるが、これに限定されるものではない。たとえば、着目画素領域を含む画像領域と参照粒子像との間の残差二乗和を変数とする算術関数（たとえば、１／Ｘ^２：Ｘは残差二乗和）を類似度とすることもできる。

ピーク検出部１２は、各部分画像Ｐ_ｊについて、類似度分布ＤＳ_ｊ，ｋ（ｘ，ｙ）から閾値以上の類似度を有するピークを検出し、この検出結果をベクトル算出部１３とピーク分布生成部１４に与える。ここで、各部分画像Ｐ_ｊについて最大の類似度を持つ単一のピークのみを検出してもよいし、あるいは、各部分画像Ｐ_ｊについて複数個のピークを検出してもよい。以下、各部分画像Ｐ_ｊについて最大の類似度を持つ単一のピークのみを検出し、当該ピークを用いて対象粒子像を検出する方法を「ＭＲＰ（Multi-Reference Pickup）法」と呼び、各部分画像Ｐ_ｊについて少なくとも１つのピークを検出し、当該ピークを用いて対象粒子像を検出する方法を「ＭＲＭＰ（Multi-Reference Multiple Pickup）法」と呼ぶこととする。

ベクトル算出部１３は、各部分画像Ｐ_ｊについて当該各部分画像Ｐ_ｊの局所領域内の所定の基準位置（たとえば、部分画像Ｐ_ｊの中心）からピークＰｋの位置に至る位置設定ベクトルを算出する。各ピークＰｋとこれに対応する位置設定ベクトルとの組がピーク分布生成部１４に与えられる。

ピーク分布生成部１４は、各部分画像Ｐ_ｊについて位置設定ベクトルにより指定された位置にピークＰｋの位置を設定して、全体画像領域に対応する２次元ピーク分布を生成する。ここで、ピークＰｋの位置は、全体画像領域における各部分画像Ｐ_ｊの位置（大域的な位置）を基点として位置設定ベクトルで指定される位置に設定される。この２次元ピーク分布は、ピーク分布画像記憶部１６３に格納される。より具体的に説明すると、ピーク分布生成部１４は、全体画像領域における各部分画像Ｐ_ｊの所定の基準位置（ｘ_ｊ，ｙ_ｊ）を、当該各部分画像Ｐ_ｊについて算出された位置設定ベクトルｖ＝（δｘ_ｊ，δｙ_ｊ）だけシフトさせ、シフトされた基準位置Ｐ（ｘ_ｊ＋δｘ_ｊ，ｙ_ｊ＋δｙ_ｊ）を位置設定ベクトルｖ＝（δｘ_ｊ，δｙ_ｊ）により指定された位置とする。

そして、粒子像探索部１５は、ピーク分布画像記憶部１６３から２次元ピーク分布を読み出し、この２次元ピーク分布に現れるピークの密度（以下「ピーク密度」と呼ぶ。）が閾値以上となる領域（対象粒子像の存在確率が高い領域）を探索して各対象粒子像の位置を検出する。ここで、ピーク密度は、単位面積当たりのピークの個数を意味する。ピーク密度が閾値以上となる領域は、２次元ピーク分布が局在化した領域である。低Ｓ／Ｎ比の全体画像に基づいて生成された２次元ピーク分布には、ノイズに埋もれた対象粒子像と同じ領域またはその周辺領域に多数のピークが現れることが期待されるので、ピーク密度が高い領域では対象粒子像の存在確率が高い。

さらに、粒子像探索部１５は、２次元ピーク分布に現れるピークの位置の類似度が閾値以上となる領域を探索し、その探索結果に基づいて対象粒子像の位置を検出する機能を有する。ここで、粒子像探索部１５は、２次元ピーク分布に現れる複数のピークの位置が重複するときは、これら複数のピークが有する類似度の加算結果（累積類似度）に基づいて対象粒子像の位置を検出する。Ｓ／Ｎ比が比較的高い画像領域では、同一位置に複数のピークが重複して現れる場合があり、その位置に対象粒子像が存在する確率は非常に高い。粒子像探索部１５は、そのような対象粒子像の位置を確実に検出することができる。

図３に示されるように、粒子像探索部１５は、ピーク密度算出部１５０と位置検出部１５１とを含む。ピーク密度算出部１５０は、２次元ピーク分布に対して探索窓を走査させつつ、この探索窓内のピークの個数を計数し、その計数値に基づいてピーク密度を順次算出する機能を有している。位置検出部１５１は、（i）ピーク密度が第１の閾値以上となったとき、もしくは、（ii）探索窓内のピーク位置の類似度または累積類似度が第２の閾値以上となったとき、当該探索窓の位置に対象粒子像が位置すると判定する。あるいは、位置検出部１５１は、（i），（ii）の双方の検出条件が満たされたときに限り、当該探索窓の位置に対象粒子像が位置すると判定することもできる。ユーザは、全体画像の画質に応じて第１の閾値および第２の閾値をそれぞれ任意の値に設定することにより検出精度を調整することが可能である。粒子像探索部１５は、このようにして検出された対象粒子像の位置を示すデータを粒子画像抽出部１７に与える。

なお、探索窓の形状としては、少なくとも１つの画素を含む円形状や矩形状が挙げられるが、特に限定されるものではない。

粒子画像抽出部１７は、粒子像探索部１５により検出された位置にある対象粒子像を含む粒子画像を全体画像から抽出し（切り出し）、これら粒子画像のデータを画像メモリ２１を介して図２の立体構造解析部３に出力する機能を有する。図２の立体構造解析部３は、粒子像検出部２から入力された粒子画像に基づいて対象粒子像の立体構造を推定する。

粒子像検出部２は、全体画像から抽出された粒子画像に基づいて、参照画像記憶部１６２に記憶されている参照粒子像の品質を向上させる機能ブロック１８〜２０を有している。

画像分類部１８は、粒子画像から派生画像を生成し分類する機能を有している。コントローラ２２の制御により、画像分類部１８は、粒子画像抽出部１７から入力された粒子画像を、参照画像記憶部１６２内の参照画像Ｒ_１〜Ｒ_Ｎと比較することにより、派生画像を１以上のグループに分類する。

具体的には、画像分類部１８は、粒子画像抽出部１７から入力された各粒子画像Ｓ（ｉ）に基づいて、参照画像記憶部１６２内の参照画像Ｒ_１〜Ｒ_Ｎとの間でピーク類似度を有する派生画像Ｓ（ｉ，０）〜Ｓ（ｉ，Ｐ_ｉ）（Ｐ_ｉは正整数）を生成する。ここで、ｉは、粒子画像を一意に特定する番号である。派生画像は、数学的には、参照画像Ｒ_１〜Ｒ_Ｎに対して各粒子画像Ｓ（ｉ）の相対位置を、Ｘ軸並進方向（水平画素方向）、Ｙ軸並進方向（垂直画素方向）および角度方向にそれぞれ移動させることで得られる画像である。類似度として相互相関係数を使用する場合、各参照画像Ｒ_ｋの画素値を要素とするベクトルと、各派生画像の画素値を要素とするベクトルとの内積すなわち相関係数に比例する量を類似度として算出することができる。生成された派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）（Ｔは２以上の正整数）は、派生画像記憶部１６４に一時的に格納される。

次に、画像分類部１８は、派生画像記憶部１６４から派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）を読み出し、これら派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）を１以上のグループに分類する。具体的には、画像分類部１８は、派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）の各々に対してベクトル化処理を行う。派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）が、それぞれ、水平画素数ｎと垂直画素数ｍを有する２次元画像であるとすれば、ベクトル化処理では、各２次元画像のｎ×ｍ個の画素値を要素とするｎ×ｍ次元の位置ベクトルが生成される。

派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）は、それぞれ対応する位置ベクトルにより、ｎ×ｍ次元のベクトル空間内に位置付けられる。画像分類部１８は、派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）を、ｎ×ｍ次元のベクトル空間（画素空間）内での近接度に応じてグループ分けすることができる。たとえば、このベクトル空間内で、互いの距離が閾値以下となるように局所的な空間内に偏って分布する所定数以上の派生画像を１つのグループとして抽出することが可能である。

あるいは、画素空間と共役なイメージ空間内で、局所的な空間内に偏って分布する派生画像を１つのグループとして抽出することもできる。イメージ空間の次元が画素空間の次元よりも低い場合には、イメージ空間で派生画像を分類し、画素空間の次元がイメージ空間の次元よりも低い場合には、画素空間で派生画像を分類することが好ましい。これにより、派生画像の分類処理を効率的に行うことができる。このような分類法は、対応分析（corresponding analysis）を用いた階層分類法（hierarchical ascendant classification）として知られている。階層分類法は、たとえば、非特許文献３（van Heel M., "Multivariate statistical classification of noisy images (randomly oriented biological macromolecules)", Ultramicroscopy, vol.13, pp.165-183 (1984)）や非特許文献４（van Heel M., "Classification of very large electron microscopical image data sets", Optik, vol. 82, pp.114-126 (1989)）に開示されている。

平均化処理部１９は、複数の派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）の画素値をグループごとに平均化することによりグループごとに２次元平均化像を生成する。この平均化は、ランダムノイズを除去し得る平均化（たとえば、算術平均化）であればよく、特に限定されるものではない。

なお、平均化処理部１９は、派生画像の画素値に当該派生画像の類似度に応じた係数を重み付けてもよい。派生画像の類似度が高いほど、当該派生画像の画素値に重み付けるべき係数を大きな値に設定し、一方、派生画像の類似度が低いほど、当該派生画像の画素値に重み付けるべき係数を小さな値に設定することが好ましい。これにより、２次元平均化像の品質向上が可能となる。類似度の値と重み付け係数の値との対応関係を示すルックアップテーブル（図示せず）を予め用意しておけば、平均化処理部１９は、このルックアップテーブルを参照して、類似度に応じた最適な重み付け係数を設定することが可能である。

画像分類部１８は、下限Ｓ_ｍａｘ−Δと上限Ｓ_ｍａｘとの幅Δを各粒子画像または全体画像のＳ／Ｎ比に応じて変動させてもよい。すなわち、各粒子画像または全体画像のＳ／Ｎ比が高いほどΔを小さくし、粒子画像のＳ／Ｎ比が低いほどΔを大きくすることが好ましい。これにより、平均化処理部１９で生成される２次元平均化像の品質の向上が可能となる。粒子画像のＳ／Ｎ比と幅Δの値との対応関係を示すルックアップテーブル（図示せず）を予め用意しておけば、画像分類部１８は、このルックアップテーブルを参照して、各粒子画像または全体画像のＳ／Ｎ比に応じた最適な幅Δを設定することが可能である。

更新処理部２０は、参照粒子像を含む参照画像Ｒ_１〜Ｒ_Ｎを２次元平均化像のデータで置き換えることで画像記憶部１６１の記憶内容を更新する機能を有する。すなわち、更新処理部２０は、平均化処理部１９で生成された２次元平均化像群の品質が参照画像Ｒ_１〜Ｒ_Ｎのそれよりも良好か否かを判定する。この判定は、両方の画像のＳ／Ｎ比や輪郭抽出結果を対比することで実行される。そして、更新処理部２０は、２次元平均化像群の品質が参照画像Ｒ_１〜Ｒ_Ｎのそれよりも良好であると判定したときに、参照画像群を２次元平均化像群で置き換える。これにより、ピーク検出部１２は、更新された参照画像Ｒ_１〜Ｒ_Ｐ（Ｐは正整数）を用いて、より高精度に類似度分布のピークを検出することが可能となる。

（粒子構造解析処理）
以上の構成を有する粒子構造解析装置１の動作を以下に説明する。図４は、粒子構造解析装置１による粒子構造解析処理の手順の一例を概略的に示すフローチャートである。図４の処理では、参照画像記憶部１６２の記憶内容を更新する処理は実行されない。

図４を参照すると、ステップＳ１では、部分画像抽出部１１は、複数の対象粒子像を含む全体画像を画像記憶部１６１から読み出し、続けて、全体画像から部分画像Ｐ_１〜Ｐ_Ｍを抽出する（ステップＳ２）。これら部分画像Ｐ_１〜Ｐ_Ｍのうち所定方向（たとえば、Ｘ軸方向）に隣り合う部分画像Ｐ_ｉ，Ｐ_ｉ＋１は、同一の対象粒子像を含み得るように互いに重複する。

図５は、全体画像３０と部分画像Ｐ_１〜Ｐ_Ｍの一部Ｐａ〜Ｐｆとを模式的に示す図である。図５に示されるように、部分画像Ｐａ〜ＰｃはＸ軸方向に互いに重複するように全体画像３０から抽出される。部分画像Ｐｄ〜Ｐｆについても同様である。よって、部分画像Ｐａ〜Ｐｃは、検出されるべき同一の対象粒子像３１を含み、部分画像Ｐｄ〜Ｐｆは、検出されるべき同一の対象粒子像３２を含む。部分画像Ｐａ，Ｐｂ，Ｐｃの位置は、それぞれ、全体画像３０の有効領域の原点Ｏを基点とするベクトルＶａ，Ｖｂ，Ｖｃで特定される。

次に、図４を参照すると、ステップＳ３で、ピーク検出部１２は、処理対象となる部分画像Ｐ_ｉを表す番号ｉを初期値「１」に設定して（ステップＳ３）、ピーク検出処理を開始する。ピーク検出部１２は、参照画像記憶部１６２から参照画像Ｒ_１〜Ｒ_Ｎを読み出し、これら参照画像Ｒ_１〜Ｒ_Ｎと部分画像Ｐ_ｉとの間の類似度分布ＤＳ_ｉ，ｋ（ｘ，ｙ，θ）を算出し（ステップＳ４）、次いで、類似度分布ＤＳ_ｉ，ｋ（ｘ，ｙ，θ）から、閾値以上の類似度を有する単数または複数のピークＰｔ（ｉ，１）〜Ｐｔ（ｉ，Ｑ）（Ｑは正整数）を検出する（ステップＳ５）。

ただし、粒子像検出部２がＭＲＰ法による処理モードで動作する場合、ピーク検出部１２は、コントローラ２２からの指示に従って類似度分布ＤＳ_ｉ，ｋ（ｘ，ｙ，θ）から最大の類似度を持つ単一のピークＰｔ（ｉ，１）のみを検出する。一方、粒子像検出部２がＭＲＭＰ法による処理モードで動作する場合、ピーク検出部１２は、コントローラ２２からの指示に従って類似度分布ＤＳ_ｉ，ｋ（ｘ，ｙ，θ）から複数のピークＰｔ（ｉ，１）〜Ｐｔ（ｉ，Ｑ）を検出することができる。

次のステップＳ６で、ベクトル算出部１３は、部分画像Ｐ_ｉの局所領域内の所定の基準位置（部分画像Ｐ_ｉの中心）からピークＰｔ（ｉ，１）〜Ｐｔ（ｉ，Ｑ）の位置に至る位置設定ベクトルＶ_ｉ（１）〜Ｖ_ｉ（Ｑ）を算出する。

次に、ベクトル算出部１３は、全ての部分画像Ｐ_１〜Ｐ_Ｍについてベクトル算出処理が終了したか否かを判定する（ステップＳ７）。全ての部分画像Ｐ_１〜Ｐ_Ｍについてベクトル算出処理が終了していない場合（ステップＳ７のＮＯ）、番号ｉがインクリメントされ（ステップＳ８）、ステップＳ４〜Ｓ７が繰り返し実行される。

一方、全ての部分画像Ｐ_１〜Ｐ_Ｍについてベクトル算出処理が終了している場合（ステップＳ７のＹＥＳ）、ピーク分布生成部１４は、全体画像の領域に対して、位置設定ベクトルＶ_ｉ（ｊ）（ｉは１〜Ｍの整数、ｊは１〜Ｑの整数）で指定される位置にピークＰｔ（ｉ，ｊ）の位置を設定することにより２次元ピーク分布を生成する（ステップＳ９）。その後、ピーク分布生成部１４は、この２次元ピーク分布をピーク分布画像記憶部１６３に記録する（ステップＳ１０）。

図６（Ａ），（Ｂ）は、それぞれ、異なる方向から視たときの同一の参照粒子像を含む参照画像Ｒ_１，Ｒ_２を例示する図である。図７（Ａ）〜（Ｃ）は、図６（Ｂ）の参照画像Ｒ_２と図５の部分画像Ｐａ，Ｐｂ，Ｐｃとに基づいて算出された位置設定ベクトルＶＳａ，ＶＳｂ，ＶＳｃを示す図である。図７（Ａ）〜（Ｃ）では、いずれも、対象粒子像３１の略中心位置にピークが検出されている。このため、位置設定ベクトルＶＳａ，ＶＳｂ，ＶＳｃは、それぞれ、部分画像Ｐａ，Ｐｂ，Ｐｃの中心から対象粒子像３１の略中心位置に至るベクトルである。

この場合、ピーク分布生成部１４は、全体画像の領域内で、部分画像Ｐａ，Ｐｂ，Ｐｃをそれぞれ位置設定ベクトルＶＳａ，ＶＳｂ，ＶＳｃだけシフトさせる。この結果、図７（Ａ）〜（Ｃ）に示されるように、部分画像Ｐａ，Ｐｂ，Ｐｃはそれぞれシフトして部分画像ＳＰａ，ＳＰｂ，ＳＰｃとなる。そして、ピーク分布生成部１４は、全体画像の領域内で、シフトされた部分画像ＳＰａ，ＳＰｂ，ＳＰｃの基準位置（たとえば、中心）にピークの位置を設定する（ステップＳ９）。

ステップＳ１０で２次元ピーク分布が記録された後は、粒子像探索部１５は、この２次元ピーク分布をピーク分布画像記憶部１６３から読み出し、この２次元ピーク分布のピーク密度およびピーク値（上記の類似度または累積類似度）を用いて対象粒子像を探索することにより各対象粒子像の位置を検出する（ステップＳ１１）。

図８は、２次元ピーク分布５０の一例を模式的に示す図である。図８に示されるように、２次元ピーク分布５０には、複数のピークＰｔ，Ｐｔ，…が現れている。粒子像探索部１５は、図８に示されるように、２次元ピーク分布５０に対して探索窓５１を走査させつつ、この探索窓内のピーク密度を順次算出する。同時に、粒子像探索部１５は、２次元ピーク分布５０に現れる複数のピークＰｔの位置が重複するときは、これら複数のピークが有する類似度を加算して累積類似度を得る。粒子像探索部１５は、（i）ピーク密度が第１の閾値以上となったとき、あるいは、（ii）探索窓内のピーク位置の類似度または累積類似度が第２の閾値以上となったときに、探索窓５１の位置に対象粒子像が存在すると判定し、その判定結果を粒子画像抽出部１７に与える。

その後、粒子画像抽出部１７は、画像記憶部１６１から全体画像を読み出し、この全体画像から、粒子像探索部１５による探索結果に基づいて対象粒子像が存在する画像領域を粒子画像として抽出する（ステップＳ１２）。粒子画像群は、画像メモリ２１を介して図２の立体構造解析部３に転送される。

立体構造解析部３は、粒子像検出部２から転送された粒子画像群に基づいて対象粒子像の立体構造を推定する（ステップＳ１３）。以上で粒子構造解析処理は終了する。

上記粒子構造解析処理が奏する効果は以下の通りである。

上述の通り、各部分画像Ｐ_ｉと参照粒子像との間の類似度分布のピークが検出されるとともに、各部分画像Ｐ_ｉ内の基準位置からピークの位置に至る位置設定ベクトルが算出される。そして、位置設定ベクトルにより指定された位置にピークの位置を設定することで２次元ピーク分布が生成される。位置設定ベクトルは、対象粒子像の存在確率が比較的高い領域を指すベクトルであるので、２次元ピーク分布では、対象粒子像の位置にピークが集中しやすい。それ故、低Ｓ／Ｎ比の全体画像からでも、２次元ピーク分布に基づいて高精度かつ短時間で対象粒子像を検出することが可能となる。

また、粒子像探索部１５は、２次元ピーク分布のピーク密度や累積類似度に基づいて対象粒子像を探索するので、高い精度で対象粒子像を検出することが可能である。さらに、粒子像探索部１５は、２次元ピーク分布に対して探索窓を走査させつつ、この探索窓内のピークの個数をピーク密度として順次算出する。それ故、少ない計算量で高速にピーク密度や累積類似度を順次算出して対象粒子像の位置を自動的に検出することができる。

特に、粒子像検出部２をＭＲＭＰ法による処理モードで動作させた場合は、低Ｓ／Ｎ比の全体画像中でノイズに埋もれた対象粒子像の存在する領域に複数のピークが現れやすい。ＭＲＭＰ法は、上述の通り、各部分画像について複数のピークを検出し得、複数のピークを用いて対象粒子像を検出する方法であるので、非常に低いＳ／Ｎ比の全体画像から対象粒子像を高い精度で検出することができる。

したがって、立体構造解析部３は、対象粒子の立体構造を高精度に推定することが可能である。

なお、対象粒子の検出精度を向上させるために、部分画像抽出部１１は、全体画像または部分画像の空間周波数成分のうち特定帯域の空間周波数成分のみを抽出するバンドパスフィルタを有していてもよい。このバンドパスフィルタを使用すれば、特定帯域の空間周波数成分から対象粒子像を検出することができる。たとえば、ハイパスフィルタを用いて、高帯域成分のみに存在するノイズを全体画像から除去したり、ローパスフィルタを用いて、個々の対象粒子像全体あるいは電子顕微鏡像全体で不均一な輝度分布を均一な輝度分布にしたりすることで、対象粒子像の検出精度を向上させることが可能である。

あるいは、部分画像抽出部１１は、複数の異なる帯域の空間周波数成分を出力するフィルタバンクを有していてもよい。離散ウェーブレット変換を実行するフィルタを使用すると、１つの全体画像から複数の異なる帯域の空間周波数成分を得ることができる。複数の異なる帯域の空間周波数成分に基づいてそれぞれ２次元ピーク分布を生成すれば、これら２次元ピーク分布の中から、対象粒子像の検出に最適な２次元ピーク分布を選択することが可能である。

対象粒子の検出精度を向上させるために、ピーク分布生成部１４は、２次元ピーク分布の空間周波数成分のうち特定帯域の空間周波数成分のみを抽出するバンドパスフィルタ、あるいは、複数の異なる帯域の空間周波数成分を出力し得るフィルタを有していてもよい。このようなフィルタを使用することで、全体画像の品質に応じて対象粒子像の検出に最適な帯域の空間周波数成分を得ることができる。

（繰り返し処理）
図９は、粒子構造解析装置１による粒子構造解析処理の別の手順を示すフローチャートである。図９のフローチャートのステップＳ１〜Ｓ１３は、図４のフローチャートのステップＳ１〜Ｓ１３と同一であるので、その詳細な説明を省略する。図９のフローチャートでは、参照画像記憶部１６２の記憶内容を更新するために、ステップＳ２０，Ｓ２１が新たに追加されている。

図９を参照すると、ステップＳ１２では、上述の通り、粒子画像抽出部１７が、画像記憶部１６１から全体画像を読み出し、この全体画像から、粒子像探索部１５による探索結果に基づいて、対象粒子像が存在する画像領域を粒子画像として抽出する。粒子画像は、画像分類部１８に出力される。

次のステップＳ２０では、コントローラ２２は、繰り返し処理を終了するか否かを判定する。繰り返し処理が終了しない場合（ステップＳ２０のＮＯ）、ステップＳ２１のフィードバック処理が実行される。図１０は、フィードバック処理の手順を概略的に示すフローチャートである。

図１０を参照すると、ステップＳ２２０では、画像分類部１８は、粒子画像抽出部１７から入力された各粒子画像Ｓ（ｉ）から複数の派生画像Ｓ（ｉ，０）〜Ｓ（ｉ，Ｐ_ｉ）を生成する。図１１（Ａ），（Ｂ）および図１２は、派生画像Ｓ（ｉ，０）〜Ｓ（ｉ，Ｐ_ｉ）の生成方法を説明するための図である。

図１１（Ａ）に示されるように、各粒子画像Ｓ（ｉ）に基づいて、Ｐ個のピーク類似度にそれぞれ対応するＰ個の派生画像Ｓ（ｉ，０）〜Ｓ（ｉ，Ｐ_ｉ＝Ｐ）が生成される。ｐ番目の派生画像Ｓ（ｉ，ｐ）は、数学的には、パラメータｘｐ，ｙｐ，θ＝ｔｈｐで定められる量だけ、それぞれ、Ｘ軸並進方向、Ｙ軸並進方向および角度θの方向に粒子画像Ｓ（ｉ）を移動させることで生成される画像である。図１１（Ｂ）に示されるように、ｘｐは、Ｘ軸並進方向の移動量、ｙｐは、Ｙ軸並進方向の移動量、ｔｈｐは、原点Ｏの周りの回転移動量、をそれぞれ表している。ここで、各粒子画像Ｓ（ｉ）について派生画像の個数は複数であるが、これに限定されるものではなく、派生画像の個数は各粒子画像Ｓ（ｉ）について１個であってもよい。

図１２は、粒子画像Ｓ（ｉ）の移動量（Ｘ軸並進移動量、Ｙ軸並進移動量および回転移動量）とピーク類似度との関係の一例を説明するための概略図である。図１２のグラフには、移動量に応じた複数のピークが現れている。画像分類部１８は、これらピークのうち、最大類似度Ｓ_ｍａｘを上限とし、Ｓ_ｍａｘ−Δを下限とする範囲内の類似度を有するピークＳＰ_ｍａｘ，ＳＰのみを選択する。画像分類部１８は、このような範囲内のピークに対応する派生画像Ｓ（ｉ，０）〜Ｓ（ｉ，Ｐ_ｉ）を生成する（ステップＳ２２０）。

全ての粒子画像Ｓ（１）〜Ｓ（Ｔ）について派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）が生成された後、ステップＳ２２１では、画像分類部１８は、派生画像記憶部１６４から派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）を読み出す。さらに画像分類部１８は、これら派生画像Ｓ（１，０）〜Ｓ（１，Ｐ_１），…，Ｓ（Ｔ，０）〜Ｓ（Ｔ，Ｐ_Ｔ）をベクトル化し、複数のグループに分類する。

次いで、平均化処理部１９は、分類された複数の派生画像の画素値をグループごとに平均化することによりグループごとに２次元平均化像を生成する（ステップＳ２２２）。そして、更新処理部２０は、参照粒子像を含む参照画像Ｒ_１〜Ｒ_Ｎのデータを２次元平均化像群のデータで置き換えることで画像記憶部１６１の記憶内容を更新する（ステップＳ２２３）。ここで、更新処理部２０は、２次元平均化像群の品質が、参照画像Ｒ_１〜Ｒ_Ｎのそれよりも良好である場合に限り、参照画像Ｒ_１〜Ｒ_Ｎを２次元平均化像群に置きかえてもよい。その後、コントローラ２２は、処理を図９のフローチャートに戻し、ステップＳ３〜Ｓ２０の手順を繰り返し実行させる。

繰り返し処理の実行回数が所定回数に到達した場合や、対象粒子像の検出精度が一定精度に到達した場合に、コントローラ２２は、繰り返し処理を終了すると判定して（ステップＳ２０のＹＥＳ）、画像分類部１８、平均化処理部１９および更新処理部２０の動作を停止させる。そして、立体構造解析部３は、粒子像検出部２から入力された粒子画像群に基づいて対象粒子像の立体構造を推定する（ステップＳ１３）。以上で粒子構造解析処理は終了する。

上記繰り返し処理が奏する効果は以下の通りである。

上記繰り返し処理では、参照画像Ｒ_１〜Ｒ_Ｎは、Ｓ／Ｎ比が向上した２次元平均化像に置き換えられる（ステップＳ２２２，Ｓ２２３）。それ故、図９のステップＳ３〜Ｓ１２，Ｓ２０〜Ｓ２１を繰り返し実行することにより、対象粒子像の検出精度を高めることが可能となる。

特に、フィードバック処理（図１０）では、粒子画像から複数の派生画像が生成され、これら派生画像が複数のグループに分類される（ステップＳ２２０〜Ｓ２２１）。それ故、２次元平均化像の品質を向上させることができ、対象粒子像の検出精度をさらに高めることが可能となる。

（テストによる評価）
次に、粒子構造解析装置１の動作テストの実行結果について説明する。粒子構造解析装置１の機能をコンピュータプログラムで実現し、このコンピュータプログラムを実行することによりテストが実行された。

図１３は、テストに使用された参照粒子像（モデル粒子像）を含む参照画像１００〜１０３を示す図である。図１３に示されるように、参照画像１００の参照粒子像を０度、９０度、１８０度および２７０度回転させることにより、それぞれ、４個の参照画像１００，１０１，１０２，１０３が生成された。これら４個の参照画像１００〜１０３を無作為に選択し、５１２０画素×５１２０画素の２次元画像内のランダムな位置に分布させた。さらに、この２次元画像にノイズを重畳させることにより１００個または２００個の対象粒子像を含む全体画像３０が生成された。

図１４（Ａ）は、ＳＮ比（＝Ｓ／Ｎ比）０．１の全体画像３０中の対象粒子像を、図１４（Ｂ）は、ＳＮ比０．０１の全体画像３０中の対象粒子像を、図１４（Ｂ）は、ＳＮ比０．００１の全体画像３０中の対象粒子像を、それぞれ示す図である。

ＭＲＰ法およびＭＲＭＰ法を用いた第１の動作モード（図４）により、このような全体画像３０から対象粒子像を検出した。図１５（Ａ）は、全体画像３０を示す図であり、図１５（Ｂ）は、この全体画像３０の一部領域Ｗｄを拡大表示した図である。図１５（Ｂ）に示されるように、白線の枠で囲まれた領域が、ＭＲＭＰ法を用いて検出され（拾い上げられ）た対象粒子像を示している。

図１６は、１００個の対象粒子像を拾い上げたときの拾い上げ効率のＳ／Ｎ比依存性を示すグラフである。このグラフの横軸は、全体画像３０のＳＮ比（＝Ｓ／Ｎ比）を示し、グラフの縦軸は、拾い上げ効率を示している。ここで、「拾い上げ効率」とは、全体画像３０中のモデル粒子像の実際の位置と、検出された対象粒子像の位置との一致率（％）である。たとえば、検出された１００個の対象粒子像の位置のうち５０個の位置がモデル粒子像の実際の位置と一致した場合、拾い上げ効率は５０％となる。図１６のグラフ中、記号「×」は、ＭＲＰ法による平均拾い上げ効率を、記号「○」は、ＭＲＭＰ法による平均拾い上げ効率を、それぞれ示している。実線は、ＭＲＰ法による平均拾い上げ効率を示す点を結んだものである。このグラフでは、判別しにくいが、ＳＮ比が、１、０．１、０．０１のそれぞれの拾い上げ効率の値は、ＭＲＰ法とＭＲＭＰ法とでほぼ同じであった。

図１７は、２００個の対象粒子像を拾い上げたときの拾い上げ効率のＳ／Ｎ比依存性を示すグラフである。このグラフの横軸は、全体画像３０のＳＮ比（＝Ｓ／Ｎ比）を示し、グラフの縦軸は、拾い上げ効率を示している。たとえば、検出された２００個の対象粒子像の位置のうち１００個の位置がモデル粒子像の実際の位置と一致した場合、拾い上げ効率は５０％となる。図１７のグラフ中、実線は、ＭＲＰ法による平均拾い上げ効率を、○の記号は、ＳＮ比０．００１の全体画像３０に対するＭＲＭＰ法による平均拾い上げ効率を、それぞれ示している。

図１６および図１７のグラフによれば、ともに、１０^−２オーダーの低Ｓ／Ｎ比の全体画像３０から１００％の効率で対象粒子像が拾い上げられていることが分かる。特に、ＭＲＭＰ法を使用して対象粒子像を拾い上げた場合は、ＭＲＰ法を使用して対象粒子像を拾い上げた場合と比べて、拾い上げ効率が向上していることが分かる。

図１８は、Ｓ／Ｎ比０．００１の全体画像３０から１００個の対象粒子像を拾い上げたときの拾い上げ効率のパラメータ依存性を示すグラフである。このグラフの横軸は、図８の探索窓５１の半径Ｒとピーク密度に対する閾値ρとの組み合わせ（パラメータ）を示している。ピーク密度が閾値ρ以上となったときに、探索窓５１の位置に対象粒子像が位置すると判定された。一方、図１９は、Ｓ／Ｎ比０．００１の全体画像３０から２００個の対象粒子像を拾い上げたときの拾い上げ効率のパラメータ依存性を示すグラフである。このグラフの横軸も、図８の探索窓５１の半径Ｒとピーク密度に対する閾値ρとの組み合わせ（パラメータ）を示している。図１８と図１９のグラフには、ともに、ＭＲＰ法による拾い上げ効率とＭＲＭＰ法による拾い上げ効率とが示されている。

図１８および図１９のグラフによれば、ＭＲＰ法とＭＲＭＰ法の双方で、パラメータを調整すれば、略１００％の拾い上げ効率を達成できることが分かる。また、ＭＲＭＰ法を使用した場合は、ＭＲＰ法を使用した場合と比べて、拾い上げ効率のパラメータの設定値に対する依存性が低いことが分かる。言い換えれば、ＭＲＭＰ法は、パラメータに対して優れた耐ノイズ性を有し、パラメータ変動に対して高いロバストネス（robustness）を有する。

図２０は、Ｓ／Ｎ比０．００１の全体画像３０から１００個の対象粒子像を拾い上げたとき、拾い上げられた対象粒子像の累積頻度分布を示すグラフである。図２１は、Ｓ／Ｎ比０．００１の全体画像３０から２００個の対象粒子像を拾い上げたとき、拾い上げられた対象粒子像の累積頻度分布を示すグラフである。図２０および図２１は、異なるパラメータ（半径Ｒと閾値ρの組み合わせ）が使用された場合の合計１６個のグラフからなり、各グラフには、使用されたパラメータの値が示されている。また、各グラフの縦軸は、拾い上げられた対象粒子像の累積頻度を示しており、各グラフの横軸は、全体画像３０中のモデル粒子像の実際の位置と、拾い上げられた対象粒子像の位置との間の距離（単位：ピクセル）を示している。各グラフの横軸は、１００個または２００個のモデル粒子像のうちｉ番目のモデル粒子像の位置Ｐ（ｉ）と、拾い上げられた１００個または２００個の対象粒子像のうち位置Ｐ（ｉ）に最近接の対象粒子像の位置Ｑ（ｉ）との間のズレを表している。各グラフ中、破線はＭＲＰ法を用いた場合の結果を、実線はＭＲＭＰ法を用いた場合の結果を、それぞれ示している。

図２０と図２１のグラフについては、累積頻度が０から１００へ収束するときの曲線の立ち上がり部分の傾きが大きいほど、拾い上げられた対象粒子像の位置ズレは小さいので、対象粒子像の検出精度は高い。図２０および図２１のグラフによれば、ＭＲＰ法とＭＲＭＰ法の双方で、パラメータを調整すれば、高い検出効率を達成できることが分かる。また、ＭＲＭＰ法を使用した場合は、ＭＲＰ法を使用した場合と比べて、パラメータの設定値に関係せずに、概ね、累積頻度の収束性が良好である。これらグラフから、ＭＲＭＰ法は、パラメータに対して優れた耐ノイズ性を有し、パラメータ変動に対して高いロバストネスを有することが分かる。

図２２および図２３は、ＴＲＰＭ２と呼ばれるタンパク質の観測粒子像を含む電子顕微鏡を全体画像３０としたときに得られたグラフである。図２２は、拾い上げ粒子数（拾い上げられた観測粒子像の個数）と相関値との関係を示すグラフであり、図２３は、相関値と拾い上げられた観測粒子像の頻度分布を示している。各グラフには、使用されたパラメータ（半径Ｒと閾値ρの組み合わせ）の値が示されている。また、各グラフ中、破線はＭＲＰ法を用いた場合の結果を、実線はＭＲＭＰ法を用いた場合の結果を、点線は、目視で観測粒子像を特定した場合の結果を、それぞれ示している。

図２２のグラフの相関値は、拾い上げられ（検出され）た観測粒子像とＴＲＰＭ２タンパク質の３次元構造の射影画像との間の相互相関係数である。図２２のグラフは、相関値が高いほど、観測粒子像の検出精度が高いことを意味している。図２２のグラフによれば、ＭＲＰ法とＭＲＭＰ法をそれぞれ使用した場合、パラメータを調整すれば、目視で観測粒子像を検出した場合よりも高い検出精度を達成できていることが分かる。たとえば、（Ｒ＝１、ρ＝１．００）、（Ｒ＝２、ρ＝０．５０）、（Ｒ＝２、ρ＝１．００）、（Ｒ＝３、ρ＝０．６７）、（Ｒ＝５、ρ＝０．４０）というパラメータを使用した場合は、目視で観測粒子像を検出した場合よりも、高い検出精度が得られている。

また、図２２のグラフによれば、ＭＲＭＰ法を使用した場合は、ＭＲＰ法を使用した場合と比べて、観測粒子像の検出精度のパラメータ依存性が低いことが分かる。よって、ＭＲＭＰ法は、パラメータに対して優れた耐ノイズ性を有し、パラメータ変動に対して高いロバストネスを有することが分かる。

図２３のグラフの相関値は、拾い上げられ（検出され）た観測粒子像とＴＲＰＭ２タンパク質の３次元構造の射影画像との間の相互相関係数である。頻度分布の偏りが相関値の高い方にシフトするほど、観測粒子像の検出精度が高い。図２３のグラフによれば、概ね、ＭＲＰ法とＭＲＭＰ法をそれぞれ使用した場合、目視で観測粒子像を検出した場合よりも高い検出精度を達成できていることが分かる。

以上、図面を参照して本発明の実施形態について述べたが、これらは本発明の例示であり、上記以外の様々な構成を採用することもできる。たとえば、上記実施形態では、粒子構造解析装置１に入力される全体画像は電子顕微鏡像であり、参照画像記憶部１６２に記憶されている参照画像Ｒ_１〜Ｒ_Ｎは、生体高分子像を模した粒子像であるが、これに限定されるものではない。全体画像が、生体高分子像以外の対象粒子像を含む低Ｓ／Ｎ比画像であってもよいし、参照画像Ｒ_１〜Ｒ_Ｎが、生体高分子像以外の対象粒子像を模した画像であってもよい。上記実施形態の粒子構造解析装置１は、たとえば、地上を撮した衛星画像から、建築物などの構造物の像を対象粒子像として検出することができる。

また、粒子像検出部２の処理を並列化することができる。たとえば、複数のプロセッサを使用して、ピーク検出部１２における類似度算出処理を、参照画像Ｒ_ｋ（ｋは１〜Ｎのうちの任意整数）別にあるいは部分画像別に並列化することができる。

単粒子解析法の概要を説明するための図である。 本発明に係る一実施形態の画像処理システムである粒子構造解析装置の概略構成を示す機能ブロック図である。 粒子像検出部の概略構成を示す機能ブロック図である。 粒子構造解析処理の手順を概略的に示すフローチャートである。 全体画像と部分画像の一部とを模式的に示す図である。 （Ａ），（Ｂ）は、それぞれ、異なる方向から視たときの同一の参照粒子像を含む参照画像を例示する図である。 （Ａ）〜（Ｃ）は、図６（Ｂ）の参照画像と図５の部分画像とに基づいて算出された位置設定ベクトルを示す図である。 ２次元ピーク分布を模式的に示す図である。 繰り返し処理を伴う粒子構造解析処理の手順を示すフローチャートである。 フィードバック処理の手順を概略的に示すフローチャートである。 派生画像の生成方法を説明するための図である。 移動量とピーク類似度との関係の一例を説明するための概略図である。 テスト用参照画像を示す図である。 （Ａ），（Ｂ），（Ｃ）は、それぞれ、０．１、０．０１、０．００１のＳＮ比（Ｓ／Ｎ比）を有するモデル粒子像（対象粒子像）を示す図である。 （Ａ）は、全体画像を示す図であり、（Ｂ）は、全体画像の一部領域Ｗｄを拡大した図である。 対象粒子像の拾い上げ効率のＳ／Ｎ比依存性を示すグラフである。 対象粒子像の拾い上げ効率のＳ／Ｎ比依存性を示すグラフである。 対象粒子像の拾い上げ効率のパラメータ依存性を示すグラフである。 対象粒子像の拾い上げ効率のパラメータ依存性を示すグラフである。 拾い上げられた対象粒子像の累積頻度分布を示すグラフである。 拾い上げられた対象粒子像の累積頻度分布を示すグラフである。 ＴＲＰＭ２タンパク質の観測粒子像を含む電子顕微鏡像から得られたグラフである。 ＴＲＰＭ２タンパク質の観測粒子像を含む電子顕微鏡像から得られたグラフである。

符号の説明

１ 粒子構造解析装置（画像処理システム）
２ 粒子像検出部
３ 立体構造解析部
１０ 画像入力部
１１ 部分画像抽出部
１２ ピーク検出部
１３ ベクトル算出部
１４ ピーク分布生成部
１５ 粒子像探索部
１５０ ピーク密度算出部
１５１ 位置検出部
１６ データ記憶部
１６１ 画像記憶部
１６２ 参照画像記憶部
１６３ ピーク分布画像記憶部
１６４ 派生画像記憶部
１７ 粒子画像抽出部
１８ 画像分類部
１９ 平均化処理部
２０ 更新処理部
２１ 画像メモリ
２２ コントローラ
３０ 全体画像
３１，３２ 対象粒子像

Claims (29)

1. 複数の対象粒子像を含む全体画像から複数の部分画像を抽出する部分画像抽出部と、
   前記対象粒子像を模した参照粒子像の画像データが格納されている画像記憶部と、
   前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するピーク検出部と、
   前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するベクトル算出部と、
   前記全体画像の領域に対して前記位置設定ベクトルにより指定された位置に当該ピークの位置を設定して２次元ピーク分布を生成するピーク分布生成部と、
   を備え、
   前記部分画像抽出部は、前記複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように前記全体画像から前記複数の部分画像を抽出する、画像処理システム。
2. 請求項１に記載の画像処理システムであって、前記ピーク分布生成部は、前記各部分画像内の当該所定の基準位置を、当該各部分画像について算出された前記位置設定ベクトルだけシフトさせ、当該シフトされた基準位置を前記位置設定ベクトルにより指定された位置とする、画像処理システム。
3. 請求項１または２に記載の画像処理システムであって、前記ピーク検出部は、前記各部分画像について、閾値以上の類似度を有する少なくとも１つのピークを検出する、画像処理システム。
4. 請求項１に記載の画像処理システムであって、前記ピーク検出部は、前記各部分画像の中から着目画素領域を順次選択するとともに、前記着目画素領域ごとに前記参照粒子像と前記着目画素領域を含む画像領域との類似度を算出することにより前記類似度分布を得る、画像処理システム。
5. 請求項１から４のうちのいずれか１項に記載の画像処理システムであって、前記２次元ピーク分布に基づいて前記対象粒子像の位置を検出する粒子像探索部をさらに備える画像処理システム。
6. 請求項５に記載の画像処理システムであって、前記粒子像探索部は、前記２次元ピーク分布に現れるピークの密度が閾値以上となる領域を探索し、その探索結果に基づいて前記対象粒子像の位置を検出する、画像処理システム。
7. 請求項６に記載の画像処理システムであって、
   前記粒子像探索部は、
   前記２次元ピーク分布に対して探索窓を走査させつつ、前記探索窓内のピークの個数に基づいて前記密度を順次算出するピーク密度算出部と、
   前記密度が閾値以上となったときの前記探索窓の位置に前記対象粒子像が位置すると判定する位置検出部と、
   を含む、画像処理システム。
8. 請求項５から７のうちのいずれか１項に記載の画像処理システムであって、
   前記粒子像探索部は、前記２次元ピーク分布に現れるピークの位置の類似度が閾値以上となる領域を探索し、その探索結果に基づいて前記対象粒子像の位置を検出し、
   前記２次元ピーク分布に現れる複数のピークの位置が重複するときは、前記複数のピークがそれぞれ有する類似度の加算結果に基づいて前記対象粒子像の位置が検出される、画像処理システム。
9. 請求項５から８のうちのいずれか１項に記載の画像処理システムであって、
   前記粒子像探索部により検出された位置にある対象粒子像を含む粒子画像を前記全体画像から抽出する粒子画像抽出部と、
   前記粒子画像または当該粒子画像から派生させた派生画像を複数のグループに分類する画像分類部と、
   当該分類された粒子画像または派生画像を前記グループごとに平均化して２次元平均化像を生成する平均化処理部と、
   前記参照粒子像の画像データを前記２次元平均化像のデータで置き換えることで前記画像記憶部の記憶内容を更新する更新処理部と、
   をさらに備え、
   前記部分画像抽出部、前記ピーク検出部、前記ベクトル算出部、前記ピーク分布生成部および前記粒子像探索部は、前記２次元平均化像中の粒子像を前記参照粒子像として使用した処理を再度実行する、画像処理システム。
10. 請求項９に記載の画像処理システムであって、前記部分画像抽出部、前記ピーク検出部、前記ベクトル算出部、前記ピーク分布生成部、前記粒子像探索部、前記粒子画像抽出部、前記平均化処理部および前記更新処理部は、繰り返し処理を実行する、画像処理システム。
11. 請求項９または１０に記載の画像処理システムであって、前記画像分類部は、粒子が向いている方向別に前記粒子画像または前記派生画像を複数のグループに分類する、画像処理システム。
12. 請求項１から１１のうちのいずれか１項に記載の画像処理システムであって、前記全体画像は、試料に含まれる複数の生体高分子を電子顕微鏡で撮像した画像である、画像処理システム。
13. 請求項１から１２のうちのいずれか１項に記載の画像処理システムであって、前記粒子画像に基づいて前記対象粒子像の立体構造を推定する立体構造推定部をさらに備える画像処理システム。
14. （ａ）複数の対象粒子像を含む全体画像から複数の部分画像を抽出するステップと、
    （ｂ）前記対象粒子像を模した参照粒子像の画像データを画像記憶部から読み出し、前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するステップと、
    （ｃ）前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するステップと、
    （ｄ）前記全体画像の領域に対して前記位置設定ベクトルを用いて指定された位置に当該ピークの位置を設定して２次元ピーク分布を生成するステップと、
    を備え、
    前記ステップ（ａ）において、前記複数の部分画像は、当該複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように前記全体画像から抽出される、画像処理方法。
15. 請求項１４に記載の画像処理方法であって、前記ステップ（ｂ）では、前記各部分画像について、閾値以上の類似度を有する少なくとも１つのピークが検出される、画像処理方法。
16. 請求項１４に記載の画像処理方法であって、
    （ｅ）前記２次元ピーク分布に基づいて前記対象粒子像の位置を検出するステップ、
    をさらに備える画像処理方法。
17. 請求項１６に記載の画像処理方法であって、前記ステップ（ｅ）は、前記２次元ピーク分布に現れるピークの密度が閾値以上となる領域を探索し、その探索結果に基づいて前記対象粒子像の位置を検出するステップを含む、画像処理方法。
18. 請求項１７に記載の画像処理方法であって、
    前記ステップ（ｅ）は、
    前記２次元ピーク分布に対して探索窓を走査させつつ、前記探索窓内のピークの個数に基づいて前記密度を順次算出するステップと、
    前記密度が閾値以上となったときに前記探索窓の位置に前記対象粒子像が位置すると判定するステップと、
    を含む、画像処理方法。
19. 請求項１６から１８のうちのいずれか１項に記載の画像処理方法であって、
    前記ステップ（ｅ）は、前記２次元ピーク分布に現れるピークの位置の類似度が閾値以上となる領域を探索し、その探索結果に基づいて前記対象粒子像の位置を検出するステップを含み、
    前記２次元ピーク分布に現れる複数のピークの位置が重複するときは、前記複数のピークがそれぞれ有する類似度の加算結果に基づいて前記対象粒子像の位置が検出される、画像処理方法。
20. 請求項１６から１９のうちのいずれか１項に記載の画像処理方法であって、
    （ｆ）前記対象粒子像を探索することにより検出された位置にある対象粒子像を含む粒子画像を前記全体画像から抽出するステップと、
    （ｇ）前記粒子画像または当該粒子画像から派生させた派生画像を複数のグループに分類し、当該分類された粒子画像または派生画像を前記グループごとに平均化して２次元平均化像を生成するステップと、
    （ｈ）前記参照粒子像の画像データを前記２次元平均化像のデータで置き換えることで前記画像記憶部の記憶内容を更新するステップと、
    をさらに備え、
    前記ステップ（ａ）〜（ｄ）は、前記２次元平均化像中の粒子像を前記参照粒子像として使用して再度実行される、画像処理方法。
21. 請求項２０に記載の画像処理方法であって、前記ステップ（ａ）〜（ｈ）は繰り返し実行される、画像処理方法。
22. コンピュータ読み取り可能な記録媒体から読み出されてコンピュータに画像処理を実行させるプログラムであって、
    前記画像処理は、
    複数の対象粒子像を含む全体画像から複数の部分画像を抽出する部分画像抽出処理と、
    前記対象粒子像を模した参照粒子像の画像データを画像記憶部から読み出し、前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するピーク検出処理と、
    前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するベクトル算出処理と、
    前記全体画像の領域に対して前記位置設定ベクトルを用いて指定された位置に当該ピークの位置を設定して２次元ピーク分布を生成するピーク分布生成処理と、
    を含み、
    前記部分画像抽出処理では、前記複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように前記全体画像から前記複数の部分画像が抽出される、プログラム。
23. 請求項２２に記載のプログラムであって、前記画像処理は、前記２次元ピーク分布に基づいて前記対象粒子像の位置を検出する粒子像探索処理をさらに含む、プログラム。
24. 請求項２３に記載のプログラムであって、
    前記画像処理は、
    前記粒子像探索処理により検出された位置にある対象粒子像を含む粒子画像を前記全体画像から抽出する粒子画像抽出処理と、
    前記粒子画像または当該粒子画像から派生させた派生画像を複数のグループに分類する画像分類処理と、
    当該分類された粒子画像または派生画像を前記グループごとに平均化して２次元平均化像を生成する平均化処理と、
    前記参照粒子像の画像データを前記２次元平均化像のデータで置き換えることで前記画像記憶部の記憶内容を更新する更新処理と、
    をさらに含み、
    前記部分画像抽出処理、前記ピーク検出処理、前記ベクトル算出処理および前記ピーク分布生成処理は、前記２次元平均化像中の粒子像を前記参照粒子像として使用して再度実行される、プログラム。
25. 請求項２４に記載のプログラムであって、前記部分画像抽出処理、前記ピーク検出処理、前記ベクトル算出処理、前記ピーク分布生成処理、前記粒子像探索処理、前記粒子画像抽出処理、前記平均化処理および前記更新処理は繰り返し実行される、プログラム。
26. コンピュータに画像処理を実行させるプログラムを記録した記録媒体であって、
    前記画像処理は、
    複数の対象粒子像を含む全体画像から複数の部分画像を抽出する部分画像抽出処理と、
    前記対象粒子像を模した参照粒子像の画像データを画像記憶部から読み出し、前記各部分画像について当該各部分画像と前記参照粒子像との間の類似度分布を算出し、前記類似度分布のピークを検出するピーク検出処理と、
    前記各部分画像について当該各部分画像内の所定の基準位置から前記ピークの位置に至る位置設定ベクトルを算出するベクトル算出処理と、
    前記全体画像の領域に対して前記位置設定ベクトルを用いて指定された位置に当該ピークの位置を設定して２次元ピーク分布を生成するピーク分布生成処理と、
    を含み、
    前記部分画像抽出処理では、前記複数の部分画像のうち所定方向に隣り合う部分画像同士が重複するように前記全体画像から前記複数の部分画像が抽出される、記録媒体。
27. 請求項２６に記載の記録媒体であって、前記画像処理は、前記２次元ピーク分布に現れるピークの密度に基づいて前記対象粒子像を探索することにより前記対象粒子像の位置を検出する粒子像探索処理をさらに含む、記録媒体。
28. 請求項２７に記載の記録媒体であって、
    前記画像処理は、
    前記粒子像探索処理により検出された位置にある対象粒子像を含む粒子画像を前記全体画像から抽出する粒子画像抽出処理と、
    前記粒子画像または当該粒子画像から派生させた派生画像を複数のグループに分類する画像分類処理と、
    当該分類された粒子画像または派生画像を前記グループごとに平均化して２次元平均化像を生成する平均化処理と、
    前記参照粒子像の画像データを前記２次元平均化像のデータで置き換えることで前記画像記憶部の記憶内容を更新する更新処理と、
    をさらに含み、
    前記部分画像抽出処理、前記ピーク検出処理、前記ベクトル算出処理および前記ピーク分布生成処理は、前記２次元平均化像中の粒子像を前記参照粒子像として使用して再度実行される、記録媒体。
29. 請求項２８に記載の記録媒体であって、前記部分画像抽出処理、前記ピーク検出処理、前記ベクトル算出処理、前記ピーク分布生成処理、前記粒子像探索処理、前記粒子画像抽出処理、前記平均化処理および前記更新処理は繰り返し実行される、記録媒体。