O Laligant

F Truchetet

A Dupasquier

Edge enhancement by local deconvolution

Keywords: edge detection, restoration, deconvolution, multiscale merging

In this paper a new approach for blurred image restoration is presented. Our algorithm is based on human vision which zooms back and forth in the image in order to identify global structures or details. Deconvolution parameters are estimated by an edge detection and correspond to the ones of a chosen edge detection model.

The segmentation is obtained by merging multiscale information provided by multiscale edge detection. The edge detection is achieved by using a derivative approach following a generalization of Canny-Deriche filtering. This multiscale analysis performs an efficient edge detection in noisy blurred images. The merging leads to the best local representation of edge information across scales. The algorithm deals with a mixed (coarse-to-fine/fine-to-coarse) approach and searches for candidate edge points through the scales. Edge characteristics are estimated by the merging algorithm for the chosen model. Scale, direction and amplitude informations allow a local deconvolution of the original image. The noise problem is not considered in this work since it does not disturb the process. Results show that this method allows non-uniformly blurred image restoration. An implementation of the whole algorithm in a intelligent camera (DSP) has been performed.

Introduction

Shooting a real world image with a camera through an optical device gives a 2D image where at least some parts are affected by a blur. This optical blurring can generally be correctly modeled by the effect of a Gaussian filtering on the image. Furthermore, the blur of real edges can have other sources than optical ones: movement of the item in the scene, coarsiness of resolution in the imaging process, effect of diffuse lighting, etc. Restoration of blurred images is one of the main topics in image processing and many methods have been described in literature. Generally these methods are based on a filtering process, the major difficulty in image restoration being the lack of knowledge about the degradation process [START_REF] You | A regularisation approach to joint blur identification and image restoration[END_REF]. In some methods the degradation is assumed to be well known and the restoration problem is reduced to the inversion of the degradation process [START_REF] Banham | Digital image restoration[END_REF][START_REF] Costello | One-dimensional comparison of wiener filtering and richardon-lucy methods for sectioned restoration of sapce-variant digital images[END_REF]. In some other approaches the degradation is supposed to be uniform and the method tends to estimate both the original image and the blur from the degraded image characteristics [START_REF] Kundur | Blind image deconvolution[END_REF][START_REF] Banham | Spatially adaptive wavelet-based multiscale image restoration[END_REF][START_REF] Guo | Blind restoration of images degraded by space-variant blurs using iterative algorithms for both blur identifi-cation and images restoration[END_REF]. Our present work lies at the frontier between these two approaches. We assume that the blur degradation process is according to a new parametric model. This model is an edge profile characterized by three parameters: scale, amplitude and orientation. The principle of the proposed restoration method (IBEP: identifying blur edge profile) consists in locally identifying the edge profile parameters Email address: email: o.laligant@iutlecreusot.u-bourgogne.fr (O.

Laligant * ,).

and then enhancing the edge by a local deconvolution. Identification of various profiles is achieved by a multiscale edge detection assisted by a fusion process to eliminate the problems of over detection and mutual influence (delocalization). This stage produces an image of edge points associated to profile parameters. These data are then used for local and stable deconvolution. The figure 2 presents the overview of the method. Note that the method operates independantly from the noise and does not affect the signal to noise ratio. This non iterative method has been developed and implemented to achieve a real-time autonomous system.

In the first section, we present the blur profile model as well as the operator used for edge detection. The second section is dedicated to the scheme of multiscale edge detection including the phase of fusion. This stage deals with the mutual influence problem and produces information used to estimate the local model. The enhancement method illustrated with practical results and a comparison with the well known and classical Wiener filtering approach are presented in the last section.

Edge detection

Edge model

The profile of an edge affected by blurring is modelled by a low-pass exponential function. The normalized edge profile model (figure 1) is defined by equation ( 1):

C sc (u) =                1 -e -sc u 2 u ≥ 0 e sc u 2 u < 0 (1)
This function can be seen as a first-order approximation of the majority of edges encountered in real image processing and it corresponds to the filtering of an ideal step edge by a low-pass filter with impulse response: h(u) = sc 2 e -sc|u| . A somewhat more realistic model should be a step function filtered by a gaussian filter (erf function) but its form is far less suitable for mathematic treatment than the proposed approximation. The parameter s c represents the scale of the blur. The smaller the parameter value is, the greater the blur is. The amplitude of the edge will be denoted by A c and its orientation with respect to image lines direction by θ c (u is the abscissa in the direction θ c ). Hence, an edge profile EP c in the image is characterized by three parameters and will be denoted as:

EP c = (s c , A c , θ c ) (2)
The estimation of these three parameters is achieved on the basis of a multiscale edge detection method.

Edge detector

The restoration process is always more or less equivalent to high pass filtering.

It is well known that such an operation is generally unstable and very sensitive to high frequency corrupting signal (noise). In our design, in order to avoid this drawback, the restoration process is applied only where a true edge is present. Therefore these edges have to be detected and localized, even if the image is noisy. This task is realized by optimal edge detection such as the one proposed by Canny and Deriche or Shen and Castan. These filters are designed to give a optimum response with respect to constraints about signal over noise ratio, edge localization, multiple responses suppression. We have shown in previous works that these filters can be seen as special cases of a familly of filters where response is optimal for the edge model proposed in the previous section. These filters are parametrized by a scale factor and they lead to a multiscale gradient estimator. The prototype impulse response of the generalized filter family [START_REF] Laligant | Generalization of shen-castan and cannyderiche filters[END_REF] is defined as following:

f n,j (x) = -sign(x) n+1 s n+1 j n! x n e -s j |x| ( 3
)
where n is the index of the filter and j the scale level. We choose the scale factor s j as:

s j = s a j
with s, a > 0 and j ∈ Z.

Shen & Castan [START_REF] Shen | An optimal linear operator for edge detection[END_REF] and Canny-Deriche [START_REF] Deriche | Using Canny's criteria to derive a recursively implemented optimal edge detector[END_REF][START_REF] Bourennane | Generalization of canny-deriche filter for detection of noisy exponential edge[END_REF] filters correspond respectively to n = 0 and n = 1. Each filter is normalized with respect to unit step response:

0 -∞ f n,j (τ )dτ = +∞ 0 f n,j (τ )dτ = 1 (4)

Responses and performances of the detectors

The filter is to be applied on the image in order to detect, localize and identify (determine the parameters) the edges.

It is therefore necessary to analyze its response and performances with respect to the blur scale s c .

Response

The response maximum y(0) of these detectors to the edge model is (see details in appendix A):

y n,j (0) = f n,j * C sc (0) (5) = 1 - s n+1 j (s j + s c ) n+1 (6)
This result will be used to design the enhancement process (section 4). This maximum decreases as the scale factor s j increases for s c fixed. It increases as the blur parameter s c increases for s j fixed.

Edge localization

The definition given by Canny [START_REF] Canny | A computational approach to edge detection[END_REF] of a criterion for localization of an edge C in presence of additive white Gaussian noise of variance unity and detected by an operator f is as following:

L(f ) = +∞ -∞ C sc (-τ )f (τ )dτ +∞ -∞ f 2 (τ )dτ 1 2 (7)
therefore (see details in appendix B), for n = 1, 2, 3, ...

L(f n,j ) n∈{1,2,...},j∈Z = s 2 c (s j + s c ) n+1 2n -1 s j (2n)! (2s j ) n n! (8)
Taking into account the 0-order discontinuity at x = 0 for n = 0

L(f 0,j ) j∈Z = s c s j + s c √ s j (9)
For each detector, the localization criterion presents a maximum according to the detection scale s j . We define the optimal scale S n S n = arg max

s j (L n (s j )) n∈{1,2,...} (10)
Value S n of s j corresponding to the maximum in localization depends on the blur scale s c and naturally on the index n of the used detector. On both sides of this extremum, localization criterion tends toward 0 when s j → 0 and when s j → ∞. The scale of the detector for this maximum in localization is given for n ∈ {1, 2, 3, ...} by (see appendix D):

S n = 1 3 (2n -1)s c n∈{1,2,...} (11)
and for n = 0:

S 0 = s c (12)

Signal to noise criterion

The definition given by Canny for signal to noise ratio in detection of an edge C in the presence of additive white Gaussian noise of variance unity and filtered by an operator f is as following:

Σ(f ) = +∞ -∞ C sc (-τ )f (τ )dτ +∞ -∞ f 2 (τ )dτ 1 2 (13)
therefore (see details in appendix C)

Σ(f n,j ) n∈{0,1,2,...},j∈Z = (s j + s c ) n+1 -s n+1 j (s j + s c ) n+1 2 n n! s j (2n)! (14)
For an edge scale s c , this criterion is maximum for s j = 0 and decreases when s j increases. The evolution of this criterion is naturally connected to the maximal response presented above.

Strategy for restoration

The principle of restoration method IBEP (identifying blur edge profile) is based on the identification of the model for edge profiles in the image. This identification consists in estimating EP c . The purpose of the restoration is to enhance the edges by local deconvolution to obtain a step profile:

EP ∞ c = (A c , θ c ) (15)
A detection on two scales at least is necessary to determine the scale s c and the amplitude A c of EP c , orientation being deduced from the gradient direction information. A study of the proposed detectors shows on one hand that the best localization of edge is obtained for a detection scale given by equations 12, 11. On the other hand, the immunity to noise increases as the scale of detection decreases. It is the same for the maximum response. However, for too small scales, localization falls very quickly.

As blur is assumed to be variable in the image, we have to choose for detection an interval of scales s close to the range of blur scales s c .

From this analysis we can conclude that for detection, localization and characterization of the edges a strategy has to be followed. This strategy in a multiresolution approach consists in choosing a set of detection scales with a geometric regular step (noted a) compatible with expected blur in the image, then from the detection scales close to the maximum of localization, in determining exactly the scale, the direction and the amplitude of the best local model.

