Method of reconstructing a spatial distribution of a characteristic of an object.

Abstract

Method for reconstructing a spatial distribution of a characteristic () in an object, comprising: acquisition of measurements () by a sensor (15), each measurement being able to be estimated a linear operator ((), ), applied to the spatial distribution of the characteristic (), forming a direct model; using a processing unit (20), reconstruction of the spatial distribution of the characteristic of the object, by iterative minimization of an error, each iteration comprising an update of the spatial distribution of the characteristic of the object; the method being characterized in that during step b), the minimized error comprises: a data attachment component () comprising a difference between the acquired measurements and the measurements estimated by the direct model; a regularization component (), comprising a sum of a norm of a spatial gradient of the characteristic, determined in different coordinates in the object.

Description

Title of the invention: Method for reconstructing a spatial distribution of a characteristic of an object. Technical field

[0001] The technical field of the invention is the reconstruction of a characteristic of an object from non-destructive measurements carried out facing this object. PREVIOUS ART

[0002] Fluorescence imaging is a technique for locating fluorescent markers in a human or animal body. One of the main applications is the localization of fluorescent markers, or fluorophores, the latter targeting cells of interest, for example cancer cells. The protocol consists of injecting these markers into the body before a fluorescence imaging examination, during which fluorescence images are formed. A reconstruction algorithm is then implemented, so as to determine the position of the fluorophores in the object examined. This algorithm makes it possible to obtain a fluorescence map, each term of which corresponds to a fluorescence intensity at different positions in the object examined. A particularity of the reconstruction is that the fluorescence map includes positive or negative terms.

[0003] Reconstruction algorithms have been described, for example in WO2010103026.

[0004] In fields other than fluorescence imaging, other reconstruction algorithms under positivity constraint have been described in the publications Daube-Witherspoon M “An iterative image space reconstruction algorithm suitable for volume ect”. IEEE transactions on medical imaging 5(2):61-66, 1986, or in D. Seung and L Lee. “Algorithms for non-negative matrix factorization. Advances in neural information processing Systems”, 13:556-562, 2001.

[0005] The inventors propose an algorithm for reconstructing a spatial distribution of a characteristic of an object, making it possible to obtain a result close to the physical reality in the object, while being low-memory consuming during its implementation. The algorithm can be applied to fluorescence imaging, but also to other analysis methods such as image deconvolution. Statement of the invention

[0006] A first object of the invention is a method for reconstructing a spatial distribution of a characteristic in an object, the object being discretized according to different spatial coordinates defined in a reference frame, the reference frame being defined according to at least one axis, the method comprising: a. acquisition of measurements by a sensor placed facing the object, each measurement can be estimated by a linear operator, applied to the spatial distribution of the characteristic, forming a direct model; b. using a processing unit, reconstruction of the spatial distribution of the characteristic of the object, in each spatial coordinate, by iterative minimization of an error, each iteration being assigned a rank, each iteration comprising an update of the spatial distribution of the characteristic of the object, the first iteration being implemented from an initialized spatial distribution,

[0007] the method being characterized in that, during step b), the minimized error comprises: - a data attachment component, comprising a gap between the acquired measurements and the measurements estimated by the direct model; - a regularization component, comprising a sum of a norm of a spatial gradient of the characteristic, determined in different coordinates in the object;

[0008] and in that each iteration comprises an update of a previous spatial distribution, corresponding either to the initial spatial distribution or to a spatial distribution resulting from a previous iteration, the update comprising a product for each spatial coordinate, of - the previous spatial distribution; - the adjoint operator of the direct model applied to the acquired measurements; - for at least one axis of the reference frame, a product, term by term: • of said previous spatial distribution translated, along the axis of the reference frame, by at least one unit, in an increasing direction; • of said previous spatial distribution translated, along the axis of the reference frame, by at least one unit, in a decreasing direction.

[0009] The update may include, in each spatial coordinate, a multiplication of the previous spatial distribution by a linear combination of products, each product being associated with an axis of the reference frame, and including a term-by-term multiplication: • of said previous spatial distribution translated, along the axis of the reference frame, by at least one unit, in an increasing direction; • of said previous spatial distribution translated, along the axis of the reference frame, by at least one unit, in a decreasing direction.

[0010] Each term-by-term multiplication includes the inverse of a norm of the spatial gradient of the prior distribution for said coordinate.

[0011] The update may include

[0012] Or : ( / (^9=^(^(^-9)+4 / 7^0^ + ^5^ ; Q (F M ) 4rOL A ^ / -'+E A y^ ] ; M corresponds to the measurements taken; FkA is the prior spatial distribution; F is a spatial distribution of the inverse of a norm of a spatial gradient of the prior spatial distribution, in each coordinate of the object; A corresponds to each axis of the reference frame; 5«_A is an operator shifting one or more units along the A axis, in the decreasing direction; 5A-„ is a shift operator of one or more units along the A axis, in the increasing direction; n corresponds to the number of axes considered; 2 is a positive real; H is the direct model operator; H' is the adjoint operator of the direct model; O denotes the Hadamard product.

[0013] Step a) may comprise forming at least one image of the object, each image forming a spatial distribution of measurements acquired by the sensor.

[0014] According to one possibility: the sensor is an image sensor, configured to form an image of the object at an emission wavelength; the object is likely to comprise a fluorophore, emitting light, at the emission wavelength, under the effect of illumination at an excitation wavelength; in step a), at least one image of the object is acquired when the object is illuminated at the excitation wavelength, different from the emission wavelength, each measurement being representative of a quantity of light emitted by the fluorophore in different coordinates in the object.

[0015] According to one possibility: the sensor is an image sensor, configured to form an image of the object at an emission wavelength; the object is likely to absorb light, at the emission wavelength, in step a), at least one image of the object is acquired when the object is illuminated at the emission wavelength, each measurement being representative of a quantity of light absorbed at different coordinates in the object.

