(57)【要約】【課題】  高いノイズレベルの画像中に埋もれたコントラストの低い解析対象物の画像位置を簡便且つ自動的に特定することを可能とする画像処理方法および画像処理プログラムを提供すること。【解決手段】  解析対象粒子が多数撮影されている画像内に、ローカルな画像を切り出すための仮想的なフレームを多数配置させ（Ｓ１００）、これらのフレームをＴＥＭ画像内でランダムに動かしながら（Ｓ１０４）、各フレームにより切り出された全画像の加算平均を求め（Ｓ１０５）、この平均画像のコントラストが高くなるように各フレームの位置移動と回転移動の双方に制約を加える（Ｓ１０７〜Ｓ１１４）。この制約は、温度パラメータＴにより制御されており、温度パラメータＴを高い値から低い値の減少させるための手順を定めるアニーリングスケジュールに従って温度を徐々に減少させる（Ｓ１１５）と、最終的には各フレームが粒子画像位置へと高い精度で配置される。【選択図】  図１
2005122284.tif 000002
【特許請求の範囲】
【請求項１】
画像内で解析対象像の位置を特定するための画像処理方法であって、
前記画像内の局所部位画像を切り出すための複数のフレームを当該画像内に配置する第１のステップと、
前記複数のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第２のステップと、
前記平均フレーム画像を基に第１のエネルギ値（Ｅ１）を算出する第３のステップと、
前記複数のフレームの中から１つのフレームを選択して当該選択されたフレームの前記画像内での位置および角度を変更してフレームを再配置する第４のステップと、
前記再配置後のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第５のステップと、
前記第５のステップで得られた平均フレーム画像を基に第２のエネルギ値（Ｅ２）を算出する第６のステップと、
前記第１および第２のエネルギ値を比較して、Ｅ１＜Ｅ２の場合には前記第４のステップのフレーム再配置を採用する一方、Ｅ１≧Ｅ２の場合には前記第４のステップのフレーム再配置を採用するか否かを判断する第７のステップと、
前記第１〜第７のステップを所定回数だけ繰り返した後に温度パラメータＴを減少させる第８のステップと、
前記第１〜第８のステップを所定のサイクル数だけ繰り返す第９のステップと、
を備えていることを特徴とする画像処理方法。
【請求項２】
前記第７のステップにおいて、Ｅ１≧Ｅ２の場合に前記第４のステップのフレーム再配置を採用するか否かの判断は、
前記第１のエネルギ値Ｅ１と前記第２のエネルギ値Ｅ２との差（ΔＥ）から次式により採用確率Ｐ（ΔＥ）を求めて乱数との比較を行い、前記採用確率Ｐ（ΔＥ）が当該乱数よりも大きい場合は前記第４のステップにおけるフレーム再配置を採用する一方、前記採用確率Ｐ（ΔＥ）が前記乱数以下であった場合は第４のステップにおけるフレーム再配置前のフレーム配置を採用することにより実行されることを特徴とする請求項１に記載の画像処理方法。
【数１】
2005122284.tif 000021
【請求項３】
前記第１および第２のエネルギ値は、前記画像の実画像に基づいて算出されるエネルギ値（実画像のエネルギ値Ereal）または実画像をＦＦＴして得られた画像に基づいて得られる周波数領域でのエネルギ値（周波数領域のエネルギ値Epow）の何れか一方または前記実画像のエネルギ値Erealと前記周波数領域のエネルギ値Epowの和（Ereal＋Epow）であることを特徴とする請求項１または２に記載の画像処理方法。
【請求項４】
前記実画像のエネルギ値Erealは、
前記平均フレーム画像内に設けられた（x,y）座標の各座標点での輝度d(x,y)を基に次式によりピクセル輝度平均値daveを求め、前記平均フレーム画像d(x,y)からピクセル輝度平均値daveを引いた値d’(x,y)を求め、さらに、当該d’(x,y)値を次式に代入して求められるものであることを特徴とする請求項３に記載の画像処理方法。
【数２】
2005122284.tif 000022
【数３】
2005122284.tif 000023
【数４】
2005122284.tif 000024
（NとMはフレームのx軸およびy軸のピクセルサイズ、ＨはErealの値が０．５となるときの輝度値）
【請求項５】
前記周波数領域のエネルギ値Epowは、前記平均フレーム画像内に設けられた（x,y）座標の各座標点での輝度d(x,y)を基に次式によりピクセル輝度平均値daveを求め、前記平均フレーム画像d(x,y)からピクセル輝度平均値daveを引いた値d’(x,y)を求め、さらに、当該d’(x,y)値からＦＦＴによりそのスペクトルの振幅｜D’(u,v)｜を求め、当該スペクトルの振幅｜D’(u,v)｜に透過周波数を規定する円形マスクＦ（u,v）を掛け、周波数領域における座標（u,v）での各ピクセルにつき次式により総和をとることにより求められるものであることを特徴とする請求項３に記載の画像処理方法。
【数５】
2005122284.tif 000025
【数６】
2005122284.tif 000026
【数７】
2005122284.tif 000027
【数８】
2005122284.tif 000028
（NとMはフレームのx軸およびy軸のピクセルサイズ）
【請求項６】
前記第８のステップにおける温度パラメータＴの減少率Ｒは、次式で与えられるものであることを特徴とする請求項１乃至５の何れかに記載の画像処理方法。
【数９】
2005122284.tif 000029
（Ｔ０：初期温度、Ｔｅ：演算を終了させる最終温度、Ｃmax：第８のステップにおける所定回数）
【請求項７】
画像内で解析対象像の位置を特定するための画像処理をコンピュータに実行させるための画像処理プログラムであって、
前記画像内の局所部位画像を切り出すための複数のフレームを当該画像内に配置する第１のステップと、
前記複数のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第２のステップと、
前記平均フレーム画像を基に第１のエネルギ値（Ｅ１）を算出する第３のステップと、
前記複数のフレームの中から１つのフレームを選択して当該選択されたフレームの前記画像内での位置および角度を変更してフレームを再配置する第４のステップと、
前記再配置後のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第５のステップと、
前記第５のステップで得られた平均フレーム画像を基に第２のエネルギ値（Ｅ２）を算出する第６のステップと、
前記第１および第２のエネルギ値を比較して、Ｅ１＜Ｅ２の場合には前記第４のステップのフレーム再配置を採用する一方、Ｅ１≧Ｅ２の場合には前記第４のステップのフレーム再配置を採用するか否かを判断する第７のステップと、
前記第１〜第７のステップを所定回数だけ繰り返した後に温度パラメータＴを減少させる第８のステップと、
前記第１〜第８のステップを所定のサイクル数だけ繰り返す第９のステップと、
を備えていることを特徴とする画像処理プログラム。
【請求項８】
前記第７のステップにおいて、Ｅ１≧Ｅ２の場合に前記第４のステップのフレーム再配置を採用するか否かの判断は、
前記第１のエネルギ値Ｅ１と前記第２のエネルギ値Ｅ２との差（ΔＥ）から次式により採用確率Ｐ（ΔＥ）を求めて乱数との比較を行い、前記採用確率Ｐ（ΔＥ）が当該乱数よりも大きい場合は前記第４のステップにおけるフレーム再配置を採用する一方、前記採用確率Ｐ（ΔＥ）が前記乱数以下であった場合は第４のステップにおけるフレーム再配置前のフレーム配置を採用することにより実行されることを特徴とする請求項７に記載の画像処理プログラム。
【数１０】
2005122284.tif 000030
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像処理方法および画像処理プログラムに関し、特に、高いノイズレベルの画像中に埋もれたコントラストの低い解析対象物の画像位置を簡便且つ自動的に特定することを可能とする画像処理方法および画像処理プログラムに関する。
【背景技術】
【０００２】
従来、タンパク質の構造を解析するに際しては、解析対象であるタンパク質を結晶化させて試料作成し、この試料をＸ線結晶構造解析法とNMR法とにより解析・評価する手法が主としてとられてきたが、近年では、極低温透過型電子顕微鏡（Cryo-EM）を用いてタンパク質の構造解析を行う方法（単粒子構造解析法）が開発されつつある（例えば、非特許文献1参照）。
【０００３】
この単粒子構造解析法は、透過型電子顕微鏡（ＴＥＭ）中で極低温に冷却した粒子状のタンパク質試料から得られたＴＥＭ画像から、個々の球状タンパク質試料の画像を切り出し、これらの試料画像を加算平均することで試料画像のノイズを低減させ、得られた加算平均画像から粒子状タンパク質試料の３次元構造を求める手法であり、近年の画像処理技術の高度化とコンピュータの高速化を背景に、有力・有用な解析手法として注目されてきている。
【０００４】
【非特許文献１】W. V. Nicholson and R. M.Glaeser, “Review:Automatic Particle Detection in Electron Microscopy”, Journal ofStructural Biology, Vol.133, pp.90-101 (2001).
【発明の開示】
【発明が解決しようとする課題】
【０００５】
従来の単粒子構造解析法で画像処理を行うに際しては、電子顕微鏡画像中に撮影されているタンパク質の粒子画像の各々を解析者自らが目視によりマニュアルで切り出し、それらの画像を加算平均するという作業が求められる。
【０００６】
しかしながら、極低温Cryo-EMで得られるＴＥＭ画像のコントラストは低くならざるを得ないことに加えノイズレベルも高いため、タンパク質粒子の画像中での位置を特定することは経験を積んだ研究者にとっても難しい作業であるために迅速な解析が困難であり、さらに、解析者の熟練度や経験によって得られる結果の個人差が大きくならざるを得ないという問題があった。このことが、単粒子構造解析法の一般的な解析手法としての普及を阻害する大きな要因の一つとなっていた。
【０００７】
本発明はこのような問題に鑑みてなされたもので、その目的とするところは、高いノイズレベルの画像中に埋もれたコントラストの低い解析対象物の画像位置を簡便且つ自動的に特定することを可能とする画像処理方法および画像処理プログラムを提供することにある。
【課題を解決するための手段】
【０００８】
本発明は、このような目的を達成するために、第１の発明は、画像内で解析対象像の位置を特定するための画像処理方法であって、前記画像内の局所部位画像を切り出すための複数のフレームを当該画像内に配置する第１のステップと、前記複数のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第２のステップと、前記平均フレーム画像を基に第１のエネルギ値（Ｅ１）を算出する第３のステップと、前記複数のフレームの中から１つのフレームを選択して当該選択されたフレームの前記画像内での位置および角度を変更してフレームを再配置する第４のステップと、前記再配置後のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第５のステップと、前記第５のステップで得られた平均フレーム画像を基に第２のエネルギ値（Ｅ２）を算出する第６のステップと、前記第１および第２のエネルギ値を比較して、Ｅ１＜Ｅ２の場合には前記第４のステップのフレーム再配置を採用する一方、Ｅ１≧Ｅ２の場合には前記第４のステップのフレーム再配置を採用するか否かを判断する第７のステップと、前記第１〜第７のステップを所定回数だけ繰り返した後に温度パラメータＴを減少させる第８のステップと、前記第１〜第８のステップを所定のサイクル数だけ繰り返す第９のステップと、を備えていることを特徴とする。
【０００９】
好ましくは、前記第７のステップにおいて、Ｅ１≧Ｅ２の場合に前記第４のステップのフレーム再配置を採用するか否かの判断は、前記第１のエネルギ値Ｅ１と前記第２のエネルギ値Ｅ２との差（ΔＥ）から次式により採用確率Ｐ（ΔＥ）を求めて乱数との比較を行い、前記採用確率Ｐ（ΔＥ）が当該乱数よりも大きい場合は前記第４のステップにおけるフレーム再配置を採用する一方、前記採用確率Ｐ（ΔＥ）が前記乱数以下であった場合は第４のステップにおけるフレーム再配置前のフレーム配置を採用することにより実行される。
【００１０】
【数１】
2005122284.tif 000003
【００１１】
また、好ましくは、前記第１および第２のエネルギ値は、前記画像の実画像に基づいて算出されるエネルギ値（実画像のエネルギ値Ereal）または実画像をＦＦＴして得られた画像に基づいて得られる周波数領域でのエネルギ値（周波数領域のエネルギ値Epow）の何れか一方または前記実画像のエネルギ値Erealと前記周波数領域のエネルギ値Epowの和（Ereal＋Epow）である。
【００１２】
また、好ましくは、前記実画像のエネルギ値Erealは、前記平均フレーム画像内に設けられた（x,y）座標の各座標点での輝度d(x,y)を基に次式によりピクセル輝度平均値daveを求め、前記平均フレーム画像d(x,y)からピクセル輝度平均値daveを引いた値d’(x,y)を求め、さらに、当該d’(x,y)値を次式に代入して求められるものである。
【００１３】
【数２】
2005122284.tif 000004
【００１４】
【数３】
2005122284.tif 000005
【００１５】
【数４】
2005122284.tif 000006
【００１６】
（NとMはフレームのx軸およびy軸のピクセルサイズ、ＨはErealの値が０．５となるときの輝度値）
【００１７】
また、好ましくは、前記周波数領域のエネルギ値Epowは、前記平均フレーム画像内に設けられた（x,y）座標の各座標点での輝度d(x,y)を基に次式によりピクセル輝度平均値daveを求め、前記平均フレーム画像d(x,y)からピクセル輝度平均値daveを引いた値d’(x,y)を求め、さらに、当該d’(x,y)値からＦＦＴによりそのスペクトルの振幅｜D’(u,v)｜を求め、当該スペクトルの振幅｜D’(u,v)｜に透過周波数を規定する円形マスクを掛け、周波数領域における座標（u,v）での各ピクセルにつき次式により総和をとることにより求められるものである。
【００１８】
【数５】
2005122284.tif 000007
【００１９】
【数６】
2005122284.tif 000008
【００２０】
【数７】
2005122284.tif 000009
【００２１】
【数８】
2005122284.tif 000010
【００２２】
（NとMはフレームのx軸およびy軸のピクセルサイズ、ｒは透過周波数を規定する円形マスクの半径）
【００２３】
さらに、好ましくは、前記第８のステップにおける温度パラメータＴの減少率Ｒは、次式で与えられるものである。
【００２４】
【数９】
2005122284.tif 000011
【００２５】
（Ｔ０：初期温度、Ｔｅ：演算を終了させる最終温度、Ｃmax：第８のステップにおける所定回数）
【００２６】
第２の発明は、画像内で解析対象像の位置を特定するための画像処理をコンピュータに実行させるための画像処理プログラムであって、前記画像内の局所部位画像を切り出すための複数のフレームを当該画像内に配置する第１のステップと、前記複数のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第２のステップと、前記平均フレーム画像を基に第１のエネルギ値（Ｅ１）を算出する第３のステップと、前記複数のフレームの中から１つのフレームを選択して当該選択されたフレームの前記画像内での位置および角度を変更してフレームを再配置する第４のステップと、前記再配置後のフレームにより切り出された画像を平均化処理して平均フレーム画像を形成する第５のステップと、前記第５のステップで得られた平均フレーム画像を基に第２のエネルギ値（Ｅ２）を算出する第６のステップと、前記第１および第２のエネルギ値を比較して、Ｅ１＜Ｅ２の場合には前記第４のステップのフレーム再配置を採用する一方、Ｅ１≧Ｅ２の場合には前記第４のステップのフレーム再配置を採用するか否かを判断する第７のステップと、前記第１〜第７のステップを所定回数だけ繰り返した後に温度パラメータＴを減少させる第８のステップと、前記第１〜第８のステップを所定のサイクル数だけ繰り返す第９のステップと、を備えていることを特徴とする。
【００２７】
好ましくは、前記第７のステップにおいて、Ｅ１≧Ｅ２の場合に前記第４のステップのフレーム再配置を採用するか否かの判断は、前記第１のエネルギ値Ｅ１と前記第２のエネルギ値Ｅ２との差（ΔＥ）から次式により採用確率Ｐ（ΔＥ）を求めて乱数との比較を行い、前記採用確率Ｐ（ΔＥ）が当該乱数よりも大きい場合は前記第４のステップにおけるフレーム再配置を採用する一方、前記採用確率Ｐ（ΔＥ）が前記乱数以下であった場合は第４のステップにおけるフレーム再配置前のフレーム配置を採用することにより実行される。
【００２８】
【数１０】
2005122284.tif 000012
【発明の効果】
【００２９】
本発明の画像処理方法では、解析対象粒子が多数撮影されている画像内にローカルな画像を切り出すための仮想的なフレームを多数配置させ、これらのフレームをＴＥＭ画像内でランダムに動かしながら各フレームにより切り出された全画像の加算平均を求め、この平均画像のコントラストが高くなるように各フレームの位置移動と回転移動の双方に制約を加える。この制約は、温度パラメータＴにより制御されており、温度パラメータＴを高い値から低い値の減少させるための手順を定めるアニーリングスケジュールに従って温度を徐々に減少させることとしたので、高いノイズレベルの画像中に埋もれたコントラストの低い解析対象物の画像位置を簡便且つ自動的（full-auto）に特定することを可能とする画像処理方法および画像処理プログラムを提供することが可能となる。また、特定された位置の画像を加算平均することにより、解析対象物の形状を高精度で再現することが可能となる。
【発明を実施するための最良の形態】
【００３０】
以下に、図面を参照して本発明の実施の形態について説明する。なお、以下の実施例においては、ノイズレベルの高い画像中に埋もれたコントラストの弱い解析対象物の画像の好適な例として、透過型電子顕微鏡画像（ＴＥＭ画像）中にノイズに埋もれた状態で撮影されているコントラストが微弱なタンパク質粒子像を取り上げ、本発明の画像処理方法および画像処理プログラムを説明するが、本発明はこのような電子顕微鏡画像の処理に限定されるものではなく、例えば、航空写真中に撮影されている車両位置の特定などの一般的な適用が可能であることはいうまでもない。
【００３１】
本発明の画像処理方法は、組み合わせ最適化問題の一つの解法として知られるアニーリング（焼きなまし）法（Simulated Annealing）に基づくアルゴリズムに従って粒子画像の自己集積化を実行する手法である。本発明の方法では、解析対象である粒子が多数撮影されている１枚のＴＥＭ画像内に、ローカルな画像を切り出すための仮想的なフレームを多数（例えば数１００個）配置させる。これらのフレームをＴＥＭ画像内でランダムに動かしながら、各フレームにより切り出された全画像の加算平均を求め、この平均画像のコントラストが高くなるように各フレームの位置移動と回転移動の双方に制約を加える。この制約は、温度パラメータ（Ｔ）により制御されており、温度パラメータＴを高い値から低い値の減少させるための手順を定めるアニーリングスケジュール（annealing
schedule）に従って温度を徐々に減少させると、最終的には各フレームが粒子画像位置へと高い精度で配置される。そして、このようにして配置された各フレームにより切り出された画像を加算平均して得られる画像は、解析対象である粒子の形状を高い精度で再現することができる。
【００３２】
本発明の画像処理方法では、ローカル画像の切り出しに用いるフレームの数と大きさ（半径）のみを予め設定するだけで画像処理が実行可能であり、従来法で要求されていた解析対象粒子に関する形状因子等の特別な情報に関しては不要となる。従って、テンプレートや学習用画像データなしでの画像処理が可能となるから、処理の完全自動化を図ることが可能となる。
【００３３】
以下に、本発明の画像処理方法の手順を具体的に説明する。
【００３４】
図１は本発明の画像処理方法を実行する際のシミュレーテッド・アニーリング・アルゴリズム（ＳＡアルゴリズム）を説明するためのフローチャートであり、図２は解析対象がタンパク質粒子である場合の初期設定状態でのＴＥＭ像（図２（ａ））、ＳＡアルゴリズムによる位置推定プロセス（図２（ｂ）））およびタンパク質粒子像位置の推定結果（図２（ｃ）））を説明するための各画像である。
【００３５】
本発明の画像処理方法では、図２に示すようにＴＥＭ画像内にローカルな部位の画像を切り出すためのフレーム（図中のＦ１〜Ｆ５）を配置し、これらのフレームにより切り出された画像を平均化して得られる平均フレーム画像の「スコア」が高くなるように、図１に示したＳＡアルゴリズムに従って各フレームのＴＥＭ画像内での位置と回転角度を変化させる。ここで、「スコア」とは、後述する手順により算出される「系のエネルギ値」であり、処理の妥当性の指標となるものである。
【００３６】
先ず、個々のフレームの大きさは単一のタンパク質粒子の像が入る大きさに設定され、各フレームの初期位置は各々のフレームが互いに重ならないように、ＴＥＭ画像中にランダムに配置される。なお、各フレームの回転角度の初期値は全て０°としている。また、後述するフレームの位置・角度の変更・非変更回数のカウンタ（変更回数カウンタIte1、および、非変更回数カウンタIte2）の最大数（Itemax1およびItemax2）を設定し、さらに、温度減少サイクルの繰り返し回数Ｃmaxを設定する（ステップＳ１００）。
【００３７】
次に、図２（ｂ）に示すように、ＴＥＭ画像内に配置されたフレームによりローカルな画像を切り出し、これらの画像を平均化して平均フレーム画像を求め、この平均フレーム画像をもとに下記の要領で系の初期エネルギ値を計算する（ステップＳ１０１）。ここで、エネルギ値としては、実画像に基づいて算出されるエネルギ値（実画像のエネルギ値Ereal）または実画像をＦＦＴして得られた画像に基づいて得られる周波数領域でのエネルギ値（周波数領域のエネルギ値Epow）の何れのエネルギ値を用いてもよいが、ここでは、これらのエネルギ値の和（Ereal＋Epow）を系のエネルギ値（E）とした。
【００３８】
実画像のエネルギ値Erealの算出は以下の手順により実行する。先ず、平均フレーム画像内に（x,y）座標を設けて各座標点での輝度d(x,y)を求め、これらの輝度d(x,y)を基に、次式（１）によりピクセル輝度平均値daveを求める。
【００３９】
【数１１】
2005122284.tif 000013
【００４０】
ここで、NとMはフレームのx軸およびy軸のピクセルサイズである。
【００４１】
次に、平均フレーム画像d(x,y)からピクセル輝度平均値daveを引いた値d’(x,y)を求め（式２）、このd’(x,y)値を飽和型関数である式（３）のMichaeals-Meaintai関数に代入してその総和を実画像のエネルギ値Erealとした。
【００４２】
【数１２】
2005122284.tif 000014
【００４３】
【数１３】
2005122284.tif 000015
【００４４】
ここで、Ｈは飽和型関数の値が０．５となるときの輝度値である。このような取扱いにより、一般的に用いられる差の２乗和のように、輝度値が高くなるにつれて急速にスコアが上昇し平均フレーム画像の一部分の輝度値のみが高くなる、という不都合を回避することが可能となる。
【００４５】
一方、周波数領域のエネルギ値Epowの算出は、上式（２）で求めた値d’(x,y)からＦＦＴによりそのスペクトルの振幅｜D’(u,v)｜を求め、このスペクトルの振幅｜D’(u,v)｜に透過周波数を規定する円形マスクを掛け、周波数領域における座標（u,v）での各ピクセルにつき次式（４）により総和をとることにより実行する。ここで、円形マスクF(u,v)は次式（５）で定義される関数であり、円形マスクの大きさ（半径ｒ）は中央位置（um,vm）から半径ｒ＝６．４（pixel）に設定した。
【００４６】
【数１４】
2005122284.tif 000016
【００４７】
【数１５】
2005122284.tif 000017
【００４８】
このような初期設定に続き、ステップＳ１０２において変更回数カウンタIte1＝０、非変更回数カウンタIte2＝０と設定し、ＴＥＭ画像上のフレームから１つをランダムに選択し（ステップＳ１０３）、このフレームの位置と角度とを変更する（ステップＳ１０４）。ここでは、変更前後での変化量を、位置についてはx軸、ｙ軸ともに７σの正規乱数とし、角度については３０σの正規乱数とした。
【００４９】
上述した計算に基づいて、位置・角度変更後の平均フレーム画像のエネルギ値Ｅ２を算出して変更前のエネルギ値Ｅ１との差（ΔE＝Ｅ１−Ｅ２）を求め（ステップＳ１０５）、変更後のエネルギ値と変更前のエネルギ値とを比較して、変更後のエネルギ値が変更前のエネルギ値よりも高いかどうかの判断を行う（ステップＳ１０６）。ここで、変更後のエネルギ値が変更前のエネルギ値よりも高い場合（ステップＳ１０６：Ｙｅｓ）は、スコアを変更後のエネルギ値に更新し（ステップＳ１０７）、ステップＳ１０４でのフレームの位置・角度変更を採用して保持して（ステップＳ１０８）、変更回数カウンタIte1を１増加させる（ステップＳ１０９）。
【００５０】
逆に、変更後のエネルギ値が変更前のエネルギ値以下の場合（ステップＳ１０６：Ｎｏ）は、変更前のスコアとの変化量（ΔE）から次式（６）により採用確率（遷移確率）P(ΔE)を求め（ステップＳ１１０）、０〜１の乱数との比較（P(ΔE)＞（乱数０〜１））を行い（ステップＳ１１１）、乱数が採用確率P(ΔE)未満の場合（ステップＳ１１１：Ｙｅｓ）はステップＳ１０４におけるフレームの変更を採用し、乱数が採用確率P(ΔE)以上であった場合（ステップＳ１１１：Ｎｏ）は、フレームの位置・角度を変更前の場所・値に戻して（ステップＳ１１２）、非変更回数カウンタIte2を１増加させる（ステップＳ１１３）。
【００５１】
【数１６】
2005122284.tif 000018
【００５２】
ステップＳ１１４では、変更回数カウンタIte1の値および非変更回数カウンタIte2の値をそれぞれの最大数（Itemax1およびItemax2）と比較し、「Ite1＜Itemax1かつIte2＜Itemax2」であれば（ステップＳ１１４：Ｙｅｓ）ステップＳ１０３に戻り以降のステップを繰り返し、「Ite1＜Itemax1かつIte2＜Itemax2」でなければ（ステップＳ１１４：Ｎｏ）ステップＳ１１５に移行して系の温度パラメータＴを減少させる。
【００５３】
温度パラメータＴの減少率Ｒは次式（７）で与えられる。ここで、Ｔ０は初期温度、Ｔｅは演算を終了させる最終温度、そして、Ｃmaxは所定のサイクル回数である。すなわち、温度パラメータＴの減少率Ｒは、所定のサイクル回数Ｃmaxに従って次式（７）により設定され、温度パラメータＴは次式（８）により各サイクル毎に指数関数的に減少する。
【００５４】
【数１７】
2005122284.tif 000019
【００５５】
【数１８】
2005122284.tif 000020
【００５６】
この温度減少を実行させる条件としては、例えば、フレームの変更回数Ite1と非変更回数Ite2をそれぞれ集計し、変更回数Ite1が全フレーム数の２０倍となるか、または、非変更回数Ite2が全フレーム数の１００倍となった場合とする。
【００５７】
１回の温度減少を行う毎にサイクルカウンタCを１増加させ（ステップＳ１１６）、サイクルがステップＳ１００で設定した所定回数Cmax以上かどうかを判断（ステップＳ１１７）して、所定回数以上であれば（ステップＳ１１７：Ｙｅｓ）終了し、所定回数に満たなければ（ステップＳ１１７：Ｎｏ）、ステップＳ１０２に移行して変更回数カウンタIte1＝０、非変更回数カウンタIte2＝０と設定し直した上で、ステップＳ１０３以降のステップを繰り返す。
【００５８】
このような一連のフレーム配置操作が終了すると、全てのフレームがタンパク質粒子像（Ｔ１〜Ｔ５）の上に配置されることとなり（図２（ｃ））、解析対象であるタンパク質粒子像の存在位置が特定されるとともに、これらのフレームにより切り出された画像を平均化して得られる平均画像はタンパク質粒子像を高い精度で再現することとなる。
【実施例１】
【００５９】
図３は、本発明の画像処理方法の処理プログラムを実行可能な電子計算機により、タンパク分子のモデル画像の画像処理を実行した場合の温度減少のサイクル回数毎のフレーム配置を説明するための図である。なお、これらの図では、タンパク質粒子の位置認識状態の理解を容易にするために、粒子が存在するフレームを白で、存在しないフレームを黒で示した。
【００６０】
サイクル回数０（図３（ａ））では、フレームはランダムに配置されており、殆どのフレーム内にはタンパク質粒子は存在していない。サイクル回数が１２０（図３（ｂ））、１５０（図３（ｃ））、３００（図３（ｄ））と多くなるにつれて、タンパク質粒子位置に配置されるサークルは徐々に増加し、サイクル回数３００では全てのフレームが粒子位置に配置されている。
【実施例２】
【００６１】
図４は、オリジナル画像にノイズを加えた場合の、粒子位置認識レベルのノイズ含有量依存性を説明するための図である。
【００６２】
図４（ａ）はオリジナル画像であり、図４（ｂ）〜（ｄ）はこのオリジナル画像にノイズを含有させた合成画像である。分散σ＝２０の正規ランダムノイズを加えた場合（図４（ｂ））は、ほぼ全てのフレーム内に粒子が存在している。ノイズ含有量を増やし分散σ＝３０の正規ランダムノイズを加えると（図４（ｃ））、フレームが粒子位置に配置される確率は約８０％に低下し、さらに、分散σ＝４０の正規ランダムノイズを加えた場合（図４（ｃ））には、フレームが粒子位置に配置される確率は５０％程度に低下する。すなわち、粒子位置認識レベルは画像に含まれるノイズレベルに依存することとなる。
【００６３】
しかしながら、図４（ｄ）のようなノイズレベルの高い画像から目視により粒子位置を特定することは殆ど不可能であるのに対して、本発明の画像処理方法によれば、５０％程度の確率で粒子位置を自動的に特定することが可能である。
【実施例３】
【００６４】
図５は、２%の劣化ウラン溶液によるNegative stain処理を施したNaチャネルのＴＥＭ画像に対する本発明の画像処理方法の適用例を説明する図である。図５（ａ）は上述したアルゴリズムに基づいてフレーム配置を実行した最終結果の画像であり、略全てのフレームが粒子位置に配置されている。図５（ｂ）は温度減少のサイクル回数毎のフレーム平均画像を示しており、初期画像（サイクル回数＝０）ではノイズのみで粒子の形状を確認することはできないが、サイクル回数７０のフレーム平均画像中には粒子形状が確認できるようになり、サイクル回数を１２０、１５０、３００と増加させるにつれて、徐々に本来のNaチャネルに近い形状が生成されていく。図５（ｃ）はサイクル回数毎の温度およびスコア（系のエネルギ値）の変化を示しており、温度が指数関数的に低下するに伴ってスコア（エネルギ）が増加し、処理の妥当性が高くなる様子が理解できる。
【実施例４】
【００６５】
図６は、アポフェリチン（apoferritin）タンパク質の極低温Cryo-EM画像に本発明の方法による画像処理を施した結果を説明するための図で、図６（ａ）は上述したアルゴリズムに基づいてフレーム配置を実行した最終結果の画像であり、略全てのフレームが粒子位置に配置されている。図６（ｂ）は温度減少のサイクル回数毎のフレーム平均画像を示しており、初期画像（サイクル回数＝０）ではノイズのみで粒子の形状を確認することはできないが、サイクル回数７０のフレーム平均画像中には構造が認識可能な程度の粒子形状が確認できるようになり、サイクル回数を１２０、１５０、３００と増加させるにつれて、徐々に本来のアポフェリチンタンパク質に近い形状が生成されていく。図６（ｃ）はサイクル回数毎の温度およびスコア（系のエネルギ値）の変化を示しており、温度が指数関数的に低下するに伴ってスコアが増加し、処理の妥当性が高くなる様子が理解できる。
【００６６】
このように、本発明の画像処理方法によれば、コントラストの低い解析対象画像が高いノイズレベルに埋もれた画像であっても、タンパク粒子の画像中での位置を的確に捉え、その形状を高い精度で再現することができる。
【００６７】
上述したように、本発明の方法をタンパク質解析に応用する場合を想定すると、本発明の画像処理方法で実行している内容は、タンパク質が結晶化してゆく状態をコンピュータ上で模擬的に再現することに類似している。従って、タンパク質の結晶化の際には、ランダム配列をしていたタンパク分子がより安定な構造（系の最小エネルギ状態）となるように規則的に集積してタンパク質結晶となるのと類似して、各フレームにより切り出された画像の平均画像は、ノイズに埋もれて形状が不明確な当初の状態から徐々に粒子状構造が再現されてゆくこととなる。
【００６８】
このように、本発明の画像処理方法では、ノイズに埋もれた状態で存在している解析対象物の形状を予め知る必要はなく、後述する自己集積する幾つかのパラメータを設定するだけで任意の形状の粒子を抽出（pick-up）することが可能である。また、本発明の画像処理方法はノイズレベルの高い画像についても適用可能であることから、目視レベルでは認識不可能な不鮮明な粒子画像であっても、高精度かつ簡便に抽出することができる。
【産業上の利用可能性】
【００６９】
本発明は、高いノイズレベルの画像中に埋もれたコントラストの低い解析対象物の画像位置を簡便且つ自動的に特定することを可能とする画像処理方法および画像処理プログラムを提供する。
【図面の簡単な説明】
【００７０】
【図１】本発明の画像処理方法を実行する際のシミュレーテッド・アニーリング・アルゴリズム（ＳＡアルゴリズム）を説明するためのフローチャートである。
【図２】解析対象がタンパク質粒子である場合の初期設定状態でのＴＥＭ像（ａ）、ＳＡアルゴリズムによる位置推定プロセス（ｂ）およびタンパク質粒子像位置の推定結果（ｃ）を説明するための各画像である。
【図３】本発明の画像処理方法の処理プログラムを実行可能な電子計算機により、タンパク分子のモデル画像の画像処理を実行した場合の繰り返し回数毎のフレーム配置を説明するための図である。
【図４】オリジナル画像にノイズを加えた場合の、粒子位置認識レベルのノイズ含有量依存性を説明するための図である。
【図５】２%の劣化ウラン溶液によるNegative stain処理を施したNaチャネルのＴＥＭ画像に対する本発明の画像処理方法の適用例を説明する図である。
【図６】アポフェリチン（apoferritin）タンパク質の極低温無染色Cryo-EM画像に本発明の方法による画像処理を施した結果を説明するための図である。
【符号の説明】
【００７１】
Ｆ１〜Ｆ５  フレーム
Ｔ１〜Ｔ５  タンパク質粒子
【図１】
2005122284.tif 000031
【図２】
2005122284.tif 000032
【図３】
2005122284.tif 000033
【図４】
2005122284.tif 000034
【図５】
2005122284.tif 000035
【図６】
2005122284.tif 000036