JP2006202115A - 画像処理方法および画像処理プログラム - Google Patents

Abstract

【課題】 立体構造が未知の対象物の３次元構造を高精度で再構築することができる画像処理方法を提供すること。
【解決手段】 ２次元の平均画像から３次元の初期構造モデルを作製し、平均画像と再投射画像との直接的な相関値を最適化の尺度として用いて３次元構造を最適化する。先ず、Ｎ枚の２次元画像を仮想球面上の初期画像位置にランダムに配置して初期３次元構造を構築し、これを初期画像位置に再投射して新たな２次元画像を得る。そして、再投射前後の２次元画像の相互相関を求めてスコア（Ｌ_１）を算出する。再投射後の１枚の２次元画像を選択してその座標パラメータを変化させて３次元構造を再構築し、仮想球面上への再投射により新たな２次元画像を得る。そして、新たなスコア（Ｌ_２）を算出して座標パラメータの変化の許容もしくは拒絶を判断する。これらのステップを繰り返して３次元構造を最適化する。
【選択図】 図１

Description

本発明は、２次元画像から３次元構造を再構築する画像処理方法に関し、より詳細には、角度情報が未知の２次元画像から３次元構造を再構築することを可能とする画像処理方法およびこれを実行するための画像処理プログラムに関する。

２次元平面に投影された立体構造を有する観察対象の像（２次元画像）を基に、当該観察対象の３次元構造を再構築するための画像処理手法は、薬学、天文学、あるいは電子顕微鏡による構造解析といった分野などで広範囲に用いられている。

たとえば、タンパク質の構造解析分野では、従来は、主に、予め結晶化されたタンパク質の試料を作製し、この結晶試料をＸ線回折法により構造解析する方法が用いられてきたが、近年、タンパク質の単粒子の電子顕微鏡画像からその３次元構造を再構築する手法（単粒子構造解析法）が提唱されるようになってきた（非特許文献１などを参照）。

単粒子構造解析法では、極低温透過型電子顕微鏡（Cryo-EM）を用いて撮影されたタンパク質試料の電子顕微鏡画像（２次元画像）から数千ないし数万の単粒子画像を切り出し、これらの２次元画像の実空間（３次元空間）での位置や角度を推定し、それぞれに適合した位置や角度で画像を加算平均してバックノイズを低減させる。そして、この加算平均画像から、それぞれの２次元画像に映し出されたタンパク質粒子の３次元角度を推定し、その３次元構造が再構築される。

このように、単粒子構造解析法では、タンパク質試料を結晶化させるといった煩雑な準備を必要とせず、精製されたタンパク質さえあればその３次元構造を決定することができるという大きな利点をもつ。また、このような３次元構造の再構築手法は、タンパク質の構造解析に限らず、２次元平面に投影された立体構造を有する観察対象の像を基に、当該観察対象の３次元構造を再構築するための画像処理手法として広範囲の分野で用いることができ、特に、近年の画像処理技術の高度化やコンピュータの高速化を背景として、注目されてきている。

I. I. Serysheva et al., "Electron cryomicroscopy and angular reconstitution used to visualizethe skeletal muscle calcium release channel", Structural Biology, volume 2 number 1 January 1995, p.18-24. Martin VAN HEEL, "ANGULAR RECONSTRUCTION: A POSTERIORI ASSIGNMENT OF PROJECTIONDIRECTIONS FOR 3D RECONSTRUCTION",Ultramicroscopy 21 (1987) p.111-124.

上述したように、単粒子構造解析では、立体構造を有する単粒子を平面上に投影した２次元画像の加算平均画像を求め、この平均画像がもつ本来の３次元空間での角度（３次元角度）を推定して３次元構造を再構築する。より詳細には、個々の平均画像の３次元角度を推定し、広域増強フィルタを掛けた後に３次元投射面に再投射させ、全平均画像の３次元投射面への再投射の総和を求めることで３次元構造が再構築される。

ここで、個々の平均画像は３次元投射面にランダムに分散しているため、個々の単粒子画像の３次元投射面への投射角度は未知である。このため、投射角度を推定するいくつかの方法が提案されている。

投影角度の推定方法は、単粒子の電子顕微鏡像を撮影する際に、観察試料を固定させるメッシュを保持するための試料ステージ（試料ホルダー）を、電子顕微鏡体内で傾斜させて撮影した単粒子画像を処理する方法と、試料傾斜させることなく撮影して得られた単粒子画像を処理する方法の、２つに大別される。

試料傾斜させて単粒子画像を撮影する第１の方法は、同一の粒子を複数の異なる角度から観察した画像を得ることができるため、単粒子構造の非対称構造（たとえば、右巻きか左巻きか）までもが特定でき、その推定精度も高い。

