METHOD FOR DETERMINING FUNCTIONAL VOLUMES FOR DETERMINING BIOKINETICS

Abstract

A method for determining functional volumes for the search for kinetics representing the trend of concentration of a radioactive tracer in an area of biological tissue, being applied to spatial components comprises the following steps applied iteratively according to a Markov chain Monte-Carlo scheme: generation of a set K0λ made up of a set of candidate kinetics associated with values of probability of appearance, depending on the concentration λ of radioactive tracer; a labeling step during which, for each spatial component of index K, the probabilities of selection of the kinetics of the set K0λ are weighted by introducing a function λk representative of the concentration of radioactive tracer in this component to obtain a set of indicative values, an indicative value Dk designating the kinetic with which the spatial component K is associated; construction of functional volumes, a functional volume VFj made up of the set of spatial components which share the same indicative value Dk.

Description

The invention relates to a method for determining functional volumes for positron emission tomography and applies notably to the fields of positron emission tomography and pharmacokinetics.

Positron emission tomography, usually referred to by the acronym PET, is a nuclear medicine imaging technology which uses molecules marked with a positron emitting radionucleide to image, in vivo, the molecular interactions of biological processes with a minimum of perturbations. PET is the modality of reference for molecular imaging because of its molecular sensitivity of the order of the femto-mole and of the quantitative character of the measurements that it makes possible. PET imaging is used in particular to determine the functional volumes hereinafter in the description designated by the acronym VF. A functional volume VF corresponds to an area consisting of tissues exhibiting identical metabolic properties.

The state of the art in the field of estimating functional volumes proceeds indirectly by reconstructing, in a first step, a temporal series of discretized images then by estimating the functional volumes from this series of images. Now, reducing the dose of radioactive tracer, also called radiopharmaceutical, contributes to degrading the signal-to-noise ratio of each of the images to the detriment, in this indirect approach, of the determination of the functional volumes.

PET is also often used in clinical procedures in a static, 2D or 3D spatial context, in which the time variable is disregarded, for want of tools allowing for an analysis of the dynamic data that is robust, automatic and validated. This use is significantly sub-optimal. In effect, PET molecular imaging is by nature dynamic inasmuch as the distribution of the radioactive tracer injected into the organism evolves over time as a function of the molecular interactions between the tracer and the tissues. The objective of the interpretation of the examination can then be summarized as a quantitative and objective characterization of the spatial extensions exhibiting similar temporal pharmacokinetic behaviors, hence the notion of functional volumes. It is quantitative inasmuch as it has to be pronounced on the uncertainty associated with the molecular concentration contained in each functional volume in terms of quantity of radioactive tracer. It is objective inasmuch as the number, the positioning and the extension of the functional volumes do not rely on the supervision of an expert having first defined one or more regions of interest.

The state of the art in the dynamic study of biokinetics by PET imaging is based on a temporal segmentation of the information in order to obtain a 3D series. This means that the acquisition is subdivided into time slices reconstructed independently of one another. Although it allows for a dynamic tracking of the distribution of the radioactive tracer from the series of reconstructed images, this approach presents the drawback that only a proportion of the observations is used to reconstruct an image. This statistical loss fatally induces a deterioration of the quality of the reconstruction which is all the greater when the time slice is short.

Most of the time, the determination of the spatial extension of the functional volume is defined a priori, from the anatomical imaging. Then, from the temporal series of images, the concentration of mean activity inside the functional volume, which in reality corresponds to an anatomical volume, is measured, and measured for each time slice. Existing methods make it possible to determine the extension of the functional volume directly from the dynamic images. Source separating techniques are also used to determine the distinct molecular processes from the reconstructed images. However, the number of kinetics is fixed in these approaches.

A. Reader et al. take up this principle in their paper entitled Joint estimation of dynamic PET images and temporal basis functions using fully 4D ML-EM, Physics in Medicine and Biology, 51(21):5455-5474, 2006, by using a dictionary of temporal functions to seek the breakdown in a functional base which best adjusts the dataset. From a limited number of temporal basis functions, typically about ten or so, it is possible a priori to describe, by their combination, an infinite number of time activity curves, usually referred to by the acronym TAC. These 4D reconstruction methods make it possible to enhance the signal-to-noise ratio of the reconstructed images. All these techniques seek to perform a classification, also referred to by the term “clustering”, of space-time data of the disintegrations, but proceed indirectly to this end since a step of reconstructions of voxelized images over the entire space is necessary. In effect, these techniques require the discretization of all the spatial and temporal dimensions, resulting in the manipulation and storage of large volumes of data and considerable computation times. Also, the latter approach does not make it possible to quantify the estimation error without approximation.

Works on reconstruction in emission tomography without recourse to the discretization of the space have been conducted in a non-parametric Bayesian framework. This approach has been extended to the space-time reconstruction as presented in the paper by D. Fall, E. Barat, C. Comtat, T. Dautremer, T. Montagu, and S. Stute entitled Continuous space-time reconstruction in 4D PET, IEEE Medical Imaging Conference, 2011. This paper is called Fall, Barat et al. hereinafter in the description. This method addresses the problem of the quantitative determination of the reconstruction error of the functional volumes and the objective of reducing the dose injected into the patient. It further proposes a space-time “clustering” in a continuous space directly from the observations. However, in this approach, the TACs are all considered as normalized kinetics. This modeling, although it allows for the automatic detection of any number of kinetics, is not optimal inasmuch as the discrimination between kinetics uses only their forms and does not take into account the concept of concentration, that is to say the amplitude of the TACs expressed in terms of Becquerel per second, this unit being denoted Bq.s⁻¹. The limitation of this modeling comes from the model retained for the clustering of the spatial components with a particular kinetic. In effect, the latter implements a hierarchical model in which the random distribution of the kinetics, represented by K₀, is itself distributed according to a Pólya tree Dirichlet process. Then, the distribution of the space-time components denoted H is chosen as a Dirichlet process whose central measurement involves K₀. No information concerning the concentration is present in this stochastic model.

One aim of the invention is notably to mitigate the abovementioned drawbacks.

To this end, the subject of the invention is a method for determining functional volumes for the search for kinetics representing the trend of the concentration of a radioactive tracer in an area of biological tissue. The method is applied to spatial components and comprises the following steps applied iteratively according to a Markov chain Monte-Carlo scheme:

• • a step of generation of a set K₀^λ made up of a set of candidate kinetics associated with values of probability of appearance of these kinetics, these values depending on the concentration λ of radioactive tracer;
  • a labeling step during which, for each spatial component of index K, the probabilities of selection of the kinetics of the set K₀^λ are weighted by introducing therein a function λ_k representative of the concentration of radioactive tracer in this component so as to obtain a set of indicative values, an indicative value D_k designating the kinetic with which the spatial component K is associated;
  • a step of construction of functional volumes, a functional volume VF_j being made up of the set of the spatial components which share the same indicative value D_k.

According to one aspect of the invention, a step computes averages of the functional volumes VF_j obtained during the iterative process.

In one embodiment, a step determines statistical estimators concerning the functional volumes VF_j obtained during the iterative process so as to quantify the estimation uncertainty. The variance of the functional volumes VF_j obtained during the iterative process can then be determined. A credible interval can also be determined.

By way of example, the set K₀^λ is obtained by the determination of two sets of scalar variables V_j and Γ_j, these variables relating to kinetic components of index j, in which V_j represents a fixed weighting coefficient independent of the concentration λ and Γ_j represents a reference concentration value associated with the kinetic of index j.