[0016] According to one possibility: - the sensor is a sensor intended to detect an ionizing ray, of type X or gamma; - in step a), the object is irradiated by an X or gamma irradiation beam, each measurement being representative of an absorption of the beam by the object.

[0017] The direct model may include a Radon transform.

[0018] The characteristic may be an emission or reflection or backscattering characteristic of an electromagnetic wave or an acoustic wave.

[0019] A second object of the invention is a processing unit, configured to implement step b) of a method according to the first object of the invention from measurements carried out by a sensor placed facing an object, so as to obtain a spatial distribution of a characteristic of the object.

[0020] A third object of the invention is a measurement system, configured to reconstruct a spatial distribution of a characteristic of an object, the object being discretized according to different spatial coordinates associated with a reference frame, and comprising: - a sensor, arranged facing the object, and configured to acquire measurements, each measurement being able to be estimated by a linear operator, applied to the spatial distribution of the characteristic, forming a direct model; - a processing unit, configured to implement step b) of a method according to the first subject of the invention.

[0021] The invention will be better understood upon reading the description of the exemplary embodiments presented in the remainder of the description, in conjunction with the figures listed below. FIGURES

[0022] [Fig.l] shows a diagram of a device enabling implementation of the invention, so as to reconstruct a fluorescence map of an object.

[0023] [Fig.2] illustrates an impulse response of the sensor described in connection with the device of [Fig.l].

[0024] [Fig.3] shows schematically an acquisition of images by the device described in connection with [Fig.l],

[0025] [Fig.4] illustrates the principle of a reduction in the increase.

[0026] [Fig.5] shows schematically the main steps of a method according to the invention.

[0027] Figures 6A to 6D are images resulting from a reconstruction of a fluorescence map of a mouse embryo. PRESENTATION OF SPECIAL EMBODIMENTS

[0028] [Fig.l] represents a device allowing an implementation of the invention. In this example, the device is configured to acquire images of a sample, so as to locate fluorescence light sources in the object. It This is one of the analysis modalities that can implement the method. Other possible modalities are described below, for example an X-ray or an X-ray tomography, or a reconstruction of light-absorbing areas in a sample.

[0029] The device comprises a light source 11, configured to illuminate a sample 10. The sample is a solid volume that is to be analyzed. Under the effect of illumination by the light source, the sample emits emission light. The light source emits illumination light 12 in an illumination spectral band. A portion of the light emitted by the light source propagates through the sample. Under the effect of illumination, the sample emits fluorescence light 13 in an emission spectral band.

[0030] The light source may be a laser source or a light-emitting diode. The light source may be fiber-based, with the light emitted by the light source being guided to the sample with an optical fiber.

[0031] The sample contains fluorophores, emitting fluorescence light, in the emission spectral band, when illuminated in an excitation spectral band. The aim is to determine the position of the fluorophores in the sample. This involves reconstructing a three-dimensional spatial distribution of light emission inside the sample. This makes it possible to identify areas of the sample with a high concentration of fluorophores.

[0032] The sample is for example a biological tissue, which one wishes to analyze in order to identify possible singularities. The objective is to identify local concentrations of fluorophores in the sample, the latter being able to be used to help determine a pathological state.

[0033] The sample 10 is discretized into voxels, called “object voxels”. Each object voxel corresponds to an elementary volume of the sample. It may for example be a volume of 100 nm x 100 nm x 100 nm to 10 pm x 10 pm x 10 pm. In the following, each object voxel is identified by a three-dimensional spatial coordinate r defined in an XYZ frame. N corresponds to the total number of object voxels discretizing the object.

[0034] The device comprises a pixelated image sensor 15, configured to form an image of the sample. The image sensor 15 is coupled to an optical system 16, the latter making it possible to conjugate an object plane, called the focal plane, with the image sensor. The image sensor and the optical system are aligned along an optical axis A. The assembly formed by the image sensor and the optical system is configured so that the focal plane can be translated parallel to the optical axis A. The optical system can combine an objective and a tube lens. The image sensor can be a CMOS type sensor. The focal plane can thus be translated according to different depths inside the sample.

[0035] Alternatively, the image sensor is associated with a confocal diaphragm, the latter allowing the successive observation of different slices of the sample.

[0036] Generally, the image sensor is configured to acquire images in different planes, extending to different depths in the sample. The sample is delimited by a surface S, forming an interface between the sample and an ambient medium in which the image sensor extends. The ambient medium is generally air. A depth in the sample corresponds to a distance from the surface S, parallel to the optical axis A. The image sensor can be coupled to a spectral filter, for example of the bandpass type, so as to detect a light wave in the emission spectral band.

