# Patent

## Title
Device and method for spatial reconstructing of absorbers mapping

## Abstract

The invention concerns a method for locating at least one absorber in a diffusing medium, using at least one excitation radiation and at least one detector (Φfluo), including: a) for at least one pair (radiation source-detector), at least one excitation by the radiation source, and at least one detection of the fluorescence signal emitted by the absorber after this excitation, b) identification of meshing of the volume into mesh elements, and c) estimation of the location of the absorber in its diffusing medium, by computing a function (Pm) of at least one of three parameters.

## Background
CROSS REFERENCE TO RELATED APPLICATION
The present divisional application claims the benefit of priority under 35 U.S.C. 120 to application Ser. No. 12/889,030, filed on Sep. 23, 2010, and claims the benefit of priority under 35 U.S.C. 119 from French Application No. 09 56610, filed on Sep. 24, 2009. The entire contents of application Ser. No. 12/889,030 is hereby incorporated herein by reference.


TECHNICAL FIELD AND PRIOR ART
The invention concerns the field of fluorescence molecular imaging of biological tissues using time-resolved optical methods.
It applies in particular to optical molecular imaging of small animals and to optical molecular imaging in man (brain, breast, other organs in which fluorophores can be injected).
The application that is notably targeted is the detection of prostate cancer.
At the current time, the locating of diseased cells can be made using biopsy techniques.
However, to find the position, even approximate position, of these cells, several samples must be taken. This invasive technique is evidently difficult to implement and is time-consuming.
The problem therefore arises of determining the location, even approximate location, of fluorophores attached to a region of a medium e.g. a biological medium.
In some cases, it is even sought to locate fluorophores with a view to subsequent surgical intervention.
Optical techniques for fluorescence molecular imaging are therefore being increasingly developed at the present time through the use of specific fluorescent markers. These preferably attach themselves to the target cells of interest (e.g. cancer cells) and offer better detection contrast than non-specific markers. The objective of these techniques is to locate fluorescence spatially but also to determine the concentration thereof.
Optical tomography systems use various light sources. Existing apparatus operates in continuous mode or in frequency-mode (which uses frequency-modulated lasers) and finally in temporal mode which uses pulsed light sources, such as pulsed lasers.
Temporal data are obtained when a pulsed source is used delivering a short signal at a determined rate. In this case, the term time-resolved diffuse imaging is used. Temporal data are data which contain the most informational content on the imaged tissue, but for which reconstruction techniques are the most complex. The measurement at each acquisition point is effectively a time-dependent function called TPSF for Temporal Spread Function.
It is sought to extract simple parameters from the TPSF, whose theoretical expression is known. Resolution of the inverse problem then allows the distribution of fluorescence to be found.
Locating techniques are known from document EP 1 884 765 for example.
However, these techniques do not give satisfaction notably in so-called “back-diffusion” geometry in which sources and detectors lie on one same side of the object or of the medium being examined.
DISCLOSURE OF THE INVENTION
The inventors have found, when applying known techniques for time-resolved diffuse optical imaging, that a certain number of problems arise. Different distributions of fluorophores may effectively lead to the same measurements, a measurement being the intensity of fluorescence for example or the mean time-of-flight, notably in reflection geometry.
For any medium, knowing its optical properties, it is possible to compute photon density for a Dirac point source located at rs and at t=0s, using Green functions.
These functions Gs(r,t) can be determined using a Monte-Carlo simulation, or by solving the radiation transfer equation, or most often, by solving the diffusion equation (1) to which boundary conditions are added:




  
    
      
        
          
            
              ∂
              
                
                  G
                  s
                
                ⁡
                
                  (
                  
                    r
                    ,
                    t
                  
                  )
                
              
            
            
              c
              ⁢
              
                ∂
                t
              
            
          
          -
          
            
              ∇
              D
            
            ⁢
            
              ∇
              
                
                  G
                  s
                
                ⁡
                
                  (
                  
                    r
                    ,
                    t
                  
                  )
                
              
            
          
          +
          
            
              μ
              a
            
            ⁢
            
              
                G
                s
              
              ⁡
              
                (
                
                  r
                  ,
                  t
                
                )
              
            
          
        
        =
        
          δ
          ⁡
          
            (
            
              
                r
                -
                
                  r
                  s
                
              
              ,
              t
            
            )
          
        
      
    
    
      
        (
        1
        )
      
    
  




The Gs moments, i.e. the quantities:




  
    
      
        
          
            G
            s
            
              (
              n
              )
            
          
          ⁡
          
            (
            r
            )
          
        
        =
        
          
            
              M
              n
            
            ⁡
            
              (
              
                G
                s
              
              )
            
          
          =
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                Gs
                ⁡
                
                  (
                  
                    r
                    ,
                    t
                  
                  )
                
              
              ⁢
              
                t
                n
              
              ⁢
              
                d
                t
              
            
          
        
      
    
    
      
        (
        
          1
          ′
        
        )
      
    
  