しかしながら、低温冷却試料からの画像を撮影する必要があるような場合には、同一の単粒子に複数回の露光を行う結果となるため、電子線により損傷を受けて構造変化を生じてしまうという問題がある。また、低温冷却用試料ステージを傾斜させるためには機械的に複雑なアタッチメントを必要とし、試料ドリフトなどが生じやすいために安定した電子顕微鏡像を得ることは困難である。このような理由により、固定ステージを用い、試料傾斜させることなく単粒子画像を撮影して単粒子画像を処理する第２の方法の方が汎用性が高く望ましい。

ところが、試料傾斜させることなく単粒子画像を撮影して単粒子画像を処理する場合には、得られた粒子画像の３次元角度（オイラー角）の推定は困難であり、高分解能の３次元構造をいきなり求めることは一般的には困難である。したがって、先ずラフな３次元モデル（初期の３次元構造）を構成しておいて、この３次元モデルの精度を徐々に高めてゆくという方法がとられる。

初期の３次元構造を構成するには、得られた２次元画像に映し出された単粒子が３次元空間内でもつオイラー角を、単粒子の３次元構造が未知の状態で推定する必要がある。２次元画像に映し出された単粒子が３次元空間内でもつオイラー角を求める方法としては、ラドンの中央断面定理によるコモンラインを同定しそれぞれの相対角を求める方法が一般的である(非特許文献２を参照)。

コモンラインの同定方法には、平均画像の回転投射データを用いて各画像間の相関係数から求めるシノグラム法や、フーリエ空間上で複数の画像のコモンラインを求める方法などが知られている。しかしながら、これらの方法でコモンラインを同定すると、仮に極めてノイズの少ない画像情報に基づいて同定したとしても、その結果に大きな誤差を生じる可能性がある。これは、コモンラインを同定する演算過程において、２次元の画像情報が１次元の投射データに変換されるために、画像情報の一部が欠落してしまうことによるものと考えられる。すなわち、コモンラインの同定を行なうことで３次元構造を再構築する方法は、２次元画像を１次元画像に潰すというプロセスを必要とし、このために大きな誤差を生じ易いという欠点がある。

ところで、初期の３次元構造が得られると、これをリファレンス（初期構造モデル）として、２次元画像に映し出された単粒子が３次元空間内でもつオイラー角の最適化が行われる。この最適化の方法としては、リファレンスの再投射画像とのマッチングにより角度を推定する方法や、オイラー角に対応する３次元ラドン変換を用いて相互相関により角度を求める方法（たとえば、非特許文献２を参照)、あるいはwavelet変換により求める方法などが知られている。

一般に、リファレンスデータを用いたマッチングは、コモンライン法に比べてロバストでノイズに対して強く精度も高い。これは、２次元の画像情報をそのまま処理していることによると考えられる。しかしながら、初期の３次元構造（すなわちリファレンス）の精度が低い場合には、角度の微調整を繰り返す毎に３次元的な形状が崩れる場合も生じ得る。このため、良好な初期の３次元構造を如何にして作製するのかが極めて重要となる。

また、従来の３次元構造の再構成手法は、基礎となる２次元画像にノイズが含まれていないこと（換言すれば、単粒子の３次元的な中心が合っていること）を前提とし、この前提の下での３次元構造の再構築を行なっている。このため、無視できないノイズが含まれている２次元画像に基づいて３次元構造の再構築を行なうと、得られた結果の精度は極めて低いものとならざるを得ないという問題がある。

本発明は、このような問題に鑑みてなされたもので、その目的とするところは、角度情報がなくても２次元画像そのままのイメージからの３次元構造の再構築を可能とし、立体構造が未知の対象物の３次元構造を高精度で再構築することができる画像処理方法およびこれを実行するためのプログラムを提供することにある。

本発明はこのような課題を解決するために、請求項１に記載の発明は、複数の２次元画像から３次元構造を再構築するための画像処理方法であって、Ｎ枚（Ｎは２以上の整数）の２次元画像を仮想球面上の初期画像位置にランダムに配置する第１のステップと、該ランダム配置されたＮ枚の２次元画像を前記仮想球面の中心点に投射して初期３次元構造を構築する第２のステップと、該初期３次元構造を前記初期画像位置に再投射してＮ枚の新たな２次元画像を得る第３のステップと、前記初期画像位置の各点において再投射前後の２次元画像の相互相関を求めて第１のスコア（Ｌ_１）を算出する第４のステップと、前記再投射後のＮ枚の２次元画像から１枚の２次元画像を任意に選択し、該選択された２次元画像の座標パラメータをランダムに変化させる第５のステップと、該座標パラメータ変化後の１枚の２次元画像とその余の（Ｎ−１）枚の２次元画像から３次元構造を再構築する第６のステップと、該再構築された３次元構造を前記仮想球面上に再投射してＮ枚の新たな２次元画像を得る第７のステップと、前記再投射前後の２次元画像の相互相関を求めて第２のスコア（Ｌ_２）を算出する第８のステップと、該第２のスコア（Ｌ_２）と前記第１のスコア（Ｌ_１）との比較に基づいて前記第５のステップにおける座標パラメータの変化の許容もしくは拒絶を判断する第９のステップと、前記座標パラメータの変化が許容された場合に前記選択された２次元画像の座標パラメータを置き換える第１０のステップと、を備えていることを特徴とする。