[0037] The device comprises a processing unit 20, arranged to process the images formed by the image sensor, so as to estimate the three-dimensional spatial distribution of fluorescence light emission in the sample.

[0038] The spatial distribution of fluorescence emission F can be linked to the measurement model based on an additive noise model:

[0039] M = P*F + B(Ï)

[0040] where: • F corresponds to a 3D tensor in which each term F(r) is the fluorescence intensity in an object voxel with three-dimensional coordinate r: this corresponds to what we are trying to estimate. • M is the set of fluorescence images formed by the image sensor: this corresponds to the set of measurements forming a 3D tensor, at the same discretization step as F. The measurements M(r') are defined in image voxels r'. The number of image voxels is N' = N; • P is an operator representing the instrument response: This is a 3D tensor, representing an impulse response of the measuring instrument (PSF: Point Spread Function). P is determined using the same discretization step as the tensors M and F. The impulse response P is established by modeling and can be recalibrated by experimental tests. The operator is of dimension (N, N'), with, in this example N' = N. P can be obtained by a ray propagation method, described in Van Toey, J "Beam-Propagation method: analysis and assessment. [Fig.2] represents the PSF of the image sensor 15 in an XZ plane. The value of the PSF, in each coordinate, corresponds to the gray level. • * denotes the discrete 3D convolution product operator. • B denotes the noise term. It is a 3-D tensor of the same dimension as M..

[0041] In this application, both the response of the instrument P and the intensity of fluo rescence are positive. The noise B can be positive and negative but it is small so that the measurement M is also positive.

[0042] The image voxels are distributed according to a regular sampling step Ax, Ay, Az respectively along axes X, Y and Z, Z corresponding to the optical axis. The object voxels are distributed according to a discretization step Ax, Ay, nAz, where n corresponds to the refractive index of the object, which is assumed to be homogeneous.

[0043] For example, the sampling of the object and the measurements may for example comprise 2000 x 2000 points along the X and Y axes and 1000 points along the Z axis.

[0044] One of the steps of the reconstruction is to obtain an estimate of the measurements M = P*F approaching M. It is a question of determining an error term s, between M and M. More generally, — prfp\, H corresponding to the model operator direct, according to which the measurements are estimated from the spatial distribution of the characteristic in the object.

[0045] A first option is to minimize an error corresponding to a com comparison of m and M: _ y -M || corresponds here to a data attachment term. t'D(F) is a scalar.

[0046] Preferably, the inventors consider that it is preferable for the error to combine a data attachment term and a regularization term. The data attachment term is a comparison of the measurements M and their estimation p^. The regularization term aims to obtain a representation of the spatial distribution of fluorescence F representative of reality. In particular, when the regularization term takes into account a spatial gradient of the fluorescence intensity, in each pixel, such that

[0047] h dF^r) '2 l <W) x2 W +\~J +\~dr)

[0048] The regularization term can be a total variation, such that 100491 F) = TV (F ) = <3>

[0050] where denotes the spatial position of the nth voxel. eR{F^ is a scalar.

[0051] The regularization term corresponds to a sum of the L2 norm of the spatial gradient of the spatial distribution of fluorescence defined in each object pixel.

[0052] Expression (3) corresponds to a regularization called “total variation”. Other regularization terms are possible, in particular of the form:

[0053] / n———:^772 \ “with o< «^2(3'). y Va. r ■'■'«m dx J Sy / ùz / I

[0054] In the following example, « = 1

[0055]

[0056]

[0057]

[0058]

[0059]

[0060]

[0061]

[0062]

[0063]

[0064]

[0065]

[0066]

[0067]

[0068]

[0069]

[0070]

[0071]

[0072]

[0073]

[0074]

[0075] The error to be minimized can be such that: e(F) =8D(F) +àcr(F) = || MM(F) || TV(F) (4) 2 is a weighting factor. 2 is a positive real. Expression (4) includes a data attachment term eD(F) = || MM(F) || 2 ct un dc æêulaHsad^^ = TV(F ), the balance between the two terms (data attachment and regularization) being measured by the weighting factor 2. Taking into account the regularization term makes it possible to denoise the fluorescence map F tends to create groups of pixels in which F(f) is homogeneous, the borders of the groups of pixels being generally sharp. Minimizing data attachment. efF)^ It ?*F - M || 2 (5) The canonical form of (5) corresponds to: ^(F) =4 II II 2(6) where R(P)f = vec(P*F) O) with the vectorization operator, which forms a vector where each term is a voxel of P*F / is a vector of dimension N, each term f(r) of which corresponds to an intensity fluorescence, assumed to be positive. R(P) is a matrix, of dimension (N,N), representing the 3-D convolution by the impulse response P. R(P)f is a vector, of dimension N, each term of which corresponds to the convolution of F by P for each image voxel. eDiF) = eD(F) = | fTR (P)TR(P)fj fR ( P ) Tm +1 m^m - jm^R(P) f m is a vector in which each term m(r) corresponds to a measure. We set = ^)^(^) (10) Qd is a symmetric semi-positive definite (SSPD) matrix, of dimension N x N. Generally speaking: a square matrix A is symmetric if Ar = A- a square matrix A is semi-positive definite if VU 0, UrAu > 0M being a vector of the same size as the matrix. Given (9), A / L if QJ+s it is a vector of dimension N.

