特許 7490185 解析対象の３次元再構成装置、３次元再構成方法及び３次元再構成プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-05-17

(45)【発行日】2024-05-27

(54)【発明の名称】解析対象の３次元再構成装置、３次元再構成方法及び３次元再構成プログラム

(51)【国際特許分類】

   H01J 37/26 20060101AFI20240520BHJP

   H01J 37/22 20060101ALI20240520BHJP

   G06T 7/55 20170101ALI20240520BHJP

【ＦＩ】

H01J37/26

H01J37/22 501Z

G06T7/55

【請求項の数】 16

(21)【出願番号】P 2021552414

(86)(22)【出願日】2020-10-14

(86)【国際出願番号】 JP2020038798

(87)【国際公開番号】W WO2021075465

(87)【国際公開日】2021-04-22

【審査請求日】2023-09-29

(31)【優先権主張番号】P 2019190733

(32)【優先日】2019-10-18

(33)【優先権主張国・地域又は機関】JP

(73)【特許権者】

【識別番号】510222006

【氏名又は名称】株式会社バイオネット研究所

(73)【特許権者】

【識別番号】515130201

【氏名又は名称】株式会社Ｐｒｅｆｅｒｒｅｄ Ｎｅｔｗｏｒｋｓ

(74)【代理人】

【識別番号】100122574

【弁理士】

【氏名又は名称】吉永 貴大

(72)【発明者】

【氏名】大橋 正隆

(72)【発明者】

【氏名】前田 新一

【審査官】右▲高▼ 孝幸

(56)【参考文献】

【文献】Sjors H. W. Scheres，A Bayesian View on Cryo-EM Structure Determination，Journal of Molecular Biology，2012年01月13日，volume 415, issue 2, pages 406-418

【文献】Michael Stolken et al.，Maximum likelihood based classification of electron tomographic data，Journal of Structural Biology，2011年01月，volume 173, issue 1, pages 77-85

【文献】Pilar Cossio et al.，BioEM: GPU-accelerated computing of Bayesian inference of electron microscopy images，Computer Physics Communications，2017年01月，volume 210, pages 163-171

【文献】Ali Punjani et al.，Building Proteins in a Day: Efficient 3DMolecular Structure Estimation withElectron Cryomicroscopy，IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE，2017年04月，vol.39, no.4, pp.706-718

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｈ０１Ｊ ３７／００

Ｇ０６Ｔ ７／５５

(57)【特許請求の範囲】

【請求項1】

電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成装置であって、
前記解析対象の３次元構造を投影することにより得られる投影像に関し、
前記解析対象の３次元構造の事後確率の下界を前記投影条件の確率分布及び前記解析対象の３次元構造に関して最大化することによって、前記解析対象の投影条件と３次元構造とを推定する、３次元再構成装置。

【請求項2】

前記投影条件の確率分布に関する前記解析対象の３次元構造の事後確率の下界の最大化は、
前記投影条件の確率分布を表現するパラメータを勾配法によって更新することにより、前記解析対象の３次元構造の事後確率を最大化させる第１の事後確率最大化手段と、
連続量である前記投影条件を離散化した投影条件ごとにその事後確率を計算することにより、前記解析対象の３次元構造の事後確率を最大化する第２の事後確率最大化手段と、
大域的最適化を使用して前記解析対象の３次元構造の事後確率を最大化する第３の事後確率最大化手段と、
を推定の経過に応じて切り替えることによって行われる、請求項１に記載の３次元再構成装置。

【請求項3】

前記投影条件は、
前記解析対象の変形を定めるパラメータからなる第１の投影条件、前記電顕画像の撮影条件であるコントラスト伝達関数（CTF）に関する補正量からなる第２の投影条件、及び、粒子画像に目的とする前記解析対象が適切に写っているか否かを示す量からなる第３の投影条件の少なくとも１つを含む、請求項１又は２に記載の３次元再構成装置。

【請求項4】

前記投影条件は、
前記解析対象がどの方向から投影されているかを表す投影方向からなる第４の投影条件、前記解析対象の中心位置を表す平行移動量からなる第５の投影条件、前記解析対象の対称性を表現する量からなる第６の投影条件、及び、複数の３次元構造のいずれに属するかを示す量からなる第７の投影条件の少なくとも１つを含む、請求項１から３のいずれか１項に記載の３次元再構成装置。

【請求項5】

前記解析対象の３次元構造の事後確率の下界の最大化において、
前記解析対象の３次元構造を複数の基底関数の線形和によって表現する機能と、
実空間における２つの前記基底関数の距離を求める機能と、
２つの前記基底関数の距離が近傍であるか判定する機能と、
前記基底関数に対応する係数の差を計算する機能と
前記基底関数に対応する係数の差がエッジかどうか判定する機能を有し、
すべての前記基底関数の組み合わせに対して求めた前記距離と、前記係数の差と、前記エッジかどうかの判定結果に基づいた３次元構造の事前確率を使って、前記３次元構造の事後確率を最大化する、請求項１から４のいずれか１項に記載の３次元再構成装置。

【請求項6】

前記投影像の生成において、
前記解析対象の３次元構造を複数の基底関数の線形和によって表現する機能と、
補助画素を作成して撮像素子により決まる解像度以上の投影像を生成する機能と、を有する、請求項１から４のいずれか１項に記載の３次元再構成装置。

【請求項7】

前記第３の投影条件の粒子画像に目的とする解析対象が適切に写っているか否かを示す量の事後確率の推定において、
粒子画像の特徴量を抽出する機能と、
粒子画像に目的の解析対象が適切に写っているか否かをラベル付けする機能と、
粒子画像に目的の解析対象が適切に写っているか否かを識別する識別器に前記特徴量と前記ラベルによって学習させる機能と、
前記識別器から出力される前記粒子画像に目的の解析対象が写っている確率を使って前記解析対象の３次元構造の事後確率の下界を前記第３の投影条件の確率分布と前記解析対象の３次元構造に関して最大化する機能と、を有する、請求項３に記載の３次元再構成装置。

【請求項8】

電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成装置であって、
Cryo-TEM法により生成された解析対象が氷包埋された試料を用いて、前記電顕画像からコントラスト伝達関数（CTF）を計算するCTF推定部と、
前記電顕画像から解析対象を抽出処理する粒子抽出部と、
前記解析対象の３次元構造を再構成する再構成部と、を備え、
前記再構成部は、
前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する下界最大化処理部と、
与えられた前記投影条件により、投影像と投影像生成行列を生成する投影像生成部と、
前記下界最大化処理部による前記解析対象の３次元構造の事後確率の下界の最大化について収束判定を行う収束判定部と、を備える、３次元再構成装置。

【請求項9】

電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成方法であって、
Cryo-TEM法により生成された解析対象が氷包埋された試料を用いて、前記電顕画像からコントラスト伝達関数（CTF）を計算するCTF推定ステップと、
前記電顕画像から解析対象を抽出処理する粒子抽出ステップと、
前記解析対象の３次元構造を再構成する再構成ステップと、を備え、
前記再構成ステップは、
前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する下界最大化処理ステップと、
与えられた前記投影条件により、投影像と投影像生成行列を生成する投影像生成ステップと、
前記下界最大化処理ステップにおける前記解析対象の３次元構造の事後確率の下界の最大化について収束判定を行う収束判定ステップと、を備える、３次元再構成方法。

【請求項10】

電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成プログラムであって、
コンピュータに、
Cryo-TEM法により生成された解析対象が氷包埋された試料を用いて、前記電顕画像からコントラスト伝達関数（CTF）を計算するCTF推定機能と、
前記電顕画像から解析対象を抽出処理する粒子抽出機能と、
前記解析対象の３次元構造を再構成する再構成機能と、を実行させ、
前記再構成機能は、
前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する下界最大化処理機能と、
与えられた前記投影条件により、投影像と投影像生成行列を生成する投影像生成機能と、
前記下界最大化処理機能における前記解析対象の３次元構造の事後確率の下界の最大化について収束判定を行う収束判定機能と、を備える、３次元再構成プログラム。

【請求項11】

前記解析対象が、分子である、請求項１から８のいずれか１項に記載の３次元再構成装置。

【請求項12】

前記分子が生体分子である、請求項１１に記載の３次元再構成装置。

【請求項13】

前記解析対象が、分子である、請求項９に記載の３次元再構成方法。

【請求項14】

前記分子が生体分子である、請求項１３に記載の３次元再構成装置。

【請求項15】

前記解析対象が、分子である、請求項１０に記載の３次元再構成プログラム。

【請求項16】

前記分子が生体分子である、請求項１５に記載の３次元再構成プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、解析対象の３次元再構成装置、３次元再構成方法及び３次元再構成プログラムに関する。

【背景技術】

【0002】

Cryo-TEM法は、試料となる生体分子を含む水溶液を液体窒素で冷やされた液化エタンで急速凍結することにより生体分子を非結晶の薄い氷に包埋し、氷包埋されたままの状態で低温ステージを持つ透過型電子顕微鏡を使用して観察する手法である。

【0003】

生体分子の３次元構造は、Cryo-TEM法により、氷包埋された多量の生体分子を透過型電子顕微鏡により撮影し、得られた画像（以降、電顕画像という場合がある）から生体分子が写っている部分を切り出した数万から数百万枚の粒子画像をもとに、推定される。２次元画像から３次元構造を推定する処理は、３次元再構成と呼ばれる。さらに、Cryo-TEM法により撮影された電顕画像をもとに生体分子の３次元構造を解析する手法は、単粒子解析と呼ばれる。

【0004】

３次元構造から２次元の投影像を生成する過程には、生体分子がどの方向から投影されているかを表す投影方向などの投影条件が必要とされる。しかしながら、生体分子は様々な向きで氷包埋されており、投影条件も粒子画像ごとに異なる。それに加えて、生体分子の向きは既知ではない。図２１は、生体分子の向きによって投影像が異なることを説明するためのイメージ図であり、投影方向（破線矢印）に対し、３つの氷包理された生体分子の向き及び対応する投影像を例示している。したがって、Cryo-TEM法による氷包埋された生体分子の３次元構造の推定には、投影方向等の未知の投影条件の推定も必要になる。そこで、従来、粒子画像ごとに、生体分子がどの方向から投影されているかを表す投影方向、生体分子の中心位置を表す平行移動量、生体分子の対称性を表現する量、及び３次元構造のいずれに属するかを示す量（クラス分類）についての投影条件の推定も行われてきた。

【0005】

ところで、透過型電子顕微鏡で生体分子のような非常に小さな物体を撮影する場合は、位相コントラストを使用して撮影を行う。位相コントラストとは、試料を通り抜けてきた電子の波の位相のずれによる干渉模様によって形成される像のことをいう。位相コントラストを撮影する場合、試料にフォーカスを合わせて撮影するとコントラストが非常に低くなるため、フォーカスを意図的にずらした状態で撮影する。フォーカスがずれた状態のことはデフォーカスと言い、フォーカスのずれ量のことをデフォーカス量という。さらに、位相コントラストによって得られる像には、収差とデフォーカスに起因する位相のズレが更に加わる。この位相のズレを＋１～－１の周波数応答として表したものをコントラスト伝達関数（以降、CTFという場合がある。）と呼ぶ。CTFは、図１２（左図：正面、右図：側面）に示すように、高周波ほど、＋１～－１間の振動が激しくなるという特性を持っている。透過型電子顕微鏡により撮影された生体分子の投影像は、位相のズレが加わることにより、図１３のようになる。図１３は、透過型電子顕微鏡により撮影された生体分子の投影像のシミュレーション画像（図１３の左右図は、後述するように、デフォーカス量の違いを示す。）である。

【0006】

生体分子の電子線による激しい損傷を防ぐため、照射できる電子線量は制限される。さらに、Cryo-TEM法では、凍結時に氷の厚さを薄く均一に作ることが困難である。これらの要因から、粒子画像の信号雑音比（以降、Ｓ／Ｎ比という場合がある。）は非常に低い。さらに、十分なコントラストを得るためにデフォーカス量を大きくすると、CTFが＋１～－１の間で激しく振動し、情報の欠落の度合いも大きくなってしまう。

【0007】

このような条件下で行われる３次元構造の推定は、推定精度を向上させるために多量の粒子画像を必要とする。多量の粒子画像を処理することは、実験にかかる労力や計算量の増大という問題を引き起こす。また、Ｓ／Ｎ比が低く、情報の欠落も多い中、高分解能を得ようとして３次元構造を推定すると、ノイズを含んだ粒子画像に３次元再構成像が過剰に適合してしまう過剰適合（overfitting）という問題が発生する。このような中、過剰適合を防ぎつつ生体分子の３次元構造を推定する方法の一つとして、非特許文献１に記載の手法がある。

【0008】

従来の単粒子解析では、撮影時のデフォーカス量を別の手法で計算し、計算したデフォーカス量が真値として３次元再構成を行う。前述のとおり、CTFは収差とデフォーカスに起因するものであり、デフォーカス量が変わるとCTFも変化する。しかしながら、前述のとおり、単粒子解析で使用される画像のＳ／Ｎ比が低いため、デフォーカス量を正しく計算することが難しく、従来は不確かなデフォーカス量を使って３次元再構成を行っていた。

【0009】

例えば、非特許文献１に記載の技術では、前述した４つの投影条件である、投影方向と、平行移動量と、生体分子の対称性と、クラス分類の投影条件（以降、隠れ変数という場合がある。）を一定間隔ごとに離散化し、離散化された条件ごとに事後確率を推定した後、推定した投影条件の事後確率を使って生体分子の３次元構造を推定するということが行われている。

【0010】

しかしながら、解析の対象となる生体分子は、非常に柔らかく変形しやすい。この点、非特許文献１の手法では、生体分子の変形が考慮されていない。また、他の手法では、らせん状の生体分子の場合、変形の影響を抑えるために、生体分子を細かいセグメントに分割すること一般的に行われるが、分割したセグメント内では変形が考慮されないため、本質的な部分が解決されているとは言えない。変形がある生体分子の粒子画像を使って３次元構造を推定すると、与えられた粒子画像を表現するように３次元構造を推定しようとするが、生体分子の変形が考慮されていないため、推定精度の低下につながる。

【0011】

生体分子の３次元再構成において使用される粒子画像は、急速凍結された試料が含まれた水溶液を透過型電子顕微鏡によって撮影することにより得られる。粒子画像にはCTFによる変調が加えられていることから、粒子画像に加えられているCTFの周波数応答を決めるデフォーカス量を別の手法で計算し、計算したデフォーカス量が真値であるとして３次元再構成を行う。しかしながら、透過型電子顕微鏡によって撮影された粒子画像は、Ｓ／Ｎ比が非常に低く、デフォーカス量を正しく計算することが難しいという問題がある。また、従来は粒子が多量に映った電顕画像ごとにデフォーカス量が計算されていたが、氷包埋された粒子は粒子ごとに氷の表面からの高さが異なるため、デフォーカス量も粒子ごとに異なるという問題もある。