請求項２に記載の発明は、請求項１に記載の画像処理方法において、前記第５ないし第１０のステップを所定の回数だけ繰り返す第１１のステップを備えていることを特徴とする。

請求項３に記載の発明は、請求項１または２に記載の画像処理方法において、前記第９のステップにおける判断は、前記第１のスコア（Ｌ_１）と第２のスコア（Ｌ_２）がＬ_２＞Ｌ_１の関係を満足する場合に前記座標パラメータの変化を許容するものであることを特徴とする。

請求項４に記載の発明は、請求項３に記載の画像処理方法において、前記第９のステップにおける判断は、Ｌ_２≦Ｌ_１の場合には、温度パラメータＴと遷移確率Ｐ（ΔＬ＝Ｌ_１−Ｌ_２）に基づいて前記座標パラメータの変化を許容するサブステップを備えていることを特徴とする。

請求項５に記載の発明は、請求項２に記載の画像処理方法において、前記第１１のステップに続き、温度パラメータＴ_１を所定値だけ減少させて新たな温度パラメータＴ_２を設定し、該温度パラメータＴ_２を用いて前記第５乃至第１１のステップを所定回数だけ反復する第１２のステップを備えていることを特徴とする。

請求項６に記載の発明は、請求項１乃至５の何れか１項に記載の画像処理方法において、前記座標パラメータは３つの角度変数と２つの位置変数を含んでいることを特徴とする。

請求項７に記載の発明は、請求項６に記載の画像処理方法において、前記３つの角度変数はオイラー角（α、β、γ）であり、前記２つの位置変数は前記２次元画像の中心位置調整用の変数（ｐ_ｘ、ｐ_ｙ）であることを特徴とする。

請求項８に記載の発明は、請求項１乃至７の何れか１項に記載の画像処理方法において、前記スコアＬは、次式で算出されるものであることを特徴とする。

ここで、ｍとｎは２次元画像中のピクセル位置を示す変数であり、

は２次元画像ｆ_ｉの平均輝度値であり、

は再投射画像ｇ_ｉの平均輝度値である。

請求項９に記載の発明は、３次元角度の推定方法であって、請求項６乃至８の何れか１項の画像処理方法により、前記２次元画像に映し出された３次元構造での角度情報を求めることを特徴とする。

請求項１０に記載の発明は、プログラムであって、コンピュータに、請求項１乃至８の何れか１項に記載された各ステップを実行させるためのものであることを特徴とする。

本発明の画像処理手法によれば、２次元の平均画像から３次元の初期構造モデルを作製し、平均画像と再投射画像との直接的な相関値を最適化の尺度として用いて３次元構造を最適化することとしたので、角度情報のない２次元画像からでも３次元構造を再構築することが可能な画像処理手法が提供される。これにより、立体構造が未知の対象物の３次元構造を高精度で再構築することができる画像処理方法が提供される。

また、本発明によれば、３次元角度が未知の２次元画像から、当該２次元画像に映し出された３次元構造での角度を簡便に推定することが可能となる。

さらに、本発明の画像処理プログラムによれば、これらの画像処理方法をコンピュータで自動的に実行させることができる。

以下に、図面を参照して本発明を実施するための最良の形態について説明する。

本発明の３次元構造の再構築方法は、先ず、２次元の平均画像から３次元の初期構造モデルを作製し、この３次元構造を最適化するものである。本発明の３次元構造の再構築方法では、平均画像の角度推定にシミュレーテッド・アニーリング（Simulated Annealing:ＳＡ）法が、３次元構造の最適化の尺度として平均画像と再投射画像との直接的な相関値が用いられ、これにより、ノイズに強くロバストな角度推定が可能となる。つまり、３次元の角度推定にシミュレーテッド・アニーリング法を適用することで、あたかもコンピュータ上で結晶化が行われるように、平均画像から３次元構造を自然に構築することが可能となる。

より具体的に説明すると、平均画像からフィルタードバックプロジェクション（ＦＢＰ）により３次元構造をつくりながら、逐次的に平均画像の３次元角度（オイラー角）と画像の位置パラメータを変化させ、その都度、平均画像と再プロジェクション画像との相関をスコアとして計算する。そして、スコア値が徐々に向上するように、シミュレーテッド・アニーリング・アルゴリズムを用いて繰り返し計算を行う。この繰り返し計算により３次元構造は最適化され、真の立体像に近い３次元構造が再構築されることになる。

先ず、本発明の理解の容易化のために、シミュレーテッド・アニーリング法について簡単な説明をする。シミュレーテッド・アニーリング法は、最適化問題を解くためのモンテカルロ法的な近似解法の一つである。最適化問題とは、特定の制約条件下で与えられた状態空間で定義される関数の最大値または最小値を与える状態空間の要素を求める問題であり、シミュレーテッド・アニーリング法は、たとえば、高温状態にある金属を徐冷することで低欠陥の金属結晶を得るというアニーリング（焼き鈍し）に着想を得て、このような物理現象をコンピュータ上で模擬的に実現して最適化問題を解こうとする手法であり、極めて強力な最適化手法の一つである。