According to one aspect of the invention, V_j is determined by using the expression:

$V_jD,T\stackrel{\mathrm{ind}}{~}\mathrm{Beta}(1+m_j,v+\text{?}m_l\text{?}\text{indicates text missing or illegible when filed}$

an expression in which:

• • for any j, m_j designates the number of components of the distribution H allocated to the j^th component of K₀^Λ, H being a random measurement distributed according to a partially dependent Dirichlet process, K₀^Λ being an innumerable collection of random probability measurements which accepts, as distribution, a local Dirichlet process;
  • J_kn designates the greatest index of the components of K₀^Λ having at least one component of H allocated;
  • Beta(,) represents the beta probability law.

In one embodiment, Γ_j is determined by using the expression:

${\Gamma }_jV_j,D,w,\tilde{Z}\stackrel{\mathrm{ind}}{~}U_j^{*}\left({\Gamma }_j\right)$ $\mathrm{with}$ $U_j^{★}\left({\Gamma }_j\right)\propto \left(a_{j,\Gamma }^{\prime }<{\Gamma }_j<b_{j,\Gamma }^{\prime }\right)\times \prod _{k\in {}_j^{+}}^{}\left(1-V_j\left(d\left({\lambda }_k,{\Gamma }_j\right)<\psi \right)\right)$

expressions in which:

• • Γ_j designates the parameter characteristic of the concentration associated with the kinetic j;
  • D designates the set of the variables D_k which indicate with which kinetic of K₀^λ the spatial component K is associated;
  • D_j={k≦k_n:D_k=j} represents the set of the indices K of the spatial components associated with the kinetic j;
  • D_j⁺={k≦k_n:D_k>j} represents the set of the indices K of the spatial components associated with a kinetic of index greater than j;
  • w designates the set of the weightings w_k of the spatial components k in the Dirichlet mixture H;
  • {tilde over (Z)} designates the set of the parameters, average and covariance matrix,
  • {tilde over (Z)}_k associated with the spatial component K of the Dirichlet mixture H;
  • here indicates an independent random generation for each parameter Γ_j;
  • ∞ here signifies “proportional to”;
  • ( ) represents the “gate” function such that (condition)=1 if condition is true, and 0 otherwise;
  • V_j designates a coefficient associated with each kinetic involved in the weighting of said kinetic in K₀^λ;
  • d( ) represents the Euclidian distance in : d(x₁,x₂)=x₁ x₂| in which x₁, x₂∈;
  • ψ is a threshold chosen a priori which makes it possible to parameterize K₀^λ.

According to one aspect of the invention, the labeling step is performed:

• • by independently generating for any k≦k_n:

$\left(D_kV,v,Q^{★},\Gamma ,w,\tilde{Z},C,T\right)\stackrel{\mathrm{ind}}{~}\sum _{j=1}^{J^{★}}p_{j,k}{\delta }_j(·)$ $\mathrm{with}$ $p_{j,k}\propto \left(p_j\left({\lambda }_k\right)>{\upsilon }_k\right)\max \left(p_j\left({\lambda }_k\right),\text{?}\right)\text{?}f_{Q_j^{★}}\left(T_i\right)$ $\mathrm{and}$ $\sum _{j=1}^{J^{★}}p_{j,k}=1$ $\text{?}\text{indicates text missing or illegible when filed}$

• • by independently generating for any K, such that k_n<k≦k*:

$\left(D_kV,v,Q^{★},\Gamma ,w,\tilde{Z}\right)\stackrel{\mathrm{ind}}{~}\sum _{j=1}^{J^{★}}p_{j,k}{\delta }_j(·)$ $\mathrm{with}$ $p_{j,k}\propto \left(p_j\left({\lambda }_k\right)>{\upsilon }_k\right)\max \left(p_j\left({\lambda }_k\right),\text{?}\right)$ $\mathrm{and}$ $\sum _{j=1}^{J^{★}}p_{j,k}=1.$ $\text{?}\text{indicates text missing or illegible when filed}$

in these expressions,

• • k_n designates the number of spatial components with which at least one coincidence event is associated;
  • k* designates the total number of spatial components of the Dirichlet mixture H for a given iteration;
  • V designates the set of the coefficients V_j involved in the weightings of the kinetics mixture K₀^λ;
  • v represents the set of the auxiliary variables v_k, used in the representation of the Dirichlet kinetics process;
  • Q* designates the set of the Pólya trees Q*_j characterizing the kinetics;
  • Γ designates the set of the parameters Γ_j characteristic of the concentration associated with the kinetic j;
  • C designates the set of the classification variables C_i which associate each coincidence event i with a spatial component of the Dirichlet mixture H;
  • T represents the set of the times of occurrences T_i of the coincidence events;
  • p_j,k designates the probability of associating the component K of the Dirichlet mixture H with the kinetic j;
  • p_j(λ) represents the weight function taken at λ associated with the probability of the kinetic j of K₀^λ;
  • δ_j( ) represents the Dirac measurement function such that δ_j(j′)=1 if j=j′ and is 0 otherwise;
  • here indicates an independent random generation for each classification variable D_k;
  • max( ) represents the “maximum of” function;
  • ζ designates is a threshold parameter chosen arbitrarily and involved in the implementation of the “slice sampling” algorithm of the local Dirichlet kinetics process;
  • {i: C_i=k} designates the set of the coincidence events i which are associated with the spatial component K;
  • f_Q*j represents the kinetic j, i.e. the density of the Pólya tree j of the local Dirichlet mixture of kinetics;
  • J* designates the total number of components of the local Dirichlet mixture of kinetics for a given iteration.

According to another aspect of the invention, the concentration λ_k is obtained by the application of a function φ_λ defined such that λ_k=φ_λ({tilde over (Z)}_k, S) with:

${\lambda }_k=\text{?}\left(\sum _{m=1}^{\infty }\text{?}\left(\text{?}\text{?}\right)\right)$ $\mathrm{with}$ $\text{?}~\left(\text{?}\text{?}\right)$ $\text{?}\text{indicates text missing or illegible when filed}$

• • in which:
  • _{{tilde over (x)} k} designates the mathematical expectation relative to the Gaussian probability law of {tilde over (x)}_k, N({tilde over (x)}_k|{tilde over (Z)}_k);
  • N( ) represents the Gaussian probability law.

Also the subject of the invention is a positron emission tomography device implementing the method described previously.

The method according to the invention offers the decisive advantage of increasing the distances between functional volumes and thereby increasing the separation quality.

Also, the method according to the invention lies in a novel approach that makes it possible to shrewdly establish the a posteriori law of the concentration of each kinetic, thus conferring on the method the possibility of providing an estimation of the uncertainty on the asymptotic metabolization parameters of a tissue.

Advantageously, the method according to the invention makes it possible to substantially reduce the dose injected into the patient for a four-dimensional PET imaging while providing functional information to a molecular scale.

A repeated use of PET imaging during therapeutic monitoring can then be envisaged without any dosimetric prejudice to the patient. A wider use in pediatrics and in translational research can also be considered.

Furthermore, a reduction of the dose of radioactive tracer injected automatically brings about a reduction of the exposure to ionizing radiations of the medical personnel involved in preparing and carrying out the examination and of the entourage of the patient.

Other features and advantages of the invention will become apparent from the following description given by way of illustration and as a nonlimiting example, in light of the attached drawings in which:

FIG. 1 schematically illustrates the various steps of the method for determining functional volumes;

FIG. 2 illustrates how the weights and auxiliary variables of the distribution H of the space-time components are generated;

FIG. 3 schematically illustrates how the spatial parameters used by the method are generated;

FIG. 4 schematically illustrates one way of determining the local and auxiliary variables of the set K₀^λ;

FIG. 5 illustrates how the kinetics are allocated to the components of the set K₀^λ;

FIG. 6 shows how the places of emission are determined and how the components of H are determined;

FIG. 7 represents one way of computing the concentration λ_k of radioactive tracer;

FIG. 1 schematically illustrates the various steps of the method for determining functional volumes.

The aim of the method is to determine functional volumes for the search for kinetics. It is applied to measurements presented in the form of coincidence events 106 and comprises a series of steps applied iteratively 111 according to a Markov chain Monte-Carlo scheme.

In a first step 108, a set K₀^λ is generated made up of a set of candidate kinetics associated with values of probability of appearance of these kinetics, these values depending on the concentration λ of radioactive tracer.

In a so-called labeling second step 109, for each spatial component of index k 106, the probabilities of selection of the kinetics of the set K_oλ are weighted 102 by introducing therein a function λ_k. This function is representative of the concentration of radioactive tracer of this component so as to obtain a set of indicative values 103. An indicative value D_k designating the kinetic with which the spatial component k is associated is determined.

In a third step 110, the functional volumes are constructed, a functional volume VF_j, being made up of the set of the spatial components which share the same indicative value D_k.

The method for determining functional volumes described hereinbelow makes it possible to mitigate the independence of the model used by Fall, Barat et al. with regard to the concentration information. To achieve this, the concentration of radioactive tracer is estimated at any point of space by marginalization, that is to say by integration over time of the distribution of space-time activity. This estimated concentration is then considered as a random co-variable endogenous to the problem. The term endogenous underscores the fact that this co-variable is not external because the latter is determined by beginning the reconstruction. Provided with this co-variable denoted λ, the proposed method implements a modeling in the distribution of the kinetics dependent on the concentration of radioactive tracer. This amounts to replacing the distribution K₀ with a collection of distributions K₀^Λ={K₀^λ: λ∈Λ}. This collection of distributions then accepts an independent Dirichlet process as distribution.

In one embodiment, a local Pólya tree Dirichlet process IDP is chosen. One consequence is that the distribution H of the space-time components is no longer generated according to a Dirichlet process but according to a distribution referred to hereinbelow in the description as partially dependent Dirichlet process and designated by the acronym pdDP.

According to one aspect of the invention, the construction of the pdDP is parameterized by a function φ_λ which makes it possible to link the concentration co-variable λ to the parameters of the marginal spatial distribution obtained by integration over time of the space-time distribution.

The general problem of reconstruction of a probability distribution G(·) for which only samples of a projection thereof denoted F(·) can be observed is formulated in the non-parametric Bayesian framework for the inverse problems by using the following expressions:

G˜

F(·, t)=∫_xP(·|x)G(dx, t)

Y _i , T _i |F F, for i=1, . . . , n   (1)

in which:

• • n is the number of samples observed;
  • G(·) is a random probability distribution distributed according to , defined over (X×T, σ(X)σ(T));
  • the distribution G(·) is estimated from observations distributed according to F(·): (Y, T)′={(Y₁, T₁), . . . , (Y_n, T_n)};
  • P(·|x) is the assumed known probability distribution, called projector, indexed by x and defined over (y, σ(y));
  • the notation σ(ε) corresponds to the sigma-algebra over ε.

The object that is intended to be reconstructed is the spatial distribution of activity in the field of view of the tomography. In PET, each coincidence event is assigned to a virtual line LOR, LOR being an acronym for “line-of-response”, joining the centers of the two detectors. A response line LOR of index l is characterized by the parameters defining its position in the space of the detectors. A correlation is constructed below between the parameters of the LORs and the index l such that the variable y_i corresponds to the index of the LOR having recorded the i^th observed coincidence. More specifically, the observation y_i designates the coordinates of the LOR which joins the two detectors having recorded the coinciding photons deriving from the annihilation of the positrons having taken place at x_i.

In this context, the terms of the equation (1) correspond to the following interpretation:

• • X⊂³ corresponds to the region of the reconstruction space;
  • T(⁺ corresponds to the acquisition time band;
  • Y_i corresponds to the index of the observed LOR;
  • T_i corresponds to the time of occurrence of the i^th event observed;
  • P(y=m|x) for m=1, . . . , M is the distribution corresponding to the projector, that is to say to the probability of recording an event on the m^th detection LOR for a localized emission at x, that is to say the place of annihilation for any x∈X,

$\sum _{m=1}^M\left(y=mx\right)={\int }_y^{}\left(yx\right)=1$

• • G(·) is the space-time distribution of the emissions corresponding to the recorded events. The set of observations is truncated because of the limited axial field of the imager.

In order to take account of these truncations in the reconstruction, an LOR of zero index is defined, the latter corresponding to the non-observable emissions such that y*=y∪{0}. It is then possible to define the probability of projection P*(y−m|x) over (y*, σ(y*)) and

$\left(yx\right)=\frac{{}^{*}\left(yx\right)}{{\int }_y^{}{}^{*}(yx)}=\frac{{}^{*}\left(yx\right)}{\sum _{m=1}^M{}^{*}\left(y-mx\right)}$

The space-time density of the places of emissions can then be written:

$G^{★}\left(x,t\right)\propto \frac{G\left(x,t\right)}{{\int }_y^{}{}^{*}(yx)}$

an expression in which ∫_yP*(dy|x) designates the probability of recording a localized emission at x by taking into account the geometry of the system and the physical properties of the detection system.

The expression (2) formally defines the proposed model which makes it possible to take into account the spatial concentration as endogenous co-variable of the reconstruction problem:

Y _i |X _i P(Y_i|X_i)

X _i , T _i |Z _i , Q _i N(X_i|Z_i)×Q_i(T_i)

Z _i, Q_i|HH

H˜pdDP(θ, G₀, K₀⁻, φ_x)

K₀^Λ={K₀^λ: λ∈Λ}˜IDP(v, PT_L(A, Q₀), γ₀, ψ)   (2)

In this expression,

• • H is a random measurement distributed according to a partially dependent Dirichlet process pdDP defined as

$\begin{array}{cc}H(·)=\sum _{k=1}^{\infty }w_k\text{?}\text{?}(·)\text{?}\text{indicates text missing or illegible when filed}& \left(3\right)\end{array}$

• • δ_ξ(·) the Dirac measurement at ξ.

N(m, Σ) designates multi-variate normal-Gaussian law in dimension p=3 parameterized by its mean m and its covariance matrix Σ, of density:

$f_{}\left(x\right)=\frac{1}{\left(2\pi \right)\text{?}{\sum }^{\frac{1}{2}}}\exp \left(-\frac{1}{2}{\left(x-m\right)}^T\sum _{}^{-1}\left(x-m\right)\right)$ $\text{?}\text{indicates text missing or illegible when filed}$

The generative model defined in the expression (2) and the concept of pdDP of the expression (3) are explained hereinbelow.

The expression of H in (3) is a statistical mixture in which w_k corresponds to the weight of the K^th component, {tilde over (Z)}_k to the spatial parameters and {tilde over (Q)}_k to the kinetic, that is to say the temporal probability distribution, of the K^th space-time component of the mixture. The distributions for the w_k and {tilde over (Z)}_k are identical to those of a Dirichlet process.

w=w₁, w₂, . . . ˜GEM(θ) applies, as indicated in the paper by J. Pitman entitled “Some developments of the Blackwell-MacQueen urn scheme”, T. F. et al., editor, Statistics Probability and Game Theory, Papers in honor of David Blackwell, volume 30 of Lecture Notes-Monograph Series, pages 245-267, Institute of Mathematical Statistics, 1996, such that:

• • w₁=v₁;
  • for any k≧2, w_k=v_k·Π_i−1^k−1(1 v_i) in which v₁, v₂, . . . are identically distributed according to a law Beta (1, θ), of which the can be obtained by using the following expression:

$f_{\mathrm{Beta}}\left(\upsilon a,b\right)=\frac{\Gamma \left(a+b\right)}{\Gamma \left(a\right)\Gamma \left(a\right)}{{\upsilon }^{a-1}\left(1-\upsilon \right)}^{b-1}$

in which the symbol Γ represents the gamma function.

{tilde over (Z)}={tilde over (Z)}₁, {tilde over (Z)}₂, . . . are identically distributed according to the a priori law G₀ which can be, for example, chosen as an inverse normal-Wishart model (NIW_{ρ, n0,μ0,Σ0}) such that {tilde over (Z)}_k=(m_k, Σ_k) with:

m _k|Σ_kN(m_k|μ₀, Σ_k/ρ)

Σ_k¹W(Σ_k¹|n₀,(n₀Σ₀)¹).

in which W designates the Wishart distribution for which the density is expressed by using the expression:

$\begin{array}{cc}f_{}\left(\sum _{}^{-1}\right)=\frac{n_0\sum _{0}\text{?}}{\text{?}\sum _{}\text{?}}\exp \left(-\frac{n_0}{2}\mathrm{Tr}\left(\underset{}{\sum _{}^{-1}}\sum _{0}\right)\right)\text{?}\text{indicates text missing or illegible when filed}& \left(4\right)\end{array}$

in which n₀ is the degree of freedom of the Wishart law and Σ₀ the associated scale matrix.

This is where the spatial concentration becomes involved as endogenous co-variable for the distribution of the kinetic parameters {tilde over (Q)}_k. First of all, the spatial distribution S(·) is introduced as marginal law of the space-time distribution H. The expression of S(·) is given by:

$S(·)={\int }_{ℳ\left(\right)}^{}H\left(·,\right)=\sum _{k=1}^{\infty }w_k\text{?}(·)$ $\text{?}\text{indicates text missing or illegible when filed}$

in which M(T) represents the set of the measurements over T.

It can be noted that, by the definitions of w_k and {tilde over (Z)}_k, S(·) is distributed according to a Dirichlet process of concentration parameter 0 and a basis measurement G₀:

S˜DP(θ, G₀)

{tilde over (Q)}={tilde over (Q)}₁, {tilde over (Q)}₂. . . are then independently distributed according to:

{tilde over (Q)} _k |{tilde over (Z)} _k , S K ₀ ^{λ k}

with λ_k=φ_λ({tilde over (Z)}_k, S) in which φ_λ(·) takes its values from Λ⊂. For example, the function φ_λ can be chosen as being the value of the expectation of the spatial density taken over the component K:

${\lambda }_k=\text{?}\left(\text{?}w_m\left(\text{?}\text{?}\right)\right)$ $\mathrm{with}$ $\text{?}~\left(\text{?}\text{?}\right)$ $\text{?}\text{indicates text missing or illegible when filed}$

It can be noted that, in the case where {tilde over (Q)}_k is independent of the marginal S such that λ_k=φ*_λ({tilde over (Z)}_k), we obtain:

{tilde over (Q)} _k |{tilde over (Z)} _k K ₀ ^{φ* k ({tilde over (Z)} k )}

{tilde over (Z)} _k , {tilde over (Q)} _k F ₀

expressions in which F₀({tilde over (Z)}_k, {tilde over (Q)}_k)=G₀({tilde over (Z)}_k)K₀^{φ*k({tilde over (Z)}k)}({tilde over (Q)}_k). This implies that H˜DP(θ, F₀).

The distribution of H is a Dirichlet process. In this situation the concentration is not involved as endogenous co-variable. However, the model then covers the possibility of having recourse to an exogenous co-variable originating from information of anatomical type for example.

The distribution of the K₀^λ remains to be established. K₀^Λ={K₀^λ: λ∈Λ} is defined as an innumerable collection of random probability measurements which accepts, as distribution, a local Dirichlet process denoted IDP, defined as follows. If:

• V=V₁, V₂, . . . are identically distributed such that V_jBeta (1, v),
• Q*=Q*₁, Q*₂, . . . are identically distributed such that Q*_jPT_{L(A, Q0}) is a Pólya tree with L-levels, of parameters A and of basis distribution Q₀,
• Γ=Γ₁, Γ₂, . . . are identically distributed such that Γ_jγ₀,then,

$K_0^{\lambda }(·)=\sum _{j=1}^{\infty }p_j\left(\lambda \right)\text{?}(·)$ $\text{?}\text{indicates text missing or illegible when filed}$

with p(λ)=p₁(λ),p₂(λ), . . . defined by the expression in which:

$p_j\left(\lambda \right)={\tilde{V}}_j\left(\lambda \right)\prod _{l<j}\left(1-{\tilde{V}}_l\left(\lambda \right)\right)$ $\mathrm{in}\mathrm{which}$ ${\tilde{V}}_j\left(\lambda \right)=V_j\left(d\left(\lambda ,{\Gamma }_j\right)<\psi \right)$

and d(λ, λ′) represents the Euclidian distance between λ and λ′.

It will be noted that if K₀^Λ−{K₀^λ: λ∈Λ}˜IDP(v, PT(A, Q₀), γ₀, ψ), then for any λ∈Λ, K₀^λ|λ˜DP(v, PT(A, Q₀)). The conditional law of K₀^λ knows that λ is a Pólya tree Dirichlet process.

It will be noted that in the limit case where ψ→∞, there is, for any λ∈Λ, K₀^λ=K₀ with K₀˜DP(v, PT(A, Q₀)), consequently

H˜DP(θ, G₀×K₀)

The definition of a Pólyatree to complete the model remains to be reviewed. We adopt the Lavine formulism presented in his paper entitled Some aspects of Pólya tree distributions for statistical modelling, Ann. Statist., 20:1222-1235, 1992, for the definition of the dyadic Pólya tree. Take E={0, 1} and the E^l the l-cartesian product E× . . . ×E with E⁰=. Let E*=∪_l=0^∞E^l. Let B_={T} and let Π={B₀, B₁, B₀₀, B₀₁, . . . } and separable binary tree of partitions of T such that sets B_{∈1 . . . ∈l+1} are obtained by binary subdivision of the sets B_{∈1 . . . ∈l} and ∪_l=0^∞B_{∈1 . . . ∈l+1} defines the measurable sets. The index 0 (respectively 1) is interpreted as 0 (respectively 1).

It is said that the random measurement Q over T accepts as distribution a Pólya tree of parameters (A, Q₀) that will be noted Q˜PT(A, Q₀), if there are positive real values A={α_∈: ∈∈E*} and random variables W={W_∈: ∈∈E*} such that:

• • the collection W consists of mutually independent random variables,
  • ∀_∈∈E*, W_∈˜Beta (α_∈0, α_∈1),
  • ∀l=1, 2, . . . et ∀_∈∈E^l,

$Q\left(B_{{\epsilon }_1\dots {\epsilon }_l}\right)=\left(\prod _{\underset{\left\{k:{\epsilon }_k=0\right\}}{k=1}}^lW_{{\epsilon }_1\dots {\epsilon }_{k-1}}\right)\left(\prod _{\underset{\left\{k:{\epsilon }_k=1\right\}}{k=1}}^lW_{{\epsilon }_1\dots {\epsilon }_{k-1}}\right)$

with factors equal to W_ in which 1−W_ if k−1.Consider the probability measurement Q₀. The Pólya tree can be centered around Q₀ by taking:

$B_{{\epsilon }_1\dots {\epsilon }_l}=\left(Q_0^{-1}\left(\frac{\sum _{k=1}^l{\epsilon }_k2^{k-1}}{2^l}\right),Q_0^{-1}\left(\frac{\sum _{k=1}^l{\epsilon }_k2^{k-1}+1}{2^l}\right)\right)$ $\mathrm{and}{\alpha }_{\epsilon 0}={\alpha }_{\epsilon 1},\forall \epsilon \in E^{*}$

In this case, we note α_∈=a_l for l=1, 2, . . . and ∀_∈∈E^l. This construction scheme defines the canonical Pólya tree. In this context, if Q˜PT(A, Q₀), then (Q)=Q₀.

The construction of the partition can be continued to infinity but, since it is not possible to consider computing an infinite tree in practice, the Pólya trees are often simplified by specifying the parameters and by constructing the partition up to a bounded level L. If we note f_Q(·), the density of probability of Q˜PT_L(A, Q₀), then:

∀t∈B_{∈1 . . . ∈L}, f_Q(t)=2^Lf_Qn(t)·Q(B_{∈1 . . . ∈L})

The justification for a partial specification of the Pólya tree can be supported by the fact that, for the continuous distributions, the parameters rapidly increase as a function of the levels of the tree. The updating of the law from the data then has little a priori influence on the law because the number of observations in the high levels become small compared to the a priori parameter.

Since the objective is to culminate with an algorithm for reconstruction of the functional volumes, the method according to the invention can be provided with auxiliary variables which will render the iterative reconstruction procedure and the implementation of an MCMC-type method, MCMC standing for “Markov chain Monte-Carlo”, effective.

We introduce C=C₁, C₂, . . . , C_n, the data classification variables, to the components of H such that (Z_i, Q_i)=({tilde over (Z)}_Ci, {tilde over (Q)}_Ci) for any i<n. Similarly, let D=D₁, D₂, . . . such that {tilde over (Q)}_k=Q*_Dk for any K.

The auxiliary variables for the sampling by slices, usually referred to by the expression “slice sampling”, are described hereinbelow.

In order to implement an effective Gibbs sampler, we follow a so-called slice sampling technique as described in the paper by M. Kalli, J. E. Griffin and S. G. Walker entitled Slice sampling mixture models, Statistics and Computing, 21(1):93-105, 2011. This approach requires the introduction of additional variables.

Let u=u₁, u₂, . . . , u_n be uniform random variables. The joint density of (X_i, T_i, u_i), for any positive sequence (ξ_k), can be determined by using the expression:

$\left(X_i,T_i,u_i|w,\tilde{Z},\tilde{Q}\right)\sim \sum _{k=1}^{\infty }\left(u_i|0,{\xi }_k\right)w_k\left(X_i|{\tilde{Z}}_k\right){\tilde{Q}}_k\left(T_i\right)$

• in which u(·|a, b) is the uniform distribution over ]a, b]. We easily check that, by marginalization on u_i of the preceding formula, we obtain the distribution of (X_i, T_i|w, {tilde over (Z)}, {tilde over (Q)}). We propose adopting a dependent slice variable (ξ_k), that is to say that, for any K, ξ_k=min (w_k, ζ)
• in which ζ∈]0, 1] and is independent of w_k. In practice, values of ζ corresponding to the expectation of the weight of the first component that does not have any datum allocated offer satisfactory mixture properties

