Lensfree method for imaging biological samples in three dimensions

Abstract

A method for three-dimensional imaging of a sample (302) comprises: receiving (102) interference patterns (208) acquired using light-detecting elements (212), wherein each interference pattern (208) is formed by scattered light from the sample (302) and non-scattered light from a light source (206; 306), wherein the interference patterns (208) are acquired using different angles between the sample (302) and the light source (206; 306); performing digital holographic reconstruction applying an iterative algorithm to change a three-dimensional scattering potential of the sample (302) to improve a difference between the received interference patterns (208) and predicted interference patterns based on the three-dimensional scattering potential; wherein the iterative algorithm reduces a sum of a data fidelity term and a non-differentiable regularization term and wherein the iterative algorithm includes a forward-backward splitting method alternating between forward gradient descent (108) on the data fidelity term and backward gradient descent (110) on the regularization term.

Claims

1. A method for three-dimensional imaging of a sample, said method comprising: receiving a plurality of interference patterns acquired using light-detecting elements for detecting incident light, wherein each interference pattern is formed by scattered light scattered from the sample in three dimensions and non-scattered light from a light source, wherein the interference patterns are acquired using different angles between the sample and the light source, wherein the plurality of interference patterns includes at least four interference patterns; determining a scattering potential of the sample measured in each of three dimensions based on the received interference patterns; applying an iterative algorithm to change the three-dimensional scattering potential of the sample by calculating a data fidelity term representing a difference between the received interference patterns represented by a measured amplitude of an optical field forming the interference patterns and predicted interference patterns based on the determined three-dimensional scattering potential formed by a light propagation model that maps the three-dimensional scattering potential of the sample to a scattered optical field in a sensor plane of the light-detecting elements and an incident optical field from the light source which is superposed with the scattered optical field to form the predicted interference patterns, wherein the iterative algorithm reduces a sum of the data fidelity term and a non-differentiable regularization term and wherein the iterative algorithm includes a forward-backward splitting method alternating between forward gradient descent on the data fidelity term and backward gradient descent on the regularization term; performing digital holographic reconstruction on the received interference patterns based on the applying; and reconstructing a three-dimensional image representation of the sample based on the digital holographic reconstruction.

2. The method according to claim 1, wherein the iterative algorithm includes a primal-dual splitting method for reducing the regularization term.

3. The method according to claim 1, wherein the forward-backward splitting method comprises a fast iterative shrinkage-thresholding algorithm, FISTA.

4. The method according to claim 1, wherein the forward gradient descent on the data fidelity term comprises solving a phase retrieval problem using Wirtinger derivatives.

5. The method according to claim 1, wherein the plurality of interference patterns comprises four to ten interference patterns.

6. The method according to claim 1, wherein the regularization terms are set for reducing artifacts and shape distortion in the holographic reconstruction caused by a limited number of illumination angles and lack of phase information of a diffracted optical field.

7. The method according to claim 6, wherein the regularization term comprises at least one of a L1 norm, a L2 norm, total variation or bound constraint.

8. The method according to claim 1, wherein the iterative algorithm is applied until a stopping criterion is met.

9. The method according to claim 1, wherein the sample is a three-dimensional organoid.

10. A computer program product comprising computer-readable instructions such that when executed on a processing unit the computer-readable instructions will cause the processing unit to perform the method according to claim 1.

11. A device for three-dimensional imaging of a sample, said device comprising: at least one light source configured to illuminate a sample from a plurality of different angles; at least one image sensor, each comprising an array of light-detecting elements for detecting incident light, wherein the at least one image sensor is configured to acquire a plurality of interference patterns, wherein each interference pattern is formed by scattered light scattered from the sample in three dimensions and non-scattered light from the light source and each interference pattern is acquired for a different angle between the sample and the light source, wherein the plurality of interference patterns includes at least four interference patterns; a processing unit configured to determine a scattering potential of the sample measured in each of three dimensions based on the received interference patterns and apply an iterative algorithm to change the three-dimensional scattering potential of the sample by calculating a data fidelity term representing a difference between the acquired interference patterns represented by a measured amplitude of an optical field forming the interference patterns and predicted interference patterns based on the determined three-dimensional scattering potential formed by a light propagation model that maps the three-dimensional scattering potential of the sample to a scattered optical field in a sensor plane of the light-detecting elements and an incident optical field from the light source which is superposed with the scattered optical field to form the predicted interference patterns, wherein the iterative algorithm reduces a sum of the data fidelity term and a non-differentiable regularization term and wherein the iterative algorithm includes a forward-backward splitting method alternating between forward gradient descent on the data fidelity term and backward gradient descent on the regularization term; a processing unit configured to perform digital holographic reconstruction on the acquired interference patterns based on the application of the iterative algorithm; and a reconstructing unit configured to reconstruct a three-dimensional image representation of the sample based on the digital holographic reconstruction.