【0012】

前述した図１３は、デフォーカス量の違いにおける投影像の変化を示したものである。左図がデフォーカス量１００Åの時の投影像で、右図がデフォーカス量５００Åの時の投影像である。デフォーカス量の違いにより投影像が大きく変化することがわかる。

【0013】

例えば、非特許文献１に記載の３次元再構成処理の流れは、デフォーカス量を計算するCTF推定、粒子画像の抽出を行う粒子抽出、粒子画像の２次元分類、粒子画像の３次元分類、精密化となっている。粒子画像の２次元分類及び３次元分類の後には、より高分解能な３次元再構成像を得るために、ユーザーによって、それ以降の工程で使う粒子画像を選別する処理が行われる。しかしながら、ユーザーによる粒子画像の選別には、ユーザーの恣意性が含まれ、最終的に推定される３次元再構成像が、ユーザーごとに異なってしまうという問題がある。

【0014】

投影条件の事後確率の計算は、投影方向や平行移動量や生体分子の対称性やクラス分類それぞれの投影条件を個別に離散化した後、全組み合わせに対して尤度を計算する必要があり、推定精度を向上させるために離散化する間隔を細かくした場合、組み合わせが膨大になるという問題がある。

【0015】

透過型電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像から行う生体分子の３次元再構成では、前述のとおり、粒子画像のＳ／Ｎ比が非常に低いことから、粒子画像に過剰に適合してしまう過剰適合が問題になる。非特許文献１に記載の技術では、３次元再構成像の事前確率（以降、事前分布又はPriorという場合がある。）を導入することで過剰適合を抑止している。事前確率による過剰適合の抑止は、再構成像が滑らかになる効果をもたらすが、抑止しすぎると分解能の低下につながるという問題がある。非特許文献１に記載の技術では、周波数空間上で３次元再構成を行っており、周波数空間上では実空間上での隣接関係の情報が失われてしまうため、隣接間の相関を利用した過剰適合の抑止を行うことは難しいという問題があった。

【0016】

非特許文献１に記載の手法では、周波数空間上で生体分子の３次元構造の推定を行っている。しかしながら、実際、センサーによる撮像は実空間上で行われており、撮像をより正確にモデル化する場合に、周波数空間では難しいという問題がある。

【先行技術文献】

【非特許文献】

【0017】

【文献】Scheres, Sjors HW. “A Bayesian view on cryo-EM structure determination.” Journal of molecular biology 415.2 (2012): 406-418.

【発明の概要】

【発明が解決しようとする課題】

【0018】

本発明は、上記従来技術の課題を解決するためになされたものであり、計算量を抑えつつ、高分解能な３次元再構成像を生成するための解析対象の３次元再構成装置、３次元再構成方法及び３次元再構成プログラムを提供することを目的とする。

【課題を解決するための手段】

【0019】

本発明は、上記目的を達成するために、以下の構成によって把握される。
（１）本発明の第１の観点は、電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成装置であって、前記解析対象の３次元構造を投影することにより得られる投影像に関し、前記解析対象の３次元構造の事後確率の下界を前記投影条件の確率分布及び前記解析対象の３次元構造に関して最大化することによって、前記解析対象の投影条件と３次元構造とを推定する。

【0020】

（２）前記投影条件の確率分布に関する前記解析対象の３次元構造の事後確率の下界の最大化は、前記投影条件の確率分布を表現するパラメータを勾配法によって更新することにより、前記解析対象の３次元構造の事後確率を最大化させる第１の事後確率最大化手段と、連続量である前記投影条件を離散化した投影条件ごとにその事後確率を計算することにより、前記解析対象の３次元構造の事後確率を最大化する第２の事後確率最大化手段と、大域的最適化を使用して前記解析対象の３次元構造の事後確率を最大化する第３の事後確率最大化手段と、を推定の経過に応じて切り替えることによって行われてもよい。

【0021】

（３）前記投影条件は、前記解析対象の変形を定めるパラメータからなる第１の投影条件、前記電顕画像の撮影条件であるコントラスト伝達関数（CTF）に関する補正量からなる第２の投影条件、及び、粒子画像に目的とする前記解析対象が適切に写っているか否かを示す量からなる第３の投影条件の少なくとも１つを含むようにしてもよい。

【0022】

（４）前記投影条件は、前記解析対象がどの方向から投影されているかを表す投影方向からなる第４の投影条件、前記解析対象の中心位置を表す平行移動量からなる第５の投影条件、前記解析対象の対称性を表現する量からなる第６の投影条件、及び、複数の３次元構造のいずれに属するかを示す量からなる第７の投影条件の少なくとも１つを含むようにしてもよい。

【0023】

（５）前記解析対象の３次元構造の事後確率の下界の最大化において、前記解析対象の３次元構造を複数の基底関数の線形和によって表現する機能と、実空間における２つの前記基底関数の距離を求める機能と、２つの前記基底関数の距離が近傍であるか判定する機能と、前記基底関数に対応する係数の差を計算する機能と前記基底関数に対応する係数の差がエッジかどうか判定する機能を有し、すべての前記基底関数の組み合わせに対して求めた前記距離と、前記係数の差と、前記エッジかどうかの判定結果に基づいた３次元構造の事前確率を使って、前記３次元構造の事後確率を最大化してもよい。

【0024】

（６）前記投影像の生成において、前記解析対象の３次元構造を複数の基底関数の線形和によって表現する機能と、補助画素を作成して撮像素子により決まる解像度以上の投影像を生成する機能と、を有してもよい。

【0025】

（７）前記第３の投影条件の粒子画像に目的とする解析対象が適切に写っているか否かを示す量の事後確率の推定において、粒子画像の特徴量を抽出する機能と、粒子画像に目的の解析対象が適切に写っているか否かをラベル付けする機能と、粒子画像に目的の解析対象が適切に写っているか否かを識別する識別器に前記特徴量と前記ラベルによって学習させる機能と、前記識別器から出力される前記粒子画像に目的の解析対象が写っている確率を使って前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する機能と、を有してもよい。

【0026】

（８）本発明の第２の観点は、電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成装置であって、Cryo-TEM法により生成された解析対象が氷包埋された試料を用いて、前記電顕画像からコントラスト伝達関数（CTF）を計算するCTF推定部と、前記電顕画像から解析対象を抽出処理する粒子抽出部と、前記解析対象の３次元構造を再構成する再構成部と、を備え、前記再構成部は、前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する下界最大化処理部と、与えられた前記投影条件により、投影像と投影像生成行列を生成する投影像生成部と、前記下界最大化処理部による前記解析対象の３次元構造の事後確率の下界の最大化について収束判定を行う収束判定部と、を備える。

【0027】

（９）本発明の第３の観点は、電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成方法であって、Cryo-TEM法により生成された解析対象が氷包埋された試料を用いて、前記電顕画像からコントラスト伝達関数（CTF）を計算するCTF推定ステップと、前記電顕画像から解析対象を抽出処理する粒子抽出ステップと、前記解析対象の３次元構造を再構成する再構成ステップと、を備え、前記再構成ステップは、前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する下界最大化処理ステップと、与えられた前記投影条件により、投影像と投影像生成行列を生成する投影像生成ステップと、前記下界最大化処理ステップにおける前記解析対象の３次元構造の事後確率の下界の最大化について収束判定を行う収束判定ステップと、を備える。

【0028】

（１０）本発明の第４の観点は、電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から解析対象の３次元構造を推定するシステムにおける３次元再構成プログラムであって、コンピュータに、Cryo-TEM法により生成された解析対象が氷包埋された試料を用いて、前記電顕画像からコントラスト伝達関数（CTF）を計算するCTF推定機能と、前記電顕画像から解析対象を抽出処理する粒子抽出機能と、前記解析対象の３次元構造を再構成する再構成機能と、を実行させ、前記再構成機能は、前記解析対象の３次元構造の事後確率の下界を投影条件の確率分布と前記解析対象の３次元構造に関して最大化する下界最大化処理機能と、与えられた前記投影条件により、投影像と投影像生成行列を生成する投影像生成機能と、前記下界最大化処理機能における前記解析対象の３次元構造の事後確率の下界の最大化について収束判定を行う収束判定機能と、を備える。

【発明の効果】

【0029】

本発明によれば、計算量を抑えつつ、高分解能な３次元再構成像を生成するための解析対象の３次元再構成装置、３次元再構成方法及び３次元再構成プログラムを提供することができる。

【図面の簡単な説明】

【0030】

【図1】本発明の第１の実施形態に係る３次元再構成装置の一例を示すブロック図である。

【図2】本発明の第１の実施形態に係る投影像生成部の処理の一例を説明するためのブロック図である。

【図3】本発明の第１の実施形態に係る投影像生成行列算出部の処理の一例を説明するためのブロック図である。

【図4】図３の本発明の第１の実施形態に係る２次元画素位置生成部の処理の一例を説明するためのイメージ図である。

【図5】図３の本発明の第１の実施形態に係る補助画素生成部の処理の一例を説明するためのイメージ図である。

【図6】図５の本発明の第１の実施形態に係る補助画素生成部の処理の他の例（画素領域内の分割）を説明するためのイメージ図である。

【図7】補間点から画素位置情報(ξ₃,η₃)を生成する処理の一例を説明するためのイメージ図である。

【図8】補間点を３つ用意して座標を変換して投影像を作った例である。

【図9】本発明の第１の実施形態に係る３次元画素位置情報(ξ₅,η₅,ζ₅)の生成方法を説明するためのイメージ図である。

【図10】本発明の第１の実施形態に係る投影像生成行列計算部の処理の一例を説明するためのブロック図である。

【図11】センサーによる撮像のモデル化の一例を説明するためのイメージ図である。

【図12】コントラスト伝達関数（CTF）を説明するイメージ図である。

【図13】透過型電子顕微鏡により撮影された生体分子の投影像のシミュレーション画像である。

【図14】３次元構造の事前確率を変えて３次元構造の推定を行った例である。

【図15】本発明の第１の実施形態に係る下界最大化処理部の処理の一例を説明するためのブロック図である。

【図16】本発明の第４の実施形態に係る下界最大化処理部の処理の一例を説明するためのブロック図である。

【図17】本発明の第６の実施形態に係る下界最大化処理部の処理の一例を説明するためのブロック図である。

【図18】本発明の第７の実施形態に係る３次元再構成装置の下界F(q(z),θ)の最大化の処理を説明するためのフロー図である。

【図19】本発明の第８の実施形態に係る３次元再構成装置の識別器の学習を説明するためのフロー図である。

【図20】第４の投影条件である生体分子がどの方向から投影されているかを表す投影方向α_i,β_i,γ_iと、第５の投影条件である生体分子の中心位置を表す平行移動量x_i,y_i と、を説明するためのイメージ図である。

【図21】生体分子の向きによって投影像が異なることを説明するためのイメージ図である。

【図22】第１の事後確率最大化手段である勾配法において、勾配をモンテカルロ近似によって算出する場合を説明するためのイメージ図である。

【図23】実際に電子顕微鏡により撮影された画像を用いて、本発明によって３次元構造を推定した結果である。

【図24】第１の投影条件の変形を定めるパラメータを用いたアルゴリズムと、第１の投影条件の変形を定めるパラメータを用いないアルゴリズムでの性能の比較である。

【図25】本発明の３次元再構成装置のハードウェア構成の一例を示すブロック図である。

【発明を実施するための形態】

【0031】

（実施形態）
以下、本発明の実施形態について、図面を参照しつつ詳細に説明する。本実施形態において、解析対象は、生体分子、合成分子などの分子や金属を挙げることができる。生体分子としては、ＤＮＡ、ＲＮＡ、タンパク質、糖類、脂質、ＡＴＰなどを挙げることができる。合成分子は、前記生体分子以外の分子であって、天然のものではなく、ヒトの指揮下又はヒトが創作もしくは指揮した制御下で合成されたものをいう。金属としては、典型的な金属だけでなく、金属化合物、Ｓｉなどの半金属や合金をも含む。なお、本発明は以下の実施形態に限定されることはなく、本発明の主旨を逸脱しない範囲で種々変形することが可能である。

【0032】

（第１の実施形態）
まず、本発明の第１の実施形態に係る３次元再構成装置が実現する３次元再構成方法の処理について、説明する。３次元再構成装置の構成は、後述する。本発明の第１の実施形態に係る３次元再構成装置は、電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件とから、生体分子の３次元構造を推定する。生体分子の３次元構造の推定は、生体分子の３次元構造の事後確率の下界の最大化によって行う。生体分子の３次元構造から２次元の投影像を生成するためには、生体分子がどの方向から投影されているかを表す投影方向などの投影条件が必要であるが、単粒子解析において投影条件は未知である。このような未知の要因を隠れ変数と定義し、隠れ変数に関して周辺化した事後確率の下界を最大化することにより、生体分子の３次元構造を推定する。以下、詳細に説明する。

【0033】

粒子画像をy_i、生体分子の３次元構造を表すパラメータをθ（以降、３次元構造パラメータθとする場合がある。）、粒子画像y_iの投影条件すなわち隠れ変数をz_iとしたとき、３次元構造の事後確率は、p(θ|y) ∝ p(y|θ)p(θ)のように、尤度p(y|θ)と事前確率p(θ)の積に比例することから、３次元構造の事後確率の対数をとり、かつ、正規化定数を無視したものは、以下の数式１のようになる

【0034】

【数1】

【0035】

ここで、KLは、カルバック・ライブラー・ダイバージェンスを表す。カルバック・ライブラー・ダイバージェンスは、２つの確率分布の距離を表し、分布が近いほど小さい値をとり、かつ、常に０以上の値となる。またq(z_i)は、隠れ変数z_iに関する任意の確率分布（以降、q(z_i)を試験分布という場合がある。）、p(z_i)は、隠れ変数z_iの事前確率、p(θ)は、３次元構造パラメータθの事前確率を、それぞれ表す。また、粒子画像y_i及び隠れ変数z_iの添字iは、粒子画像の通し番号を表し、Mは、粒子画像枚数を表す。さらに、添え字iがある場合は、i番目の粒子画像を表し、添え字がない場合は、粒子画像の集合を表す。隠れ変数z_iの添え字に関しても同様である。

【0036】

事前確率p(z_i)は、投影条件すなわち隠れ変数z_iに関して事前に何らかの知識がある場合には、その事前知識に沿うように事前確率を定義する。例えば、隠れ変数z_iがある値をとりやすいというような場合には、その値を中心にしたガウス分布を定義するなどが考えられる。また、隠れ変数z_iに関する事前知識がない場合には、一様分布を定義する。

【0037】