シミュレーテッド・アニーリング法においては、対象となる最適化問題での変数｛ｘ_ｉ｝と最適化の程度を判断するためのコスト関数Ｅがそれぞれ、物理系における粒子とその系のエネルギに例えられる。一般に、系のエネルギに相当するコスト関数Ｅがｎ個の変数（すなわち、ｉ＝１〜ｎ）の関数であるとすると、コスト関数Ｅは、Ｅ＝Ｅ（ｘ_１，ｘ_２，・・・，ｘ_ｎ）と表記される。３次元構造の再構築のような画像復元の問題では、各ピクセル値を変数｛ｘ_ｉ｝とし、コスト関数Ｅを最小にする変数｛ｘ_ｉ｝の組を見つけることで最適化問題が解かれる。

現実の物理系においては、最小エネルギ状態とは結晶構造を構成する各原子の熱揺らぎがない状態に対応するが、系の温度をＴ（≠０Ｋ）とすると、各原子は熱的に揺らいでおり、それぞれの自由エネルギが最も低い状態で熱平衡となる。系が熱平衡状態にある場合に基底状態よりΔＥだけ高いエネルギ状態をとる確率は、ボルツマン分布Ｐに従って、Ｐ（ΔＥ）＝ｅｘｐ（−ΔＥ／ｋＴ）である。ここで、ｋはボルツマン定数であり、系の温度Ｔを徐々に下げてゆくと原子の熱揺らぎは抑制され、系のエネルギが低下することになる。

シミュレーテッド・アニーリング法による最適化のプロセスでは、先ず、初期変数｛ｘ_ｉ｝として適当な値を与えておき、各パラメータにランダムな微小変化Δｘ_ｉを順番に加えてコスト関数の変化量ΔＥを求める（アニーリング）。この変化量ΔＥが負であれば、微小変化Δｘ_ｉを最適解に近づけるものとして採用する一方、変化量ΔＥが正であれば、微小変化Δｘ_ｉをボルツマンの確率に従って受け入れる。このような演算を全ての変数ｘ_ｉに対して繰り返し実行することで、コンピュータの中で擬似的な熱平衡状態を作り出すのである。

次に、系の温度Ｔを僅かに下げる（クーリング）と、系は非平衡状態となり、上記と同様の演算を行なうことで新たな熱平衡状態が実現される。このような演算を系のエネルギが充分に低くなるまで実行すると、最終的には最適解を得ることができる。

図１は、本発明の３次元構造の再構築方法の手順を説明するためのフローチャートであり、図２は、本発明の３次元構造の再構築方法の各ステップで得られる３次元構造の様子を説明するための図である。また、図３は、本発明の３次元構造の再構築方法を実行するに際して選択される５つの変数を説明するための図である。なお、以下の説明では、３次元構造が再構築される対象物の例としてタンパク質の単粒子を採用する。また、以下に説明する本発明の３次元構造の再構築方法の各手順をプログラム化すれば、必要な画像処理をコンピュータで実行可能なことは明らかである。

本発明の３次元構造の再構築方法においては、図３に示すように、立体形状を有する試料から得た複数の平均画像を仮想的な球体の表面上に配置させ、これらの平均画像を基にして、球体の中心部に３次元構造を再構築する。ここで、それぞれの平均画像は本来、実際の試料を特定の方向から撮影した２次元画像である。したがって、それぞれの平均画像の撮影方向を特定するための３つの角度変数と、被写体である試料の中心位置を２次元平面に投射させた位置と平均画像の中心位置とのズレを２次元平面内で修正するための２つの位置変数の、計５つの変数を選択する。そして、各平均画像ごとにこれらの変数を最適化すれば、本来の試料形状である３次元構造が再構築されるはずである。

ここでは、図３に示すように、地球儀様の仮想的な透明球体の表面に平均画像をランダムに配置させてこの球体の中心部に３次元構造を再構築するに際し、座標パラメータとして、以下の変数を用いている。すなわち、３つの角度変数であるオイラー角として、回転軸であるＸ軸、Ｙ軸、およびＺ軸それぞれに対応させてα、β、γとした。なお、回転させる軸の順番は、Ｚ軸、Ｙ軸、Ｘ軸の順とした。これにより、地球儀の経度がγに、緯度がβに、そしてαが平均画像中の回転角に対応する。なお、以下では、βとγを「オイラー角」、αを「回転角」と呼ぶこともある。この角度変数に加えて、平均画像の縦と横の位置変数（ｐ_ｘ、ｐ_ｙ）も考慮することとして、中心位置のズレ調整を可能とした。このように、本発明の３次元構造の再構築方法においては、各平均画像の変数（移動パラメータ）の数は５となる。