12. The device according to claim 11, wherein the device comprises a plurality of fixed light sources for illuminating the sample from the plurality of different angles.

13. The device according to claim 11, wherein the device comprises a movable light source, which is movable between a plurality of positions for illuminating the sample from the plurality of different angles.

14. The device according to claim 11, wherein the device comprises a single image sensor for acquiring the plurality of interference patterns.

15. A method for three-dimensional imaging of a sample, said method comprising: receiving a plurality of interference patterns acquired using light-detecting elements for detecting incident light, wherein each interference pattern is formed by scattered light scattered from the sample in three dimensions and non-scattered light from a light source interfering at the light-detecting elements for acquisition by the light-detecting elements, and wherein the interference patterns are acquired using different angles between the sample and the light source, wherein the plurality of interference patterns includes at least four interference patterns; determining a scattering potential of the sample measured in each of three dimensions based on the received interference patterns; applying an iterative algorithm to change the three-dimensional scattering potential of the sample by calculating a data fidelity term representing a difference between the received interference patterns represented by a measured amplitude of an optical field forming the interference patterns and predicted interference patterns based on the determined three-dimensional scattering potential formed by a light propagation model that maps the three-dimensional scattering potential of the sample to a scattered optical field in a sensor plane of the light-detecting elements and an incident optical field from the light source which is superposed with the scattered optical field to form the predicted interference patterns, wherein the iterative algorithm reduces a sum of the data fidelity term and a non-differentiable regularization term and wherein the iterative algorithm includes a forward-backward splitting method alternating between forward gradient descent on the data fidelity term and backward gradient descent on the regularization term; and performing digital holographic reconstruction on the received interference patterns based on applying the iterative algorithm.

16. A device for three-dimensional imaging of a sample, said device comprising: at least one light source configured to illuminate a sample from a plurality of different angles; at least one image sensor, each comprising an array of light-detecting elements for detecting incident light, wherein the at least one image sensor is configured to acquire a plurality of interference patterns formed at the at least one image sensor, wherein each interference pattern is formed by scattered light scattered from the sample in three dimensions and non-scattered light from a light source interfering at the light-detecting elements for acquisition by the at least one image sensor, wherein the plurality of interference patterns includes at least four interference patterns; and a processing unit configured to determine a scattering potential of the sample measured in at least three dimensions based on the received interference patterns and apply an iterative algorithm to change the three-dimensional scattering potential of the sample by calculating a data fidelity term representing a difference between the acquired interference patterns represented by a measured amplitude of an optical field forming the interference patterns and predicted interference patterns based on the determined three-dimensional scattering potential formed by a light propagation model that maps the three-dimensional scattering potential of the sample to a scattered optical field in a sensor plane of the light-detecting elements and an incident optical field from the light source which is superposed with the scattered optical field to form the predicted interference patterns, wherein the iterative algorithm reduces a sum of the data fidelity term and a non-differentiable regularization term and wherein the iterative algorithm includes a forward-backward splitting method alternating between forward gradient descent on the data fidelity term and backward gradient descent on the regularization term; and a processing unit configured to perform digital holographic reconstruction on the acquired interference patterns based on the application of the iterative algorithm.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The above, as well as additional objects, features and advantages of the present inventive concept, will be better understood through the following illustrative and non-limiting detailed description, with reference to the appended drawings. In the drawings like reference numerals will be used for like elements unless stated otherwise.

(2) FIG. 1 is a flow chart of a method according to an embodiment.