[0076]

[0077]

[0078]

[0079]

[0080]

[0081]

[0082]

[0083]

[0084]

[0085]

[0086]

[0087]

[0088]

[0089]

[0090]

[0091]

[0092]

[0093]

[0094] c= -R(P)Tm^c^ ^mTR(P)f = yTR(P)Tm cte is a constant which does not depend on / : (12') The constant not involved in the minimization of it can be neglected thereafter. The QD matrix can be decomposed into two symmetric matrices, one of positive elements, and one of QD, each representing the positive and negative terms of the QD matrix respectively. Q d = QJ-Q d ^ Similarly, vector 6 can be decomposed into two vectors of positive elements c'+ etc, each representing the respectively positive and negative terms of c. c = c+-c(14) The minimization of the data attachment error is performed iteratively. Each iteration is assigned a rank k. During each iteration, the vector / is updated according to the expression: fk^_ fk-lÇ\ QDfk'l+c (15) JJ The symbol 0 corresponds to the Hadamard product. In this example, which concerns a deconvolution of a fluorescence measurement, Qd is a matrix of positive elements, because R(P) only contains positive terms, so D=() and 2D+= Qd Similarly, c is a vector containing only negative terms (cf. (12)), so c' = - c Expression (15) becomes: Expression (16) can be expressed as: pk pk-1Q PpM Q 7 ) Pr corresponds to the returned tensor P, according to a central symmetry. It is an adjoint operator of the previously defined direct model, applied to the measures M. More generally, if H denotes the operator corresponding to the direct model, H denotes the adjoint operator. In this example, H( F) = P*F- The adjoint operator H is such that H (M} = PR*M or H (P*FkA ) = Pp*P*Fkl- Thus, (17) can be expressed by: pk pk'^QH (17 )

[0095]

[0096]

[0097]

[0098]

[0099]

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106]

[0107]

[0108]

[0109] Regularization term. According to (3), the regularization term SR can be such that: ^F^TViF ) = || GIF) || F£n\G(F(r^^ (20) Dx is an "X-axis derivative" operator; DY is a "Y-axis derivative" operator; Dz is a “Z-axis derivative” operator The term | G {F ( r ) ) | ^ = corresponds to the L2 norm of the spatial gradient G(F( r )) of the spatial distribution F at the object voxel r. We can show that whatever F(f), i lGW<2 (21) (r) ) l2 “ 2 |G(F0)|, ^2^(^())12 Expression (21) reflects the fact that | G (F ( r ) ) | is bounded above by a parabola, tangent at a point G ( Fo ) • i' avec G0- G(fS 2—i^i— + 2iGol 07 We apply the MM method for the upper bound of a lower bound: when a real-valued function (in this case \G(F(r) ) | ), is upper bounded and tangent at Fo by a function M(F(r) ), (in this case J |G(F(r)) j ), if a point F(r) 2 |Goj decreases the value of M(F(r) ) then this point then decreases the value of |G(F(r) ) See figure 3. More generally, if a point with abscissa u minimizes M(u), this same point minimizes |G(m) I- The fluorescence map F is updated iteratively at each iteration being assigned a rank k. K is a non-zero positive integer. In the first iteration, k = 1. The first iteration is performed from an initialized fluorescence map F°. Subsequently, we consider that j Go | corresponds to the value of the gradient G ( FA'! ( r ) ) during iteration k - 1. ■me ) The term 1£ | G (F*^ r„ ) ) | being constant, it does not contribute to the minimization. Generally speaking, = aT£a (24), where - a is a vector comprising the terms a(rn); - B is a diagonal matrix formed from the terms ... / 3( / ^)-

[0110] The term j „ ) can be expressed in canonical form according to the expression: [Un yTRr(Dx)rR(Dx)f+yTRT(Dr)r^ (25)where - R(Dx ), R(Dy) and 7?(Dz ) are matrix representations of the respective convolutions by Dx, DY and Dz. - r is a diagonal matrix of dimension N, / / of elements ) l2

[0112] Which amounts to (27)

[0114] Each term of the diagonal of L is the inverse of an L1 norm of a spatial gradient of at2 each coordinate rn. Global Criterion

[0115] Taking into account (4), (9) and (25), we can write:

[0116] = + (30)

[0117] We now pose

[0118] R(p)TR(p) + a[Rt(Dx) FR(Dx) + RT(DY) VR(DY) + Rr(Dz) FR(DZ)] (31) [0H9] 2= Qd+ ^[Q^Qy + Q^ <32)

[0120] RT(DX) rR(Dx) (33)

[0121] Q y ^R t {D y ) VR(D y ) (34)

[0122] Qz~ Rr(Dz) rR(Dz) (35)

[0123] Now, in the case where we choose the forward derivative,

[0124] R(DX) = '1 -1 0 1 0 ■■ -I ■■ 1 -1 '■ 0 1 . 36 and R(dx)t^= 1 -1 0 0 1 . ■ -1 '■ 0 0 1 0 -J 1 1 (37)

