The present invention relates to a solution for solving an ill-posed inverse problem in image analysis, e.g. in an electron tomography application in order to recover a structure of a sample. The solution is provided for instance as a method comprising steps of determining reliable prior knowledge about the solution, determining initial guess for the solution and determining the corresponding forward operator, deciding upon model of stochasticity, deciding on suitable regularization method, deciding on updating scheme, and producing a sequence using the set configuration.
|
1. A method for providing reliable solutions to multicomponent inverse problems comprising the steps of:
carrying out each of the following steps in a processing unit of a computer;
given a multicomponent inverse problem, determining reliable prior knowledge about its solution, determining an initial guess for a solution, and determining a corresponding forward operator;
given measured data, deciding upon a model for the stochasticity of the data;
deciding upon which component-wise reconstruction/regularization method to use for recovering each component separately based on said prior knowledge of the solution, said model for the stochasticity of the data, and said forward operator;
deciding upon an updating scheme for the component-wise prior knowledge of the solution based on said component-wise reconstruction/regularization method;
deciding upon a component-wise updating scheme for parameters of said component-wise regularization/reconstruction method based on said model for the stochasticity of the data, said pre-processed data, and said updated component-wise prior knowledge of the solution; and
producing a sequence which is obtained by intertwining two sequences, the first of the sequences from an iterated reconstruction/regularization method and the second of the sequences from a component-wise reconstruction/regularization method using said reconstruction/regularization method for each component separately and applying said updating schemes for the component-wise prior knowledge of the solution and reconstruction/regularization parameters.
2. The method according to
3. The method according to
4. The method according to
5. The method according to
6. The method according to
7. The method according to
8. The method according to
9. The method according to
10. A data acquisition device comprising the processing unit, a non-transient computer readable memory unit, and an interface unit, wherein the processing unit is arranged to run the reconstruction method according to
11. A computational device with a non-transient computer readable memory medium comprising data representing a solution of the inverse problem obtained from the method according to
|
The present invention relates to solving ill-posed inverse problems (typically those that occur in image analysis and image enhancement) where there are several different quantities that are among the unknowns and where the error level in the data is not known beforehand. Such an inverse problem occurs in particular when one uses the transmission electron microscope in a tomographic setting (electron tomography) in order to recover the structure of the sample.
Here we provide the background of the invention in a less formalized language, stressing ideas rather than mathematical rigour. The purpose is to provide the big picture and the mathematical formalization of all concepts introduced here is given in subsection 5.1.
A computational approach to many scientific problems is often based on mathematical models idealizing how input data is transformed into output data. In this context, a forward problem can be characterized as a problem where the goal is to determine (or predict) the outcome when the models are applied to the input data, whereas in an inverse problem the goal is to determine the causes that, using the models, yields the output data.
Most inverse problems, especially the ones that we consider, can be formulated as the problem of solving an (operator) equation. More formally, assume that the object of interest in the inverse problem can be represented, or modeled, by an elements in a suitable set . Assume further that one cannot directly determine fε
by measurements; instead one has a model of an experiment involving f describing how input data is transformed into output data. This model, derived using information applicable to the experiment, can be represented by an operator T that maps elements in
into another set
, called the data space. In reality an experiment only yields finitely many, say m, noisy data points, so
is in this case an m-dimensional vector space. Hence, formally we have T:
and the inverse problem can be stated as the problem of solving the equation
T(f)=g where gεis given.
The corresponding forward problem is to calculate g when f is given. To account for the stochasticity in the measurement process, one also needs to view data as a sample of some random variable and not as a fixed point. In this context, the inverse problem would be the problem to estimate fε from the data, which now is a single sample of some random element with values in the data space
. A method that claims to (approximately) solve the inverse problem will be referred to as a reconstruction method.
A vast majority of inverse problems that are of interest are ill-posed. Intuitively, an inverse problem is well-posed (i.e. not ill-posed) if it has a unique solution for all (admissible) data and the solution depends continuously on the data. The last criterion simply ensures that the problem is stable, i.e. small errors in the data are not amplified. If one has non-attainable data then there are no solutions, i.e. we have non-existence. Next, there can be solutions but they are not unique, i.e. we have non-uniqueness. Finally, even if there is a unique solution, the process of inversion is not continuous, so the solution does not depend continuously on the data and we have instability. If a problem is ill-posed, non-existence is usually not the main concern. In fact, existence of a solution is an important requirement that can many times be achieved by modifying the problem as shown in subsection 5.1.2. Non-uniqueness is considered to be much more serious. The inverse problem with exact data has a unique solution whenever T is injective. In general, there is not much more to say about this case, although for a particular operator T it can of course be a difficult problem to determine if it is injective. The inverse problems occurring in practical applications, however, in general have finite data (so the data space is finite dimensional), while the space
is infinite dimensional, and the forward operator can impossibly be injective. In this case one either has to choose one solution by some criteria or introduce uniqueness by adding additional information. This process is formalized in subsection 5.1.2. Finally, we have the case when one has instability, i.e. when the solution does not depend continuously on the data. This is the main reason for failure of many reconstruction methods that are based on numerically evaluating a discretized version of the inverse of the forward operator. Such an approach works well for well-posed problems, but if the inverse is not continuous, then one will experience numerical instabilities even when the data are only perturbed by small errors. This can partly be dealt with by the use of regularization methods which in general terms replaces an ill-posed problem by a family of neighboring well-posed problems. One has to keep in mind though that no mathematical trick can make an inherently ill-posed problem well-posed. All that a regularization method can do is to recover partial information about the solution in a stable manner. Thus, the “art” of applying regularization methods will always be to find the right compromise between accuracy and stability.
An important class of inverse problems are multicomponent inverse problems where the space naturally splits into a number of components. This splitting of
can be used in deriving efficient reconstruction/regularization methods and regularization methods that take advantage of such a splitting is referred to here as component-wise regularization. An important example of a multicomponent inverse problem is blind deconvolution i.e. one seeks to de-convolve when the convolution kernel is unknown, so both the kernel and the function that is convolved is to be recovered. Another example is the identification problem in emission tomography where both the attenuation map and the activity map are unknown and needs to be reconstructed from measured data [8]. Finally, the inverse problem in electron tomography is an example of a multicomponent inverse problem where a number of experimental parameters must be recovered besides the scattering potential [5]. See also [7] for further examples from electrical engineering, medical and biological imaging, chemistry, robotics, vision, and environmental sciences.
Component-wise regularization/reconstruction, which is a method for dealing with multicomponent inverse problems (see subsection 5.1.6), has been described in the literature, see e.g. the survey article [7] that describes a component-wise non-linear least squares method. Iterated regularization methods have also been descried in the literature, see e.g. [12, 9, 2] for an analysis of iterated Tikhonov regularization. The combination as described here of component-wise regularization/reconstruction and iterated regularization/reconstruction is however not previously described.
It is an object of the present invention to remedy at least some of the problems with the existing technologies and provide an algorithm and method that yields a stable reconstruction of an ill-posed multicomponent inverse problem.
This invention provides a new method to combine iterated regularization with a component-wise regularization. When using a tolerance based variational regularization in any of the components, the invention provides the possibility to update the prior, which in turn redefines the regularization functional, for that component. The invention also provides a new method to iteratively update the regularization parameter for that component and this can be done without knowledge of the total error in the data. Finally, the algorithm is adaptable to a number of multicomponent inverse problem, especially to problems where the space of unknowns splits into two components and the inverse problem is severely ill-posed only in one component.
The present invention is realized in a number of aspects in which a first is a method for providing a stable solution to a multicomponent inverse problem, comprising the steps of:
A second aspect of the present invention, the prior updating within the sequence may be combined with a smoothing operation.
Yet another aspect of the present invention, a data ensemble is provided, comprising the formulation of the inverse problem, defining the forward operator modeling the experiment wherein the inverse problem and forward operator has been derived using the scientific knowledge about how input data is related to output data in the experiment, and a model of the probabilistic distribution of the data;
Yet another aspect of the present invention, a computer readable storage media containing the data ensemble may be provided.
Still another aspect of the present invention, a signal relating to the data ensemble for transportation in a communication network is provided.
In the following the invention will be described in a non-limiting way and in more detail with reference to exemplary embodiments illustrated in the enclosed drawings, in which:
In
The processing device 202, 300 is shown in detail in
Actual experimental data may be pre-processed as to represent data originating from the multicomponent inverse problem in question. The type of pre-processing is determined by the inverse problem and the experimental setup. In a formal sense, unperturbed data are assumed to lie in the range of the forward operator. Other types of pre processing may also be done on the experimental data. The present invention is a new method to process experimental data to get better solutions of multicomponent inverse problems. This algorithm is useful for many problems, including electron microscopy tomography. The algorithm can simultaneously reconstruct several unknown components. Furthermore, the algorithm does not require an accurate estimate of the error level in the data as these can be adaptively updated within an iterative scheme.
The present invention is an algorithm and virtual machine to provide high-quality solutions (reconstructions) of data from a multicomponent inverse problem. Such data are acquired from electron microscope when used in electron tomography (ET) [5] and SPECT when there are no attenuation maps available [8].
The method applies to a large range of multicomponent inverse problems including different data acquisition geometries and stochastic models for the errors in the data. It can be adapt to each of these situations by altering the choices of the component-wise regularization/reconstruction methods and the choices for selecting and updating the corresponding priors and regularization/reconstruction parameters.
The algorithm substantially includes three steps:
In this section we provide the mathematical background necessary for providing an in depth description of the details of the algorithm.
Inverse Problems and Concept of Ill-Posedness
We begin by formally introducing the concept of an inverse problem, the forward problem, and forward operator.
Definition 5.1.
Let be Banach spaces and define data[f] as a random element in
whose distribution depends on fε
. The (measured) data, gε
is defined as a sample of data[f] and the forward operator is the (not necessarily linear) mapping T:
defined by
T(f):=E[data[f]]
provided that the expectation value of data[f] exists. The inverse problem is the problem to estimate fε from the data, i.e. from a single sample of data[f]. The forward problem is to generate a sample in
of data[f] when fε
is given.
Typically, data[f] will be the sum of a random element cT[f], whose distribution depends on T(f), and an independent random element E with a fixed distribution. Another common approach is to consider the measured data g itself as an estimator of T(f), which is different from the setting in Definition 5.1 where the data g is assumed to be a sample of the random variable data[f] with expectation value T(f). The inverse problem is then reduced to the problem of solving the operator equation
T(f)=g. (1)
Definition 5.2.
A reconstruction method will in this paper refer to a method that claims to (approximately) solve the inverse problem in Definition 5.1 (or (1)). Formally it is defined by a reconstruction operator Rλ: and the vector λε
is the parameters of the reconstruction method (it could e.g. be some stepsize or number of iterations).
We intentionally do not put additional requirements on the reconstruction operator. It is however clear that if a reconstruction method is to be of any use, the reconstruction operator should in some sense approximately provide a least squares solution of the inverse problem in question.
Definition 5.3.
We say that the data g in the inverse problem in Definition 5.1 is attainable whenever there exists fε such that T(f)=g. Moreover, we say that we have finite data if
is finite dimensional and otherwise we have infinite data. We have data with only additive error when cT(f)=T(f) (i.e. it is deterministic). If we in addition assume that E≡0, then we have exact data. Finally, we say that the problem is linear in f when the forward operator is a linear function of f.
The formal definition of ill- and well-posedness is due to Hadamard and it is based on the formulation in (1) of the inverse problem.
Definition 5.4.
The inverse problem in (1) is well-posed (or properly posed) if the forward operator is surjective (existence) and injective (uniqueness) when acting on and the inverse of the forward operator is continuous (stability). Otherwise the problem is called ill-posed (or improperly posed).
Remark 5.5.
As stated in [3], in order to apply the above criteria for well-posedness in a specific case, one needs to specify the notion of solution, the range of T (in order to determine which data are considered admissible), and the topology that defines the concept of continuity. Mathematically, one can always make a problem well-posed by choosing specifications in an artificial way, e.g. one can artificially choose the topology in a way that makes T−1 continuous. This is of course irrelevant in an applied problem, where the specifications have to be appropriate for the problem in question.
An important case is when the space in Definition 5.1 naturally splits into a direct sum:
=
V1x . . .
(2)
There are two main reasons for the splitting in (2) of . The first is that the elements in the various components represent different physical entities. A typical example with N=2 could be that
V1 denotes an infinite dimensional Banach space, typically some set of real valued positive functions, that represent an unknown density that is to be recovered whereas
V2 is a finite dimensional vector space defining some experimental parameters that also needs to be recovered. The second reason is that the inverse problem might posses very different stability characteristics when one considers it as an inverse problem only in the
-variable. As an example, it might, be severely ill-posed w.r.t. the elements in
for j=1, . . . , N−1 and well-posed w.r.t. the elements in
We will use the term multicomponent inverse problem to refer to the inverse problem in Definition 5.1 together with the splitting (2).
As already indicated, ill-posed inverse problems require special treatment since an acceptable reconstruction method for such problems must include a regularization. Thus, in our terminology, a regularization method is a reconstruction method whereas there are reconstruction methods that are not regularization methods. In order to develop a reasonable mathematical theory for regularization methods we need of course to restrict the operators T we choose to consider. A basic assumption for a reasonable theory is to assume that
are Banach spaces and T:
is continuous and weakly sequentially closed [4] [3, p. 241].
Existence and Uniqueness of Solutions
Our first task is to relax the notion of a solution of the inverse problem so that existence and uniqueness is achieved in a wider class of problems.
Definition 5.6 (Least Squares Solutions).
Consider the inverse problem given in Definition 5.1 (or (1)). We say that f†εX is a least squares solution whenever
Since an inverse problem can have a least squares solution, even if it does not have an exact solution, this gives us the existence of solutions in a wider class of problems. Note that when data gε is attainable, then the condition in (3) reads as T(f†)=g. Next, we consider the problem of uniqueness. In most cases there are infinitely many solutions to (3), so we need a way to select one. Following [14], this is done by introducing a functional which performs this selection and thereby enforces uniqueness.
Definition 5.7.
Consider the inverse problem given in Definition 5.1 (or (1)) and assume that it has least squares solutions. Also, let ρε be a fixed element which we call the prior and we are given a fixed functional Sρ:
(that depends on the prior ρ). Then f†ε
is a Sρ-minimizing least squares solution if it is a least squares solution that minimizes Sρ, i.e.
Sρ(f†)=inf{Sρ(f):fεfulfills (3)}.
The most common case is when Sρ(f):=∥f−ρ∥, in which case one talks about ρ-minimum norm least squares solutions.
The functional Sρis a measure of the distance to the prior ρ, even though it is not required to define a metric. The idea is that Sρ(f)≧0 with equality when f equals ρ. Another common requirement is that Sρis convex. Remark 5.8. When we have a linear inverse problem, i.e. T is linear and when Sρ(f)=S(f−ρ) for some functional S, then without loss of generality one can choose ρ≡0, since f0† is a S0-minimizing least squares solution of T(f)=g−T(ρ) if and only f0†+ρ is a Sρ-minimizing least squares solution of T(f)=g.
It can happen that the equation Tf=g does not have a least squares solution, but if it does it also has a Sρ-minimizing least squares solution in a wide class of problems and one can prove the following theorem [15].
Theorem 5.9.
Suppose that X and H are Banach spaces, and T:. Finally, assume that the following holds for some ρε
:
Theorem 5.9 provides us with the much needed concept of a unique solution to the inverse problem in Definition 5.1 (or (1)). We are now ready to introduce the notion of a regularization method in a rigorous way and discuss what is commonly included in a mathematical analysis of such a method.
Definition 5.10 (Regularization Method).
A regularization (of the inverse problem in Definition 5.1) is defined as a family {Rλ}λ of continuous maps
Rλ where λε]0,λ01[x . . . x]0,λ0k[,
for some fixed λ0j ε]0, ∞] such that for all gε there exists a parameter choice rule
λ:→]0,λ01[x . . . x]0,λ0k[
such that
for some least squares solution f†ε. The parameter λ(ε,gε) is called the regularization parameter.
We see that each Rλ is a reconstruction operator and the above definition formalizes the fact that a reconstruction method that is a regularization must be stable (Rλ are continuous) and it must have a convergence property in the sense that the regularized solutions converge (in the norm) to a least squares solution of the inverse problem when the error ε in the data goes to zero and the regularization parameter is chosen according to the associated parameter choice rule.
Remark 5.11.
The most common case is when one has a single regularization parameter, i.e. the case k=1 in Definition 5.10. One can also extend the definition of a regularization given in Definition 5.10 so that it incorporates perturbations in the forward operator (modeling errors). In this case T is replaced with some approximation Tδ which is in a δ-neighborhood of T in some appropriate metric. The parameter choice rule must now also include limiting as δ→0. We refer to [1, Section 1.2] for more details on this.
In the inverse problems literature one categorizes the parameter choice rules into two distinct types, a-priori and a-posteriori parameter choice rules.
Definition 5.12.
Let λ:→]0,λ01[x . . . x]0, λ0k[ be a parameter choice rule as in Definition 5.10. If λ does not depend on gε but only on ε, then we say that λ is an a-priori parameter choice rule. Otherwise, we say that λ is an a-posteriori parameter choice rule.
Hence, with an a-priori parameter choice rule the choice of the regularization parameter depends only on the error level ε and not on the actual data gε. As shown in [3, Theorem 3.3], one cannot have a parameter choice rule that depends only on the data gε, since the resulting regularization method is either not stable or not convergent for ill-posed inverse problems.
Variational Regularization Methods
One important class of regularization methods are the variational regularization methods. Even though one can consider cases that makes use of several regularization parameters, we shall see that variational regularization methods always have a natural single regularization parameter. We therefore restrict our description to this case.
Variational regularization methods are reconstruction methods that can be formulated as solving an optimization problem which is defined by a regularization functional1 1Note that the concept of regularization operator which occurs in inverse problems literature is synonymous with the concept of regularization and is therefore not the same as regularization functional.
Sρ:
and a data discrepancy functional
D:
Furthermore, variational regularization methods can be subdivided into tolerance and penalization based methods. The difference lies in the way the data discrepancy functional enters the optimization problem. Penalization based regularization methods are defined as {RS
with λ being the regularization parameter. Tolerance based regularization methods are defined as {RS
with ε being the regularization parameter. In most of the cases D is given as the square of some norm on , but it can also be a log-likelihood function that depends on assumptions regarding the probabilistic distribution of data[f].
Iterated Reconstruction/Regularization Methods
Iterated reconstruction/regularization methods are based on defining a sequence {fj}j⊂ that is obtained by applying a reconstruction/regularization method (which depends on the iterate j) on the inverse problem. Thus, if f0ε
is given, then one forms the iterates {fj}j ⊂
as
fj:=Rλ
The most common case is when Rλ
Component-Wise Reconstruction/Regularization
Let us now consider the inverse problem in Definition 5.1 with the splitting in (2) of . One can of course regularize the inverse problem in Definition 5.1 directly without considering the structure obtained from the splitting (2). This splitting can however be useful in designing a reconstruction/regularization method. Since the various components represent different physical entities, one natural approach is to regularize independently in each
Vj-component. It is however non-trivial to formally show that the resulting reconstruction method is in fact a regularization of the multicomponent inverse problem given by Definition 5.1 and the splitting (2).
Let us more formally describe how to reconstruct/regularize independently in each Vj-component. Let us begin with some notation. Elements in
are denoted by f:=(f1, . . . , fN) where fj ε
Vj. For a given l=1, . . . , N we introduce the following notation:
f|l|:=(f1, . . . ,fl−1,fl+1,fN)εx . . . x
×
x . . . x
and (f|l|;h):=(f1, . . . , fl−1h, fl+1, fN) for hεNow, for fixed f|l|, assume that the family {Rl,λ
Rl,λ→
defines a reconstruction (and if necessary a regularization) method of the corresponding inverse problem where the forward operator is given by
h→T((f|l|;h)) for hεl.
We now define a sequence {fj}j ⊂ in the following way: Assume that f0ε
is given and let σ:
→{1, . . . ,N} denote an index ordering that determines the ordering of the component-wise regularizations. For j≧1 we now define
fj,σ(j):=Rσ(j),λ
fj:=(fj−1[σ(j)];fj,σ(j)). (6)
Mathematical results of reconstruction methods that define a sequence {fj}j ⊂ are scarce and difficult to obtain. It is a highly non-trivial task to derive conditions for when such a reconstruction method is a regularization method, even though each of the maps Rl,λ
-component. Important issues are the question of convergence, i.e. to determine whether the sequence converges and to what it converges. Another issue is to determine the influence of the choice of index ordering σ.
An example of a component-wise reconstruction method is separable non-linear least squares which is nicely reviewed in [7]. This method is applicable to multicomponent inverse problems where the forward operator is a linear combination of nonlinear functions in certain variables. Thus, in this case the space splits into two parts (so N=2 in (2)) where the first part contains the linear variables and the second part contains the nonlinear variables. The idea behind separable nonlinear least squares is to explicitly eliminate the linear variables and then end up with a non-linear least squares problem for recovering the remaining non-linear variables. See [7] for further details.
The Algorithm
The algorithm we are about to describe is designed to solve a multicomponent inverse problem which we assume is given by Definition 5.1 (or in (1)) together with the splitting in (2) of the space . The algorithm yields a sequence which is obtained by intertwining two sequences, the first obtained from an iterated reconstruction/regularization method and the second from a component-wise reconstruction/regularization method. Within this sequence one can also update the component-wise priors and the component-wise regularization/reconstruction parameters. To simplify the notation, we describe the algorithm for the case when the space
splits into two parts (i.e. N=2 in (2)), so we assume that
=
(7)
The extension to the general case is straightforward and simply makes appropriate use of (6).
Elements in will be denoted by (f, h) ε
and the (measured) data is denoted by gε
(which is a sample of data[f, h]). Also, let D:
denote a fixed Hilbert space norm on
. Next, consider the associated inverse problem of recovering the
-component when the
-component is assumed to be known. For fixed hε
and fixed prior ρε
(see Definition 5.7) we therefore let the associated reconstruction/regularization method be denoted by
Rρ,λ→
where λ1(·;ρ,h) is the associated parameter selection rule (which now depends on hεand the prior ρε
), i.e.
λ1(·;ρ,h)→]0,λ01,1[x . . . x]0,λ01;k
for some fixed λ01,j ε[0,∞] with j=1, . . . , k1. Similarly, consider the inverse problem of recovering the -component when the
-component is assumed to be known. For fixed fε
and fixed prior τε
(see Definition 5.7) we therefore let the associated reconstruction/regularization method be denoted by
Rτ,λ→
where λ2(·;τ,f) is the associated parameter selection rule, (which now depends on fεand the prior τε
), i.e.
λ2(·;τ,f):→[0,λ02,1[x . . . x]0,λ02,k
for some fixed λ02,j ε]0,∞] with j=1, . . . , k2. Given a starting point (f0,h0) ε a prior ρε
in the case when the
-prior is not updated and a prior τε
in the case when the
-prior is not updated, the iterates generate a sequence {(fj,hj}j ⊂
defined recursively as follows:
where Fj:for j=1,2 is an operator that acts a smoothing.
Example 5.13.
A typical example is when ⊂
(Ω,
) (for some domain Ω⊂
) denotes a fixed subset representing the feasible
-components. In particular, it is assumed that all functions in
are non-negative and
can e.g. be a finite dimensional vector space representing unknown parameters that need to be recovered as well but that are not the primary interest. If the multicomponent inverse problem is not ill-posed w.r.t. the
-component, then we recover that component by a simple least squares approach. In this case the reconstruction method Rτ,λ
In the case with finite data (so for some m) one often defines the Hilbert space norm D as
D(c,g):=√{square root over ((c−g)t·Σ−1·(c−g))}{square root over ((c−g)t·Σ−1·(c−g))} for c,gε
where we assume that the positive definite m×m-matrix Σ is known (often it represents the covariance matrix of the stochastic variable data[f, h]). This concludes our example of how to recover the -component when the
-component is assumed to be known. Let us now consider the converse, namely to recover the
-component when the
-component is assumed to be known. To make matters interesting, assume that the multicomponent inverse problem is ill-posed w.r.t. the
-component. In our example we choose to use a tolerance based variational regularization method to recover elements in
(when the
V2-component is known), so our regularization method Rρ,λ
) and the Hilbert space norm D in the following way:
Common choices of Sρwhen ⊂
(Ω,
) are
Sρ(f):=∥L(f−ρ)∥2 where L is a linear operator,
The former assignment above yields the class of Tikhonov regularization methods, and the latter yields relative entropy regularization. To conclude our example we need to specify how steps (8) and (9) are defined. In the latter case one defines λ1 (ε, g; ρ, h) and for fixed 0<δ<1, one possibility is to define it as
Note that λ1 (ε,g;ρ,h) does not explicitly depend on ε so we do not assume any prior knowledge of the error level ε in the data. Finally, the smoothing operator F1 in step (8) is given as a smoothing convolution, i.e. F1(f):=ρ*φb where the kernel φb εV1 could e.g. be a fixed low-pass filter with resolution given by the cut-off threshold b>0. Under these assumptions, our iterates defined by steps (8)-(13) take the form
The parameters δ,b>0 are fixed and act as regularization parameters to the method. Setting δ=0 implies that there is no regularization while a choice of δ=1 implies that the data discrepancy is ignored and the solution will be equal to the prior. Increasing b simply means that we smooth the prior more. Roughly speaking, b is chosen so that spatial frequencies up to this threshold are expected to be reconstructed with reasonable stability without much regularization. Some details at higher frequencies can also be recovered, but require regularization to preserve stability. In implementations of the algorithm, the value of b is sometimes decreased as the number of iterations increases. When we keep the prior fixed, then the above iterates are a special case of a component-wise reconstruction method of the inverse problem in Definition 5.1 with the split (7) where we employ tolerance based entropy regularization to the -component and non-linear least squares reconstruction to the
-component. However, in our method one also combines this with an iterated variational regularization scheme that provides has the option of modifying the regularization functional by updating the prior.
Flowchart of Invention
As we shall see, the reconstruction problem in electron tomography (ET) provides an example of a multicomponent inverse problem. We begin with a very brief recollection of the actual experimental setting in using the transmission electron microscope (TEM) in ET as means for structure determination. Next, we state the most common model for the image formation in a TEM which in turn forms the basis for the forward operator. Finally, we formally state the inverse and forward problems in ET.
The Experimental Setting
A specimen that is to be imaged in a TEM must first be physically immobilized since it is imaged in vacuum. Moreover, such specimens also need to be thin (about 100 nm) if enough electrons are to pass through to form an image. The purpose of sample preparation is to achieve this while preserving the structure of the specimen. Sample preparation techniques are rather elaborate and depend on the type of specimen, see [11] (or [6, section 2.2]) and tire references therein for more details. However, from a simplistic point of view one can say that in-vitro specimens are flash-frozen in a millisecond (cryo-EM). In-situ specimens are chemically fixed, cryo-sectioned and immuno-labeled, and can also be treated in a similar way to in-vitro specimens, [13].
Data collection in ET is done by mounting the specimen on a holder (goniometer) that allows one to change its positioning relative to the optical axis of the TEM. For a fixed position, the specimen is radiated with an electron beam and the resulting data, referred to as a micrograph, is recorded by a detector. Hence, each fixed orientation of the specimen yields one micrograph and the procedure is then repeated for a set of different positions. The most common data acquisition geometry is single axis tilting where the specimen plane is only allowed to rotate around a fixed single axis, called the tilt axis, which is orthogonal to the optical axis of the TEM. The rotation angle is called the tilt angle and the angular range is commonly in [−60°, 60°].
The Forward and Inverse Problems in ET
Before dealing with the inverse problem in ET one first needs to properly define the forward problem and the associated forward operator that models the image formation. A proper derivation of the forward operator is outside the scope of this paper and we merely provide a very brief outline on how one arrives at the expression, for the forward operator that occurs in the standard phase contrast model used by the ET community. The interested reader is referred to [5, 10] for a more detailed exposition.
The idea that data in a micrograph can be interpreted as a kind of “projection” of the specimen underlies most models for image formation used in ET. The starting point is to assume that we have perfect coherent imaging, i.e. the incoming electron wave is a monochromatic plane wave (coherent illumination) and electrons only scatter elastically. The scattering properties of the specimen are in this case given by the Coulomb potential and the electron-specimen interaction is modeled by the scalar Schrödinger equation. The picture is completed by adding a description of the effects of the optics and the detector of the TEM, both modeled as convolution operators. However, the basic assumption of perfect coherent imaging must be relaxed. Inelastic scattering introduces partial incoherence which is accounted for within a coherent framework by introducing an imaginary part to the scattering potential, called the absorption potential. The incoherence that stems from incoherent illumination is modeled by modifying the convolution kernel that describes the effect of the optics. Next, as shown in [5], taking the first order Born approximation and linearizing the intensity enables one to explicitly express the measured intensity in terms of the propagation operator acting on the scattering potential of the specimen convolved with point spread functions describing the optics and detector. The standard phase contrast model used by the ET community for the image formation in TEM is based on replacing the propagation operator by its high energy limit as the wave number tends to infinity. This yields a model for the image formation that is based on the parallel beam transform (see (20) for a definition). As we shall see, in many cases one makes one further assumption, namely that the optics and detector are perfect. The resulting model for the image formation can only account for the amplitude contrast, so the main contrast mechanism, namely phase contrast, is not adequately captured by this model for image formation.
The Scattering Potential and Intensity
In order to precisely state the forward operator we introduce the scattering potential f:3→C that defines the structure of the specimen. Following [5], the scattering potential is given as
where m denotes the electron mass at rest, V: is the potential energy2 that models elastic interaction, and Vabs:
is the absorption potential that models the decrease in the flux, due to inelastic scattering, of the non-scattered and elastically scattered electrons. In scattering theory one usually wants a potential that fulfills the Rollnick condition which would be the case when fεε
(Ω,
∩L2(Ω,
). Under the assumptions and approximations outlined in the previous paragraph leading to the standard phase contrast model, the expression for the intensity generated by a single electron is given as 2The potential energy is related to the Coulomb (electrostatic) potential U:
→
by V=−eU where e is the charge of the electron.
for ωτS2 and zεω⊥ where ω⊥:={χε:χ·w=0}, k is the particle wave number3 w.r.t. the homogeneous background medium (which in our case is vacuum), and M denotes the magnification, Moreover, the functions fre,fim:
denote the real and imaginary parts of f, respectively and P denotes the parallel beam transform (X-ray transform) which is defined as 3We use the convention that the relation between the wave number k and the wavelength λ is given by k=2π/λ.
Finally, the point spread functions PSFre and PSFim in (19) model the effect of the optics and incoherent illumination of the TEM. A precise expression of these can e.g. be found in [5] or [10, chapter 65].
The Actual Measured Data and Forward Operator
The expression for the actual data measured on a micrograph needs to account for the counting stochasticity and the detector (usually a slow-scan CCD camera). Following [5], the actual data on a micrograph delivered by the detector from pixel (i,j) is given as a sample of the random variable data[f](ω)i,j which is defined as
where
where Dose(ω) is the incoming dose (in number of electrons hitting the specimen per area unit).
The forward operator in ET, denoted by T, is defined as the expected value of data[f](ω)i,j, i.e.
T(f)(ω)i,j:=E[data[f](ω)i,j] for ωεS2 and pixel (i,j).
From the definition of the forward operator, it is easy to see that
where εi,j:=E[Ei,j]. In many cases it is customary to further simplify the above expression by the following approximation:
where |Δi,j| is the area of (i, j):th pixel Δi,j and zi,j εΔi,j is some suitably chosen point (typically the midpoint). Then, the resulting expression for the forward operator reads as
T(f)(ω)i,j=gaini,j|Δi,j|Dose(ω){PSFdetI(f)(ω,·)}(zi,j)+εi,j. (22)
Definition 5.15.
We have a fixed finite set S0 of directions on a smooth curve S⊂S2 that defines our data collection geometry. The scattering properties of the specimen are assumed to be fully described by the complex valued scattering potential f defined in (18). For each direction ωεS0, the specimen is probed by a monochromatic wave and the resulting intensity is measured in a finite set of pixels (i,j) defined by fixed subsets Δi,j ⊂ω⊥. The forward problem is to generate a sample of data[f](ω)i,j for ωεS0 and the pixels (i,j) when f is given. The inverse problem is to determine f when a sample of data[f](ω)i,j is known for ωεS0 and finitely many pixels (i, j).
Additional Difficulties
The real experimental situation in ET is unfortunately more complicated and a proper formulation of the inverse problem shows that one needs to recover the values of a number of parameters along with the scattering potential. The short exposition here closely follows the one given in [5].
First, there are parameters that have no corresponding interpretation in the actual real-world ET experiment. The optical setup used in deriving the forward operator does not correspond to the actual setup in the TEM and some of the parameters defining this setup enters into the definition of PSFre and PSFim. Formally, different choices of these parameters can define the an optical setup that has the same imaging properties as the actual setup in the TEM. However, there is a natural way to choose these parameters prior to reconstruction and therefore we considered them as known to sufficient degree of accuracy. Next, there are parameters, such as detector related parameters (e.g. gaini,j and the parameters that enters into PSFdet), that can be determined by separate calibration experiments independently of the ET data collection. These are also considered as known prior to reconstruction.
A more difficult class of parameters are those that are independent of the specimen but unique for each ET data collection. Some of these, such as the wavenumber k and magnification M, are in most cases known to sufficient degree of accuracy. However, the defocus Δz, spherical aberration Cs, and some of the envelope parameters have nominal values that must be adjusted, either by an analysis of the recorded micrographs and/or by performing additional measurements after the ET the data collection. The problem of determining these parameters is usually referred to as CTF estimation.
Finally, there are parameters that depend on the specimen. There are two examples of such parameters. The first relates to the fact that in ET one is dealing with a region of interest problem (local tomography), so one has a predefined region of interest (ROI) Ω⊂ and one seeks f in a predefined region of interest (ROI) Ω⊂
which is a strict subset of the support of f. Since Ω is subset of the support of f, when reconstructing f one would have to make assumptions regarding the values of f outside the Ω. A typical such assumption is that f equals some constant average value outside Ω which then becomes a specimen dependent parameter. More formally, define flocal:f|Ω and set Δl(ω,χ) to denote the path length of the part of the line t
tω+χ that lies outside the region of interest (ROI) Ω, but inside the support of f. Let the constant γre+iγim ε
represent the average value of f outside the ROI, so in (19) we need to introduce the following approximation:
P(fre)(ω,χ)≈P(flocalre)(ω,χ)+γreΔl(ω,χ)
P(fim)(ω,χ)≈P(flocalim)(ω,χ)+γimΔl(ω,χ).
Note that equality holds when the line is entirely in Ω. The γre and γim parameters are our specimen dependent parameters. The other example of a specimen dependent parameter is when one assumes that the imaginary part of the scattering potential is proportional to the real part. More formally, assume that
f(χ)=(Qph+iQamp)F(χ) (23)
for constants Qph, Qamp>0 and some function F:. The ratio Qamp/Qph is called the amplitude contrast ratio. Under this assumption, the intensity generated by a single electron becomes
which in turn yields the following expression for the forward operator:
The point spread function PSF depends on Qamp and Qph and is given as
PSF(ω,z):={QphPSFre(ω,·)+Qamp PSFim(ω,·)}(z).
With this assumption one needs to recover the scaling factor (amplitude contrast ratio), which is specimen dependent.
The true inverse problem for ET is a multicomponent inverse problem where we have to recover not only the scattering potential f but also a number of parameters. The space of unknowns in the inverse problem in ET therefore naturally splits in the following way:
=
×Z.
Here, ⊂
(Ω,
) is the set of complex valued functions that have positive real and imaginary parts and in ET we require that
⊂
(Ω,
)∩
(Ω,
). If one assumes that (23) holds at the expense of introducing the amplitude contrast ratio as an unknown parameter, then
is the set of positive real valued functions. Z is a finite dimensional vector space whose elements are the various parameters that are considered unknown prior to reconstruction.
Finally we show an example of how the of the invention on the inverse problem in ET. We create ET data using a simulator based on (21) with the forward operator in (22). The phantom, i.e. the object to be reconstructed, is the RNA polymerase II, and actual scattering potential is calculated from the protein data base information of the RNA polymerase II. The size of the RNA polymerase II is about 400 kDa and it is assumed to be embedded in an amorphous ice slab that is 50 nm thick. The simulator allows generation of images similar to those recorded in an electron microscope. In the simulator tests, an underfocus of 1 μm was used and a full single axis tilt series was produced comprising 121 tilts between −60° and 60°.
The three-dimensional reconstruction from the present innovation is compared both to the original phantom (i.e. the true answer) and to the filtered backprojection reconstruction (FBP) that has been regularized by an additional low-pass filtering (low-pass FBP). This latter filtering, which in our case reduces the resolution to 7 nm, is necessary in order to gain stability and the above value for the low-pass filtering represents the best trade off between stability and resolution if FBP is to be used on this particular example.
We see clearly that the present invention, as compared to the FBP approach, produces reconstructions with higher resolution and where the background clutter is severely suppressed.
Edin, Anders, Öktem, Ozan, Rullgàrd, Hans, Lagerros, Johan, Öfverstedt, Lars-Göran
Patent | Priority | Assignee | Title |
10229092, | Aug 14 2017 | City University of Hong Kong | Systems and methods for robust low-rank matrix approximation |
8921787, | Feb 10 2010 | Mochii, Inc. | Incoherent transmission electron microscopy |
Patent | Priority | Assignee | Title |
7831097, | Mar 16 2006 | Siemens Medical Solutions USA, Inc | System and method for image reconstruction |
20080063247, | |||
WO9733588, |
Executed on | Assignor | Assignee | Conveyance | Frame | Reel | Doc |
Nov 08 2007 | Okinawa Institute of Science and Technology Promotion Corporation | (assignment on the face of the patent) | / | |||
Jun 28 2009 | LAGERROS, JOHAN | Sidec Technologies AB | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 023243 | /0493 | |
Jul 01 2009 | EDIN, ANDERS | Sidec Technologies AB | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 023243 | /0493 | |
Jul 04 2009 | OKTEM, OZAN | Sidec Technologies AB | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 023243 | /0493 | |
Jul 05 2009 | OFVERSTEDT, LARS-GORAN | Sidec Technologies AB | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 023243 | /0493 | |
Jul 13 2009 | RULLGARD, HANS | Sidec Technologies AB | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 023243 | /0493 | |
May 20 2011 | Sidec Technologies AB | Okinawa Institute of Science and Technology Promotion Corporation | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 026358 | /0687 |
Date | Maintenance Fee Events |
Oct 31 2016 | LTOS: Pat Holder Claims Small Entity Status. |
Oct 31 2016 | M2551: Payment of Maintenance Fee, 4th Yr, Small Entity. |
Dec 02 2020 | M2552: Payment of Maintenance Fee, 8th Yr, Small Entity. |
Date | Maintenance Schedule |
Jun 18 2016 | 4 years fee payment window open |
Dec 18 2016 | 6 months grace period start (w surcharge) |
Jun 18 2017 | patent expiry (for year 4) |
Jun 18 2019 | 2 years to revive unintentionally abandoned end. (for year 4) |
Jun 18 2020 | 8 years fee payment window open |
Dec 18 2020 | 6 months grace period start (w surcharge) |
Jun 18 2021 | patent expiry (for year 8) |
Jun 18 2023 | 2 years to revive unintentionally abandoned end. (for year 8) |
Jun 18 2024 | 12 years fee payment window open |
Dec 18 2024 | 6 months grace period start (w surcharge) |
Jun 18 2025 | patent expiry (for year 12) |
Jun 18 2027 | 2 years to revive unintentionally abandoned end. (for year 12) |