follow the differential equations given by (2):




  
    
      
        
          
            
              -
              
                ∇
                D
              
            
            ⁢
            
              ∇
              
                
                  G
                  s
                  
                    (
                    n
                    )
                  
                
                ⁡
                
                  (
                  r
                  )
                
              
            
          
          +
          
            
              μ
              a
            
            ⁢
            
              
                G
                s
                
                  (
                  n
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
          
        
        =
        
          
            δ
            ⁡
            
              (
              
                
                  r
                  -
                  
                    r
                    s
                  
                
                ,
                n
              
              )
            
          
          +
          
            
              n
              c
            
            ⁢
            
              G
              s
              
                (
                
                  n
                  -
                  1
                
                )
              
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




Therefore, Gs(0) is the solution of the differential equation (2) with only the Dirac as source term, and Gs(n) is the solution when the source term is proportional to Gs(n-1). The sequence of Gs(n) functions is therefore found by successive steps.
For an infinite homogeneous medium, an analytical solution of (1) is known. This is:




  
    
      
        
          
            G
            s
            ∞
          
          ⁡
          
            (
            
              r
              ,
              t
            
            )
          
        
        =
        
          
            c
            
              
                (
                
                  4
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  Dct
                
                )
              
              
                3
                /
                2
              
            
          
          ⁢
          
            exp
            ⁡
            
              (
              
                
                  -
                  
                    
                      r
                      2
                    
                    
                      4
                      ⁢
                      Dct
                    
                  
                
                -
                
                  
                    μ
                    a
                  
                  ⁢
                  ct
                
              
              )
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




The moments of this function are also determined analytically;
The 2 first are:

Gs(0)(r)=exp(−k·|r−rs|)/(|r−rs|D)  (4)

Gs(1)(r)=exp(−k·|r−rs|)/(2c√{square root over (μaD)})  (5)

These formulas are of interest when simulating measurements for which it is possible to insert the sources (and detectors) in the medium (for example by means of optical fibres at the tip of a biopsy needle).
In the remainder of this application, some embodiments have recourse to moments 0 and 1 of the measurement Msd(t) prompted by the source s and measured by the detector d.
They are useful since they give the total detected energy Isd and the time-of-flight Tsd.
Indeed:




  
    
      
        
          I
          sd
        
        =
        
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  
                    M
                    sd
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
                ·
                
                  t
                  0
                
              
              ⁢
              
                d
                t
              
            
          
          =
          
            M
            sd
            
              (
              0
              )
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




and




  
    
      
        
          T
          sd
        
        =
        
          
            
              
                ∫
                0
                ∞
              
              ⁢
              
                
                  
                    M
                    sd
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
                ·
                
                  t
                  1
                
                ·
                
                  d
                  t
                
              
            
            
              
                ∫
                0
                ∞
              
              ⁢
              
                
                  
                    M
                    sd
                  
                  ⁡
                  
                    (
                    t
                    )
                  
                
                ·
                
                  t
                  0
                
                ·
                
                  d
                  t
                
              
            
          
          =
          
            
              M
              sd
              
                (
                1
                )
              
            
            /
            
              M
              sd
              
                (
                0
                )
              
            
          
        
      
    
    
      
        (
        7
        )
      
    
  




The description of the chain of physical processes is the following; each photon is:

    
    
        emitted by the source as per a probability density described by the time response of the source,
        and then propagates, following the probability described by the function Gs(r,t), as far as a point r of the mesh,
        it is then absorbed according to the concentration and effective cross-section of the fluorophore at r,
        and then emitted by the fluorophore at the fluorescence wavelength at r,
        it then propagates as far as the detector d as per a probability Gd(r,t)
        and it is finally detected by the detector, which has a time response.
    
    


Owing to the inverse return properties of light, Gd(r,t) is the solution of (1) by replacing rs by rd(rd being the position of the detector under consideration).
As a result, for a “Dirac” fluorophore arranged at r, the measurement Msd performed by the detector d for a medium illuminated by the source s is written:

Msd(t)=S(t)*Gs(r,t)*F(r,t)*Gd(r,t)*D(t)  (8)

The sign * represents temporal convolution, F(t) is the time response of the fluorophore, S(t) is the temporal function of the source and D(t) is the time response of the detector.
In practice, one is capable of determining the source signal detected by the detector, making it possible to arrive at a convolution of the respective time responses of the source and detector. In other words, it is not necessary to estimate the time response of the source and the time response of the detector separately i.e. the convolution S(t)*D(t).
Here, we have assumed that the temporal functions of the sources and detectors are independent of their position; however, if this were not the case, this problem could easily be treated by writing Ss and Dd, the indices s and d being the indices of the different sources and different detectors.
For an assembly of fluorophores in volume Ω, we obtain:




  
    
      
        
          
            M
            sd
          
          ⁡
          
            (
            t
            )
          
        
        =
        
          
            
              ∫
              
                ∫
                ∫
              
            
            Ω
          
          ⁢
          
            S
            ⁡
            
              (
              t
              )
            
          
          *
          
            
              G
              s
            
            ⁡
            
              (
              
                r
                ,
                t
              
              )
            
          
          *
          
            F
            ⁡
            
              (
              
                r
                ,
                t
              
              )
            
          
          *
          
            
              G
              d
            
            ⁡
            
              (
              
                r
                ,
                t
              
              )
            
          
          *
          
            D
            ⁡
            
              (
              t
              )
            
          
          ⁢
          
            d
            r
          
        
      
    
    
      
        (
        9
        )
      
    
  




Owing to the computing properties of moments on the products of convolutions, for the first two moments of Msd(t) we obtain:




  
    
      
        
          M
          sd
          0
        
        =
        
          
            
              ∫
              
                ∫
                ∫
              
            
            Ω
          
          ⁢
          
            
              S
              
                (
                0
                )
              
            
            ·
            
              
                G
                s
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
            ·
            
              
                F
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
            ·
            
              
                G
                d
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
            ·
            
              D
              
                (
                0
                )
              
            
          
          ⁢
          
            d
            r
          
        
      
    
    
      
        (
        10
        )
      
    
  
  
    
      
        
          M
          sd
          
            (
            1
            )
          
        
        =
        
          
            
              ∫
              
                ∫
                ∫
              
            
            Ω
          
          ⁢
          
            d
            r
          
          ⁢
          
            {
            
              
                
                  
                    
                      
                        S
                        
                          (
                          1
                          )
                        
                      
                      ·
                      
                        
                          G
                          s
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          F
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          G
                          d
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        D
                        
                          (
                          0
                          )
                        
                      
                    
                    +
                  
                
              
              
                
                  
                    
                      
                        S
                        
                          (
                          0
                          )
                        
                      
                      ·
                      
                        
                          G
                          s
                          
                            (
                            1
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          F
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          G
                          d
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        D
                        
                          (
                          0
                          )
                        
                      
                    
                    +
                  
                
              
              
                
                  
                    
                      
                        S
                        
                          (
                          0
                          )
                        
                      
                      ·
                      
                        
                          G
                          s
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          F
                          
                            (
                            1
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          G
                          d
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        D
                        
                          (
                          0
                          )
                        
                      
                    
                    +
                  
                
              
              
                
                  
                    
                      
                        S
                        
                          (
                          0
                          )
                        
                      
                      ·
                      
                        
                          G
                          s
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          F
                          
                            (
                            0
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        
                          G
                          d
                          
                            (
                            1
                            )
                          
                        
                        ⁡
                        
                          (
                          r
                          )
                        
                      
                      ·
                      
                        D
                        
                          (
                          0
                          )
                        
                      
                    
                    +
                  
                
              
              
                
                  
                    
                      S
                      
                        (
                        0
                        )
                      
                    
                    ·
                    
                      
                        G
                        s
                        
                          (
                          0
                          )
                        
                      
                      ⁡
                      
                        (
                        r
                        )
                      
                    
                    ·
                    
                      
                        F
                        
                          (
                          0
                          )
                        
                      
                      ⁡
                      
                        (
                        r
                        )
                      
                    
                    ·
                    
                      
                        G
                        d
                        
                          (
                          0
                          )
                        
                      
                      ⁡
                      
                        (
                        r
                        )
                      
                    
                    ·
                    
                      D
                      
                        (
                        1
                        )
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




Taking inspiration from (7) we obtain:




  
    
      
        
          M
          sd
          
            (
            1
            )
          
        
        =
        
          
            
              ∫
              
                ∫
                ∫
              
            
            Ω
          
          ⁢
          
            
              d
              r
            
            ·
            
              S
              
                (
                0
                )
              
            
            ·
            
              
                G
                s
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
            ·
            
              
                F
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
            ·
            
              
                G
                d
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                r
                )
              
            
            ·
            
              
                D
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                
                  
                    T
                    s
                  
                  +
                  
                    
                      T
                      s
                    
                    ⁡
                    
                      (
                      r
                      )
                    
                  
                  +
                  
                    T
                    ⁡
                    
                      (
                      r
                      )
                    
                  
                  +
                  
                    
                      T
                      d
                    
                    ⁡
                    
                      (
                      r
                      )
                    
                  
                  +
                  
                    T
                    d
                  
                
                )
              
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




In this formula:

    
    
        Ts is the mean time of the source function,
        Ts(r) is the mean time-of-flight between s and r,
        T(r) is the lifetime τ of the fluorophore, assumed to be invariant as per position,
        Td(r) is the mean time-of-flight between r and the detector d,
        Td is the mean response time of the detector.
    
    


In infinite, homogeneous cases, according to formulas (4) and (5) it is seen that Ts(r), the mean time-of-flight between s and r, is just equal to:

|r−rs|/(2c√{square root over (μaD)}).

Therefore, v=2c√{square root over (μaD)} can be interpreted as the mean velocity of the photons in this medium.
With regard to the prostate, typically v=5 cm/ns is obtained. The same applies to Td(r).
If the fluorophore is a single, point fluorophore, the formulas (10) and (12) become simpler.
If the fluorophore lies at position r0, the formulas become:

Isd=Msd(0)=S(0)·Gs(0)(r0)·F(0)(r0)·Gd(0)(r0)·D(0)  (13)

and




  
    
      
        
          T
          sd
        
        =
        
          
            
              M
              sd
              
                (
                1
                )
              
            
            
              M
              sd
              
                (
                0
                )
              
            
          
          =
          
            
              T
              s
            
            +
            
              
                T
                s
              
              ⁡
              
                (
                
                  r
                  0
                
                )
              
            
            +
            τ
            +
            
              
                T
                d
              
              ⁡
              
                (
                
                  r
                  0
                
                )
              
            
            +
            
              T
              d
            
          
        
      
    
    
      
        (
        14
        )
      
    
  




S(0)D(0) and Ts+Td can be measured during a calibration operation, during which the source is placed facing the detector. This gives a calibration signal whose zero-order moment is equal to S(0)D(0) and the first-order normalized moment equals to Ts+Td
Therefore, the measured time-of-flight is equal to the mean lag of the source to which are added the time-of-flight between the source and the fluorophore, the lifetime of the fluorophore, the time-of-flight between the fluorophore and the detector and the response time of the detector.
Let us now move on to the problem of fluorescence reconstruction, first with regard to an assembly of fluorophores. It is ascertained, in back-diffusion, that the reconstruction of fluorescence maps in a diffuse medium is very difficult.
The inventors have found, notably in reflection geometry, that different configurations lead to identical measurements. For example, a fluorophore that is broadly distributed can produce a measurement equivalent to a point fluorophore. It therefore appears that during the reconstruction of fluorophores in reflection, there is a strong need for prior hypotheses. One type of hypothesis that can be used is to state that the fluorophores are point fluorophores located at given sites. Said hypothesis can be called a supporting constraint hypothesis.
One first solution is to consider that the fluorophore is single and located as a point in the medium being examined.
Another solution is to consider a plurality of point fluorophores distributed in the examined medium.
The above formula (14) can be used to reconstruct fluorescence by “triangulation”, for example using the technique described in document EP 1 884 765.
In particular, when the medium is homogeneous infinite:

Tsd−Ts−τ−Td=(|r−rs|+|r−rd|)/(2c√{square root over (μaD)}).

Therefore, by combining several source-detector pairs, it is possible to proceed back to position r. The disadvantage of this method is that it provides as many solutions as there are source-detector pairs.
For this purpose, the invention proposes a method to locate at least one fluorophore, in a diffusing medium, using at least one pulsed radiation source capable of emitting radiation to excite this fluorophore and at least one detector capable of measuring a fluorescence signal emitted by this fluorophore, which comprises:

    
    
        illuminating the medium by a radiation source,
        detecting, by at least one detector, the signal produced by the medium at the fluorescence wavelength,
        for at least one source-detector pair, performing temporal distribution Msd(t) of the signal received by the detector,
        the diffusing medium being discretized into M voxels,
    
    


characterized in that the method also comprises the computing of at least one basic parameter χNi,m, which, for at least one of said temporal distributions Msd(t), combines at least one magnitude obtained from at least one moment of said distribution, with at least one estimation of this magnitude, this estimation being made by considering that there are N fluorophores, each then occupying a voxel of the medium as per a distribution m.
Said method may also comprise the determination of a combined parameter PmN, combining at least one or two basic parameters χNi,m.
At least one of the basic parameters may be:

    
    
        a first basic parameter χN1,m, which is the sum, for all pairs (source, detector), of the differences between the value of the measured intensity Msd0 for each source-detector pair, and an estimation of zero-order moment for each source-detector pair obtained using Green functions Gsmn and Gmnd, for the source and detector of each pair, this estimation being made by considering that the N fluorophores are distributed in the voxels of the medium as per configuration m in which the fluorophores are distributed in the voxels mn,
        or a second basic parameter χN1,m, which is the sum, for all source-detector pairs, of the differences for each source-detector pair between the mean measured source-detector time-of-flight (first-order normalized moment of M(t)) corrected by known temporal magnitudes relating to the source, to the detector and to the fluorophore, and the estimation of this corrected time-of-flight, this estimation being made by considering that the N fluorophores are distributed in the voxels of the medium as per configuration m.
    
    


For example, the basic coefficient χ1mN can be such that:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...α
          N
        
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                
                  (
                  0
                  )
                
              
              -
              
                
                  M
                  
                    sd
                    ,
                    m
                  
                  
                    theo
                    ⁢
                    
                      
                    
                    ⁢
                    0
                  
                
                ⁡
                
                  (
                  
                    
                      α
                      1
                    
                    ,
                    
                      ...α
                      N
                    
                  
                  )
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  




where:




  
    
      M
      
        sd
        ,
        m
      
      
        theo
        ⁢
        
          
        
        ⁢
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...α
          N
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          sm
          n
        
      
      ⁢
      
        G
        
          
            m
            n
          
          ⁢
          d
        
      
    
  




    
    
        a configuration m corresponding to N fluorophores distributed in the voxels mn, 1≦n≦N,
        Gsmn (respectively Gmnd) representing the energy transfer functions between the source and the voxel mn (respectively the voxel mn and detector d). The coefficients αn, which can be considered as intensities, are obtained by the minimization operation, as is χ1m;
        each element (Msd0−Msd,mtheo0 (α1, . . . αN)) of the sum then representing the difference between:
        the value of the measured intensity Msd0 for the selected source-detector pair,
        and an estimation of the zero-order moment Msd0, obtaining using firstly the Green functions Gsmn and Gmnd, and secondly a set of coefficients αn each of which is assigned to a fluorophore and which represents the intensity of fluorescence emission by the latter.
    
    


The second basic parameter χ2mN can be such that:




  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
          ′
        
        ,
        
          ...α
          N
          ′
        
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                (
                
                  
                    T
                    sd
                  
                  -
                  
                    T
                    s
                  
                  -
                  σ
                  -
                  
                    T
                    d
                  
                
                )
              
              -
              
                (
                
                  
                    T
                    
                      sd
                      ,
                      m
                    
                    theo
                  
                  ⁡
                  
                    (
                    
                      
                        α
                        1
                        ′
                      
                      ,
                      
                        ...α
                        N
                        ′
                      
                    
                    )
                  
                
                )
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              T
              sd
            
            )
          
        
      
    
  




where:





  
    
      T
      
        sd
        ,
        m
      
      theo
    
    ⁡
    
      (
      
        
          α
          1
          ′
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            a
            N
            ′
          
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  




Ts, Td and τ respectively representing the response time of the source, of the detector and the duration of fluorescence of the fluorophore, Tsmn and Tmnd representing the respective times-of-flight between the source and the fluorophore located in the voxel mn, and between the fluorophore located in voxel mn and the detector, the coefficient σ2 (Tsd) corresponding to an estimation of the variance of distribution Tsd, Ttheosd (α′1 . . . α′N) then representing an estimation of the first-order normalized moment of function Msd(t), the coefficients α′N being obtained with the minimization operation.
With said method, it is also possible to determine a combined parameter PmN, equal to the sum or to the product of the basic coefficients χN1,m, and χN2,m, the combined parameters PmN of lowest value then corresponding to the fluorophore distributions the closest to the effective fluorophore distribution.
The number of fluorophores can equal 1, each basic parameter χN=11,m, and χN=12,m, then being determined by considering that the fluorophore is single and is located in voxel m. χN=11,m, and χN=12,m can then be determined according to the following equations:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      min
      
        α
        m
      
    
    ⁢
    
      
        ∑
        
          s
          ,
          d
        
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                0
              
              -
              
                
                  
                    α
                    m
                  
                  ·
                  
                    G
                    sm
                  
                
                ⁢
                
                  G
                  md
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  




and




  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      ∑
      
        s
        ,
        d
      
    
    ⁢
    
      
        
          (
          
            
              (
              
                
                  T
                  sd
                
                -
                
                  T
                  s
                
                -
                τ
                -
                
                  T
                  d
                
              
              )
            
            -
            
              (
              
                
                  T
                  sm
                
                +
                
                  T
                  
                    d
                    ⁢
                    
                      
                    
                    ⁢
                    m
                  
                
              
              )
            
          
          )
        
        2
      
      
        
          σ
          2
        
        ⁡
        
          (
          
            T
            sd
          
          )
        
      
    
  




Msd0 representing the zero-order moment of the function Msd(t),
Gsm and Gmd being energy transfer functions, for each source-detector pair sd, between respectively the source s and the voxel m and between voxel m and the detector d,
Tsm and Tdm representing the respective times-of-flight between the source and voxel m, and voxel m and the detector d,
Ts, Td and τ respectively representing the response times of the source, of the detector and the duration of fluorescence of the fluorophore, the coefficient σ2(Tsd) corresponding to an estimation of the variance of distribution Tsd.
With this method, it is also possible to determine a combined parameter PmN, equal to the sum or to the product of the basic coefficients χN=11,m, and χN=12,m, the combined parameters PmN=1 of lowest value then corresponding to the distributions of the fluorophores the closest to the effective distribution of the fluorophores.
The invention also proposes a method to locate at least one fluorophore or at least one absorber in a diffusing medium, using at least one pulsed radiation source capable of emitting radiation to excite this fluorophore or this molecule, and at least one detector capable of measuring a fluorescence signal (Φfluo) emitted by this fluorophore or an emission signal emitted by this absorber, comprising:
a) for at least one pair (radiation source, detector), at least one excitation by radiation derived from the radiation source of this pair, and at least one detection by the detector of this pair of fluorescence signal emitted by this fluorophore after this excitation, or of the emission signal emitted by this absorber,
b) identification of meshing (M) of the volume into mesh elements (m);
c) an estimation of the location of the fluorophore or of the absorber in its diffusing medium, by computing a function (χ1m, χ′1m, χ2m, χ′2m, Pm) of at least one of the three following parameters, for each element of the mesh:

    
    
        χ1m (or χ1mN) which is the sum, for all source-detector pairs, of the differences between the value of measured intensity Msd0 for each source-detector pair, optionally corrected by estimation of intensity contribution, measured by the detector, of already located fluorophores, and an estimation of the zero-order moment for each source-detector pair obtained using the Green functions Green Gsm and Gmd (or Gsmn and Gmnd) for the source and detector of each pair and the mesh element;
        χ2m (or χ2mN) which is the sum, for all source-detector pairs, of the differences for each source-detector pair between the mean measured source-detector time-of-flight (first-order normalized moment of M(t)), optionally corrected by estimation of the time-of-flight contribution by already located fluorophores, and this same, modelled, first-order normalized moment,
        χ3m which is the sum, over all the source-detectors pairs, of the correlations for each source-detector pair between firstly a value of the measured intensity (Isd=M(0) (t)) determined by the zero-order moment of M(t) denoted M(0) (t), optionally corrected by estimation of the intensity contribution measured by the detector of already located fluorophores, and secondly a value of this zero-order moment modelled using the product Gsm·Gdm of the Green functions Gsm and Gmd, for the source and detector of each pair and the mesh element.
    
    


