Markovian models for one dimensional structure estimation on heavily noisy imagery.
Abstract
Radar (SAR) images often exhibit profound appearance variations due to a variety of factors including clutter noise produced by the coherent nature of the illumination. Ultrasound images and infrared images have similar cluttered appearance, that make 1 dimensional structures, as edges and object boundaries difficult to locate. Structure information is usually extracted in two steps: firstly, building and edge strength mask classifying pixels as edge points by hypothesis testing, and secondly computing on pixel wide connected edges. With constant false alarm rate (CFAR) edge strength detectors for speckle clutter, the image needs to be scanned by a sliding window composed of several differently oriented splitting subwindows. The accuracy of edge location for these ratio detectors depends strongly on the orientation of the subwindows. In this work we propose to transform the edge strength detection problem into a binary segmentation problem in the undecimated wavelet domain, solvable using parallel 1d Hidden Markov Models. For general dependency models, exact estimation of the state map becomes computationally complex, but in our model, exact MAP is feasible. The effectiveness of our approach is demonstrated on simulated noisy reallife natural images with available ground truth, while the strength of our output edge map is measured with Pratt’s, Baddeley an Kappa proficiency measures. Finally, analysis and experiments on three different types of SAR images, with different polarizations, resolutions and textures, illustrate that the proposed method can detect structure on SAR images effectively, providing a very good start point for active contour methods.
Keywords:Haar wavelet edge detection; Graphical models, SAR images; undecimated wavelet transform
1 Introduction
Satellite born Synthetic Aperture Radar (SAR) images represent the complex reflectivity map of a scene, and unlike optical images, the interpretation is not consistent with common visual perception. Furthermore, the direct application of conventional image processing tools conceived from an optical point of view usually gives suboptimum results on SAR data. Among the most important differences between optical and SAR images are the appropriate noise models. In special, a complex SAR image may be represented as the convolution of the local complex reflectivity of the observed area with the impulse response of the SAR system. Hence, the SAR data pixels are roughly the low pass filtered version of the complex local scattering properties of the observed scene. A usual model for the complex image is obtained considering its real and imaginary parts as independent, zero mean and equally distributed random variables. Consequently, the intensity of a SAR image follows a negative exponential distribution, and can be further modeled as a product of the radar cross section and a noise term called speckle, which is exponentially distributed with mean equal to one, Tello Alonso et al. (2011). In the framework of SAR processing, it is useful to take the logarithm of the original signal in order to manage the multiplicative speckle; thus, after this transformation, the speckle is not only additive and signal independent, but also approximately Gaussian. Moreover, the logarithm operation is helpful in reducing the large dynamic range of SAR data.
A crucial point for SAR image automatic interpretation is the low level step of scene segmentation, i. e. the decomposition of the image in a tessellation of uniform areas. Edgebased segmentation schemes aim to find out the transitions between uniform areas, rather than directly identifying them. The related algorithms generally work in two stages: they firstly compute an edge strength map of the scene and finally extract the local maxima of this map. The first stage is usually achieved using an statistical edge detection filter. It consists in scanning the image with a sliding tworegion window and evaluating for each position if there is a change between the two regions by hypothesis testing. The test are dependant of the true underlying conditional distributions that are emitted by the states (classes) of the segmentation field.
Our algorithm for edge strength mask detection in heavily noisy images, like logamplitude satellite born Synthetic Aperture Radar (SAR) images, which is proposed in this paper, relies on the difference of statistical behavior along the wavelet scales of the Gaussian noise in front of the edges. The discontinuities are highlighted by the wavelet transform, and they tend to persist over scales. These two key properties, very well known in the optical image processing literature, were described by Romberg et al. (2001) as:

Two Populations: Smooth signal/image regions are represented by small wavelet coefficients, while edges, ridges, and other singularities are represented by large coefficients.