(3) FIG. 2 is a schematic view illustrating acquisition of an interference pattern.

(4) FIG. 3 is a schematic view of a device according to an embodiment.

DETAILED DESCRIPTION

(5) Referring now to FIG. 1, a method 100 for three-dimensional imaging of a sample according to an embodiment will be described. The method 100 is configured to provide a reconstruction of a three-dimensional scattering potential of the sample.

(6) According to the method, a set of interference patterns are received 102 as input for three-dimensional reconstruction. The interference patterns may be formed by scattered light from the sample and non-scattered light from a light source, wherein the interference patterns are acquired using different angles between the sample and the light source based on a set of interference patterns. The non-scattered light from the light source may be passed along a common optical path with the light being scattered by the sample. Thus, the interference pattern may be formed within a wavefront passing the sample in a so-called in-line holography set-up.

(7) The interference patterns may be detected using one or more image sensors. An image sensor may comprise a plurality of light-detecting elements in an array. The plurality of light-detecting elements may detect incident light so as that the array of light-detecting elements may detect the interference pattern.

(8) Thanks to the use of different angles between the light source and the sample in the acquisition of the interference patterns, the combined information in the interference patterns may be used for reconstructing the three-dimensional scattering potential of the sample.

(9) The lack of phase information in the acquired interference patterns and using a limited number of illumination angles implies that the reconstructed image may comprise considerable artifacts and shape distortion.

(10) The method 100 is configured to perform digital holographic reconstruction in order to provide a three-dimensional scattering potential of the sample while avoiding or at least reducing artifacts and shape distortion in the reconstruction. The method may thus comprise formulating an inverse problem with appropriate regularization terms, wherein the inverse problem is defined for finding the three-dimensional scattering potential of the sample that would provide the interference patterns as detected by the one or more image sensors. Thus, the method may comprise formulating 104 the inverse problem and initialize parameters.

(11) The inverse problem may be generally described as:

(12) $\begin{array}{cc}{f}^{*}=\arg \underset{f}{\min }\left\{𝒟\left(f\right)+R\left(f\right)\right\},& \#\left(1\right)\end{array}$
wherein f is the three-dimensional scattering potential of the sample, (f) is a data fidelity term which evaluates the difference between the measured interference pattern and the interference pattern calculated by a light propagation model from the three-dimensional scattering potential, and (f) is a regularization term.

(13) The goal of (Equation 1) is to find the scattering potential of the object f that minimizes the sum of the data fidelity term (f) and the regularization term (f).

(14) In an embodiment, the data fidelity term may be expressed as follows:

(15) $\begin{array}{cc}𝒟\left(f\right)=\underset{p=1}{\overset{P}{.Math.}}\frac{1}{2}{.Math..Math.S_p\left(f\right)+u_p^{\mathrm{inc}}.Math.-y_p.Math.}^2& \#\left(2\right)\end{array}$
where, with reference to FIG. 2, S.sub.p(.) is a light propagation model that maps f to a scattered optical field in a sensor plane 202, u.sub.p.sup.inc denotes a corresponding incident optical field 204 from a light source 206 which is superposed with the scattered optical field to form the interference pattern, and y.sub.p represents the measured amplitude of the field 208 (detected interference pattern by an image sensor 210 comprising a plurality of light-detecting elements 212). The subscript p indicates the p.sup.th illumination angle.

(16) The light propagation model may be linear or non-linear. The image sensor may be configured to detect a large field of view and therefore, a linear model may be chosen to save time and memory in computations.

(17) As for (f), various regularization terms can be used, which may for example depend on the specific object to be imaged. However, at least some of the regularization terms may be non-differentiable (e.g., L.sub.1 norm and total variation).

(18) According to an embodiment, the regularization terms may comprise at least one of a L1 norm, a L2 norm, total variation or bound constraint. In such case, the inverse problem may be phrased as