先ず、変数の初期化を行い（ステップＳ１０２）、３次元構造を再構築するための複数の２次元画像（平均画像）を球体表面上にランダムに配置する。具体的には、Ｎ枚の平均画像のそれぞれのオイラー角をランダムに設定し、かつその回転角と位置変数（ｐ_ｘ、ｐ_ｙ）を何れも０とする（ステップＳ１０３、図２（ａ））。

次に、初期設定された角度変数でフィルタードバックプロジェクション（ＦＢＰ）を行い、初期３次元構造を求め（ステップＳ１０４、図２（ｂ））、この初期３次元構造から各平均画像位置に再プロジェクションを行う（ステップＳ１０５）。そして、この再投射画像と元の平均画像との相互相関を計算して各平均画像のスコアの総和（初期スコア）を求める（ステップＳ１０６）。なお、後の取り扱いを簡便化するために、各平均画像のスコアの総和を全平均画像枚数で割って正規化した値を「初期スコア」とし、この値に１００を掛けることで初期スコアの最大値が１００になるようにしている。

このような初期設定に続き、サイクル数Ｃを０に設定して（ステップＳ１０７）、Ｎ枚の平均画像の中から変数変化させる平均画像をランダムに１枚選択し、この平均画像のオイラー角と回転角・位置変数をランダムに変化させる。

平均画像のオイラー角の変化は、平均画像の全域移動と変数変化前の平均画像の位置を中心とする制限移動が、それぞれ１／２の確率となるように行なう。対称性をもたない粒子（たとえば、タンパク質ＲＦ２）の全域移動では、角度変数γが−１８０°〜＋１８０°、角度変数βが−９０°〜＋９０°、の範囲内でランダムに位置を決定する。

一方、対称性をもつ粒子の場合には、その対称性に応じて角度変数が取り得る範囲に制限を加えることが可能である。例えば、４回対称性をもつ粒子の場合は、角度変数γのみを０〜＋９０°の範囲とする。ここで、角度変数βが＋９０°あるいは−９０°の方向に近づくにつれて、角度変数γにおける１°刻みの球体上の間隔は狭くなる。これは、ちょうどメルカトル世界地図で緯度が高くなるほど面積が拡大されることに対応する。もし、一様乱数で平均画像の位置を決定した場合は、上部と下部付近（すなわち、角度変数β＝±９０°）に移動が片寄ることになる。そこで、仮想球体をＺ軸方向から眺めたＸＹ平面座標を用いて、この範囲でランダムに移動点を決める。

制限移動では、現在の位置を中心として正規乱数σ＝１０°で設定する。平均画像の回転角度（オイラー角のα）についても、０°〜＋３６０°（＝±０°）の範囲での一様分布乱数による設定と、現在の角度変数を中心とした正規乱数σ＝１０°の設定とが、それぞれ１／２の確率となるように角度変数を変化させる。また、平均画像の位置変数の変化は、σ＝１pixelの正規分布乱数により決定した。

初期設定後は、以下に説明するステップを所定のサイクル数Ｃだけ繰り返すこととなるが、平均画像の変数変化処理が受け入れられなかったサイクルの数（失敗回数）をIte_１、平均画像の変数変化処理が受け入れられたサイクルの数（成功回数）をIte_２、としてカウントする。したがって、サイクルの反復数の初期値を、予めIte_１＝０、Ite_２＝０と設定する（ステップＳ１０８）。

そして、Ｎ枚の平均画像の中から平均画像をランダムに１枚選択して、この平均画像のオイラー角と回転角・位置変数をランダムに変化させ（ステップＳ１０９）、この平均画像の移動後の位置で当該平均画像のＦＢＰを行い、再び３次元構造を構築する（ステップＳ１１０、図２（ｃ））。

ここで、ｉ番目の平均粒子画像をｆ_ｉ(px_i,py_i)とし、その高周波領域の増強フィルタの画像をｆ_ｉ(px_i,py_i)’、オイラー角α、β、γのバックプロジェクション演算子をＨ(α_i,β_i,γ_i)とすると、３次元構造の密度データＶは、次式（１）で与えられる。ここで、Ｎは計算に使用した平均画像の枚数であり、px_iとpy_iは平均画像を縦方向および横方向（ｘ方向およびｙ方向）にずらしたピクセル（pixel）の数である。

このようにして再構築された３次元構造Ｖを基に、全ての平均画像の対応する角度に対して再プロジェクションを行う（ステップＳ１１１、図２（ｄ））。この再プロジェクション画像をｇ_ｉとすると、再プロジェクション演算はバックプロジェクション演算の逆演算であるから、再プロジェクション演算子はバックプロジェクション演算子Ｈ(α_i,β_i,γ_i)を逆変換関数として次式（２）で与えられる。

次に、この再投射画像ｇ_ｉと平均画像ｆ_ｉとの相互相関の２乗値を計算し、この総和を新しいスコアＬ_ｎｅｗとする（ステップＳ１１２）。ここでのＬ_ｎｅｗは請求項１のＬ_２に対応するものである。具体的なスコアＬの計算は、次式（３）および（４）で実行される。