Edge image representation

Multiscale edge detection

For each point I(x, y) of the image, the gradient vector is obtained by a 2Dseparable convolution. Following J. Canny [START_REF] Canny | A computational approach to edge detection[END_REF], let h n,j be the regularization operator associated with f n,j :

dh n,j (x) dx = f n,j (x). A 2D separable operation is defined as: h n,j (x, y) = h n,j (x).h n,j (y). Therefore the gradient vector ----→ G(x, y)
is defined as:

----→ G(x, y) =         G x (x, y) G y (x, y)         =         d dx (h n,j (x, y) * I(x, y)) d dy (h n,j (x, y) * I(x, y))         n∈{0,1,...},j∈[0,1,...,p-1] (16) =         f n,j (x)h n,j (y) * I(x, y) h n,j (x)f n,j (y) * I(x, y)         n∈{0,1,...},j∈[0,1,...,p-1]
It is easy to compute the edge orientation:

θ = Arg[ ----→ G(x, y)]. Multiscale anal-
ysis is achieved by a multiscale set of p operators {f n,j } j∈[0,1,...,p-1] of scales

s j = s a j with j ∈ [0, 1, ..., p -1]
, s being initial scale and a > 0. In this way, a set of p gradient images is obtained:

{ - → G j } j∈[0,1,...,p-1]
. Fine scales and coarse scales correspond respectively to small and large values of s j .

Candidate edge points {L j (x, y)} j∈[0,1,...,p-1] in the gradient images are localized by extraction of local maxima. After detection and localization, the next step (see 2) consists in estimating s c and A c in the original image. However, one observes variations of localization and amplitude due not only to noise but also to influence of neighbor edges. This phenomenon has not been really analyzed in the literature dedicated to edge detection. In the next section, we propose a simple and pratical study leading to some rules allowing these perturbations to be avoided or taken into account in the restoration process.

Influence of a neighbor edge

It is the distortion induced by a neighbor edge on the detection of the current edge that is under study. To simplify the discussion of this influence we con-sider only infinite scale edge profiles. However, the observations will remain valid whatever may be the scales. Two kinds of influence can be distinguished depending on the respective configuration of the two edges. Opposite edges lead to a substractive influence while two edges in a staircase-like configuration give an additive influence.

• Firstly we have a configuration of two parallel opposite edges (figure 3), cur-

rent profile EP ∞ c = (∞, A c > 0, θ c ) and neighbor profile EP ∞ v = (∞, A v < 0, θ v = θ c ). h n is the step response to the detector f n . The response to EP ∞ c
is the superposition of two responses (to EP ∞ c alone and EP ∞ v alone):

D n,j (x) = A c h n,j (x) + A v h n,j (d -x) (17)
d being the distance between the two edges and