(19) $f^{*}=\arg \underset{f}{\min }\left[\underset{p=1}{\overset{P}{.Math.}}\frac{1}{2}{.Math..Math.S_p\left(f\right)+u_p^{\mathrm{inc}}.Math.-y_p.Math.}^2+{\mu }_{L_1}{.Math.f.Math.}_{L_1}+{\mu }_{\mathrm{TV}}{.Math.f.Math.}_{\mathrm{TV}}\right],$
where ∥f∥.sub.L1, ∥f∥.sub.TV denote the L.sub.1 norm and the total variation of f, respectively and μ.sub.L1 and μ.sub.TV are weights for the corresponding regularization terms.

(20) The L1 norm corresponds to a sparsity constraint which may be suitable to apply to sparse samples such as isolated particles. The larger the weight of this regularization term, the sparser the sample should be.

(21) The total variation imposes piece-wise constant constraint to the sample which may be suitable for samples whose refractive index does not change too significantly within a small region.

(22) The bound constraint is applicable if the range of the value of the refractive index of the sample is known so that the method may be forced to find the solution within this range.

(23) The method further comprises performing digital holographic reconstruction based on the formulated problem. In view of the regularization term being non-differentiable, the method may use an iterative algorithm to reduce or minimize the sum in (Equation 1).

(24) The method may use a forward-backward splitting (FBS) strategy to handle the non-differentiability. The FBS alternates between forward gradient descent 108 on (f) (Equation 3) and backward gradient descent 110 on (f) (Equation 4) until a stopping criterion is met. The iterative algorithm may thus also include determining 112 whether a stopping criterion is met. The stopping criterion may be that a minimum has been found, that a maximum number of iterations has been reached or that the sum of the data fidelity term and the regularization term has been reduced below a set threshold.

(25) The iterative algorithm by FBS may be defined as

(26) $\begin{array}{cc}{\hat{f}}^{k+1}=f^k-{\tau }^k\nabla 𝒟\left(f^k\right),& \#\left(3\right)\\ f^{k+1}=pro{x}_{ℛ}\left({\hat{f}}^{k+1},{\tau }^k\right)=\arg \underset{f}{\min }\left\{{\tau }^kℛ\left(f\right)+\frac{1}{2}{.Math.f-{\hat{f}}^{k+1}.Math.}^2\right\},& \#\left(4\right)\end{array}$
where f.sup.k is the three-dimensional scattering potential at the k.sup.th iteration, τ.sup.k is a step size that decides how far the algorithm should move the solution in the direction of the gradient at the k.sup.th iteration and prox denotes the proximal operator whose objective is to find a scattering potential close to both {circumflex over (f)}.sup.k+1 and the minimizer of the regularization term R.

(27) The minimizing of the sum in (Equation 1) may thus be reduced to two subproblems related to phase retrieval in (Equation 3) and regularizations in (Equation 4).

(28) Solving of (Equation 3) may be performed by calculating the gradient of the data fidelity term with the help of Wirtinger derivatives:

(29) $\begin{array}{cc}\nabla 𝒟\left(f^k\right)=\underset{p=1}{\overset{P}{.Math.}}{S}_p^H\left\{S_p\left(f^k\right)+u_p^{\mathrm{inc}}-y_p\odot \frac{S_p\left(f^k\right)+u_p^{\mathrm{inc}}}{.Math.S_p\left(f^k\right)+u_p^{inc}.Math.}\right\}& \#\left(5\right)\end{array}$
where S.sub.p.sup.H is the Hermitian conjugate of S.sub.p and A⊙B denotes the Hadamard product of A and B.

(30) Next, (Equation 4) may be solved by applying a primal-dual splitting method, as further explained in (Equation 11)-(Equation 14) below. The benefit of this approach is a composite of common regularization terms including non-differentiable ones may be used. Therefore, using primal-dual splitting is very useful in real-world applications because biological samples with different optical properties often require different a priori constraints in the regularization.

(31) Once the stopping criterion is met, the method may output 114 a result in the form of a three-dimensional scattering potential of the sample.

(32) An accelerated variant of FBS may be used, namely a fast iterative shrinkage-thresholding algorithm (FISTA) to reduce reconstruction time. FISTA may be defined as