ここで、

は平均画像の平均輝度値であり、

は再プロジェクション画像の平均輝度値である。また、ｍとｎはピクセル（pixel）の縦と横の位置を示す変数（ｍ_ｃおよびｎ_ｃのそれぞれを画像の縦と横の中心を示す定数）としている。なお、ここでの平均画像ｆ_ｉは、平均画像画位置を修正した後のものｆ_ｉ(px_i,py_i)を用いることとしている。

このような各平均画像の移動（角度変数および位置変数の変化）が選択されるか拒絶されるかは、シミュレーテッド・アニーリングのアルゴリズムにより決定される。もし、新しいスコアＬ_ｎｅｗが以前のスコアＬ_ｏｌｄよりも高い場合（Ｌ_ｎｅｗ＞Ｌ_ｏｌｄ、ステップＳ１１３：Ｙｅｓ）は、この角度・位置移動を確定して新しいスコアＬ_ｎｅｗを以前のスコアＬ_ｏｌｄと置き換え（ステップＳ１１４）、成功回数Ite_２の値を１付加してステップＳ１１６へと移行する。ここでのＬ_ｏｌｄは請求項１のＬ_１に対応するものである。

一方、新しいスコアＬ_ｎｅｗが以前のスコアＬ_ｏｌｄ以下の場合（Ｌ_ｎｅｗ≦Ｌ_ｏｌｄ、ステップＳ１１３：Ｎｏ）は、以前のスコアＬ_ｏｌｄと新しいスコアＬ_ｎｅｗの差（ΔＬ＝Ｌ_ｏｌｄ−Ｌ_ｎｅｗ）と現在の温度Ｔとから、ボルツマンの確率分布（式５）にしたがって遷移確率Ｐ（ΔＬ）を求め（ステップＳ１１７）、０〜１の範囲にある乱数値と比較して、遷移確率Ｐ（ΔＬ）が乱数値よりも大きければ（ステップＳ１１８：Ｙｅｓ）、移動が確定されてステップＳ１１４に移行する。一方、遷移確率Ｐ（ΔＬ）が乱数値以下であれば（ステップＳ１１８：Ｎｏ）、以前の角度・位置へと戻されて（ステップＳ１１９）、失敗回数Ite_１の値を１付加してステップＳ１１６へと移行する。

遷移確率Ｐ（ΔＬ）は温度パラメータＴにより制御されており、比較的温度が高い場合は多少スコアＬが低下しても移動が確定される確率は高くなる。一方、比較的温度が低い場合は、スコアＬが高くなる移動のみが確定されるようになる。なお、計算開始の初期段階では、温度パラメータＴは比較的高い値をとるため、スコアＬには依存せずにほぼランダムに平均画像が移動し、温度パラメータＴが低下する（クーリング）につれて徐々にスコアＬが高くなる方向にのみ平均画像が移動するようになる。

ステップＳ１１６では、失敗回数Ite_１および成功回数Ite_２が、予め与えられた最大失敗回数Ite_ｍａｘ１および最大成功回数Ite_ｍａｘ２と比較され、失敗回数Ite_１または成功回数Ite_２の少なくとも一方が最大回数以上であれば（ステップＳ１１６：Ｙｅｓ）、ステップＳ１２１に移行して温度パラメータＴを低く設定し直す（クーリング）。一方、失敗回数Ite_１と成功回数Ite_２の何れもが最大回数に満たなければ（ステップＳ１１６：Ｎｏ）、ステップＳ１０９に戻って以降の各ステップを繰り返す。

このような一連の処理を繰り返し、系の温度を徐々に下げてゆく。ここでの繰り返しプロセスは、たとえば、平均画像枚数Ｎの１０倍の画像移動が成功した場合（Ite_２ｍａｘ＝１０Ｎ）、あるいは平均画像枚数Ｎの３０倍の画像移動を試行した場合（Ite_１ｍａｘ＝３０Ｎ）、を１つのサイクルとし、１つのサイクルが終了するとステップＳ１２２でサイクル数Ｃに１が付加される。

温度パラメータＴは、次式（６）および（７）にしたがって指数関数的に低下させる。ここで、Ｔ_０は初期温度、Ｔ_ｅは終了温度である。また、Ｃ_ｍａｘは繰り返し演算の最大サイクル数であり、サイクル数ＣはステップＳ１２３で最大サイクル数Ｃ_ｍａｘと比較され、サイクル数Ｃが最大サイクル数Ｃ_ｍａｘよりも大きいと判断（ステップＳ１２３：Ｙｅｓ）された場合には演算が終了し（ステップＳ１２４）、サイクル数Ｃが最大サイクル数Ｃ_ｍａｘ以下であればステップＳ１０８に移行してその後の各ステップを繰り返す。そして、最終的に図２（ｅ）に示すような、最適化された３次元構造が再構築される。

たとえば、最大サイクル数Ｃ_ｍａｘは１００〜１０００程度に設定され、初期温度はＴ_０＝１、終了温度はＴ_ｅ＝０．０１などとされるが、これらのパラメータは処理するデータに応じて適宜変更可能であることはいうまでもない。