カルバック・ライブラー・ダイバージェンスは、常に０以上の値をとることから、生体分子の３次元構造の事後確率の下界F(q(z),θ)は、以下の数式２のように書くことができる（以降、単に下界とする場合がある）。本発明の第１の実施形態に係る３次元再構成装置は、数式２の３次元構造の事後確率の下界F(q(z),θ)を、投影条件の確率分布q(z)と生体分子の３次元構造パラメータθに関して最大化することにより、生体分子の３次元構造を推定する。ここで、投影条件の確率分布q(z)は、投影像が与えられた下での投影条件の尤もらしさを示す確率分布を示す。

【0038】

【数2】

【0039】

３次元構造の事後確率は、以下の数式３のように、下界F(q(z),θ)と、任意の確率分布q(z)と隠れ変数zの事後確率p(z|y,θ)のカルバック・ライブラー・ダイバージェンスとの和で書くことができる。

【0040】

【数3】

【0041】

さらに、３次元構造の事後確率、すなわち、数式３の左辺は、任意の確率分布q(z)に依存しないので、下界F(q(z),θ)のq(z)に関する最大化と、カルバック・ライブラー・ダイバージェンスのq(z)に関する最小化は等価である。したがって、３次元構造の推定は、下界F(q(z),θ)の最大化、又は、カルバック・ライブラー・ダイバージェンスKL[q(z)|p(z|y,θ)]の最小化のどちらで行ってもよい。

【0042】

以上が、本発明の第１の実施形態に係る３次元再構成装置による、３次元構造推定（以降、３次元構造推定を単に最適化という場合がある）の枠組みである。つまり、本発明の第１の実施形態に係る３次元再構成装置は、３次元構造の事後確率の下界F(q(z),θ)をq(z)とθに関して最大化する機能を有する。

【0043】

本発明の第１の実施形態に係る３次元再構成装置では、隠れ変数すなわち投影条件について、生体分子の変形を定めるパラメータd_iを第１の投影条件とし、CTF（コントラスト伝達関数）に関する補正量ρ_iを第２の投影条件とし、粒子画像に目的とする生体分子が適切に写っているか否かを示す量g_iを第３の投影条件とする。また、生体分子がどの方向から投影されているかを表す投影方向α_i,β_i,γ_iを第４の投影条件とし、生体分子の中心位置を表す平行移動量x_i,y_iを第５の投影条件とし、生体分子の対称性を表現する量φ_iを第６の投影条件とし、複数の３次元構造のいずれに属するかを示す量k_i（クラス分類）を第７の投影条件とする。そうすると、隠れ変数z_iは、以下の数式４に示すように、各投影条件を要素とするベクトル量で表される（なお、数式４では、第４から第７及び第１から第３の投影条件の順で記載）。当然のことながら、必ずしもすべての隠れ変数を推定する必要はなく、上記投影条件のいくつかの組み合わせを隠れ変数として推定して、残りを固定値としてもよい。例えば、第１から第３の投影条件の少なくとも１つ、及び／又は、第４から第７の投影条件の少なくとも１つを推定するようにしてもよい。

【0044】

【数4】

【0045】

図２０は、第４の投影条件である生体分子がどの方向から投影されているかを表す投影方向α_i,β_i,γ_iと、第５の投影条件である生体分子の中心位置を表す平行移動量x_i,y_iと、を説明するためのイメージ図である。生体分子がどの方向から投影されているかを表す投影方向は、生体分子の向きで表現する。３次元物体の向きを表すためには３変数を必要とし、α_i,β_i,γ_iとする。生体分子の中心位置を表す平行移動量は、投影方向に並行な軸の移動量は無視し、投影方向に垂直な面内の移動量をx_i,y_iとする。

【0046】

投影像は、３次元構造パラメータθからの３次元物体の生成、生成された３次元物体の投影方向に対する積分、投影面に対する平行移動、生体分子の変形、CTFの適応、センサーによる撮像によって決まるものとして、生成過程を定義する。隠れ変数zすなわち投影条件によって決まる投影像を生成するための行列（以降、投影像生成行列とする場合がある。）をW(z_i)とし、３次元構造をaとしたとき、投影像y_iは、以下の数式５で生成されるものとする。ここで、ε_iは投影像に含まれるノイズを表す。

【0047】

【数5】

【0048】

さらに、粒子画像のノイズε_iが正規分布に従うと仮定すると、尤度p(y_i|z_i,θ)は、以下の数式６のように定義できる。

【0049】

【数6】

【0050】

ここで、Nは、粒子画像y_iの画素数を表し、σ_i ²は、ノイズの分散を表す。３次元構造a、分散σ_i ²は、前述の３次元構造パラメータθ={a,σ_i ²}の構成要素とする。また、Ｔは、行列の転置演算を表す。ここで、３次元構造aはベクトル量である。σ_i ²は、θ={a,σ_i ²}のように便宜上１つのみ書いたが、実際には粒子画像ごとに定義する。当然のことながら、粒子画像同士で共通のσ²としてもよい。数式６では、σ_i ²をスカラー値とし、ノイズの分散は画素位置によらず一定としたが、画素間の相関を考慮して行列Σ_iとし、尤度p(y_i|z_i,θ)を定義してもよい。

【0051】

また、粒子画像のノイズε_iを正規分布以外としてもよい。例えば、ノイズをガウス分布とポアソン分布の混合分布と仮定して最適化を行う手法を採用してもよい。

【0052】

また、数式６では、３次元構造パラメータθ中の３次元構造aは１つとして書いたが、クラス分類の数をＫ個としたときは、３次元構造aをＫ個用意し、クラス分類k_iの値に応じて切り替えて計算する。この場合、３次元構造パラメータθは、θ={a_k,σ_i ²}のようになる。a_kの添え字kはクラス番号を表し、１～Ｋの値をとる。a_kもθ={a_k,σ_i ²}のように便宜上１つのみ書いたが、実際にはクラス分類の数だけ定義する。クラス分類を表す変数をk_iとしたとき、変数k_iは、第７の投影条件である複数の３次元構造のいずれに属するかを示す量を指す。また、分散σ_i ²はクラス間では同一としたが、クラスごとに別々の分散σ_ik ²を定義してもよい。

【0053】

電顕画像には、目的とする生体分子のほかに、別の生体分子や、氷の結晶などが写っている場合がある。また、粒子画像においても、目的とする生体分子が重なっている場合や、生体分子が粒子画像の中心から大きく外れて欠けてしまっている場合もある。このような粒子画像があると、最終的な３次元再構成像の推定精度が悪くなる。

【0054】

目的とする生体分子が通常の生成過程によって投影された粒子画像y_iを正常粒子画像、通常の生成過程によって投影されたとすることはできない粒子画像y_iを異常粒子画像とする。通常の生成過程によって投影されたとはできないとは、前述の別の生体分子又は氷の結晶などが写っている場合や、目的とする生体分子が重なっている場合のほかに、生体分子が粒子画像の中心から大きく外れて欠けてしまっている場合などのことを指す。

【0055】

粒子画像y_iが正常粒子画像であるか異常粒子画像であるかを示す変数（以降、正常粒子画像変数とする場合がある。）をg_i（g_i=1のとき正常粒子画像，g_i=0のとき異常粒子画像）とし、尤度p(y_i|z_i,θ)を以下の数式７のようにおく。

【0056】

【数7】

【0057】

ここで、p_reg(y_i|z_i,θ)は正常粒子画像の尤度で、前述の数式６で計算できる。p_irreg(y_i|z_i,θ)は、異常粒子画像の尤度で、正常粒子画像の生成過程とは別の生成過程で生成されたとして尤度を定義する。正常粒子画像同様にノイズが正規分布に従うと仮定して、数式６のように定義してもよい。その際、通常の生成過程から生成される投影像生成行列W(z_i)とは異なる投影像生成行列W’(z_i)から投影像が生成される。正常粒子画像変数は、第３の投影条件である粒子画像に目的とする生体分子が適切に写っているか否かを示す量g_iを指す。

【0058】

さらに、３次元構造パラメータθは、異常粒子画像の３次元構造も含めて、θ={a_k,a’,σ_i ²}とする。a_kは正常粒子画像の３次元構造で、a’は異常粒子画像の３次元構造である。クラス分類k_i及び正常粒子画像変数g_iの値に応じて、３次元構造aを切り替えて尤度の計算を行う。例えば、正常粒子画像変数g_iがg_i=1で、かつ、クラス分類k_iがkのときは、正常粒子画像の３次元構造a_kを使用し、正常粒子画像変数g_iがg_i=0のときは、異常粒子画像の３次元構造a’を使用し、尤度の計算を行う。

【0059】

また、異常粒子画像の尤度p_irreg(y_i|z_i,θ)を平均０のガウス分布のように定義してもよい。この場合は、投影像は投影条件や３次元構造に依存せず、ノイズだけが加えられた状態を表す。

【0060】

事後確率を隠れ変数g_iに関してのみ書くと、前述の数式（１）は、以下の数式８のように書ける。このとき、隠れ変数g_iは、０か１かの離散値なので積分を積算にする。

【0061】

【数8】

【0062】

事後確率のg_iに関する最大化によって、以下の数式９が得られる。

【0063】

【数9】

【0064】

３次元構造の事後確率の下界の３次元構造パラメータθに関する最大化は、数式９で得られた＾q(g_i)（＾q：キュー・ハット＜推定値＞を示す）を用いると、以下の数式１０を最大化することになる。

【0065】

【数10】

【0066】

このとき、正常粒子画像の３次元構造パラメータθに依存しない項は無視する。数式１０は、正常粒子画像の対数尤度log p_reg(y_i|g_i,θ)に、粒子画像y_iが正常粒子画像の時の事後確率＾q(g_i=1)を掛け、事前確率log p(θ)を足したものを最大化することを表している。本発明の第１の実施形態に係る３次元再構成装置では、数式１０を最大化することにより、３次元構造を推定する。

【0067】

次に、３次元再構成装置の構成について、説明する（ここでは処理の実行単位を示す各部を基に構成を説明し、ハードウェア面での構成は図２５を参照して後に説明する）。図１は、本発明の第１の実施形態に係る３次元再構成装置１０の一例を示すブロック図である。３次元再構成装置１０は、例えば、CTF推定部１１と、粒子抽出部１２と、再構成部１６（下界最大化処理部１３と、投影像生成部１４と、収束判定部１５とを含む）とを備える。３次元再構成装置１０は、各部に対応して、CTF推定ステップと、粒子抽出ステップと、再構成ステップ（下界最大化処理ステップと、投影像生成ステップと、収束判定ステップとを含む）とを備える３次元再構成方法を実現する。そして、３次元再構成方法は、３次元再構成装置１０の内部又は外部に備えられた、図示しないメモリやハードディスク等の記憶装置にデータを記録し、図示しないCPU（Central Processing Unit）やGPU（Graphic Processing Unit）等の演算装置（コンピュータ）が、各部に対応して、CTF推定機能と、粒子抽出機能と、再構成機能（下界最大化処理機能と、投影像生成機能と、収束判定機能とを含む）とに係る各種プログラムを読み込んで実行することにより、達成される。ここで、プログラムとは、例えば、透過型電子顕微鏡により撮影された電顕画像より生体分子の３次元構造を推定するための３次元再構成プログラムである。

【0068】

なお、透過型電子顕微鏡に搭載されている撮像装置には、動画像を撮影できるものがある。その場合、撮影された動画像のフレームの位置合わせを行った後に足し合わせて１枚の画像としたものを使用し、３次元再構成処理を行ってもよい。

【0069】

以下、図１に示す第１の実施形態に係る３次元再構成装置１０の構成を例にして説明しつつ、３次元再構成装置１０における処理の一例について説明する。

【0070】

図１内の本発明の第１の実施形態に係るCTF推定部１１は、Cryo-TEM法により生成された生体分子が氷包埋された試料を用いて、透過型電子顕微鏡により撮影された電顕画像２３からCTF２５を計算して推定する。CTF２５は、例えばデフォーカス量と球面収差係数の値から計算によって求められる。装置固有の値である球面収差係数はユーザーによって指定され、CTF推定部１１では電顕画像２３ごとのデフォーカス量を計算する。

【0071】

この際、CTFは周波数が高周波になるほど減衰することが知られており、周波数ごとの減衰量を表す関数は包絡関数と呼ばれている（前述の図１２の右図における包絡線envelopeを参照）。推定するCTF２５に、周波数ごとの減衰量を表す包絡関数を適応させてもよい。

【0072】

また、デフォーカス量は、透過型電子顕微鏡により撮影された画像のアモルファス部分のフーリエ変換像と、計算により求められたCTF２５が一致するように推定されることが一般的である。アモルファス部分のフーリエ変換像自体をCTF２５としてもよい。

【0073】

図１内の本発明の第１の実施形態に係る粒子抽出部１２は、電顕画像２３から生体分子（以下、粒子とする場合がある。）を拾い上げる処理（以下、抽出処理とする場合がある。）を行う。抽出処理は、ユーザーによって粒子の場所を指示させる方法（マニュアル抽出処理と呼ばれる）と、ユーザーがいくつかの粒子を抽出処理により電顕画像２３から拾い上げ、拾い上げた粒子画像をテンプレートとして残りの粒子を自動で拾い上げる方法（自動抽出処理と呼ばれる）などがある。第１の実施形態に係る３次元再構成装置１０は、特定の抽出処理に限定されるものではない。粒子抽出部１２は、入力された複数枚の電顕画像２３から粒子を拾い上げ、切り出された粒子画像y_i２１を出力する。

【0074】

なお、先に、図１内の本発明の第１の実施形態に係るCTF推定部１１では、電顕画像２３からCTF２５を推定すると説明したが、図１内の本発明の第１の実施形態に係る粒子抽出部１２の後に、CTF推定部１１の処理を行い、粒子画像ごとにCTF２５を推定するようにしてもよい。

【0075】

図１内の本発明の第１の実施形態に係る再構成部１６は、下界最大化処理部１３と、投影像生成部１４と、収束判定部１５を備える。下界最大化処理部１３では、３次元構造の事後確率の下界F(q(z),θ)をq(z)とθに関して最大化する。投影像生成部１４では、詳しくは後述するように、与えられた投影条件による、投影像１１２と投影像生成行列W(z_i)１１３を生成する。収束判定部１５では、収束判定を行う。収束判定は、前回の３次元構造の事後確率の下界F(q(z),θ)と現在の３次元構造の事後確率の下界F(q(z),θ)の差分絶対値が閾値以下であるかで判定を行い、閾値以下であれば収束したとして処理を終了する。収束していない場合は、下界最大化処理部１３で３次元構造の事後確率の下界の最大化を継続する。下界F(q(z),θ)の計算は、前述の数式２に従って行う。

【0076】

なお、前回の３次元構造の事後確率の下界F(q(z),θ)と現在の３次元構造の事後確率の下界F(q(z),θ)の差分絶対値をローパスフィルターに通した出力を、閾値と比較して、収束判定をしてもよい。これにより、突発的な大きな変動による処理終了を防ぐことができ、安定して処理終了の判断が行える。