Persistence: Large and small wavelet coefficient values cascade along the scale/space representation.
Interpreting these properties statistically, Two Populations imply that wavelet coefficients have nonGaussian marginal statistics, with a large peak at zero (many small coefficients) and heavy tails (a few large coefficients), while Persistence implies that wavelet coefficient values are statistically dependent along the scale/space representation. Nevertheless, the critical subsampling in the Discrete Wavelet transform (DWT), which discards every second wavelet coefficient at each decomposition level, results in wavelet coefficients that are highly dependent on their location in the subsampling lattice. Small shifts in the input waveform may cause large changes in the wavelet coefficients, large variations in the distribution of energy at different scales, and possibly large changes in reconstructed waveforms.
In Tello Alonso et al. (2011), the authors advocate for the use of the undecimated wavelet transform for detecting and enhancing edges on SAR images. They propose a deterministic multiscale thresholding scheme to produce an edge strength mask that single out edge candidates to apply snake algorithms in order to get one pixel wide continuous edges. This combination algorithms, in the radiometric domain, are popular in the literature of SAR edge detection. Radiometric edge strength masks are usually obtained by hypothesis testing, using the likelihood ratio filter, the robust rank order or the Wilcoxon test, or costume made maximum likelihood ratio tests (see Lim et al (2006), Germain et al (2001), Giron et al (2012) and references therein).
Our approach is more automatic than the one introduced by Tello Alonso et al (2011), since we estimate all parameters from the data using a Bayesian model in the undecimated wavelet domain. We model each string of pixel related wavelet coefficients as independent 1d Hidden Markov Model with Laplacian and Gaussian marginal distributions, related to “edge” and “noedge” states, respectively, as in Wang et al (2009). We derive our own MAP equations to compute the model parameters in order to predict intrascale and interscale edge strength maps. Then we combine these maps into one edge strength map that can be processed with standard Gradient based edge detection algorithms, watershed algorithms or active contour algorithms.
In Sections 2 and 3 we introduce wavelet transforms and the Bayesian modeling of wavelet coefficients, leaving the proof of the theorems to Appendix A. In Section 4 we discuss processing of the obtained esge strenght mask, and develop some synthetic experiments, adding noise to good quality optical images with available ground truth. Edge quality is measured with three different supervised edge measures, Pratt’s FOM, Baddely and Kappa. We also provide experiments with benchmark SAR images, that compare with other results from the literature. The conclusions and prospects are left for for Section 5.
2 Undecimated Haar wavelet transform
The undecimated discrete wavelet transform (UDWT) has been independently discovered several times for different purposes and under different names, shift/translation invariant wavelet transform, stationary wavelet transform or redundant wavelet transform. The key point is that it is redundant, shift invariant, and gives an approximation to the continuous wavelet transform denser than the approximation given by the discrete wavelet transform (DWT), Guo et al. (1995). The most cited implementations of this idea are Mallat’s à trous algorithm and Beylkin’s algorithm; for details, see Starck et al. (2007) and references therein.
In our implementation, we use an algorithm based on the idea of no decimation. The size of the coefficients array do not diminish from level to level, so by using all coefficients at each level, we get very well allocated highfrequency information. In the UDWT algorithm, the number of pixels involved in computing a given coefficient grows slower and so the relation between the frequency and spatial information is more precise, Gyaourova et al. (2001).
The undecimated wavelet transform can be defined with all types of wavelet filters, but we use the Haar wavelet since its analyzing filter (the shortest possible) gives good localization and poor suppression of discontinuities. Tello Alonso et al. (2011) consider Beylkin’s stationary Haar wavelet transform for edge detection on SAR data. Thresholding Beylkin’s algorithm produces very smooth results, which makes the images to look a little blurred. On the other hand, it makes the regions of constant intensity look almost noise clean. In the case of SAR imagery, we preferred to denoise first with a Wiener filter, applying the undecimated wavelet transform next, with the Haar filter. In Figure 1 we can see plots of 3 level Haar UDWT coefficients of a natural image taken form the Berkeley Database, Martin et al. (2001). The undecimated wavelet transform allows to link the pixel with the coefficient, beside giving better edge localization than the DWT.
3 Wavelet Reduced Bayesian Net Model
The 2d undecimated wavelet transform consists of three bands of detail coefficients and one approximation band per level (see Figure 1), each of which has the same size as the original image. To statistically match the twopopulation property of the wavelet transform, we consider the collection of small/large wavelet coefficients from an edge detection perspective. Romberg et al. (2001), in a denoising context, modeled the wavelet coefficients as outcomes emitted by a Gaussian mixture that approximates the marginal pdf. In our case, “small” and “large” mean differently. We continue to expect most of the nodes to be “small”, since smooth regions produce chains of small coefficients, while singularities are chains of “large” coefficients. But we also expect large coefficients in close regions near the edges, so Gaussian distributions will give poor localization, making the algorithm return thick edges instead of thin ones. A more heavytailed distribution, like the Laplacian, will be better suited for modeling the state “edge”, while the Gaussian will be kept as a marginal for the “noedge” state. Thus, the conditional densities of the coefficients will be
The Reduced Bayesian Net of our detector involves two independent Hidden Markov models (HMM) for the vertical and horizontal bands, each of maximum decomposition level . Thus, for each pixel we have two sequences of hidden states , and , and represent the horizontal and vertical bands respectively, and where each and take the “edge” or “noedge” values. The hidden state sequences have a 1d Markov chain property that models the Persistence property of the wavelets transforms. In this model, the observed wavelet coefficients sequence has been emitted by the hidden state sequence. The conditional distribution of the coefficients, as we said, is a LaplacianGaussian mixture.
The edge detector finds the two hidden state sequences for each pixel given the observed sequences of vertical and horizontal coefficients and combine them in two different ways:

The edge maps from each directional band are combined with the majority vote rule, generating Horizontal and Vertical edge maps.

The Horizontal and Vertical maps are combined with the boolean “or” operator, see Figure LABEL:flow.
Edge maps may be obtained with another choices of deterministic combinations. The Horizontal and Vertical edge maps may be obtained combining the levels of the undecimated wavelet transform with the “and” operator, decision that labels a pixel as edge point only if it is labeled as an edge point in all levels of the wavelet decomposition. Applying the median operator inside each band allows us to repair sequences broke due to noise, if there is enough evidence that they are edge related. The final combination with the “or” operator allows us to keep independent horizontal and vertical edges, as well as corner points. The redundance of the undecimated wavelet transform keeps the information of the diagonal edges in the vertical and horizontal coefficients, suppressing the need of a third HMM model for the Diagonal band.
3.1 Model’s parameters
Each of the HMM models has a set of parameters , which has to be determined from the data. The parameters are the transition probabilities between the two hidden states , the initial probabilities in each state , and the variance parameters and of the emission densities, Laplacian and Gaussian, respectively. The means of the emission densities are set to zero, following the sparsity hypothesis of the model. We will tie the model for each pixel in the same band, considering all sequences of coefficients in each band as independent observations of the same model. Thus for the sample of sequences of random variables (where model the noisy wavelet coefficients), our estimated parameters are the values that maximize :
3.1.1 EM algorithm for HMM parameter estimation
We will compute the parameters using the well known Expectation Maximization paradigm, thus we will iterate until convergence the steps

Step E: compute the conditional expectation ,