(33) $\begin{array}{cc}{\hat{f}}^k=f^k-\gamma \nabla 𝒟\left(f^k\right),& \#\left(6\right)\\ x^k=pro{x}_{ℛ}\left({\hat{f}}^k,\gamma \right)=\arg \underset{f}{\min }\left\{\gamma ℛ\left(f\right)+\frac{1}{2}{.Math.f-{\hat{f}}^k.Math.}^2\right\},& \#\left(7\right)\\ {\alpha }^{k+1}=\frac{1+\sqrt{1+4{\left({\alpha }^k\right)}^2}}{2},& \#\left(8\right)\\ f^{k+1}=x^k+\frac{{\alpha }^k-1}{{\alpha }^{k+1}}(x^k-x^{k-1}),& \#\left(9\right)\end{array}$
wherein γ∈(0, 1/Lip(∇)), where Lip(∇) denotes the Lipschitz constant of ∇D and f.sup.1=x.sup.0∈.sup.L.sup.x.sup.×L.sup.y.sup.×L.sup.z and L.sub.x, L.sub.y, L.sub.z represent a size of a region of interest in the sample. This iteration may be performed for a set number of iterations.

(34) In an embodiment, (Equation 7) and (Equation 6) may be formulated as:

(35) $\begin{array}{cc}{x}^{*}=\arg \underset{x}{\min }\left\{\mathrm{\gamma ℛ}\left(x\right)+\frac{1}{2}{.Math.x-\left(f^k-\gamma \nabla 𝒟\left(f^k\right)\right).Math.}^2\right\},& \#\left(10\right)\end{array}$
wherein the gradient of the data fidelity term may be solved using Wirtinger derivatives as defined in (Equation 5). Further, (Equation 10) may be solved by using a function P(x), which may be set as:
P(x)=1/2∥x−(f.sup.k−γ∇(f.sup.k))∥.sup.2  #(11)
Further, γ(x) may be rewritten as G(x)+Σ.sub.i.sup.MF.sub.i(K.sub.ix), where K.sub.i is an appropriate linear operator. Then, (Equation 10) may further be rewritten as

(36) $\begin{array}{cc}{x}^{*}=\arg \underset{x}{\min }\left\{P\left(x\right)+G\left(x\right)+{.Math.}_q^M{F}_q\left(K_{q}x\right)\right\}.& \#\left(12\right)\end{array}$

(37) (Equation 12) may further be solved by an iteration in a set number of iterations i, by calculating
x.sup.i+1=prox.sub.τG(x.sup.i−τ(∇P(x.sup.i)+Σ.sub.q.sup.MK*.sub.qy.sub.q.sup.i))  #(13)
where K.sub.q* is the Hermitian conjugate of K.sub.q, and calculating
y.sub.q.sup.i+1=prox.sub.σF.sub.q*(y.sub.q.sup.i+σK.sub.q(2x.sup.i+1−x.sup.i))  #(14)
where F*.sub.q is the Fenchel conjugate of F.sub.q and τ and σ are parameters which are set in initialization to τ>0 and σ>0.

(38) Thanks to the method according to the present inventive concept, the digital holographic reconstruction may determine a three-dimensional scattering potential of a sample based on only a few acquired interference patterns. In an embodiment, only four interference patterns may be used for determination of the three-dimensional scattering potential of the sample.

(39) Since the digital holographic reconstruction may be based on very few acquired interference patterns, the interference patterns may be acquired in a very short time, enabling very fast three-dimensional imaging. This implies that a time between two different three-dimensional images of the sample in a sequence of three-dimensional images may be very short. Hence, the method allows imaging of quickly changing processes, such as imaging fast biological processes.

(40) Also or alternatively, a high throughput of the three-dimensional imaging may be provided since a short time is required between subsequent three-dimensional images. Further, the imaging may be performed with a large field of view, which is also beneficial for providing a high throughput. For example, digital holographic imaging devices may have a larger field of view than conventional imaging/microscopy devices. A large field of view may mean that a large lateral extension of the sample may be simultaneously imaged.

(41) The method may be used for determining the three-dimensional scattering potential of a biological sample. In particular, thanks to the method enabling fast determination of the three-dimensional scattering potential, the method may be used for imaging fast biological processes, e.g. providing an acquisition rate of more than one image per second.

(42) Referring now to FIG. 3, a device 300 for three-dimensional imaging of a sample 302 according to an embodiment is described. The device 300 is able to acquire interference patterns for enabling reconstruction of a three-dimensional scattering potential of a sample 302.