[0125] The matrices R(DX) and R(Dx) 7 allow to obtain, for each voxel, a rate of variation along the X axis. These are matrix operators allowing to obtain, for each term of the fluorescence map, expressed here in a vector form, a difference between two terms along the X axis. More precisely, each term of the fluorescence map is assigned a coordinate along the X axis. The matrices R(Dx) and rÇ]) )T allow to make a difference between two terms of the fluorescence map adjacent along the X axis or spaced by an increment along the X axis.

[0126] Similarly, the matrices R(DY) and R(E)Y^T make it possible to obtain, for each voxel,

[0127] a rate of variation along the Y axis. This is a matrix operator used to obtain, for each term of the fluorescence map, expressed here in vector form, a difference between two terms along the Y axis. Each term of the fluorescence map is assigned a coordinate along the Y axis. The matrices R(DY) and r(Dy^t allow a difference to be made between two terms of the fluorescence map that are adjacent along the Y axis or spaced by an increment along the Y axis. Similarly, the matrices R[D7) and r^d^t allow us to obtain, for each voxel, a rate of variation along the Z axis. This is a matrix operator allowing to obtain, for each term of the fluorescence map, expressed here in a vector form, a difference between two terms along the Z axis. Each term of the fluorescence map is assigned a coordinate along the Z axis. The matrices R{D7) and R^D^7 allow to make a difference between two terms of the fluorescence map adjacent along the Z axis or spaced one increment apart along the Z axis.

[0128]

[0129] We set j+ = '10- * 1 0 ■ 0 1 ■ 0 . ■ '■ 0 ■■ 10 '' 0 f 9 38 / R^ ,R(DX '>xf = 0 c 0 0 0 ■. 0 ■. • 1 0 1 o°J (41) R(DAy+ = 0 1 ■ 0 . 0 ■■ JÛ (40 ) 1 ( 1 . ■ o ■■ 0° ■■ 0 *. ■■■ 10]

[0130] r( f) j+ and RÇ[)}T+ include the positive terms of R(Dx) and r(r res-

[0131]

[0132]

[0133]

[0134] respectively. R(Dx) and R^JJX^T~ include the opposites of the negative terms of R() and R(Dx)T respectively. 0 " a: R(D x )=R(D A .) + -R(D x )' Similarly; we can define: R(Dy ) + and R( Dy ) r+' 9 which include the positive terms of R(DY ) and r(£)^ ) T

[0135]

[0136]

[0137]

[0138]

[0139]

[0140]

[0141]

[0142]

[0143]

[0144]

[0145]

[0146]

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153]

[0154]

[0155]

[0156]

[0157]

[0158]

[0159] respectively R( Dy ) and 2? ( Dj ) T ' 9ui comprise the opposites of the negative terms of R (DY ) and R(DY )T respectively; R(DZ)+e^ R( Dzf+- 9^ include the positive terms of R(D:) and D ) respectively; R( D7 ) and Z)z)7' 9^ comprise the opposites of the negative terms of R^DZ) and R(Dz)T respectively; R(DZ) = R(Dzf-R(Dz) According to (33), Qx= RT(DX) rR(Dx) (42) <4= WDxD- (KD. ¥ ))<43 ' Qx= R(DX)T+VR(DX)*+ R(Dx)tTR(Dx)'-R(Dx)T+rR(Dx) - R(Dx)TTR(Dxf (44) We ask: Qx = R(Dx) “ru(Dx) + + R(Dx)TTR(Dx)’ <45) g + is a symmetric matrix of positive elements. X Q x = R(D x ) T+ VR(D x )\R(D^^ Qx is a symmetric matrix of positive elements. q x = e / -e;< 47 ) Similarly, we pose: Qy = R(DY)T+VR(Dy) + +R(Dy)TTR(DYy (5°) Qy= R(DY)T+VR(DYy + R{Dy)TVR(DyŸ ^^+ and Qy are symmetric matrices of positive elements. Qt = Qy-Qr^ Qz= R(DzDrR(Dzy + R(DzfTR(Dzy^ Qz= R(DZ) TXR (Dz)+R(DzyYR(Dz)*^ Qz= Qz-Qzw ôz+ and Qz are symmetric matrices of positive elements. According to (31) to (35) 101601 Q= Q / qo-q- Q y -Q y + Qz-Q z<56)

[0161] Knowing that QD = 0.

[0162] Q is an SSPD matrix.

[0163] Let+= qO1(Q \ Q / + qA <57>

[0164] And 2 = 42+ Q + Q_) (58)

[0165] 2+ and Q are symmetric matrices of positive elements.

[0166] According to (30),

[0167] = + = 1fT Qf + cf f (59)

[0168] We find a formalism expressed in relation to expression (11). As indicated in (15), according to this formalism the vector / is updated according to the expression:

[0169] A ; 4-lfT) Qfk"]+c (60), knowing that c = c (61) JJ

[0170] Thus,

[0171] or -q (62) 'f Q'T'

[0172] Expression (61) can be formulated as:

[0173] k (63) FF ° cW

[0174] with

[0175] 6(^)= ^0^(^010^ / + ÿD^*rO(Drs'C-') + 00^00(0 / )0^ D / *(rO.o / )) (64)

[0176] Q:(FtA) = + oO / rGii / O') / + Dr^rODr' • / *'))+ i(i) / OoOoO*C') + o / OrOo / 'O')) (65)

[0177] F includes the terms of the diagonal of the matrix F spatialized so that at each coordinate rn, p ( - 1 / | G ( F^' ( r„ ) ) | •

[0178] g (FÆ4) and 2+( F'"') are operators applied to the fluorescence map Fk~\ respectively defined by (64) and (65).

[0179] Expression (63) corresponds to a formula for updating the fluorescence map F between two successive iterations FkA and Fk- Using a more general notation, involving the adjoint model H of the direct model, applied to the measurements,

[0180] (66) Q {F^

[0181] and

[0182] ^(^))+^0^^01^^)+ i)F(rOox*Ft-')')++ i)r'>-f(rODY,'Ft-'))+ 1(0 / 0^000^)+ o^rOoO*!^')) (67)

