Hough transform and astronomical data analysis

Vol. 40, No. 4, pp. 479485, 1996 Copyright @ 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0083-6656/96 $32.00 + 0.00

Vistas in Astronomy

PII: SOO83-6656(96)00032-3

HOUGH TRANSFORM AND ASTRONOMICAL DATA ANALYSIS PASCAL BALLESTER * European Southern Observatory, Karl-Schwarzschildstr.

2, D-85748 Garching, Germany

The Radon and the Hough transforms are backprojection methods used for shape analysis in particular in the domains of automatic control, medicine, scene analysis, and remote sensing. Recent applications of these techniques in astronomy have demonstrated their value for the determination of astrophysical parameters and for image analysis. This paper reviews recent developments of the Radon and Hough transform work, and analyzes Hough transform application issues for the purpose of astronomical data analysis. Copyright @ 1996 Elsevier Science Ltd. Abstract-

1. INTRODUCTION The Radon transform has been introduced in 1917 as a method to reconstruct an object from its projected views. The transform received much attention from the mathematical physics community and is well-known in the domains of tomography, synthetic aperture radar (SAR), and seismology. Paul Hough patented in 1962 the Hough transform originally as an algorithm for the detection of particle tracks in bubble chambers. Princen et al. (1992) have shown that a formal equivalence between the two transforms could be obtained under particular assumptions and considered more generally the Hough transform as a maximum likelihood method. The Hough transform is a robust method for the detection of multi-dimensional features in images. Computation and storage requirements however grow exponential with the dimension of the parameter space. Most of the literature about the Hough transform attempts to bring solutions to this problem, as well as to the formation of artefacts in the accumulator, also known as correlated noise. Although recent survey articles as in Leavers (1993) indicate no application of the Hough transform in astronomy, publications by Ragazzoni and Barbieri (1993, 1994, 1996) and Ballester (1991, 1994) have illustrated the interest of the Hough transform in this field. We discuss in this article the present applications of the Radon and Hough transforms in astrophysics and review applications of the HT for the analysis of natural images as well as new varieties of Hough transform. *E-mail:[email protected]. 479

P Ballester

480

2. RADON TRANSFORM The Radon transforms is an integral transform written as:

Hb,d)

=

1(x, y) 6 {x cos(p) + y sin(p) - d} dx dy

where the 6 symbol denotes the Dirac delta function. By the Radon transform, a point in the image corresponds to a sinusoid in the parameter space H(p, d) and a point in H(p, d) corresponds to a straight line in the image space. The intensity of peaks in the parameter space depends on the length of the corresponding straight line in the image (Jain, 1989). The Radon transform or derived methods are involved whenever the model of an object needs to be reconstructed from its projections. In space-based astronomy recent applications have been developed in particular for the BATSE experiment of the Compton Gamma Ray Observatory (Zhang et al., 1993) and for the Rotating Slit Aperture Telescope (Touma et al., 1992).

3. HOUGH TRANSFORM The HT of a list of pixels (x;,yi) is represented by: N

Hh,

.... a,) =

2

A

{f(xi, y/l, (al, .... a,,))

(2)

i=l

where A is the Kronecker delta function:

A(t) =

1 ift=O 0 t otherwise

and f((x, y), (al, . . . . a,,)) = 0 the characteristic relation of the sought-for feature. Hough used the slope-intercept parametrisation of the straight line (yi - bxi - a = 0) which presents the disadvantages that both the slope and the intercept are unbounded, also that the method is sensitive to the choice of co-ordinate axes in the picture plane. Duda and Hart (1972) proposed the normal parametrization: f((x, Y), (P, 4) = x cm(p) + y sin(p) - d = 0

(4)

and indicated the method could be generalised to any kind of analytical curve, for example circles. In practice the detection of features in the accumulator array constructed with equation (2) is limited by noise on the sample data values. Thrift and Dunn (1983) suggested to alleviate the problem by transforming each data point into a sinusoidal band:

H(al,

. . . . a,,) = ic

{D(h,yl),

(al...., a,))]

i=l

where D((x,y), (al, .. . . a,))) represents the distance of the locus (x,y) to the curve of parameter (al, . . .. a,) and c(r) is a voting kernel, a positive monotonously decreasing