h n,j (x) = e -s j |x| n i=0 s n-i j (n -i)! |x| n-i (18)
An edge is assumed to be detected where response is maximum. Therefore an edge of amplitude A c is detected at location x c close to 0 if there exists a zero-crossing in the range ]x c ; x v [ where x v is the location of the detection of a neighbor profile; x c is negative or null (x c 0). Absolute value of

x c increases and the amplitude of the response maximum decreases as the influence of the neighbor edge increases.

Proposition 1 Subtractive influence is the influence of a neighborhood of gradient sign opposite to the current point one. It tends to decrease the detection amplitude as the detector scale decreases and it displaces the highest response maximum toward the fine scales.

• Secondly, let us consider the other configuration where edges profiles form a stair (figure 4):

EP ∞ c = (∞, A c > 0, θ c ) and EP ∞ v = (∞, A v > 0, θ v = θ c ).
In this case, the profile EP ∞ c of amplitude A c will be detected if there exists a minimum between two maximum responses to two profiles EP ∞ c and EP ∞ v (so that they can be separated).

The influence of the neighbor edge tends to displace the response to current edge: in that case x c is positive. The amplitude of the response maximum increases as the influence of the neighbor edge increases (detection scale decreases). This remark leads to a second proposition:

Proposition 2 Additive or amplification influence is the influence due to a neighborhood whose gradient sign is the same as the current point one. It tends to increase the detection amplitude as the detector scale decreases and it displaces the highest response maximum toward the coarse scales.

Remark 1 Due to its infinite derivative at 0, Shen & Castan filter does not delocalize the responses whatever the influence. Nevertheless, the amplitude increases or decreases following the above propositions.

Merging

Merging, as considered in IBEP (see figure 2) is defined as the operation consisting in the choosing of the scale providing the best local edge representation (at a given pixel) among all available scales [START_REF] Dupasquier | Restoration of blurred images by edge detection[END_REF].

We have seen, in a previous section, how in a first step, a set of gradient{

- → G j } j∈[0,...,p-1]
and local maximum {L j } j∈[0,...,p-1] images are obtained from the multiscale analysis method. In the next step one has to localize and to characterize local edge profiles (EP c ) by merging these data.

Gradient orientations are coded into eight sectors of 45 degrees according to the Freeman code (see figure 8):

{θ 0 = {[-22.5 • , 22.5 • [}, θ 1 = {[22.5 • , 67.5 • [}, θ 2 = {[67.5 • , 112.5 • [}, ...θ 7 = {[292.5 • , 337.5 • [}}.
According to the discussion on the influence of neighbor edges, it is obvious that subtractive or additive influence will delocalize the current detected edge.

This delocalization has to be taken into account. But the initial scale a is assumed to be small enough to avoid any delocalization on the first scale.

The merging algorithm joins the coarse-to-fine approach and the fine-to-coarse one. Coarse-to-fine approach allows to follow the evolution of the tight edge profiles across the scales. The best localization of these edges is performed by this approach. Conversely, fine-to-coarse approach leads to the best detection for wide edge profiles. A good edge characterization of the various profiles occuring in a real world image needs the cooperation of these two approaches.

Local maxima of edge gradient are searched in a neighborhood cone through the scales (figures 6,7). During this search, local maximums under influence of the neighborhood are picked out and eliminated. The greatest amplitudes give the best localization (see figure 5).

The merging algorithm is summarized in the following lines.

• for each current point of the original image, from multiscale gradient information:

• search for the closest maximums in the coarse-to-fine approach. Ordered maximums lead to a first list of candidates points.

• Search for the closest maximums in the fine-to-coarse approach. Ordered maximums lead to a second list of candidates points.

• Discard points under positive influence (see below), points under soustrac-tive influence eliminate themselves.

• Through the two ordered lists: a point is valid if each list represents a maximum at the same scale, same direction and same localization.

A valid point produces an edge if this localization is exactly the position of the current point.

Elimination of points under positive influence is carried out as follows: maximum are sorted in decreasing scales (fine-to-coarse). Two consecutive maximums which do not present the same localization are eliminated.

Following the merging stage, a segmentation method is needed to obtain a binary edge points image. Classic methods like a simple thresholding or the coïncidence assumption of Marr and Hildreth [START_REF] Witkin | Scale-space filtering[END_REF] [12] can be used. This last one selects points according to the persistence of the detection across the scales and consequently it is well adapted to our merging method.

Finally, we obtain a list P c of ordered points at the current position (x c , y c ):

P c = {(Gx j A , Gy j A , s j A ), (Gx j B , Gy j B , s j B ), ...}). These points are the best edge information available and represent gradient information.

Restoration

The IBEP restoration principle is based on a local evaluation in image of scale and amplitude equivalent to the model. The deconvolution based on the model allows to restore the contour. Deconvolution parameters are those of the model EP c , they can be deduced from the multiscale edge analysis described in the previous section.

Model parameters estimate

This part is developped on the basis of the Shen and Castan [START_REF] Shen | An optimal linear operator for edge detection[END_REF] operator regarding its efficency and simplicity.

Following the notation introduced in the equation 3, the definition of Shen-Castan operator is:

f 0,j (x) = -sign(x)s j e -s j |x| (19)
characterized by:

                                       y 0,j (0) = 1 - s j s j +sc L(f 0,j ) = sc s j +sc √ s j s 0 [max(L)] = s c Σ(f 0,j ) = sc s j +sc 1 √ s j (20)
The above relations are defined in the continuous domain, to improve the precision of the results, discrete versions can be introduced. The normalized filter ( +∞ 0 f 0,j [k] = 1) becomes:

f 0,j [k] = -N s j sign[k]e -s j |k| (21)
with N s j = e s j -1. The maximum response y 0,j [0] to a simple edge profile is:

y 0,j [0] = ∞ k=-∞ C sc [l -k]f 0,j [k] (22) = A c 1 - N s j N s j +sc
with N s j +sc =e s j +sc -1.

Four detection scales are processed with a geometrical progression step a = 2.

A minimum of two scales of detection is necessary to estimate the parameters of the current edge profile EP c . As presented in section 3.3 edge points are selected by a segmentation method which is the coincidence assumption [START_REF] Witkin | Scale-space filtering[END_REF] assisted by a simple threshold. The merging stage produces the list P c with all needful information leading to the profile estimation. θ c is directly deduced from the gradient information of the first point in the list P c . The resolution of a system of equations gives a solution for A c and s c .

Deconvolution

The edge profile is a monodimensional object characterized by its amplitude A c , its equivalent scale s c and its orientation θ c ∈ {θ 0 , θ 1 , ..., θ 7 }. In the direction θ c , the model profile can be seen as the convolution of a step edge with the following Point Spread Function:

T (u) = s c e -sc|u| 2 ( 23
)
where u is the abscissa in direction θ c . Let I 0 be the original image not affected by blur and I be the degraded image. In direction θ c :

I(u) = I 0 * T (u) (24)
By direct inversion of the blurring operator one obtains a deconvolution function:

T -1 (u) = δ(u) - 1 s 2 c δ (u) (25)
The enhanced image could be deduced as follows:

I 0 (u) = I(u) - 1 s 2 c I (u) (26)
However, as it is well known, the straightforward application of this equation leads to an unstable result whatever the choice of the second derivative operator. We solved this problem by simply substituting the second derivative of the edge model (EP c ) for the second derivative of the image. Along an edge profile EP c = (s c , A c , θ c ), the enhanced image is given by:

I 0 (u) = I(u) - 1 s 2 c A c C sc (u) (27)
In practice, the deconvolution is separable and it can be applied along the two directions. Depending on its orientation, an edge is enhanced as following:

• θ c = θ 0 or θ 4 : I 0 (x) = I(x) -1 s 2 c A c C sc (±x) (C sc (+x) for θ 0 , C sc (-x) for θ 4 ).
• θ c = θ 2 or θ 6 : I 0 (y) = I(y) -1

s 2 c A c C sc (±y) • θ c = θ 2k+1 k∈[0,3] (odd index) : I 0 (x, y) =                I(x) -1 s 2 x A c C sx (±x) I(y) -1 s 2 y A c C sy (±y) with s x = s y = 1 √ 2 s c .
The deconvolution has to be a local process, we choose to limit it so that 95 % of the amplitude of the theoretical profile is corrected. This limit corresponds to a spatial range around the current point of [x min , x max ] with -s c x min = s c x max = 3, x = 0 being the position of current edge point in the image I. Two patterns with two different blurring scales for edge profiles have been designed in figure 13 (original image in figure 12). The enhanced image in figure 14 confirms the ability of the IBEP method to restore different blurs if the equivalent scales are included in the analysis range. Furthermore, this edge enhancement method does not increase the noise which is simply let unchanged.

Results

Blanks of restoration are the results of no edge detection or wrong estimation in the fusion stage. The effects inherent to edge detection by the derivative approach, as the rounding of corners (see figure 14), affect the restoration process.

To measure the efficiency of the method, results have been compared with those given by the standard Wiener restoration method. In this method there is always a trade off between restoration and noise or blur. The synthetic image (figure 12) has been locally convolved by a P SF function with two different scale parameters s 1 et s 2 to obtain two edge profiles in the image (figure 13).

Then gaussian noise of variance v is added. The original image (sharp and clean) is used as reference for quality assessment of the restoration methods.

Table 1 presents the results with respect to the degradation parameters for the two methods. Note that only the best result obtained with Wiener restoration has been consigned in the table.

This table shows that in absence of noise, the presented enhancement method gives better results than Wiener filter, this last being not adaptive. As the variance of noise increases, the performances of the enhancement method measured by mean square error decrease. IBEP method does not process in any way the noise and MSE table takes into account the degradation due to the noise as well as the one due to blurring. However, these results prove that the method is stable for images corrupted by noise.

Conclusion

We have presented a new edge enhancement adaptive method based on the modelization of the edge profiles. Local estimation of the edge model is achieved by a multiscale edge detection process. Edge detection process includes mul-tiscale edge detection, multiscale merging and segmentation, this last stage being performed by coincidence assumption and simple thresholding.

This very simple method is efficient especially when the degradation process is not space-invariant and gives better results than Wiener restoration whatever its adjustment. The edge enhancement can be used in cooperation with a denoising method to build a complete restoration system. Our method does not need any parameter adjustment: scales of detection and number of levels have only to be chosen according to scale range in the image to restore.

IBEP has been implemented in an autonomous system composed by a camera and a DSP. The figure 15 shows an example of obtained result. The rate of the system is 0.5 image per second.

The mutual influence of the edges has been taken into account to minimize estimation errors of the model; lack of restoration are mainly due to these errors. Some improvements are expected from another study still in progress in which we try to interpolate the correct values of the neighborhood. An other way for improving could be based for example on ridge detector. The IBEP method does not take into account the noise and these results show that it is stable against this noise.

M SE v = 0 v = 2 v =

A Response of the detectors to the edge model

The response maximum y(0) of the detectors (eq. 3) to the edge model (eq. 1)

is :

y n,j (0) = f n,j * C sc (0) (A.1) = +∞ -∞ f n,j (τ )C sc (-τ )dτ = 0 -∞ f n,j (τ )dτ - 0 -∞ f n,j (τ ) e scτ 2 dτ + ∞ 0 f n,j (τ ) e -scτ 2 dτ
Symmetry considerations of the filters give:

y n,j (0) = 0 -∞ f n,j (τ )dτ - 0 -∞ f n,j (τ )e scτ dτ
For x < 0, the filters are defined by:

f n,j (x) = -sign(x) n+1 s n+1 j n! x n e s j x (A.2)
Hence the second term can be interpreted as the integral of a new filter f n,j+c

with a new scale s j + s c and can be written:

y n,j (0) = 0 -∞ f n,j (τ )dτ - s n+1 j (s j + s c ) n+1 0 -∞ f n,j+c (τ )dτ
Finally, introducing the regularization (integral) operator, we obtain:

y n,j (0) = h(0) - s n+1 j (s j + s c ) n+1 h n,j+c (0) = 1 - s n+1 j (s j + s c ) n+1 (A.3) B Localization L(f ) = +∞ -∞ C sc (-τ )f (τ )dτ +∞ -∞ f 2 (τ )dτ 1 2
We write: This scale is denoted S n . The demonstration of S 0 is trivial.

f n,j (

Figure 9

 9 Figure 9 shows an example of enhancement on an image section. The sharpness of the edge profiles is clearly improved. Figures 10 and 11 present a restored real image for which there exists no original image. This example shows that strip borders are not restored. Indeed, strip borders present a configuration of two opposite and close edge profiles. This case of subtractive influence whatever the detection scales does not allow to estimate the profile parameters.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Edge model defined as an exponential function with s c > 0

s 1

 1 = 0.2, s 2 = 0.4

I 0 2 .

 02 (x, y) -I 0 (x, y) Degradation parameters: v noise variance, s 1 , s 2 blur scales. The best results obtained with Wiener restoration have been consigned.

Fig. 4 .

 4 Fig. 4. Additive influence of a neighbor edge: the detection amplitude increases and the highest response maximum moves toward the coarse scales as the detector scale decreases.

Fig. 5 .

 5 Fig. 5. (a) Edge profiles (s c = 3 and s c = 0.7) and (b), (c), (d) multiscale analysis (s 0 = 1.6, s 1 = 0.8, s 2 = 0.4). The best detection is achieved by the first scale for "s c = 3" profile while the second scale performs the best detection for "s c = 0.7" profile.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Fine to coarse approach.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Gradient orientations according to the Freeman code θ i .

Fig. 10 .

 10 Fig. 10. Original real degraded image.

Fig. 11 .

 11 Fig. 11. Restoration of a real degraded image. Lines not restored are under subtractive influence whatever the detection scales.

Fig. 12 .

 12 Fig. 12. "Square image": original image.

Fig. 13 .

 13 Fig. 13. "Square image": blurred and noisy test synthetic image.

Fig. 14 .

 14 Fig. 14. Restored image by the edge enhancement method.

Fig. 15 .

 15 Fig. 15. Restoration by an autonomous system: camera (740 × 580 pixels) and DSP (ADSP 2181 -40MHz). This figure shows the blurred image and a piece (256 × 256 pixels) of restoration (inside the square outline).

2 ) 4 ) 2 c 2 j 0 n 2 2 j 1 = a 2 n 2 ( 2 c 8 ) 2 . 2 ) 3 = d ds j s 2 c= n!2 n s 2 c

 24220221222822322 x) x 0 = -a n x n e -s j |x| (B.1) with a n = sign(x) n+1 s n+1 jWe have to take into account the 0-order discontinuity (at the abscissa x = 0) for n = 0. Then for n = 0: j τ 0 e -s j τ dτ= s 2 j s c s j + s c (B.3)And more generally for n ∈ {1, 2, ...}: nτ n-1 + s j τ n e -s j τ dτ= a n s c n +∞ 0 τ n-1 e -(s j +sc)τ dτa n s c s j +∞ 0 τ n e -(s j +sc)τ dτ = a n n! s c (s j + s c ) na n n!s j s c (s j + s c ) n+1 = a n n! s (s j + s c ) n+1 (B.5) e -2s j τ s j τ nnτ n-1 2 dτ τ 2n e -2s j τ dτ -+∞ 0 2s j nτ 2n-1 e -2s j τ dτ + +∞ 2τ 2n-2 e -2s j τ dτ = 2a 2 n (2n)! s (2s j ) 2n+1 -(2n -1)! 2s j n (2s j ) 2n + (2n -2)! n 2 (2s j ) 2n-2s j ) 2n-1 (2n)! -4n(2n -1)! + 4n 2 (2n -0,j ) j∈Z = s c s j + s c √ s j (B.7) L(f n,j ) n∈{1,2,...},j∈Z = s (s j + s c ) n+1 2n -1 s j (2n)! (2s j ) n n! (B.Remark Replacing n! (2n!) by the Γ n+1 (Γ 2n+1 ) function extends the relation B.8 to real n.C Signal over noise ratioΣ(f ) = +∞ -∞ C sc (-τ )f (τ )dτ +∞ -∞ f 2 (τ )dτ the result A.3 of A.1 we obtain: Σ(f n,j ) n∈{0,1,2,...},j∈Z = (s j + s c ) n+1s n+1 j (s j + s c ) n+1 2 n n! s j (2n)! (CRemark Replacing n! (2n!) by the Γ n+1 (Γ 2n+1 ) function extends the relation C.2 to real n.D Scale of the detectors for maximum in localizationFor a given edge scale s c , we search for the scale s j leading to the maximum in localization. This particular scale is denoted S n .The maximum in localization corresponds to :dL(f n,j ) ds j s j =Sn = 0 Therefore for n ∈ {1, 2, ...} dL(f n,j ) ds j (s j + s c ) n+1 2n -1 s j (2n)! (2sj ) n n! + s c ) n+1 -(n + 1) (s j + s c ) n s j + s c ) -(n + 1) = 0 (D.2) Finally ∀s c > 0, n ∈ {1, 2, ...} dL(f n,j ) ds j = 0 ⇒ s j = 1 3 (2n -1) s c (D.3)

  PhD in electronics at the same university in 1977. He was with Thomson-CSF for two years as a research engineer and he is currently full professor and the head of the image processing group at Le2i, Laboratory of Electronic, Computer, and Imaging Sciences, at Université de Bourgogne, France. His research interests are focused on image processing for artificial vision inspection and particularly on wavelet transforms, multiresolution edge detection, and image compression. He is member SPIE, IEEE, and ISIS (a research group of CNRS).

	Olivier Laligant received his PhD degree in 1995 from the Université de Bour-gogne, France. He is an assistant professor in the computing, electronic, imaging department (Le2i) at Université de Bourgogne, France. His research interests are focused on multiscale edge detection, merging of data, wavelet transform, and re-cently motion estimation. Frédéric Truchetet received the master degree in physics at Dijon University, France, in 1973 and a Alain Dupasquier is engineer from CNAM (Conservatoire National des Arts et Métiers) in an industrial group of manufacture of plastic toys.

Table 1

 1 Performances of IBEP and Wiener restoration on the "square image". M SE =

	4


Antoine Ayache
email: ayache@cict.fr

Jacques L Evy V Ehel

23) Probabilit es 8) Analyse math ematique Processus a r egularit e locale prescrite

R esum e

La fonction de H older associ ee a une fonction permet de mesurer les variations locales de son irr egularit e. Dans 4], il est montr e que H est la fonction de H older d'une fonction continue f si et seulement si H s' ecrit comme une limite inf erieure de fonctions continues. Une construction possible de f repose sur une g en eralisation de la c el ebre fonction de Weierstrass. Dans cette note nous transposons ce r esultat au cadre stochastique.

Processes with prescribed local regularity

Abridged english version

Let X = fX(t)g t2 0;1] be a continuous and nowhere di erentiable stochastic process. The H older process of X is de ned for all t as : X (t) = sup lim sup h!0 jX(t + h) X(t)j jhj = 0 :