$\left(\zeta \approx \frac{0}{\left(n+0\right)\left(1+0\right)}\right).$

•  Provided with this choice of ξ_k, the following expression is obtained:

$\left(X_i,T_i,u_i|w,\tilde{Z},\tilde{Q}\right)\sim {\zeta }^{-1}\sum _{u_i<\zeta <w_k}w_k\left(X_i|{\tilde{Z}}_k\right){\tilde{Q}}_k\left(T_i\right)+\sum _{u_i<w_k\le \zeta }\left(X_i|{\tilde{Z}}_k\right){\tilde{Q}}_k\left(T_i\right)$

• in which the two sums are finite by virtue of #{j: w_j>ε}<∞ for any ε>0.

Similarly, let v=v₁, v₂, . . . be uniform auxiliary variables and ζ∈]0, 1]. Then, for any K, the following applies:

$\left({\tilde{Q}}_k,{\upsilon }_k|V,Q^{*},\Gamma ,w,{\tilde{Z}}_k,S\right)\sim {\varsigma }^{-1}\sum _{{\upsilon }_k<\varsigma <p_j\left({\lambda }_k\right)}p_j\left({\lambda }_k\right){\delta }_{Q_j^{*}}\left({\tilde{Q}}_k\right)+\sum _{{\upsilon }_k<p_j\left({\lambda }_k\right)\le \varsigma }{\delta }_{Q_j^{*}}\left({\tilde{Q}}_k\right)$