χ1mN, χ2mN, χ3m can have the form already indicated above, or as in the remainder of this application.
Said function of at least one of the three parameters can be a linear function of at least one of said parameters and/or a non-linear function of at least one of said parameters.
Therefore, for example, another parameter is χN=14,m:




  
    χ
    
      4
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            s
          
        
      
    
    ⁢
    
      [
      
        
          
            ∑
            sd
          
          ⁢
          
            
              
                (
                
                  
                    M
                    sd
                    
                      (
                      0
                      )
                    
                  
                  -
                  
                    
                      M
                      sd
                      
                        theo
                        ⁢
                        
                          
                        
                        ⁢
                        0
                      
                    
                    ⁡
                    
                      (
                      
                        
                          α
                          1
                        
                        ,
                        
                          ...
                          ⁢
                          
                            
                          
                          ⁢
                          
                            α
                            N
                          
                        
                      
                      )
                    
                  
                
                )
              
              2
            
            
              
                σ
                2
              
              ⁡
              
                (
                
                  M
                  sd
                  0
                
                )
              
            
          
        
        +
        
          
            ∑
            sd
          
          ⁢
          
            
              
                (
                
                  
                    (
                    
                      
                        T
                        sd
                      
                      -
                      
                        T
                        s
                      
                      -
                      σ
                      -
                      
                        T
                        d
                      
                    
                    )
                  
                  -
                  
                    (
                    
                      
                        T
                        sd
                        theo
                      
                      ⁡
                      
                        (
                        
                          
                            α
                            1
                          
                          ,
                          
                            ...
                            ⁢
                            
                              
                            
                            ⁢
                            
                              α
                              N
                            
                          
                        
                        )
                      
                    
                    )
                  
                
                )
              
              2
            
            
              
                σ
                2
              
              ⁡
              
                (
                
                  T
                  sd
                
                )
              
            
          
        
      
      ]
    
  




where:




  
    
      M
      sd
      
        theo
        ⁢
        
          
        
        ⁢
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          sm
          n
        
      
      ⁢
      
        G
        
          
            m
            n
          
          ⁢
          d
        
      
    
  




and:




  
    
      T
      sd
      theo
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  




An operator selects one or more of these parameters and/or one or more functions of these parameters, and computes the same to obtain initial information on the location of at least one fluorophore, optionally after comparing the computed results with one or more threshold values, or one or more minimum values, or one or more maximum values.
Irrespective of the embodiment of the invention, it is additionally possible, for the locating of several fluorophores, to carry out a step to adjust all the intensity contributions of already located fluorophores to the measured intensity.
It is also possible to perform a step to compute a correlation between the measured intensity and all the intensity contributions of all the fluorophores.
If there are several fluorophores of unknown number, the method may further comprise a step (S′1) to increment the number of fluorophores if the computing of a correlation is not satisfactory.
The medium surrounding the fluorophore or the absorber may be of infinite type or semi-infinite type or in slab form limited by two parallel surfaces, or of any shape its outer surface being discretized into a series of planes.
One method according to the invention is particularly well suited to the case in which the sources of radiation and the detectors have reflection geometry.
The invention also concerns a device to locate at least one fluorophore, in a diffusing medium, comprising at least one pulsed radiation source capable of emitting radiation to excite this fluorophore, and at least one detector capable of measuring a fluorescence signal emitted by this fluorophore, comprising:

    
    
        means, for at least one source-detector pair, to perform temporal distribution Msd(t) of the signal received by the detector,
        means to produce a mesh (M) of the volume in mesh elements m;
        means to compute at least one basic parameter χNi,m, which, for at least one of said temporal distributions Msd(t), combines at least one magnitude obtained from at least one moment of said distribution, and at least one estimation of this magnitude, this estimation being made by considering that there are N fluorophores, each then occupying a voxel of the medium as per a distribution m.
    
    


The basic coefficient χ1mN can be such that:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            a
            N
          
        
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                
                  (
                  0
                  )
                
              
              -
              
                
                  M
                  
                    sd
                    ,
                    m
                  
                  
                    theo
                    ⁢
                    
                      
                    
                    ⁢
                    0
                  
                
                ⁡
                
                  (
                  
                    
                      α
                      1
                    
                    ,
                    
                      ...
                      ⁢
                      
                        
                      
                      ⁢
                      
                        α
                        N
                      
                    
                  
                  )
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  




    
    
        where:
    
    






  
    
      M
      
        sd
        ,
        m
      
      
        theo
        ⁢
        
          
        
        ⁢
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          sm
          n
        
      
      ⁢
      
        G
        
          
            m
            n
          
          ⁢
          d
        
      
    
  




    
    
        
        
            a configuration m corresponding to N fluorophores distributed in the voxels mn, 1≦n≦N,
            Gsmn(respectively Gmnd) representing the energy transfer functions between the source and the voxel mn (respectively the voxel mn and the detector d). The coefficients αn, which can be considered as intensities, are obtained with the minimization operation, as is χ1m;
            Each element (Msd0−Msd,mtheo0 (α1, . . . αN)) of the sum then representing the difference between:
            the value of the measured intensity Msd0 for the selected source-detector pair,
            and an estimation of the zero-order moment Msd0, obtained using firstly the Green functions Gsmn and Gmnd, and secondly a set of coefficients αn of which each one is assigned to a fluorophore and represents the intensity of the fluorescence emission thereof.
        
        
    
    


And the second basic parameter χ2mN can be such that:




  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
          ′
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
            ′
          
        
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                (
                
                  
                    T
                    sd
                  
                  -
                  
                    T
                    s
                  
                  -
                  σ
                  -
                  
                    T
                    d
                  
                
                )
              
              -
              
                (
                
                  
                    T
                    
                      sd
                      ,
                      m
                    
                    theo
                  
                  ⁡
                  
                    (
                    
                      
                        α
                        1
                        ′
                      
                      ,
                      
                        ...
                        ⁢
                        
                          
                        
                        ⁢
                        
                          α
                          N
                          ′
                        
                      
                    
                    )
                  
                
                )
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              T
              sd
            
            )
          
        
      
    
  




where:




  
    
      T
      
        sd
        ,
        m
      
      theo
    
    ⁡
    
      (
      
        
          α
          1
          ′
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
            ′
          
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  




Ts, Td and τ respectively represent the response time of the source, of the detector and the duration of fluorescence of the fluorophore, Tsmn and Tmnd representing the respective times-of-flight between the source and the fluorophore located in the voxel mn, and between the fluorophore located in voxel mn and the detector, the coefficient σ2(Tsd) corresponding to an estimation of the variance of distribution Tsd, Ttheosd,m (α′1 . . . α′N) then representing an estimation of the first-order normalized moment of function Msd(t), the coefficients α′N being obtained with the minimization operation.
Said device may comprise means to determine a combined parameter PmN, equal to the sum or to the product of the basic coefficients χN1,m, and χN2,m. 
If the number of fluorophores is equal to 1, each basic parameter χN=11,m, and χN=12,m, can be determined by considering that the fluorophore is single and located in voxel m.
χN=11,m, and χN=12,m, can be determined according to the following equations:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      min
      
        α
        m
      
    
    ⁢
    
      
        ∑
        
          s
          ,
          d
        
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                0
              
              -
              
                
                  
                    α
                    m
                  
                  ·
                  
                    G
                    sm
                  
                
                ⁢
                
                  G
                  md
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  




and




  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      ∑
      
        s
        ,
        d
      
    
    ⁢
    
      
        
          (
          
            
              (
              
                
                  T
                  sd
                
                -
                
                  T
                  s
                
                -
                τ
                -
                
                  T
                  d
                
              
              )
            
            -
            
              (
              
                
                  T
                  sm
                
                +
                
                  T
                  
                    d
                    ⁢
                    
                      
                    
                    ⁢
                    m
                  
                
              
              )
            
          
          )
        
        2
      
      
        
          σ
          2
        
        ⁡
        
          (
          
            T
            sd
          
          )
        
      
    
  




Msd0 representing the zero-order moment of the function Msd(t),
Gsm and Gmd being energy transfer functions, for each source-detector pair sd, between the source s and voxel m and between voxel m and detector d respectively,
Tsm and Tdm representing the respective times-of-flight between the source and voxel m, and voxel m and the detector d,
Ts, Td and τ respectively representing the response times of the source, of the detector and the duration of fluorescence of the fluorophore,
the coefficient σ2(Tsd) corresponding to an estimation of the variance of distribution Tsd.
The invention also concerns a device to locate at least one fluorophore or at least one absorber in a diffusing medium, comprising at least one radiation source capable of emitting radiation to excite this fluorophore or this molecule, and at least one detector capable of measuring a fluorescence signal (Φfluo) emitted by this fluorophore or an emission signal emitted by this absorber, comprising:
a) means for producing a mesh (M) of the volume in mesh elements m;
b) means to estimate the location of the fluorophore or of the absorber in its diffusing medium, by computing a function (Pm) of at least one of the three following parameters:

    
    
        χ1m (or χ1mN) which is the sum over all source-detector pairs, of the differences between the value of the measured intensity Msd0 for each selected source-detector pair, optionally corrected by the estimation of the contribution to the intensity made by already located fluorophores, and an estimation of the zero-order moment obtained using the Green functions Gsm and Gmd (or Gsmn and Gmnd),
        χ2m (or χ2mN) which is the sum over all source-detector pairs of the differences, for at least one measurement, between the mean measured source-detector time-of-flight (first-order normalized moment of M(t)), optionally corrected by the estimation of the contribution of already located fluorophores to the time-of-flight, and this same modelled, first-order normalized moment,
        χ3m which is the sum over all source-detector pairs of the correlations between firstly a value of the measured intensity (Isd=M(0) (t)) determined by the zero-order moment of M(t), denoted M(0) (t), optionally corrected by the estimation the contribution, to the intensity, by already located fluorophores, and secondly a value of this modelled zero-order moment (the product Gsm·Gdm).
    
    


Said function of at least one of the three parameters can be a linear function of at least one of said parameters and/or a non-linear function of at least one of said parameters.
Whatever its embodiment, said device may comprise means, if there are a plurality of fluorophores, to perform a step to adjust all intensity contributions by already located fluorophores to the measured intensity.
It may additionally comprise means to compute a correlation between the measured intensity and all the intensity contributions by all the fluorophores.
Means can be provided to increment the number of fluorophores if the computing of a correlation is not satisfactory.
In said device, the sources of radiation and the detectors preferably have reflection geometry.



BRIEF DESCRIPTION OF THE DRAWINGS
FIGS. 1A-1B schematically illustrate an experimental configuration of fluorescence tomography in reflection geometry, which can be implemented under the invention.
FIGS. 1C-1E respectively show a series of laser pulses and single emitted photons, and a fluorescence curve obtained from data on the single photons, and curves of measured fluorescence.
FIG. 2 is a schematic of a configuration of fluorophores, sources and detectors used for tests.
FIGS. 3A-3C, 4A-4C illustrate the application of non-combined criteria, according to the invention.
FIGS. 5A-5B illustrate the application of a combined criterion according to the invention.
FIG. 6 schematically illustrates a configuration of 2 fluorophores, 2 sources and 2 detectors.
FIGS. 7 and 8 are examples of a method according to the invention.
FIGS. 9A-9C and 10A-10H illustrate examples of application of a method according to the invention.
FIGS. 11A-11C and 12A-12C illustrate the results obtained with a method according to the invention.



DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION
An experimental device in reflection geometry (in which the sources and detectors lie on the same face of the object or of the examined medium) is illustrated in FIG. 1A.
The radiation source and the detection means respectively may be the tip of an optical fibre conveying excitation light into the medium and the tip of an optical fibre which collects part of the emission light.
The diffusing medium 2 can be seen in the figure, containing fluorescent inclusions Fi (here one inclusion F1 is shown) and a set of sources Si (here each of the sources e.g. a laser source Li, is connected to a fibre, therefore the source can be considered to be the tip of the fibre closest to the medium 2) sequentially bringing the light into the medium.
A set of detection points is designated Di: here each of the detectors PMi is connected to a fibre, therefore the detector can be considered to be the tip of the fibre closest to the medium 2. These detectors are used to sample the density of fluorescence photons.
In general, the number ns of sources can be any number. For example, there may be only one source, or 2 or 3, or any number ns>3.
The number nd of detectors can be any number. For example, there may be a single detector, or 2 or 3 or any number nd>3.
Preferably ns>2 or >3 or >5 and nd>2 or >3 or >5.
In applications in which great compactness is required, in endoscopy for example, the number ns of sources can lie between 1 and 10.
Each source Li can be a pulsed laser source or an optical fibre connected to a remote pulsed laser (which is the case for the configuration in FIG. 1A). It is also possible to have a single source connected to several fibres (configuration in FIG. 1B in which there is a single laser source L, the other aspects of the system being the same as those in FIG. 1A), and means e.g. an optical switch or multiplexer to select the fibre in which the beam is to be sent.
Each detector Di may comprise an array of camera pixels of CCD or CMOS type, or may comprise a photomultiplier, and it may be connected to an optical fibre which transmits the signal from the examined medium towards the detector, which is the case for the configuration in FIGS. 1A and 1B.
According to one embodiment, at the output of the detector(s), temporal analysis of the signals is performed using a dedicated measuring instrument (TCSPC—Time Correlated Single Photon Counting to count time correlated photons or intensified camera with temporal gate).
One of the pulse sources of radiation Li may also be used as means to trigger the TCSPC card (for example via a link 91, 92, 9 between each source Li and the means 24).
It is preferable to work with pulses in the femtosecond domain, provided the radiation source is adequate i.e. it comprises one or more laser sources Li each pulse of which has a temporal width also lying in the femtosecond domain. The use of pulsed sources of picosecond laser diode type can also be envisaged, the advantage of this solution being lower cost.
According to one particular embodiment, one or more sources Li are a pulsed laser diode, for example at a wavelength close to 630 nm and with a repeat rate of about 50 MHz.
The laser light preferably passes through an interference filter to remove any secondary modes. An interference filter and/or a coloured filter absorbing low wavelengths can be placed in front of each detector PMi or the tip Di of the corresponding fibre to select the fluorescence light (for example: λ>650 nm, the source being at a wavelength of 631 nm for example) of a fluorophore Fi arranged in the medium 2 and to optimize elimination of the excitation light.
In the text of the present application, irrespective of the embodiment of the invention, one particular issue relates to the position of the radiation source and/or the position of the detection means.
When fibres are used, these positions are most often meant to be those of the tips of the fibres which convey the radiation into the diffusing medium 2, or to the interface of this medium and/or those which collect the diffused radiation, the latter being placed in the medium or at the interface with this medium, and in this case they are not to be understood as the positions of the source Li properly so-called or of the detector PMi properly so-called. The fibres conveying the radiation into the medium can be called excitation fibres, and the fibres collecting the diffused radiation can be called emission fibres or collection fibres.
In applications of “endoscopy” type, use is preferably made of a fibre system, the fibres possibly being integrated in an endo-rectal probe for example, at least one fibre transmitting the source of pulsed excitation light towards the region to be examined.
This source of pulsed light may be external to the probe.
At least one other fibre transmits the optical emission signal of the region to be examined towards a photodetector which may also be external to the probe. With the Time Correlated Single Photon Counting technique (TCSPC), a photon is detected using a photomultiplier that is emitted by the fluorophore after a pulse of the radiation source.
The system therefore allows time-resolved detection of the fluorescence pulses.
It allows collection of the fluorescence photons.
FIG. 1C shows a series of laser pulses Ii (i=1-4) and a series of corresponding single photons pi (i=1-4), these being detected by a TCSPC-type system (Time Correlated Single Photon Counting).
Each photon is in fact detected relative to the output of the corresponding pulse: in FIG. 1C, ti represents the time elapsed between each laser pulse Ii and the instant of detection of each photon pi.
It is then possible to establish a statistical distribution or histogram of the photon arrival time, as illustrated in FIG. 1D, which shows the number of detected fluorescence photons, as a function of elapsed time t relative to each laser pulse. This histogram represents an experimental measurement of the function Msd(t) mentioned previously, this function Msd(t) being the theoretical histogram of photon arrival time. From said histogram, it is possible to determine statistical parameters, notably the intensity or mean arrival time, or mean time (in fact, the weighted mean of the abscissas with the ordinates of the histogram). Therefore, for each pair grouping together a source and a detector, called source-detector pair and denoted sd, this histogram can be obtained which is likened to the function Msd(t), after a large number of measurements or acquisitions.
Said curve Msd(t) which, as can be seen (see also the example in FIG. 1E), allows all the information to be used over a wide time window either side of the point of maximum intensity (and not only in the rising part of the signal) can therefore be subsequently processed to obtain characteristic data such as arrival time, or mean time as will be explained below.
Electronic means 24 such as a microcomputer or computer are programmed to store in memory and process data derived from the detectors. More precisely, a central unit 26 is programmed to implement a processing method according to the invention: to compute one or more parameters χNi,m or χ3,m or one or more functions of these parameters and/or to apply a method such as described below for example in connection with one of FIG. 7 or 8. These means 26 can optionally be used to make one or more comparisons of the value(s) found for the parameter(s) with one or more threshold values.
An operator selects one or more of these parameters and/or one or more functions of these parameters, and computes the same to obtain initial data on a location of at least one fluorophore, optionally after comparing the computed results with one or more threshold values or with one or more minimum values or with one or more maximum values.
The electronic means 24 can optionally be used to control the triggering of the radiation source(s) and/or of the detector(s), for example for synchronization of all these means.
A monitoring screen or display or viewing means 27 can be used, after processing, to show the positioning of the fluorophore(s) in the examined medium 2. The means 24 and the display means 27 also enable an operator to determine a mesh such as used in the present application, as explained below.
If a locating result is not suitable or is not satisfactory, for example for lack of precision, at least one other parameter can be chosen or at least one iteration can be performed.
Other detection techniques can be used, for example with an intensified camera e.g. an ultra-rapid intensified camera of “gated camera” type; in this case, the camera opens onto a time gate of width about 200 ps for example, then this gate is shifted e.g. in steps of 25 ps.
The excitation light, at wavelength λx excites the fluorophores Fi which re-emit a so-called emission light at wavelength λm>λx with a lifetime τ. This lifetime corresponds to the mean duration during which the excited electrons remain in this state before returning to their initial state. For cases not concerning fluorescence imaging, the method can be used to locate absorbers in the medium i.e. sites at which the absorption coefficient of the medium is higher than the mean absorption coefficient of the medium. In the remainder of this application, said absorber can be likened to a fluorophore, by considering that its fluorescence wavelength, or emission wavelength, is equal to the excitation wavelength, and also by considering that the duration of fluorescence is zero, i.e τ=0. In this case, the filter which may optionally be placed upstream of the detector is to be adapted to the excitation wavelength.
Irrespective of the medium and geometry under consideration one or more fluorescent inclusions are located in the medium, which has a light propagation velocity c (=3*108 m/s/n, n=optical index), a diffusion coefficient D and an absorption coefficient μa.
The medium 2 surrounding the fluorophore F1 may have any geometric shape; it may for example be of infinite type.
It may also be of semi-infinite type, limited by a wall. In this case, the excitation and emission optical fibres can be placed at the interface of the diffusing medium, or at a short depth of the order of a few mm.
A semi-infinite medium, in which each of the fibre tips conveying an excitation signal and collecting a fluorescence signal is arranged at more than 1 cm, or more than 1.5 cm or more than 2 cm from the wall limiting this medium can be treated, according to the method of the invention, as an infinite medium.
The surrounding medium may also form a slab limited by two parallel surfaces.
The invention also applies to a medium of any shape, whose outer surface is broken down into a series of planar facets.
In all cases, it is particularly well suited to the case of reflection geometry in which the source(s) and the detector(s) lie on the same face of the object or of the examined medium.
The Green functions Gs(r,t) of a device according to the invention can be determined using a Monte-Carlo simulation, for example to solve the radiation transfer equation, or most often by solving the diffusion equation (1) to which boundary conditions are added (partial current boundary condition [Robin] or extrapolated boundary condition [Dirichlet]), as explained in the article by S. R. ARRIDGE: “Optical tomography in medical imaging”, Inverse Problems, April 1999, vol. 15(2), R41-R93. These determinations or simulations can be carried out using computing means such as the means 26, 27 in FIGS. 1A and 1B. The moments of the Green functions Gs(r,t), such as defined previously (1′), can be obtained by iteratively solving equation (2).
According to one embodiment of the invention, a fluorophore is a single fluorophore F1, or a plurality of fluorophores substantially grouped together in one same place.
According to another embodiment, the fluorophore comprises N fluorophores distributed in different places.
Using time-resolved optical measurements, it is sought to locate these fluorescent inclusions i.e. to reconstruct the fluorescence map.
According to the invention, the fluorescence map is reconstructed using an efficient method from a computing viewpoint, and which has the property of concentrating the fluorescence on small supports. The term “small supports” represents the fact that the fluorophore(s) are discretely distributed in certain places of the examined medium. The inventors considered that this type of prior information is well suited to the reconstruction of fluorophores in diffusing tissues, notably in some applications related to the early detection of cancerous tumours in organs such as the prostate, breast, testicles, brain, etc.
In one first embodiment, the case of a single fluorophore is considered.
Depending on the method used, the data on M0sd and M1sd is applied (already defined above, see for example Equations 11-14) to reconstruct the position of the fluorophores which emit fluorescence radiation.
In all cases, the fluorophore search volume is meshed into M mesh elements (or voxels) m, whose size depends on desired accuracy: the higher the precision the smaller the individual volume of a mesh element, and the greater the number of mesh elements.
This mesh can be decided by an operator of the computer system 24.
According to the invention, one or more first parameters χim are used called basic parameters (or criteria, but in the remainder hereof the expression “parameters” will be used) which are defined below, or a combination thereof, this combination then possibly being called a combined parameter Pm. According to this embodiment, the index m designates a voxel.
Each of these basic parameters depends on the variable m, and hence on the selected mesh element in the volume. All the mesh elements are selected one by one, and for each mesh element the parameter is computed which is a sum of the contributions of all the possible pairs (source-detector) of the system or part of all these pairs. In other words, each parameter represents the sum for different acquisitions, each acquisition corresponding to a source-detector-pair, of a combination between a magnitude (here and in the rest of this document this expression is equivalent to “function” or “quantity” or “parameter”) determined from at least one moment of Msd(t), and the estimation of this magnitude, this estimation being made by considering a single fluorophore positioned in the voxel m.
By combination is meant any arithmetic operation, whether linear or not. It will be seen in the remainder of the description that this combination may assume the form of a subtraction or of a product for example.
It will also be seen that the magnitude determined from at least one moment of the function Msd(t) can be the zero-order moment of this function, or the first-order normalized moment, from which known and assumed constant temporal magnitudes (or “functions” or “parameters”) are subtracted.
According to one preferred embodiment, this magnitude can be the zero-order moment of Msd(t), denoted M(0)sd.
Therefore, a first basic parameter χ1m is defined:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      min
      
        α
        m
      
    
    ⁢
    
      
        ∑
        
          s
          ,
          d
        
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                0
              
              -
              
                
                  
                    α
                    m
                  
                  ·
                  
                    G
                    sm
                  
                
                ⁢
                
                  G
                  md
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  




χ1m, can be considered as a residual between the measurement of moment M(0)sd, and its estimation, given by the product αmGsmGmd. χ1m is minimal for the mesh elements m giving the best positioning of the fluorophore. Its value results from computing of the minimum value of the sum indicated above, for a set of coefficients αm. In other words, the set of coefficients αm is sought which minimizes the above sum.
The term σ2(Msd0) shown in the denominator part is optional. It is an estimation of the uncertainty associated with the determination of Msd0. If Msd0 is the number of photons detected by the pair of detectors (s, d), and if the noise is Poissonian, then σ2(Msd0)=Msd0 
Therefore, each element (Msd0−αmGsmGmd)2 of the sum represents the difference between:

    
    
        the value of the measured intensity Msd0 for the selected source-detector pair, when the source emits a pulse and the detector detects a fluorescence signal emitted by the fluorophore,
        and an estimation of the zero-order moment Msd0, obtained using the Green functions Gsm and Gmd, which can be estimated in the manner already described above, for example using means 24.
    
    


Again according to one preferred embodiment, the magnitude determined from at least one moment of Msd(t) is the first-order normalized moment, denoted Tsd, from which are subtracted the time of the function of the source Ts, the response time of the detector Td, and the fluorescence time τ. In this case, a second parameter χ2m is used:




  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      ∑
      
        s
        ,
        d
      
    
    ⁢
    
      
        
          (
          
            
              (
              
                
                  T
                  sd
                
                -
                
                  T
                  s
                
                -
                τ
                -
                
                  T
                  d
                
              
              )
            
            -
            
              (
              
                
                  T
                  sm
                
                +
                
                  T
                  
                    d
                    ⁢
                    
                      
                    
                    ⁢
                    m
                  
                
              
              )
            
          
          )
        
        2
      
      
        
          σ
          2
        
        ⁡
        
          (
          
            T
            sd
          
          )
        
      
    
  




χ2m, also called the criterion of time-of-flight measurement residual, is also small for mesh elements verifying good temporal agreement, and hence contributing towards better positioning of the fluorophore. More specifically, χ2m corresponds to a difference (or residual), for at least one acquisition, between the mean measured source-detector time-of-flight (first-order normalized moment of M(t)) and this same, modelled, first-order normalized moment, the fluorophore then being assumed to correspond to the voxel m.
In the above expression of this parameter χ2m:

    
    
        Tsd is the measured source-detector time-of-flight, it is the first-order normalized moment of the measurement,
        Ts is the mean time of the source function, previously mentioned;
        Td is the response time of the detector, previously mentioned;
        τ is the lifetime of the fluorophore;
        Tsm is the mean duration of the source-fluorophore trajectory (or source−voxel (m));
        Tdm is the mean duration of the fluorophore-detector trajectory (or detector−voxel (m)).
    
    