This process allows to measure the local variations of the irregularity of X. It is shown in 4] in the deterministic context, that a function H is the H older function of a continuous function f if and only if H is of the form H(t) = lim inf H l (t), for all t, where the H l are continuous functions (In fact, it is not restrictive to suppose that the H l are H older functions whose H older constants increase slowly when l goes to in nity i.e c l = O(l)). The aim of this note is to construct a continuous random process W extending the well-known Weierstrass function and whose H older process may be, with probability 1, any lower limit of continuous functions with values in a compact of ]0; 1 . More precisely, let IL = (l p ) p2I N be the sequence de ned as l 1 = 1 and for all p 2, l p = 1 a 0 1 b 0 l p 1 ] + 1 where a, b, a 0 , b 0 are xed reals satisfying 1 > b 0 > b > a > a 0 > 0. Let (H l ) l2I L denote a sequence of ( ; c l )-H older functions with values in a; b] and such that c l = O(l) (we say that f is a ( ; c)-H older function if for all t 1 , t 2 we have jf(t 1 ) f(t 2 )j cjt 1 t 2 j ). At last, let ( l ) l2I L denote a sequence of real random variables such that P 1 p=1 P(j lp j = 2 l 2 p ; l p ]) < +1: This last condition is satis ed when, for instance, the l are square integrable, have the same distribution and a bounded probability density. It follows from the Borel-Cantelli Lemma that there exists an event , with probability 1, such that for all ! 2 there exists p 0 satisfying for all p p 0 , l 2 p j lp (!)j l p : Set = 4k, where k 2 IN . We will show that the random, lacunary Fourier series, fW(t)g t2 0;1] de ned for all t 2 0; 1] and ! 2 as W(t; !) = P l2I L lHl(t) sin( l t) l (!) and W(t; !) = 0 else, satis es almost surely, for all t, W (t) = lim inf H l (t). Thus, this construction allows to build stochastic processes with an arbitrary singularities spectrum. These processes are not obtained via a multiplicative cascade and yet they can be multifractal.

Introduction

Pour simpli er, nous supposons que tous les processus dont il est question dans cette note sont index es par 0; 1] et a valeurs r eelles. L'extension de nos r esultats a des processus a valeurs r eelles index es par un pav e compact de IR d , ne posent pas de probl emes essentiels. Par ailleurs, nous nous int eressons a des processus dont les trajectoires sont presque sûrement des fonctions continues et nulle-part d erivables. Soit X un processus de ce type, le processus de H older associ e a X, est d e ni pour tout t par X (t) = sup n lim sup h!0 jX(t+h) X(t)j jhj = 0 o : Il permet de mesurer les variations locales de l'irr egularit e de X. Dans 4], il a et e montr e qu'une fonction H est la fonction de H older d'une fonction f d eterministe et continue si et seulement si H s' ecrit comme une limite inf erieure de fonctions continues i.e pour tout t, H(t) = lim inf H l (t) o u les H l sont continues (on peut supposer, sans perte de g en eralit e, que les H l sont des fonctions H olderiennes avec des constantes de H older croissant lentement i.e en O(l), voir 4]). Une construction possible repose sur une g en eralisation de la c el ebre fonction de Weierstrass. L'objectif de cette note est de transposer ce r esultat au cadre stochastique. Plus pr ecis ement, nous construisons un processus al eatoire continu, g en eralisant la fonction de Weierstrass et dont le processus de H older peut être, avec probabilit e 1, n'importe quelle limite inf erieure de fonctions continues, a valeurs dans un compact de ]0; 1 .

2 Construction de W Dans toute cette partie a, b, a 0 , b 0 sont des r eels tels que 1 > b 0 > b > a > a 0 > 0 et IL = (l p ) p2I N est la suite d'entiers d e ni par l 1 = 1 et pour tout p 1 par l p = 1 a 0 1 b 0 l p 1 ] + 1. (H l ) l2I L d esigne une suite de fonctions ( ; c l )-H olderiennes a valeurs dans a; b] et telles que c l = O(l) (on dit qu'une fonction f est ( ; c)-H olderienne si pour tout t 1 , t 2 on a jf(t 1 ) f(t 2 )j cjt 1 t 2 j ). f" l g l2I L d esigne une suite de variables al eatoires r eelles telles que P 1 p=1 P(j lp j = 2 l 2 p ; l p ]) < +1: Cette derni ere condition est par exemple v eri ee lorsque les " l sont identiquement distribu ees, admettent un deuxi eme moment et une densit e born ee. Ainsi d'apr es le Lemme de Borel-Cantelli, il existe un ev enement , de probabilit e 1, tel que : pour tout ! 2 , il existe p 0 , v eri ant pour tout p p 0 l 2 p j" lp (!)j l p :

(1) Ces in egalit es joueront un rôle essentiel dans tout ce qui va suivre. Tous les processus consid er ees ci-dessous seront suppos es nuls en tout ! = 2 . Notre objectif sera d' etablir le Th eor eme suivant.

Th eor eme 1 Soit un entier naturel non nul divisible par 4 et soit fW(t)g t2 0;1] , la s erie de Fourier al eatoire lacunaire, d e nie lorsque ! 2 par W(t; !) = X l2I L lHl(t) sin( l t)" l (!);

(2) alors pour tout ! 2 et tout t 2 0; 1], on a W (t; !) = lim inf p!1 H lp (t). Fixons un r eel " > 0 arbitrairement petit, t 2 0; 1] et ! 2 . Il existe p 1 p 0 tel que pour tout p p 1 , H lp (t) H(t) . Posons pour tout s 2 0; 1], W p1 (s; !) = P p1 p=1 Hl p (s)lp sin( lp s)" lp (!) et W(s; !) = W(s; !) W p1 (s; !). Il r esulte du Th eor eme des accroissements nis qu'il existe C 1 > 0 une constante telle que pour tout s 0 ; s 00 2 0; 1], on a jW p1 (s 0 ; !) W p1 (s 00 ; !)j C 1 js 0 s 00 j :

Etape 1 Montrons qu'il existe une constante C 2 > 0, telle que pour jhj 1. Pour tout N 2 e " 4 posons h N = ( 1) (N) 2 N o u (N) 2 f0; 1g est tel que j sin( N (t + h N )) sin( N t)j 1:

j f W(t + h; !) f W(t; !)j C 2 jhj
(5) On a, en utilisant (1) et le fait que h N = ( 1

) (N) 2 N , j f W(t + h N ; !) f W(t; !) NHN (t) (sin( N (t + h N )) sin( N t))" N (!)j A(h N ) + A 0 (h N ) + A 00 (h N ) o u A(h) = X l2e I Lne "
lHl(t) j sin( l (t + h)) sin( l t)jl A 0 (h) = X l 2 e " l < N lHl(t) j sin( l (t + h)) sin( l t)jl et A 00 (h) = X l2I L j lHl(t+h) lHl(t) jl:

On en d eduit que A(h) +1 X n=0 n(H(t)+") nj sin( n (t + h)) sin( n t)j C 7 +1 X n=0 n(H(t)+ " puisque sup n2I N n"
2 n < +1. Ensuite en raisonnant comme dans l' etape 1, on peut montrer que A(h) C 8 jhj H(t)+ " 2 : (6) Soit L N = supfl 2 e " ; l < Ng: On a A 0 (h) C 8 LN X n=0 n(H(t) 3 2 ") j sin( n (t+h)) sin( n t)j C 8

(1 H(t)+ [START_REF] Barral | Continuity of the multifractal spectrum of a random statistically self-similar measure[END_REF] 2 ")(LN +1) (1 H(t)+ 3 2 ") 1 jhj:

(7) Nous supposerons que [START_REF] Barral | Continuity of the multifractal spectrum of a random statistically self-similar measure[END_REF] 2 " < min(a a 0 ; b 0 b). On a donc 1 a 0 > 1 (H(t) [START_REF] Barral | Continuity of the multifractal spectrum of a random statistically self-similar measure[END_REF] 2 ") > 1 (H(t) + 3 2 ") > 1 b 0 : Par suite, (1 (H(t)

+ 3 2 "))N (1 (H(t) 3 2 "))L N , et (7) entra^ ne alors A 0 (h N ) C 0 9 N(1 (H(t)+ 3 2 ")) jh N j C 9 jh N j H(t)+ 3 2 " :
(8) En raisonnant comme dans l' etape 1, on peut montrer que jA 00 (h)j C 10 jhj:

(9) Finalement il r esulte de (1), ( 5), ( 6), ( 8) et (9) que, pour tout N 2 e " 4 , assez grand

j f W(t + h N ; !) f W(t; !)j N 2 N(H(t)+ " 4 ) j sin( N (t + h N )) sin( N t)j A(h N ) A 0 (h N ) A 00 (h N )
C 11 jh N j (H(t)+ " 3 ) C 8 jh N j H(t)+ " 3 C 9 jh N j H(t)+ 3 2 " C 10 jh N j > C 12 jh N j (H(t)+ " , pour t = P 1 i=1 t i 2 i ; t 1 = 0 ou 1 et t 2 I k n ). ( n ) n converge faiblement vers une mesure appel ee mesure binômiale. et sa fonction de r epartition ont un spectre multifractal en forme de \, et la fonction de H older associ ee est partout discontinue. Elle est donn ee par H(t) = ' 0 (t) log 2 m 0 ' 1 (t) log 2 m 1 , o u ' 0 (t) = lim inf n!1 P n i=1 (1 ti) n et ' 1 (t) = 1 ' 0 (t). On peut g en eraliser cette construction en rendant les \poids" m 0 et m 1 al eatoires, de fa con a obtenir un processus stochastique multifractal, c'est-a-dire dont le spectre n'est pas r eduit a un point. On parle alors de cascade al eatoire, et ce type de processus est presque syst ematiquement invoqu e quand il s'agit de mod eliser des ph enom enes stochastiques multifractals.

De ce point de vue, l'avantage de la construction pr esent ee dans cette note est qu'elle permet d'obtenir, de fa con el ementaire, un processus dont le spectre multifractal est le plus g en eral possible, puisque le processus de H older de W est de la forme la plus g en erale. Ainsi, pour obtenir un spectre du même type que celui d'une mesure binômiale, c'est-a-dire en \, il su t de prendre une suite de fonctions continues, dont la limite inf erieure est de la forme H(t) = ' 0 (t) log 2 m 0 ' 1 (t) log 2 m 1 . Mais des spectres tr es di erents sont aussi possibles. Notons que l'obtention d'un spectre prescrit est crucialement li ee au fait que notre construction garantit que, presque surement, pour tout t, l'exposant de H older est bien H(t), et pas seulement le r esultat moins fort que, quel que soit t, presque surement l'exposant de H older est H(t). En particulier, les ensembles E sont, dans notre cas, \d eterministes", alors qu'ils sont al eatoires pour la mesure binômiale stochastique. On pourrait rendre les E pour W al eatoires en choisissant un H al eatoire. D'ailleurs, la d etermination du spectre presque sûr des cascades al eatoires est beaucoup plus di cile que pour W 2].

En r esum e, notre construction fournit un proc ed e simple pour obtenir un processus al eatoire dont le processus de H older peut être prescrit dans le cadre le plus g en eral, ce qui autorise en principe l'obtention de n'importe quel spectre multifractal presque sûr. Cette construction n'est pas a base de cascade, et elle pourrait donc permettre de mod eliser des ph enom enes multifractals dans lesquels aucun m ecanisme de cascade ne semble a l' uvre (comme par exemple certains tra cs sur Internet).

Remeciements

Les auteurs tiennent a remercier vivement Julien Barral et Jacques Peyri ere pour leurs importants commentaires sur les versions pr ec edentes de cette note.

R ef erences

  remarques concernent l'analyse multifractale (voir par exemple 3]). En e et, dans le cas o u H est tr es irr eguli ere, ce qui est possible avec la construction ci-dessus, on pr ef ere souvent d ecrire la r egularit e de W de fa con globale, a travers son spectre de singularit e. Rappelons que ce spectre est la fonction f qui a tout associe f( ) = dim H E , o u dim H d esigne la dimension de Hausdor et E = ft : H(t) = g. Pour une fonction de Weierstrass standard ou un mouvement Brownien fractionnaire, le spectre est r eduit a un point et ces fonctions sont parfois quali ees de \monofractales". L'exemple paradigmatique de fonction multifractale est la fonction de r epartition d'une mesure binômiale: on se xe un r eel m 0 dans 0; 1=2] et on d e nit la suite de mesures ( n ) par n n ; (k +1)2 n ; ' 0 (k; n) est la proportion de 0 dans le d eveloppement en base 2 d'un r eel de I k n jusqu' a l'ordre n, et ' 1 (k; n) = 1 ' 0 (k; n) (i.e. ' 0 (k; n) = P n i=1 (1 ti) n

) j sin( n (t + h)) sin( n t)j;


Development of innovative pH sensor to evaluate phagocytosis of nanoparticles

Introduction: Inhaled nanoparticles (NP) exhibit variable toxicity levels which mainly depend on their physicochemical characteristics (size, morphology, crystallinity, chemical surface composition…). Biological effects monitoring thanks to usual tests (ROS, TNFα, LDH) are performed on alveolar macrophages collected from the respiratory system. In this frame, evaluation of NP uptake by macrophages appears as complementary data useful for NP toxicity assessment. The aim of this work deals with the development of pH sensible NP allowing the quantification of NP phagocytosed by macrophages, specially the step of fusion between phagosomes and lysosomes.

Material and methods: Two types of fluorescent NP with variable and well-characterized sizes and surface coatings have been synthesized. One type with FITC fluorescence and one other with two fluorochromes (FITC within the core and pHrodo in the porous polysiloxane shell). Red fluorescence of pHrodo probe increases as pH acidifies allowing the distinction of intra-lysosomal engulfed NP in macrophages (cell line RAW 264.7). Observations were realized by confocal microscopy performing fluorescence spectral analysis. Confocal acquisitions were realized with the two types of NP including fluorescent labeling of cell nuclei with hoechst as well as actin labeling with phalloidin for cellular limit detection.

Results: Incubation of pHrodo NP in different acellular mediums with decreasing pH values showed the emission peak of red fluorescence in acidic conditions. Moreover confocal microscopic observations led to distinguish entirely engulfed NP (yellow labeling due to colocalisation of NP in acidic vesicles) from those which are just adherent to the cell membrane (FITC green labeling).

Conclusion/Perspectives:

This study allows the validation of pHrodo NP model and the confocal observation of NP in contact with macrophages. We already developed cytometric quantification methods for micrometric particles, which can now be applied to NP.

Key words: Fluorescent nanoparticles, phagocytosis pH sensor, confocal microscopy.


François Buret

Monique Dauge

Patrick Dular

Laurent Krähenbühl
email: laurent.krahenbuhl@ec-lyon.fr

Victor Péron

Ronan Perrussel

Clair Poignard

Damien Voyer

2D electrostatic problems with rounded corners

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

High-voltage applications require a precise knowledge of the electric field in the area where the geometry of the structure is sharp. On a real device, the geometry is not "exactly sharp" but present rounded edges or corners. The accurate description of these rounded shapes, especially when the geometry involves several corners, can be cumbersome in a numerical model. In addition most of the time only a rough (statistical) description of the rounded shape is available due to the manufacturing tolerance. Dauge et al. have proposed in [START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF] a theoretical approach to tackle the "rounded shape" problems by an accurate asymptotic analysis. Define the potentials v in Ω and v 0 in Ω (see Fig. 1) by

Γ N Γ 1 Γ 1 Γ 0 Ω (b) Domain with a sharp corner Ω. O 1 θ R ω Γ 0 ∞ Ω∞ (c) Unbounded profile domain Ω∞.
         v = 0, on Ω v = 0, on Γ 0 v = 1, on Γ 1 ∂ n v = 0, on Γ N ,          v 0 = 0, on Ω v 0 = 0, on Γ 0 v 0 = 1, on Γ 1 ∂ n v 0 = 0, on Γ N , (1)
where characterizes the "size" of the rounded corner and ∂ n denotes n • ∇, n being the unitary outward normal on the boundary of the domain. Throughout the paper ω denotes the angle of the sharp corner and α = π/ω. The main idea of [START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF] consists in expanding v into two sums in power of α . Using a smooth radial cut-off function ϕ defined by

ϕ(ρ) = 1, if ρ d 1 0, if ρ d 0 , with d 0 < d 1 , (2)
d 0 , d 1 being fixed corner distances, the expansion writes for any integer n 1 (ϕ(./ ) is the function t → ϕ(t/ ))

v =ϕ . v 0 + ϕ . n p=1 b p pα v pα + (1 -ϕ) n p=1 B p pα V pα . + r nα , (3)
where b p and B p are real parameters and r nα is such that

∃˜ > 0, ∃C n > 0, ∀ < ˜ , Ω ∇r nα 2 dx < C n (n+1)α , (4) i
.e. the energy norm of the error converges as (n+1)α to 0, which is written r nα = O H1 ( (n+1)α ). The functions v pα satisfy a boundary value problem in the sharp domain Ω (see Fig 1(b)), whereas the functions V qα are the so-called profile terms that satisfy the Laplace equation in the infinite domain

Ω ∞ of R 2 , which is a localization of the rounded corner (see Fig 1(c)).
The paper [START_REF] Krähenbühl | Numerical treatment of rounded and sharp corners in the modeling of 2D electrostatic fields[END_REF] was devoted to provide numerically the first order terms of the theoretical expansion of [START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF]. The following remarks have been observed:

• the exact solutions close to the corner, computed for several values of the curvature radius , are quasi-similar, up to a "scaling factor" (related to ). It is also noticed that the "shape" of the solutions close to the corner (their "shape" but not their amplitude) weakly depend on other elements of the studied structure, such as the distance to the boundaries: if the corner geometry is self-similar 1 , it is also said that the dominant term of the solutions close to the corner is self-similar. • the exact solutions far from the corner are weakly influenced by the change of the curvature radius , and they converge to the solution on the domain with the sharp corner when goes to zero.

The aim of the present paper is to push forward numerically the expansion (3) on the numerical model studied in [START_REF] Krähenbühl | Numerical treatment of rounded and sharp corners in the modeling of 2D electrostatic fields[END_REF], especially when non-symmetric structures are involved.

II. SECOND ORDER EXPANSION

The heuristics for the construction of the two first orders are detailed in [1, subsection 4.1]. The roughest approximation v 0 of v , far from the corner, that is defined in Ω by [START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF], writes in the neighborhood of the corner

v 0 (x) ρ→0 ∞ p=1 a p ρ pα sin(pαθ) = ∞ p=1 a p s pα (ρ, θ), (5)
where (ρ, θ) are the polar coordinates. As the behaviors of v 0 and v are different in the corner, v 0 should be truncated in the corner. Expansion (5) enforces the coefficient B 1 of (3) to equal the first singular coefficient a 1 , while the profile term V α satisfies, for p = 1,

       -X V pα = [ X ; ϕ]s pα , in Ω ∞ , V pα | Γ∞ = 0, lim R→+∞ V pα = 0. (6)
For any couple (ν, u), we remind that [∆; ν]u = ∆(νu)-ν∆u.

In [START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF], it is especially shown that in the neighborhood of +∞

V α (X) ρ→+∞ +∞ p=1 A p s -pα (X). (7)
Then, necessarily, b 2 and B 2 respectively equal a 1 A 1 and a 2 given by ( 5) and ( 7), leading to the second order:

v = ϕ • v 0 + α (1 -ϕ)a 1 V α • + 2α a 1 A 1 v 2α + (1 -ϕ)a 2 V 2α • + O H 1 ( 3α ), (8)
where V 2α is the profile term that satisfies (6) for p = 2 and v 2α is the correction far from the corner defined by

     -v 2α = [ ; (1 -ϕ)]s -α , on Ω, v 2α = 0, on Γ 0 , and v = 0, on Γ 1 , ∂ n v 2α = 0, on Γ N . (9)
In [2, page 9], the profile problem (25) satisfied by v ∞ did not involve the cut-off function ϕ as in (6). Actually, the profiles V α and v ∞ can be linked through the equality V α = v ∞ -ϕs α . However, for profiles of order p > 2, the relationship is more complex and the use of the cut-off function seems to be more convenient [START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF].

III. NUMERICAL RESULTS

Two L-shape geometries with a rounded corner are considered; their dimensions are specified on Fig. 2. For both cases, α equals 2/3. The finite element method has been used for solving (1), (6), and (9) as it was done in [START_REF] Krähenbühl | Numerical treatment of rounded and sharp corners in the modeling of 2D electrostatic fields[END_REF].

The errors in the energy norm of the first and second order expansions, plotted on Fig. 3, behave respectively as 2α and 3α , independently of the geometry. This is remarkable since a 2 = 0 in the symmetric configuration (see Fig. 2(a)) and we could think that the first order approximation would have lead to a better approximation than in the non symmetric configuration (see Fig. 2(b)). Nonetheless, this intuition is obviously not correct and the correction v 2α far from the corner plays an equivalent role than the correction close to the corner regarding the energy norm in the whole domain.

The normal electric field along the electrodes in the non symmetric configuration (from the bottom right to the top left of the electrode, see Fig. 2(b)) for two values of are presented in Fig. 4. This normal field has been computed respectively from v , from first order expansion and from second order expansion given by (8). The behavior of the normal field is "closer" for the order 2, in particular the location of the maximum is roughly equivalent to the correct solution. However the maximum electric field is overestimated by both approximations, requiring to go further in the expansion. Curvilinear abscissa on Γ 0 . Normal electric field (V/m).

Exact. Order 1. Order 2.

(a) Rounded corner with = 0.4. Origin for the curvilinear abscissa is at the center of the rounded corner.

Fig. 1 .

 1 Fig. 1. Considered domains Ω , Ω and Ω∞.

Fig. 2 .

 2 Fig. 2. Two considered problems and their dimensions. = 0.4.

Fig. 3 .

 3 Fig. 3. Convergence of the error in the energy norm.

  Rounded corner with = 0.2.

Fig. 4 .

 4 Fig.4. Normal electric field along Γ 0 for the exact solution and the two first order approximations. Non symmetric configuration is considered (Fig.2(b)). Origin for the curvilinear abscissa is at the center of the rounded corner.

O θ ρ ω Γ N Γ N Γ 1 Γ 1 Γ 0 Ω (a) Domain with a rounded corner Ω . O θ ρ ω Γ N

Roughly, it means that a single parameter, here , and a basic geometry are sufficient to describe the corners for any value of . For a precise definition, refer to[START_REF] Dauge | Around the Research of Vladimir Maz'ya II: Partial Differential Equations[END_REF].


J Wolff

C M Marques

F Thalmann
email: fabrice.thalmann@ics-cnrs.unistra.fr

A thermodynamic approach to phase coexistence in ternary phospholipid-cholesterol mixtures

Keywords: numbers: 87.16.dt, 64.75.-g, 87.14.Cc

We introduce a simple and predictive model for determining the phase stability of ternary phospholipid-cholesterol mixtures. Assuming that competition between the liquid and gel order of the phospholipids is the main driving force behind lipid segregation, we derive a Gibbs free-energy of mixing, based on the thermodynamic properties of the lipids main transition. A numerical approach wax devised that enable the fast and efficient determination of the ternary diagrams associated with our Gibbs free-energy. The computed phase coexistence diagram of DOPC/DPPC/cholesterol reproduces well known features for this system at 10 • C, as well as its evolution with temperature.

Lipid membranes have captured the attention of physicists, chemists and biologists alike owing to their prominent structural and functional role in living cells, and to the pervading use of lipid based products in the pharmaceutical, food and cosmetic industries [START_REF] Mouritsen | Life-as a matter of fat: the emerging science of lipidomics[END_REF][START_REF] Heimburg | Thermal Biophysics of Membranes[END_REF]. In the recent years, ternary lipid systems comprising phospholipids and cholesterol were established as important model systems in relation with lateral lipid segregation phenomena in membranes. At the root of this rich phenomenology lies a complex phase behavior of certain lipid mixtures, in which cholesterol molecules play a unique role. Ternary systems comprising a phospholipid in the liquid state, a second phospholipid in the gel state and cholesterol display several unique characteristics, the two most important of them being arguably a wide binary coexistence gap between the so-called liquid disordered (L d ) and liquid ordered (L o ) phases, and a region of triple coexistence [START_REF] Zhao | [END_REF][4][5].

A prominent feature of phospholipid membranes is their main thermodynamic transition, from a low temperature gel or rippled gel phase (usually termed L β ,P β ), to a high temperature liquid phase (L α ) [START_REF] Heimburg | Thermal Biophysics of Membranes[END_REF]. A sharp change in the statistical conformations of the hydrophobic chains occurs at the transition, associated with membrane thinning and an increase of the area per phospholipid. The transition takes place at a melting temperature T m whose precise value depends on the chemical structure of the chains. It is a weakly first order transition for membrane of pure composition. The thermotropic behavior of phospholipid binary mixtures have been studied since the seminal work of Mabrey and Sturtevant [6] and supports the view that liquid-prone and gel-prone lipids tend to phase separate. This suggests that mixtures composed of phospholipids with same head group have their miscibility behavior governed by the competition between their liquid and gel conformations. Hence, a description of ternary lipid systems must necessarily, and to some extent, account for the internal state of the phospholipids.

We introduce in this Letter a simple, robust and pre-dictive model of phospholipid-cholesterol mixtures phase coexistence. Our model aims at predicting the miscibility behavior of phospholipids, based upon the thermodynamic properties of their main transition [7][8][9][10]. The thermotropic behavior of the phospholipid tails is modeled with an internal order parameter, in a way that schematically reproduces the thermodynamic features of the liquid-gel transition. The peculiar features of cholesterol-phospholipid interactions [11,12] are accounted for by means of a non-ideal mixing entropy generalizing the one proposed by Komura et al. [13]. In addition, we outline an original numerical procedure for the derivation of ternary phase diagrams.

The stability of a mixture is governed by its Gibbs freeenergy of mixing G mix . This thermodynamic function expresses the free-energy difference between the homogeneous mixed state and the system comprising the separated pure components with same global composition. Ternary mixtures in membranes are therefore described by G mix (T, n 1 , n 2 , n 3 ), where the n i are the molar content of each lipid, and T the temperature. We assume that the hydrostatic pressure remains constant, close to the standard value, and that lipid membranes have no tension, so that these two intensive parameters are absent from G mix . Membranes are assumed to be well hydrated with a large excess of water, and free from any other external fields [14].

The simplest realization of our model, described in this Letter, assumes that the areas per lipid a i do not change with temperature. This allows us to express G mix in terms of the surface fractions φ 1 , φ 2 and φ C of each lipid component, φ C standing for the cholesterol. We claim that it is acceptable to neglect all changes in the area per lipid at T = T m , in a first approximation to the description of the mixture phase coexistence, especially in the absence of membrane tension. This restriction is not a fundamental one, and versions of our model accounting for membrane expansion at the transition will be presented elsewhere [15].

A frequent approximation to non-ideal solutions is obtained by combining an ideal mixing entropy with a quadratic enthalpic term. For a membrane of area A, the free-energy of mixing per unit of area, g mix = βG mix /A, reads:

g mix = φ 1 a 1 ln(φ 1 ) + φ 2 a 2 ln(φ 2 ) + φ C a C ln(φ C ) +χ 12 φ 1 φ 2 + χ 1C φ 1 φ C + χ 2C φ 2 φ C , (1)
with φ 1 +φ 2 +φ C = 1 and β the Boltzmann factor 1/k B T , and a 1 , a 2 and a C are the areas of the three different lipid species. One observes that expression (1) alone does not account for the main transition, and that the effective interaction parameters χ ij are expected to display a non trivial and unknown temperature dependence. We could not obtain, based on (1), any diagram comparable with the experimental ones. Moreover, our attempts to do so lead to unexpectedly large values of a C , conflicting with the view that cholesterol molecules are smaller than phospholipids.

Experiments and simulations either suggest that cholesterol may form transient bound complexes with phospholipid moieties surrounding them [11,12]. The presence of such complexes prevents cholesterol molecules from occupying too many neighboring positions, thus giving the semblance of a large excluded volume. To account for the statistical behavior of cholesterol molecules, we derived the following entropic contribution, where all lipid areas have been set to unit:

G (E) (φ 1 , φ 2 ) = φ 1 ln(φ 1 ) + φ 2 ln(φ 2 ) + φ C ln(φ C ) -(φ 1 + φ 2 ) log(φ 1 + φ 2 ) + (δ -1 -φ C ) log(1 -δφ C ).
(

) 2
This expression originates from a lattice enumeration of configurations, following Flory, where each cholesterol molecule deprives an average number δ of neighboring sites from the possibility of being occupied by another cholesterol molecule. Expression (2) reduces to ideal mixing when δ → 1, and is similar to the entropy of Komura et al. for binary systems and δ = 2 [13]. It accounts for both the mixing term and the excluded volume between cholesterol molecules, consistent with the presence of cholesterol complexes. It finally sets a maximal value φ C ≤ 1/δ to the mixture content in cholesterol. Several Landau order parameters have been proposed for describing the liquid-gel chain melting [7,8,10,13,16]. Theories based on a single scalar order parameter can pretend to capture only rough features of the main transition. This is why, for simplicity and generality, we restrict ourselves to a schematic approach, a two-state model based on a scalar order parameter m restricted to the interval [-1, 1]. The mean-field Ising model under finite external field serves as a guide [START_REF] Chaikin | Principle of condensed matter physics[END_REF].

G (I) (T, φ 1 , φ 2 , m) = -m{h 1 (T )φ 1 + h 2 (T )φ 2 } -2 Jm 2 + ( 1 + m 2 ) ln( 1 + m 2 ) + ( 1 -m 2 ) ln( 1 -m 2 ), (3)
with

h 1 (T ) = ∆H 1 2RT 2 1 (T -T 1 ) ; h 2 (T ) = ∆H 2 2RT 2 2 (T -T 2 ), ( 4
)
and R is the gas constant. Eq. ( 3) represents a meanfield Ising model in the presence of an effective external field ) is shaped as a double well potential, with a linear bias induced by h eff . The order parameter m fluctuates freely around an optimal value m ⋆ located at the global minimum of G (I) , with (T, φ 1 , φ 2 ) constant.

h eff = h 1 φ 1 + h 2 φ 2 . For J ≥ Jc = 1/4, the m dependence of G (I
In practice, as soon as J is larger than 0.3, the two local minima are close to ±1. We conventionally ascribe m g ≃ -1 to represent the L β gel phase, while m l ≃ 1 is set to represent the L α liquid phase. The coupling J expresses that neighboring lipid tails tend to be in the same state, irrespective of the chemical species to which they belong. This is the main cooperative effect that induces the liquid-gel separation.

The term h eff acts as a switch that drives the order parameter from its m ⋆ = m g gel state to its m ⋆ = m l liquid state, upon changing its sign from negative to positive. One sees from (4) that if lipid 1 is pure, h 1 (T ) changes its sign at T = T 1 while m ⋆ (T ) suddenly jumps from m g to m l , with an amplitude ∆m ⋆ ≃ 2. The enthalpic change is close to ∆H 1 , so that ∆H 1 can be identified with the molar latent heat at the main transition. The difference in Gibbs free-energy between the liquid and gel phases at T = T 1 is then extrapolated linearly with T at the vicinity of the transition. The contribution G (I) is thus entirely specified by the temperature and the enthalpy difference at the transition of the lipid species composing the mixture. This term induces a liquid-gel separation in the whole temperature range [T 1 , T 2 ].

Cholesterol-phospholipid interactions are further accounted for with a specific interaction term

G (C) (φ 1 , φ 2 , m) = χφ 1 φ C -ξm(φ 1 + φ 2 )φ C (5)
Eq. ( 5) comprises a m-dependent contribution, proportional to ξ, which marks the preference of cholesterol for either the m = 1 liquid order, or the m = -1 gel order, according to the sign of ξ. This interaction does not discriminate between saturated and unsaturated lipids, since it describes a coupling to the chain order. Other cholesterol-phospholipid interactions are accounted for by a m-independent contribution in χ. This order-insensitive interaction introduces a coupling between cholesterol and the unsaturated component. We found that interactions between cholesterol and the saturated component are not required to describe the experimental phase diagrams.

The sum G(T,

φ 1 , φ 2 , m) = G (E) + G (I) + G (C
) is our main thermodynamic potential. In the absence of surface tension, and assuming that changes in internal state occur at constant volume, m fluctuates around its optimal value m ⋆ , with m = m ⋆ in the thermodynamic limit. Therefore, the free-energy of mixing g mix (T, φ 1 , φ 2 ) is obtained by minimizing G(T, φ 1 , φ 2 , m) with respect to m at fixed temperature and surface fractions [10].

g mix (T, φ 1 , φ 2 ) = min m [G(T, φ 1 , φ 2 , m)]| T,φ1,φ2 = G(T, φ 1 , φ 2 , m ⋆ (T, φ 1 , φ 2 )). (6)
The Gibbs free-energy surface g mix arises from projecting out m from the four-dimensional potential G. Equivalently, g mix can be viewed as the restriction of G to the submanifold ∂ m G = 0. For J > Jc , the constraint ∂ m G = 0 defines two stable submanifolds, a negative m g (T, φ 1 , φ 2 ) and a positive m l (T, φ 1 , φ 2 ) branch, linked respectively with the gel and liquid phases. At fixed (T, φ 1 , φ 2 ), only one of the two branches is the stable one, while the other branch remains metastable.

In ternary mixtures it is known that phase coexistence is dictated by the convexity properties of g mix (T, φ 1 , φ 2 ) relative to φ 1 , φ 2 . Minimization of g mix is equivalent to finding gmix , the lower convex hull of g mix , that displays regions of three kinds. Points of the convex hull in contact with the original surface are stable monophasic regions. Regions lying above triple tangent planes corresponds to triple coexistence. The remaining part of the convex hull correspond to developable patches of surface, wrapped around double tangent lines that connects pairs of points with distinct composition (tie-lines) [START_REF] Gibbs | The Scientific Papers of J. Willard Gibbs[END_REF]. In practice, we used a public domain routine, qhull, to compute a discrete approximation of gmix [START_REF] Barber | [END_REF]. The starting point is a fine mesh discretization of g mix over the relevant domain of composition {φ 1 , φ 2 ≥ 0, φ 1 + φ 2 ≤ 1 -1/δ}. Then qhull computes a triangulated surface approaching gmix . Triple coexistence correspond to facets with all sides much larger than the initial mesh size. Double coexistence is associated with elongated triangles, with a shortest side much smaller than the two longer sides, the latter being oriented parallel with the tie-lines. Finally, small facets of the convex minimization are linked with stable monophasic regions. The projection of the triangulated gmix surface onto the composition plane provides directly a fast and accurate approximation of the Gibbs phase diagram shown in fig. 1(a). This discrete solution can then be refined according to the usual rules for the determination of ternary phase diagrams. Details of these calculations will be given in a forthcoming publication [15].

Arguably, the most widely studied ternary lipid phase diagram for membranes containing cholesterol is the ternary mixture DOPC/DPPC/Chol. DOPC is a double tail molecule with 18 carbons per tail and a cisunsaturated bond on the ninth carbon. It has a low main transition temperature T 1 at -21 • C and an associated enthalpy change ∆H 1 evaluated at 7.7 kcal/mol [START_REF] Cevc | Phospholipids Handbook[END_REF]. DPPC has also a phosphatidylcholine head but here the two 16-carbon tails have no unsaturated bonds, leading to a higher transition temperature T 2 at 42 • C, and associated enthalpy change ∆H 2 = 8.7 kcal/mol [START_REF] Cevc | Phospholipids Handbook[END_REF]. Within our approach, the above temperature and enthalpy values and a single coupling parameter J = 0.35 -see eq. ( 3) -describe well published data [START_REF] Schmidt | [END_REF] for the experimental binary DOPC/DPPC mixing behavior. For a complete description of the ternary phase diagram our model further requires three parameters: δ, ξ and χ. The value of δ is bound by the maximum amount of cholesterol that a membrane can hold, estimated here at φ C ≃ 0.45 or equivalently δ = 1/φ C = 2.2 [START_REF] Zhao | [END_REF]. The parameter ξ, as defined in eq.( 5), quantifies the preference of cholesterol for the liquid disordered phase with respect to the gel phase [START_REF] Mouritsen | Life-as a matter of fat: the emerging science of lipidomics[END_REF]22]. Binary diagrams of DPPC/Chol computed with ξ = 1.8 reproduce well the temperature behavior of the miscibility gaps occurring between the L β -gel phase and the L o liquid ordered phase as determined from Giant Unilamellar Vesicles in [23] and from NMR [24] -see fig. [START_REF] Heimburg | Thermal Biophysics of Membranes[END_REF]. Eq.( 5) also includes a repulsive m-independent interaction between cholesterol and the unsaturated component, given by χ. At 10 • C a value of χ = 6.0 promotes both a triple coexistence region and a large biphasic domain, with oblique L o /L d tie lines of positive slope and a plait point on the left side of the Gibbs triangle, in agreement with experimental findings [5] [4]. A second critical point, initially hidden by the triple coexistence region, is predicted to emerge at the binary detachment temperature. Note also that within our approach the extension of the L β domain, on the right bottom part of the diagram, gradually decreases up to T 2 .

As a summary we presented a model for phase coexistence in ternary phospholipid-cholesterol mixtures. Our approach is based on a mean-field two-state description of the phospholipids and a non-ideal mixing entropy for cholesterol. All but one parameter in our system are bound by known thermodynamic or physical properties. We combined analytical theory and numerical convex minimization to show the predictive power of our model for the ternary mixture DOPC/DPPC/Chol. Resulting phase coexistence diagrams reproduce well known features for this system at 10 • C as well as the temperature variation of the diagram. We are also encouraged by preliminary results from related ternary systems that further indicate that our model successfully describes a large variety of lipid/lipid/cholesterol mixtures.

FIG. 1 :

 1 FIG. 1: (Color online) (a) Discrete approximation to the ternary phase diagram obtained from the lower convex hull method explained in the text. (b) Final phase diagram determined from the analytical refinement of the discrete approximation.

  FIG.2:(Color online) Computed evolution of the DPPC/DOPC/Chol triangular diagram with temperature. The computed diagrams display the main expected features for this lipid mixtures.

  . A larger value of χ would open there a miscibility gap, unseen for DOPC/Chol, while a smaller value would drastically reduce the biphasic L o /L d domain. In our model a continuous path connects the L o and L d regions around the critical point, where the order parameter m remains positive. The L β phase on the lower right corner of the phase diagram is associated with a negative m-value and separated from the rest of the diagram by a discontinuity in m across the coexistence gap. The evolution of the Gibbs diagram with respect to temperature is displayed in fig. (2). The diagram was drawn here under the assumption of a constant ξ parameter, and a constant value of the coupling contribution k B T • χ denoting its pure enthalpic nature. Our diagram at the lower temperature in fig. (2), displays the main expected features for a DOPC/DPPC/Chol mixture. As the temperature increases, the triple coexistence triangle shrinks and vanishes between 15 • C and 20 • C. Also, the L o /L d biphasic region, resting initially on the left side of the triangle, detaches and assumes a closed shape of decreasing surface. A first critical point, on the left side of the Gibbs diagram, evolves slowly away from the edge


We consider the Block Cimmino projection method for the solution of the linear system of equations Ax = b, where A is a nonsingular large sparse unsymmetric matrix of order n. The blocks are obtained by partitioning the system into p strips of rows, with p ≤ n, as in:

   A 1 . . . A p    x =    b 1 . . . b p    . (1)
The Block Cimmino method projects the current iterate simultaneously onto the manifolds corresponding to the strips, and takes a convex combination of all the resulting vectors. Our aim is to focus on preprocessing and partitioning strategies of the original system to bring it to the form (1).

If the matrix A is ill conditioned then there are some linear combinations of rows that are almost equal to zero and, as mentioned by [START_REF] Bramley | Row projection methods for large nonsymmetric linear systems[END_REF], these linear combinations may occur inside the blocks or across the blocks after row partitionings of the form [START_REF] Arioli | A Block Projection Method For Sparse Matrices[END_REF]. Assuming that the projections in the Block Cimmino algorithm are computed exactly on the subspaces, then the rate of convergence of the Block Cimmino algorithm depends only on the conditioning across the blocks. If we consider additionally Conjugate Gradient acceleration of the Block Cimmino method, as in [START_REF] Arioli | A Block Projection Method For Sparse Matrices[END_REF][START_REF] Bramley | Row projection methods for large nonsymmetric linear systems[END_REF], the convergence behaviour of the resulting method is directly linked to the spectrum of the n × n matrix E RJ formed as the sum of the previously mentioned projections, and given by

E RJ = p i=1 A i T A i A i T -1 A i , (2)
assuming, for simplicity, that the block-rows A i have full row rank. An efficient implementation of Block Cimmino requires a combination of a robust method for computing the projections and a partitioning strategy that minimizes the ill-conditioning across the blocks.

If the block rows A i are nearly mutually orthogonal, i.e. AA T is strongly block-diagonally dominant, we can expect that the method will converge very quickly, if the projections are computed accurately. Conversely, the structure of AA T tells us about the orthogonality of the subspaces represented by block partitions of A. In all our experiments we first normalize the matrix AA T using MC77 [START_REF] Ruiz | A scaling algorithm to equilibrate both row and column norms in matrices[END_REF] so that the 2-norms of the rows and columns are close to one. Then the entries in AA T correspond to the cosine of the principal angle between every pair of rows, and we may expect that if such a cosine is relatively small the corresponding pair of rows are almost orthogonal. That is, if the block (i, j)th entry of the scaled AA T is close to zero then the subspaces corresponding to the blocks A i and A j can be considered to be orthogonal. If AA T is block tridiagonal, the blocks of A are such that the even numbered blocks are orthogonal as are also the odd-numbered blocks. Thus if we solve the projected subproblems accurately (using say a direct method) then we also solve the subproblems corresponding to the odd and even numbered blocks accurately. The partition of A is thus of the form

[A] = B 1 T B 2 T T ,
where

B 1 = i A i / i odd and B 2 = i A i / i even .
Such a partitioning is called a twoblock partitioning. We denote by m 1 and m 2 the number of rows in B 1 and B 2 respectively. It was shown by [START_REF] Elfving | Block-iterative methods for consistent and inconsistent linear equations[END_REF] that, with such a two-block partitioning, the spectrum of matrix E RJ is

λ k = 1 + cos Ψ k λ k = 1 -cos Ψ k-m 2 λ k = 1 k = 1, . . . , m 2 k = m 2 + 1, . . . , 2m 2 k = 2m 2 + 1, . . . , n
where {Ψ k } m 2 1 are the principal angles between R(B 1 T ) and R(B 2 T ).

The main idea behind the so called two-block partitioning strategy is to exploit structural orthogonality between the subspaces R(A i T ) defined by the partitioning (1). This structural orthogonality can be analysed on the basis of the sparsity pattern of the normal equations matrix AA T . The idea is to find a permutation to transform the normal equations matrix AA T into block tridiagonal form. That is we determine a permutation matrix P such that B = PAA T P T is in block tridiagonal form using an implementation of algorithm of [START_REF] Cuthill | Reducing the bandwidth of sparse symmetric matrices[END_REF]. We then use this permutation to solve the row-wise permuted system of equations

PAx = Pb (3)
using the block Cimmino algorithm.

The main problem with this first preprocessing strategy for general sparse matrices, is that the block tridiagonal structure has diagonal blocks with very differing sizes leading to a partitioning with a bad degree of parallelism because of unbalanced tasks in the block-row projections.

In our second preprocessing strategy, we take into account the numerical values in the matrix AA T , and do not enforce strict orthogonality between every second block of rows. We drop small entries of AA T and by using the relationship to principal angles can relax the orthogonality condition in a controlled way in an attempt to obtain much better parallelism compared to the previous strategy. We thus only keep nonzero entries above a given tolerance and permute the resulting filtered matrix into block tridiagonal form using the Cuthill-McKee algorithm and use the row partitioning defined by this block tridiagonal structure to solve (3).

Of course, since numerical values have been dropped from the normal equations matrix, the resulting partition will not provide two subsets of structurally orthogonal blocks of rows, but, if the values dropped are sufficiently small, we may expect that the numerical properties of the resulting iteration matrix will be relatively close to that of the "strict" two-block partitioning case. Additionally, since the filtered matrix has less entries than the original normal equations matrix, the resulting block tridiagonal permuted matrix will normally have a smaller bandwidth than B and this may help to define a partitioning on A with more blocks, better balanced projections, and a higher degree of parallelism.

To illustrate this preprocessing strategy, we show preliminary results obtained with the bmw3 2 matrix (of size N = 227362) from the Tim Davis collection. The original matrix, before preprocessing, has the symmetric pattern shown in the left subplot above, in which there is no evident structural orthogonality between sets of rows.

We perform some runs on a 24 (Intel i7) cores parallel machine, with a generic partitioning into 24 equal sized blocks. The Block Cimmino method needs then 915 iterations to reach a scaled residual value of the order 10 -7 . It requires 5.46 sec for the analysis, 2.41 sec for the factorization of the augmented systems and 63.25 sec for the iterations.

Following the second preprocessing strategy, we filter AA T matrix under the value 0.1, which corresponds to the cosine of an angle of 84 • , and we obtain the pattern in the right subplot after Cuthill-McKee and a reduction of nonzero entries in AA T of 88%. With this block tridiagonal structure, we can define 18 almost balanced partitions, the biggest one of size 16129 rows and the smallest of 11062 rows. The Block Cimmino method needs then 296 iterations to converge to the same level, with 7.19 sec for the analysis, 3.29 sec for the factorization, 24.62 sec for the iterations, and with 2.48 sec for the extra work in the preprocessing.

We intend to develop and compare the benefits of such strategies on parallel architectures, different test examples and also with variants where we may drop entries in A in combination with entries in AA T

Partitioning

  Strategies for the Block Cimmino algorithm List of authors: M. Zenadi * T. Drummond † I. S. Duff ‡ R. Guivarch * D. Ruiz *