expression in which λ_k=φ_λ({tilde over (Z)}_k, S) and p_j(λ_k) are defined above.

The model (2) is then rewritten by adding to it the hidden and auxiliary variables:

$\begin{array}{cc}Y_i|X_i\stackrel{\mathrm{ind}}{\sim }\left(Y_i|X_i\right)X_i,T_i|C_i,D\stackrel{\mathrm{ind}}{\sim }\left(X_i|{\tilde{Z}}_{C_i}\right)\times Q_{D_{C_i}}^{*}\left(T_i\right)C_i,u_i|w\stackrel{\mathrm{iid}}{\sim }\sum _{k=1}^{\infty }\left(u_i|0,\min \left(\zeta ,w_k\right)\right)w_k{\delta }_k\left(C_i\right)D_k,{\upsilon }_k|V,\Gamma ,w,\tilde{Z}\stackrel{\mathrm{ind}}{\sim }\sum _{j=1}^{\infty }\left({\upsilon }_k|0,\min \left(\varsigma ,p_j\left({\lambda }_k\right)\right)\right)p_j\left({\lambda }_k\right){\delta }_j\left(D_k\right)w\sim \mathrm{GEM}\left(\theta \right){\tilde{Z}}_k\stackrel{\mathrm{iid}}{\sim }{\mathrm{ℐ}}_{\varrho ,n_0,{\mu }_0,{\Sigma }_0}V_j\stackrel{\mathrm{iid}}{\sim }\mathrm{Beta}\left(1,\nu \right){\Gamma }_j\stackrel{\mathrm{iid}}{\sim }{\Upsilon }_0Q_j^{*}\stackrel{\mathrm{iid}}{\sim }{\mathrm{PT}}_L\left(,Q_0\right)\mathrm{with}p_j\left({\lambda }_k\right)={\tilde{V}}_j\left({\lambda }_k\right){\Pi }_{l<j}\left(1-{\tilde{V}}_l\left({\lambda }_k\right)\right),{\tilde{V}}_j\left({\lambda }_k\right)-V_j\left(d\left({\lambda }_k,{\Gamma }_j\right)<\psi \right)\mathrm{and}{\lambda }_k={}_{{\tilde{x}}_k}\left(\sum _{l=1}^{\infty }w_l\left({\tilde{x}}_k|{\tilde{Z}}_l\right)\right)& \left(5\right)\end{array}$