Step M: maximize the above expression on .
Theorem 3.1
The EM estimation scheme for our wavelet model consists on the following steps: in the th step, the forward and backward variables and , are computed as
where edge, noedge, respectively. The state probability variables, and , became:
where for any .
Thus the new parameters are :
where is the size of the image.
The proof of the theorem is left to Appendix A. The EM algorithm is a local descent algorithm that can be trapped in local minimum, but as we use it to construct an edge strength mask. We need to verify that the mask contains the true edges of the image, even if it provides, as all statistical edge detectors do, thick edges.
3.2 Sequence of hidden states
Given the sequence of observed coefficients on a given image, and the estimated parameters of the model over the image, we are interested in finding the optimal sequence of hidden states that explains , for each pixel location . That is to say, we want to find the sequence for which .
(1) 
Viterbi algorithm solves this problem, (see Appendix B). We apply this algorithm for each pixelwise sequence of coefficients, in the vertical and horizontal band, using their respective estimated parameters, obtaining binary maps of edge positions, maps for each band, horizontal and vertical.
3.3 Edge strenght Mask
In Sun et al. (2004), strings of labeled coefficients produce a label “edge” only if they have all the evidence in favor of being an edge, that is, the string produced by a branch in the wavelet tree is an unbroken sequence of states. In our case, we compute the most likely state maps that may have emitted the observed image, and decide that a pixel is a candidate to edge point if it is flagged as one in more than one half of the state maps of any of the directional bands considered.
The final edge strength map has a constant rate of false alarms but a low misclassification rate, so the second stage in edge detection, that is made usually with active contours, can successfully locate the real edges.
4 Experimental Results
The appropriate noise model for speckled images is the multiplicative one. We apply the logarithmic transformation to obtain an image with the appropriate additive model, in order to consider Laplacian Gaussian mixtures as the correct emission distributions. Natural optical images with high variance Gaussian noise added are examples that fit the additive noise model and have edges on all scales, providing good challenge for our algorithm. Also, there are several classical databases of images for edge detection which provide benchmark images with ground truth. In Figure 3 we show three different Ground Truth images, made by different people, that show the different scales that edges of natural images have. The database was made adding seven levels of Gaussian noise, , to the test images, downloaded form the Berkeley Database [13].
We will compare our edge detector, WBND, with a simple but effective wavelet based thresholding scheme (HTHW), also based on the undecimated Haar wavelet transform, and Canny edge detector, the usual choice among practitioners. The deterministic algorithm HTHW was included to show the implicit power for edge detection that the projection in wavelet space possesses, and to show the reducing performance effect that noise produces on its detection capability.
In a second experiment, we will work with several SAR images obtained from different sensors, showing different kinds of textures at different resolutions. We will work with an ESAR PRI (high resolution SAR image) of ice land, produced to study seasonal polynyas, regions of open water inside the Artic casquet, an ESAR VV three looks (moderate resolution SAR image) of a region outside the city of Munich, produced to study automatic generation of road maps, and finally, an ESAR full polarimetric complex image without any processing, with several kinds of textures, from a countryside region of Germany, which has been extensively studied in other recent edge detection proposals, Giron et al. (2012).
We will begin detailing the HTHW algorithm we use for making the comparison’s simulation.
4.1 Hard Thresholding Haar wavelet detector (HTHW)
Many wellknown edge detectors, like Prewit, Roberts, Sobel and Canny, are based on the principle of matching local image regions with specific edge patterns. These edge patterns are in fact impulse responses of highpass filters that suppress homogeneous and smooth regions while emphasizing discontinuities. The edge detection is performed by convolving the image with a set of directional patterns, followed by thresholding, Kitanovski et al. (2009). Wavelet based edge detectors also involve these ideas.
We implemented a simple edge detector thresholding the undecimated Haar wavelet decomposition at the first level, which is roughly like computing the first derivative in three directions, horizontal, vertical and diagonal. We label a pixel location as edge if the maximum of the sizes of its related three coefficients surpass a given threshold. We also keep the information about the orientation of the biggest coefficient to apply nonmaximal suppression as Canny’s method does.
Our statistical algorithm, WBND, select the coefficients as being in the state “edge” or “noedge” not only when their amplitude is large (as HTHW does), but also when it is moderately large but still has high probability of being produced by the “edge” emission distribution.
4.2 Initial value selection
The thresholds for Canny and HTHW method are image dependent, and chosen using the Prattquality curves method discussed on FernandezGarcia et al. (2004). We moved the parameters in a one dimensional fashion, 100 steps, and let the quality measures choose the best map. In the case of Canny, we fixed the noise smoothing parameter to the default, and moved only the threshold, using the Matlab function “edge”. In the case of HTHW, there is only one parameter to move, since it does not have a builtin smoothing procedure. Nevertheless, we applied a median filter to the image before computing the edge map, in order to compete in a fair way with Canny.
For our method, WBND, we choose noninformative priors for the state values, and transition probabilities that strongly penalize the change from one state to the other.
The initial variance parameters of the conditional emission distributions, and , were carefully chosen from histograms of the coefficients of each image.
4.3 Quality measures
In this section we define three quality measures well known in the computer vision literature. These measures need ground truth maps to produce the comparison output, which are always available when working with artificial images.
4.3.1 Pratt’s quality measure
Let us consider the ground truth of an image , and another computed edge map over the same image, and the set of edge points of both maps and be the Euclidean distance. Thus, Pratt’s measure is defined by
with .
4.3.2 Baddeley’s error measure
Baddeley’s distance is defined by
with .
4.3.3 Kappa’s quality measure
The index Kappa is defined by
where is the observed proportion of edges and is the expected proportion of edges in the image, one computed over the estimated edge map and the other computed over the ground truth map .
4.4 Results on optical images
We compute 50 edge maps per image sampling the parameter space of Canny and HTWD. We add 7 levels of noise to each of these images, and compute the distance from them to the Ground Truth, generating a 100 point quality curve for each image. The edge map that makes the maximum value over such curve is the map chose to compete with WBND. In Figure 3 and 2 we plot the maximum value of the measure (over the 100 edge maps) as a function of the noise level.
In Figure 3 we can see plots of the three quality measures as functions of seven levels of noise, computed over three images: a simple synthetic image, Star, and two natural images, Giant and Building, with detailed hand made ground truth. The images has very different scaling behavior in their structural edge system. Star, the first panel of Figure 3, is a very simple artificial image, with very thin edges and no texture. The other two images were natural images that depict man made constructions, with edges in several resolutions.
The positive quality measures, Pratt and Kappa, score maximum quality with the 1 value. Baddeley is an error measure, scoring maximum quality with zero. Keeping in mind these ideas, the measure curves have different interpretations from image to image. In general, with low level of noise, the optimized Canny found the boundaries of the objects almost perfectly, but with high level of noise, errors add up and quality worsens. The edges outside the boundary, as well as the texture, are not well resolved by Canny’s method in any image.
This edge strength maps are the first step in boundary and structure detection. The next step would be to apply an active contour method to deliver a one pixel wide, continuous edge structure. In Figure 2 we can see that WBND gives a strong set of lines to follow. The other two methods give too many false alarms outside the boundary of the object. We should stress again that the maps shown are the best maps chosen by the quality measures contrasting them against the Ground Truth.
Pratt  Baddeley  Kappa  

