François Rousseau
email: rousseau@unistra.fr

Piotr A Habas Bicg

Colin Studholme
email: colin.studholme@ieee.org

A supervised patch-based approach for human brain labeling

We propose in this work a patch-based image labeling method relying on a label propagation framework. Based on image intensity similarities between the input image and an anatomy textbook, an original strategy which does not require any non-rigid registration is presented. Following recent developments in non-local image denoising, the similarity between images is represented by a weighted graph computed from an intensity-based distance between patches. Experiments on simulated and in-vivo MR images show that the proposed method is very successful in providing automated human brain labeling.

Introduction

Automated brain labeling [START_REF] Miller | Mathematical textbook of deformable neuroanatomies[END_REF][START_REF] Collins | Automatic 3-D model-based neuroanatomical segmentation[END_REF][START_REF] Lancaster | Automated labeling of the human brain: A preliminary report on the development and evaluation of a forward-transform method[END_REF][START_REF] Klein | Mindboggle: a scatterbrained approach to automate brain labeling[END_REF] is a crucial step in neuroimaging studies since it provides a segmentation map of contiguous meaningful brain regions. Label propagation (also called label fusion) is a way to segment an image and it is usually performed by registering one or several labeled images to the image to be segmented. The labels are then mapped and fused using classification rules to produce a new label map of the input image [START_REF] Rohlfing | Performance-based classifier combination in atlasbased image segmentation using expectation-maximization parameter estimation[END_REF][START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF][START_REF] Artaechevarria | Combination strategies in multi-atlas image segmentation: application to brain MR data[END_REF][START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF][START_REF] Isgum | Multi-atlas-based segmentation with local decision fusion-application to cardiac and aortic segmentation in CT scans[END_REF][START_REF] Collins | Towards accurate, automatic segmentation of the hippocampus and amygdala from MRI by augmenting ANIMAL with a template library and label fusion[END_REF][START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF][START_REF] Van Rikxoort | Adaptive local multi-atlas segmentation: application to the heart and the caudate nucleus[END_REF][START_REF] Lötjönen | Fast and robust multi-atlas segmentation of brain magnetic resonance images[END_REF][START_REF] Yushkevich | Nearly automatic segmentation of hippocampal subfields in in vivo focal t2-weighted MRI[END_REF][START_REF] Sdika | Combining atlas based segmentation and intensity classification with nearest neighbor transform and accuracy weighted vote[END_REF][START_REF] Langerak | Label fusion in Atlas-Based segmentation using a selective and iterative method for performance level estimation (SIMPLE)[END_REF][START_REF] Khan | Optimal weights for local multi-atlas fusion using supervised learning and dynamic information (SuperDyn): Validation on hippocampus segmentation[END_REF]. The principle of such registration-based labeling approach is shown in Figure 1. This segmentation approach is highly versatile since the main prerequisite is an anatomy textbook, i.e. a set of measured images (such as Magnetic Resonance (MR) imaging or Computerized Tomography (CT)) and the corresponding label maps. The use of label propagation has been extensively investigated for automatic brain MR segmentation, especially for structures such as hippocampus, caudate, putamen, amygdala, etc., and cortex areas. The keypoints of registration-based label propagation approaches concern the accuracy of the non-rigid registration, the fusion rules [START_REF] Duda | Pattern Classification[END_REF][START_REF] Rohlfing | Performance-based classifier combination in atlasbased image segmentation using expectation-maximization parameter estimation[END_REF][START_REF] Warfield | Simultaneous truth and performance level estimation (STA-PLE): an algorithm for the validation of image segmentation[END_REF][START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF][START_REF] Artaechevarria | Combination strategies in multi-atlas image segmentation: application to brain MR data[END_REF], the selection of the labeled images [START_REF] Wu | Optimum template selection for atlas-based segmentation[END_REF][START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF][START_REF] Van Rikxoort | Adaptive local multi-atlas segmentation: application to the heart and the caudate nucleus[END_REF] and the labeling errors in primary manual segmentation.

Because of the use of a registration algorithm, one makes the implicit (and strong) assumption that there exists a one-to-one mapping between the input image and all the anatomical images of the textbook. In the case where no one-to-one mapping is possible, a registration-based labeling framework propagates incorrect labels. Moreover, local incorrect matching due to inherent registration errors can also lead to segmentation errors. Finally, as shown in a recent evaluation study [START_REF] Klein | Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration[END_REF], there is a non-negligible discrepancy in term of quality of matching between non-rigid registration techniques. Thus, even if several algorithms are now freely available, the non-rigid registration procedure remains a complex step to setup and it is usually quite CPU time consuming.

In this work, following recent developments in non-local image denoising [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Katkovnik | From local kernel to nonlocal Multiple-Model image denoising[END_REF], we propose an alternate strategy for label propagation which does not require any non-rigid registration. The proposed algorithm makes use of local similarities between the image to be labeled and the images contained in the anatomy textbook. The key idea is similar to a fuzzy block matching approach which avoids the constraint of a strict one-to-one mapping. The method described in this article has been developed irrespective of a similar patch-based approach recently proposed by Coupé et al. [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF]. The contributions of our work can be summarized as follows: 1) a patch-based framework for automated image labeling, 2) investigation of several patch aggregation strategies (pair-wise vs group-wise labeling and pointwise vs multipoint estimation), 3) comparison of fusion rules (majority voting rule vs STAPLE [START_REF] Warfield | Simultaneous truth and performance level estimation (STA-PLE): an algorithm for the validation of image segmentation[END_REF]), 4) application to human brain labeling using publicly available simulated (Brainweb) and in-vivo (IBSR and NA0-NIREP) datasets, 5) comparison with a non-rigid registration-based technique.

2 Patch-based Label Propagation

Patch-based Principle

Recently, Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] have proposed a very efficient denoising algorithm relying on a non-local framework. Since then, this non-local strategy has been studied and applied in several image processing applications such as non-local regularization functionals in the context of inverse problems [START_REF] Kindermann | Deblurring and denoising of images by nonlocal functionals[END_REF][START_REF] Mignotte | A non-local regularization strategy for image deconvolution[END_REF][START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Peyré | Non-local regularization of inverse problems[END_REF][START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF] or medical image synthesis [START_REF] Roy | MR contrast synthesis for lesion segmentation[END_REF].

Let us consider, over the image domain Ω, a weighted graph w that links together the voxels of the input image I with a weight w(x, y), (x, y) ∈ Ω 2 . This weighted graph w is a representation of non-local similarities in the input image I.

In [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], the non-local graph w is used for denoising purpose using a neighborhood averaging strategy (called non-local means (NLM)):

∀x ∈ Ω, I nlm (x) = y∈Ω w(x, y)I(y) y∈Ω w(x, y) (1)
where w is the graph of self-similarity computed on the noisy image I, I(y) is the gray level value of the image I at the voxel y and I nlm is a denoised version of I.

The weighted graph reflects the similarities between voxels of the same image. It can be computed using a intensity-based distance between patches [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF]:

w(x, y) = f x ′ ∈PI (x),y ′ ∈PI (y) (I(x ′ ) -I(y ′ )) 2 2N β σ2 (2)
where P I (x) is a 3D patch of the image I centered at voxel x; f is a kernel function (f (x) = e -x in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]), N is the number of voxels of a 3D patch; σ is the standard deviation of the noise and β is a smoothing parameter.

With the assumption of Gaussian noise in images, β can be set to 1 (see [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] for theoretical justifications) and the standard deviation of noise σ is estimated via pseudo-residuals as defined in [START_REF] Gasser | Residual variance and residual pattern in nonlinear regression[END_REF].

Label Propagation

In this work, we propose to investigate the use of such a non-local patch-based approach for label propagation.

Let T be an anatomy textbook containing a set of T1-weighted MR images I and the corresponding label maps

L: T = {(I i , L i ), i = 1, • • • , n}.

Weighted Graph and Label Propagation

Let us consider, over the image domain Ω, a weighted graph w i that links together voxels x of the input image I and voxels y of the image I i with a weight w i (x, y), (x, y) ∈ Ω 2 . w i is computed as follows:

w i (x, y) = f x ′ ∈PI (x),y ′ ∈PI i (y) (I(x ′ ) -I i (y ′ )) 2 2N β σ2 . (3)
This weighted graph w i is a representation of non-local similarities between the input image I and the image I i of the textbook (see Figure 2).

The assumption on which the proposed method relies is the following one: if patches of the input image I are locally similar to patches of the anatomy textbook, they should have a similar label. The label propagation procedure is then performed using the set of graphs {w i } i=1,••• ,n which reflects the local similarities between I and