From this generative model, it is possible to express the different a posteriori conditional laws which will be involved in the Gibbs sampler.

The sampling of the a posteriori conditional laws by Markov chain Monte-Carlo MCMC algorithm is described hereinbelow with the support of FIGS. 2 to 7.

From the expression (5), the MCMC sampler generates, in succession, samples originating from the following conditional distributions:

• Places of emission/allocation to the components of H: (X, C|w,u,{hacek over (Z)}, Q*, D, Y, T)
• Spatial parameters: ({hacek over (Z)}|C, X)
• Weights and auxiliary variables of H: (w, u|C)
• Allocation of the kinetics to the components of K₀^Λ: (D|V, v, Q*, Γ, w, {hacek over (Z)}, C, T)
• Pólya trees of the kinetics: (Q*|D, C, T)

Weightings, local variables and auxiliaries of K₀^Λ: (V, Γ, v|D, w, {hacek over (Z)}, T)

FIG. 2 illustrates how the weights and auxiliary variables of the distribution H of the space-time components are generated, that is to say how samples of conditional distribution (w,u|C) are generated.

We denote k_n the latent number of components of the space-time mixture for n observations and n_k=#{i: C_i=k} for k≦k_n. We adopt the approach presented in the paper by Kalli et al., previously cited, and proceed with the joint generation of (w, u|C) by first of all sampling w₁, w₂, . . . , w_kn|C, then u|w₁, w₂, . . . , w_kn, C, and finally w_kn+1, w_kn+2, . . . |u, w₁, w₂, . . . , w_kn, C.

