Document 1911686

M. van Ginkel, L.J. van Vliet, P.W. Verbeek, M.A. Kraaijveld, E.P. Reding and H.J. Lammers, Robust Curve
Detection using a Radon Transform in Orientation Space Applied to Fracture Detection in Borehole Images, in
R.L. Lagendijk et al. (eds), ASCI’01, Proc. 7th Annual Conference of the Advanced School for Computing and
Imaging (Heijen, The Netherlands, May 30 - June 1), ASCI, Delft, 2000, 299-306.
Robust Curve Detection using a Radon Transform in Orientation
Space Applied to Fracture Detection in Borehole Images
M. van Ginkel
L.J. van Vliet
P.W. Verbeek
Pattern Recognition Group
Delft University of Technology
Delft, The Netherlands
M.A. Kraaijveld
E.P. Reding
H.J. Lammers
Shell International E&P
Technology Applications and Research
Rijswijk, The Netherlands
Abstract
subjects are likely to be unknown to a large part of
the image processing community. We would to emphasize though that the technique is a general method for
detecting parameterised curves, not only fractures in
borehole images.
In our study two image types were used: micro resistivity and acoustic images, but in the current paper we
will restrict ourselves to micro resistivity images. Micro
resistivity images are obtained via many small buttons
(0.1 inch), placed on four pads, that are pressed to the
borehole wall. The amount of current that flows into
the formation via each button is a function of the formation resistivity. These tools are fitted with either 64
or 192 buttons. The pads do not cover the entire wall.
The recorded image contains large vertical stripes for
which no data is available. Micro resistivity imaging
yields a high resolution image of the subsurface.
If the borehole is intersected with a planar event,
the unwrapping of the cylindrical image to a plane will
show a sinusoid. Fractures in the formation, therefore,
can be detected as sinusoids in the image.
A complication is that a fracture is often not an isolated event. Typically, fractures are hidden in a prominent background of bedding planes. The bedding generally manifests itself as stacked sinusoids with a common amplitude (resulting from the bedding dip) and
phase (resulting from the bedding azimuth). The automatic estimation of bedding dip and azimuth can be
performed fairly well as robust estimates of the dip and
azimuth can be obtained by spatial averaging [11]. The
detection of fractures, however, is much more difficult
as they are relatively rare, they have a small coverage
in the image and their orientation deviates significantly
from that of the background. In other words: fractures
are the exceptions where the bedding is the rule. Also,
multiple (intersecting) fractures may be present.
The standard approach for fracture detection, i.e.
computing an edge map and applying a Radon trans-
We present a novel approach to parameterised curve
detection. The method is based on the generalised
Radon transform, which is traditionally applied to a
2D edge/line map. The novelty of our method is the
mapping of the original 2D image to a 3D orientation
space, which then forms the input for the Radon transform. The orientation space representation offers important advantages. It can represent multiple intersecting structures, contains local orientation information,
and offers a better SNR. A Radon transform applied to
this space is automatically constrained by local orientation information, thus reducing the amount of false
detections. The method is applied to fracture detection
in borehole images. The challenge is to detect intersecting (partial) sinusoidal curves in a heterogeneous
and noisy background. The orientation space representation has proven to be instrumental for good detection
of the fractures.
1. Introduction
The detection of curves in digital images is an important step in many image processing problems. In this
paper we propose a novel approach to curve detection.
We restrict ourselves to parameterised curves. The
method is based on a generalised Radon transform, or
equivalently, a Hough transform [1].
The images that we study in this paper are borehole images. Borehole images are acquired by lowering a tool and measuring a physical quantity during
the ascent back up, see for example [12]. Some of the
characteristics of these images are specific to borehole
images, but most of the problematic features (distortions, background noise) are of general importance for
curve detection. Below we use ample space to introduce
borehole imaging and fracture detection, since these
84
form to it, has been widely applied [7, 6, 14]. These
existing methods commonly suffer from two problems.
First, the edge maps used are obtained by crude segmentation approaches. Second, these methods do not
take a number of important properties of fractures in
the images into account. The results of the existing
methods typically show just a few, isolated, detected
fractures, indicating that these methods indeed fail to
perform well.
There are two problems with the detection of these
sinusoids. First, the borehole images suffer from a number of distortions, most of which can be identified in
figure 3 in section 5:
If two oriented structures are present at (ψ, z), with orientations φ1 and φ2 , I [φ] (ψ, z, φ1 ) and I [φ] (ψ, z, φ2 ) will
be large, whereas I [φ] (ψ, z, φ) will be small for other
values of φ.
Contrast this orientation representation with a conventional orientation estimator, which can only represent a single orientation at any given point (ψ, z) [9,
10, 8]:
I(ψ, z) → φ(ψ, z)
(2)
Orientation space offers the following advantages over
the conventional approach:
• Representation of intersecting structures
• Control over the orientation selectivity (this is important because in section 2 we will see that there
is a trade-off between selectivity and localisation.
• Imperfect sinusoids due to:
– Variations in the speed of the tool
• Improved SNR. The orientation space transform
spreads noise over the whole orientation space,
whereas the signal is mapped to very specific locations.
– The borehole is not perfectly cylindrical
– The fracture is not perfectly planar
• Partial sinusoids due to:
• An additional (orientation) constraint for the detection of curves. This is property is specific to
the Radon technique described below.
– Missing data due to the partial coverage
– True partial fractures (either through cementation or drilling induced fractures)
We now come to the key point of the paper. Curves
in the two-dimensional image space are transformed
into other curves in the three-dimensional orientation
space. In this space each point not only encodes a spatial position, but also the local orientation: the normal
of the curve. These new curves have the same parameters as the original curves, since they are obtained
by a deterministic mapping. We can therefore use the
Radon transform to compute a parameter space. The
local maxima in parameter space correspond to the
curves in orientation space and correspond to the sinusoids in the original image. The notion of constraining
a Radon transform is not new [1], but the elegant, geometric, way it is done here for intersecting curves here
is.
A simple example: consider a circle (r cos φ, r sin φ).
The orientation space representation of this circle is
(r cos φ, r sin φ, φ), as depicted in figure 1.
A Radon transform can give a lot of false detections,
because it is free to join curve segments that do not in
reality belong to the same curve. A common way to reduce errors is to impose some constraints on the Radon
transform, for instance making use of local orientation
information. By applying the Radon transform in orientation space, the Radon transform is constrained by
default, since points in this space represent both the
spatial position and local orientation.
• Structures other than fractures introduced by
some of the effects above, especially the partial
coverage and the speed variations
Second, the presence of multiple, intersecting, curves
is a problem for the standard Radon technique. Even
though the Radon transform itself is capable of dealing
with intersecting curves, we must first obtain an edge
map. This edge map must be able to represent intersecting curves, because if we loose information at this
stage we cannot retrieve at at a later stage.
This problem can be solved by adopting an orientation space approach [15, 5]. An orientation space
is a representation for multiple, intersecting, oriented
structures. An orientation space is the result of some
transformation:
I(ψ, z)
→
I [φ] (ψ, z, φ)
(1)
The transformation consists of applying rotated copies
of an orientation selective filter. The exact transformation is given in the next section, but note that there
is no unique orientation transformation. How do we
interpret the orientation space I [φ] ? For a given point
(ψ, z) in the original image, orientation space gives the
strength of a particular orientation φ at that location.
85
beek [15]. Chen uses the orientation of the longest line
in sight to construct his orientation space. In [15] a
local Hough transform is proposed to construct orientation space. The approach we use here, based on a
linear filter bank, was proposed in [5]. A similar approach using Gabor filters was recently used in [2].
A linear orientation space is defined as follows:
800
700
600
500
400
300
200
100
I [φ] (ψ, z, φ) = I(ψ, z) ∗ h(ψ, z; φ)
0
(3)
The orientation selective template filter h(ψ, z) is rotated over an angle φ to obtain h(ψ, z, φ). Convolution
is denoted by ”∗”.
The choice of h(ψ, z) is largely free. One constraint
on h(ψ, z) follows from the observation that we should
deal separately with orientation and scale. This is
achieved by imposing that the Fourier transform of
h(ψ, z) has the following structure:
Figure 1: Orientation space representation of a circle.
Several periods of the curve along the φ-axis are shown.
The first period is displayed as a stem plot. Notice that
the circle is transformed into a double helix.
The approach has two other important properties.
First, it is segmentation-free. It is neither necessary nor
desirable to segment orientation space prior to computing the Radon transform. In fact, finding a good
segmentation procedure would be extremely difficult.
Second, the vertical distortions in the borehole image give a very strong response in orientation space.
Fortunately, we know that the orientation of these distortions is vertical and can simply mask out the corresponding regions in orientation space.
F{h}(f, θ) = Φ[f ] (f )Φ[θ] (θ),
(4)
with f and θ the polar coordinates in Fourier space,
Φ[f ] the radial component of the filter and Φ[θ] the
angular part.
There is one major constraint on h, which follows naturally when we investigate sampled orientation spaces. The Nyquist theorem [13] states that in
order to sample a signal correctly, it must be bandlimited. In order to compute the convolution both the
input image I(ψ, z) and the filter h(ψ, z) must be bandlimited. This is the first minor constraint on h. The
convolution result, i.e. I [φ] (ψ, z, φ) for a fixed value of
φ is bandlimited. Sampling along the x and y axes is
therefore allowed, but what about the φ axis?
In [5] we derived the conditions under which I [φ]
may be sampled along the φ axis. These conditions are
in the form of constraints on the filter h(ψ, z). These
constraints are, in fact, the same as those imposed on
a steerable filter [4]. A steerable filter is a filter for
which the filter response under an arbitrary angle can
be synthesized from the filter response under a fixed
set of orientations. Although the viewpoint of orientation space and steerable filters differs, the latter having
been introduced for adaptive filtering, the mathematical structure is the same. If we sample along the φ
axis, we are in fact computing the filter response for
a finite number of orientations. If we sample properly
we can compute I [φ] for arbitrary φ from the sampled
version, or, equivalently, for an arbitrary orientation of
the filter, which is exactly what a steerable filter is.
The filter we use has the same radial component as
the one in [5], but we use a slightly different angular filter. The angular filter proposed in [5] is the filter with
the optimal orientation selectivity for a given number of
The orientation constraint reduces the amount of
false detections, but not sufficiently. To reduce the
number of false detections and the select the most
salient events, a number of post-processing steps are
performed. In these steps the original image data is inspected to see whether a detected curve is a true curve,
due to noise or is the result of lining up curve segments
which do not belong together. These methods work
well, but are rather ”ad-hoc”.
The structure of the paper is as follows; in section 2
we describe orientation space in more detail. Section 3
is devoted to the Radon transform in orientation space
and its discretisation. The post-processing steps are introduced in section 5, followed by experimental results
in section 6.
Due to space limitations we only present an outline
of our method, with the emphasis on the two key components: the orientation space and the Radon transform applied in orientation space. An extended version
is in preparation.
2. Orientation space
Orientation space was first introduced by Chen and
Hsu [3] and later, independently, by van Vliet and Ver86
Orientation space
Borehole image
Fracture
Response to single
bedding plane
Response to
bedding planes
3
2.5
Response
to fracture
2
1.5
φ
z
1
25
20
0.5
ψ
15
10
0
Bedding planes
ψ
0
Fracture touching
bedding plane
1
2
5
3
4
5
6
z
0
Projection of fracture response
onto bottom plane
Figure 2: On the left: a schematic representation of a borehole image containing bedding and a single fracture. On
the right the corresponding orientation space representation. Note that in orientation space the fracture intersects
the bedding information at just two locations.
that we will detect both even (lines) and odd (edges)
structures.
samples, but the filter response suffers from some secondary responses. Here we use a Gaussian profile with
approximately the same orientation selectivity. This
filter is not truly bandlimited, but the aliasing effects
are negligible. The angular filter is as follows:
(N θ)2
Φ[θ] (θ; N ) = exp(−
),
2π 2
3. The Radon transform
As stated in the introduction, a generalised Radon
transform is used to transform the edge map (orientation space) into a parameter space. In this parameter
space local maxima correspond to sinusoidal curves in
the original image space. The generalised Radon transform I [R] for an arbitrary parameterised curve ~c(s; p~)
is defined as
Z
[R]
I (~
p) = I(~c(s; p~))ds,
(7)
(5)
where N is the amount of filters (i.e. the number of
samples along the φ axis), and is related to the orientation selectivity sa (=the standard deviation of the
Gaussian kernel) as sa = π/N . Note that there is
trade-off between selectivity and localisation. As sa
becomes smaller, smaller orientation differences can be
detected, but the filter becomes more and more elongated in the spatial domain, thus loosing localisation.
The radial component of the filter is given by:
Φ[f ] (f ; fc , bf ) =
|f |
fc
f2c2
b
f
exp(−
f 2 − fc2
)
2b2f
s
with s the distance along the curve, and p~ the parameters of the curve. There are two popular methods of
discretising this equation. The first, which we refer
to as the Radon algorithm, visits all the points p~ in
parameter space and computes the corresponding integrals. The second or Hough algorithm visits all points
in input space and determines for which points p~ the
corresponding curve passes through the selected point
in input space. It then adds the selected point’s contribution to those points p~. The Radon algorithm is the
most suitable of the two for our application.
Instead of finding sinusoidal curves in the original
image space I, we are looking for curves, corresponding
to the sinusoidal curves in I, in orientation space I [φ] .
(6)
which is a bandpass filter with a nearly Gaussian profile
with (band)width bf . It attains its maximum for fc :
Φ[f ] (fc ; fc , bf ) = 1 and goes to zero as f goes to zero.
The resulting filter F{h} = Φ[f ] (f ; fc , bf )Φ[θ] (θ; N )
is a quadrature filter, because Φ[θ] (θ; N ) guarantees
that F{h}(f, θ) is effectively zero for θ outside the interval (−π/2; π/2). Using a quadrature filter ensures
87
Figure 2 shows the schematic structure of a borehole
image and its representation in orientation space. In
this space intersecting fractures and bedding planes lie
apart and the position of points on the curve encodes
the curve’s tangent (local orientation).
A sinusoidal curve
~c(ψ; φd , d, A) = (ψ, z(ψ))T
(8)
z(ψ) = A sin(ψ − φd ) + d
(9)
Another important topic is sampling the parameter space. In our application the sampling distances
along the three parameter axes are determined by the
required accuracy in the parameters of the sinusoids.
All three parameters, amplitude, phase and depth, are
expressed in the physical units of the input image. The
desired accuracy is in the order of one to four pixels.
Tests on a synthetic image (not shown in the current
paper) show that if the parameter space is sampled using a sampling spacing that correspond to two pixels in
the original image, the resulting errors in the parameters found are dominated by the sampling distance.
The resulting error is therefore in the order of two pixels, and within the desired range.
with
where we parameterise ~c as a function of the ψ coordinate, rather than the distance along the curve. The
parameters of the curve are the amplitude A, the phase
φd , and the depth d. The curve ~c is mapped onto the
curve ~c[φ] using the orientation space transformation:


ψ
 (10)
A sin(ψ − φd ) + d
~c[φ] (ψ; φd , d, A) = 
arg(A cos(ψ − φd ) + j)
4. Post-processing
The parameter space computed by the Radon transform of the previous section, contains a large number
of local maxima which do not correspond to actual sinusoidal events in the image, despite the orientation
constraints.
Rather than trying to find criteria which can distinguish between good and bad local maxima in parameter
space directly, we go back to the edge map to decide
whether a particular curve is a good or a bad pick.
We distinguish between three types of bad picks:
Note that the curve is given as a function of ψ rather
than the distance along the curve. The model is given
with respect to the following scaling of the dimensions:
ψ The periodic dimension along the borehole wall,
rescaled to the interval [0; 2π).
z The depth dimension in the original pixel units.
1. Picks due to background noise. Integration along a
curve through a noisy region results in a significant
value in parameter space.
φ In radians.
The result of a curve integral depends on the relative
scaling between the axes. We must therefore be careful
in our choice of coordinate system. The ψ and z dimensions are both physically meaningful dimensions in
the original image space. The integral should therefore
be evaluated with respect to the original scaling, i.e.
the original ψ dimension in pixel units, not rescaled to
[0; 2π). There is no comparable solution for the scaling
along the φ axis. We have therefore chosen to exclude
the ψ axis from contributing to ds in the curve integral:
2. Picks due to connecting up segments that belong
to different fractures.
3. Picks due to detecting a single fracture multiple
times. Given a correctly picked fracture with amplitude A, curves with a slight larger or small amplitude will still give a good fit on the flanks of the
sinusoid.
I [R] (φd , d, A) =
s
Z
(11)
∂z(ψ) 2
[φ]
[φ]
|I |(~c (ψ; φd , d, A)) 1 + (
) dψ
∂ψ
ψ
These erroneous picks would all have less support
than the correct one if the data were well-behaved.
Unfortunately, because of the missing data and vertical distortions, this assumption is no longer true. It
still holds locally, especially for the third type of wrong
pick, but over the image as a whole, good or bad picks
cannot be distinguished on the basis of the integrated
support alone.
To reduce the amount of incorrect picks we have
adopted the following scheme:
The second pitfall lies in the fact that sinusoids with a
larger amplitude have a longer curve length. If the edge
strength at each point along two sinusoids is the same,
the one with the larger amplitude will accumulate more
edge strength in I [R] . This bias for larger amplitudes
can be counteracted by dividing the accumulated edge
strength by the length of (the 2D projection of) the
curve.
For all the local maxima in parameter space in decreasing order with respect to their amplitude, do:
88
Borehole image Geologist’s picks Algorithm’s picks
5 meter
Distortion
Partial
support
Figure 3: Left: a 5 m interval of a micro resistivity image in a heavily fractured formation. Middle: the events
picked by a geologist. Right: events picked by the algorithm.
1. Recompute the corresponding Radon integral.
This value can have changed due to the operations
performed in a previous execution of step 2. If the
value is less than 95% of the original value, reject
this maximum, else
In the second post-processing step, we do not reject
any picks, but do reorder them. This is done by recomputing the Radon integrals for the remaining picks, but
only support that is concentrated within a sufficiently
narrow orientation band around the integration curve
is accumulated. This is done on the original edge map,
not the one in which the support has been erased. The
idea is that support due to a noisy background has no
preferred orientation and is therefore spread along the
φ axis in orientation space, whereas ”true” support is
concentrated along the integration curve.
After this reordering, the first N sinusoids are displayed, with N typically in the range 10 to 25.
2. accept it, and erase the support for this curve and
a band around it from orientation space. The
width of the band corresponds to the response
width of the orientation space filter.
3. Go on to the next maximum (step 1).
This step was designed to remove errors of type 2
and 3 by preventing the same patch of support from
being used for multiple sinusoids. The erasing is performed in orientation space to minimise accidental removal of support for other sinusoids. The magic number 95% performed for all tested images, indicating
that the sensitivity to its setting is low.
5. Experiments
The scheme has been applied to four test images. The
most interesting of these, a microresistivity image, is
depicted in figure 3. This particular image contains
89
very few bedding planes, but has been selected on the
basis of the large number of intersecting fractures. The
distortions mentioned in the introduction are clearly
visible. The image’s dimensions are 292 by 2048 pixels. The orientation space parameters were set to:
N = 17, fc = 0.095 and bf = 0.080. The dimensions of I [R] (φd , d, A) were (73, 512, 50) pixels, with a
sampling distance corresponding to four pixels in the
original image. The maximum amplitude that can be
detected was therefore 4 × 50 = 200 pixels. The three
steps, computing orientation space, computing parameter space, and the post-processing all take roughly 15
minutes on a 300 MHz UltraSparc IIi system, resulting
in a total processing time of about 45 minutes.
Rolling Grants program 94RG12 of the Netherlands
Organization for Fundamental Research of Matter
(FOM).
References
[1] D.H. Ballard. Generalizing the Hough transform
to detect arbitrary shapes. Pattern Recognition,
13(2):111–122, 1981.
[2] J. Chen, Y. Sato, and S. Tamura. Orientation
space filtering for multiple orientation line segmentation. Pattern Analysis and Machine Intelligence, 22(5):417–429, May 2000.
[3] Y.S. Chen and W.H. Hsu. An interpretive model
of line continuation in human visual perception.
Pattern Recognition, 22(5):619–639, 1989.
6. Conclusions
We have presented a new scheme for robust curve detection. There are two essential steps in the analysis
procedure. The first is choosing a correct edge(/line)
map representation, which can cope with multiple intersecting structures. The second is that the local orientation information embedded in this representation
can be used to constrain the Radon transform used
to detect and find the parameters of the curves, in an
elegant geometric way.
The main drawback of the method is that it requires
a large amount of computation, and this becomes worse
when we increase the dimensionality of the problem.
We do not yet not how to properly construct an orientation space for D-dimensional space, but it would have
to be 2D − 1 dimensional, comprised of the original D
spatial dimensions and D − 1 independent orientation
parameters (the Dth parameter would correspond to a
magnitude). The method also scales badly with respect
to the complexity of the curve, because the dimensionality of the parameter space is equal to the number of
parameters of the curve. This is a well known property
of Radon and Hough transforms.
With respect to the application of the scheme to
the detection of fractures it is clear that the detection
performance is very good and close to that of an experienced interpreter. Although an interpreter still has to
QC the detected fractures, a significant reduction of the
interpretation time can be achieved. Moreover, there is
a clear advantage in the application of the automated
techniques as they provide repeatable and auditable
results.
[4] W.T. Freeman and E.H. Adelson. The design
and use of steerable filters. IEEE transactions
on Pattern Analysis and Machine Intelligence,
13(9):891–906, September 1991.
[5] M. van Ginkel, P.W. Verbeek, and L.J. van Vliet.
Orientation selectivity for orientation estimation.
In Proceedings of the 10th Scandinavian Conference on Image Analysis (Lappeenranta, Finland),
volume I, pages 533–537, June 9-11 1997.
[6] K. Glossop, P.J.G. Lisboa, P.C. Russell, A. Siddans, and G.R. Jones. An implementation of the
hough transformation for the identification and labelling of fixed period sinusoidal curves. Computer
Vision and Image Understanding, 74(1):96–100,
April 1999.
[7] J. Hall, M. Ponzi, M. Gonfalini, and G. Maletti.
Automatic extraction and characterisation of geological features and textures from borehole images
and core photographs. In SPWLA 37th Annual
Logging Symposium, June 16-19 1996.
[8] M. Kass and A. Witkin. Analyzing oriented patterns. Computer Vision, Graphics and Image Processing, 37:362–385, 1987.
[9] H. Knutsson. Producing a continuous and distance
preserving 5-d vector representation of 3-d orientation. In IEEE Computer Society Workshop on
Computer Architecture for Pattern Analysis and
Image Database Management, pages 175–182. Miami Beach, Florida, November 18-20 1985.
Acknowledgements
[10] H. Knutsson. Representing local structure using
tensors. In The 6th Scandinavian Conference in
Image Analysis, pages 244–251, June 19-22 1989.
This work was partially supported by Shell International Exploration and Production B.V. and the
90
[11] M.A. Kraaijveld and W.J.M. Epping. Harnessing advanced image analysis technology for quantitative core and borehole image interpretation.
In proceedings of the 41st SPWLA annual symposium, Dallas, page paper VV, June 4-7 2000.
[12] M. Lovell, G. Williamson, and P. Harvey. Borehole
imaging: Applications and case histories, Geological Society special publication, no. 159. Geological
Society, London, 1999.
[13] A.V. Oppenheim, A.S. Willsky, and I.T. Young.
Signals and Systems. Prentice-Hall, 1983.
[14] B.B. Thapa, P. Hughett, and K. Karasaki. Semiautomatic analysis of rock fracture orientations
from borehole wall images. Geophysics, 62(1):129–
137, January/February 1997.
[15] L.J. van. Vliet and P.W. Verbeek. Segmentation of overlapping objects. in L.J. van Vliet and
I.T. Young (eds), Abstracts of the ASCI Imaging
Workshop, Venray, The Netherlands, 5-6, October
1995.
91