{I i } i=1,••• ,n .
In the context of denoising [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], Buades et al. assume that every patch in a natural image has many similar patches in the same image. Thus, the graph w is computed between each voxel of the input image (which has led to the term "non-local" algorithm):w = {w(x, y), ∀(x, y) ∈ Ω 2 }. This is because of the assumption that similarities can be found in the entire image I. However, in the context of human brain labeling, the location of the brain structures is not highly variable and it is then not desirable to go through the entire image domain Ω to find good examples for the segmentation of a specific part of a brain. In this work, the graph is computed using a limited number of neighbors: w = {w(x, y), ∀x ∈ Ω, y ∈ N (x)}, where N (x) is the neighborhood of the voxel x. The size of the considered neighborhood N is directly related to the brain variability. The influence of this parameter is evaluated in the Section 3.2.

Pair-wise Label Propagation

For clarity sake, let start with an input image I and a textbook T which contains only one anatomical image and the corresponding labeled image: T = (I, L) (this is the basis of pair-wise label propagation techniques). Using the patch-based approach described previously, the image I can be labeled using the following equation:

∀x ∈ Ω, L(x) = y∈N (x) w(x, y)L(y) y∈N (x) w(x, y) (4)
where L(y) is a vector of [0, 1] M (M is the total number of labels) representing the proportions for each label at the voxel y in the image L (this notation unifies the case where the anatomy textbook contains hard or fuzzy label maps). Thus,

L(x) = (l 1 (x), l 2 (x), • • • , l M (x)
) are the membership ratios of the voxel x with respect to the M labels, such that k l k (x) = 1 and

l k (x) ∈ [0, 1], ∀k ∈ [[1, M ]].
Then, for each voxel x of I, the equation ( 4) leads to a fuzzy labeling of the input image I since L(x) ∈ [0, 1] M . A hard segmentation H of the image I can be obtained by taking the component of L(x) with the highest value. The hardening of the fuzzy label vector L(x) to get a binary label vector H(x) is done as follows:

h k (x) = 1 if l k (x) = L(x) ∞ 0 otherwise. (5)
Now, let us consider a textbook T containing n pairs of images:

T = {(I i , L i ), i = 1, • • • , n}.
The most straightforward approach is to perform n times the pair-wise procedure previously described, leading to n

fuzzy segmentations {L i (x)} i=1,••• ,n (or n hard segmentations {H i (x)} i=1,••• ,n ).
Then, these n label maps can be fused using a classifier combination strategy [START_REF] Kittler | On combining classifiers[END_REF]. Thus, for instance, under the assumption of equal priors and by hardening the fuzzy label maps L i to get a set of hard segmentations {H i } i=1,••• ,n , the final labeling Λ(x) at the voxel x can be obtained by applying the majority voting (MV) rule:

λ k (x) = 1 if n i=1 h ik (x) = max M j=1 n i=1 h ij (x) 0 otherwise. ( 6
)
where h ij is a binary value corresponding to the label j of the image i at the voxel x and Λ

(x) = (λ 1 (x), λ 2 (x), • • • , λ M (x)).
The pair-wise labeling approach using the majority voting rule to fuse the labels is described in Algorithm 1. A version of this pair-wise technique using STAPLE as the final label fusion is also evaluated in the Section 3.

Algorithm 1 Pair-wise labeling method using the majority voting rule inputs: an image I and an anatomy textbook 4) Compute H i (x) using Equation ( 5) end for Compute Λ(x) by aggregating the set of labels {H i } i=1,••• ,n (Equation ( 6)) end for

T = {(I i , L i ), i = 1, • • • , n} ouput: a label image Λ for all x ∈ Ω do for i = 1 to n do Compute L i (x) using Equation (

Group-wise Label Propagation

We propose to study a group-wise combination strategy which takes all the images of the textbook to produce a fuzzy labeling. Indeed, the Equation ( 4) can be applied for an arbitrary number of labeled images:

∀x ∈ Ω, L(x) = n i=1 y∈N (x) w i (x, y)L i (y) n i=1 y∈N (x) w i (x, y) . (7)
Again, this leads to a fuzzy labeling which can be thresholded to obtain a hard labeling by taking the maximum of each vector L(x):

h k (x) = 1 if l k (x) = L(x) ∞ 0 otherwise. (8)
In this case, the final segmentation Λ can be set to L (resp. H) if a final fuzzy (resp. hard) labeling is desired. Unlike the pair-wise approach, the weight of each label map L i is automatically set by using the local patch-based similarity measure (see Equation ( 3)) with the input image I. There is no need to apply a classifier combination strategy. This is described in Algorithm 2.

Algorithm 2 Group-wise labeling method inputs: an image I and an anatomy textbook

T = {(I i , L i ), i = 1, • • • , n} ouput: a label image Λ for all x ∈ Ω do Compute L(x) using Equation (7) Compute H(x) using Equation (8) Λ(x) ← H(x) end for

Pointwise and Multipoint Estimation

In the previous section, according to the classification proposed by Katkovnik et al. in [17], the patchbased labeling techniques (both pair-wise and group-wise) provide pointwise estimates (Equations ( 4) and ( 7)). This means that the pair-wise and group-wise methods estimate one label vector L for every voxel x. However, since the patch-based similarity measure is the core of the proposed labeling methods, one can obtain a label patch estimate at each considered voxel. This is called a multipoint label estimator. In contrast to the pointwise estimator, a multipoint estimator gives the estimate for a set of points (in our case, a patch).

For instance, a multipoint estimate for the group-wise labeling method is given by:

∀x ∈ Ω, P L (x) = n i=1 y∈N (x) w i (x, y)P Li (y). ( 9
)
where P Li (y) is a 3D patch of the label map L i centered at the voxel y, and P L (x) is one multipoint label estimate centered at the voxel x. Thus, by going through the entire image domain Ω, one label patch estimate is obtain for each voxel. Since there is an overlap between all these label patches (each containing N voxels), for every voxel x ∈ Ω, we have N label estimates. Indeed, the label patch estimate of each voxel y which belongs to the patch P (x) contributes to the final label estimate of the voxel x. These N estimates can then be aggregated using a combination classifier. In this work, we have used the majority voting rule to fuse these N estimates. As suggested in [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF], it is possible to speed-up the algorithm by considering a subset Ω ⋆ of the image domain Ω with the constraint that there is at least one label estimate for each voxel. For instance, by using patches of size 3 × 3 × 3 voxels, one can consider only voxels with even spatial index (in each dimension) which leads to a speed-up of 8. Such an approach (denoted as fast multipoint in this article) is evaluated in the Section 3. Moreover, for pointwise and multipoint estimation, based on the mean and variance of patches [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF], a voxel preselection can be used to avoid useless computation.

Finally, several patch-based label propagation algorithms can be derived from the proposed framework depending of the following choices: pair-wise (PW) vs group-wise (GW), point-wise vs (possibly fast) multipoint estimates, majority voting (MV) rule vs STAPLE.

Experiments

Evaluation Framework

In this work, experiments have been carried out on three publicly available image datasets: Brainweb [START_REF] Aubert-Broche | A new improved version of the realistic digital brain phantom[END_REF], the Internet Brain Segmentation Repository (IBSR) database and the NA0 database developed for the Non-rigid Image Registration Evaluation Project (NIREP). These three complementary image datasets provide different challenges for label propagation techniques: Brainweb dataset is used to mainly evaluate the separation power between 2 principal brain structures (white matter and gray matter), IBSR is a well known dataset for segmentation algorithms evaluation since it contains 32 brain structures (white matter, cortex, internal gray structures such as hippocampus, caudate, thalamus, putamen, etc.), and the NA0 dataset which has been originally created for registration evaluation purpose, provide also a good evaluation framework for cortical parcellation algorithms.

The Brainweb dataset1 is a set of 20 T1-weighted simulated data with these specific parameters: SFLASH (spoiled FLASH) sequence with TR=22ms, TE=9.2ms, flip angle=30 deg and 1 mm isotropic voxel size. Each anatomical model consists of a set of 3-dimensional tissue membership volumes, one for each tissue class: background, cerebro-spinal fluid (CSF), gray matter (GM), white matter (WM), fat, muscle, muscle/skin, skull, blood vessels, connective (region around fat), dura matter and bone marrow. Example images from the Brainweb dataset are shown in Figure 3 (first column). For the IBSR dataset2 , the MR brain data sets and their manual segmentations are provided by the Center for Morphometric Analysis at Massachusetts General Hospital. It contains 18 images of healthy brains and the corresponding segmentation of the whole brain into 32 structures. Example images from the IBSR dataset are shown in Figure 3 (second column). The following pre-processing has been applied on IBSR images: 1) N3-based bias field correction using MIPAV3 , 2) ITK-based histogram matching using 3DSlicer 4 , 3) affine registration using ANTs5 [START_REF] Avants | Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[END_REF]. For each image to segment, all the images of the anatomy textbook are registered to this current image. Then, the label maps are also transformed (using a nearest neighbor interpolation). A total of 18*17 affine registrations have been performed.