対称性を有するタンパク質粒子の画像を処理対象とする場合は、その対称性に応じてオイラー角の範囲を狭めることができる。ここで使用したＲＦ２は、完全に非対称な粒子のため、オイラー角の範囲は全域となる（γ＝−１８０°〜＋１８０°、β＝−９０°〜＋９０°）。

以下に、実施例により本発明をより詳細に説明する。

本実施例は、本発明の手法によりＲＦ２タンパクのモデルデータからプロジェクション画像を作製した例である。

ＲＦ２は完全な非対称性を有するタンパク質であり、その構造は極めて複雑である。テストデータとして、ＲＦ２の単粒子の３次元モデルデータからオイラー角や回転角をランダムに変化させた６０枚のプロジェクション画像を作製した。

図４（ａ）は、テストデータ用に作製した６０枚のプロジェクション画像の中から、１０枚の画像を例として示す図である。また、図４（ｂ）は、これらの１０枚のプロジェクション画像の位置を±２pixelの範囲でランダムに変化させ、さらに正規分布1０σのノイズを加えた画像である。本実施例では、図４（ｂ）に示したような、テストデータ用に作製した６０枚のプロジェクション画像に位置変化とノイズ付与を施した６０枚のプロジェクション画像をテストデータとして画像処理を実行している。

先ず、初期値として、各ＲＦ２の６０枚のテストデータをオイラー角上にランダムに配置する。この様子を図４（ｃ）に示してある。この画像配置からランダムに選択された画像の位置を変化させ、シミュレーテッド・アニーリングのプロセスによりその画像位置変化の許容もしくは拒絶の決定を行う。

画像位置変化に伴う一連のサイクルを繰り返しながら、徐々に温度を低下させてゆくと、所定のサイクル数の計算を繰り返した後には、図４（ｄ）に示すように、各平均画像のオイラー角、回転角、位置がほぼ妥当な配置をとることが確認できる。

図５は、一連のサイクルを繰り返して再構築される３次元構造（図５（ａ））、スコア変化の様子（図５（ｂ））、および再構築された３次元構造とモデルデータとの比較結果（図５（ｃ））、を説明するための図である。

計算初期段階では、系の温度パラメータＴが高いために、各平均画像はスコアに依存せず移動し、得られるスコアは低い（図５（ｂ））。また、再構築される３次元構造もいびつな形状を有しておりノイズも極めて高い(図５（ａ）のサイクル０〜２８０)。

温度パラメータが０．３以下にまで低下すると、スコアは徐々に高くなり、再構築される３次元構造もその輪郭がはっきりしてくると同時にその周囲のノイズも低下する(図５（ａ）のサイクル２８０〜３８０)。

さらに、温度パラメータが０．０５以下まで低下すると、スコアの上昇は徐々に鈍り、３次元構造の細部の構造が現れてくる(図５（ａ）のサイクル４８０)。この段階では、最終的な位置や角度の微調整が行われ、最終的に再構築された３次元構造は、図５（ｃ）に示すように、極めてモデルデータと類似している。つまり、本発明の画像処理手法により、高い精度で３次元構造が再構築されることがわかる。

本発明は、単粒子構造解析における３次元角度推定問題に、シミュレーテッド・アニーリングのプロセスを適応したものとして把握することも可能である。すなわち、従来の画像処理方法では極めて困難であった非対称粒子や対称性の低い粒子のオイラー角（３次元角度）を推定し、この３次元角度に基づいて自動的に３次元構造が再構築される。

また、本発明の画像処理方法は、２次元画像を用いたコンピュータ上での「結晶化プロセス」として捉えることができ、繰り返し計算の回数を増やして、温度の低下をスローにすることで、得られる３次元構造の再構築精度を向上させることができる。したがって、今後のコンピュータの計算速度の向上や並列化などによって、ほぼ完全な３次元構造を自動的に再構築することが可能である。

さらに、本発明の画像処理方法につき、タンパク質の単粒子を電子顕微鏡で撮影した画像の処理を例として説明してきたが、本発明はこのような画像処理に限定されるものではなく、透過型電子顕微鏡（ＴＥＭ）や走査型電子顕微鏡（ＳＥＭ）で撮影された種々の試料の２次元画像は勿論のこと、原子間力顕微鏡（ＡＦＭ）や光学顕微鏡（ＯＭ）により撮影された画像やＸ線トポグラフ（ＸＲＴ）像、あるいは航空写真など、極めて広範囲に適用することができる。

以上、本発明の３次元構造の再構築に関する画像処理の手法について説明したが、上記実施例は本発明を実施するための例にすぎず、本発明はこれらに限定されるものではない。これらの実施例を種々変形することは本発明の範囲内にあり、更に本発明の範囲内において他の様々な実施例が可能であることは上記記載から自明である。