【0077】

投影像生成部１４について、詳しく説明する。図２は、本発明の第１の実施形態に係る投影像生成部１４の処理の一例を説明するためのブロック図である。隠れ変数z_i２０とCTF２５を投影像生成行列算出部１１０に入力し、投影像を生成するための行列である投影像生成行列W(z_i)１１３を生成する。行列乗算処理部１１１では、投影像生成行列W(z_i)１１３と３次元構造パラメータθ２２に含まれる３次元構造aの行列の乗算を行って投影像１１２を生成する。

【0078】

投影像生成行列算出部１１０では、図３に示した、本発明の第１の実施形態に係る投影像生成行列算出部１１０の処理の一例を説明するためのブロック図のように、２次元画素位置生成部１４１に入力された隠れ変数z_i２０をもとに、一連の生成過程を経て、投影像生成行列計算部１４８において投影像生成行列W(z_i)１１３を計算する。

【0079】

投影像の生成過程は、３次元構造パラメータθからの３次元物体の生成、生成された３次元物体の投影方向に対する積分、投影面に対する平行移動、生体分子の変形、CTFの適応、センサーによる撮像によって決まるものとして、撮像モデルを定義する。投影像生成行列W(z_i)は、これらひとつひとつの処理の合成によって表されているものとすると、以下の数式１１のように分解することができる。投影像生成行列W(z_i)の(z_i)は、投影像生成行列W(z_i)１１３が隠れ変数z_iの関数となることを表しており、投影像生成行列W(z_i)１１３は、隠れ変数z_iに依存して変化する。以降では行列ひとつひとつの算出方法を説明する。

【0080】

【数11】

【0081】

まず、２次元画素位置生成部１４１では、図４に示した、図３の本発明の第１の実施形態に係る２次元画素位置生成部１４１の処理の一例を説明するためのイメージ図のように、指定された画像サイズW_i,H_iに合わせてGrid状の画素位置情報(ξ₁,η₁)１４９を生成する。画素位置情報(ξ₁,η₁)１４９は投影画像の元画素１６０の位置を表しており、画素位置情報(ξ₁,η₁)１４９はすべての画素位置情報の集合でベクトルである。例えば、W_i=100,H_i=50の場合、画素位置情報(ξ₁,η₁)１４９は、(0,0),(1,0)…(99,0),(0,1)…(99,49)というデータになる。図４の点線矩形領域は画素領域１６１を表す。

【0082】

次に、補助画素生成部１４２では、図５に示した、図３の本発明の第１の実施形態に係る補助画素生成部１４２の処理の一例を説明するためのイメージ図のように、２次元画素位置生成部１４１で生成された画素位置情報(ξ₁,η₁)１４９の画素領域１６１内に補助画素１６２を生成し、元画素１６０と補助画素１６２を合わせたものを画素位置情報(ξ₂,η₂)１５０として出力する。例えば、W_i=100,H_i=50で間隔を０．５とし、図５のように補助画素を置いた場合、画素位置情報(ξ₂,η₂)１５０は、(-0.5,-0.5),(0.5,-0.5)…(0,0),(1,0)…(99,0),(-0.5,0.5),(0.5,0.5)…(99,49)…(99.5,49.5)というデータになる。図５の補助画素１６２の配置は一例であり、図６のように補助画素１６２をより細かく置いてもよい。図６の左図（図５の例）では、画素領域１６１が区間Ｓ_１からＳ_４に４分割された例を、図６の右図では、画素領域１６１が区間Ｓ_１１からＳ_１４、区間Ｓ_２１からＳ_２４、区間Ｓ_３１からＳ_３４及び区間Ｓ_４１からＳ_４４に１６分割された例を、それぞれ示している。

【0083】

図３に戻り、次に、変形処理部１４３では、変形処理を行う。変形処理は、画素位置情報(ξ₂,η₂)１５０を新しい画素位置情報(ξ₃,η₃)１５１に変換することによって行われる。変形処理用のパラメータは、第１の投影条件である変形を定めるパラメータを指し、隠れ変数z_i２０に含まれており（前述の数式４中のd_i）、変形処理部１４３では隠れ変数z_i２０を用いて座標変換後の画素位置情報(ξ₃,η₃)１５１を出力する。変形処理の方法は問わないが、例えば、画素位置情報(ξ₂,η₂)１５０一点一点の変換先を隠れ変数z_i２０としてもよい。

【0084】

また、変形処理は、図７に示した、補間点１６３から画素位置情報(ξ₃,η₃)１５１を生成する処理の一例を説明するためのイメージ図のように、いくつかの補間点から画素位置情報(ξ₂,η₂)１５０を画素位置情報(ξ₃,η₃)１５１に変換する処理によって行ってもよい。補間点１６３を通る曲線１６５を求め、画素位置情報(ξ₂,η₂)１５０から曲線１６５に直交する垂線１６６を下ろし、画素位置情報(ξ₂,η₂)１５０から垂線１６６の交点までの距離をη₃とし、あらかじめ補間点１６３の中から決めておいた始点１６４から垂線１６６の交点までの曲線１６５に沿った距離をξ₃とする。このとき、(ξ₃,η₃)には、画素位置情報(ξ₂,η₂)１５０と曲線１６５及び始点１６４の位置関係１６７により符号を付与する。すなわち、位置関係１６７は、画素位置情報(ξ₂,η₂)１５０が始点１６４の垂線１６６の左かつ曲線１６５の上では(ξ₃,η₃)の符号が（－，＋）、始点１６４の垂線１６６の左かつ曲線１６５の下では（－，－）、始点１６４の垂線１６６の右かつ曲線１６５の上では（＋，＋）、始点１６４の垂線１６６の右かつ曲線１６５の下では（＋，－）となる。曲線１６５はどのような曲線でもよく、単純に補間点１６３間を直線でつないでもよいし、スプライン曲線で表現してもよい。図８は、補間点を３つ用意して座標を変換して投影像を作った例である。画像中の点が補間点を表し、左図が変形無し、右図が変形ありの場合である。

【0085】

図３に戻り、次に、平行移動処理部１４４では、座標の平行移動処理を行う。平行移動用パラメータは、第５の投影条件である生体分子の中心位置を表す平行移動量を指し、隠れ変数z_i２０に含まれているとものする（前述の数式４中のx_i,y_i）。平行移動後の座標位置を画素位置情報(ξ₄,η₄)１５２、平行移動用パラメータをx_i,y_iとすると、以下の数式１２のように変換される。

【0086】

【数12】

【0087】

次に、３次元画素位置生成部１４５では、図９に示した、本発明の第１の実施形態に係る３次元画素位置情報(ξ₅,η₅,ζ₅)１５３の生成方法を説明するためのイメージ図のように、画素位置情報(ξ₄,η₄)１５２から３次元画素位置情報(ξ₅,η₅,ζ₅)１５３を生成する。３次元画素位置情報(ξ₅,η₅,ζ₅)１５３のうちξ₅,η₅は画素位置情報(ξ₄,η₄)１５２のξ₄,η₄と同じとし、ζ₅方向、すなわちξ₄η₄平面に垂直な方向の画素位置情報を生成する。このとき、ζ₅方向に生成する画素位置情報の間隔と個数はあらかじめユーザーにより指定されているものとする。例えば、(ξ₄,η₄)の画素位置情報の数が１００個でζ₅方向に１００個の画素位置情報を生成したとすると、最終的に３次元画素位置情報(ξ₅,η₅,ζ₅)１５３は、１００００点の画素位置情報の集合になる。

【0088】

図３に戻り、次に、回転処理部１４６では、座標の３次元回転処理を行う。３次元回転処理用パラメータは、第４の投影条件である生体分子がどの方向から投影されているかを表す投影方向を指し、隠れ変数z_i２０に含まれているものとし（前述の数式４中のα_i,β_i,γ_i）、パラメータに従って回転処理を行い、回転後の座標位置、画素位置情報(ξ₆,η₆,ζ₆)１５４を生成する。この座標の３次元回転処理は、３次元物体の姿勢を回転させるために行う。３次元物体の姿勢の表現方法としては、例えばオイラー角を使ってもよい。オイラー角は３次元物体の姿勢を３つの角度で表現し、あらかじめ決められた軸の順番で軸回りに回転させることにより、姿勢を制御する。

【0089】

次に、対称性処理部１４７では、生体分子の構造の対称性に基づいて座標の変換を行う。対称性のパラメータは、第６の投影条件である生体分子の対称性を表現する量を指し、隠れ変数z_i２０に含まれているものとする（前述の数式４中のφ_i）。例えば、生体分子にζ軸を中心に１２０度ごとの対称性があるとした場合は、変換後の座標は、以下の数式１３によって求められる。この場合、１２０度が隠れ変数z_iの構成成分φ_iの値になる。

【0090】

【数13】

【0091】

最後に、投影像生成行列計算部１４８では、前段で求めた画素位置情報と隠れ変数z_i２０とCTF２５をもとに投影像生成行列W(z_i)１１３を求める。

【0092】

図１０は、本発明の第１の実施形態に係る投影像生成行列計算部１４８の処理の一例を説明するためのブロック図である。行列W_3d計算部１７１では、３次元構造パラメータθ２２に含まれる３次元構造aから３次元物体を作る行列W_3dを生成する。以降では、３次元構造を基底関数で表現する方法を説明する。これは、生体分子の３次元構造を複数の基底関数の線形和によって表現する機能を指す。

【0093】

画素位置情報(ξ₇,η₇,ζ₇)１５５、基底関数をf_bとしたとき、３次元物体は、以下の数式１４によって計算される。

【0094】

【数14】

【0095】

数式１４中のlは、画素位置情報(ξ₇,η₇,ζ₇)１５５であって３次元空間上の位置を表すベクトルであり、μ_i及びΔd_iは、基底関数f_bのパラメータであって、μ_iは基底の３次元空間上の位置を表すベクトル、Δd_iは基底関数の広がりを表すスカラー値であり、Bは、基底関数の数を表す。μ_iは、基底関数の位置を表しており、実空間における２つの基底関数の距離を計算する際に使用される。例えば、２つの基底関数のパラメータをμ_iとμ_jとしたとき、基底関数の距離は、μ_iとμ_jの距離||μ_i－μ_j||₂を指す。ここで、||…||₂はＬ２ノルムを表す。

【0096】

さらに、a_iは、３次元構造パラメータθの構成要素である３次元構造aの要素を表す。前述のとおり、３次元構造aはベクトル量であり、ベクトルのサイズは基底関数の数であるBになる。ここで、基底関数f_b(l|μ_i,Δd_i)に対応する係数とはa_iを指し、基底関数f_b(l|μ_i,Δd_i)と基底関数f_b(l|μ_j,Δd_j)の係数の差とはa_i－a_jを指す。

【0097】

基底関数は上述の例以外でもよく、例えば、基底関数をデルタ関数とすると単純なボクセル表現と同じになる。数式１４中のV(l)は、係数a_iと基底関数f_bとの線形和で表現されているので、行列に置き換えることができる。

【0098】

行列W_prj計算部１７２では、３次元物体を投影方向に投影するため、行列W_prjを画素位置情報(ξ₄,η₄)１５２と画素位置情報(ξ₅,η₅,ζ₅)１５３と画素位置情報(ξ₆,η₆,ζ₆)１５４から生成する。投影方向に投影するため、行列W_prjは、以下の数式１５のように３次元物体を回転させるための行列W_rotと投影像を生成する行列W_intに分解することができる。

【0099】

【数15】

【0100】

３次元物体を回転させるための行列W_rotは、画素位置情報(ξ₅,η₅,ζ₅)１５３と画素位置情報(ξ₆,η₆,ζ₆)１５４から生成できる。さらに、投影処理は以下の数式１６のように線形式で得られるので、行列に書き直すことが可能である。V’(ξ₅,η₅,ζ₅)は、３次元構造aから生成された３次元物体W_3daをW_rot(z_i)で回転させたものを表す。投影処理では、V’(ξ₅,η₅,ζ₅)をη₅方向に積算を行うことを表している。

【0101】

【数16】

【0102】

ある画素位置の値が別の画素位置に移動する行列の求め方を説明する。以下の数式１７の場合、Vの１番目の要素は、V’の２番目の要素になり、Vの２番目の要素は、V’の３番目の要素になり、Vの３番目の要素は、V’の１番目の要素になる。移動先の要素番号を行、移動元の要素番号を列とし、その場所の値を１にし、これをすべての画素位置に対して行うことにより、変換行列を作成する。すなわち、移動元と移動先の対応関係が決まれば変換行列が生成できる。この例では対応関係が１対１の場合を示したが、複数の移動元の値の線形和から移動先の値が決まるようにすることも可能である。

【0103】

【数17】

【0104】

行列W_tr計算部１７３では、３次元物体を平行移動させるための行列W_trを画素位置情報(ξ₃,η₃)１５１と画素位置情報(ξ₄,η₄)１５２から生成する。

【0105】

行列W_def計算部１７４では、投影画像の変形処理を行う行列W_defを画素位置情報(ξ₂,η₂)１５０と画素位置情報(ξ₃,η₃)１５１から生成する。

【0106】

行列W_ctf計算部１７５では、画像にCTF２５を適応する行列を生成する。ここで、CTF２５は、CTF推定部１１によって推定される。CTF２５は周波数ごとの応答を表した周波数空間上の関数であり、CTF２５の適応を周波数空間上で行う場合は要素ごとの掛け算になる。

【0107】

周波数空間上で行う場合は、フーリエ変換を行う行列、CTF２５とフーリエ変換像との要素ごとの掛け算を行う行列、フーリエ逆変換を行う行列を一つにまとめて行列W_ctfとすることにより、画像にCTF２５を適応する行列を生成することができる。フーリエ変換、フーリエ逆変換は、線形演算のため行列で表現することが可能である。

【0108】

また、実空間上でCTF２５を適応する場合はCTF２５をフーリエ逆変換し、畳み込み演算を行うことになる。この場合も畳み込み演算は、線形演算のため行列に変換することが可能である。

【0109】

さらに、CTF２５の補正項をW_cmp(z_i)とした場合、行列W_ctf(z_i)は、以下の数式１８のようになる。補正項W_cmp(z_i)は、第２の投影条件である電顕画像の撮影条件であるCTF２５に関する補正量を指し、隠れ変数z_i２０に含まれている（前述の数式４中のρ_i）。ここで、W_ctfは、CTF推定部１１によって推定されたCTF２５によって決まる行列である。

【0110】

【数18】

【0111】

例えば、補正項W_cmp(z_i)をスカラー値ρ_iとした場合、数式１８は、以下の数式１９のように書ける。

【0112】

【数19】

【0113】