BD Id  Image  HTHW  Canny  WBND  HTHW  Canny  WBND  HTHW  Canny  WBND 
86000  building  0,427  0,529  0,629  1,908  1,734  1,757  0,098  0,173  0,222 
119082  giant  0,563  0,537  0,649  1,846  1,721  2,257  0,190  0,269  0,230 
69015  koala  0,454  0,494  0,606  2,298  2,177  1,985  0,077  0,172  0,281 
187071  mountain  0,391  0,488  0,655  1,740  1,911  1,459  0,109  0,204  0,307 
24077  monk  0,508  0,473  0,558  2,031  2,153  2,605  0,133  0,139  0,177 
42049  bird  0,517  0,473  0,928  1,235  1,182  0,795  0,342  0,391  0,606 
106024  penguin  0,262  0,327  0,476  1,584  1,682  2,051  0,094  0,155  0,229 
300091  surf  0,311  0,436  0,745  1,141  1,347  0,921  0,153  0,227  0,423 
108082  tiger  0,156  0,193  0,208  1,493  2,261  1,916  0,017  0,046  0,064 
Table 1 shows values of the three measures computed over edge maps of 9 images selected from the Berkeley database with the maximum level of noise added, to simulate signal to noise ratios similar to a SAR image. Quality values are computed referring to the most detailed segmentation that the Berkeley database can provide for the selected image, which is usually the boundary of the main objects, with few details. We selected pictures of animals with different textures, fur and feathers, pictures of landscapes and sea waves, and pictures of man made constructions, as buildings, statues and cars. WBND consistently shows better values than the other two methods.
5 Sar image data
5.1 ERS SAR PRI image data of Artic Ice
In this section, we will apply our edge detector to a portion of a high resolution ERS SAR PRI image (25 m spatial resolution), VV polarization, taken on February 15, 1994, on the southern part of the coastal Artic, the Franz Josef Land. This image of ice land includes a polynya visible diagonally from the top left to the bottom right of the image, see Dokken et al.(2002) for details of the imagery and this important environmental problem.
A polynya is a geographical term for an area of open water surrounded by sea ice, i.e. unfrozen sea within the ice pack. Knowledge of the distribution and frequency of coastal polynyas and leads is important in understanding largescale climate processes in the Arctic Basin. Traditionally, satellite passive microwave sensors are used for polynyas studies due to their temporal and spatial coverage. Due to the typical size of polynya features (on average 6 km in width) a relatively high spatial resolution of the satellite sensor is required. The active synthetic aperture radar (SAR) enables high resolution imaging (25 m spatial resolution) of the geographic region of interest during night and day and through the oftenoccurring cloud cover above these features. In Dokken et al. (2002), several algorithms for segmentation of Artic ice images were compared with a new segmenter, called the SAR polynya algorithm, to delineate open water, new ice, and young ice, and to define the size and shape of polynyas. This algorithm uses wavelet projection in two different ways, to recover texture signatures and to extract edges. Such a complex type of algorithm is very common in SAR processing studies, since region detection in the images is very difficult.
Dokken et al. (2002) suggest that enhancements of the SAR polynya algorithm could be to combine different image representations (e.g. curvelets and ridgelets) in order to improve the delineation of young ice types inside the polynya region. We may suggest replacing the last loop of the algorithm, the deterministic wavelet based edge detection, by our Hidden Markov Model approach based on our simulations, where we clearly outperform wavelet thresholding techniques using statistical models. In Figure 4 we can see one of the images studied in Dokken et al. (2002), with a the polynya in it. In that image we selected a small region and computed Canny’s edge map and WBND edge map. We can report that our algorithm delineate the region of thin ice much better than Canny’s algorithm, allowing a segmentation of the region where the polynia should lay.
5.2 ERS SAR image data from city
In this section we apply our algorithm to an ERS SAR image containing urban and woodland areas. We selected a subimage, taken from a 1500 x 4000 ERS1, L band, HH polarization image with an estimated number of looks equal to 2.95, corresponding to Munich, Germany. In Figure 5 we show the image, and the resulting edge map estimating only one model for the whole piece. The image shows several highways and roads, which are reasonably well detected by the algorithm. Pixel resolution is low, around 100 square meters, since it was averaged three times to reduce speckle, so an accurate road map could only be made postprocessing the edge map with subpixel accuracy. Canny’s and Sobel edge maps have thinner lines than the WBND edge map, but do not isolate the regions accurately. In Figure 5 we can see a small portion of woodland which was only well detected by our proposed algorithm.
5.3 Full Polarimetric Complex ESAR imagery
In the previous examples we worked with images with only one polarization band available. In many cases, combination of polarization signatures is important to reveal differences in texture, Mortin et al. (2012).
In this section we will work on a full polarimetric ESAR image that shows three different textures, agricultural land, cities and woodland. The image considered is a raw complex, E band SAR image from a zone near the city of Webling, Bayern, Germany. The image has pixels, and shows a variety of textures and manmade constructions, as we can see in Figure 6. We divided the image in 16 pieces, in order to compute the map in a parallel implementation with the same initial parameters. In Figure 6 we can see the whole map computed over the HH polarization band of the image, and as we found out in the previous subsection, the city roads and woodland regions, corresponding to heterogeneous and extremely heterogenous textures, were well detected by the algorithm. Canny’s algorithm produces thinner edges than our’s, but misses a large part of the boundary of heterogeneous regions.
In Figure 7 we can see another region selected in the original image, a piece of a long runway, and the same region remarked on Canny’s and WBND edge maps. Non of the algorithms could detect that long runway. An explanation for this is that the logarithmic transformation considerably reduces the natural contrast of the HH image. In the two other examples we have shown, the images had a better resolution, so the logarithmic transformation did not diminish the detection power of the algorithm. This image has three other polarizations available, but it has not been processed to reduce speckle, has no number of looks, so it is the noisiest product available.
We selected the region of the runway we showed in Figure7, in all polarizations available, and rerun the algorithm. In the second column of the Figure 7, we can see the runway found, in all images. Nevertheless, the runway is much better defined in the edge map corresponding to VV polarization. Nevertheless, keeping edges from all the different maps made with different polarizations images increases noise in the final map, so combination should be made with a penalization. Manmade constructions were well detected with HH polarization images, low contrast regions (or natural homogeneous areas as the polynyas), were better detected fine tuning VV polarization regions, and VH have not well defined edges.
The proposal of a full polarimetric edge detector is beyond the scope of this paper, but the use of a multidimensional complex wavelet transform could be the way to avoid the logarithmic transformation, maintaining the contrast range of the image in all polarizations. To our best knowledge, no statistical model has been fitted to the coefficients of a multidimensional complex wavelet transform, relating edge and noedge states, and we are currently studying different transformations and algorithmic implementations, searching for the best suited for edge detection.
6 Conclusion and Discussion
In this work we described a new approach to edge map construction in SAR images using hidden Markov Models on the undecimated wavelet domain. Most Hidden Markov models on the wavelet domain have been defined for denoising applications, and they traditionally favor orthogonal transformations. For edge detection, shift invariant transforms are better suited for the task than orthogonal ones, since redundance allows better edge localization and tracking along the wavelet representation. Also, modeling coefficients with a LaplacianGaussian mixture implies that an implicit denoising is taking place at the same time edge related coefficients are selected. This explains the reason why our algorithm maintains a good performance in our simulations when Gaussian noise with high variance was added.
Our algorithm is based on modeling the coefficients of three consecutive levels of the horizontal and vertical detail bands of the undecimated wavelet transform with a three step Markov chain with two possible hidden states, “edge” and “noedge”, which jointly fit a LaplacianGaussian mixture with zero mean. This two independent models conformed a reduced Bayesian Net on the undecimated wavelet domain, which is simple but effective in capturing the two main properties of edge related wavelet coefficients, persistence and non gaussianity. The final edge map is produced combining the most likely sequence of states giving the estimated models, keeping all vertical, horizontal and corner edges.
We investigated the performance of a deterministic wavelet thresholding detector, and the wellknown Canny edge detector along with our WBND algorithm for edge detection in both synthetic images and natural images with available ground truth. From the experimental results on the artificial image, it was observed that Canny detector performed slightly better than WBND and HTHW when moderated Gaussian noise was added, and that our method had an increasingly better performance when higher levels of noise were added. The performance was measured using Pratt’s, Baddeley and Kappa classical quality measures. On natural images with considerably more edges at different resolution sizes than synthetic images, our WBND detector appeared to be the most robust to variations in noise, performing well in all noise configurations tested.
We also report resulting edge maps on three different resolution and polarization SAR images: a VV polarization high resolution ERS SAR PRI image of iced land, an HH polarization L frequency ERS image with three looks of a region outside the city of Munich, and a low resolution SRC raw complex image, the noisiest kind of product, with the HH, VV and VH bands available. Edge detection was successful in all extremely heterogeneous textured regions, but the logarithmic transformation flatten the contrast on the SRC image, introducing detection errors in the homogeneous regions. Increasing contrast on homogeneous regions, and considering all polarizations, we can report that edge maps were better formed on the VV band than on the HH band, and combinations of the two edge maps resulted in very noisy outputs.
The proposal of a full polarimetric edge detector is beyond the scope of this paper, but the use of the undecimated or stationary complex wavelet transform could be the way to avoid the logarithmic transformation, maintaining the contrast range of SRC images. We are aware of Choi et al. (2002) work on Hidden Markov Trees model for Kingsbury’s dual complex wavelet transform, but to our best knowledge, no statistical model has been fitted to the coefficients of the complex wavelet transform of a multidimensional complex image, relating edge and noedge states, and we are currently studying different transformations and algorithmic implementations, searching for the best suited transform for edge detection.
References
 [1] Argenti, F., Bianchi, T., Lapini, A.; Alparone, L., 2012. Fast MAP Despeckling Based on Laplacian Gaussian Modeling of Wavelet Coefficients. IEEE Geoscience and Remote Sensing Letters, 9(1), 13–17.
 [2] Bianchi, T.; Argenti, F.; Alparone, L., 2008. SegmentationBased MAP Despeckling of SAR Images in the Undecimated Wavelet Domain, Geoscience and Remote Sensing, IEEE Transactions on, 46(9), 2728–2742.
 [3] Choi, H., Romberg, J.; Baraniuk, R.; Kingsbury, N., 2002. Hidden Markov tree modeling of complex wavelet transforms. Acoustics, Speech, and Signal Processing, 2000. ICASSP ’00. Proceedings. 2000 IEEE International Conference on, 1, 133–136.
 [4] Crouse, M., Nowak, R. and Baraniuk, R., 1998. Waveletbased statistical signal processing using hidden Markov models, IEEE Trans. Signal Proc., 46(4), 886–902.
 [5] Dokken, S., Winsor, P., Markus, M.,Askne, J., Bjork,G., 2002. ERS SAR characterization of coastal polynyas in the Arctic and comparison with SSM/I and numerical model investigations. Remote Sensing of Environment, 80, 321  335.
 [6] FernandezGarcia, N., MedinaCarnicer,R., CarmonaPoyato, A., MadridCuevas, F., PrietoVillegas, M., 2004. Characterization of empirical discrepancy evaluation measures, Pattern Recognition Letters, 25, 35–47.
 [7] Giron, E., Frery, A.C., CribariNeto, F., 2012. Nonparametric edge detection in speckled imagery, Mathematics and Computers in Simulation, in press.
 [8] Guang, B., Sun, H., 2005. Turbo Iterative Estimation of State Parameters in WaveletDomain Hidden Markov Models of SAR Image. Acta Electronica Sinica ,33(6), 1039–1043.
 [9] Gyaourova, A., C. Kamath, and I. K. Fodor, 2001. Undecimated Wavelet Transforms for Image DeNoising, LLNL Technical report, UCRLID150931.
 [10] Guo, H., 1995. Theory and Applications of the ShiftInvariant, Timevarying and undecimated wavelet transforms. Master Thesis. Rice University.
 [11] Joshi, D.; Jia Li; Wang, J.Z., 2006. A computationally efficient approach to the estimation of two and threedimensional hidden Markov models. Image Processing, IEEE Transactions on, 15(7), 1871– 1886.
 [12] Kitanovski, V., Taskovski, D., Panovski, L., 2008. Multiscale Edge Detection Using Undecimated Wavelet Transform. ISSPIT 2008, 385  389.
 [13] Martin, L., Fowlkes, C., Tal, D., Malik, J., 2001. A Database of Human Segmented Natural Images and its Application to Evaluating Segmentation Algorithms and Measuring Ecological Statistics, Proc. 8th Int. Conf. Computer Vision, 2, 416.
 [14] Mortin, J., Schroder, T., Hansen, A., Holt, B. and McDonald, K., 2012. Mapping of seasonal freezethaw transitions across the panArctic land and sea ice domains with satellite radar, J. Geophysical Res. Oceans, in press.
 [15] Rabiner, L., 1989. A tutorial on hidden Markov models, Proc. IEEE, 77(2), 257–286.
 [16] Romberg,J., Choi, H., Baraniuk, R., 2001. Bayesian TreeStructured Image Modeling Using WaveletDomain Hidden Markov Models IEEE transactions on image processing, 10(7), 1056–1068.
 [17] Starck, J., Fadili,J., and Murtagh, F., 2007. The Undecimated Wavelet Decomposition and its Reconstruction. IEEE transactions on image processing, 16 (2), 297–309.
 [18] Sun, J., Gu, D., Chen,Y., Zhang,S., 2004. A multiscale edge detection algorithm based on wavelet domain vector hidden Markov tree model. Pattern Recognition, 37, 1315 – 1324.
 [19] Tello Alonso, M.; LopezMartinez, C.; Mallorqui, J.J.; Salembier, P., 2011. Edge Enhancement Algorithm Based on the Wavelet Transform for Automatic Edge Detection in SAR Images. Geoscience and Remote Sensing, IEEE Transactions on, 49(1),222–235.
 [20] Zhang, R., Ouyang, W., Cham, W.K., 2009. Image multiscale edge detection using 3D hidden Markov model based on the nondecimated wavelet. Proceeding ICIP’09, 2149–2152.