Astronomical Data Analysis

0

2

4

6

8

0.0

10

481

0.5

1.0

1.5

2.0

Slope

X

Fig. 1. Two sets of colinear points (y = 4.5 + 0.7 x 4 and y = x - 4) and the construction of the HT for straight lines detection. Two accumulation points appear at (0.7,4.5) and (I .o,-4.0).

function of the distance r to the curve. Depending on the criterion to optimise, several voting kernels can be used. For instance,

1-

c(r) = i

2

r

(> -

rln

0

if jr\ 5 r,n if IrI > r,n

(6)

is maximised when the term I:, ri2 reaches its minimum, which is the least-square criterion limited to a domain of influence of radius r,,. 3.1. HT and M-estimates Kyriati (1991) and Princen (1992) have pointed out the formal equivalence between the HT and M-estimates defined by Huber (1981) as any estimate expressed by a minimum problem of the form:

P((Xi9yiJ9

(al,

....

a,))

= min!

(7)

i=l

Although the HT is actually a maximum problem, the equivalence is obtained for p(x) = -D(x) and the choice p(x) = - logf(x) gives the ordinary ML estimate. In this formulation f is called a probability density function. The least-square is an M-estimate for f a Gaussian density. Practical drawbacks of LS regression are well known: 0 it allows one to fit only one feature at a time, it provides accurate parameter estimates in the presence of Gaussian noise, but parameter estimates are not reliable in the presence of outliers, and l errors are only considered on the dependent variable.

l

Hunt et al. (1990) analysed the performance of the HT compared to statistical detection theory processors and indicated that the performance of the HT is not dramatically less than the one of the optimal detection theory processor. Princen et al. (1992) derive a number of conclusions concerning the accuracy of the HT as an M-estimate:

482

I? Ballester

Accuracy is maximised if the HT kernel function is roughly like the ML kernel near its centre. l Accuracy increases with the line-segment length. . In the normal representation, worst-case results are obtained for segments lying far from the foot of the normal. Therefore the choice of a co-ordinate system can significantly improve accuracy. l

3.2. Applications

in Astrophysics

The HT can be compared to M-robust regression methods when it is applied to data points. Ragazzoni and Barbieri (1993, 1994) used the HT for time-series analysis, in particular for astronomical light curves. The light curve of an eclipsing binary star (GW Cephei) has been determined with an improved accuracy by Hough Transform. The new period value was confirmed by recent observations (Ragazzoni and Barbie& 1996). Ballester (1994) applied the HT for automated arc line identification by cross-matching an arc spectrum and a line catalogue, allowing initial inaccuracies of initial guess values (average dispersion, central wavelength) of the order of +30X. The method is fast and reliable with a 2D accumulator for quasi-linear dispersion relations and is robust to missing lines and coverage differences. The method presents some similarities with the projection method described in Zuiderwijk (1995).

4. IMAGE ANALYSIS BY HT As the HT is based on analytical representation of features, one can expect it to work best for the detection of artificial objects presenting regular shapes. Indeed an important number of applications of the HT can be found in industrial control or scene analysis. In this respect the HT can be used to analyse instrument calibration images, such as for the detection of echelle orders (Ballester, 1991). The HT is also used for the analysis of natural images in particular in geology and biology. Detection of geological alignments by HT has been studied by Saether et af. (1994) based on a Moving Window Hough algorithm performing directional filtering on LANDSAT images. Wadge and Cross (1989) use a windowed Hough Transform to identify alignments of volcanic cones. Gradient information allows multi-stage detection of circles: all vectors normal to the circle boundary must intersect at the circle centre. This constraint enables to perform a two dimensional Hough transform to determine the centre of the circle, foliowed by a onedimensional transform (or histogram) to determine the radius. The method is known as 21HT which together with the Gering HT with gradients appear to be efficient both in terms of storage requirements and speed with respect to alternative methods like Standard HT, Generalised Hough transform or Modified Fast Hough Transform (Yuen et al., 1990) These methods have been used in geology for the detection of circular geological features by Cross (1988) Taud and Parrot (1992) and compared to photointerpretation by Capellini et al. (1992). The detection of ellipses would in principle require a five-dimensional parameter space for (X0, Yo. a, 6, 0). Several methods have been developed for the multi-stage detection of ellipses, most of them involving gradient information. Thomas et al. (1992) used the HT for

Astronomical Data Analysis

483

the detection of cell nuclei, with good detection rates. The method does not require image segmentation and is found robust and providing few mistaken identifications.

5. IMPLEMENTATIONS OF THE HOUGH TRANSFORM 5.1. Computational Complexity

The definition of the HT as an efficient signal theory detector indicates that computational load is unavoidable to achieve optimal performance. Storage requirements and computational complexity are however two standard problems that most of the HT literature tries to alleviate. We can list a number of possible strategies: l

l l l l

use physical assumptions to restrict the parameter space or the sub-window of the image on which the HT is performed, develop a parallel implementation, use gradient information to perform a multi-stage transformation, adaptively quantise the parameter space, randomly sample the image features.

An active field of research aims at developing new varieties of HT based on the two latter approaches. 5.2. Noise Another standard problem of the Hough transform is inherent bias and noise in the accumulator. The detection power of the HT as well as the achievable accuracy on the parameter estimates are mostly limited by the following factors: l

l l

Uncorrelated noise: Noise in the observed positions generates systematic errors by spreading and displacing the accumulation peaks. Correlated noise: Coincidental intersections in the parameter space generate artefacts. Quantization errors.. Discretization of the pixel and the parameter spaces introduce quantization errors.

The choice of the quantization steps has been studied by a number of authors and is still an open problem in HT theory. Heuristic determination by means of simulation studies for a given type of problem is often the retained approach. Recently Lam et al. (1994) have shown that the optimal quantization is a function of the length of the lines and therefore that no optimal quantization is possible if lines of different lengths have to be detected. 5.3. Varieties of Hough Transform Computational load, adverse effects of the noise and choice of the quantization steps are the object of most of the research concerning the HT. We will indicate in this article only the most active areas of development and refer to the survey article by Leavers (1993) and the comparison by Kalvanien et al. (1995) for a more complete presentation. Varieties of HT can be separated in two groups, respectively, as non-randomised (also called deterministic) and randomised.

484

I? Ballester

Deterministic transforms involve a one-to-many mapping or a one-to-one mapping using gradient information. An extensive search for a maximum is performed over an accumulator array. One-to-one mapping assumes each boundary point is member of only one feature. Deterministic transforms include multiresolution approaches in which the quantization steps are adapted in a coarse-to-fine strategy and statistical approaches consisting of a refinement of the Hough Transform in the sense of ML theory. Randomised transform, also called Probabilistic HTs are based on randomised sampling of the image points and are the most recent development in HT research. They involve a many-to-one mapping and usually alleviate the problem of choosing quantization steps by using a dynamically linked list structure, virtually providing an accumulator of infinite size and accuracy. These techniques are promising as they can easily be extended to problems of higher dimensionality without exponential growth of the computational cost.

6. CONCLUSIONS Applied to data sets, the HT appears to be a powerful method for the analysis of linear or quasi-linear problems. It performs efficient approximate M-estimates and does not rely strongly on initial guess values. Signal detection theory operators can provide better accuracy at the price of greater computational complexity. The HT has been successfully used in astrophysics for the determination of astrophysical relations and in the domain of instrumental calibrations. For the analysis of astronomical images the transform could be used for directional filtering and for objects detection, provided that the method can be made robust to noise and shape irregularities characteristic of astronomical objects. The possibility to parallelise the method would make it appropriate for processing large data sets.

7. QUESTIONS Ql: G. L. Gimel’Farb: In practice, due to computational errors and noise, the hit peaks in parameter space will be distributed between the neighboring “voxels”. How can we search for such “blurred” peaks (for instance how to choose the proper windows in the parameter space) to avoid artefacts?

Answer: Voting kernels are used to avoid artefacts due to the noise. The choice of the influence domain of the kernel depends on the expected noise. Q2: Y Di Gesit: Which is the best approach to analyze astronomical data via Hough transform? Classical one? Circular? Generalised? Does it depend on the kind of data strongly? Answer. SHT often appears as a precise method. Other approaches like RHT are faster but can show limitations, e.g. for the detection of small segments. Q3: M. Kurtz: Could these methods be used to improve photometry a geometric shape.

of stars? The PSF is

Answer. HT is mostly a geometrical transform. Features must be resolved to be detected. Q4: G. Paturef: Is it possible to estimate the chance that the detected feature is real or not?

Astronomical Data Analysis Answer: The intensity

posteriori features.

probability

485

of a maximum is proportional to the length of the feature. A (based on the HT expression) provides the likelihood ratio of

References Ballester P. (1991) In Proc. 3rd ESOIST-ECF Data Analysis Workshop (Edited by Grosbol P. J. and Warmels R. H.), p. 23. Ballester P. (1994) Astron. Astrophys. 286, 101 l-1018. Capellini V., Fini S., Harrigan E. and Mecocci A. (1992) Visual Forum (Edited by C. Arcelli et al.), pp. 119-126. Plenum Press, New York. Cross A. M. (1988) Znt. J; Remote Sensing 9, 1519-1528. Duda R. 0. and Hart P. E. (1972) Graphics and image processing. Commun. ACM 15, 11-15. Huber P. J. (1981) Robust Statistics. Wiley, New York. Hunt D. J., Nolte L. W., Reibman A. R. and Ruedger W. H. (1990) CVGZP 52,38&401. Jain A. K. (1989) Fundamentals of Digitaf Image Processing. Prentice-Hall, New York. Kalvanien H., Hirvonen P., Xu L. and Oja E. (1995) Image and Vision Computing 13,239252.

Kyriati N. and Bruckstein A. M. (1991) IEEE Trans. Pattern Analysis and Machine Zntelligence 13,602-606.

Lam W. C. Y., Lam L. T. S., Yuen K. S. Y. and Leung D. N. K. (1994) Pattern Recog. Lett. 15, 1127-l 135.

Leavers V F. (1993) CVGZP Image Understanding 58, 250. Princen J., Illingworth J. and Kittler J. (1992) JMZV 1, 153-168. Ragazzoni R. and Barbieri C. (1993) In Conf Applic. Time Series in Astron. Meteorol. (Edited by 0. Lessi), p. 233. Ragazzoni R. and Barbieri C. (1994) PASP 106,683-687. Ragazzoni R. and Barbieri C. (1996) Commissions 27 and 42 of the IAU Information Bulletin on Variable Stars, IBVS#4293. Saether B., Rueslatten H. and Gronlie A. (1994) In Proc. Znt. Geosci. Remote Sensing Symp. (IGARSS’94) 1Pasadena, Calif. Taud H. and Parrot J. F. (1992) Znt. 1 Remote Sensing 13, 319-335. Thrift P. R. and Dunn S. M. (1983) CVGZP 21,383. Thomas A. D. H., Davies T. and Luxmoore A. R. (1992) In Proc. ZEE Coffoq. Applic. Images Processing in Mass Health Screening Digest, London. Touma H., Kadiri S., Martin F., Aime C., Darcourt J. and Lanteri H. (1992) ESA, Targets for Space-Based Znterferometry, pp. 197-200. Wadge G. and Cross A. M. (1989) Int. J Remote Sensing 10,455474. Yuen H. K., Princen J., Illingworth J. and Kittler J. (1990) Image and Vision Computing 8, 71-77.

Zhang S. N., Fishman G. J., Harmon B. A., Paciesas W. S., Rubin B. C., Meegan C. A., Wilson C. A. and Finger M. H. (1993) Bull. Am. Astron. Sot. 182, H71.05. Zuiderwijk E. (1995) Astron. Astrophys. Suppl. Ser. 112, 537-549.