The term (Tsd−Ts−Td−□) represents the “net” time-of-flight, i.e. the measured time-of-flight, from which are subtracted the mean times of the source function, of detector response and of fluorescence lifetime. This net time-of-flight represents the duration of the light trajectory between the source, the fluorophore and the detector.
Tsm+Tdm represents an estimation of this net time-of-flight when the fluorophore is present in the voxel m.
The magnitudes Ts, Td, τ are data which can be known: τ depends on the fluorophore, whilst the sum Ts+Td can be obtained by calibration testing as previously described. Tsm and Tdm are modelled either analytically or by solving equation (2), by considering that the fluorophore is single and located in the voxel m.
Preferably, and as is the case in the above formula, the parameter χ2m can be normalized by a statistical magnitude (again, here and in the rest of the document, this expression is similar to “function” or “parameter”) relative to the distribution of Tsd. This statistical magnitude acts as confidence indicator. Here σ2 (Tsd) has been chosen, but another statistical magnitude (e.g. standard deviation or variance) could be chosen, to indicate the confidence assigned to a distribution, but normalization by variance is optimal. This statistical magnitude is used to weight the measurements corresponding to a source-detector pair, relative to other measurements corresponding to another source-detector pair: the smaller the standard deviation of the distribution of Tsd, the more the measurements derived from this source-detector pair can be considered to be reliable and the greater the importance they can be given in determining the second coefficient, which is explained below. The magnitude σ2 (Tsd) can easily be estimated from a width ΔTsd of function Msd(t) according to the following equations:



  
    
      Δ
      ⁢
      
        
      
      ⁢
      
        T
        sd
      
    
    =
    
      
        
          
            ∑
            i
          
          ⁢
          
            
              
                (
                
                  
                    t
                    i
                  
                  -
                  
                    T
                    sd
                  
                
                )
              
              2
            
            ⁢
            
              
                M
                sd
              
              ⁡
              
                (
                
                  t
                  i
                
                )
              
            
          
        
        
          M
          sd
          
            (
            0
            )
          
        
      
    
  




  
    
      
        σ
        2
      
      ⁡
      
        (
        
          T
          sd
        
        )
      
    
    =
    
      
        Δ
        ⁢
        
          
        
        ⁢
        
          T
          sd
          2
        
      
      
        M
        sd
        
          (
          0
          )
        
      
    
  



where:

    
    
        ti designates the abscissa i of function Msd(t)
    
    


As for the first basic parameter χ1m, χ2m is minimal for the mesh elements m in which the fluorophore is effectively positioned.
Therefore, when these two first parameters χ1m and χ2m are used, the mesh elements of the volume for which they are both minimum, are mesh elements in which the fluorophore may most probably be located.
Additionally, the inventors have evidenced the particular advantage of jointly using these two basic parameters χ1m and χ2m. When the sources and detectors lie in a plane P (which is typically the case for a reflection configuration), the first basic parameter χ1m gives particularly pertinent information on the coordinates of the fluorophore in this plane. In other words, the voxels m for which the coefficient χ1m is minimum are generally distributed in a direction perpendicular to this plane, but their coordinates in this plane are little dispersed.
In complementary manner, and unexpectedly, the second basic parameter χ2m gives particularly pertinent information on the position of the fluorophore in the direction perpendicular to this plane. For example, a parameter χ2m, combining a measurement and an estimation of the time-of-flight, gives an indication on the depth of the fluorophore in relation to a plane formed by the sources and detectors. In other words, the voxels m for which the coefficient χ2m is minimum are generally distributed substantially parallel to the plane formed by the sources and detectors, and their coordinates in the direction perpendicular to this plane are little dispersed. These voxels are then located at one same depth in relation to the point at which the sources and detectors are located.
The advantage can therefore be appreciated of combining these two first basic parameters: the probability of the presence of the fluorophore in a voxel m will be higher the more the first basic parameter χ1m and second basic parameter χ2m are minimum for this voxel. Owing to the complementarity existing between these two first basic parameters χ1m and χ2m, pertinent information will be obtained regarding the coordinates of the fluorophore in this plane defined by the sources and the detectors (by the first basic parameter χ1m), but also pertinent information regarding the coordinates of the fluorophore in the direction orthogonal to this plane i.e. to the depth (by the second basic parameter χ2m). The coordinates of the voxel likely to contain the fluorophore are thus defined.
A third basic parameter χ3m can be formed by considering that the magnitude determined from the distribution Msd(t) is the zero-order moment, as for the determination of the first basic parameter. This magnitude is then combined with the estimation of this magnitude, said estimation being made by considering that the fluorophore is single and located in voxel m. But here this combination is not a difference but a correlation, in other words a product. Therefore, this third basic parameter corresponds to the correlation, for at least one source-detector pair, between a value of the measured intensity (Isd=M(0)(t)) determined by the zero-order moment of M(t), and a value of this estimated zero-order moment (the product Gsm·Gdm). This coefficient of correlation corresponds to the sum, for each modelled source-detector pair, of the product between the measured value of lsd (intensity for one excitation by the source S and for one detection by the detector D) and the modelled intensity:




  
    
      χ
      
        3
        ⁢
        m
      
    
    =
    
      
        
          ∑
          sd
        
        ⁢
        
          
            I
            sd
          
          ⁢
          
            G
            sm
          
          ⁢
          
            G
            
              d
              ⁢
              
                
              
              ⁢
              m
            
          
        
      
      
        
          
            ∑
            
              s
              ,
              d
            
          
          ⁢
          
            
              G
              sm
              2
            
            ·
            
              G
              
                d
                ⁢
                
                  
                
                ⁢
                m
              
              2
            
          
        
      
    
  
  ,




m being the index corresponding to the voxel under consideration, s and d representing the source and detector under consideration.
It is also possible to normalize this coefficient:




  
    
      ∑
      sd
    
    ⁢
    
      
        I
        sd
      
      ⁢
      
        G
        sm
      
      ⁢
      
        G
        
          d
          ⁢
          
            
          
          ⁢
          m
        
      
    
  
  
    
      
        ∑
        
          s
          ,
          d
        
      
      ⁢
      
        
          G
          sm
          2
        
        ·
        
          G
          
            d
            ⁢
            
              
            
            ⁢
            m
          
          2
        
      
    
  




using:




  
    
      χ
      
        3
        ⁢
        m
      
    
    =
    
      
        
          ∑
          sd
        
        ⁢
        
          
            I
            sd
          
          ·
          
            G
            sm
          
          ·
          
            G
            
              d
              ⁢
              
                
              
              ⁢
              m
            
          
        
      
      
        
          
            
              ∑
              
                s
                ,
                d
              
            
            ⁢
            
              I
              sd
              2
            
          
        
        ·
        
          
            
              ∑
              
                s
                ,
                d
              
            
            ⁢
            
              
                G
                sm
                2
              
              ·
              
                G
                
                  d
                  ⁢
                  
                    
                  
                  ⁢
                  m
                
                2
              
            
          
        
      
    
  
  ,




this normalized coefficient then lying between 0 and 1.
Unlike the two preceding parameters, χ3m is maximal for the mesh elements m giving the best positioning of the fluorophore, notably in a plane substantially parallel to the plane formed by the sources and detectors.
To implement a method according to the invention, at least one of the above basic parameters is chosen, and a combination Pm is formed of these parameters. For example, Pm can be called a combined parameter, resulting from a combination of at least one first parameter (for only one parameter the term “function” of this parameter would best be used, but in the present application for reasons of simplification, the general expression “combination” will be used). Said combination Pm is also a function of m, i.e. of the voxel under consideration.
Said combination Pm of these parameters, according to one preferred embodiment, can be the sum of the first and of the second basic parameters:

χ1m+χ2m 

Therefore, when using this combined parameter Pm=χ1m+χ2m, the mesh elements of the volume for which Pm is minimum are mesh elements in which the fluorophore may most probably be located.
Another example of a combination of the basic parameters can be the product of the first and second parameter for example:

χ1m·χ2m 

Here again, this combined parameter Pm is minimum for the mesh element(s) in which the fluorophore is most probably located.
Yet another example of a combination of the basic parameters can be the ratio of the third and second basic parameters for example:

χ3/χ2m 

This combined parameter Pm is maximum for the mesh element(s) in which the fluorophore is most probably located, since the third parameter is then minimum, whilst the second is maximum. On the other hand, it will be sought to maximize the inverse combined parameter χ2m/χ3m, if this is the one chosen as combined parameter.
It is also possible to choose only one parameter, a function of one of the above basic parameters, for example the square of one of these parameters, or its inverse.
To form a combined parameter, it is also possible to use at least two of the basic parameters, at least one thereof being the argument of a non-linear function. For example, it is possible to use the product χ1m2.χ2m, the product of the square of the first by the second. Another example is the ratio χ2m2/χ3m.
It is also possible to combine the three above parameters χ1m, χ2m and χ3m by taking into account the fact that the two first have their value increased, for a certain voxel, with the probability of presence of the fluorophore in this voxel, whilst the last one has its value decreased.
In this case of a single fluorophore, each value of the chosen basic parameter or combined parameter, for a given voxel, can be arranged in a vector X (M,1), each coordinate of the vector representing the value Pm of the voxel m (1<i<m). Said vector can be called a fluorescence map.
Irrespective of the parameter chosen, whether a first basic parameter χim (i designating the index of the basic parameter) or a combined parameter Pm, it will be sought to determine whether it has a value higher or lower than a given value Plimite, to determine whether or not the corresponding voxel m is the one in which the fluorophore is probably or even most probably located. Plimite for example is one half of the maximum value assumed by this parameter over all the voxels under consideration and, if the fluorophore is a single fluorophore, it is possible for example to choose the voxel for which the parameter has the maximum value.
As a variant, Plimite is at least equal to one half of the maximum value assumed by this parameter over all the voxels under consideration and, if the fluorophore is a single fluorophore, it is possible for example to choose the voxel for which the parameter has the minimum value.
In relation to the value that the chosen parameter assumes for a given voxel m0, this voxel may appear in an image representing the medium under consideration, for example on the screen 27, with a certain colour scale in relation to the value of Pm0. As a variant, it is only the voxels for which the chosen parameter is greater or lower than the selected threshold value which will be shown in this image.
According to one example of embodiment, the parameter χ3m/χ2m, is used, the ratio of the third to the second.
The example is taken of the configuration in FIG. 2 which is a schematic of a configuration used for tests; this configuration comprises 3×3 sources S1-S9, positioned in plane z=0, and 3×3 detectors D1-D9, positioned in the same plane as the sources. The fluorophore F is positioned in the plane z=1.6 cm.
Illumination is carried out with the 9 sources and a fluorescence signal is formed by the 9 detectors.
The method presented above is applied to this configuration in FIG. 2, only the basic parameter χ3m being chosen, the map of correlations in FIGS. 3A and 3B is obtained.
It can notably be seen in FIG. 3B, that the position is relatively well determined along axes x and y. However, the position along z is not at all precise. FIG. 3A shows slices of FIG. 3B at different altitudes of z: the identified region for location of the fluorophore appears as a lighter region: this region extends along direction z, which is why its trace can be seen at different altitudes z in FIG. 3A (the principle is the same in FIGS. 4A, 5A, 5B, 11C and 12C described below).
FIG. 3C illustrates the value of χ3m in relation to depth, along a vertical straight line (parallel to z) and passing through the fluorophore. The shape of the resulting curve, very flat, confirms this lack of precision: the parameter χ3m has almost the same value everywhere, this value is close to 1, therefore the fluorophore can be everywhere.
The map of criterion −χ2m, with a factor σm taken to be constant and equal to 1 for the needs of the simulation, is illustrated in FIGS. 4A and 4B. More exactly, it is the criterion max(−χ2m)/6−χ2m which is shown, limited to the positive values. FIG. 4A shows slices of FIG. 4B at different altitudes of z.
It can be seen in FIG. 4B that this time a more precise result is obtained along z than along x and y.
FIG. 4C illustrates the value of the combined parameter Pm1=max(−2m)/6−χ2m in relation to depth, along a vertical straight line (parallel to z) and passing through the fluorophore. The shape of the resulting curve, minimum in the vicinity of a particular value (here close to z=1.5 cm), shows that precise location along z can be obtained by using this combined parameter.
It is also possible to use a parameter Pm which combines the two above parameters. For example:

Pm2=χ3m/χ2m,