In this context, by following the new definition and the corollary twenty of the paper by J. Pitman entitled Some developments of the Blackwell-MacQueen urn scheme in T. F. et al., editor, Statistics, Probability and Game Theory; Papers in honor of David Blackwell, volume 30 of Lecture Notes—Monograph Series, pages 245-267, Institute of Mathematical Statistics, 1996, we are able to express the a posteriori law of the Dirichlet process as the sum of a finite sum of components whose weights are distributed by a Dirichlet distribution and of a data-independent Dirichlet process. The presence of the auxiliary variable u makes it possible to generate only a finite number k* of weights and of places of emission in the generation of the independent Dirichlet process.

In a first stage 200, the w_k are generated for k≦k_n.

(w₁, w₂, . . . , w_kn, r_kn|C)˜Dirichlet (n₁, n₂, . . . , n_kn, θ)

In a second stage 201, the samples u_i|w₁, w₂, . . . , w_kn, C are generated.

u _i |C u(u_i|0, min (w_Ci, ζ))

u*=min(u)

In a third stage 202, the samples w_k are generated for k>k_n.As long as r_k−1>u*,

U _k˜Beta (1, θ)

w _k =U _k r _k−1

r _k =r _k−1(1−U_k)

In a fourth stage 203, k*=min({k: r_k<u*}) is posited.Clearly, w_k<u* for any k>k*. For this reason, a finite set of w_k is sufficient for the allocation of the data to the components. It should be added that we also abandon in the computations all the w_k such that w_k<u* for k≦k_n.

FIG. 3 schematically illustrates how the spatial parameters are generated, that is to say how samples of conditional distribution {tilde over (Z)} C, X) are generated.

In a first stage 300, ({tilde over (Z)}_k|C, X) are generated for any k≦k_n, from an inverse normal-Wishart law whose parameters are updated from the observations, through the conjugation property:

$m_k|{\Sigma }_k,C,X\stackrel{\mathrm{ind}}{\sim }\left({\mu }_k,\frac{{\Sigma }_k}{{\varrho }_k}\right)$ $\mathrm{with}$ ${\varrho }_k=\varrho +n_k$ ${\mu }_k=\frac{\varrho {\mu }_0+n_k}{\varrho +n_k}$ $=\frac{1}{n_k}\sum _{\left\{i:C_i=k\right\}}X_i$ $\mathrm{and}$ ${\Sigma }_k^{-1}|C,X\stackrel{\mathrm{ind}}{\sim }\left({\eta }_k,{\left({\eta }_k{\tilde{\Sigma }}_k\right)}^{-1}\right)$ $\mathrm{with}$ ${\eta }_k=n_0-n_k$ ${\tilde{\Sigma }}_k=\frac{1}{{\eta }_k}\left[n_0{\Sigma }_0+n_k+\frac{\left({\mu }_0-\right){\left({\mu }_0-\right)}^{\top }}{\frac{1}{\varrho }+\frac{1}{n_k}}\right]$ $=\frac{1}{n_k}\sum _{\left\{i:C_i=k\right\}}\left(X_i-\right){\left(X_i-\right)}^{\top }$

In a second stage 301, the samples {tilde over (Z)}_k are generated:

{tilde over (Z)}_kNIW_{ρ, n0,μ0,Σ0}, for k_n<k≦k*

FIG. 4 schematically illustrates one way of determining the local and auxiliary variables of the set K₀^λ.

For this, (V, Γ, v|D, w, {tilde over (Z)}, T) must be generated.

We consider as a priori law γ₀(·)=u(·|a_Γ, b_Γ). For any j, we define

D _j ={k≦k _n : D _k =j}

D _j ⁺ ={k≦k _n : D _k >j}

and designate, for any j, m_j=#D_j as the number of components of H allocated to the j^th component of K₀^Λ.We also define the greatest index of the components of K₀^Λ having at least one component of H allocated:

J _{k n} =max({j: m_j>0})

We adopt the same generation scheme as for the sampling of (w, u|C), that is to say we first generate (V_j, Γ_j|D, w, {tilde over (Z)}, T) for j≦J_kn 400, then

$\left(vD,V_1,\cdots ,V_{J_{{\kappa }_n}},{\Gamma }_1,\cdots ,{\Gamma }_{J_{{\kappa }_n}},w,\tilde{Z}\right)$

401 and, finally, V_j and Γ_j for j>j_kn 402.For any j≦J_kn, the following are generated independently:

$V_j|D,T\stackrel{\mathrm{ind}}{\sim }\mathrm{Beta}\left(1+m_j,\nu \sum _{l=j+1}^{J_{{\kappa }_n}}m_l\right)$

For any j≦J_kn, we posit:

$a_{j,\Gamma }^{\prime }=\max \left(\underset{k\in {}_j}{\max }\left({\lambda }_k-\psi \right),a_{\Gamma }\right)$ $b_{j,\Gamma }^{\prime }=\min \left(\underset{k\in {}_j}{\min }\left({\lambda }_k+\psi \right),b_{\Gamma }\right)$

then generate:

${\Gamma }_j|V_j,D,w,\tilde{Z}\stackrel{\mathrm{ind}}{\sim }{}_j^{*}\left({\Gamma }_j\right)$ $\mathrm{in}\mathrm{which}$ ${}_j^{*}\left({\Gamma }_j\right)\propto \left(a_{j,\Gamma }^{\prime }<{\Gamma }_j<b_{j,\Gamma }^{\prime }\right)\times \prod _{k\in {}_j^{+}}\left(1-V_j\left(d\left({\lambda }_k,{\Gamma }_j\right)<\psi \right)\right)$

The sampling of u*_j requires a Metropolis-Hastings (MH) step. We propose an independent MH sampler with, as candidate distribution:

g(Γ_j)=u(Γ_j|a′_j,Γ, b′_j,Γ)

Thus, on the iteration \|⊥, for any j≦J_kn,

• • generate Γ*_ju(Γ_j|a′_j,Γ, b′_j,Γ).
  • Compute the acceptance ratio w(Γ*_j, Γ_j^(t))

$\omega \left({\Gamma }_j^{*}·{\Gamma }_j^{\left(t\right)}\right)=\begin{array}{c}\prod _{k\in {}_j^{+}}\left(1-V_j\left(d\left({\lambda }_k,{\Gamma }_j^{*}\right)<\psi \right)\right)\\ \prod _{k\in {}_j^{+}}\left(1-V_j\left(d\left({\lambda }_k,{\Gamma }_j^{\left(t\right)}\right)<\psi \right)\right)\end{array}$

in which Γ_j^(t) is the value of the Markov chain on the iteration t.

• • Test the acceptance of the proposition Γ*_j, i.e. assign Γ_j^(t+1)=Γ*_j with a probability

min(1, ω(Γ*_j, Γ_j^(t))),

otherwise assign Γ_j^(t+1)=Γ_j^(t).

• • For any k≦k_n, posit ρ_j,k=ρ_j−1,k(1−V_j(d(λ_k, Γ_j)<ψ)) and ρ*_j=min(ρ_j,1, . . . , ρ_j,kn).
  • For any k≦k_n, independently generate

${\upsilon }_k|D,V_1,\dots ,V_{J_{{\kappa }_n}},{\Gamma }_1,\dots ,{\Gamma }_{J_{{\kappa }_n}},w,\tilde{Z}\stackrel{\mathrm{ind}}{\sim }\left({\upsilon }_k|0,\min \left(p_{D_k}\left({\lambda }_k\right),\varsigma \right)\right)$

and posit

v*=min(v)

• • Generate V_j et Γ_j for j>J_kn. As long as ρ*_j−1>v*,

V _j Beta(1, v)

Γ_ju(Γ_j|a_Γ, b_Γ)

• • For any k≦k_n, posit ρ_j,k=ρ_j−1,k(1−V_j(d(λ_k, Γ_j)<ψ)) and ρ*_j=min(ρ_j,1, . . . , ρ_j,kn).
  • Posit J*=min({j: ρ*_j<v*}).

FIG. 5 illustrates how the kinetics are allocated to the components of K₀^λ, that is to say how (D|V, v, Q*, Γ, w, {tilde over (Z)}, T) is generated.

In a first stage 500, the following is generated independently for any k≦k_n:

$\left(D_k|V,v,Q^{*},\Gamma ,w,\tilde{Z},C,T\right)\stackrel{\mathrm{ind}}{\sim }\sum _{j=1}^{J^{*}}p_{j,k}{\delta }_j(·)$ $\mathrm{with}$ $p_{j,k}\propto \left(p_j\left({\lambda }_k\right)>{\upsilon }_k\right)\max \left(p_j\left({\lambda }_{\kappa }\right),\varsigma \right)\prod _{\left\{i:C_i=k\right\}}f_{Q_j^{*}}\left(T_i\right)$

and Σ_j=1^J*P_j,k−1.In a second stage 501, the following is generated for any K, such that k_n<k<k*:

$D_k|V,v,\Gamma ,w,\tilde{Z}\stackrel{\mathrm{ind}}{\sim }\sum _{j=1}^{J^{*}}p_{j,k}{\delta }_j(·)$ $\mathrm{with}$ $p_{j,k}\propto \left(p_j\left({\lambda }_k\right)>{\upsilon }_k\right)\max \left(p_j\left({\lambda }_{\kappa }\right),\varsigma \right)$ $\mathrm{and}\sum _{j=1}^{J^{*}}p_{j,k}-1.$

The Pólya trees of the generated kinetics must then be generated. According to the conjugation property of the Pólya tree, for any j≦J_kn, (Q*_j|D, C, T) is generated independently according to a Pólya tree distribution with updated parameters A′_j−{α′_j,τ: ε∈E*} such that

Q* _j |D, C, T PT _L(A′_j, Q₀)

with

α′_j,τ−α_τ+n_j,τ

and

n _j,∈ =#{i∈{1, . . . , n}: (D_Ci=j)̂(T_i∈B_∈)}

Then, for any J_kn<j≦J*, independently generate

Q* _j |D, C, T PT _L(A, Q₀)

FIG. 6 shows how the places of emission are determined and how the components of H are determined. This is by generating (C, X|w, u, {tilde over (Z)}, Q*, D, Y, T).

The conditional law has a form which does not allow for a direct generation. In effect, it involves a projector P(Y_i|X_i) which incorporates geometrical and physical characteristics of the detection device, which in practice precludes finding an analytical expression for the distribution that offers a hope of simple generation. We therefore have recourse to a Metropolis-Hastings step. We begin by expressing the conditional distribution by using the Bayes rule:

$C_i,X_i|Y_i,T_i,u_i,w,\tilde{Z},Q^{*},D\stackrel{\mathrm{ind}}{\sim }\left(Y_i|X_i\right)\times {\tilde{G}}_i\left(C_i,X_i\right)$ $\mathrm{with}$ $f_{{\tilde{G}}_i}\left(C_i,X_i\right)\propto \sum _{k=1}^{{\kappa }^{*}}\left(w_k>u_i\right)\max \left(w_k,\zeta \right)f_{}\left(X_i|{\tilde{Z}}_k\right)f_{Q_{D_k}^{*}}\left(T_i\right){\delta }_k\left(C_i\right)$

and P(Y_i|X_i) corresponds to the likelihood of the projector for the sets of observed events (truncated). Since this likelihood does not accept any standard form allowing for direct sampling, we use an independent Metropolis-Hastings algorithm and have to define a candidate law that we denote G*_i for i=1, . . . , n. On each iteration t+1 we have to generate (C*_i, X*_i)˜G*_i and compute the acceptance ratio. This step is the most costly in terms of computations.

We thus propose a candidate law of the form:

f _{G* i} (C*_i, X*_i)∝f_N(X*_i|m*_Yi, Σ*_Yi)×f_Gi(C*_i, X*_i)

We replace the projector P(Y_i|X_i) with N(X_i|m*_Yi, Σ*_Yi) in the conditional law to construct the candidate distribution 600, in which m*_Yi and Σ*_Yi are chosen in such a way as to suitably approximate P(Y_i|X_i). We can then express f*_Gi(·) by using the expression:

$f_{G_i^{*}}\left(C_i^{*},X_i^{*}\right)\propto \sum _{k=1}^{{\kappa }^{*}}\left(w_k>u_i\right)\max \left(w_k,\zeta \right){\stackrel{\_}{\eta }}_{k,Y_i}f_{}\left(X_i^{*}|{\stackrel{\_}{m}}_{k,Y_i},{\stackrel{\_}{\Sigma }}_{k,Y_i}\right)f_{Q_{D_k}^{*}}\left(T_i\right){\delta }_k\left(C_i^{*}\right)$

with, for any k≦k* and m≦M,

η _k,m =f _N(m*_m|m_k, Σ*_m−Σ_k)

m _k,m=((Σ*_m)⁻¹+Σ_k⁻¹)⁻¹((Σ*_m)⁻¹m*_m+Σ_k⁻¹m_k)

Σ _k,m=((Σ*_m)⁻¹+Σ_k⁻¹)⁻¹