The evaluation database NA06 consists of a population of 16 annotated 3D MR image volumes corresponding to 8 normal adult males and 8 females. The 16 MR data sets have been segmented into 32 gray matter regions of interest (ROIs) (see Figure 3 (third column). MR images were obtained in a General Electric Signa scanner operating at 1.5 Tesla, using the following protocol: SPGR/50, TR 24, TE 7, NEX 1 matrix 256 × 192, FOV 24 cm. 124 contiguous coronal slices were obtained, with 1.5 or 1.6 mm thick, and with an interpixel distance of 0.94 mm. Three data sets were obtained for each brain during each imaging session. These were coregistered and averaged post hoc using Automated Image Registration (AIR 3.03, UCLA) 7 . The final data volumes had anisotropic voxels with an interpixel spacing of 0.7 mm and interslice spacing of 1.5 -1.6mm. All brains were reconstructed in three dimensions using Brainvox8 . Before tracing ROIs, brains were realigned along a plane running through the anterior and posterior commissures (i.e., the AC-PC line); this ensured that coronal slices in all subjects were perpendicular to a uniformly and anatomically defined axis of the brain. The following pre-processing has been applied on NA0 images: 1) N3-based bias field correction using MIPAV, 2) ITK-based histogram matching using 3DSlicer.

As it is usually used for label propagation method evaluation, for each aforementioned dataset, a leaveone-out study is performed. Each image is separately selected as the image to be segmented. Using the provided segmentations of the remaining images as the anatomy textbook, label propagation is performed to obtain a segmentation of the considered image. In all cases, the Dice Index (DI) overlap is used as a segmentation quality measure:

DI = 2.T P 2.T P + F P + F N ( 10
)
where T P is the number of true positives, F P is the number of false positives and F N , the number of false negatives. The Dice index is computed for hard segmentations.

Neighborhood N , Patch P and β

Experiments were carried out on the Brainweb dataset to determine the influence of the size of the neighborhood N (which can be viewed as a local search area for similar patches), the size of the patches P and the smoothing parameter β (Equation 2). A pair-wise label propagation approach with the MV rule to fuse label maps has been used to determine the optimal parameter set. We found that the highest overlap value for both gray matter and white matter is obtained with the smallest size of patches (i.e. 3 × 3 × 3 voxels), the lower value of β (0.5) and large neighborhood N (the dice index stabilizes at a neighborhood size of 11 × 11 × 11 voxels). These findings are similar to some results obtained for image denoising [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF] or image reconstruction [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF]. β is a noise level dependent parameter. While the lowest value of β provides the best overlap results on Brainweb images (which are almost noiseless), we have found that a value of 1 is more appropriate for in-vivo brain MR images. In the next sections, the following setting has been used: a patch size of 3 × 3 × 3 voxels, a neighborhood size of 11 × 11 × 11 voxels, and β = 1.

Influence of the size of the anatomy textbook

The evolution of the mean Dice index with respect to the size of the anatomy textbook has also been studied. This experiment has been carried out on the IBSR dataset, using up to 30 image permutations for each image to label. The behavior of the fast multipoint patch-based approach is compared with a non-rigid registrationbased method (using ANTs [START_REF] Avants | Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[END_REF], with the following command line: ANTS 3 -m PR[target.nii, source.nii, 1, 2] -i 100x100x10 -o output.nii -t SyN[0.25] -r Gauss[3,0] and the majority voting rule as label fusion strategy). In both cases, it has been found that increasing the size of the anatomy textbook provides higher overlap values (which stabilize quickly around 8 images). Figure 4 shows the evolution of the mean Dice index with respect to the number of images contained in the anatomy textbook. For all the results shown in the next sections, the maximum number of images has been used.

Pointwise vs blockwise

Comparison of pair-wise (using the majority voting rule for label fusion) vs group-wise approaches is shown in Table 1 for the Brainweb images, in Table 2 for IBSR and in Table 6 for NA0-NIREP. We have also compared three types of estimators: pointwise vs multipoint vs fast multipoint. As mentioned in [START_REF] Katkovnik | From local kernel to nonlocal Multiple-Model image denoising[END_REF] ("in multipoint image estimation a weighted average of a few sparse estimates is better than single sparse estimate alone"), the results obtained on the three datasets suggest that the use of a multipoint estimator leads to the best segmentation (i.e. the highest overlap measures). Finally, the pair-wise approach and the group-wise approach lead to very similar results.

Aggregation strategy

For the label fusion step, the majority voting rule is compared to the simultaneous truth and performance level estimation (STAPLE) [START_REF] Warfield | Simultaneous truth and performance level estimation (STA-PLE): an algorithm for the validation of image segmentation[END_REF]. STAPLE estimates the performance of each classifier iteratively and weights it accordingly, relying on an expectation-maximization (EM) optimization approach. We used the implementation provided by T. Rohlfing 9 . As already reported in [START_REF] Artaechevarria | Combination strategies in multi-atlas image segmentation: application to brain MR data[END_REF], STAPLE based fusion rule does not necessary lead to higher Dice coefficients compared to the majority voting rule.

Comparison with previously reported results

The results obtained on the IBSR dataset are reported in Table 5 (highest Dice coefficients are in bold). In this table are also shown results reported in the literature, including non-rigid registration-based approaches and segmentation techniques. This experiment clearly shows that the proposed patch-based framework is very competitive with respect to recently published methods.

Computational Time

Table 7 shows the comparison of the computational time for a non-rigid registration-based approach (using ANTs [START_REF] Avants | Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[END_REF]) and the proposed patch-based techniques. We choose ANTs as it has been ranked as one of the best among fourteen non-rigid registration algorithms [START_REF] Klein | Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration[END_REF]. The use of pointwise or multipoint estimator leads to higher computational time than using non-rigid registration. However, the fast multipoint approach provides a very time efficient algorithm. Moreover, the multi-threading of the proposed label propagation method is particularly adapted due its blockwise nature. As eight processors were used for our experiments, the computational time with multi-threading is about 8 times smaller. In this case, the propagation of one label image with the fast multipoint estimator can be performed in about 1mn.

Visual Evaluation

Figures 5 and6 show visual comparison between the proposed patch-based approach with a non-rigid registration-based approach (ANTs + majority voting). The first row of this figure corresponds to the mean image obtained using the T1-weighted MR images of the anatomy textbook. Such mean image is one way of visual evaluation used for non-rigid registration algorithms. The proposed patch-based method can also provided such mean image by simply using intensities (instead of labels) of the images contained in the textbook. The second row corresponds to obtained brain segmentation. This figure illustrates well the fact that the proposed technique allows one to get a better delineation of brain structures (and especially the cortex for instance).

Combining patch-based strategies and registration-based techniques

The main purpose of the proposed patch-based method is to avoid long time computation due to nonrigid registration. However, one can wonder if the two techniques should cooperate to reach better results.

To one side, one of the advantages of the patch-based approach is to possibly consider multiple examples within the same image of the anatomy textbook. The consequence of the local search windows use is to relax the one-to-one constraint usually involved in non-rigid registration-based label propagation approaches. Experiments on the IBSR dataset have clearly shown that the use of the patch-based approach leads to very satisfactory segmentation results for brain structures with sharp contrast. On the other side, the use of one-to-one constraint leads to a topologically regularized segmentation. We have experimentally observed that this aspect is important for cortex parcellation. Indeed, the delineation of cortical areas does not rely on intensity contrast within the cortex. It appears then that the intensity similarity assumption which the patch-based approach relies on is not sufficient to provide the highest quality cortical parcellation. Table 8 shows the mean Dice index obtained on the NA0-NIREP dataset using the multipoint pairwise patch-based approach, a non-rigid registration-based technique using ANTs (using majority voting or STAPLE as fusion rule) and the mixed approach. This latter technique simply consists in fusing the labeling obtained with a non-rigid registration algorithm, with the hard cortical mask obtained with a patch-based technique. The mixed approach clearly leads to higher Dice index compared to the two other techniques (which provide similar overlap scores separately). It is also interesting to note that in the combining approach, STAPLE seems to provide a better label fusion framework than the majority voting rule. Figure 7 presents a cortex parcellation result using the combined approach with STAPLE.