なお、CTF推定部１１で推定されたデフォーカス量に直接補正を行ってもよい。その場合は、補正されたデフォーカス量から行列W_ctf(z_i)を生成する。

【0114】

行列W_img計算部１７６では、センサーによる撮像によって生成される撮影画像を計算する行列を画素位置情報(ξ₁,η₁)１４９と画素位置情報(ξ₂,η₂)１５０から求め出力する。

【0115】

図１１は、センサーによる撮像のモデル化の一例を説明するためのイメージ図である。画素位置情報(ξ₂,η₂)１５０の３点と、それぞれの位置のW_3dとW_prj(z_i)とW_tr(z_i)とW_def(z_i)とW_ctf(z_i)を適応した値I₁,I₂,I₃１８０、すなわち、以下の数式２０によって求められた画素位置情報(ξ₂,η₂)１５０の値、とから求められる平面１８１を計算し、画素位置情報(ξ₂,η₂)１５０の３点の領域内で平面を積分する。積分により求められた値を画素位置情報(ξ₂,η₂)１５０の３点の領域内の画素値とする。これを図６の本発明の第１の実施形態に係る補助画素１６２の配置例で示した画素領域１６１内の分割されたすべての区間（区間Ｓ_１１からＳ_１４、区間Ｓ_２１からＳ_２４、区間Ｓ_３１からＳ_３４及び区間Ｓ_４１からＳ_４４）に対して行い、積分で求めた画素値の総和を画素領域１６１内の画素値とする。同様の処理をすべての画素領域１６１で行うことで最終的な撮像画像を得る。この処理は、補助画素１６２を作成し投影像を得る処理を指す。

【0116】

画素位置情報(ξ₂,η₂)１５０の３点の領域内での平面１８１の積分は、同領域の体積を求めることと等しいため、線形演算で記述することができ、撮像画像を得るための演算を行列で表現することが可能である。

【0117】

【数20】

【0118】

以上、前述した図３の本発明の第１の実施形態に係る投影像生成行列算出部１１０の一例では、図１０に示したように、３次元構造パラメータθからの３次元物体の生成W_3d、生成された３次元物体の投影方向に対する積分W_prj(z_i)、投影面に対する平行移動W_tr(z_i)、生体分子の変形W_def(z_i)、CTF２５の適応W_ctf(z_i)、センサーによる撮像W_imgの順番で観測モデルを定義した。しかしながらこの適応順番は上記に限られない。例えば、生体分子の変形W_def(z_i)をW_3dの後に定義してもよい。この場合、変形処理は３次元空間に対して行われる。

【0119】

また、CTF２５の適応W_ctf(z_i)をW_prj(z_i)の中に含めてもよい。本来、デフォーカス量は投影方向に対して徐々に変化しているため、投影方向に垂直な平面ごとにCTF２５を変えながら適応することにより、より正確にCTF２５の適応が行える。

【0120】

また、本発明の第１の実施形態に係る投影像生成部１４では、正常粒子画像の投影像の生成過程と、異常粒子画像の投影像の生成過程を定義し、隠れ変数z_iに含まれるg_iの値に応じて投影像生成行列W(z_i)を生成するが、投影像生成部１４の処理について前述した説明では、正常粒子画像／異常粒子画像の別を問わず投影像の生成過程を示した。これに対し、例えば、正常粒子画像の投影像の生成過程では、基底関数の位置μ_iを球状に配置し、異常粒子画像の投影像の生成過程では、基底関数の位置μ_iを円筒状に配置することにより、別々の生成過程を定義することもできる。

【0121】

その他、前述のとおり、異常粒子画像の尤度を平均０のガウス分布と定義した場合は、投影像は投影条件や３次元構造に依存せず、一定値を生成する生成過程を定義すればよい。

【0122】

ところで、図１に示す本発明の第１の実施形態に関する３次元再構成装置１０は、３次元構造の事後確率の下界F(q(z),θ)の最大化によって３次元構造を推定する。下界F(q(z),θ)の最大化処理は、下界最大化処理部１３によって行われる。下界の最大化の手法は問わないが、以降では、ＥＭアルゴリズムの枠組みを利用した下界の最大化の例を示す。ＥＭアルゴリズムの枠組みを利用した下界の最大化では、下界F(q(z),θ)のq(z)に関する最大化と、θに関する下界の最大化を交互に行うことにより、下界F(q(z),θ)の最大化を行う。

【0123】

図１５は、本発明の第１の実施形態に係る下界最大化処理部１３の処理の一例を説明するためのブロック図である。勾配法用事後確率推定部２００では、下界F(q(z),θ)のq(z)に関する最大化を行う。まず初めに、事後確率＾q(q_i=1)の計算を以下の数式２１を用いて行う。この計算は全粒子画像に対して行う。

【0124】

【数21】

【0125】

このとき、尤度p_reg(y_i|θ)と尤度p_irreg(y_i|θ)は、以下の数式２２により、それぞれの生成過程で得られる下界F(q(z),θ)を推定値にする。すなわち、尤度p_reg(y_i|θ)は数式２２の上の段の式の右辺を、尤度p_irreg(y_i|θ)は数式２２の下の段の式の右辺を推定値とする。p(g_i)は事前確率で、例えば一様分布や前回の尤度p_reg(y_i|θ)、p_irreg(y_i|θ)でもよい。

【0126】

【数22】

【0127】

次に、残りの隠れ変数zの確率分布に関して最大化を行う。隠れ変数zの任意の確率分布q(z)をあるパラメータで特徴づけられる分布族から選び、パラメータを最適化することにより下界F(q(z),θ)をq(z)に関して最大化する。確率分布のパラメータとは、例えば、選んだ確率分布がガウス分布の場合、平均と分散のことを指す。

【0128】

確率分布のパラメータの最適化には勾配法を使用する。以下、この勾配法を使用した最大化を、第１の事後確率最大化手段という。確率分布のパラメータをτとしたとき、勾配は、下界F(q(z|τ),θ)をτに関して微分することにより得られる。さらに、勾配は、q(z)による期待値をとる形をしているので、モンテカルロ近似を使用すると、勾配は、以下の数式２３で得られる。

【0129】

【数23】

【0130】

ここで、z_i~q(z_i)は、任意の確率分布q(z_i)からz_iを確率分布に従うようにサンプリングすることを表している。さらに、^－τ_ik（^－τ：タウ・バー＜平均値＞を示す）は、現在の確率分布のパラメータの値を表している。ここで、添え字iは粒子画像番号を表し、添え字kはクラス分類番号を表す。すなわち、パラメータτは、画像ごと、クラス分類ごと、隠れ変数の要素ごとに定義し、すべてのパラメータτに対して勾配を計算する。また、数式２３の尤度p(y_i|z_i..,θ)は、前述の数式６を用いて計算する。

【0131】

さらに、計算された勾配を用いて確率分布のパラメータを以下の数式２４のように更新する。ここで、εはユーザーによってあらかじめ設定されたパラメータである。数式２３と数式２４によるパラメータの更新処理は、定義したすべての確率分布のパラメータに対して行う。

【0132】

【数24】

【0133】

次に、勾配法用３次元構造推定部２０１では、下界F(q(z),θ)の３次元構造パラメータθに関する最大化を行う。下界F(q(z|τ),θ)の３次元構造パラメータθに関する最大化は、尤度p(y_i|z_i,θ)を前述の数式６としたことから解析的に解け、さらに、q(z)に関する期待値をモンテカルロ近似すると、下界F(q(z|τ),θ)のθに関する最大化によって得られるθの更新式は、以下の数式２５と数式２６のようになる。

【0134】

【数25】

【0135】

【数26】

【0136】

このとき、θに関する事前確率p(θ)を以下の数式２７とする。σ_i ²に関する事前確率は一様分布とし、最大化に寄与しないため、p(θ)=p(a)とする。勾配法用３次元構造推定部２０１では、数式２５によって３次元構造aを更新し、数式２６によって分散σ_i ²を更新する。このとき分散σ_i ²の更新はすべての粒子画像に対して行う。

【0137】

【数27】

【0138】

今回の例では下界F(q(z),θ)のθに関する最大化を解析的に解いて行ったが、解析的に解けない場合は勾配法を使ってもよい。勾配法で行う場合は、下界F(q(z),θ)をθに関して偏微分して勾配を求め、勾配の方向にθを更新し、これを収束するまで繰り返すことにより、下界F(q(z),θ)のθに関する最大化が行える。

【0139】

以上が、本発明の第１の実施形態に関する３次元再構成装置の下界最大化処理部１３の処理になる。本発明の第１の実施形態に係る下界最大化処理部１３の処理の一例では、勾配法用事後確率推定部２００と勾配法用３次元構造推定部２０１の処理を交互に行うと説明したが、例えば、勾配法用事後確率推定部２００による下界F(q(z),θ)のq(z)に関する最大化を複数回繰り返した後に、勾配法用３次元構造推定部２０１による下界F(q(z),θ)のθに関する最大化を行ってもよい。

【0140】

また、本実施形態では実空間上の処理として説明を行ったが、同等の値が得られる場合は周波数空間で処理を行ってもよい。CTF２５の処理のように、実空間から周波数空間に変換して処理を行い、処理後実空間に逆変換するという方法でもよい。

【0141】

（第２の実施形態）
例えば、図１に示す第１の実施形態に係る３次元再構成装置１０（又は、第１の実施形態の変形例に係る３次元再構成装置。以下、同様とする。）の構成によって、本発明の実施形態に係る３次元再構成方法に係る処理を実現し、電顕画像から生体分子の３次元再構成像を得ることができる。しかしながら、本発明の実施形態に係る３次元再構成方法に係る処理を実現することが可能な構成は、第１の実施形態に係る３次元再構成装置１０に示す構成に限られない。

【0142】

本発明の第２の実施形態に係る３次元再構成装置は、本発明の第１の実施形態に係る３次元再構成装置１０に対し、図１５の下界最大化処理部１３中の勾配法用事後確率推定部２００で計算される勾配の算出方法の点において、異なる。

【0143】

隠れ変数z_iに関する勾配は、以下の数式２８のようになる。すなわち、第１の実施形態に係る３次元再構成装置１０では、勾配の計算はモンテカルロ近似を使って算出したが、本発明の第２の実施形態に係る３次元再構成装置では、数式２８内の∇τ_ikf(q(z_i),θ)を近似し、近似した値をもとに勾配を求める。∇τ_ikf(q(z_i),θ)は、以下の数式２９のようになる。

【0144】

【数28】

【0145】

【数29】

【0146】

例えば、∇τ_ikf(q(z_i),θ)を一次式で近似する場合、数式２９は、以下の数式３０のように書ける。

【0147】

【数30】

【0148】

ある隠れ変数z_iのときの∇τ_ikf(q(z_i),θ)を数式２９を使って計算し、これを複数の隠れ変数z_iに対して行う。隠れ変数z_iと計算した∇τ_ikf(q(z_i),θ)を用いて数式３０の近似式の係数Aτ_ikとbτ_ikを求める。ここで、係数Aτ_ikはベクトルで、係数bτ_ikはスカラー値である。数式２９を使って計算する∇τ_ikf(q(z_i),θ)の数は、最低限、近似式の係数が求まる数以上とする。例えば、係数Aτ_ikが５次元のベクトルの場合、求める係数の数が６個になるので、６個以上の∇τ_ikf(q(z_i),θ)を計算する。近似式の係数はどのように求めてもよいが、例えば最小二乗法で求めることができる。

【0149】

本発明の第１の実施形態に係る勾配の計算では、モンテカルロ近似を行って勾配を求めるが、勾配を求める式は前述の数式２３なので、Ｓ回尤度を計算する必要がある。尤度の計算には前述の数式６より投影像の生成が必要になる。一般的に近似精度を良くしようとした場合、尤度の計算回数Ｓを大きくする必要がある。一方、本発明の第２の実施形態に係る勾配の計算では、最低限、近似式の係数が求まる数だけ尤度計算を行えばよいので計算量の削減が可能になる。

【0150】

次に、求めた近似式である数式３０を使って以下の数式３１の計算を行い、勾配を求める。このとき、数式３１の積分は確率分布q(z_i|τ_i)の関数形によっては解析的に解けることがある。その場合は、解析的に解いた式を使って勾配を求める。解析的に解けない場合は、区分求積法等の数値積分を行ってもよい。勾配の更新式である前述の数式２４及び事後確率＾q(q_i=1)の計算は、本発明の第１の実施形態に係る３次元再構成装置の下界最大化処理部１３と同様に計算する。

【0151】

【数31】

【0152】

（第３の実施形態）
例えば、図１に示す第１の実施形態の係る３次元再構成装置１０（又は、第１の実施形態の変形例に係る３次元再構成装置。以下、同様とする。）の構成によって、本発明の実施形態に係る３次元再構成方法に係る処理を実現し、電顕画像から生体分子の３次元再構成像を得ることができる。しかしながら、本発明の実施形態に係る３次元再構成方法に係る処理を実現することが可能な構成は、第１の実施形態に係る３次元再構成装置１０に示す構成に限られない。

【0153】

本発明の第３の実施形態に係る３次元再構成装置は、本発明の第１の実施形態に係る３次元再構成装置１０に対し、図１５の下界最大化処理部１３中の勾配法用３次元構造推定部２０１で使用される３次元構造の事前確率p(θ)が実空間上での隣接関係を考慮して定義されている点において、異なる。

【0154】

３次元構造aの事前確率として、実空間上での隣接間の関係性を考慮した事前確率を使用することができる。以下の数式３２は、空間的に近傍の係数aは近い値をとるという３次元構造の事前確率の一例である。

【0155】

【数32】

【0156】

数式３２内のs,tは、係数aの番号であり、s~tは、s,tが空間的に近傍であることを意味している。すなわち、数式３２は、近傍の係数aの差分の２乗和を表している。ここで、係数a_sと係数a_tが空間的に近傍であるかどうかは、係数a_sと係数a_tの基底関数f_b(l|μ_s,Δd_s)とf_b(l|μ_t,Δd_t)の３次元空間上の距離||μ_s－μ_t||₂がユーザーにより指定された閾値以内というような条件で決めてもよい。この処理は、２つの基底関数の距離が近傍であるか判定する機能を指す。また、数式３２内のconstは、係数aとは無関係の項をまとめたものである。

【0157】

また、空間的に近傍の係数aが近い値をとり、エッジを保存するような３次元構造の事前確率として以下の数式３３のようなものを使ってもよい。

【0158】

【数33】

【0159】

数式３３内のη_stは、係数a_sと係数a_tの間にエッジがあるか否かを表すエッジ変数で、０又は１の値をとる。λは、係数a_sと係数a_tの値の差をエッジとするか否かを決めるための閾値で、係数a_sと係数a_tの値の差の２乗がλより小さい場合η_stは１の値をとり、大きい場合は０の値をとる。これは、基底関数に対応する係数aの差がエッジかどうか判定する機能を指す。