[0183] In this example, / / '(M) — PR*P (68) and / / P*FkA} — P^P^F^1 (69)

[0184] The notation DXR+ means DF returned, that is to say after application of a transformation, by central symmetry, of the matrix

[0185] By observing the effect of matrices (38) to (41), we see that: - The convolution of F by leaves F unchanged. - The convolution of F by Dx'+ leaves F unchanged. - The convolution of F by Dy shifts F in the negative A direction, preferably of a unit. - The convolution of F by [)SR~ shifts F in the positive A direction, preferably by one unit

[0186] In the case of a 3-dimensional fluorescence map, A is successively equal to X, Y and Z.

[0187] In the end we calculate in a very simple way

[0188] c,(^)-(^)*^+4^®^+M,<®w*']+Mr®s^l+Mr®w^]]-HM^,))+43r0^'+Mr®5^]+M^®!^]+Mr®s^*,Jl = (70)

[0189] QAF^} = + + +S^[rOFÀM] ] (71)

[0190] 5^ is a shift operator in the direction of increasing x (same for 5—y; shift towards decreasing y and: shift towards decreasing z).

[0191] is a shift operator in the direction of increasing x (same for Sy^; shift towards increasing y and; shift towards increasing z).

[0192] Whatever the axis chosen, and whatever the direction, the shift is carried out according to a number of units equal to 1 or close to 1, for example between 1 and 10 or preferably between 1 and 5.

[0193] Expressions (70) and (71) can be written:

[0194] + (72)

[0195] Q (fm) =4rOEAS^AFw + EASA^ (73)

[0196] Where A corresponds to each axis of the frame considered. One or more axes can be taken into account, most generally two (two-dimensional spatial distribution) or three axes (3D spatial distribution). In (71), n ​​corresponds to the number of axes considered.

[0197] An advantage of the update according to the update expressions (62) or (63) or (66) is that they are simple to implement, from a memory point of view. Indeed, the update is obtained by shifts of the fluorescence map PkF according to the three directions X, Y and Z.

[0198] A particularity of the invention is that taking into account the regularization component results, during the update, in a successive shift of the fluorescence map along the three axes of the reference frame, towards the positive or negative directions. These are inexpensive operations in terms of computing power.