Discussion

Label propagation is a versatile image segmentation technique which can be applied to a large variety of images. In this work, we focus on the development of a new label propagation framework applied to automated human brain labeling. A patch-based label propagation has been proposed whose purpose is to relax the one-to-one constraint existing in non-rigid registration-based techniques. Indeed, this image similarity-based approach can be seen as a one-to-many block matching technique. It allows the use of several good candidates (i.e. the most similar patches) to estimate label patches. Several patch-based algorithms have been derived depending on the patch aggregation strategies and the label fusion rules (if needed). Comparison with previously proposed non-rigid registration-based methods on publicly available in-vivo MR brain images has shown a great potential of this method.

The proposed approach is also related to supervised learning methods which try to learn the link between image intensity patches and label patches. In our work, instead of learning the link between image intensity and labels, the label estimation relies on a weighted graph which represents directly the similarity between intensity of the input image and intensity of the anatomy textbook images. As mentioned in the Introduction, the work presented in this article shares similarities with the one done by Coupé et al. [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF]. Indeed, they have proposed a patch-based strategy relying on a group-wise technique for label fusion. The main methodological differences concern the smoothing parameter in the definition of patch weights and the voxel pre-selection. About validation, they have investigated one aggregation strategy (group-wise) on two regions of interest: hippocampus and lateral ventricles (on a non-publicly available dataset).

The experiments on three freely available T1-weighted MR image datasets (Brainweb, IBSR, NA0-NIREP) have shown that the proposed framework can lead to high quality segmentation of a large number of brain structures. It has been shown that CPU time consuming non-rigid registration steps can be avoided. The proposed patch-based technique does not require an accurate correspondence between the input image and the anatomy textbook. The usual assumption of one-to-one mapping is relaxed by using local search windows. Further work may investigate the use of more complex similar patch search steps. For instance, using a rotation-invariant patch match is one of the possible extensions of the proposed method (in particular when considering large patches), while keeping in mind that it might increase the CPU time. However, the results obtained with translation-based patch match are already very competitive with respect to existing methods. Some possible methodological extensions are local adaptivity of patch size and search window size, choice of kernel function f and the use of graph operators.

As the algorithm relies on the assumption of similar intensities between the input image and the anatomy textbook, the segmentation accuracy may depend on the contrast of the structure to segment or intensity variations (such as intensity bias). Intensity variations might disturb the search process of similar patches, which would lead to a final labeling obtained with less relevant patch examples. The sensitivity of intensity variations can be captured by the smoothing factor beta. However, using a fixed global value of β in conjunction with common correction techniques such as intensity bias correction and histogram matching, the proposed algorithm provides very satisfactory results. Also for the IBSR data, it can be noticed that the segmentation of the pallidum is not as good as the segmentation of other brain structures. The reason for failure is that this structure has no clear contrasted boundaries in MR T1-weighted images. Thus, contrary to registration-based techniques which can propagate spatial relations of the structures, a pure intensity-based method such as the one proposed in this work cannot provide a satisfactory delineation of the pallidum. This assertion is corroborated by the experiments performed on the NA0-NIREP dataset. Because of the key assumption of intensity similarity, a non-rigid registration-based label propagation approach outperforms the proposed patch-based method. In this context of cortex parcellation, we have shown that these two techniques can be complementary. In particular, the joint use of non-rigid registration and patch-based strategy can significantly improve the segmentation result. Thus, further work may explore the incorporation of a regularization term and prior information such as spatial relations or a topological atlas into the patchbased framework.

Contrary to registration-based techniques, especially for the group-wise approach using a patch selection, if there is no correspondence between the input and the anatomy textbook (i.e. the image similarity is lower than the threshold used in the patch selection step), no label is propagated. This prevents an incorrect label being introduced in the estimation of the final label image. Thus, further work is to evaluate the robustness of this label propagation method to images which contain lesions or tumors, and to use this patch-based technique to potentially detect these pathological patterns. Another research direction concerns the use of a database of images as IBSR to segment pathological images. Such a "cross-site" study will analyze the robustness of the method in the presence of large anatomical variability, including the initial affine registration, intensity and shape variations, and the research of examples in the anatomy textbook.

Finally, as mentioned in other works [START_REF] Wu | Optimum template selection for atlas-based segmentation[END_REF][START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF][START_REF] Van Rikxoort | Adaptive local multi-atlas segmentation: application to the heart and the caudate nucleus[END_REF] and as shown in the experiments performed on the Brainweb dataset, the label image can be estimated using a subset of the anatomy textbook (the observed convergence rate is about 10 example images). This remark leads to the issue of atlas selection, i.e. identifying a subset of representative examples in the textbook with respect to the input image. This can be particularly important when using a large anatomy textbook to avoid useless heavy computational burden.

Figure 1 :

 1 Figure 1: Principle of registration-based label propagation methods. The input data (shown with green borders) are an anatomy textbook (i.e. a set of n anatomical images with the corresponding label maps), and one anatomical image I. The set of anatomical images of the textbook is (non-linearly) registered to the input image I, and each label map is deformed with respect to the estimated transformation h i . The final image segmentation (shown with red borders) is then obtained by fusing all the deformed label maps (shown with blue borders).

Figure 2 :

 2 Figure 2: Weighted graph building. The set of graphs {w i } i=1,••• ,n is a representation of non-local interactions between the input image I and the images {I i } i=1,••• ,n of the textbook. N i (x) is the neighborhood of the voxel x in the image I i .

Figure 3 :

 3 Figure 3: MR image datasets used for the evaluation and the corresponding segmentation. First column: Brainweb database, second column: IBSR database, third column: NA0-NIREP database.

Figure 4 :

 4 Figure 4: Mean dice index with respect to the number of images used for label propagation Left: patchbased approach, right: registration-based approach using ANTS (using, for each case, majority voting to fuse labels). (a) : left white matter, (b): left cortex, (c): left thalamus, (d): left hyppocampus.