【0160】

３次元構造の事前確率p(a)の対数をとったものが、数式３２又は数式３３になるように前述の数式２７のΣ_a ^-1を決めることにより、下界F(q(z),θ)の３次元構造パラメータθに関する最大化は、本発明の第１の実施形態に係る３次元再構成装置の勾配法用３次元構造推定部２０１の前述の数式２５を用いて行うことが可能となる。

【0161】

（第４の実施形態）
例えば、図１に示す第１の実施形態の係る３次元再構成装置１０（又は、第１の実施形態の変形例に係る３次元再構成装置。以下、同様とする。）の構成によって、本発明の実施形態に係る３次元再構成方法に係る処理を実現し、電顕画像から生体分子の３次元再構成像を得ることができる。しかしながら、本発明の実施形態に係る３次元再構成方法に係る処理を実現することが可能な構成は、第１の実施形態に係る３次元再構成装置１０に示す構成に限られない。

【0162】

本発明の第４の実施形態に係る３次元再構成装置は、本発明の第１の実施形態に係る３次元再構成装置１０に対し、図１の下界最大化処理部１３で行われる下界の最大化手法に関して、連続量である投影条件、例えば投影方向等を離散化した投影条件ごとに投影条件の事後確率を計算することにより生体分子の３次元構造の事後確率の最大化を行う点において、異なる。以下、この離散化を考慮した最大化を、第２の事後確率最大化手段という。

【0163】

連続量である隠れ変数z_i（例えば、投影方向等）を一定規則に従って離散的にサンプリングし、サンプリングした条件でのみ任意の確率分布である試験分布q(z)が値を持つように試験分布を定義する。

【0164】

３次元構造の事後確率の下界F(q(z),θ)のq(z)に関する最大化は、試験分布q(z)と隠れ変数z_iの事後確率p(z_i|y_i,θ)に関するカルバック・ライブラー・ダイバージェンスのq(z)に関する最小化と等価なため、以下の数式３４のように、q(z)がp(z_i|y_i,θ)と等しいときに、下界F(q(z),θ)はq(z)に関して最大となる。このとき、試験分布q(z)を上述のように、サンプリングした条件でのみ値を持つように定義したことから、数式３４の分母の積分演算は積算演算になる。

【0165】

【数34】

【0166】

すなわち、これは、連続量である隠れ変数z_i（例えば、投影方向等）を一定規則に従って離散的にサンプリングし、隠れ変数z_iの全組み合わせに対して尤度p(y_i|z_i,θ)と事前確率p(z_i|θ)の積の積算を行って分母を計算することになる。

【0167】

隠れ変数z_iの離散化は、例えば、以下の数式３５のようにすることにより、隠れ変数z_iを、範囲κ_i(n).startからκ_i(n).endの間でκ_i(n).num個生成することができる。数式３５中のnは、隠れ変数z_i２０の要素の番号を表し1~Nzの値をとり、jは、サンプルの番号を表し1~κ_i(n).numの値をとるものとする。さらに、全組み合わせＳは、以下の数式３６によって求められる。例えば、隠れ変数の要素が前述の数式４の場合、ｎ＝１はα_iになり、Nzは１０になる。

【0168】

【数35】

【0169】

【数36】

【0170】

図１６は、本発明の第４の実施形態に係る下界最大化処理部１３の処理の一例を説明するためのブロック図である。離散用事後確率推定部３００では、前述の数式３４を全粒子画像及び全隠れ変数の要素に対して行う。事後確率＾q(q_i=1)の計算は、前述の数式２１を用いて行う。

【0171】

離散用３次元構造推定部３０１では、下界F(q(z),θ)の３次元構造パラメータθに関する最大化を行う。離散用３次元構造推定部３０１では、以下の数式３７によって３次元構造aを更新し、以下の数式３８によって分散σ_i ²を更新する。このとき、分散σ_i ²の更新はすべての粒子画像に対して行う。

【0172】

【数37】

【0173】

【数38】

【0174】

（第５の実施形態）
本発明の第５の実施形態に係る３次元再構成装置は、本発明の第４の実施形態に係る３次元再構成装置に対し、図１６の離散用３次元構造推定部３０１で行われる下界F(q(z),θ)の３次元構造パラメータθに関する最大化手法に関して、事後確率の高い投影条件を選択する機能と、選択されなかった残りの投影条件に関してモンテカルロ近似を行う機能を有する点において、異なる。

【0175】

本発明の第４の実施形態に係る３次元再構成装置では、前述のとおり、３次元構造aの更新式である前述の数式３７内の隠れ変数z_iの組み合わせ数Ｓは前述の数式３８で決まるが、隠れ変数の要素数が増える、又は、分割数κ_i(n).numが増えると組み合わせ数が膨大になる。

【0176】

これに対し、本発明の第５の実施形態に係る３次元再構成装置は、隠れ変数z_iの全組み合わせの中から事後確率p(z_i|y_i,θ)の高い順に一定数を選択し、残りの組み合わせをモンテカルロ近似することにより、下界F(q(z),θ)の３次元構造パラメータθに関する最大化を行う。

【0177】

まず、図１６の離散用事後確率推定部３００で求めた事後確率p(z_i|y_i,θ)を確率の高い値から低い値へ降順にソートし、確率がある一定基準以内の隠れ変数z_iを選択する。基準は、例えば、上位Ｋ番目までや、確率の総和が何％以内までなどがあげられる。

【0178】

次に、総サンプル数をＳ、基準を満たす隠れ変数z_iの数をＫとしたときに、残りのＳ－Ｋ個の隠れ変数z_iから事後確率p(z_i|y_i,θ)に応じてＮ－Ｋ個の隠れ変数z_iをランダムにサンプリングする。

【0179】

３次元構造aと分散σ_i ²の更新式は、以下の数式３９、数式４０のようになる。

【0180】

【数39】

【0181】

【数40】

【0182】

ここで、Ω_iは、基準以内の隠れ変数z_iの集合を表し、z_ij~p(z_i|y_i,θ)＼Ω_iは、基準以外の隠れ変数z_iの集合の中から事後確率p(z_i|y_i,θ)に応じて隠れ変数z_ijをサンプリングすることを表している。また、T_iは、基準以外の隠れ変数z_iの事後確率p(z_i|y_i,θ)の総和を表す。

【0183】

基準を満たす隠れ変数z_iの数がＫで、サンプリングした隠れ変数z_ijの数がＮ－Ｋなので、全体でＮ個の隠れ変数z_iを使って投影像を計算することになる。ＮがＳよりも小さければ計算量の削減が行える。一般的に推定の初めの段階では、隠れ変数すなわち投影条件はいろいろな可能性があるため確率分布の形状は広くなるが、後半になれば投影条件の推定が進み可能性が絞られてくるので確率分布の形状が急峻になる。この場合、ほとんどの条件で確率が０に近くなる。例えば基準を上位Ｋ番目までとした場合に、基準以外の隠れ変数z_iの事後確率p(z_i|y_i,θ)の総和を表すT_iは非常に小さな値になるため、モンテカルロ近似を行っても十分な近似精度が得られ、かつ、計算量の削減ができる。

【0184】

（第６の実施形態）
本発明の第６の実施形態に係る３次元再構成装置は、本発明の第１の実施形態に係る３次元再構成装置１０に対し、図１の下界最大化処理部１３で行われる下界の最大化手法に関して、第１の実施形態における第１の事後確率最大化手段と第４の実施形態における第２の事後確率最大化手段を推定の経過に応じて切り替える機能を有する点において、異なる。

【0185】

本発明の第４の実施形態に係る３次元再構成装置の図１６の下界最大化処理部１３では、前述のとおり、隠れ変数z_iを一定規則に従って離散的にサンプリングし、下界F(q(z),θ)のq(z)に関する最大化を行っている。隠れ変数z_iの組み合わせ数Ｓは前述の数式３８で決まるが、分割数κ_i(n).numが増えると組み合わせ数が膨大になる。

【0186】

また、本発明の第１の実施形態に係る３次元再構成装置１０の図１５の下界最大化処理部１３では、前述のとおり、勾配法によって確率分布のパラメータを更新している。しかしながら、現在の推定値が実際の解から遠い場合に勾配を検出することが難しくなる。また、勾配法は局所的な最適化を行うため局所解に収束してしまう可能性もある。

【0187】

これに対し、本発明の第６の実施形態に係る３次元再構成装置の下界最大化処理部１３では、推定の初期段階では図１６の本発明の第４の実施形態に係る３次元再構成装置の下界最大化処理部１３の処理を使って最適化をし（第４の実施形態における第２の事後確率最大化手段）、推定の後半では図１５の本発明の第１の実施形態に係る３次元再構成装置１０の下界最大化処理部１３の処理（第１の実施形態における第１の事後確率最大化手段）を行う。

【0188】

また、図１６の本発明の第４の実施形態に係る３次元再構成装置の下界最大化処理部１３の処理を行う際は、組み合わせ数が膨大にならないように分割数κ_i(n).numを少なくし粗く推定を行い、推定の後半で行われる図１５の本発明の第１の実施形態に係る３次元再構成装置の下界最大化処理部１３の処理で詳細な推定を行う。

【0189】

図１７は、本発明の第６の実施形態に係る下界最大化処理部１３の処理の一例を説明するためのブロック図である。勾配法用事後確率推定部２００及び勾配法用３次元構造推定部２０１は、図１５の本発明の第１の実施形態に係る３次元再構成装置の下界最大化処理部１３の勾配法用事後確率推定部２００及び勾配法用３次元構造推定部２０１の処理と同じで、離散用事後確率推定部３００及び離散用３次元構造推定部３０１は、図１６の本発明の第４の実施形態に係る３次元再構成装置の下界最大化処理部１３の離散用事後確率推定部３００及び離散用３次元構造推定部３０１の処理と同じである。

【0190】

制御部４００では、制御の切り替えを行う。制御の切り替えは例えば、下界F(q(z),θ)の値が閾値以下だった場合、離散用事後確率推定部３００及び離散用３次元構造推定部３０１の処理を行い、閾値を超えたら勾配法用事後確率推定部２００及び勾配法用３次元構造推定部２０１の処理を行うという方法がある。

【0191】

（第７の実施形態）
本発明の第７の実施形態に係る３次元再構成装置は、本発明の第１の実施形態に係る３次元再構成装置１０に対し、図１の下界最大化処理部１３で行われる下界の最大化手法に関して、第１の実施形態における第１の事後確率最大化手段と第４の実施形態における第２の事後確率最大化手段に加え、以下に説明する第３の事後確率最大化手段を推定の経過に応じて切り替える機能を有する点において、異なる。

【0192】

下界F(q(z),θ)をq(z)とθに関して最大化して３次元構造パラメータθを推定する問題は、以下の数式４１のように書ける。この問題は一般的に大域最適化と言われる。一方、初期値の近傍の最適解を見つける問題は局所最適化と呼ばれる。大域最適化に比べて局所最適化は計算コストが少ないが、推定値局所解に収束する可能性がある。一方、大域最適化は局所解に陥ることはないが、計算量が多くなるという問題がある。また、下界F(q(z),θ)の計算コストが大きい場合、計算量の多さはより問題になる。

【0193】

【数41】

【0194】

本発明の第７の実施形態に係る３次元再構成装置は、第３の事後確率最大化手段である大域最適化と、第１の実施形態における第１の事後確率最大化手段と第４の実施形態における第２の事後確率最大化手段のいずれか又は両方を組み合わせて３次元構造の推定を行う。

【0195】

本発明の第７の実施形態に係る３次元再構成装置では、隠れ変数z_iの中の大域最適化を行う変数をφ_iとする。φ_iは、前述した第６の投影条件である生体分子の対称性を表現する量である。本発明の第７の実施形態に係る３次元再構成装置では、大域最適化を行う変数を１つとしたが、当然ながら大域最適化を行う変数を複数にしてもよい。

【0196】

また、大域最適化の手法は問わないが、本発明の第７の実施形態に係る３次元再構成装置ではベイズ最適化を使って説明する。

【0197】

ベイズ最適化では、獲得関数の値が最大となる変数を次の候補値として使用する。獲得関数は、それまでの候補値と目的関数の値によって更新される。本発明の第７の実施形態に係る３次元再構成装置では、目的関数は３次元構造の事後確率の下界F(q(z),θ)になる。

【0198】

任意の確率分布q(z)を大域最適化で行う変数φに関する確率分布とそれ以外の変数z＼φに関する確率分布に分けると、以下の数式４２のように書ける。ここで、z＼φはφ以外の隠れ変数を表す。

【0199】

【数42】

【0200】

本発明の第７の実施形態に係る３次元再構成装置では、獲得関数α(φ)より取得した候補値^－φ（^－φ：ファイ・バー＜平均値＞を示す）を固定し残りの隠れ変数z＼φの確率分布q(z＼φ)とθに関して下界F(q(z＼φ),θ)を最大化する。最大化した下界の値（目的関数の値）をy_nとすると、以下の数式４３のように書ける。下界の最大化は本発明の第１～６の実施形態又は変形例のいずれかを使って行う。

【0201】

【数43】

【0202】

それまでに候補となったすべての隠れ変数^－φ_nとその下界の値y_nを使って、以下の数式４４のように、次の候補値^－φ_n+1を得る。

【0203】

【数44】

【0204】

図１８は、本発明の第７の実施形態に係る３次元再構成装置の下界F(q(z),θ)の最大化の処理を説明するためのフロー図である。まず、獲得関数の値が最大となる隠れ変数φの候補値^－φ_nを取得する（Ｓ１００）。候補値^－φ_n固定して、下界F(q(z|^－φ_n,τ),θ)のq(z)とθに関する最大化を行う（Ｓ１０１）。次に収束判定を行う（Ｓ１０２）。収束判定は、本発明の第１の実施形態に係る３次元再構成装置１０の収束判定部１５の処理で判定する。収束していない場合は下界最大化処理（Ｓ１０１）に戻る。収束している場合は、候補値^－φ_nと最大化した下界の値y_nを保存する（Ｓ１０３）。次に処理を継続するか判定する（Ｓ１０４）。処理の終了判定はどのようなものでもよいが、例えば一定回数処理を行ったかどうかで判定してもよい。処理を継続する場合は候補値^－φ_nの取得（Ｓ１００）に戻る。

【0205】

（第８の実施形態）
本発明の第８の実施形態に係る３次元再構成装置は、本発明の第１の実施形態に係る３次元再構成装置１０に対し、第３の投影条件である粒子画像に目的とする生体分子が適切に写っているか否かを示す量の事後確率の推定手法において、異なる。

【0206】

本発明の第１の実施形態に係る３次元再構成装置１０では、前述のとおり、第３の投影条件である粒子画像に目的とする生体分子が適切に写っているか否かを示す量の事後確率＾q(g_i=1)について、数式２１を用いて行っている。

