Restoration of SAR Images using Recovery of ... - CiteSeerX
Restoration of SAR Images using Recovery of Discontinuities and non-linear Optimization Florence Tupin, Marc Sigelle, Ammar Chkeif and Jean-Pierre Veran Ecole Nationale Superieure des Telecommunications Departement Images, 46 rue Barrault 75634 Paris Cedex 13 France E-mail: [email protected] Phone: 33 1 45 81 76 27 Fax: 33 1 45 81 37 94
Abstract. In this paper, we study the behaviour of contour recovery
when ltering radar images. We start from recent methods lying on an equivalence scheme between implicit and explicit boundary processes in image restoration [1, 2]. Here we extend them to the processing of synthetic aperture radar (SAR) images. First we set up a general bayesian frame enabling recovery of discontinuities in such restoration methods. Then we exhibit an extension of the Geman-Reynolds-Charbonnier theorem allowing convenient ltering of SAR images. Due to the high dynamics of radar ERS-1 images, a deterministic algorithm is proposed integrating di erent statistical hypotheses for observation and regularization parts. Besides, we use a well-adapted SAR edge detector instead of the usual gradient in the boundary estimation step of an iterative boundary/intensity restoration algorithm. Intensities are then estimated with a deterministic non-linear method. Finally, the particular behaviour or radar statistics ( law) lead us to de ne a new potential function adapted to speckle regularization while respecting region discontinuities.
1 Introduction Image reconstruction or restoration can often be formulated as minimizing some energy functional: U (f j p) (1) where f = ff g 2 is the image to be reconstructed, p = fp g 2 a given observation and S = fsg the (primary) lattice of sites. It leads generally to non-convex optimization with respect to intensity variables f . This is the more so when a regularization term implicitly accounting for discontinuities is used [3, 4]. In the past few years, several methods have been drawn out in order to overcome the intrinsic complexity of these optimization problems: Graduated Non-Convexity [3] consists in approaching original energy functional (1) by a sequence of convex functions whose minima are by some means guaranteed to converge to the desired optimum. Mean-Field approximation was rst designed for the truncated quadratic also called weak-string model [5, 6]. First, average estimation of the boundary value between neighbor sites is obtained by mean- eld computations. Then, a gradient-type algorithm allows to estimate intensities. s
s
S
s
s
s
S
Anisotropic di usion processes have also been deeply investigated (see [7]). More recently, there has been renewed interest for restoring an intensity
image by simultaneously looking for the explicit line process associated to the primal energy [1, 2, 4]. This approach can be viewed as solving the dual problem of problem (1). It can be viewed as introducing discrete or continuous boundary variables on the dual lattice of the original one [1, 8]. It generally leads to a sequence of convex optimization problems. Moreover, the line process itself can contain some regularization part enabling boundary continuity [9]. Here we shall consider only non-interacting either explicit or implicit line processes. In this paper we follow this last line of algorithms by extending the method originally designed by P. Charbonnier et al. [2]. It was rst aimed to reconstruct images blurred with white gaussian noise. We are interested here in the ltering of synthetic aperture radar (SAR) images which su er from speci c speckle e ects. Many lters have been proposed in the last twenty years. Most of them are based on the statistical properties of speckle images [10, 11]; maximum homogeneous region detection [12] or structure (strong scatterers, edges and lines) detection [13] improve ltering results. We propose in this paper an adaptation of bayesian regularization methods which preserve discontinuities. More speci cally one should possibly process and remove speckle noise while trying to preserve spatial boundaries between regions of interest such as elds or cities, or even to restore thin separating objects such as roads.
2 A Bayesian Model of the Recovery of Boundaries
A bayesian formulation of explicit/implicit boundary processes can be derived as follows. Let us note C = f(r; s)g the set of second order cliques, i.e. of neighbor sites pairs. We also note B the boundary process and b = fb g( )2C any realization of it. Up to now we make no hypotheses on the discrete/continuous nature of this process. The joint posterior probability of f and b writes as: Pr(F = f; B = b j P = p) / Pr(P = p j F = f; B = b) Pr(F = f; B = b) = Pr(P = p j F = f; B = b) Pr(F = f j B = b) Pr(B = b) (2) From now on we make the following assumptions: First, likelihood of data does not depend on boundary realization b, due to a classical site-site independance assumption for data formation. Thus we obtain: Pr(P = p j F = f; B = b) = Pr(P = p j F = f ) which can always be put in an exponential form: Y Pr(P = p j F = f ) = Pr(P = p j F = f ) = Z1 exp f? U1 (p j f )g rs
s
2
s
S
s
s
s
r;s
1
where Z1 is a normalization constant which does not depend on f . For example in the classic additive white gaussian noise assumption we have: 2 Pr(P = p j F = f ) = p 1 exp f? (p 2?2f ) g 2 s
s
s
s
s
s
leading to the quadratic energies: 2 X (p ? f )2 U (p j f ) = (p 2?2f ) and U1 (p j f ) = 22 2 s
s
s
s
s
s
s
s
S
Then, we write the second term of nal member of eq.2 as:
Pr(F = f j B = b) = Z 1(b) exp f? U2 (f j b)g = Z 1(b) exp f? 2
2
X
2C
b U (f ; f )g rs
rs
r
s
(r;s)
(3) meaning that regularization between neighbor sites is accounted for when no boundary exists between them. For example one can assume a gaussian-quadratic expression for U : U (f ; f ) = (f ? f )2 or a more general expression of the type: g( G ), where is a regularization coecient, G = f ? f the local intensity gradient between sites r and s and a suitable gradient normalization coecient. Notice also that normalization constant Z2 does depend here of the conditional line process realization b. Last, third term of eq.2 is the prior distribution for line process. Assuming now boundary variables b to be binary and independant, we write it as: X X (1 ? b )g = Z1 exp f? (b )g Pr(B = b) = Z1 exp f?U3 (b)g = Z1 exp f? 3 3 3 ( )2C ( )2C rs
rs
r
s
r
s
rs
rs
s
r
rs
rs
rs
r;s
r;s
i.e. a stationary binomial model for non-interacting boundaries. Here the value of allows to master the rate of contours in posterior image f . Let us now assume that normalization constant Z2 (b) in equation (3) is nearly independant of b. This consists indeed in neglecting the contour regularity information contained in Z2 (b). Joint MAP estimates of f and b result then from: max Pr(F = f; B = b j P = p) i :e : min f U (f; b j p) = U1 (p j f ) + U2 (f j b) + U3(b) g ( ) ( ) f;b
f;b
The MAP estimate for intensity image itself is thus: f^ = arg min f min U (f; b j p) = U1 (p j f ) + min [ U2 (f j b) + U3(b) ] g f
b
b
which enables to de ne a global e ective regularization energy : U~2 (f ) = min [ U2 (f j b) + U3 (b) ] b
and also locally, due to the independance assumption between b variables: U~ (f ; f ) = min [ b U (f ; f ) + (b ) ] (4) rs
rs
r
s
brs
rs
rs
r
s
rs
In the case where all pixel energies are quadratic, this leads to minimize [3] : X X X (p ? f )2 2 (1 ? b ) (5) b ( f ? f ) + + U (f; b j p) = 2 2 2 ( )2C ( )2C s
s
rs
s
S
r;s
r
rs
s
r;s
This allows estimation of boundary variables once realization f is known: (
2 ^b = 0 if (f ? f ) (there is boundary between r and s) 1 if (f ? f )2 < (no boundary between r and s) r
s
r
s
rs
and to write the e ective interaction potential also called weak-string model [3]: U~ (f ; f ) = min [ (f ? f )2; ]. Using this implicit model is then equivalent to solve problem (5) with explicit boundary variables. Nevertheless smoother, derivable e ective potential functions of the local intensity gradient have been proposed, among them the famous example of Geman-Reynolds-Mc Clure: [1]: rs
r
s
r
s
2 U~ (f ; f ) = ( ( G ) ) with (u) = 1 +u u rs
rs
r
s
Notice that this potential and weak-string model possess the same total variation and curvature at origin characteristics when chosing = and = 2 . We shall assume these relations ful lled in the following. More general requirements for e ective potentials have been explicited in [1, 2]. In the next section we present a new, extended theoretical framework for such potentials functions.
3 Extending the Geman-Reynolds-Charbonnier theorem Here we generalize the Geman-Reynolds-Charbonnier theorem [1, 2] which allows passing from a general implicit to an explicit boundary model for deterministic restoration purposes. Our extension in view of SAR image processing is the following:
Theorem 1 (an extension of Geman-Reynolds-Charbonnier). Let and H two functions de ned on R (H 0 and monotonously increasing) such that +
is strictly concave increasing, and let L and M be de ned as: L = !lim 0(u) and M = lim 0 (u) +1 !0+ u
Then
u
{ there exists a strictly decreasing convex function : [L; M ] 7! [ ; ], such that:
(H(u)) = inf (bH(u) + (b)) with = lim
(v) where (v) = (v) ? v 0 (v) !1 (v) , = lim !0+ L
v
b
M
v
{ for any value of u, the value of b reaching the in mum is b = 0 (H(u)) 0. u
It results from Legendre transform applied to (x) and substituting x = H(u) in it, due to the inversibility of H : R+ 7! R+ . In the classic case where H(u) = u2 we re-obtain the conditions of Geman-Reynolds-Charbonnier Theorem, which states that: ? (u2) = inf bu2 + (b) with = !lim ( (v) ? v 0 (v)) and = lim (v) +1 !0+ L
b
v
M
v
and where the \boundary value" b for which in mum is reached is b = 0(u2 ) (6) The theorem is applied to each local interaction energy considered as a function of local intensity gradient: U~ (f ; f ) = ( ( f ? f )2 ). It allows then to transform a non-convex, implicit model into a convex, explicit one for optimization purposes. For a gaussian observation model one has: 2 X X ( ( f ? f )2 ) U (f j p) = (p 2?2f ) + 2 ( )2C so that total energy can be put in the convex form: 2 X X X U (f; b j p) = (p 2?2f ) + b (f ? f )2 + (b ) (7) 2 ( )2C ( )2C with = , = 2 0 and: b = 0( ( f ? f )2 ) ) 0 (8) In order to optimize an energy function such as (7) we use the ARTUR variant [2] (see g.1). At some step (N ) a new intensity image is computed. Then, boundary variables b are estimated from eq.8, allowing in turn to minimize previous energy by a Newton-Raphson-like algorithm. For this purpose, the zeros of every partial derivative @U (f;@f b j p) are estimated using the following gradient-descent iterative scheme: @U (f; b j p) ( ) 2 @ U (f; b j p) ( ) 8s 2 S @f where ! = f ( +1) ? f ( ) = ? ! @f 2 (9) pixels being updated in a sequential mode. After some iterations of this kind a new intensity image f ( +1) is obtained, so that boundary values b can in turn be re-estimated from eq.8 and so on. Convergence is attained when both boundary and intensity variables are stable [2]. However, the gaussian observation model as well as the truncated gaussian regularization one are not always adapted to the processing of SAR images. Our general aim is to apply extended Theorem 1 when attachment to data functions are no more quadratic and to regularization models which do not necessarily depend on the local intensity gradient. u
u
r
rs
s
s
s
r
r
s
S
s
r;s
s
rs
s
s
s
S
r
s
rs
r;s
r
r;s
s
rs
rs
s
n
n
n
s
s
n
s
s
N
rs
Initialization (0)
f
=p
Step (N) Chose boundary variables using eq.8
step (1), (2) ..(n)..
Optimize pixel image using eq.9
Stop test Final result
Fig.1. Modi ed Artur Algorithm
4 Statistics of SAR images We shortly describe the well-known statistics of radar images. Denoting a exp(i ) the complex eld received by the radar p antenna, the intensity is de ned by I = jaj2 and the amplitude by A = I = jaj. L-looks images are obtained by averaging L independent images (either by dividing the available bandwith of the SAR system in L parts, or by spatially averaging L pixels [14]). The ERS-1 images, on which results are presented in this paper, are PRI Products (Precision Images) i:e: 3-look amplitude images with a resolution of 12.5 meters. Under some assumptions implying a rough surface compared to the radar wavelength (5.6cm for ERS-1) statistics of an area with constant re ectivity are given by [15]: 2 2 ?1 Pr(A j ) = ?2L(L) A 2 exp (? LA2 ) (10) denoting by ? the Gamma function [16] and by the square root of the mean re ected intensity of the area. This probability density function is called in the following and implies the well-known speckle phenomenon which a ects radar images. Other probability density functions could be used to take into account textured regions [17] [18]. L
L
L
5 Edge detection on SAR images The use of the radiometric di erence G is ill-adapted to SAR image statistics [20]. On radar images, such a detector gives di erent detection results depending on the radiometric mean of a region since the standard deviation increases proportionally to the mean. An edge detector widely used in radar imagery is based on the ratio between radiometric means [19] [20]. Denoting Re(r) and Re(s) two adjacent regions associated to each site of a clique (r; s) (horizontal, vertical or diagonal) as shown on gure 2, the edge detector response d is obtained by: rs
rs
d = 1 ? min A ; A A A
r
s
s
r
rs
with A the radiometric mean computed on the neighborhood Re(r) around r. To preserve thin structures, the neighborhood size can be chosen of 5 2 pixels, although this implies more false alarms than with larger regions. This detector has proved to have a good behaviour on SAR images [20]. Let us also de ne a ponctual edge detector: r
d = 1 ? min AA ; AA
r
s
s
r
rs
site r
diagonal clique site s vertical clique
horizontal clique pixels of region Re(s) pixels of region Re(r)
Fig.2. Used regions to compute drs In the following, we propose to use these results to restore the original scene image while preserving the discontinuities.
6 Methodology of the study We propose in this paper to do some experiments to analyse the behaviour of restoration techniques with recovery of discontinuities on SAR images. We test several con gurations and analyse the obtained results, particularly the restoration e ects on bright and dark regions. The following choices can be made:
{ Choice of the probability density function:
We have tested two possibilities, gaussian and distribution. The second choice should be more adapted to SAR images as claimed in section 4.
{ Choice of the regularization function H:
The usual choice is a quadratic function. We will show in this paper that other functions can be used in the case of SAR images.
{ Choice of the boundary process de nition: The right expression is b = 0(H( G )). Using the remarks made on edge rs
rs
detection in SAR images, we propose to test the boundary process de ned by b = 0(H( d )) where is an appropriate control parameter of d . Of course, in the case of this modi cation, we are not any more in the theoretical framework de ned by the generalized Geman-Reynolds-Charbonnier theorem, and all the convergence prooves are not valid any more. Therefore, in this case, the results proposed are surely not the MAP estimates. Let us also note that the a new choice of the regularization function H provides a new boundary process. Besides, we propose to introduce a new way to estimate the boundary process, by using average values instead of amplitudes in the gradient or ratio. The boundary process which is estimated using average values is noted b (or b respectively). Once again, in this case the theoretical framework is not veri ed any more. The following table summarizes the di erent tests we have made (with G = f ? f and R = ff ). rs
rs
rs
rs
rs
rs
r
r
s
rs
s
distribution gaussian reg. function x = Grs , H(x) = x2 x = Grs , H(x) = x2 x = Rrs , H(x) = x + x1 ? 2 boundary brs brs brs brs brs section 7.1 7.2 8.1.1 8.1.2 8.2
Table 1. Summary of the tests and the corresponding sections. For all the tested con gurations, estimations of the boundary process with ponctual and average values are tried. As for the function, the following func-
tion is used in all tests: (x) = 1 +x x . For each con guration, the regularized result of a synthetic 3-look amplitude image ( g.3a) is given. Besides, the remarks made on the regularization process are also veri ed on the boundary process image which shows where edges are detected. Then a comparison of all the obtained results is given in section 9. The best method is applied on a real (ERS-1) SAR image ( g.5a) showing well de ned elds and roads.
7 Gaussian probability and quadratic regularization function Recall that in the case of a gaussian probability and a quadratic H function, the energy U (f; b j p) is expressed by:
U (f; b j p) =
X
2
s
(p ? f )2 + s
X
s
S
2C
b (f ? f )2 + rs
s
X
r
c
With (x) = 1 +x x [1], the boundary process is given by:
2C
(b ) rs
c
1 b = 0 (H( G )) = 2 G 1 + ( )2 rs
rs
rs
Since the rst and second derivatives are: @U = 2(f ? p ) + 2 X b (f ? f ) and @ 2 U = 2(1 + X b ) @f @f 2 2Ns 2Ns the iterative scheme is given by: s
s
rs
s
f( s
n+1)
s
r
rs
s
r
? f( ) = ? n
s
1+
1X r
2Ns
b
f( ) ? p + n
s
rs
X
s
r
2Ns
r
!
b (f ( ) ? f ( ) ) rs
n
s
n
r
(11)
The regularization method has been applied on a synthetic 3-look amplitude image ( g.3a), in both cases: right b and modi ed b (which is equivalent to replace the gradient edge detector by a ratio edge detector). For each case, estimation with the amplitude values and local average values have been tested. rs
rs
7.1 Using brs The results obtained are not very statisfying. A di erent behaviour in dark and bright regions is observed: the dark areas are strongly regularized, while the bright areas show many variations. Note that these e ects, although present in both cases (for b or b ), are reduced when the average values are used to estimate b ( g.3b). Besides, the result contains many isolated bright points in the case of the ponctual b estimator. rs
rs
rs
rs
7.2 Using brs
Since the boundary process is now well adapted to the detection of discontinuities in SAR images, the previous e ects have disapeared. Resulting images are satisfying with a good edge preservation and a strong ltering of homogeneous regions ( g.3c, see also section 9). Let us also note that the homogeneous areas are much more regularized in the case of the average estimation b , than in ponctual estimation. These satisfying results seem at rst in contradiction with speckle behaviour since a SAR image does not follow a gaussian pdf. An explanation of this result is the behaviour of the ltering process. Indeed, equation (11) can be written as: ! X 1 f ( +1) = p + b f( ) X b 1+ 2Ns rs
n
s
s
r
rs
rs
2Ns
n
r
r
showing that the ltered value is a weighted mean of the observation and the neighboring values, with weights depending of the edge presence. Therefore, after a few iterations, the probability density function turns out to be almost gaussian.
8 probability
8.1 Quadratic regularization function
To improve the restoration method a distribution, more adapted to SAR images (x4), has been used [21]: 2 2 ?1 Pr(p jf ) = k pf 2 exp(? Lp2 ) f k being a normalization constant depending on the number of looks used to create the SAR image (eq.10). Therefore, the following data attachment energy ! is deduced: 2 X p U1 (pjf ) = L 2 log f + f 2 s
s
s
L
s
s
L
s
2
s
s
s
s
S
providing the expression of sequence f ? f ( ): ! ?1 2L (f ( )2 ? p2) + 2 X b (f ( ) ? f ( ) ) 2L (3p2 ? f ( )2 ) + 2 X b f ( )3 2Ns f ( )4 2Ns Many problems arise using this expression. First, in order to assure the convergence of the Newton-Raphson method, coecient ! in eq.9 has to be positive, which implies: 2L (3p2 ? f ( )2 ) + 2 X b > 0 f ( )4 2Ns In practice, starting with the original image p for f initialization, the condition 3p2 > f 2 is always veri ed and no convergence problem accurs. (n+1)
s
s
n
n
n
n
s
s
n
rs
s
s
s
r
r
n
n
s
s
rs
s
s
s
rs
s
r
n
s
n
r
Using brs The previous phenomenon (stronger ltering of the dark areas) is now reversed, which is of course still not satisfying. In the case of a ponctual boundary process estimation, there are many isolated dark points in the result, which is not the case with the average estimation b . rs
Using brs The same phenomenon of completely di erent ltering behaviours
according to the local radiometric average is observed and the same remark about isolated dark points holds. Of course, the aim of the regularization process is to smooth every homogeneous region of the image in the same way, whatever its radiometric mean. But it is not the case using a distribution and a quadratic regularization function. Let us consider two extreme cases: jf ? p j ) ; 1 K 10) : when f ! 1 (precisely for f 2 >> KL max(1 ; jf ? f j s
s
X
s
n+1)
s 2NX r
r
s
r
b f( rs
f(
s
s
2Ns
r
n)
b
rs
therefore, only the regularization term counts on bright regions, leading to very smoothed areas; jf ? p j when f ! 0 (precisely for f 2
Abstract. In this paper, we study the behaviour of contour recovery
when ltering radar images. We start from recent methods lying on an equivalence scheme between implicit and explicit boundary processes in image restoration [1, 2]. Here we extend them to the processing of synthetic aperture radar (SAR) images. First we set up a general bayesian frame enabling recovery of discontinuities in such restoration methods. Then we exhibit an extension of the Geman-Reynolds-Charbonnier theorem allowing convenient ltering of SAR images. Due to the high dynamics of radar ERS-1 images, a deterministic algorithm is proposed integrating di erent statistical hypotheses for observation and regularization parts. Besides, we use a well-adapted SAR edge detector instead of the usual gradient in the boundary estimation step of an iterative boundary/intensity restoration algorithm. Intensities are then estimated with a deterministic non-linear method. Finally, the particular behaviour or radar statistics ( law) lead us to de ne a new potential function adapted to speckle regularization while respecting region discontinuities.
1 Introduction Image reconstruction or restoration can often be formulated as minimizing some energy functional: U (f j p) (1) where f = ff g 2 is the image to be reconstructed, p = fp g 2 a given observation and S = fsg the (primary) lattice of sites. It leads generally to non-convex optimization with respect to intensity variables f . This is the more so when a regularization term implicitly accounting for discontinuities is used [3, 4]. In the past few years, several methods have been drawn out in order to overcome the intrinsic complexity of these optimization problems: Graduated Non-Convexity [3] consists in approaching original energy functional (1) by a sequence of convex functions whose minima are by some means guaranteed to converge to the desired optimum. Mean-Field approximation was rst designed for the truncated quadratic also called weak-string model [5, 6]. First, average estimation of the boundary value between neighbor sites is obtained by mean- eld computations. Then, a gradient-type algorithm allows to estimate intensities. s
s
S
s
s
s
S
Anisotropic di usion processes have also been deeply investigated (see [7]). More recently, there has been renewed interest for restoring an intensity
image by simultaneously looking for the explicit line process associated to the primal energy [1, 2, 4]. This approach can be viewed as solving the dual problem of problem (1). It can be viewed as introducing discrete or continuous boundary variables on the dual lattice of the original one [1, 8]. It generally leads to a sequence of convex optimization problems. Moreover, the line process itself can contain some regularization part enabling boundary continuity [9]. Here we shall consider only non-interacting either explicit or implicit line processes. In this paper we follow this last line of algorithms by extending the method originally designed by P. Charbonnier et al. [2]. It was rst aimed to reconstruct images blurred with white gaussian noise. We are interested here in the ltering of synthetic aperture radar (SAR) images which su er from speci c speckle e ects. Many lters have been proposed in the last twenty years. Most of them are based on the statistical properties of speckle images [10, 11]; maximum homogeneous region detection [12] or structure (strong scatterers, edges and lines) detection [13] improve ltering results. We propose in this paper an adaptation of bayesian regularization methods which preserve discontinuities. More speci cally one should possibly process and remove speckle noise while trying to preserve spatial boundaries between regions of interest such as elds or cities, or even to restore thin separating objects such as roads.
2 A Bayesian Model of the Recovery of Boundaries
A bayesian formulation of explicit/implicit boundary processes can be derived as follows. Let us note C = f(r; s)g the set of second order cliques, i.e. of neighbor sites pairs. We also note B the boundary process and b = fb g( )2C any realization of it. Up to now we make no hypotheses on the discrete/continuous nature of this process. The joint posterior probability of f and b writes as: Pr(F = f; B = b j P = p) / Pr(P = p j F = f; B = b) Pr(F = f; B = b) = Pr(P = p j F = f; B = b) Pr(F = f j B = b) Pr(B = b) (2) From now on we make the following assumptions: First, likelihood of data does not depend on boundary realization b, due to a classical site-site independance assumption for data formation. Thus we obtain: Pr(P = p j F = f; B = b) = Pr(P = p j F = f ) which can always be put in an exponential form: Y Pr(P = p j F = f ) = Pr(P = p j F = f ) = Z1 exp f? U1 (p j f )g rs
s
2
s
S
s
s
s
r;s
1
where Z1 is a normalization constant which does not depend on f . For example in the classic additive white gaussian noise assumption we have: 2 Pr(P = p j F = f ) = p 1 exp f? (p 2?2f ) g 2 s
s
s
s
s
s
leading to the quadratic energies: 2 X (p ? f )2 U (p j f ) = (p 2?2f ) and U1 (p j f ) = 22 2 s
s
s
s
s
s
s
s
S
Then, we write the second term of nal member of eq.2 as:
Pr(F = f j B = b) = Z 1(b) exp f? U2 (f j b)g = Z 1(b) exp f? 2
2
X
2C
b U (f ; f )g rs
rs
r
s
(r;s)
(3) meaning that regularization between neighbor sites is accounted for when no boundary exists between them. For example one can assume a gaussian-quadratic expression for U : U (f ; f ) = (f ? f )2 or a more general expression of the type: g( G ), where is a regularization coecient, G = f ? f the local intensity gradient between sites r and s and a suitable gradient normalization coecient. Notice also that normalization constant Z2 does depend here of the conditional line process realization b. Last, third term of eq.2 is the prior distribution for line process. Assuming now boundary variables b to be binary and independant, we write it as: X X (1 ? b )g = Z1 exp f? (b )g Pr(B = b) = Z1 exp f?U3 (b)g = Z1 exp f? 3 3 3 ( )2C ( )2C rs
rs
r
s
r
s
rs
rs
s
r
rs
rs
rs
r;s
r;s
i.e. a stationary binomial model for non-interacting boundaries. Here the value of allows to master the rate of contours in posterior image f . Let us now assume that normalization constant Z2 (b) in equation (3) is nearly independant of b. This consists indeed in neglecting the contour regularity information contained in Z2 (b). Joint MAP estimates of f and b result then from: max Pr(F = f; B = b j P = p) i :e : min f U (f; b j p) = U1 (p j f ) + U2 (f j b) + U3(b) g ( ) ( ) f;b
f;b
The MAP estimate for intensity image itself is thus: f^ = arg min f min U (f; b j p) = U1 (p j f ) + min [ U2 (f j b) + U3(b) ] g f
b
b
which enables to de ne a global e ective regularization energy : U~2 (f ) = min [ U2 (f j b) + U3 (b) ] b
and also locally, due to the independance assumption between b variables: U~ (f ; f ) = min [ b U (f ; f ) + (b ) ] (4) rs
rs
r
s
brs
rs
rs
r
s
rs
In the case where all pixel energies are quadratic, this leads to minimize [3] : X X X (p ? f )2 2 (1 ? b ) (5) b ( f ? f ) + + U (f; b j p) = 2 2 2 ( )2C ( )2C s
s
rs
s
S
r;s
r
rs
s
r;s
This allows estimation of boundary variables once realization f is known: (
2 ^b = 0 if (f ? f ) (there is boundary between r and s) 1 if (f ? f )2 < (no boundary between r and s) r
s
r
s
rs
and to write the e ective interaction potential also called weak-string model [3]: U~ (f ; f ) = min [ (f ? f )2; ]. Using this implicit model is then equivalent to solve problem (5) with explicit boundary variables. Nevertheless smoother, derivable e ective potential functions of the local intensity gradient have been proposed, among them the famous example of Geman-Reynolds-Mc Clure: [1]: rs
r
s
r
s
2 U~ (f ; f ) = ( ( G ) ) with (u) = 1 +u u rs
rs
r
s
Notice that this potential and weak-string model possess the same total variation and curvature at origin characteristics when chosing = and = 2 . We shall assume these relations ful lled in the following. More general requirements for e ective potentials have been explicited in [1, 2]. In the next section we present a new, extended theoretical framework for such potentials functions.
3 Extending the Geman-Reynolds-Charbonnier theorem Here we generalize the Geman-Reynolds-Charbonnier theorem [1, 2] which allows passing from a general implicit to an explicit boundary model for deterministic restoration purposes. Our extension in view of SAR image processing is the following:
Theorem 1 (an extension of Geman-Reynolds-Charbonnier). Let and H two functions de ned on R (H 0 and monotonously increasing) such that +
is strictly concave increasing, and let L and M be de ned as: L = !lim 0(u) and M = lim 0 (u) +1 !0+ u
Then
u
{ there exists a strictly decreasing convex function : [L; M ] 7! [ ; ], such that:
(H(u)) = inf (bH(u) + (b)) with = lim
(v) where (v) = (v) ? v 0 (v) !1 (v) , = lim !0+ L
v
b
M
v
{ for any value of u, the value of b reaching the in mum is b = 0 (H(u)) 0. u
It results from Legendre transform applied to (x) and substituting x = H(u) in it, due to the inversibility of H : R+ 7! R+ . In the classic case where H(u) = u2 we re-obtain the conditions of Geman-Reynolds-Charbonnier Theorem, which states that: ? (u2) = inf bu2 + (b) with = !lim ( (v) ? v 0 (v)) and = lim (v) +1 !0+ L
b
v
M
v
and where the \boundary value" b for which in mum is reached is b = 0(u2 ) (6) The theorem is applied to each local interaction energy considered as a function of local intensity gradient: U~ (f ; f ) = ( ( f ? f )2 ). It allows then to transform a non-convex, implicit model into a convex, explicit one for optimization purposes. For a gaussian observation model one has: 2 X X ( ( f ? f )2 ) U (f j p) = (p 2?2f ) + 2 ( )2C so that total energy can be put in the convex form: 2 X X X U (f; b j p) = (p 2?2f ) + b (f ? f )2 + (b ) (7) 2 ( )2C ( )2C with = , = 2 0 and: b = 0( ( f ? f )2 ) ) 0 (8) In order to optimize an energy function such as (7) we use the ARTUR variant [2] (see g.1). At some step (N ) a new intensity image is computed. Then, boundary variables b are estimated from eq.8, allowing in turn to minimize previous energy by a Newton-Raphson-like algorithm. For this purpose, the zeros of every partial derivative @U (f;@f b j p) are estimated using the following gradient-descent iterative scheme: @U (f; b j p) ( ) 2 @ U (f; b j p) ( ) 8s 2 S @f where ! = f ( +1) ? f ( ) = ? ! @f 2 (9) pixels being updated in a sequential mode. After some iterations of this kind a new intensity image f ( +1) is obtained, so that boundary values b can in turn be re-estimated from eq.8 and so on. Convergence is attained when both boundary and intensity variables are stable [2]. However, the gaussian observation model as well as the truncated gaussian regularization one are not always adapted to the processing of SAR images. Our general aim is to apply extended Theorem 1 when attachment to data functions are no more quadratic and to regularization models which do not necessarily depend on the local intensity gradient. u
u
r
rs
s
s
s
r
r
s
S
s
r;s
s
rs
s
s
s
S
r
s
rs
r;s
r
r;s
s
rs
rs
s
n
n
n
s
s
n
s
s
N
rs
Initialization (0)
f
=p
Step (N) Chose boundary variables using eq.8
step (1), (2) ..(n)..
Optimize pixel image using eq.9
Stop test Final result
Fig.1. Modi ed Artur Algorithm
4 Statistics of SAR images We shortly describe the well-known statistics of radar images. Denoting a exp(i ) the complex eld received by the radar p antenna, the intensity is de ned by I = jaj2 and the amplitude by A = I = jaj. L-looks images are obtained by averaging L independent images (either by dividing the available bandwith of the SAR system in L parts, or by spatially averaging L pixels [14]). The ERS-1 images, on which results are presented in this paper, are PRI Products (Precision Images) i:e: 3-look amplitude images with a resolution of 12.5 meters. Under some assumptions implying a rough surface compared to the radar wavelength (5.6cm for ERS-1) statistics of an area with constant re ectivity are given by [15]: 2 2 ?1 Pr(A j ) = ?2L(L) A 2 exp (? LA2 ) (10) denoting by ? the Gamma function [16] and by the square root of the mean re ected intensity of the area. This probability density function is called in the following and implies the well-known speckle phenomenon which a ects radar images. Other probability density functions could be used to take into account textured regions [17] [18]. L
L
L
5 Edge detection on SAR images The use of the radiometric di erence G is ill-adapted to SAR image statistics [20]. On radar images, such a detector gives di erent detection results depending on the radiometric mean of a region since the standard deviation increases proportionally to the mean. An edge detector widely used in radar imagery is based on the ratio between radiometric means [19] [20]. Denoting Re(r) and Re(s) two adjacent regions associated to each site of a clique (r; s) (horizontal, vertical or diagonal) as shown on gure 2, the edge detector response d is obtained by: rs
rs
d = 1 ? min A ; A A A
r
s
s
r
rs
with A the radiometric mean computed on the neighborhood Re(r) around r. To preserve thin structures, the neighborhood size can be chosen of 5 2 pixels, although this implies more false alarms than with larger regions. This detector has proved to have a good behaviour on SAR images [20]. Let us also de ne a ponctual edge detector: r
d = 1 ? min AA ; AA
r
s
s
r
rs
site r
diagonal clique site s vertical clique
horizontal clique pixels of region Re(s) pixels of region Re(r)
Fig.2. Used regions to compute drs In the following, we propose to use these results to restore the original scene image while preserving the discontinuities.
6 Methodology of the study We propose in this paper to do some experiments to analyse the behaviour of restoration techniques with recovery of discontinuities on SAR images. We test several con gurations and analyse the obtained results, particularly the restoration e ects on bright and dark regions. The following choices can be made:
{ Choice of the probability density function:
We have tested two possibilities, gaussian and distribution. The second choice should be more adapted to SAR images as claimed in section 4.
{ Choice of the regularization function H:
The usual choice is a quadratic function. We will show in this paper that other functions can be used in the case of SAR images.
{ Choice of the boundary process de nition: The right expression is b = 0(H( G )). Using the remarks made on edge rs
rs
detection in SAR images, we propose to test the boundary process de ned by b = 0(H( d )) where is an appropriate control parameter of d . Of course, in the case of this modi cation, we are not any more in the theoretical framework de ned by the generalized Geman-Reynolds-Charbonnier theorem, and all the convergence prooves are not valid any more. Therefore, in this case, the results proposed are surely not the MAP estimates. Let us also note that the a new choice of the regularization function H provides a new boundary process. Besides, we propose to introduce a new way to estimate the boundary process, by using average values instead of amplitudes in the gradient or ratio. The boundary process which is estimated using average values is noted b (or b respectively). Once again, in this case the theoretical framework is not veri ed any more. The following table summarizes the di erent tests we have made (with G = f ? f and R = ff ). rs
rs
rs
rs
rs
rs
r
r
s
rs
s
distribution gaussian reg. function x = Grs , H(x) = x2 x = Grs , H(x) = x2 x = Rrs , H(x) = x + x1 ? 2 boundary brs brs brs brs brs section 7.1 7.2 8.1.1 8.1.2 8.2
Table 1. Summary of the tests and the corresponding sections. For all the tested con gurations, estimations of the boundary process with ponctual and average values are tried. As for the function, the following func-
tion is used in all tests: (x) = 1 +x x . For each con guration, the regularized result of a synthetic 3-look amplitude image ( g.3a) is given. Besides, the remarks made on the regularization process are also veri ed on the boundary process image which shows where edges are detected. Then a comparison of all the obtained results is given in section 9. The best method is applied on a real (ERS-1) SAR image ( g.5a) showing well de ned elds and roads.
7 Gaussian probability and quadratic regularization function Recall that in the case of a gaussian probability and a quadratic H function, the energy U (f; b j p) is expressed by:
U (f; b j p) =
X
2
s
(p ? f )2 + s
X
s
S
2C
b (f ? f )2 + rs
s
X
r
c
With (x) = 1 +x x [1], the boundary process is given by:
2C
(b ) rs
c
1 b = 0 (H( G )) = 2 G 1 + ( )2 rs
rs
rs
Since the rst and second derivatives are: @U = 2(f ? p ) + 2 X b (f ? f ) and @ 2 U = 2(1 + X b ) @f @f 2 2Ns 2Ns the iterative scheme is given by: s
s
rs
s
f( s
n+1)
s
r
rs
s
r
? f( ) = ? n
s
1+
1X r
2Ns
b
f( ) ? p + n
s
rs
X
s
r
2Ns
r
!
b (f ( ) ? f ( ) ) rs
n
s
n
r
(11)
The regularization method has been applied on a synthetic 3-look amplitude image ( g.3a), in both cases: right b and modi ed b (which is equivalent to replace the gradient edge detector by a ratio edge detector). For each case, estimation with the amplitude values and local average values have been tested. rs
rs
7.1 Using brs The results obtained are not very statisfying. A di erent behaviour in dark and bright regions is observed: the dark areas are strongly regularized, while the bright areas show many variations. Note that these e ects, although present in both cases (for b or b ), are reduced when the average values are used to estimate b ( g.3b). Besides, the result contains many isolated bright points in the case of the ponctual b estimator. rs
rs
rs
rs
7.2 Using brs
Since the boundary process is now well adapted to the detection of discontinuities in SAR images, the previous e ects have disapeared. Resulting images are satisfying with a good edge preservation and a strong ltering of homogeneous regions ( g.3c, see also section 9). Let us also note that the homogeneous areas are much more regularized in the case of the average estimation b , than in ponctual estimation. These satisfying results seem at rst in contradiction with speckle behaviour since a SAR image does not follow a gaussian pdf. An explanation of this result is the behaviour of the ltering process. Indeed, equation (11) can be written as: ! X 1 f ( +1) = p + b f( ) X b 1+ 2Ns rs
n
s
s
r
rs
rs
2Ns
n
r
r
showing that the ltered value is a weighted mean of the observation and the neighboring values, with weights depending of the edge presence. Therefore, after a few iterations, the probability density function turns out to be almost gaussian.
8 probability
8.1 Quadratic regularization function
To improve the restoration method a distribution, more adapted to SAR images (x4), has been used [21]: 2 2 ?1 Pr(p jf ) = k pf 2 exp(? Lp2 ) f k being a normalization constant depending on the number of looks used to create the SAR image (eq.10). Therefore, the following data attachment energy ! is deduced: 2 X p U1 (pjf ) = L 2 log f + f 2 s
s
s
L
s
s
L
s
2
s
s
s
s
S
providing the expression of sequence f ? f ( ): ! ?1 2L (f ( )2 ? p2) + 2 X b (f ( ) ? f ( ) ) 2L (3p2 ? f ( )2 ) + 2 X b f ( )3 2Ns f ( )4 2Ns Many problems arise using this expression. First, in order to assure the convergence of the Newton-Raphson method, coecient ! in eq.9 has to be positive, which implies: 2L (3p2 ? f ( )2 ) + 2 X b > 0 f ( )4 2Ns In practice, starting with the original image p for f initialization, the condition 3p2 > f 2 is always veri ed and no convergence problem accurs. (n+1)
s
s
n
n
n
n
s
s
n
rs
s
s
s
r
r
n
n
s
s
rs
s
s
s
rs
s
r
n
s
n
r
Using brs The previous phenomenon (stronger ltering of the dark areas) is now reversed, which is of course still not satisfying. In the case of a ponctual boundary process estimation, there are many isolated dark points in the result, which is not the case with the average estimation b . rs
Using brs The same phenomenon of completely di erent ltering behaviours
according to the local radiometric average is observed and the same remark about isolated dark points holds. Of course, the aim of the regularization process is to smooth every homogeneous region of the image in the same way, whatever its radiometric mean. But it is not the case using a distribution and a quadratic regularization function. Let us consider two extreme cases: jf ? p j ) ; 1 K 10) : when f ! 1 (precisely for f 2 >> KL max(1 ; jf ? f j s
s
X
s
n+1)
s 2NX r
r
s
r
b f( rs
f(
s
s
2Ns
r
n)
b
rs
therefore, only the regularization term counts on bright regions, leading to very smoothed areas; jf ? p j when f ! 0 (precisely for f 2