(43) The device 300 may comprise at least one light source 306a-d. The at least one light source 306a-d may be configured to illuminate a sample 302 from a plurality of different angles;

(44) As shown in FIG. 3, fixed light sources 306a-d may be used. In such case, the relation between the light source 306a-d and the sample 302 may be accurately defined on manufacture or set-up of the device 300. However, it should further be realized that one or more movable light sources may alternatively or additionally be used, whereby a flexibility of the angle between the light source and the sample 302 may be used.

(45) The light sources 306a-d may be configured to illuminate the sample 302 from angles that are far apart. This implies that the interference patterns may acquire different information in relation to the sample 302 to facilitate accurately reconstructing the three-dimensional scattering potential of the sample 302.

(46) The light sources 306a-d may be evenly distributed over angles in relation to the sample 302. The relation between the light source 306a-d and the sample 302 may be defined by a polar angle and an azimuthal angle. In FIG. 3, the light sources 306a-d are shown including a centrally placed light source 306a that illuminates the sample 302 from straight above, i.e. having a 0° polar angle and a 0° azimuthal angle to the sample 302. The remaining three light sources 306b-d are arranged with a 45° polar angle in relation to the sample 302 and with different azimuthal angles of 0°, 120°, and 240°, respectively.

(47) The light sources 306a-d may illuminate the sample 302 using an at least partially coherent light source 306a-d. In this regard, the light source 306a-d may be laser sources. According to an alternative, the light sources 306a-d may be light emitting diodes, which may be combined with a pinhole in order to generate at least partially coherent light.

(48) The device 300 may further comprise at least one image sensor 310 for detecting the interference pattern that may be formed by scattered light from the sample 302 and non-scattered light from the light source 306a-d. The image sensor 310 may comprise an array of light-detecting elements for detecting intensity of incident light on the respective light-detecting elements.

(49) The image sensor 310 may for example be a charge-coupled device (CCD) or a complementary metal-oxide-semiconductor (CMOS) image sensor.

(50) The light sources 306a-d and the image sensor 310 may be used for sequentially detecting the interference patterns. However, according to an alternative, the interference patterns may be simultaneously detected. In such case, the interference patterns may be acquired by different regions of the image sensor 310 or even by different image sensors such that the interference patterns may be separately determined. However, the interference patterns may overlap and still allow the digital holographic reconstruction to determine the three-dimensional scattering potential of the sample 302.

(51) The device 300 may further comprise a processing unit 330. The processing unit 330 may be configured to perform the digital holographic reconstruction as described above with reference to FIG. 1.

(52) The processing unit 330 may be arranged in a common housing with the light sources 306 a-d and the image sensor 310. Hence, the device 300 may be formed as a self-contained unit, which is able to acquire interference patterns of a sample 302 and determine the digital holographic reconstruction.

(53) Alternatively, the processing unit may be distributed in two or more units which perform different parts of calculations. A part of the distributed processing may then be performed within the housing of the device 300 before processed information is transmitted to an external unit. As yet another alternative, the processing unit may be external to the housing of the device 300, such that the acquired interference patterns are communicated to an external unit for digital holographic reconstruction. In an embodiment, the processing for performing digital holographic reconstruction may even be done “in the cloud”.

(54) The processing unit 330 may for instance comprise a general-purpose processing unit, which may be provided with instructions for performing digital holographic reconstruction. Alternatively, the processing unit 330 may be implemented as firmware arranged e.g. in an embedded system, or as a specifically designed processing unit, such as an Application-Specific Integrated Circuit (ASIC) or a Field-Programmable Gate Array (FPGA), which may be configured to implement functionality of the processing unit 330.

(55) The instructions for performing digital holographic reconstruction may be provided in form of software, which may be separately delivered from the device 300, and which may e.g. be loaded to a processing unit 330 of an existing device 300 for improving functionality of the device 300.

(56) In the above the inventive concept has mainly been described with reference to a limited number of examples. However, as is readily appreciated by a person skilled in the art, other examples than the ones disclosed above are equally possible within the scope of the inventive concept, as defined by the appended claims.