本発明は、角度情報がなくても２次元画像そのままのイメージからの３次元構造の再構築を可能とし、立体構造が未知の対象物の３次元構造を高精度で再構築することができる画像処理方法およびこれを実行するためのプログラムを提供する。

本発明の３次元構造の再構築方法の手順を説明するためのフローチャートである。 本発明の３次元構造の再構築方法の各ステップで得られる３次元構造の様子を説明するための図である。 本発明の３次元構造の再構築方法を実行するに際して選択される５つの変数を説明するための図である。 図（ａ）はテストデータ用に作製した６０枚のプロジェクション画像の中から、１０枚の画像を例として示す図、図（ｂ）はこれらの１０枚のプロジェクション画像の位置を±２pixelの範囲でランダムに変化させ、さらに正規分布1０σのノイズを加えた画像を示す図、図（ｃ）は初期値として各ＲＦ２の６０枚のテストデータをオイラー角上にランダムに配置した様子を示す図、そして図（ｄ）は最終的に得られる３次元構造をモデルデータと対比して示す図である。 図（ａ）は一連のサイクルを繰り返して再構築される３次元構造を示す図、図（ｂ）はスコア変化の様子を示す図、そして図（ｃ）は再構築された３次元構造とモデルデータとの比較結果を説明するための図である。

Claims (10)

1. 複数の２次元画像から３次元構造を再構築するための画像処理方法であって、
   Ｎ枚（Ｎは２以上の整数）の２次元画像を仮想球面上の初期画像位置にランダムに配置する第１のステップと、
   該ランダム配置されたＮ枚の２次元画像を前記仮想球面の中心点に投射して初期３次元構造を構築する第２のステップと、
   該初期３次元構造を前記初期画像位置に再投射してＮ枚の新たな２次元画像を得る第３のステップと、
   前記初期画像位置の各点において再投射前後の２次元画像の相互相関を求めて第１のスコア（Ｌ_１）を算出する第４のステップと、
   前記再投射後のＮ枚の２次元画像から１枚の２次元画像を任意に選択し、該選択された２次元画像の座標パラメータをランダムに変化させる第５のステップと、
   該座標パラメータ変化後の１枚の２次元画像とその余の（Ｎ−１）枚の２次元画像から３次元構造を再構築する第６のステップと、
   該再構築された３次元構造を前記仮想球面上に再投射してＮ枚の新たな２次元画像を得る第７のステップと、
   前記再投射前後の２次元画像の相互相関を求めて第２のスコア（Ｌ_２）を算出する第８のステップと、
   該第２のスコア（Ｌ_２）と前記第１のスコア（Ｌ_１）との比較に基づいて前記第５のステップにおける座標パラメータの変化の許容もしくは拒絶を判断する第９のステップと、
   前記座標パラメータの変化が許容された場合に前記選択された２次元画像の座標パラメータを置き換える第１０のステップと、
   を備えていることを特徴とする画像処理方法。
2. 前記第５ないし第１０のステップを所定の回数だけ繰り返す第１１のステップを備えていることを特徴とする請求項１に記載の画像処理方法。
3. 前記第９のステップにおける判断は、前記第１のスコア（Ｌ_１）と第２のスコア（Ｌ_２）がＬ_２＞Ｌ_１の関係を満足する場合に前記座標パラメータの変化を許容するものであることを特徴とする請求項１または２に記載の画像処理方法。
4. 前記第９のステップにおける判断は、Ｌ_２≦Ｌ_１の場合には、温度パラメータＴと遷移確率Ｐ（ΔＬ＝Ｌ_１−Ｌ_２）に基づいて前記座標パラメータの変化を許容するサブステップを備えていることを特徴とする請求項３に記載の画像処理方法。
   
5. 前記第１１のステップに続き、温度パラメータＴ_１を所定値だけ減少させて新たな温度パラメータＴ_２を設定し、該温度パラメータＴ_２を用いて前記第５乃至第１１のステップを所定回数だけ反復する第１２のステップを備えていることを特徴とする請求項２に記載の画像処理方法。
6. 前記座標パラメータは３つの角度変数と２つの位置変数を含んでいることを特徴とする請求項１乃至５の何れか１項に記載の画像処理方法。
7. 前記３つの角度変数はオイラー角（α、β、γ）であり、前記２つの位置変数は前記２次元画像の中心位置調整用の変数（ｐ_ｘ、ｐ_ｙ）であることを特徴とする請求項６に記載の画像処理方法。
8. 前記スコアＬは、次式で算出されるものであることを特徴とする請求項１乃至７の何れか１項に記載の画像処理方法。
   ここで、ｍとｎは２次元画像中のピクセル位置を示す変数であり、
   は２次元画像ｆ_ｉの平均輝度値であり、
   は再投射画像ｇ_ｉの平均輝度値である。
9. 請求項６乃至８の何れか１項の画像処理方法により、前記２次元画像に映し出された３次元構造での角度情報を求めることを特徴とする３次元角度推定方法。
10. コンピュータに、請求項１乃至８の何れか１項に記載された各ステップを実行させるためのプログラム。