【0207】

これに対し、本発明の第８の実施形態に係る３次元再構成装置では、事前に学習させた識別器を使って粒子画像が正常粒子画像であるか識別し、識別器が出力した正常粒子画像である確率を事後確率＾q(g_i=1)の推定値とするものである。

【0208】

図１９は、本発明の第８の実施形態に係る３次元再構成装置の識別器の学習を説明するためのフロー図である。識別器の学習は、全粒子画像からあらかじめ設定した枚数を取り出す（Ｓ２００）。次に、取り出した粒子画像に正常粒子画像又は異常粒子画像のラベル（タグ）を付与する（Ｓ２０１）。さらに、取り出した粒子画像の特徴量を抽出し（Ｓ２０２）、抽出した特徴量を用いて識別器の学習を行う（Ｓ２０３）。

【0209】

ここで、特徴量はどのようなものでもよいが、例えば、スケール不変特徴変換によるSIFT特徴量やフィッシャーベクトル等が考えられる。また、識別器についてもどのようなものでもよいが、混合ガウスモデルGMMやランダムフォレストのようなものがある。

【0210】

また、図１９の本発明の第８の実施形態に係る３次元再構成装置の識別器の学習を説明するためのフロー図では、取り出した粒子画像に正常粒子画像又は異常粒子画像のラベルを付与したが、例えば解析対象の生体分子は未知であるが、異常粒子画像の特徴は既知の場合（例えば、別の生体分子の構造が既知など）は、異常粒子画像の特徴量を事前に抽出しておいてもよい。その場合は、正常粒子画像の特徴量のみを抽出し、異常粒子画像の特徴量は事前に取得しておいたものを使用して識別器を学習させる。

【0211】

又は、正常粒子画像の識別器と、異常粒子画像の識別器を別々に学習させてもよい。正常粒子画像の識別器と、異常粒子画像の識別器を別々に学習させた場合の事後確率は、p_regを正常粒子画像の識別器の出力、p_irregを異常粒子画像用の識別器の出力とすると、以下の数式４５の結果を＾q(g_i=1)の推定値とする。

【0212】

【数45】

【0213】

３次元構造の推定は、学習した識別器を用いて全粒子画像を識別し、識別した結果を＾q(g_i=1)として、本発明の第１～７の実施形態に係る３次元再構成装置の処理を行う。

【0214】

（実施形態の効果）
本発明によれば、計算量を抑えつつ、高分解能な３次元再構成像を生成する３次元再構成装置、３次元再構成方法及び３次元再構成プログラムを提供することができる。具体的には、以下のとおりである。

【0215】

上記した実施形態では、透過型電子顕微鏡により撮影された電顕画像から抽出した複数枚の粒子画像と、電顕画像の撮影条件から生体分子の３次元構造を推定するシステムにおける３次元再構成装置であって、３次元構造の事後確率の下界を投影条件の確率分布及び３次元構造に関して最大化することによって、投影条件と生体分子の３次元構造とを推定する。

【0216】

その際、生体分子の３次元構造を投影することにより得られる投影像に関し、生体分子の変形を定めるパラメータからなる第１の投影条件、電顕画像の撮影条件であるCTFに関する補正量からなる第２の投影条件、及び、粒子画像に目的とする生体分子が適切に写っているか否かを示す量からなる第３の投影条件の少なくとも１つを含む投影条件を与えた３次元再構成装置により、生体分子の３次元構造を推定する。また、生体分子がどの方向から投影されているかを表す投影方向からなる第４の投影条件、生体分子の中心位置を表す平行移動量からなる第５の投影条件、生体分子の対称性を表現する量からなる第６の投影条件、及び、複数の３次元構造のいずれに属するかを示す量（クラス分類）からなる第７の投影条件の少なくとも１つを含む投影条件を与えた３次元再構成装置により、生体分子の３次元構造を推定する。

【0217】

従来手法では変形を定めるパラメータが考慮されていないため、変形がある生体分子の３次元構造を適切に推定することができないが、本実施形態では、第１の投影条件である生体分子の変形を定めるパラメータを推定する投影条件に加えることにより、変形がある生体分子でも精度良く推定することが可能となる。さらに、構造変化のあるような生体分子に対しても、同一生体分子の構造多形として推定が可能になる。

【0218】

また、推定する投影条件に第２の投影条件であるCTFに関する補正量を加えることにより、Ｓ／Ｎ比が非常に低い電顕画像から計算されたCTFを補正しながら生体分子の３次元構造を推定することができ、より高精度な３次元再構成像を得ることが可能になる。また、電顕画像から抽出した粒子画像ごとのCTFは、氷の厚さが均一でないことから粒子画像ごとに異なるため、CTFに関する補正量を加えることにより、粒子画像ごとのCTFのズレも吸収できる。

【0219】

また、推定する投影条件に第３の投影条件である粒子画像に目的とする生体分子が適切に写っているか否かを示す量を加えることにより、投影像の写り方の適切さによって粒子画像を最適化の中で取捨選択でき、これにより、従来手法にあったユーザーによる恣意性を排除できる。

【0220】

本実施形態の一態様として、第１から第３の投影条件の確率分布に関する３次元構造の事後確率の下界の最大化は、投影条件の確率分布を表現するパラメータを、例えば、モンテカルロ近似によって算出した勾配を用いた勾配法によって更新することにより、生体分子の３次元構造の事後確率を最大化させる第１の事後確率最大化手段と、本来連続量である投影条件、例えば投影方向等を離散化した投影条件ごとにその事後確率を計算することにより、生体分子の３次元構造の事後確率を最大化する第２の事後確率最大化手段と、大域的最適化を使用して生体分子の３次元構造の事後確率を最大化する第３の事後確率最大化手段と、を推定の経過に応じて切り替えることによって生体分子の３次元構造の事後確率の最大化を行ってもよい。

【0221】

非特許文献１に記載の技術では推定精度を向上させるために、離散化を細かく行うと計算コストが膨大になるという問題があるが、本実施形態のように、第１の事後確率最大化手段である勾配法では、勾配を例えばモンテカルロ近似によって算出することにより確率の低い条件は少なく評価し、確率の高い条件を多く評価するといった効率のよい評価ができ、結果的に計算量の削減になる。例えば、図２２のように投影方向を水平、垂直１度間隔で離散的に分割した場合、分割数は６４８００となる。図２２中の矢印は、投影方向を表している。ある方向からの投影像と、反対方向からの投影像は同じになることを考慮すると、分割数は半分でよく３２４００となる。これは１粒子あたり３２４００個の投影像を生成する必要があることを表している。仮に、モンテカルロ近似の際のサンプリング数を１０００とすると、パラメータの更新１回あたりで比較すると約３０分の１の計算量になる。勾配法によるパラメータの更新回数を１０回とした場合でも、１０００×１０＝１００００個の投影像を計算すればよく、約３分の１の計算量になる。また、先の例のように離散的に分割した場合は１度以下の精度で推定することはできないが、勾配法の場合はパラメータの更新量を細かくすることにより、１度以下の精度で推定することが可能となる。

【0222】

さらに、一般的に勾配法では局所解に収束する可能性があるが、第１の事後確率最大化手段である勾配法に加えて、第２の事後確率最大化手段である離散化された条件ごと事後確率を推定する手法と、第３の事後確率最大化手段である大域的最適化手法とを組み合わせることにより、推定精度の低下を抑えつつ計算量の削減が行える。

【0223】

本実施形態の一態様として、生体分子の３次元構造の事後確率の下界の最大化において、生体分子の３次元構造を複数の基底関数の線形和によって表現する機能と、実空間における２つの基底関数の距離を求める機能と、２つの基底関数の距離が近傍であるか判定する機能と、基底関数に対応する係数の差を計算する機能と基底関数に対応する係数の差がエッジかどうか判定する機能を有し、すべての基底関数の組み合わせに対して求めた距離と、係数の差と、エッジかどうかの判定結果に基づいた３次元構造の事前確率を使って、３次元構造の事後確率を最大化してもよい。

【0224】

粒子画像のＳ／Ｎ比が低いことに起因する過剰適合の問題を、空間的な隣接関係に基づいた３次元構造の事前確率を導入することにより、過剰適合を防ぎながら過度な抑止にならないように効率的に解決することができる。

【0225】

図１４は、３次元構造の事前確率を変えて３次元構造の推定を行った例である。左図は隣接間を考慮しない場合、中図は隣接間が滑らかになるようにした場合、右図は隣接間を滑らかにしつつエッジを保存した場合である。実線の円内部を見ると、中図と右図の例では隣接間が滑らかになっていることがわかる。さらに、点線の円内部を見ると、中図の例ではエッジがぼやけているのに対し、右図の例ではエッジが保存されていることがわかる。

【0226】

本実施形態の一態様として、投影像の生成において、生体分子の３次元構造を複数の基底関数の線形和によって表現する機能と、補助画素を作成して撮像素子により決まる解像度以上の投影像を生成する機能を有し、生成した投影像を用いて生体分子の３次元構造を推定してもよい。

【0227】

生体分子の３次元構造を基底関数の線形和によって表現することにより、３次元空間の任意の位置での３次元構造の値を取得することができ、投影像を任意の分解能で作成することができる。実際の画素に加え、より細かな補助画素の値も計算することにより、より高解像度な投影像を生成することができ、それにより推定精度を向上させることができる。

【0228】

本実施形態の一態様として、第３の投影条件の粒子画像に目的とする生体分子が適切に写っているか否かを示す量の事後確率の推定において、粒子画像の特徴量を抽出する機能と、粒子画像に目的の生体分子が適切に写っているか否かをラベル付けする機能と、粒子画像に目的の生体分子が適切に写っているか否かを識別する識別器に特徴量とラベルによって学習させる機能と、識別器から出力される目的の生体分子が写っている確率を使って３次元構造の事後確率の下界を投影条件の確率分布と３次元構造に関して最大化する機能とを有してもよい。

【0229】

特徴量とラベルによって学習させ識別器が出力する確率を目的とする生体分子が適切に写っているか否かを示す量の事後確率の推定値として使うことにより、目的とする生体分子が適切に写っていない場合の生成過程を定義することが困難な場合でも、生体分子が適切に写っているか否かを考慮して最適化することができ、従来手法にあったユーザーによる恣意性を排除できる。

【0230】

図２３は、実際に電子顕微鏡により撮影された画像を用いて、本発明によって３次元構造を推定した結果である。図２３の左図が推定した３次元構造で、右図が３次元構造の円筒軸に垂直に切った平面上の断面図である。

【0231】

図２４は、第１の投影条件の変形を定めるパラメータを用いたアルゴリズムと、第１の投影条件の変形を定めるパラメータを用いないアルゴリズムでの性能の比較である。横軸は、複数枚の粒子画像のデータセットに含まれる変形量の度合いを表し、縦軸は、推定誤差をｄＢで表している。測定は、横軸のデータセットごとに１０回行った。図２４は、測定結果を箱ひげ図を使ってグラフに描画したものである。さらに、１０回の平均値も折れ線グラフとして重ねて描画している。図２４中の薄いグレーが第１の投影条件の変形を定めるパラメータを用いたアルゴリズムによる結果で、濃いグレーが第１の投影条件の変形を定めるパラメータを用いないアルゴリズムでの結果になる。図２４のグラフの横軸は、右に行くほど変形度合いが大きくなる。薄いグレーの折れ線グラフは水平になっている。これは粒子画像に含まれる変形が大きくなっても性能が保たれていることを示している。一方、濃いグレーの折れ線グラフは、右下がりになっている。これは、逆に、粒子画像に含まれる変形が大きくなるほど、性能が劣化していることを示している。図２４の一番右の最も大きい変形が含まれているデータセットでは、第１の投影条件の変形を定めるパラメータを用いたアルゴリズムは、第１の投影条件の変形を定めるパラメータを用いないアルゴリズムに比べて、３．８６ｄＢ、性能が向上している。

【0232】

（３次元再構成装置１０のハードウェア構成）
前述した実施形態における３次元再構成装置１０及び各処理の一部又は全部は、ハードウェアで構成されていてもよいし、ＣＰＵ（Central Processing Unit）、又はＧＰＵ（Graphics Processing Unit）等が実行するソフトウェア（プログラム）の情報処理で構成されてもよい。ソフトウェアの情報処理で構成される場合には、前述した実施形態における３次元再構成装置１０の少なくとも一部の機能を実現するソフトウェアを、フレキシブルディスク、ＣＤ－ＲＯＭ（Compact Disc-Read Only Memory）、又はＵＳＢ（Universal Serial Bus）メモリ等の非一時的な記憶媒体（非一時的なコンピュータ可読媒体）に収納し、コンピュータに読み込ませることにより、ソフトウェアの情報処理を実行してもよい。また、通信ネットワークを介して当該ソフトウェアがダウンロードされてもよい。さらに、ソフトウェアがＡＳＩＣ（Application Specific Integrated Circuit）、又はＦＰＧＡ（Field Programmable Gate Array）等の回路に実装されることにより、情報処理がハードウェアにより実行されてもよい。

【0233】

ソフトウェアを収納する記憶媒体の種類は限定されるものではない。記憶媒体は、磁気ディスク、又は光ディスク等の着脱可能なものに限定されず、ハードディスク、又はメモリ等の固定型の記憶媒体であってもよい。また、記憶媒体は、コンピュータ内部に備えられてもよいし、コンピュータ外部に備えられてもよい。

【0234】

図２５は、前述した実施形態における３次元再構成装置１０のハードウェア構成の一例を示すブロック図である。３次元再構成装置１０は、一例として、プロセッサ７１と、主記憶装置７２（メモリ）と、補助記憶装置７３（メモリ）と、ネットワークインタフェース７４と、デバイスインタフェース７５と、を備え、これらがバス７６を介して接続されたコンピュータ７として実現されてもよい。

【0235】

図２５のコンピュータ７は、各構成要素を一つ備えているが、同じ構成要素を複数備えていてもよい。また、図２５では、１台のコンピュータ７が示されているが、ソフトウェアが複数台のコンピュータにインストールされて、当該複数台のコンピュータそれぞれがソフトウェアの同一の又は異なる一部の処理を実行してもよい。この場合、コンピュータそれぞれがネットワークインタフェース７４等を介して通信して処理を実行する分散コンピューティングの形態であってもよい。つまり、前述した実施形態における３次元再構成装置１０は、１又は複数の記憶装置に記憶された命令を１台又は複数台のコンピュータが実行することで機能を実現するシステムとして構成されてもよい。また、端末から送信された情報をクラウド上に設けられた１台又は複数台のコンピュータで処理し、この処理結果を端末に送信するような構成であってもよい。

【0236】