It is then possible to proceed with the sampling of C_i, X_i), for any i≦n on the iteration t+1.C*_i is generated independently 601 such that:

$C_i^{*}\stackrel{\mathrm{ind}}{\sim }\sum _{k=1}^{{\kappa }^{*}}{\tilde{w}}_{k,i}{\delta }_k(·)$ $\mathrm{with}$ ${\tilde{w}}_{k,i}\propto \left(w_k>u_i\right)\max \left(w_k,\zeta \right){\tilde{\eta }}_{k,Y,}f_{Q_{D_k}^{*}}\left(T_i\right)$ $\mathrm{and}$ $\sum _{k=1}^{{\kappa }^{*}}{\stackrel{\_}{w}}_{k,i}=1.$

Then X*_i is generated independently 602 such that:

X* _i N(X*_i| m_C*i,Yi, Σ_C*i,Yi)

The acceptance ratios ω((C*_i, X*_i), (C_i^(t), X_i^(t))) are then computed 603:

$\begin{array}{c}\omega \left(\left(C_i^{*},X_i^{*}\right),\left(C_i^{\left(t\right)},X_i^{\left(t\right)}\right)\right)=\frac{\left(Y_i|X_i^{*}\right)}{\left(Y_i|X_i^{\left(t\right)}\right)}·\frac{f_{}\left(X_i^{\left(t\right)}|m_{Y_i}^{*},{\Sigma }_{Y_i}^{*}\right)}{f_{}\left(X_i^{*}|m_{Y_i}^{*},{\Sigma }_{Y_i}^{*}\right)}\\ =\omega \left(X_i^{*},X_i^{\left(t\right)}\right)\end{array}$

It will be able to be noted that {tilde over (G)}_i(C_i, X_i) and C_i disappear from the acceptance ratio.

It then remains to test 604 the acceptance of (C*_i, X*_i), that is to say assign (C_i^(t+1), X_i^(t+1))=(C*_i, X*_i) with the probability

min(1,ω(X*_i, X_i^(t)))

otherwise assign (C_i^(t+1), X_i^(t+1))=(C_i^(t), X_i^(t)))FIG. 7 presents a way of computing the concentration λ_k. The evaluation of λ_k is necessary for the generation of (V, Γ, v|D, w, {tilde over (Z)}, T) and of (D|V, v, Q*, Γ, w, {tilde over (Z)}, C, T).The expression of λ_k for the particular choice of φ_λ({tilde over (Z)}_k, S) given previously, can be determined according to the equation (5) by using the expression:

${\lambda }_k={}_{{\tilde{x}}_k}\left(\sum _{m=1}^{\infty }w_m\left({\tilde{x}}_k|{\tilde{Z}}_m\right)\right)$

This expression can be evaluated by a significant sampling technique. For this, for a given number of particles T_λ and for any k, the following is generated independently 700 for any l≦T_λ:

{tilde over (x)} _kl N({tilde over (x)}_kl|{tilde over (Z)}_k)

Then, the following is computed 701 for any l≦T_λ:

${\tilde{\lambda }}_{\mathrm{kl}}\approx \sum _{m=1}^{{\kappa }^{*}}w_m\left({\tilde{x}}_{\mathrm{kl}}|{\tilde{Z}}_m\right)$

In a third stage 702, λ_k is evaluated:

${\lambda }_k\approx \frac{1}{T_{\lambda }}\sum _{l=1}^{T_{\lambda }}{\tilde{\lambda }}_{\mathrm{kl}}$

The functional volumes VF then remain to be determined 105.

Once the MCMC algorithm has reached balance equilibrium, the sampler provides draws of the a posteriori joint law of (X, C, D, w, u, {tilde over (Z)}, Q*|Y, T). For each draw of parameters retained, it is possible to construct a sample of the space-time activity distribution of the events recorded on the iteration (t):

$\begin{array}{cc}{f}_G^{\left(t\right)}\left(x,t\right)=\sum _{k=1}^{{\kappa }^{*}}w_kf_{}\left(x|{\tilde{Z}}_k\right)f_{{\tilde{Q}}_k}\left(t\right),& \left(6\right)\end{array}$

The space-time density corresponding to all the events, which are not necessarily observable, is obtained by normalization.Similarly, the entire distribution of the functional volumes is accessible. We marginalize over time and define the spatial distribution of a functional volume j associated with the kinetic Q*_j on the iteration (1) by using the expression:

$\begin{array}{cc}{f}_{{\mathrm{VF}}_j}^{\left(t\right)}\left(x\right)=\sum _{k:{\tilde{Q}}_k=Q_j^{*}}^{{\kappa }^{*}}w_kf_{}\left(x|{\tilde{Z}}_k\right).& \left(7\right)\end{array}$

All the elements are then available to determine the estimators of space-time activity distribution and functional volumes VF for T iterations of the algorithm that are retained. The following expressions can then be used for that:

$\begin{array}{cc}\left(f_G\left(x,t\right)|Y,T\right)\approx \frac{1}{T}\sum _{t=1}^Tf_G^{\left(t\right)}\left(x,t\right)\mathrm{and}& \left(8\right)\\ \left(f_{{\mathrm{VF}}_j}\left(x\right)|Y,T\right)\approx \frac{1}{T}\sum _{t=1}^Tf_{{\mathrm{VF}}_j}^{\left(t\right)}\left(x\right)& \left(9\right)\end{array}$

Similarly, any other functional of the a posteriori distributions can be computed, like a conditional variance or a credibility interval.

Claims

1. A method for determining functional volumes for the search for kinetics representing the trend of the concentration of a radioactive tracer in an area of biological tissue, the method being applied to spatial components and comprising the following steps applied iteratively according to a Markov chain Monte-Carlo scheme: a step of generation of a set K0λ made up of a set of candidate kinetics associated with values of probability of appearance of these kinetics, these values depending on the concentration λ of radioactive tracer;a labeling step during which, for each spatial component of index K, the probabilities of selection of the kinetics of the set K0λ are weighted by introducing therein a function λk representative of the concentration of radioactive tracer in this component so as to obtain a set of indicative values, an indicative value Dk designating the kinetic with which the spatial component K is associated;a step of construction of functional volumes, a functional volume VFj being made up of the set of the spatial components which share the same indicative value Dk.

2. The method as claimed in claim 1, further comprising a step of computation of the averages of the functional volumes VFj obtained during the iterative process.

3. The method as claimed in claim 1, comprising a step of determination of statistical estimators concerning the functional volumes VFj obtained during the iterative process so as to quantify the estimation uncertainty.

4. The method as claimed in claim 3, in which the variance of the functional volumes VFj obtained during the iterative process is determined.

5. The method as claimed in claim 3, in which a credible interval is determined.

6. The method as claimed in claim 1, in which the set K0λ is obtained by the determination (100) of two sets of scalar variables Vj and Γj, these variables relating to kinetic components of index j, in which Vj represents a fixed weighting coefficient independent of the concentration λ and Γj represents a reference concentration value associated with the kinetic of index j.

7. The method as claimed in claim 6, in which Vj is determined by using the expression:

8. The method as claimed in claim 6, in which Γj is determined using the expression:

9. The method as claimed in claim 1, in which the labeling step (109) is performed: by independently generating (500) for any k≦kn:

10. The method as claimed in claim 1, in which the concentration λk is obtained by the application of a function φλ defined such that λk=φλ({tilde over (Z)}k, S) with:

11. A positron emission tomography device implementing the method as claimed in claim 1.