Figure 5 :

 5 Figure 5: Visual evaluation of brain segmentation results (IBSR, image #7). First row: T1-weighted images, second row: corresponding segmentation, third row: misclassified voxels. Left: ground-truth, middle: patchbased technique, right: non-rigid registration-based approach (using ANTs).

Figure 6 :

 6 Figure 6: Surface rendering of segmentation results (IBSR, image #7). Left: ground-truth, middle: patchbased technique, right: non-rigid registration-based approach (using ANTs).

Figure 7 :

 7 Figure 7: Visual evaluation of cortex parcellation results (NA0-NIREP, image #6). First row: ground truth, second row: combination technique (ANTs + patch + STAPLE).

Figure 8 :

 8 Figure 8: Surface rendering of cortex parcellation results (NA0-NIREP, image #6). First row: ground truth, second row: combination technique (ANTs + patch + STAPLE).

Table 1 :

 1 Comparison (using the mean Dice index) of several possible combinations of patch-based approaches (using optimal global parameter setting) on the Brainweb dataset. For the pair-wise results, labels have been fused using the majority voting rule. Highest overlap rates are shown in bold.

	brain structure		pair-wise			group-wise
		pointwise multipoint fast multipoint	pointwise multipoint fast multipoint
	CSF	0.87	0.88	0.87	0.87	0.86	0.86
	Gray Matter	0.95	0.95	0.94	0.95	0.93	0.94
	White Matter	0.95	0.95	0.94	0.95	0.94	0.94
	Fat	0.83	0.87	0.85	0.83	0.85	0.84
	Muscle	0.66	0.70	0.67	0.66	0.67	0.66
	Muscle/Skin	0.93	0.94	0.93	0.93	0.93	0.93
	Skull	0.92	0.93	0.91	0.92	0.92	0.92
	vessels	0.61	0.59	0.56	0.61	0.56	0.57
	around fat	0.37	0.36	0.30	0.37	0.34	0.36
	dura matter	0.51	0.53	0.47	0.51	0.49	0.50
	bone marrow	0.81	0.83	0.79	0.81	0.81	0.80

Table 3 :

 3 Comparison (using the mean Dice index) of aggregation strategies (majority voting vs STAPLE) on the IBSR dataset, using a multipoint estimation approach. When structures are separated into left/right, two overlap scores are reported. Highest overlap rates are shown in bold.

http://mouldy.bic.mni.mcgill.ca/brainweb

http://www.cma.mgh.harvard.edu/ibsr

http://mipav.cit.nih.gov

http://www.slicer.org

http://www.picsl.upenn.edu/ANTS

http://www.nirep.org

http://bishopw.loni.ucla.edu/air5

http://www.nitrc.org/projects/brainvox

http://www.stanford.edu/ rohlfing/software/index.html

Acknowledgment

The research leading to these results has received funding from the European Research Council under the European Communitys Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 207667). This work is also funded by NIH/NINDS Grant R01 NS 055064.

2009 [START_REF] Artaechevarria | Combination strategies in multi-atlas image segmentation: application to brain MR data[END_REF] 2002 [START_REF] Fischl | Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain[END_REF] 2009 [START_REF] Khan | Robust atlas-based brain segmentation using multi-structure confidence-weighted registration[END_REF] 2010 [START_REF] Lötjönen | Fast and robust multi-atlas segmentation of brain magnetic resonance images[END_REF] 2010 [START_REF] Sdika | Combining atlas based segmentation and intensity classification with nearest neighbor transform and accuracy weighted vote[END_REF] 2007 [


Samy Abbes

Distributed Markov Processes

come   L'archive ouverte pluridisciplinaire

The Algebraic Model

Sites. Local and Global States. Transitions. We consider n ≥ 2 sites, indexed 1, . . . , n. Each site i has a set of local states denoted X i . The X i 's are not supposed to be pairwise disjoint. We denote by x i a local state that belongs to X i . We assume that "∅" is a special symbol that does not belong to any X i .

We define a global state of the system as any n-tuple (x 1 , . . . , x n ) ∈ X 1 × • • • × X n . It will be convenient to adopt a special notation for the set of global states, so we put

X = X 1 × • • • × X n .
To each local state x ∈ n i=1 X i , we associate a n-tuple t = (t 1 , . . . , t n ) as follows:

for j = 1, . . . , n, t j =

x, if x ∈ X j , ∅, the empty word, otherwise.

(1)

So for instance, if x belongs to S 1 and S 2 only, we have t = (x, x, ∅, . . . , ∅).

Definition 1.1 (transitions and ressources) We call transition any n-tuple t = (t 1 , . . . , t n ) defined as in Eq. ( 1), and we denote by E the set of transitions. The resources of transition t are defined as those indices i ∈ {1, . . . , n} such that t i = ∅. The set of resources of t is denoted by ρ(t), so that

ρ(t) = {i : 1 ≤ i ≤ n, t i = ∅} = {i : 1 ≤ i ≤ n, t i ∈ X i }. ( 2
)
Definition 1.2 (independence relation) We define the independence relation on transitions by declaring two transitions independent if their resources are disjoint. The independence relation is denoted with symbol " / / ", so we have:

∀u, v ∈ E u / / v ⇐⇒ ρ(u) ∩ ρ(v) = ∅.

Sequences of Transitions and Finite

Trajectories. Let S = E * be the set of finite sequences of transitions. Then S is equipped with two structures: the prefix ordering, denoted ≤, and the concatenation on sequences, denoted by the dot •. The latter identifies S with the free semigroup generated by E.

The two structures are related as follows:

∀u, v ∈ S u ≤ v ⇐⇒ ∃z ∈ S v = u • z.
(3) Furthermore, they satisfy the following compatibility relation:

∀u, v, w ∈ S u ≤ v ⇒ w • u ≤ w • v. ( 4
)
Definition 1.3 (Adjacency relation) Say that two sequences u, v ∈ S are immediately adjacent if one can deduce one from the other by inverting two neighbor transitions t, t ′ such that t / / t ′ . The symmetric and transitive closure of the immediate adjacency relation thus defined is called the adjacency relation, denoted by ∼.

Let S denote the quotient set S = S/ ∼. Then the partial order and the semigroup structures on S induce similar structures on the quotient S, which satisfy both properties (3) and ( 4), in S. We denote by ǫ the neutral element in S, which is also the least element of S, and is the image in S of the empty word.

Definition 1.4 Elements of S are called finite trajectories.

Projection of Trajectories on Local Sites.

For each site i = 1, . . . , n, let θ i : E → (X i ) * denote the i th natural projection. An inspection of Eq. ( 2) defining resources of transitions shows that θ i satisfies the following:

∀i = 1, . . . , n θ i (t) = ∅ ⇐⇒ i / ∈ ρ(t), θ i (t) ∈ X i ⇐⇒ i ∈ ρ(t). (5) Lemma 1.1 For each i = 1, . . . , n, the i th projection θ i : E → (X i ) * satisfies t / / t ′ ⇒ θ i (t) • θ i (t ′ ) = θ i (t ′ ) • θ i (t) . ( 6
)
It follows from (6) that the semigroup extension θ i : S → (X i ) * induces an application that we still denote by

θ i : S → (X i ) * .
Intuitively, if u is a finite trajectory, then θ i (u) represents the local action of u on site i.

Proposition 1.1 1.-The application θ i is a morphism of semigroups. 2.-For u, v ∈ S, we have θ 1 (u), . . . , θ n (u) = θ 1 (v), . . . , θ n (v) ⇒ u = v.
(7) It is to be noted that the mapping

θ 1 × • • • × θ n : S → (X 1 ) * × • • • × (X n ) * ,
which is injective by virtue of Proposition 1.1 above, is not surjective in general.

Action of Finite Trajectories on Global States.

If A is any set, and if s ∈ A * is a nonempty word on A, denote by γ(s) the last letter of s. The semigroup A * comes equipped with a natural right action on A, defined as follows:

∀a ∈ A, ∀s ∈ A * a • s = a, if s is empty, γ(s), if s non empty. ( 8
)
This is indeed a right action, that is, we have

∀a ∈ A ∀s, s ′ ∈ A * a • (s • s ′ ) = (a • s) • s ′ .
Next we extend the above action to a right semigroup action of S on the set X of global states, as follows. Consider first the extension componentwise of the above action to an action of S on X, and observe that the following holds:

∀x ∈ X ∀t, t ′ ∈ E t / / t ′ ⇒ x • (t • t ′ ) = x • (t ′ • t).
It follows that the right action of S on X induces a right action of the quotient semigroup S on X. By construction, the action commutes with the projections π i : X → X i in the following way: for x ∈ X and u ∈ S and putting y = x • u and y i = x i • θ i (u), then:

x u ----→ y π i     π i x i θ i (u) ----→ y i
Histories as Limiting Trajectories. Poset S lacks lubs (least upper bounds). But it has a natural completion which is constructed as follows.

Consider two increasing sequences u

= (u k ) k≥0 and v = (v k ) k≥0 in S. Define u v whenever ∀k ≥ 0 ∃p ≥ 0 u k ≤ v p . ( 9
)
This is the so-called Egli-Miller order on sequences. Consider also the equivalence relation on increasing sequences in S defined by

u ≡ v ⇐⇒ (u v) ∧ (v u). ( 10
)
A limit element is represented by some increasing sequence u of elements in S ; any two increasing sequences u and v are identified provided they satisfy u ≡ v. Denote by (S, ≤) the quotient poset, the partial order of which is inherited from the order on sequences defined in (9).

There is an obvious embedding of posets S → S. Furthermore, S has the property that any increasing sequence in S has a lub in S. And for any element x ∈ S there is a maximal element ω ∈ S such that x ≤ ω. Definition 1.5 We say that a limiting element ω ∈ S is a history if ω is a maximal element of S. The set of histories is denoted by Ω.

Local Projections of Histories. For each i = 1, . . . , n, the application θ i : S → (X i ) * is increasing since it is a semigroup morphism according to Prop. 1.1. Denote by (X i ) * the set whose elements are either finite or infinite sequences of elements of X i . For u = (u k ) k≥0 and v = (v k ) k≥0 two increasing sequences of elements of S satisfying u ≡ v as in (10), the sequences θ i (u k ) k≥0 and θ i (v k ) k≥0 are both increasing in (X i ) * since θ i is increasing, and satisfy furthermore

θ i (u k ) k≥0 ≡ θ i (v k ) k≥0 ,
where the ≡ equivalence is now taken as the equivalence relation relative to the Egli-Milner order on increasing sequences in (X i ) * . They define thus a unique limiting sequence of elements in X i , finite or infinite, that is to say, a unique element in (X i ) * . We will denote this element by θ i (u), for u the corresponding class in S. Since the data

θ 1 (u), . . . , θ n (u) ∈ (X 1 ) * × • • • × (X n ) * (11)
uniquely determines u ∈ S, we will most often represent an element u ∈ S by the corresponding tuple (11). In particular the set Ω of maximal elements in S identifies with a subset of (X

1 ) * × • • • × (X n ) * .

Concatenation of Trajectories with Histories. The concatenation

S × S → S,

seen as a left action of S on itself, extends to a left action of S on S. Indeed, if u = (u k ) k≥0 and v = (v k ) k≥0 are two increasing sequences in S such that u ≡ v (in the sense of (10)), and if x is some fixed element in S, then the sequences x • u and x • v defined by

x • u = (x • u k ) k≥0 , x • v = (x • v k ) k≥0 , are increasing and satisfy x • u ≡ x • v. It is readily seen that • : S × S → S is indeed a left action of S on S, that is to say, that x • (y • u) = (x • y) • u for
x, y ∈ S and u ∈ S; and that it restricts to a left action of S on Ω.

Furthermore, if x ∈ S is finite and if u ∈ S is such that x ≤ u there is a unique element v ∈ S such that u = x • v. This element v is • finite if and only if u is finite;
• maximal if and only if u is maximal.

Shadows.

For any finite trajectory u ∈ S, we define the shadow of u as the subset of Ω denoted by ↑ u and defined by

↑ u = {ω ∈ Ω : u ≤ ω} .
By analogy with the theory of infinite product spaces in measure theory, ↑ u can be seen as the elementary cylinder of base u.

For u ∈ S a fixed finite trajectory, the mapping

Φ u : Ω →↑ u , ω → u • ω , ( 12
)
where the concatenation u • ω is defined as above, is a bijection.

The Probabilistic Model

General Notions. As a subset of (X

1 ) * × • • • × (X n ) * , the set Ω is nat- urally equipped with a σ-algebra, that we denote by F, inherited from the Borel σ-algebra on (X 1 ) * × • • • × (X n ) * .
The σ-algebra F is generated by the collection of shadows ↑ u, for u ranging over S . It is a result of classical measure theory that if two probabilities P and P ′ on (Ω, F) satisfy P (↑ u) = P ′ (↑ u) for all u ∈ S, then P = P ′ .

Assume P is some probability on (Ω, F), and let u ∈ S be some finite trajectory. It is straightforward to check that the bijection Φ u : Ω →↑ u defined in ( 12) is not only a bijection, but is also bi-measurable. Therefore, if P (↑ u) > 0, we shall consider the probability P u on (Ω, F), image under Φ -1 u :↑ u → Ω of the conditional probability P ( • | ↑ u). The probability measure P u is characterized by its values on shadows ↑ w, for w ranging over S, and given by

∀w ∈ S, P u (↑ w) = 1 P (↑ u) P ↑ (u • w) . ( 13
)
Probabilistic Processes. Instead of merely considering a probability measure on Ω, we will consider a finite family (P x ) x∈X 0 of probability measures, indexed by some subset X 0 of X, and all defined on Ω. Intuitively, P x describes the behavior of the system starting from the global state x.

We will denote by P such a family of probabilities, and call it a probabilistic process.

Reachable and Locally Reachable States.

Definition 2.1 (reachable states) Let P = (P x ) x∈X be a probabilistic process defined on n sites. A global state y ∈ X is reachable from some

global state x ∈ X if ∃u ∈ S, y = x • u, P x (↑ u) > 0.
A weaker notion, often more useful, is that of locally reachable state.

Definition 2.2 (locally reachable state) Let P = (P x ) x∈X be a probabilistic process defined on n sites. Let i be a given site, and let s ∈ X i be a local state of site i. We say that s is locally reachable from some global state x ∈ X if there exists a global state y ∈ X such that:

1. y is reachable from x, in the sense of Def. 2.1; and 2. s = θ i (y).

Remark 2.1 If n = 1, the notion of local state coincides with that of local state; and reachable and locally reachable states are two identical notions.

Distributed Markov Processes.

Definition 2.3 (DMP) A probabilistic process P = (P x ) x∈X 0 is said to be a distributed Markov process (DMP) if for all x ∈ X 0 , and for all u ∈ S such that P x (↑ u) > 0, the two following assertions are true:

1. x • u ∈ X 0 , and 2. (P x ) u = P x•u , where (P x ) u is the probability measure on (Ω, F) defined as in ( 13), with respect to probability P x .

The following proposition gives a criterion to insure that a probabilistic process is indeed a DMP.

Proposition 2.1 A probabilistic process P = (P x ) x∈X is a DMP if and only if the value of the following quantities:

1 P x ↑ u P x ↑ (u • t) , for (u, x, t) ∈ S × X × E, ( 14
)
only depend on the transition t, and on the resulting state x • u.

Characteristic Numbers of a Distributed Markov Process. Consider a DMP P-provided that it exists, to be proved later. Then define the following family of real numbers:

P x (↑ t) , for x ∈ X and t ∈ E. ( 15
)
We call it the family of characteristic numbers associated with the DMP.

Proposition 2.2 If P and Q are two DMPs that induce the same families of characteristic numbers (15), then P = Q.

Concurrency Equations for Characteristic Numbers. Not any family of numbers p x (t), x ∈ X, t ∈ E identifies with the characteristic numbers of some DMP. Sadly, given a family p x (t), x ∈ X, t ∈ E of nonnegative numbers, it is difficult to find a convenient necessary and sufficient condition for the existence of a DMP P =

(P x ) x∈X such that ∀x ∈ X, ∀t ∈ E, p x (t) = P x (↑ t) .
However, a first condition is given in the following proposition.

Proposition 2.3 Let p x (t), x ∈ X, t ∈ E be the family of characteristic numbers of some DMP. Then the following equations are satisfied:

∀x ∈ X 0 , ∀a, b ∈ E, a / / b ⇒ p x•a (b) = p x•b (a). ( 16
)
The equations ( 16), that reflect the partial order structure of trajectories, we call the concurrency equations. What is missing are normalization equations that would reflect the structure of the synchronizations.

Stopping Times and the Strong Markov Property for DMPs

Stopping Times

Let N = N ∪ {∞} n . This set will describe the "instants" of a system distributed over n sites. The analogy with the usual time must be used carefully, since in particular not all instants make sense; more precisely, the instants which make sense for a given history are random, that is to say, depend on the given history.

Recall that S denotes the set of finite or infinite trajectories. If u ∈ S is any trajectory, we define the length of u as follows:

|u| = |θ 1 (u)|, . . . , |θ n (u)| ∈ N , ( 17
)
where | • | denotes the usual length of words, finite or infinite. We say that u ∈ S has finite length whenever u is a finite trajectory, denoted by |u| < ∞.

Notation 3.1 Let ω ∈ Ω and t ∈ N with t = (t 1 , . . . , t n ). We denote by ω t the n-tuple

ω t = (ω 1 t , . . . , ω n t ) ∈ (X 1 ) * × • • • × (X n ) * such that: ∀i = 1, . . . , n ω i t ≤ θ i (ω), ω i t = t i .
It must be noted that, in general, ω t is not a prefix of ω.

Definition 3.1 (Random times and stopping times) Let T : Ω → N be an arbitrary mapping. We denote by the short notation ω T the n-tuple ω

T (ω) ∈ (X 1 ) * × • • • × (X n ) * defined for ω ∈ Ω.
1. We say that T is a random time if ω T ≤ ω for all ω ∈ Ω.

2. If T is a random time we say that T is a stopping time if furthermore the following property holds:

∀ω, ω ′ ∈ Ω ω ′ ≥ ω T ⇒ ω T = ω ′ T .
Proposition 3.1 Let T : Ω → N be a stopping time.

1. Then T and ω T are two measurable mappings, N being equipped with the discrete σ-algebra.

2. The σ-algebras generated by T and by ω T coincide. This σ-algebra F T is also characterized as follows:

∀A ∈ F A ∈ F T ⇐⇒ ∀ω, ω ′ ∈ Ω ω ∈ A, ω ′ ≥ ω T ⇒ ω ′ ∈ A.
Proof. 1. We show that ζ = ω t is measurable. It is enough to show that for every finite trajectory u ∈ S, we have that ζ -1 {u} is a measurable subset of Ω. But, by the property of stopping times, ζ -1 {u} is either empty or ↑ u, so it is measurable in both cases. It follows that T is measurable as well, since the length mapping | • | : S → N is clearly measurable.

Shift Operator and Iteration of Stopping Times

Definition 3.2 (Shift operator) Let T : Ω → N be a stopping time. The shift operator associated with T is the mapping θ T : Ω → Ω , which is only partially defined; if T (ω) < ∞, then θ T (ω) is defined as the unique element of Ω such that ω = ω T • θ T (ω), and θ T (ω) is undefined otherwise.

Definition 3.3 (iterated stopping times) Let T : Ω → N be a stopping time, and let θ T be the associated shift operator. We define a sequence (R k ) k≥0 of mappings Ω → N , called iterated stopping times associated with T , as follows:

R 0 = 0 ∀k ≥ 0 R k+1 = R k + T • θ R k , whenever R k and T • θ R k are finite, ∞,
elsewhere.

Remark 3.1 θ 0 = Id Ω , and R 1 = T .

Proposition 3.2 Let T : Ω → N be a stopping time. Then the sequence (R k ) k≥0 of iterated stopping times is a sequence of stopping times.

Lemma 3.1 Let T : Ω → N be a stopping time, and let (R k ) k≥0 be the associated sequence of iterated stopping times. Let U : Ω → N be another stopping time. Then the mapping V : Ω → N defined as follows:

∀ω ∈ Ω, ω V = {ω T k : ω T k ≥ ω U , k ≥ 0}, if ∃k ≥ 0, ω U ≤ ω T k , ∞, otherwise.
is a stopping time.

The (Strong) Markov Property

Since constant times do not make sense in our context, a weak Markov property is excluded. Fortunately, stopping times have a meaning for us, and therefore a strong Markov property can be expected, as shown by the following theorem.

Theorem 3.1 Let T : Ω → N be a stopping time and let P = (P x ) x∈X 0 be a DMP. Then the following equality hods for all x ∈ X 0 and for every non-negative and measurable function h : Ω → R:

E x h • θ T |F T = E γ(ω T ) (h) , P x -almost surely, (18)
where both members are random variables only defined on {T < ∞}.

Construction of a Class of Distributed Markov Processes

Synchronization of Trajectories

For u a trajectory on n ≥ 2 sites, we denote by θ [1,n-1] (u) the projection of u on the n -1 first sites. In other words, w = θ [1,n-1] (u) is the unique trajectory on n -1 sites such that ∀i = 1, . . . , n -1 θ i (w) = θ i (u).

Definition 4.1 Let w be a finite trajectory on n -1 sites with local state spaces X 1 , . . . , X n-1 , and z be a finite trajectory on a n th site with local state space X n .

1. We say that w and z are compatible if there exists a trajectory u on n sites such that θ [1,n-1] (u) = w, and θ n (u) = z.

2. We say that w and z are sub-compatible if there exists a sub-trajectory w ′ ≤ w such that w ′ and z are compatible. In this case, we put v = inf{w ′ ≤ w : w ′ and z are compatible} and we define the synchronization of w and z as the unique trajectory on n sites, denoted by w z such that

θ [1,n-1] (w z) = v, and θ n (w z) = z.
Definition 4.2 Consider a process defined on n -1 sites with state spaces X 1 , . . . X n-1 , and let X n be another set. Put

S = X 1,n ∪ • • • ∪ X n-1,n ,
to be understood as the set of shared states of site n with some other site.

Let Ω [1,n-1] be the space of histories defined on the n -1 first sites, and let Ω n be the set of trajectories on the n th site. Finally, put

Λ = Ω [1,n-1] × Ω n .
For z ∈ Ω n , we consider T (z) ∈ N ∪ {∞}, defined as the first hitting time of S for z. We keep our usual notation z T to denote the sub-trajectory of z until time T (z). We define the first synchronization time as σ : Λ → N such that:

for λ = (w, z) ∈ Λ, λ σ = ∞,
if w and z T are not sub-compatible, w z, if w and z T are sub-compatible.

Proposition 4.1 The first synchronization time satisfies the following property:

∀λ, λ ′ ∈ Ω [1,n-1] × Ω n λ σ ≤ λ ′ ⇒ λ σ = λ ′ σ . ( 19
)
Of course, similarity with Def. 3.1 is to be noted. Def. 3.1 does not directly apply here since Λ = Ω [n-1] ×Ω n is not our usual space of trajectories on n sites, precisely because we are in the process of building trajectories on n sites.

Synchronization Product: Informal Description

Let Q and M be two DMP, Q over n -1 sites and M over a n th site with set of states X n . Observe that M is just a homogeneous Markov chain on X n .

We will denote by

X [1,n-1] the set X 1 × • • • × X n-1
of global states of process Q. We will use the notation X to denote

X = X [1,n-1] × X n .
The notation X is purposely identical with the one introduced to denote the set of global states on n sites.

We assume that the following assumptions are in force:

1. For any s, s ′ ∈ X n , s ′ is reachable from s with respect to M; 2. For any x ∈ X [1,n-1] , for any site i ≤ n -1 and for any local state s ∈ X i , s is a locally reachable state from x;

3. X i,n = ∅ for all i = 1, . . . , n -1.

Informally, the construction we perform intends to force the synchronization of the two processes described by Q and M, and proceeds as follows. Consider an element λ ∈ Λ = Ω [n-1] × Ω n . Consider the first synchronization time σ described in § 4.1, provided it is finite, yielding λ σ . To insure λ σ is well defined, introduce the event

∆ = {λ ∈ Λ : |λ σ | < ∞} , ( 20
)
and equip λ σ with the conditional law

K (x,y) (λ σ ) = Q x ⊗ M y (λ σ |∆), (21)
where x and y are respectively the initial global states of histories ω and ζ, and where Q x and M y are the probability measures that define the processes Q and M respectively. By construction, the random variable λ σ identifies with a finite trajectory on n sites, that brings the system with n sites to a new global state (x ′ , y ′ ); then repeat this construction, starting from (x ′ , y ′ ) this time, and then concatenate the resulting new random trajectory with λ σ previously obtained. The concatenation of infinitely many of these random trajectories eventually leads to a random history ω ∈ Ω and the law of this random history is precisely the probability measure on Ω that we are seeking.

Construction of the Synchronization Product

We now proceed to the rigorous construction of the synchronization product, informally described above. The first step consists in describing the transition matrix and initial measure of the auxiliary Markov chain that has been described in the previous subsection.

Transition Matrix of the Auxiliary Markov Chain. Consider the countable set C of values of λ σ : Λ → S, equipped with the family of probabilities K (x,y) given by Eq. ( 21); and then the product C ′ = C × X . Elements of C ′ will be noted as (w, z, x, y, ) with w a finite trajectory on n -1 sites, z a finite path on X n , and (x,

y) ∈ X [1,n-1] × X n .
Then we introduce the transition matrix L on C ′ given by:

L(w, z, x, y → w ′ , z ′ , x ′ , y ′ ) = 1 γ(w z)=(x ′ ,y ′ ) × K (x ′ ,y ′ ) λ σ = (w ′ , z ′ ) , (22)
where 1 {•} stands for the characteristic function.

Initial Measures for the Auxiliary Markov Chain. We consider a family (µ (x,y) ) (x,y) of initial measures on C ′ , for (x, y) ranging over the set of global states X 1 × • • • × X n . For (x 0 , y 0 ) ∈ X being given, µ (x 0 ,y 0 ) is defined as follows:

∀(w, z, x, y) ∈ C ′ , µ (x 0 ,y 0 ) (w, z, x, y) = 1 {x=x 0 , y=y 0 } K (x 0 ,y 0 ) (λ σ = w z). ( 23
)
Definition of the Synchronization Product.

Lemma 4.1 If assumptions 1 and 3 above are in force, and if

(w n , z n , x n , y n ) n≥0
is a sample path of the Markov chain with transition matrix L given above, then the concatenation

(w 1 z 1 ) • (w 2 z 2 ) • . . . ( 24
)
is almost surely a history, i.e., a maximal trajectory in the space Ω of trajectories over n sites.

Proof. Assume that (w n , z n ) n≥0 is a sample such that the concatenation ( 24) is not maximal. Put

∀n ≥ 1 r n = (w 1 z 1 ) • (w 2 z 2 ) • . . . • (w n z n ),
and r = (w 1 z 1 ) • (w 2 z 2 ) • . . .

= sup n≥1 r n .
Since r is assumed to be non maximal, there is a site i such that the θ i (r) is of finite length. Since θ i (r) = sup n≥1 θ i (r n ), and since θ i (r) has length i at least, it follows that i < n. Hence, infinitely often, site n has not synchronized with site i although it has non zero probability according to Assumption 3 (stated at the beginning of § 4.2). By the Borel-Cantelli Lemma, all such r together form a subset of zero probability, what was to be proved.

Let Ξ be the natural sample space for Markov chains on C ′ . For (x, y) ∈ X , let π (x,y) be the probability measure on Ξ associated with initial measure µ (x,y) defined in Eq. ( 23) and with transition matrix L defined in Eq. ( 22). According to Lemma 4.1, the concatenation of elements of a sample path ξ ∈ Ξ defines an element ω ∈ Ω. The concatenation defines thus a mapping Φ : Ξ → Ω, which is easily checked to be measurable. Definition 4.3 (synchronization product) Let Q and M be two DMP respectively defined on n -1 site for Q and on a n th site for M. We assume the two assumptions 1 and 2 are in force. Then the synchronization product of Q and M is the probabilistic process on n sites P = P (x,y) (x,y) , with (x, y) ranging over the set X of global states over n sites, and with P (x,y) the probability on (Ω, F) given as the image measure

P (x,y) = Φπ (x,y) ,
where π (x,y) is the probability on the sample space Ξ defined as before, with initial measure µ (x,y) and transition matrix L.

Properties of the Synchronization Product

Let ω ∈ Ω be a history on n sites. Then ω naturally decomposes as ω = (w, z) with w a history on the n -1 first sites and z a sample path on the n th site, whence a natural injection Ω → Λ. This injection allows to transpose on Ω whatever mappings previously defined on Λ, such as σ : Λ → N . Corollary 4.1 Let u be a given finite trajectory on n sites. Consider the sequence(S k ) k≥0 of iterated stopping times associated with S previously defined in Lemma 4.2. The mapping V : Ω → N defined as follows:

∀ω ∈ Ω, ω V = {ω S k : ω S k ≥ u, k ≥ 0}, if ∃k ≥ 0, ω S k ≥ u, ∞, otherwise,
is a stopping time.

Proof. Consider the stopping time U : Ω → N defined by:

∀ω ∈ Ω, ω U = u, if u ≤ ω, ∞, otherwise. (25)
By Lemma 4.2, S : Ω → N is a stopping time. Applying Lemma 3.1 to stopping time S and its iterated sequence (S k ) k≥0 on the one hand, and to stopping time U defined by Eq. ( 25) on the other hand, yields the result of the corollary.

Proposition 4.2 The synchronization product P described in Def. 4.3 is a DMP over n sites.

Proof. We denote by x a generic global state of the n -1 first sites, that is, an element of X [1,n-1] , and by y an element of X n . What needs to be proved is that the quantity

P (x,y) ↑ (t • t ′ )| ↑ t (26)
only depends on γ(t) and t ′ , for t, t ′ ranging over the finite trajectories on n sites and for (x, y) ranging over X = X [1,n-1] × X n . For this, we proceed by steps. In the remaining of the proof, and accordingly with the notations of ( 26), (x, y) denotes the initial global state of the system distributed over n sites, and t denotes some finite trajectory over n sites. We also put (x ′ , y ′ ) = γ(t).

Step 1. Without loss of generality, one may assume that t ′ has the form t ′ = ω S , where S is the stopping time defined in Lemma 4.2. Indeed, assume that the result is true for any t ′ of this form, then by the chain rule it is clearly true also for t ′ any finite concatenation of such trajectories. Hence the result is true for t ′ of the form t ′ = ω S k for any k ≥ 0, with (S k ) k≥0 the sequence of iterated stopping times associated with S. Let now t, t ′ be any finite trajectories over n sites.

By Corollary 4.1, we have that V : Ω → N defined by

∀ω ∈ Ω, ω V = {ω S k : ω S k ≥ t ′ , k ≥ 0}, if ∃k ≥ 0, ω S k ≥ t • t ′ , ∞, otherwise,
is a stopping time. We decompose over the possible values v of V to get:

P (x,y) ↑ (t • t ′ ) ↑ t = v P (x,y) ↑ (t • t ′ ) ∩ {t • ω V = t • v} ↑ t = v P (x,y) ↑ (t • v)| ↑ t since V is a stopping time.

Now by assumption of

Step 1 each term of the sum only depends on γ(t) and v. And by construction, the range of v only depends on t ′ . It follows that the sum itself only depends on γ(t) and t ′ .

Step 2. Without loss of generality, one may assume that t decomposes as t = w • z with w over the n -1 first sites, z over X n , and w / / z. In other words, one may assume that θ n (t) ∩ X i = ∅ for i = 1, . . . , n -1. For, if t is any finite trajectory, let z = θ n (t). Let k be the number of synchronizations of z with the n -1 first sites. Then if (w 1 z 1 ), . . . , (w k , z k ) denote the k successive synchronization times of t, then t can be written as follows:

Inductive Construction of Synchronization Products

Lemma 4.3 Let Q and M be two DMP on n -1 and 1 site respectively, with global states X [1,n-1] = X 1 × • • • × X n-1 and X n respectively. We assume that the transition matrix M associated with Markov chain M has all its coefficients positive, and that Assumptions 1-3 from § 4.2 are in force.

Then the synchronization product of Q and M, with global states in X = X [n-1] × X n satisfies the following property: for every global state x ∈ X and for any site i = 1, . . . , n, any local state s ∈ X i is locally reachable from x (see Def. 2.2).

Proof. We keep the notations of § 4.3- § 4.4. Let (x, y) ∈ X [1,n-1] × X n be any initial global state on n sites, and let s be some local state of some site i. We show the existence of some finite trajectory u ∈ S on n sites such that P (x,y) > 0 and θ i • γ(u) = s.

1. Assume that i = n.

(a) If s ∈ X j for all j = 1, . . . , n -1. Put u = (ǫ, . . . , ǫ, s). Consider any value (w, z) of λ σ defined in § 4.1. Then u • (w, z) is also a value for λ σ . Therefore:

P (x,y) (↑ u) ≥ P (x,y) (  u • (w, z) = K (x,y) λ σ = u • (w, z) > 0.
(b) Is s ∈ X j for some j = 1, . . . , n -1. By Assumption 2, s is locally reachable from x. Hence there is a finite trajectory v ∈ S [1,n-1] on n -1 sites such that θ j • γ(v) = s, and Q x (↑ v) > 0. Next extract from v all states that belong to X n , and sequentialize them into a path on X n with s as its last element (such a path exists).Then, by construction, v and z are compatible and by putting u = v z we have θ [1,n-1] (u) = v and θ n (u) = z. Furthermore, if k = |z|, and if (S k ) k≥0 denotes as before the sequence of iterated stopping times associated with stopping time S, then any ω ∈ Ω satisfying ω ≥ u satisfies S k (ω) = u. By construction of the synchronization product, and since z has positive M y -probability since by assumption the transition matrix only has positive coefficients, it follows that P (x,y) (S k = u) > 0. Since θ n (u) = z, we also have that γ • θ n (u) = s since s is the last element of z.

2. Assume that i < n. The construction of the element u is similar to the case 1b above.

Theorem 4.1 For any n ≥ 1 and for finite space states X 1 , . . . , X n of arbitrary finite size, there exists DMP inductively obtained by synchronization product with global states X 1 × • • • × X n .

  Lemma 4.2 The mapping S : Ω → N obtained by composition of the natural injection Ω → Λ with σ : Λ → N defined in Def. 4.2, is a stopping time in the sense of Def. 3.1.

t

  = (w 1 z 1 ) • . . . • (w k z k ) • r ,with r of the following form:r = u • s ,

with u over the n -1 first sites, s over X n and u / / s. Finally, putting

we have:

Since t 0 is a concatenation of finite trajectories of the form w i z i , and by the construction of the synchronization product, it follows by the chain rule that

Reinserting the right members of Eqs. ( 28) and ( 29) in ( 27), we get:

Now, by the assumption made in this step, and since r has the requested form, we find that the right member of Eq. ( 30) only depends on γ(r) = γ(t) and t ′ , what was to be shown.

Step 3. The basic case. According to the reduction established in Steps 1-2, let t be a finite trajectory over n sites satisfying θ n (t) ∩ X i = ∅ for i = 1, . . . , n -1, and let t ′ = w z. Then, according to the construction of the synchronization product P, we have:

On the last line, since both Q x and M y are Markovian, it is clear that the expression only depends on γ(t) and t ′ .

Proof. For each i = 1, . . . , n, consider any transition matrix with positive coefficients on X i . Let M i denote the probabilistic process associated on X i associated with the given transition matrix. Rename the elements of X i so that

Then proceed inductively to the construction of the DMP P i on X [1,i] = X 1 × • • • × X i , synchronization product of P i-1 with M i . By Lemma 4.3, P i has the property that any local states are locally reachable. Hence Assumption 2 is fulfilled. Assumption 1 is also fulfilled since the transition matrices all have positive coefficients. Finally, Assumption 3 is also fulfilled thanks to Eq. ( 31). Hence the construction of P i+1 is allowed, completing the proof.