This gives a precise result both in planes x, y but also in plane z, as illustrated in FIGS. 5A and 5B, which again are successive slices at different altitudes z.
It can be seen that it was therefore possible to position the fluorophore correctly.
Although a description has been given here of the determination of parameters χ1m, χ2m and χ3m, from the zero-order moment or from the first-order normalized moment of at least one distribution Msd(t), it will be appreciated that other basic parameters could be determined, combining firstly a magnitude determined from at least one moment of Msd(t), and secondly an estimation of this magnitude, this estimation being performed by considering that the fluorophore is a single fluorophore and located in voxel m.
In another embodiment, the case is treated of N (>1) fluorophores, N being a known number.
This case is illustrated in FIG. 6, in which a medium 2 contains 2 fluorophores F1, F2 and is provided with two sources S1, S2 and 2 detectors D1, D2 to detect the fluorescence emitted by the fluorophores.
When illumination is produced by the source S1, excitation of the 2 fluorophores occurs, which each re-emit fluorescence radiation captured by the detectors D1, D2.
The same applies when illumination is produced by the source S2: the signal of each of the two detectors contains a contribution of each of the two fluorophores.
The purpose is therefore to be able to identify the positions of the two fluorophores F1 and F2 using the signals produced by the two detectors D1, D2.
The starting points are the intensities and times-of-flight determined from the signals provided by the two detectors.
For each fluorophore, it is then proceeded with computing or estimating one or more basic or combined parameters having the same form as those already described above. But consideration is given to the values of these same parameters which have already been computed for one or more other fluorophores.
In other words, in the formulas which give the values of the basic parameters, the values of moments M(0)sd and M(1)sd for each source-detector pair will be replaced by a corrected value M′(0)sd and M′(1)sd which takes into account the contributions of already identified fluorophores. The corrected values are equal to the difference between the measured moments and the estimated contribution to these moments by already identified or located fluorophores. Procedure is therefore via iteration.
In this variant of the method according to the invention, applied to the case of N fluorophores, a matrix X(M,N) is reconstructed, N being the number of fluorophores, M being the number of voxels. Each element of the matrix Xmn represents the combined coefficient Pm corresponding to the voxel m (1<m<M) and to the fluorophore n (1<n≦N). This matrix is an association of N column vectors X(M,1), i.e. of N fluorescence maps.
As above, the measurements Msd(t) corresponding to different source-detector pairs are performed previously, which allows determination of the moments M(0)sd and M(1)sd for each source-detector pair.
The method is now explained in connection with FIG. 7.
Initially (S0), all the elements of the matrix are initialized to 0. N is known (number of fluorophores present in the diffusing medium based on prior information) and M is also known (number of voxels discretizing the medium).
For each value of n:

    
    
        first the coefficients (1≦m≦M) of column n are set at 0 (step S1);
        the remainder of the measurements is computed (step S2) from the above formulas 10 and 12, discretized onto a mesh whose elements are indexed by m):
    
    






  
    
      
        
          M
          sd
          
            ′
            ⁡
            
              (
              0
              )
            
          
        
        =
        
        ⁢
        
          
            M
            sd
            
              (
              0
              )
            
          
          -
          
            
              S
              
                (
                0
                )
              
            
            ⁢
            
              D
              
                (
                0
                )
              
            
            ⁢
            
              
                ∑
                m
              
              ⁢
              
                
                  ∑
                  n
                
                ⁢
                
                  ·
                  
                    G
                    sm
                  
                  ·
                  
                    X
                    mn
                  
                  ·
                  
                    
                      G
                      md
                    
                    .
                  
                
              
            
          
        
      
    
  
  
    
      
        
          M
          sd
          
            ′
            ⁡
            
              (
              1
              )
            
          
        
        =
        
        ⁢
        
          
            M
            sd
            
              (
              1
              )
            
          
          -
          
            
              S
              
                (
                0
                )
              
            
            ⁢
            
              D
              
                (
                0
                )
              
            
            ⁢
            
              
                ∑
                m
              
              ⁢
              
                
                  ∑
                  n
                
                ⁢
                
                  ·
                  
                    G
                    sm
                  
                  ·
                  
                    X
                    mn
                  
                  ·
                  
                    G
                    md
                  
                  ·
                  
                    (
                    
                      
                        T
                        s
                      
                      +
                      
                        T
                        sm
                      
                      +
                      τ
                      +
                      
                        T
                        md
                      
                      +
                      
                        T
                        d
                      
                    
                    )
                  
                
              
            
          
        
      
    
  




In these formulas Gsm and Gmd are the simulated Green functions (already mentioned above) and Xmn is the element (m,n) of the matrix X, previously initialized to 0, but comprising terms already estimated if n>1; M′(0)sd and M′(1)sd are forced to be strictly positive, by zeroing strictly negative values (step S3). S(0) and D(0) have already been defined.
Each value of M′(0)sd is the difference between the measured value of M(0)sd and the estimated value using the Gsm and Gmd values applied to the current coefficient Xmn, such as results from the preceding iteration (optionally from the initial zeroing).
It can also be said:

    
    
        that M′(0)sd is the measured intensity Msd0 for each selected source-detector pair, corrected by the estimation of the contribution to intensity by already located fluorophores for each of which at least one coefficient of the corresponding column of the matrix X is nonzero;
        whilst M′(1)sd is the mean measured source-detector first-order moment (first-order normalized moment of M(t)), corrected by the estimation of time-of-flight contribution by already located fluorophores for each of which, here again, at least one coefficient of the corresponding column of the matrix X is nonzero.
    
    


It is then possible (step S4), for each voxel m and using M′(0)sd and M′(1)sd, to determine at least one of the basic parameters χim (for example χ1m, χ2m, χ3m) and optionally at least one combined parameter Pm, Pm having the form already defined previously for the case when there is only a single fluorophore.
Therefore, χ1m is determined by replacing M(0)sd by M′(0)sd and χ2m is determined by replacing Tsd by M′(1)sd/M′(0)sd, whilst χ3m is determined by replacing Isd by M′(0)sd 
Therefore n being fixed, for every m (step S5):

X(m,n)=Pm or χim or χ′im(i=1-3).

It is optionally possible to proceed with adjusting the terms of the matrix X (step S6) by assigning a multiplicative coefficient αi to each of the N columns of the matrix X being formed, by endeavouring to cause the following sum to converge towards M(0)sd:




  
    
      ∑
      m
    
    ⁢
    
      
        ∑
        n
      
      ⁢
      
        
          G
          sm
        
        ⁢
        
          X
          mn
        
        ⁢
        
          G
          md
        
        ⁢
        
          α
          n
        
      
    
  
  ≈
  
    M
    sd
    
      (
      0
      )
    
  




Convergence can be obtained using a least square method. This step allows a significant improvement in the quality of results.
Matrix X is adjusted by multiplying each term of a column n by the coefficient αn. It is specified that these coefficients αn can be considered to be the intensities of the fluorescence signal.
Next, n is incremented (step S7) and a return is made to step S1 (in which all the elements of the column n are initialized to 0). If there is no convergence test step S6, S5 is followed by S7.
When the matrix X has been formed (n=N, step S7) by the preceding loop, it can be verified (step S8) that




  
    ∑
    m
  
  ⁢
  
    
      ∑
      n
    
    ⁢
    
      
        G
        sm
      
      ⁢
      
        X
        mn
      
      ⁢
      
        G
        md
      
      ⁢
      
        α
        n
      
    
  




is close to Msd(0).

If the correlation is not satisfactory, the above loop is reiterated from step S2 using X, which has just been computed, as the initial matrix. Correlation, for example, may comprise a so-called residual difference test with respect to data, corresponding to the determination of a difference between an estimated magnitude and an effectively measured magnitude.
Here, the estimated magnitude is:




  
    ∑
    m
  
  ⁢
  
    
      ∑
      n
    
    ⁢
    
      
        G
        sm
      
      ⁢
      
        X
        mn
      
      ⁢
      
        G
        md
      
      ⁢
      
        α
        n
      
    
  




for the measured magnitude Msd(0).
The estimated magnitude is:




  
    
      ∑
      m
    
    ⁢
    
      
        ∑
        n
      
      ⁢
      
        
          G
          sm
        
        ⁢
        
          X
          mn
        
        ⁢
        
          
            G
            md
          
          ⁡
          
            (
            
              Ts
              +
              Tsm
              +
              Tm
              +
              Tmd
              +
              Td
            
            )
          
        
      
    
  
  
    
      ∑
      m
    
    ⁢
    
      
        ∑
        n
      
      ⁢
      
        
          G
          sm
        
        ⁢
        
          X
          mn
        
        ⁢
        
          G
          md
        
        ⁢
        
          α
          n
        
      
    
  




for the measured magnitude Msd(1)/Msd(0) 
Finally, here again, in each column i.e. for each fluorophore, the row m is selected for which the basic parameter χim or combined parameter Pm, gives the minimum or maximum value, depending on the type of this parameter: as already explained above, for each of the parameters χ1m, χ2m the maximum value is taken, and for parameter χ3m the minimum value; these principles are applied to every combined parameter which results from these basic parameters.
The selected row m corresponds to a voxel m in which the fluorophore n is probably or most probably located.
Another embodiment treats the case in which N (>1) fluorophores, N being an unknown integer.
Physically, the situation is substantially the same as the one already described above in connection with FIG. 6 and with the preceding case.
However, the number of fluorophores is unknown, and hence the number of columns which the matrix X must have.
Iteration is also performed, using the same parameters already presented χ1m, χ2m, χ3m and optionally at least one combined parameter Pm.
This variant of the method is now explained with reference to FIG. 8.
N, the number of fluorophores, is initialized to 0 (step S′0), then it is incremented to 1 (step S′1).
Next, (step S′2), all the elements of the matrix are initialized to 0. Only M is known (number of voxels discretizing the medium).
Then Niter, the number of iterations, is initialized to 0 (step S′3), after which it is incremented to 1 (step S′4). During the initialization of Niter, it is also possible to determine the maximum number of iterations it is desired to perform.
Steps S′5-S′11 are identical to steps S2-S8. Reference is therefore to be made to the above description.
If the correlation is satisfactory, the number of fluorophores is equal to the current number N and the matrix X is formed of the previously determined elements (step S′12).
This leads to the end of the method (step S′13).
If the correlation is not satisfactory, it is checked whether the number of iterations is the previously fixed maximum number Niter (step S′11).
If this is not the case, it is returned to step S′3 by incrementing Niter by one unit.
If it is the case, insofar as the correlation test has failed, the number of fluorophores is incremented by one unit (return to step S′1).
The method is continued until a matrix has been determined for which the correlation test is satisfactory (S′11-S′13).
Finally here again, in each column i.e. for each fluorophore, the row m is selected for which the basic or combined parameter gives the minimum or maximum value, depending on the nature of this parameter: as already explained above, for each of the parameters χ1m, χ2m the maximum value is taken and for parameter χ3m the minimum value; and these principles are applied to every combined parameter which results from these basic parameters. The selected row m corresponds to a voxel m in which the fluorophore n is probably or most probably located.
One variant of this method for an unknown number N of fluorophores can be the following:
1) N=1 is used,
2) the reconstruction described above with reference to FIG. 7 is applied for N fluorophores,
3) testing to determine whether data attachment is satisfactory,
4) if yes, end; If not N=N+1, and it is returned to step 2.
An example is now given with two fluorophores.
The method in FIG. 7 is applied.
The two fluorophores are designated F1 and F2 in FIG. 9A. The sources S and detectors D are arranged on a cylindrical surface as can be seen in this figure.
During iteration 1, it can be seen (FIG. 9B) that one of the fluorophores is correctly reconstructed. FIGS. 10A-10D show representations of M(0)exp (measurement of M(0), FIG. 10A), M(1)exp (measurement de M(1), FIG. 10B), M(0)est (computing of M(O), FIG. 10C), M(1)est (computing of M(1), FIG. 10D). These figures show that the convergence of estimations towards measured data is not satisfactory.
During iteration 2, it can be seen (FIG. 9C) that both fluorophores are correctly reconstructed. FIGS. 10E-10H show representations of M(0)exp (measurement of M(0), FIG. 10E), M(1)exp (measurement of M(1), FIG. 10F), M(O)est (computing of M(0), FIG. 10G), M(1)est (computing of M(1), FIG. 10H). These figures show that the convergence of estimations towards measured data is satisfactory.
A fourth embodiment corresponds to a generalization of the first embodiment. In the latter, it was seen that it is possible to determine basic coefficients χim performing a sum for different acquisitions, each acquisition corresponding to a source-detector pair, of a combination between a magnitude determined from at least one moment of Msd(t), and the estimation of this magnitude, this estimation being made by considering a single fluorophore positioned in the voxel m.
Therefore, according to this first embodiment, during the determination of each basic coefficient χim, it was considered that a configuration is assumed a priori in which the fluorophore is single and located in the voxel m.
According to this fourth embodiment, basic coefficients χim are also determined similar to those described in the first embodiment, but when determining the basic coefficients χimN, an a priori configuration is taken in which there are N fluorophores located in N different voxels. The index m then corresponds to a configuration i.e. a particular distribution of these N fluorophores among the M voxels in which the examined medium has been discretized.
It can therefore be understood that there are MN possible configurations (maximum). Therefore it is possible to determine MN n-th basic parameters χimN.
The first embodiment may correspond to a particular case of this fourth embodiment in which N=1: there are therefore M parameters χimN=1 which can be determined, the index m corresponding to a particular distribution of the single fluorophore among the M voxels of the medium: the distribution in which the fluorophore is located in the voxel m.
Therefore, if N fluorophores are considered, N being a natural strictly positive integer, it is possible to determine at least one basic parameter χimN, which, for at least one acquisition Msd(t) corresponding to a source-detector pair, is equal to the combination of a magnitude obtained from at least one moment of Msd(t) with an estimation of said magnitude, this estimation being made by considering that the N fluorophores are distributed over the M voxels of the medium in a given distribution, said distribution being indexed by the parameter m.
Therefore, each i-th basic coefficient χim is obtained by considering a given distribution m of the N fluorophores in the M voxels. As recalled above, there are therefore MN possible distributions, corresponding to as many i-th possible basic coefficients. In other words, 1≦m≦MN 
If it is considered that the magnitude obtained from at least one moment of Msd(t) is a zero-order moment of Msd(t), denoted Msd(0), it is therefore possible to determine a first basic coefficient χ1mN such that:




  
    χ
    
      1
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                
                  (
                  0
                  )
                
              
              -
              
                
                  M
                  
                    sd
                    ,
                    m
                  
                  
                    theo
                    ⁢
                    
                      
                    
                    ⁢
                    0
                  
                
                ⁡
                
                  (
                  
                    
                      α
                      1
                    
                    ,
                    
                      ...
                      ⁢
                      
                        
                      
                      ⁢
                      
                        α
                        N
                      
                    
                  
                  )
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  