[0199] For example, the calculation of is a simple calculation, obtained by: - shift of all the terms of FkA, associated with a coordinate (FF z) by one rank towards the negative x: we obtain the shifted fluorescence map - multiplication, term by term, of each term of S by P; - shift the set of terms of ^POS^F^1 ], associated with a co ordinate (FF *9 of a rank towards the positive x

[0200] Thus, the update is carried out: - by calculating Q&, from the direct model operator and the adjoint operator of the direct model, with = P^P, knowing that this quantity is identical for all iterations; - by calculating c = - c — PR*M , from the adjoint operator, knowing that this quantity is identical for all iterations; - and, during each iteration, by calculating Q[Fk~^ and Q^ F'^1) linked to the regularization, knowing that these components are obtained, without requiring great computing power, by successive shifts of the pkA fluorescence map along the three axes X, Y and Z, generally by one increment, either in the positive direction or in the negative direction.

[0201] The formalism described in connection with expressions (61) to (73) makes it possible to determine, iteratively, the fluorescence map using an update formula that is simple to implement and inexpensive in terms of resources, as seen previously.

[0202] [Fig.4] shows the main steps of implementing a method according to the invention, in an example corresponding to a reconstruction of a fluorescence map (spatial distribution of fluorescence) of an object.

[0203] Step 100: illumination of the object.

[0204] During step 100, the object 10 is illuminated by the light source 11 in an illumination spectral band. In the example described, the illumination spectral band is an excitation spectral band of a fluorophore potentially present in the object.

[0205] The object can be interposed between the light source 11 and the image sensor 15, which corresponds to an acquisition configuration usually referred to by the term “in transmission”. Alternatively, the acquisition can be carried out in reflection, or backscattering, the light source 11 and the image sensor 15 being arranged facing the object.

[0206] Step 110: Acquisition of images of the object

[0207] During this step, one or more elementary images . .Mi ...Mj of the object are acquired, at different depths. Figure 5 shows a diagram of an acquisition configuration according to which the focal plane of the image sensor is successively moved to different depths in the object, which corresponds to a preferred embodiment. The light source remains fixed. The image sensor acquires as many images as there are focal lengths, each focal length being associated with a depth di in the object. This gives elementary images ..M^ I corresponding to the number of elementary images. The elementary images MI are respectively associated with the depths dA....dh. .dp In the following, the elementary images are combined to form an acquired image noted M; it is a three-dimensional image of the object. The acquired image is discretized into voxels, called image voxels M(r '), each image voxel being a pixel of an elementary image Mn forming the acquired image M.

[0208] Alternatively, the image sensor has a fixed focal length and is moved relative to the object, so that the focal plane is successively moved to different depths d^....d^.. ,d1 in the object.

[0209] Alternatively, the image sensor remains fixed and a light source is used that emits light in a narrow light blade. The light source is then translated relative to the object, preferably parallel to the optical axis. Thus, during each illumination, a layer of the object parallel to the focal plane of the image sensor is illuminated, the thickness of which is less than the difference between two successive depths d,, diA. The sensor acquires an image at each position of the light source. Different images are thus successively formed, corresponding respectively to different depths of the object.

[0210] Generally speaking, step 110 aims to obtain different images ...M,...Mj representative of different depths of the object d\....di...dI. Acquiring images by having focal planes at different depths makes it possible to obtain images in which the fluorescent sources, in the object, appear more clearly. This makes it possible to obtain a more precise reconstruction. The measurements correspond to a direct model applied to the fluorescence map.

[0211] The following steps are implemented by the processing unit 20. The latter comprises a microprocessor connected to a memory. The memory comprises instructions allowing the implementation of the processing described below. The memory of the processing unit comprises the images acquired during the step 110.

[0212] Step 120: Initialization

[0213] During this step, the spatial distribution of fluorescence in the object is initialized. Each term F(r) takes an initial value determined randomly or predefined. Preferably, the initialization is carried out such that F0 = pR*M: this is the application of the adjoint operator of the direct model to the measurements. More generally, p° —

[0214] Step 130: formation of the vector f®. During this step, each term of F°(r) forms a vector f°.

[0215] Step 140: formation of the matrices Q+, Q+, and the vectors C" and c: these matrices and vectors form input data for the algorithm. They are determined as a function of the PSF, the measurements M of the sensor and the norm of the spatial gradient taken into account in the regularization term.

[0216] Step 150: updating the vector fk according to the vector fkA resulting from a previous iteration or, during the first iteration, from the initialization (step 130). The update formula is:

[0217] fk^ (62) Q t

[0218] Expression (62) corresponds to taking into account a fluorescence map forming a vector.

[0219] The vectorization phase (step 130) is optional. The fluorescence map can be expressed in the form of a matrix, 2D or 3D, the update formula being:

[0220] kn Q(F^)+Pr*MFFU

[0221] More generally, the update formula is (66).

[0222] Step 160: repeating step 150 until an iteration stopping criterion is reached. The iteration stopping criterion is a predetermined number of iterations or a comparison between two vectors fk and fk^ resulting from two successive iterations (or between two successive fluorescence maps pk and ).

[0223] The invention was implemented to reconstruct a fluorescence map of a mouse embryo at the blastocyst stage. The nucleus of each cell was labeled by fluorescence (marker: Draq5 - registered trademark) excited at 625 nm. An acquisition of 50 image planes spaced 2 pm from each other was carried out. The diameter of the embryo was 90 pm. The size of each pixel, reduced to object space, was 108 nm.

[0224] The fluorescence map of the sample was reconstructed successively without and with implementation of the regularization. Figures 6A and 6C represent views along the XY plane (parallel to the plane of the image sensor) and XZ (perpendicular to the image sensor plane), without implementing the regularization. The update formula is as explained in (16). Figures 6B and 6D represent views along the XY plane (parallel to the image sensor plane) and XZ (perpendicular to the image sensor plane), with implementation of the regularization. It can be observed that the regularization makes it possible to obtain a reconstruction more representative of reality.

[0225] An example of reconstruction of a spatial distribution of fluorescence has been described above: the characteristic of the object is an emission intensity of a fluorescence light.

[0226] The invention can however be applied to other reconstructions of spatial distributions. It can for example be a reconstruction of the attenuation of an object to ionizing radiation, in a determined energy range, when measurements by X-ray tomography are carried out around the object. In this case, the direct model implements a Radon transform, which is also a linear operator applied to the unknowns to be sought, i.e. an absorbance of voxels of the sample. In this case, the map F to be reconstructed is the absorbance map in different voxels of the sample. The measurements M correspond to the integral of the absorbances following the linear trajectories between an X-ray source and the pixels of the detector. The measurements can be estimated according to a linear prediction model, forming the direct model: where R is a Radon transform. The attachment error to the data is calculated using the Radon transform adjoint operator R'(Af), applied to the measurements.

[0227] The invention can also be applied to an estimation of the optical absorption of an object, a difference, compared to the fluorescence modality, being that the emission wavelength of the light source corresponds to the wavelength of the photons detected by the image sensor. The characteristic is then an optical absorption in different coordinates of the object.

[0228] The feature may be a part of an image, the objective being to reconstruct a sharp image from a blurred image.

[0229] Thus, the characteristic is an optical characteristic of an object, or a characteristic of attenuation or reflection or backscattering or attenuation of a wave, in the object. It may be an electromagnetic wave (light wave, X-ray or gamma ray) or an acoustic wave, the intended applications being ultrasound. More generally, the characteristic may be a response of the object to a wave to which it is exposed.

[0230] The invention can also be applied to an operation of deconvolution of an image, for example the application of a filter whose PSF is known.

[0231] In all cases, during each iteration, taking into account the regularization results in a combination of terms, each term corresponding to a shift in the spatial distribution along at least one axis of the reference frame in which the coordinates of the object are defined.

Claims

Claims

1. Method for reconstructing a spatial distribution of a characteristic (F) in an object, the object being discretized according to different spatial coordinates (r«) defined in a reference frame, the reference frame being defined according to at least one axis (X,Y,Z), the method comprising: a. acquisition of measurements (M) by a sensor (15) placed facing the object, each measurement being able to be estimated by a linear operator (H(F ), P'^F, R(F)), applied to the spatial distribution of the characteristic (F\ forming a direct model; b. using a processing unit (20), reconstruction of the spatial distribution of the characteristic of the object, in each spatial coordinate, by iterative minimization of an error, each iteration being assigned a rank (k\ each iteration comprising an update of the spatial distribution of the characteristic of the object, the first iteration being implemented from an initialized spatial distribution, the method being characterized in that during step b), the minimized error comprises: - a data attachment component (sD ( / ), e^F}) comprising a difference between the acquired measurements and the measurements estimated by the direct model; - a regularization component (gR( s^F)), comprising a sum of a norm of a spatial gradient of the characteristic, determined in different coordinates in the object; and in that each iteration comprises an update of a previous spatial distribution, (Fa*\ fkA) corresponding either to the initial spatial distribution (f°, y0) or to a spatial distribution resulting from a previous iteration, the update comprising a product for each spatial coordinate, of - the previous spatial distribution - the adjoint operator of the direct model applied to the measurements (PR*M,R\M)y, - for at least one axis of the reference frame, a product, term by term: • of said previous spatial distribution (F^'1, fk~l) translated, along the axis of the reference frame, by at least one unit, in an increasing direction; • of said previous spatial distribution (Fm, / *-i) translated, along the axis of the reference frame, by at least one unit, in a decreasing direction.

2. Method according to claim 1, in which the update comprises, in each spatial coordinate, a multiplication of the previous spatial distribution by a linear combination of products, each product being associated with an axis of the reference frame (X,Y,Z), and comprising a term-by-term multiplication: • of said previous spatial distribution translated, along the axis of the reference frame, by at least one unit, in an increasing direction; • of said previous spatial distribution translated, along the axis of the reference frame, by at least one unit, in a decreasing direction.

3. A method according to claim 1 or claim 2, wherein each term-by-term multiplication comprises the inverse of a norm of the spatial gradient ( 1 / | G ( F1^ ( r ) ) | ) of the prior distribution for said coordinate.

4. Method according to any one of the preceding claims, in which the update comprises where: - F^ is the prior spatial distribution; - O denotes the Hadamard product; = h'^f^ ))++ EAV[r'O ; Q(Fk'1) =A[rG£As^Fk-l+^^ ]; - M corresponds to the measurements carried out; - r is a spatial distribution of the inverse of a norm of a spatial gradient of the prior spatial distribution, in each coordinate of the object; - A corresponds to each axis of the reference frame; - 5^ is a shift operator of one or more units along the A axis, in the decreasing direction; - ^a— is a shift operator of one or more units along the A axis, in the increasing direction; - n corresponds to the number of axes considered; - >1 is a positive real number; - H is the operator of the direct model; - H' is the adjoint operator of the direct model.

5. Method according to any one of the preceding claims, in which step a) comprises forming at least one image of the object, each image forming a spatial distribution of measurements acquired by the sensor.

6. Method according to claim 5, in which - the sensor is an image sensor, configured to form an image of the object at an emission wavelength; - the object is capable of comprising a fluorophore, emitting light, at the emission wavelength, under the effect of illumination at an excitation wavelength; - in step a), at least one image of the object is acquired when the object is illuminated at the excitation wavelength, different from the emission wavelength, each measurement being representative of a quantity of light emitted by the fluorophore in different coordinates in the object.

7. Method according to claim 5, in which - the sensor is an image sensor, configured to form an image of the object at an emission wavelength; - the object is capable of absorbing light, at the emission wavelength, - in step a), at least one image of the object is acquired when the object is illuminated at the emission wavelength, each measurement being representative of a quantity of light absorbed in different coordinates in the object.

8. Method according to claim 5, in which - the sensor is a sensor intended to detect an ionizing ray, of the X or gamma type; - in step a), the object is irradiated by an X or gamma irradiation beam, each measurement being representative of an absorption of the beam by the object.

9. A method according to any preceding claim, wherein the characteristic is an emission or reflection or backscattering or absorption characteristic of an electromagnetic wave or an acoustic wave.

10. Processing unit, configured to implement step b) of a method according to any one of the preceding claims from measurements carried out by a sensor arranged facing an object, so as to obtain a spatial distribution of a characteristic of the object.

11. Measurement system, configured to reconstruct a spatial distribution of a characteristic of an object, the object being discretized according to different spatial coordinates associated with a reference frame, - a sensor, arranged facing the object, and configured to acquire measurements, each measurement being able to be estimated by a linear operator (H(F), P*F, R(F)), applied to the spatial distribution of the characteristic (F\ forming a direct model; - a processing unit, configured to implement step b) of a method according to any one of claims 1 to 9.

Images