前述した実施形態における３次元再構成装置１０の各種演算は、１又は複数のプロセッサを用いて、又は、ネットワークを介した複数台のコンピュータを用いて、並列処理で実行されてもよい。また、各種演算が、プロセッサ内に複数ある演算コアに振り分けられて、並列処理で実行されてもよい。また、本開示の処理、手段等の一部又は全部は、ネットワークを介してコンピュータ７と通信可能なクラウド上に設けられたプロセッサ及び記憶装置の少なくとも一方により実行されてもよい。このように、前述した実施形態における３次元再構成装置１０は、１台又は複数台のコンピュータによる並列コンピューティングの形態であってもよい。

【0237】

プロセッサ７１は、コンピュータの制御装置及び演算装置を含む電子回路（処理回路、Processing circuit、Processing circuitry、ＣＰＵ、ＧＰＵ、ＦＰＧＡ、又はＡＳＩＣ等）であってもよい。また、プロセッサ７１は、専用の処理回路を含む半導体装置等であってもよい。プロセッサ７１は、電子論理素子を用いた電子回路に限定されるものではなく、光論理素子を用いた光回路により実現されてもよい。また、プロセッサ７１は、量子コンピューティングに基づく演算機能を含むものであってもよい。

【0238】

プロセッサ７１は、コンピュータ７の内部構成の各部又は装置等から入力されたデータやソフトウェア（プログラム）に基づいて演算処理を行い、演算結果や制御信号を各部又は装置等に出力することができる。プロセッサ７１は、コンピュータ７のＯＳ（Operating System）や、アプリケーション等を実行することにより、コンピュータ７を構成する各構成要素を制御してもよい。

【0239】

前述した実施形態における３次元再構成装置１０は、１又は複数のプロセッサ７１により実現されてもよい。ここで、プロセッサ７１は、１チップ上に配置された１又は複数の電子回路を指してもよいし、２つ以上のチップあるいは２つ以上のデバイス上に配置された１又は複数の電子回路を指してもよい。複数の電子回路を用いる場合、各電子回路は有線又は無線により通信してもよい。

【0240】

主記憶装置７２は、プロセッサ７１が実行する命令及び各種データ等を記憶する記憶装置であり、主記憶装置７２に記憶された情報がプロセッサ７１により読み出される。補助記憶装置７３は、主記憶装置７２以外の記憶装置である。なお、これらの記憶装置は、電子情報を格納可能な任意の電子部品を意味するものとし、半導体のメモリでもよい。半導体のメモリは、揮発性メモリ、不揮発性メモリのいずれでもよい。前述した実施形態における３次元再構成装置１０において各種データを保存するための記憶装置は、主記憶装置７２又は補助記憶装置７３により実現されてもよく、プロセッサ７１に内蔵される内蔵メモリにより実現されてもよい。

【0241】

記憶装置（メモリ）１つに対して、複数のプロセッサが接続（結合）されてもよいし、単数のプロセッサが接続されてもよい。プロセッサ１つに対して、複数の記憶装置（メモリ）が接続（結合）されてもよい。前述した実施形態における３次元再構成装置１０が、少なくとも１つの記憶装置（メモリ）とこの少なくとも１つの記憶装置（メモリ）に接続（結合）される複数のプロセッサで構成される場合、複数のプロセッサのうち少なくとも１つのプロセッサが、少なくとも１つの記憶装置（メモリ）に接続（結合）される構成を含んでもよい。また、複数台のコンピュータに含まれる記憶装置（メモリ））とプロセッサによって、この構成が実現されてもよい。さらに、記憶装置（メモリ）がプロセッサと一体になっている構成（例えば、Ｌ１キャッシュ、Ｌ２キャッシュを含むキャッシュメモリ）を含んでもよい。

【0242】

ネットワークインタフェース７４は、無線又は有線により、通信ネットワーク８に接続するためのインタフェースである。ネットワークインタフェース７４は、既存の通信規格に適合したもの等、適切なインタフェースを用いればよい。ネットワークインタフェース７４により、通信ネットワーク８を介して接続された外部装置９Ａと情報のやり取りが行われてもよい。なお、通信ネットワーク８は、ＷＡＮ（Wide Area Network）、ＬＡＮ（Local Area Network）、ＰＡＮ（Personal Area Network）等の何れか、又は、それらの組み合わせであってよく、コンピュータ７と外部装置９Ａとの間で情報のやり取りが行われるものであればよい。ＷＡＮの一例としてインターネット等があり、ＬＡＮの一例としてＩＥＥＥ８０２．１１やイーサネット（登録商標）等があり、ＰＡＮの一例としてＢｌｕｅｔｏｏｔｈ（登録商標）やＮＦＣ（Near Field Communication）等がある。

【0243】

デバイスインタフェース７５は、外部装置９Ｂと直接接続するＵＳＢ等のインタフェースである。

【0244】

外部装置９Ａはコンピュータ７とネットワークを介して接続されている装置である。外部装置９Ｂはコンピュータ７と直接接続されている装置である。

【0245】

外部装置９Ａ又は外部装置９Ｂは、一例として、入力装置であってもよい。入力装置は、例えば、カメラ、マイクロフォン、モーションキャプチャ、各種センサー、キーボード、マウス、又はタッチパネル等のデバイスであり、取得した情報をコンピュータ７に与える。また、パーソナルコンピュータ、タブレット端末、又はスマートフォン等の入力部とメモリとプロセッサを備えるデバイスであってもよい。

【0246】

また、外部装置９Ａ又は外部装置Ｂは、一例として、出力装置でもよい。出力装置は、例えば、ＬＣＤ（Liquid Crystal Display）、ＣＲＴ（Cathode Ray Tube）、ＰＤＰ（Plasma Display Panel）、又は有機ＥＬ（Electro Luminescence）パネル等の表示装置であってもよいし、音声等を出力するスピーカ等であってもよい。また、パーソナルコンピュータ、タブレット端末、又はスマートフォン等の出力部とメモリとプロセッサを備えるデバイスであってもよい。

【0247】

また、外部装置９Ａまた外部装置９Ｂは、記憶装置（メモリ）であってもよい。例えば、外部装置９Ａはネットワークストレージ等であってもよく、外部装置９ＢはＨＤＤ等のストレージであってもよい。

【0248】

また、外部装置９Ａ又は外部装置９Ｂは、前述した実施形態における３次元再構成装置１０の構成要素の一部の機能を有する装置でもよい。つまり、コンピュータ７は、外部装置９Ａ又は外部装置９Ｂの処理結果の一部又は全部を送信又は受信してもよい。

【0249】

（本明細書等の解釈）
本明細書及び特許請求の範囲の解釈にあたっては、以下の定義が、使用される。本明細書（請求項を含む）において、「a、b及びcの少なくとも1つ（一方）」又は「a、b又はcの少なくとも1つ（一方）」の表現（同様な表現を含む）が用いられる場合は、a、b、c、a-b、a-c、b-c、又はa-b-cのいずれかを含む。また、a-a、a-b-b、a-a-b-b-c-c等のように、いずれかの要素について複数のインスタンスを含んでもよい。さらに、a-b-c-dのようにdを有する等、列挙された要素（a、b及びc）以外の他の要素を加えることも含む。

【0250】

本明細書（請求項を含む）において、「データを入力として／データに基づいて／に従って／に応じて」等の表現（同様な表現を含む）が用いられる場合は、特に断りがない場合、各種データそのものを入力として用いる場合や、各種データに何らかの処理を行ったもの（例えば、ノイズ加算したもの、正規化したもの、各種データの中間表現等）を入力として用いる場合を含む。また「データに基づいて／に従って／に応じて」何らかの結果が得られる旨が記載されている場合、当該データのみに基づいて当該結果が得られる場合を含むとともに、当該データ以外の他のデータ、要因、条件、及び／又は状態等にも影響を受けて当該結果が得られる場合をも含み得る。また、「データを出力する」旨が記載されている場合、特に断りがない場合、各種データそのものを出力として用いる場合や、各種データに何らかの処理を行ったもの（例えば、ノイズ加算したもの、正規化したもの、各種データの中間表現等）を出力とする場合も含む。

【0251】

本明細書（請求項を含む）において、「接続される（connected）」及び「結合される（coupled）」との用語が用いられる場合は、直接的な接続／結合、間接的な接続／結合、電気的（electrically）な接続／結合、通信的（communicatively）な接続／結合、機能的（operatively）な接続／結合、物理的（physically）な接続／結合等のいずれをも含む非限定的な用語として意図される。当該用語は、当該用語が用いられた文脈に応じて適宜解釈されるべきであるが、意図的に或いは当然に排除されるのではない接続／結合形態は、当該用語に含まれるものして非限定的に解釈されるべきである。

【0252】

本明細書（請求項を含む）において、「ＡがＢするよう構成される（A configured to B）」との表現が用いられる場合は、要素Ａの物理的構造が、動作Ｂを実行可能な構成を有するとともに、要素Ａの恒常的（permanent）又は一時的（temporary）な設定（setting/configuration）が、動作Ｂを実際に実行するように設定（configured/set）されていることを含んでよい。例えば、要素Ａが汎用プロセッサである場合、当該プロセッサが動作Ｂを実行可能なハードウェア構成を有するとともに、恒常的（permanent）又は一時的（temporary）なプログラム（命令）の設定により、動作Ｂを実際に実行するように設定（configured）されていればよい。また、要素Ａが専用プロセッサ又は専用演算回路等である場合、制御用命令及びデータが実際に付属しているか否かとは無関係に、当該プロセッサの回路的構造が動作Ｂを実際に実行するように構築（implemented）されていればよい。

【0253】

本明細書（請求項を含む）において、含有又は所有を意味する用語（例えば、「含む（comprising/including）」及び有する「（having）等）」が用いられる場合は、当該用語の目的語により示される対象物以外の物を含有又は所有する場合を含む、open-endedな用語として意図される。これらの含有又は所有を意味する用語の目的語が数量を指定しない又は単数を示唆する表現（a又はanを冠詞とする表現）である場合は、当該表現は特定の数に限定されないものとして解釈されるべきである。

【0254】

本明細書（請求項を含む）において、ある箇所において「１つ又は複数（one or more）」又は「少なくとも１つ（at least one）」等の表現が用いられ、他の箇所において数量を指定しない又は単数を示唆する表現（a又はanを冠詞とする表現）が用いられているとしても、後者の表現が「１つ」を意味することを意図しない。一般に、数量を指定しない又は単数を示唆する表現（a又はanを冠詞とする表現）は、必ずしも特定の数に限定されないものとして解釈されるべきである。

【0255】

本明細書において、ある実施例の有する特定の構成について特定の効果（advantage/result）が得られる旨が記載されている場合、別段の理由がない限り、当該構成を有する他の１つ又は複数の実施例についても当該効果が得られると理解されるべきである。但し当該効果の有無は、一般に種々の要因、条件、及び／又は状態等に依存し、当該構成により必ず当該効果が得られるものではないと理解されるべきである。当該効果は、種々の要因、条件、及び／又は状態等が満たされたときに実施例に記載の当該構成により得られるものに過ぎず、当該構成又は類似の構成を規定したクレームに係る発明において、当該効果が必ずしも得られるものではない。

【0256】

本明細書（請求項を含む）において、「最大化（maximize）」等の用語が用いられる場合は、グローバルな最大値を求めること、グローバルな最大値の近似値を求めること、ローカルな最大値を求めること、及びローカルな最大値の近似値を求めることを含み、当該用語が用いられた文脈に応じて適宜解釈されるべきである。また、これら最大値の近似値を確率的又はヒューリスティックに求めることを含む。同様に、「最小化（minimize）」等の用語が用いられる場合は、グローバルな最小値を求めること、グローバルな最小値の近似値を求めること、ローカルな最小値を求めること、及びローカルな最小値の近似値を求めることを含み、当該用語が用いられた文脈に応じて適宜解釈されるべきである。また、これら最小値の近似値を確率的又はヒューリスティックに求めることを含む。同様に、「最適化（optimize）」等の用語が用いられる場合は、グローバルな最適値を求めること、グローバルな最適値の近似値を求めること、ローカルな最適値を求めること、及びローカルな最適値の近似値を求めることを含み、当該用語が用いられた文脈に応じて適宜解釈されるべきである。また、これら最適値の近似値を確率的又はヒューリスティックに求めることを含む。

【0257】

本明細書（請求項を含む）において、複数のハードウェアが所定の処理を行う場合、各ハードウェアが協働して所定の処理を行ってもよいし、一部のハードウェアが所定の処理の全てを行ってもよい。また、一部のハードウェアが所定の処理の一部を行い、別のハードウェアが所定の処理の残りを行ってもよい。本明細書（請求項を含む）において、「１又は複数のハードウェアが第１の処理を行い、前記１又は複数のハードウェアが第２の処理を行う」等の表現が用いられている場合、第１の処理を行うハードウェアと第２の処理を行うハードウェアは同じものであってもよいし、異なるものであってもよい。つまり、第１の処理を行うハードウェア及び第２の処理を行うハードウェアが、前記１又は複数のハードウェアに含まれていればよい。なお、ハードウェアは、電子回路、又は電子回路を含む装置等を含んでよい。

【0258】

本明細書（請求項を含む）において、複数の記憶装置（メモリ）がデータの記憶を行う場合、複数の記憶装置（メモリ）のうち個々の記憶装置（メモリ）は、データの一部のみを記憶してもよいし、データの全体を記憶してもよい。

【0259】

以上、本開示の実施形態について詳述したが、本開示は上記した個々の実施形態に限定されるものではない。特許請求の範囲に規定された内容及びその均等物から導き出される本発明の概念的な思想と趣旨を逸脱しない範囲において種々の追加、変更、置き換え及び部分的削除等が可能である。例えば、前述した全ての実施形態において、数値又は数式を説明に用いている場合は、一例として示したものであり、これらに限られるものではない。また、実施形態における各動作の順序は、一例として示したものであり、これらに限られるものではない

【符号の説明】

【0260】

７…コンピュータ
７１…プロセッサ
７２…主記憶装置（メモリ）
７３…補助記憶装置（メモリ）
７４…ネットワークインタフェース
７５…デバイスインタフェース
７６…バス
８…通信ネットワーク
９Ａ，９Ｂ…外部装置
１０…３次元再構成装置
１１…ＣＴＦ推定部
１２…粒子抽出部
１３…下界最大化処理部
１４…投影像生成部
１５…収束判定部
１６…再構成部
１１０…投影像生成行列算出部
１１１…行列乗算処理部
１４１…２次元画素位置生成部
１４２…補助画素生成部
１４３…変形処理部
１４４…平行移動処理部
１４５…３次元画素位置生成部
１４６…回転処理部
１４７…対称性処理部
１４８…投影像生成行列計算部
２００…勾配法用事後確率推定部
２０１…勾配法用３次元構造推定部
３００…離散用事後確率推定部
３０１…離散用３次元構造推定部
４００…制御部

【図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】