where:




  
    
      M
      
        sd
        ,
        m
      
      
        theo
        ⁢
        
          
        
        ⁢
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          sm
          n
        
      
      ⁢
      
        G
        
          
            m
            n
          
          ⁢
          d
        
      
    
  




A configuration m is considered here which corresponds to N fluorophores distributed in the voxels mn, 1≦n≦N.
Gsmn (respectively Gmnd) represent the energy transfer functions between the source and the voxel mn (respectively the voxel mn and the detector d). The coefficients αn, which can be considered as intensities, are obtained by the minimization operation, as is χ1m.
Each element (Msd0−Msd,mtheo0 (α1, . . . αN)) of the sum then represents the difference or more generally the comparison (it could well be a ratio) between:

    
    
        the value of the measured intensity Msd0 for the selected source-detector pair,
        and an estimation of the zero-order moment Msd0, obtained using firstly the Green functions Gsmn and Gmnd which can be estimated in the manner already described above, and secondly a set of coefficients αn each of which is assigned to a fluorophore and which represents the fluorescence emission intensity thereof.
    
    


This first basic parameter χ1mN is minimum for the configurations m the closest to the real distribution of the N fluorophores in the examined medium.
If it is considered that the magnitude obtained from at least one moment of Msd(t) is a first-order normalized moment of Msd(t), denoted Msd(1)/Msd(0), possibly also being denoted Tsd, from which the known temporal magnitudes Ts (mean time of the source), Td (response time of the detector) and σ (duration of fluorescence) may optionally be subtracted, it is then possible to determine a second basic parameter χ2m such that:




  
    χ
    
      2
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
          ′
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                (
                
                  
                    T
                    sd
                  
                  -
                  
                    T
                    s
                  
                  -
                  σ
                  -
                  
                    T
                    d
                  
                
                )
              
              -
              
                (
                
                  
                    T
                    
                      sd
                      ,
                      m
                    
                    theo
                  
                  ⁡
                  
                    (
                    
                      
                        α
                        1
                        ′
                      
                      ,
                      
                        ...
                        ⁢
                        
                          
                        
                        ⁢
                        
                          α
                          n
                          ′
                        
                      
                    
                    )
                  
                
                )
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              T
              sd
            
            )
          
        
      
    
  




where:




  
    
      T
      
        sd
        ,
        m
      
      theo
    
    ⁡
    
      (
      
        
          α
          1
          ′
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
            ′
          
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  




A configuration m is considered here which corresponds to N fluorophores distributed in the voxels mn, 1≦n≦N.
Gsmn(respectively Gmnd) represent the energy transfer functions between the source and the voxel mn (respectively the voxel mn and the detector d). Tthesd (α1, . . . αN) thus defined, effectively corresponds to an estimation of Tsd−Ts−Td−τ.
It can be seen that χ2mN is the difference (or a residual, and more generally a comparison, this comparison possibly also being formed by a ratio) for at least one measurement, between a magnitude formed by the mean measured source-detector time-of-flight (first-order normalized moment of M(t)) preferably but not necessarily reduced by known temporal magnitudes (Ts, Td, τ), and this same magnitude modelled using a contribution of all the fluorophores distributed as per configuration m, via coefficients α′i each of which is assigned to a fluorophore and which represents the fluorescence emission intensity thereof.
The α′n are obtained by the minimization operation, as is χ2mN. The latter, also called the criterion of time-of-flight measurement residual, is also small for configurations m corresponding to the real distribution of the N fluorophores in the medium.
Preferably and as is the case in the above formula, the parameter χ2mN can be normalized by a statistical magnitude relating to the distribution of Tsd. This statistical magnitude acts as confidence indicator, as already explained above for χ2m. Here σ2 (Tsd) has been chosen, but another statistical magnitude (for example standard variation or variance) could be chosen to indicate the confidence assigned to a distribution, but as seen previously, normalization with the coefficient σ2 (Tsd) is preferred.
Other basic parameters obtained from other moments of Msd(t) can be formed. However, the defined basic parameters χ1mN and χ2mN are particularly pertinent.
From one or more formed basic parameters χimN, it is possible to define a combined parameter PmN, combining said basic parameters χimN. By combination, it is recalled that any arithmetic operation is meant whether linear or not, the simplest operations being multiplication, addition or ratios.
According to one preferred embodiment, using the parameters χ1mN and χ2mN previously obtained, it is possible to form a parameter PmN such that PmN=χ1mN+χ2mN. This combined coefficient PmN will be minimal when the distribution of the N fluorophores of index m is the closest to reality.
When PmN is a combination of the basic parameters χ1mN and χ2mN, a set of coefficients αn (1<n<N) or α′n (1<n<N) can be chosen minimizing either χ1mN, or χ2mN, or PmN.
It is then possible to obtain a combined coefficient PmN for all or part of the MN possible distributions of the N fluorophores in the M voxels, a configuration m corresponding to a positioning of the N fluorophores in voxels m1 . . . mN. The distributions are then determined whose combined coefficient PmN is lower than a limit value PmNlimite determined by the operator. This value PmNlimite may for example be equal to the highest obtained value PmN divided by two.
It is also possible to determine a basic parameter χ4mN, such that:




  
    χ
    
      4
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
    
    ⁢
    
      [
      
        
          
            ∑
            sd
          
          ⁢
          
            
              
                (
                
                  
                    M
                    sd
                    
                      (
                      0
                      )
                    
                  
                  -
                  
                    
                      M
                      
                        sd
                        ,
                        m
                      
                      
                        theo
                        ⁢
                        
                          
                        
                        ⁢
                        0
                      
                    
                    ⁡
                    
                      (
                      
                        
                          α
                          1
                        
                        ,
                        
                          ...
                          ⁢
                          
                            
                          
                          ⁢
                          
                            α
                            N
                          
                        
                      
                      )
                    
                  
                
                )
              
              2
            
            
              
                σ
                2
              
              ⁡
              
                (
                
                  M
                  sd
                  0
                
                )
              
            
          
        
        +
        
          
            ∑
            sd
          
          ⁢
          
            
              
                (
                
                  
                    (
                    
                      
                        T
                        sd
                      
                      -
                      
                        T
                        s
                      
                      -
                      σ
                      -
                      
                        T
                        d
                      
                    
                    )
                  
                  -
                  
                    (
                    
                      
                        T
                        
                          sd
                          ,
                          m
                        
                        theo
                      
                      ⁡
                      
                        (
                        
                          
                            α
                            1
                          
                          ,
                          
                            ...
                            ⁢
                            
                              
                            
                            ⁢
                            
                              α
                              N
                            
                          
                        
                        )
                      
                    
                    )
                  
                
                )
              
              2
            
            
              
                σ
                2
              
              ⁡
              
                (
                
                  T
                  sd
                
                )
              
            
          
        
      
      ]
    
  




where:




  
    
      M
      
        sd
        ,
        m
      
      
        theo
        ⁢
        
          
        
        ⁢
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          
            sm
            n
          
          ⁢
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
        
      
    
  




and:




  
    
      T
      
        sd
        ,
        m
      
      theo
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...
          ⁢
          
            
          
          ⁢
          
            α
            N
          
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
          α
          n
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  




This again concerns the computing of a basic parameter χN4,m, which, for at least one of said temporal distributions Msd(t), combines at least one magnitude (the two magnitudes concerned here are M0sd and Tsd) obtained from at least one moment of said distribution, and at least one estimation of these magnitudes (these here being Mtheo0sd,m and Ttheo0sd,m), these estimations being made by considering that there are N fluorophores, each occupying a voxel of the medium according to a distribution m.
It is possible to determine PmN from only one basic parameter χimN, the parameter PmN even possibly being equal to a basic criterion χimN. For example, PmN can be determined so that:

PmN=χN4m 

When N=1, χ1mN=1=χ1m and χ2mN=1= 2m are effectively found, χ1m and χ2m being such as defined in the description of the first embodiment.
The invention is advantageously applied to the locating of tumours with respect to the prostate, but can also be applied to other organs, notably testicles, breast, brain, etc.
Therefore by means of the invention it was possible to reconstruct the fluorescence of capillaries buried at 1 cm or 2 cm.
In this example, by taking into account the accuracy of the temporal measurement, it was possible to reject certain aberrant measurements; since the mean times are obtained by the division of M1 by M0, this ratio becomes most imprecise if there is very little signal.
FIGS. 11A-11C et 12A-12C illustrate the results obtained.
FIG. 11A shows a configuration with a capillary at a depth of about 1 cm, with 9 sources S and 9 detectors D arranged in the plane z=0. The value of τ is 0.95 ns. The criterion χ2m is used here with the approximation Ts=Td=0. The depth found is z=1.1 cm, for x=0.7 cm and y=0 cm.
The graph in FIG. 11B shows that measurements based on intensities alone (curve I) give imprecise results along z, whilst those based on temporal results give precise results along z (curve II), which confirms FIG. 11C.
FIG. 12A shows a configuration with a capillary at a depth of about 2 cm, with 9 sources S′ and 9 detectors D′ arranged in the plane z=0. The value of τ is 0.95 ns. The criterion χ2m is used here with the approximation Ts=Td=0. The depth found is z=1.9 cm, for x=0.2 cm and y=0 cm.
The graph in FIG. 12B shows that measurements based on intensities alone (curve I′) give imprecise results along z, whilst those based on temporal results give precise results along z (curve II′), which confirms FIG. 12C.

## Claims

1. A method for locating at least N (N>1) absorber in a diffusing medium using at least one source of pulsed radiation capable of emitting radiation to excite this absorber, and at least one detector capable of measuring an emission signal emitted by this absorber comprising:
illuminating the medium by a radiation source,
detecting, by at least one detector, the emission signal emitted by this absorber,
for at least one source-detector pair, performing a temporal distribution Msd(t) of the emission signal received by the detector,
the diffusing medium being discretized into M voxels, each absorber occupying one voxel among the M voxels of the medium, according to a distribution m, and
computing at least one basic parameter χNi,m, which, for at least one of said temporal distributions Msd(t), combines at least one magnitude obtained from at least one moment of said distribution, and at least one estimation of this magnitude.

2. The method according to claim 1, also comprising determining a combined parameter PmN, combining at least one or two basic parameters χNi,m.

3. The method according to claim 1, one of the basic parameters χN2,m comprising:
a comparison between the first-order normalized moment of said temporal distribution and the modeling of this first-order normalized moment;
or a comparison between the zero-order moment of said temporal distribution and the modeling of this zero-order moment.

4. The method according claim 1, comprising determining at least one of the following basic parameters:
a first basic parameter χN1,m, which is the sum, for all source-detector pairs, of the differences between the value of the measured intensity Msd0 for each source-detector pair and an estimation of the zero-order moment, for each source-detector pair, obtained using the Green functions Gsmn and Gmnd, for the source and the detector of each pair, this estimation being made by considering that the N absorbers are distributed in the voxels of the medium as per configuration m, in which the absorbers are distributed in the voxels mn; or
a second basic parameter χN2,m which is the sum for all source-detector pairs of the differences, for each source-detector pair, between the mean measured source-detector time-of-flight (first-order normalized moment of M(t)), corrected by known temporal magnitudes relating to the source, to the detector and to the absorbers, and the estimation of this corrected time-of-flight, this estimation being made by considering that the N absorbers are distributed in the voxels of the medium as per configuration m.

5. The method according to claim 4, wherein the basic coefficient χ1mN is such that:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...α
          N
        
      
    
    ⁢
    
      
        ∑
        sd
        
          
        
      
      ⁢
      
        
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                
                  (
                  0
                  )
                
              
              -
              
                
                  M
                  
                    sd
                    ,
                    m
                  
                  
                    theo
                    0
                  
                
                ⁡
                
                  (
                  
                    
                      α
                      1
                    
                    ,
                    
                      ...α
                      N
                    
                  
                  )
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  



where:





  
    
      M
      
        sd
        ,
        m
      
      
        theo
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...α
          N
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          sm
          n
        
      
      ⁢
      
        G
        
          
            m
            n
          
          ⁢
          d
        
      
    
  




a configuration m corresponding to N absorbers distributed in the voxels mn, 1≦n≦N,
Gsmn (respectively Gmnd) representing the energy transfer functions between the source and the voxel mn (respectively the voxel mn and the detector d), the coefficients αn, which can be considered as intensities, being obtained by the minimization operation, as is χ1m;
each element (Msd0−Msd,mtheo0 (α1, . . . αN)) of the sum then representing the difference between:
the value of the measured intensity Msd0 for the selected source-detector pair,
and an estimation of the zero-order moment Msd0, obtained using firstly the Green functions Gsmn and Gmnd, and secondly a set of coefficients αn each of which is assigned to an absorber and which represents the emission intensity thereof,

and the second basic parameter χ2mN possibly being:





  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
          ′
        
        ,
        
          ...α
          N
          ′
        
      
    
    ⁢
    
      
        ∑
        sd
        
          
        
      
      ⁢
      
        
      
      ⁢
      
        
          
            (
            
              
                (
                
                  
                    T
                    sd
                  
                  -
                  
                    T
                    s
                  
                  -
                  
                    T
                    d
                  
                
                )
              
              -
              
                (
                
                  
                    T
                    
                      sd
                      ,
                      m
                    
                    theo
                  
                  ⁡
                  
                    (
                    
                      
                        α
                        1
                        ′
                      
                      ,
                      
                        ...α
                        N
                        ′
                      
                    
                    )
                  
                
                )
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              T
              sd
            
            )
          
        
      
    
  



where:





  
    
      T
      
        sd
        ,
        m
      
      theo
    
    ⁡
    
      (
      
        
          α
          1
          ′
        
        ,
        
          ...α
          N
          ′
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  



Ts, Td respectively representing the response times of the source and of the detector, Tsmn and Tmnd representing the respective times-of-flight between the source and the absorber located in voxel mn, and between the absorber located in voxel mn and the detector, the coefficient σ2 (Tsd) corresponding to an estimation of the distribution variance Tsd, Ttheosd (α′1 . . . α′N) then representing an estimation of the first-order normalized moment of function Msd(t), the coefficients α′N being obtained by the minimization operation.

6. The method according to claim 1, wherein a combined parameter PmN, is also determined equal to the sum or product of the basic coefficients χN1,m, and χN2,m, the combined parameters PmN of lowest value then corresponding to the distributions of the absorbers the closest to the effective distribution of the absorbers.

7. The method according to claim 1, the number of absorbers being equal to 1, each basic parameter χN=11,m, and χN=12,m, then being determined by considering that the absorber is single and located in the voxel m.

8. The method according to claim 7, wherein χN=11,m, and χN=12,m are determined according to the following equations:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      min
      
        α
        m
      
    
    ⁢
    
      
        ∑
        sd
        
          
        
      
      ⁢
      
        
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                0
              
              -
              
                
                  
                    α
                    m
                  
                  ·
                  
                    G
                    sm
                  
                
                ⁢
                
                  G
                  md
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  



and





  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      ∑
      
        s
        ,
        d
      
      
        
      
    
    ⁢
    
      
    
    ⁢
    
      
        
          (
          
            
              (
              
                
                  T
                  sd
                
                -
                
                  T
                  s
                
                -
                
                  T
                  d
                
              
              )
            
            -
            
              (
              
                
                  T
                  sm
                
                +
                
                  T
                  
                    d
                    ⁢
                    
                      
                    
                    ⁢
                    m
                  
                
              
              )
            
          
          )
        
        2
      
      
        
          σ
          2
        
        ⁡
        
          (
          
            T
            sd
          
          )
        
      
    
  



Msd0 representing the zero-order moment of the function Msd(t),
Gsm and Gmd being energy transfer functions, for each source-detector pair sd, between the source s and the voxel m and between the voxel m and the detector d respectively,
Tsm and Tdm representing the respective times-of-flight between the source and the voxel m and the voxel m and the detector d,
Ts, Td respectively representing the response times of the source and of the detector,
the coefficient σ2(Tsd) corresponding to an estimation of distribution variance Tsd.

9. The method according to claim 7, wherein a combined parameter PmN is also determined equal to the sum or to the product of the basic coefficients χN=11,m, and χN=12,m, the combined parameters PmN=1 of lowest value then corresponding to the distributions of absorbers the closest to reality.

10. The method according to claim 7, to locate several absorbers, further comprising at least one of the following steps:
an adjustment step (S6, S′9) adjusting all the intensity contributions of already located absorbers to the measured intensity;
a computing step (S8, S′11) computing a correlation between the measured intensity and all the intensity contributions of all the absorbers; and
a step (S′1) incrementing the number of absorbers if the computing of a correlation is not satisfactory.

11. The method according to claim 1, wherein the radiation sources and the detectors have reflection geometry.

12. A device for locating at least N (N>1) absorbers, in a diffusing medium, comprising at least one pulsed radiation source capable of emitting radiation to excite this absorber and at least one detector capable of measuring a signal emitted by this absorber, comprising:
means, for at least one source-detector pair, for performing temporal distribution Msd(t) of the signal received by the detector,
means for producing meshing of the volume into mesh elements m or voxels, each of the N absorbers being distributed in one of the M voxels of the medium as per a distribution m,
means for computing at least one basic parameter χNi,m, which, for at least one of said temporal distributions Msd(t), combines at least one magnitude obtained from at least one moment of said distribution, and at least one estimation of this magnitude.

13. The device according to claim 12, further comprising means for determining a combined parameter PmN, combining at least one or two basic parameters χNi,m.

14. Device according to claim 12, the means to compute at least one parameter comprising means for determining:
a parameter χN2,m comprising the comparison between the first-order normalized moment of said temporal distribution and the modeling of this first-order normalized moment;
or a parameter χN1,m comprising the comparison between the zero-order moment of said temporal distribution and the modeling of this zero-order moment.

15. The device according to claim 12, the means for computing at least one parameter comprising means for:
determining a first basic parameter χN1,m, which is the sum, for all source-detector pairs, of the differences between the value of the measured intensity Msd0 for each source-detector pair, and an estimation of the zero-order moment for each source-detector pair, obtained using the Green functions Gsmn and Gmnd, for the source and the detector of each pair, this estimation being made by considering that the N absorbers are distributed in the voxels of the medium as per configuration m, in which the absorbers are distributed in the voxels mn; and
determining a second basic parameter χN2,m, which is the sum, for all source-detector pairs, of the differences, for each source-detector pair, between the mean measured source-detector time-of-flight (first-order normalized moment of M(t)), corrected by known temporal magnitudes relating to the source, to the detector and to the absorber, and the estimation of this corrected time-of-flight, this estimation being made by considering that the N absorbers are distributed in the voxels of the medium as per configuration m.

16. The device according to claim 15,
wherein the basic coefficient χ1mN is such that:





  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
        
        ,
        
          ...α
          N
        
      
    
    ⁢
    
      
        ∑
        sd
        
          
        
      
      ⁢
      
        
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                
                  (
                  0
                  )
                
              
              -
              
                
                  M
                  sd
                  
                    theo
                    
                      
                    
                    0
                  
                
                ⁡
                
                  (
                  
                    
                      α
                      1
                    
                    ,
                    
                      ...α
                      N
                    
                  
                  )
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  



where:





  
    
      M
      sd
      
        theo
        0
      
    
    ⁡
    
      (
      
        
          α
          1
        
        ,
        
          ...α
          N
        
      
      )
    
  
  =
  
    
      ∑
      
        n
        =
        1
      
      N
    
    ⁢
    
      
    
    ⁢
    
      
        α
        n
      
      ⁢
      
        G
        
          sm
          n
        
      
      ⁢
      
        G
        
          
            m
            n
          
          ⁢
          d
        
      
    
  




a configuration m corresponding to N absorbers distributed in the voxels mn, 1<n<N,
Gsmn (respectively Gmnd) representing the energy transfer functions between the source and the voxel mn (respectively the voxel mn and the detector d), the coefficients αn being obtained by the minimization operation, as is χ1m;
each element (Msd0−Msdtheo0 (α1 . . . αN)) Of the sum then representing the difference between:
the value of the measured intensity Msd0 for the selected (source-detector) pair,
and an estimation of the zero-order moment Msd0, obtained using firstly the Green functions Green Gsmn and Gmnd, and secondly a set of coefficients αn each of which is assigned to an absorber and which represents the emission intensity thereof,

and the second basic parameter χ2mN possibly being such that:





  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
    N
  
  =
  
    
      min
      
        
          α
          1
          ′
        
        ,
        
          ...α
          N
          ′
        
      
    
    ⁢
    
      
        ∑
        sd
        
          
        
      
      ⁢
      
        
      
      ⁢
      
        
          
            (
            
              
                (
                
                  
                    T
                    sd
                  
                  -
                  
                    T
                    s
                  
                  -
                  
                    T
                    d
                  
                
                )
              
              -
              
                (
                
                  
                    T
                    sd
                    theo
                  
                  ⁡
                  
                    (
                    
                      
                        α
                        1
                        ′
                      
                      ,
                      
                        ...α
                        N
                        ′
                      
                    
                    )
                  
                
                )
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              T
              sd
            
            )
          
        
      
    
  



where:





  
    
      T
      sd
      theo
    
    ⁡
    
      (
      
        
          α
          1
          ′
        
        ,
        
          ...α
          N
          ′
        
      
      )
    
  
  =
  
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          
            G
            
              
                m
                n
              
              ⁢
              d
            
          
          ⁡
          
            (
            
              
                T
                
                  sm
                  n
                
              
              +
              
                T
                
                  
                    m
                    n
                  
                  ⁢
                  d
                
              
            
            )
          
        
      
    
    
      
        ∑
        
          n
          =
          1
        
        N
      
      ⁢
      
        
      
      ⁢
      
        
          α
          n
          ′
        
        ⁢
        
          G
          
            sm
            n
          
        
        ⁢
        
          G
          
            
              m
              n
            
            ⁢
            d
          
        
      
    
  



Ts, Td respectively representing the response times of the source and of the detector, Tsmn and Tmnd representing the respective times-of-flight between the source and the absorber located in the voxel mn, and between the absorber located in the voxel mn and the detector, the coefficient σ2(Tsd) corresponding to an estimation of distribution variance Tsd, Ttheosd (α′1 . . . α′N) then representing an estimation of the first-order normalized moment of function Msd(t), the coefficients α′N being obtained by the minimization operation.

17. The device according to claim 15, comprising means for determining a combined parameter PmN, equal to the sum or to the product of the basic coefficients χN1,m, and χN2,m.

18. The device according to claim 12, the number of fluorophores being equal to 1, each basic parameter χN=11,m, and χN=12,m, then being determined by considering that the absorber is single and located in the voxel m.

19. The device according to claim 18, wherein χN=11,m, and χN=12,m are determined using the following equations:




  
    χ
    
      1
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      min
      
        α
        m
      
    
    ⁢
    
      
        ∑
        sd
      
      ⁢
      
        
          
            (
            
              
                M
                sd
                0
              
              -
              
                
                  
                    α
                    m
                  
                  ·
                  
                    G
                    sm
                  
                
                ⁢
                
                  G
                  md
                
              
            
            )
          
          2
        
        
          
            σ
            2
          
          ⁡
          
            (
            
              M
              sd
              0
            
            )
          
        
      
    
  



and





  
    χ
    
      2
      ⁢
      
        
      
      ⁢
      m
    
  
  =
  
    
      ∑
      
        s
        ,
        d
      
    
    ⁢
    
      
        
          (
          
            
              (
              
                
                  T
                  sd
                
                -
                
                  T
                  s
                
                -
                
                  T
                  d
                
              
              )
            
            -
            
              (
              
                
                  T
                  sm
                
                +
                
                  T
                  dm
                
              
              )
            
          
          )
        
        2
      
      
        
          σ
          2
        
        ⁡
        
          (
          
            T
            sd
          
          )
        
      
    
  



Msd0 representing the zero-order moment of function Msd(t),
Gsm and Gmd being energy transfer functions, for each source-detector pair sd, between the source s and the voxel m and between the voxel m and the detector d respectively,
Tsm and Tdm representing the respective times-of-flight between the source and the voxel m and the voxel m and the detector d,
Ts, Td respectively representing the response times of the source and of the detector,
the coefficient σ2(Tsd) corresponding to an estimation of distribution variance Tsd.

20. The device according to claim 12, to locate several absorbers, further comprising, for cases in which there are several absorbers,
means for adjusting all the intensity contributions of already located absorbers to the measured intensity;
means for computing a correlation between the measured intensity and all intensity contributions of all the absorbers; and
if the number of the several absorbers is unknown, means for incrementing the number of absorbers if the computing of a correlation is not satisfactory.