An Unbiased Detector of Curvilinear Structures.
AbstractThe extraction of curvilinear structures is an important low-level operation in computer vision that has many applications. Most existing operators use a simple model for the line that is to be extracted, i.e., they do not take into account the surroundings of a line. This leads to the undesired consequence that the line will be extracted in the wrong position whenever a line with different lateral contrast is extracted. In contrast, the algorithm proposed in this paper uses an explicit model for lines and their surroundings. By analyzing the scale-space behavior of a model line profile, it is shown how the bias that is induced by asymmetrical lines can be removed. Furthermore, the algorithm not only returns the precise subpixel line position, but also the width of the line for each line point, also with subpixel accuracy.
Introduction
Extracting curvilinear structures, often simply called lines, in digital images is an important low-level
operation in computer vision that has many applications. In photogrammetric and remote
sensing tasks it can be used to extract linear features, including roads, railroads, or rivers, from
satellite or low resolution aerial imagery, which can be used for the capture or update of data
for geographic information systems [1, 2]. In addition it is useful in medical imaging for the
extraction of anatomical features, e.g., blood vessels from an X-ray angiogram [3] or the bones
in the skull from a CT or MR image [4].
The published schemes for line detection can be classified into three categories. The first
approach detects lines by considering the gray values of the image only [5, 6, 7] and uses purely
local criteria, e.g., local gray value differences. Since this will generate many false hypotheses
for line points, elaborate and computationally expensive perceptual grouping schemes have to
be used to select salient lines in the image [8, 9, 10, 7]. Furthermore, lines cannot be extracted
with sub-pixel accuracy.
The second approach is to regard lines as objects having parallel edges [11, 12, 13]. In a
first step, the local direction of a line is determined for each pixel. Then two edge detection
filters are applied in the direction perpendicular to the line, where each filter is tuned to detect
either the left or right edge of the line. The responses of each filter are combined in a non-linear
way to yield the final response of the operator [11]. The advantage of this approach is that
since the edge detection filters are based on the derivatives of Gaussian kernels, the procedure
can be iterated over the scale-space parameter oe to detect lines of arbitrary widths. However,
because special directional edge detection filters that are not separable have to be constructed,
the approach is computationally expensive.
The final approach is to regard the image as a function z(x; y) and extract lines from it by
using various differential geometric properties of this function. The basic idea behind these
algorithms is to locate the positions of ridges and ravines in the image function. These methods
can be further divided according to which property they use.
The first sub-category defines ridges as the point on a contour line of the image, often also
called isohypse or isophote, where the curvature of the contour line has a maximum [4, 14, 15].
One way to do this is to extract the contour lines explicitly, to find the points of maximum
curvature on them, and to then link the extracted points into ridges [14]. However, this scheme
suffers from two main drawbacks. Firstly, since no contour lines will be found for perfectly flat
ridges, such ridges will be labeled as an extended peak. Furthermore, for ridges that have a very
low gradient the contour lines will become widely separated, and thus hard to link. Another
way to extract the maxima of curvature on the contour lines is to give an explicit formula for
that curvature and its direction, and to search for maxima in a curvature image [4, 15]. However,
this procedure will also fail for perfectly flat ridges. While Sard's theorem [16] tells us that for
generic functions such points will be isolated, they occur quite often in real images and lead
to fragmented lines without a semantic reason. Furthermore, the ridge positions found by this
operator will often be in wrong positions due to the nature of the differential geometric property
used, even for images without noise [4, 17].
In the second sub-category, ridges are found at points where one of the principal curvatures
of the image assumes a local maximum [18, 15], which is analogous to the approach taken
to define ridges in advanced differential geometry [19]. For lines with a flat profile it has the
problem that two separate points of maximum curvature symmetric to the true line position will
be found [15]. This is clearly undesirable.
In the third sub-category, ridges and ravines are detected by locally approximating the image
function by its second or third order Taylor polynomial. The coefficients of this polynomial are
usually determined by using the facet model, i.e., by a least squares fit of the polynomial to the
image data over a window of a certain size [20, 21, 22, 23, 24, 25]. The direction of the line is
determined from the Hessian matrix of the Taylor polynomial. Line points are then found by
selecting pixels that have a high second directional derivative perpendicular to the line direction.
The advantage of this approach is that lines can be detected with sub-pixel accuracy without
having to construct specialized directional filters. However, because the convolution masks that
are used to determine the coefficients of the Taylor polynomial are rather poor estimators for
the first and second partial derivatives, this approach usually leads to multiple responses to a
single line, especially when masks larger than 5 \Theta 5 are used to suppress noise. Therefore, the
approach does not scale well and cannot be used to detect lines that are wider than the mask
size. For these reasons, a number of line detectors have been proposed that use Gaussian masks
to detect the ridge points [26, 27, 15]. These have the advantage that they can be tuned for a
certain line width by selecting an appropriate oe. It is also possible to select the appropriate oe
for each image point by iterating through scale space [26]. However, since the surroundings
of the line are not modeled the extracted line position becomes progressively inaccurate as oe
increases.
Evidently, very few approaches to line detection consider the task of extracting the line
width along with the line position. Most of them do this by an iteration through scale-space
while selecting the scale, i.e., the oe, that yields the maximum value to a certain scale-normalized
response as the line width [11, 26]. However, this is computationally very expensive, especially
if one is only interested in lines in a certain range of widths. Furthermore, these approaches
will only yield a relatively rough estimate of the line width since, by necessity, the scale-space
is quantized in rather rough intervals. A different approach is given in [28], where, lines and
edges are extracted in one simultaneous operation. For each line point two corresponding edge
points are matched from the resulting description. This approach has the advantage that lines
and their corresponding edges can in principle be extracted with sub-pixel accuracy. However,
since a third order facet model is used, the same problems that are mentioned above apply.
Furthermore, since the approach does not use an explicit model for a line, the location of the
corresponding edge of a line is often not meaningful because the interaction between a line and
its corresponding edges is neglected.
In this paper, an approach to line detection is presented that uses an explicit model for lines,
and various types of line profile models of increasing sophistication are discussed. A scale-space
analysis is carried out for each of the models. This analysis is used to derive an algorithm
in which lines and their widths can be extracted with sub-pixel accuracy. The algorithm uses
a modification of the differential geometric approach described above to detect lines and their
corresponding edges. Because Gaussian masks are used to estimate the derivatives of the im-
age, the algorithm scales to lines of arbitrary widths while always yielding a single response.
Furthermore, since the interaction between lines and their corresponding edges is explicitly
modeled, the bias in the extracted line and edge position can be predicted analytically, and can
thus be removed. Therefore, line position and width will always correspond to a semantically
meaningful location in the image.
The outline of the paper is as follows. In Section 2 models for lines in 1D and 2D images
are presented and algorithms to extract individual line points are discussed. Section 3 presents
an algorithm to link the individual line points into lines and junctions. The extraction of the
width of the line is discussed in Section 4. Section 5 describes an algorithm to correct the line
position and width to their true values. Finally, Section 6 concludes the paper.
2 Detection of Line Points
2.1 Models for Line Profiles in 1D
Many approaches to line detection consider lines in 1D to be bar-shaped, i.e., the ideal line of
width 2w and height h is assumed to have a profile given by
(1)
However, due to sampling effects of the sensor lines often do not have this profile. Figure 1
shows a typical profile of a line in an aerial image, where no flat bar profile is apparent. There-
fore, let us first consider lines with a parabolically shaped profile because it will make the
derivation of the algorithm clearer and provide us with criteria that should be fulfilled for arbitrary
line profiles. The ideal line of width 2w and height h is given by
(2)
The line detection algorithm will be developed for this type of profile, but the implications of
applying it to bar-shaped lines will be considered later on.
Image data
Approximating profile

Figure

1: Profile of a line in an aerial image and approximating parabolic line profile.
2.2 Detection of Lines in 1D
In order to detect lines with a profile given by (2) in an image z(x) without noise, it is sufficient
to determine the points where z 0 (x) vanishes. However, it is usually convenient to select only
salient lines. A useful criterion for salient lines is the magnitude of the second derivative z 00 (x)
in the point where z 0 Bright lines on a dark background will have z 00 (x) - 0 while
dark lines on a bright background will have z 00 (x) AE 0. Please note that for the ideal line profile
Real images will contain a significant amount of noise, and thus the scheme described above
is not sufficient. In this case, the first and second derivatives of z(x) should be estimated by
convolving the image with the derivatives of the Gaussian smoothing kernel since, under certain,
very general, assumtions, it is the only kernel that makes the inherently ill-posed problem of
estimating the derivatives of a noisy function well-posed [29, 30]. The Gaussian kernels are
given by:
oe (x) =p
2-oe
oe
The responses, i.e., the estimated derivatives, will then be:
oe
oe
oe (x)   f p (x)
s
s
Figure

2: Scale-space behaviour of the parabolic line f p when convolved with the derivatives of
Gaussian kernels for x 2 [\Gamma3; 3] and oe 2 [0:2; 2].
\Gamma2x(OE oe
oe
oe
oe 4
r 00
oe
\Gamma2(OE oe
oe
oe
oe 4 (g 000
oe
oe
where
x
Z
Equations (6)-(8) give a complete scale-space description of how the parabolic line profile
will look like when it is convolved with the derivatives of Gaussian kernels. Figure 2 shows
the responses for an ideal line with bright line on a dark background,
for x 2 [\Gamma3; 3] and oe 2 [0:2; 2]. As can be seen from this figure, r 0
for all oe. Furthermore, r 00
takes on its maximum negative value at
Hence it is possible to determine the precise location of the line for all oe. In addition it can be
seen that the ideal line will be flattened out as oe increases as a result of smoothing. This means
that if large values for oe are used, the threshold to select salient lines will have to be set to an
accordingly smaller value.
Let us now consider the more common case of a bar-shaped profile. For this type of profile
without noise no simple criterion that depends only on z 0 (x) and z 00 (x) can be given since z 0 (x)
and z 00 (x) vanish in the interval [\Gammaw; w]. However, if the bar profile is convolved with the
derivatives of the Gaussian kernel, a smooth function is obtained in each case. The responses
will be:
r- b (x,s,1,1)
s
r- b (x,s,1,1)
s
-22
Figure

3: Scale-space behaviour of the bar-shaped line f b when convolved with the derivatives
of Gaussian kernels for x 2 [\Gamma3; 3] and oe 2 [0:2; 2].
oe
r 00
oe
oe
oe

Figure

3 shows the scale-space behaviour of a bar profile with
is convolved with the derivatives of a Gaussian. It can be seen that the bar profile gradually
becomes "round" at its corners. The first derivative will vanish only at
because of the infinite support of g oe (x). However, the second derivative r 00
b (x; oe; w; h) will not
take on its maximum negative value for small oe. In fact, for oe - 0:2w it will be very close
to zero. Furthermore, there will be two distinct minima in the interval [\Gammaw; w]. It is, however,
desirable for r 00
to exhibit a clearly defined minimum at salient lines are
detected by this value. After some lengthy calculations it can be shown that
has to hold for this. Furthermore, it can be shown that r 00
b (x; oe; w; h) will have its maximum
negative response in scale-space for
3. This means that the same scheme as described
above can be used to detect bar-shaped lines as well. However, the restriction on oe must be
observed.
In addition to this, (11) and (12) can be used to derive how the edges of a line will behave
in scale-space. Since this analysis involves equations which cannot be solved analytically, the
calculations must be done using a root finding algorithm [31]. Figure 4 shows the location of
the line and its corresponding edges for w 2 [0; 4] and oe = 1. Note that the ideal edge positions
are given by From (12) it is apparent that the edges of a line can never move closer
than oe to the real line, and thus the width of the line will be estimated significantly too large for
narrow lines. However, since it is possible to invert the map that describes the edge position,
the edges can be localized very precisely once they are extracted from an image.
The discussion so far has assumed that lines have the same contrast on both sides, which is
rarely the case for real images. For simplicity, only asymetrical bar-shaped lines
f a (x) =? !
x
Line position
Edge position
True line width

Figure

4: Location of a line with width w 2 [0; 4] and its edges for oe = 1.
are considered (a 2 [0; 1]). General lines of height h can be obtained by considering a scaled
asymmetrical profile, i.e., hf a (x). However, this changes nothing in the discussion that follows
since h cancels out in every calculation. The corresponding responses are given by:
r a (x; oe; w;
a
r 00
a
The location where r 0
a
the position of the line, is given by
This means that the line will be estimated in a wrong position whenever the contrast is significantly
different on both sides of the line. The estimated position of the line will be within the
actual boundaries of the line as long as
a
The location of the corresponding edges can again only be computed numerically. Figure 5
gives an example of the line and edge positions for It can be
seen that the position of the line and the edges is greatly influenced by line asymmetry. As a
gets larger the line and edge positions are pushed to the weak side, i.e., the side that posseses
the smaller edge gradient.
Note that (18) gives an explicit formula for the bias of the line extractor. Suppose that we
knew w and a for each line point. Then it would be possible to remove the bias from the line
detection algorithm by shifting the line back into its proper position. Section 5 will describe the
solution to this problem.
Because the asymmetrical line case is by far the most likely case in any given image it is
adopted as the basic model for a line in an image. It is apparent from the analysis above that
failure to model the surroundings of a line, i.e., the asymmetry of its edges, can result in large
errors of the estimated line position and width. Algorithms that fail to take this into account
will fail to return very meaningful results.
a
x
Line Position
Edge Position
True line width
True line position

Figure

5: Location of an asymmetrical line and its corresponding edges with width
and a 2 [0; 1].
2.3 Lines in 1D, Discrete Case
The analysis so far has been carried out for analytical functions z(x). For discrete signals only
two modifications have to be made. The first is the choice of how to implement the convolution
in discrete space. Integrated Gaussian kernels were chosen as convolutions masks, mainly
because the scale-space analysis of Section 2.2 directly carries over to the discrete case. An additional
advantage is that they give automatic normalization of the masks and a direct criterion
on how many coefficients are needed for a given approximation error. The integrated Gaussian
is obtained if one regards the discrete image z n as a piecewise constant function
In this case, the convolution masks will be given by:
n;oe
For the implementation the approximation error is set to 10 \Gamma4 in each case because for images
that contain gray values in the range [0; 255] this precision is sufficient. Of course, other
schemes, like discrete analogon of the Gaussian [32] or a recursive computation [33], are suitable
for the implementation as well. However, for small oe the scale-space analysis will have to
be slightly modified because these filters have different coefficients compared to the integrated
Gaussian.
The second problem that must be solved is the determination of line location in the discrete
case. In principle, one could use a zero crossing detector for this task. However, this would
yield the position of the line only with pixel accuracy. In order to overcome this, the second
order Taylor polynomial of z n
is examined. Let r, r 0 , and r 00 be the locally estimated derivatives
at point n of the image that are obtained by convolving the image with g n
, and g 00
. Then the
Taylor polynomial is given by
The position of the line, i.e., the point where
The point n is declared a line point if this position falls within the pixel's boundaries, i.e., if
and the second derivative r 00 is larger than a user-specified threshold. Please note
that in order to extract lines, the response r is unnecessary and therefore does not need to be
computed. The discussion of how to extract the edges corresponding to a line point will be
deferred to Section 4.
2.4 Detection of Lines in 2D
Curvilinear structures in 2D can be modeled as curves s(t) that exhibit a characteristic 1D line
profile, i.e., f a
, in the direction perpendicular to the line, i.e., perpendicular to s 0 (t). Let this
direction be n(t). This means that the first directional derivative in the direction n(t) should
vanish and the second directional derivative should be of large absolute value. No assumption
can be made about the derivatives in the direction of s 0 (t). For example, let z(x; y) be an image
that results from sweeping the profile f a
along a circle s(t) of radius r. When this image is convolved
with the derivatives of a Gaussian kernel, the second directional derivative perpendicular
to s 0 (t) will have a large negative value, as desired. However, the second directional derivative
along s 0 (t) will also be non-zero.
The only remaining problem is to compute the direction of the line locally for each image
point. In order to do this, the partial derivatives r x , r y , r xx , r xy , and r yy of the image will have
to be estimated, and this can be done by convolving the image with the following kernels
x;oe (x;
oe
oe (y)g oe (x) (26)
xx;oe (x;
oe
xy;oe (x;
oe (x) (28)
oe
The direction in which the second directional derivative of z(x; y) takes on its maximum absolute
value will be used as the direction n(t). This direction can be determined by calculating
the eigenvalues and eigenvectors of the Hessian matrix
r xy r yy
The calculation can be done in a numerically stable and efficient way by using one Jacobi
rotation to annihilate the r xy term [31]. Let the eigenvector corresponding to the eigenvalue of
maximum absolute value, i.e., the direction perpendicular to the line, be given by (n x
1. As in the 1D case, a quadratic polynomial will be used to determine whether
the first directional derivative along (n x ; n y ) vanishes within the current pixel. This point will
be given by
(a) Input image (b) Line points and their response (c) Line points and their direction

Figure

points detected in an aerial image (a) of ground resolution 2m. In (b) the line
points and directions of (c) are superimposed onto the magnitude of the response.
where
r xx n 2
y
Again, (p x is required in order for a point to be declared a line point. As
in the 1D case, the second directional derivative along (n x ; n y ), i.e., the maximum eigenvalue,
can be used to select salient lines.
2.5 Example

Figures

6(b) and (c) give an example of the results obtainable with the presented approach.
Here, bright line points were extracted from the input image given in Fig. 6(a) with
This image is part of an aerial image with a ground resolution of 2 m. The sub-pixel location
of the line points and the direction (n x ; n y ) perpendicular to the line are symbolized
by vectors. The strength of the line, i.e., the absolute value of the second directional derivative
along (n x ; n y ) is symbolized by gray values. Line points with high saliency have dark gray
values.
From figure 6 it might appear, if an 8-neighborhood is used, that the proposed approach
returns multiple responses to each line. However, when the sub-pixel location of each line point
is taken into account it can be seen that there is always a single response to a given line since
all line point locations line up prefectly. Therefore, linking will be considerably easier than in
approaches that yield multiple responses, e.g., [27, 21, 22], and no thinning operation is needed
[34].
3 Linking Line Points into Lines
After individual line pixels have been extracted, they need to be linked into lines. It is necessary
to do this right after the extraction of the line points because the later stages of determining line
width and removing the bias will require a data structure that uses the notion of a left and right
side of an entire line. Therefore, the normals to the line have to be oriented in the same manner
as the line is traversed. As is evident from Fig. 6, the procedure so far cannot do this since line
points are regarded in isolation, and thus preference between two valid directions n(t) is not
made.
3.1 Linking Algorithm
In order to facilitate later mid-level vision processes, e.g., perceptual grouping, the data structure
that results from the linking process should contain explicit information about the lines as
well as the junctions between them. This data structure should be topologically sound in the
sense that junctions are represented by points and not by extended areas as in [21] or [23]. Fur-
thermore, since the presented approach yields only single responses to each line, no thinning
operation needs to be performed prior to linking. This assures that the maximum information
about the line points will be present in the data structure.
Since there is no suitable criterion to classify the line points into junctions and normal line
points in advance without having to resort to extended junction areas another approach has been
adopted. From the algorithm in Section 2 the following data are obtained for each pixel: the
orientation of the line (n x sin ff), a measure of strength of the line (the second
directional derivative in the direction of ff), and the sub-pixel location of the line (p x
Starting from the pixel with maximum second derivative, lines will be constructed by adding
the appropriate neighbor to the current line. Since it can be assumed that the line point detection
algorithm will yield a fairly accurate estimate for the local direction of the line, only three
neighboring pixels that are compatible with this direction are examined. For example, if the
current pixel is and the current orientation of the line is in the interval [\Gamma22:5
only the points (c x +1; c y
\Gamma1), are examined. The choice regarding
the appropriate neighbor to add to the line is based on the distance between the respective sub-pixel
line locations and the angle difference of the two points. Let
be the
distance between the two points and
j, such that fi 2 [0; -=2], be the angle
difference between those points. The neighbor that is added to the line is the one that minimizes
In the current implementation, used. This algorithm will select each line point
in the correct order. At junction points, it will select one branch to follow without detecting the
junction, which will be detected later on. The algorithm of adding line points is continued until
no more line points are found in the current neighborhood or until the best matching candidate
is a point that has already been added to another line. If this happens, the point is marked as a
junction, and the line that contains the point is split into two lines at the junction point.
New lines will be created as long as the starting point has a second directional derivative that
lies above a certain, user-selectable upper threshold. Points are added to the current line as long
as their second directional derivative is greater than another user-selectable lower threshold.
This is similar to a hysteresis threshold operation [35].
The problem of orienting the normals n(t) of the line is solved by the following procedure.
Firstly, at the starting point of the line the normal is oriented such that it is turned \Gamma90 ffi to the
direction the line is traversed, i.e., it will point to the right of the starting point. Then at each line
point there are two possible normals whose angles differ by 180 ffi . The angle that minimizes the
difference between the angle of the normal of the previous point and the current point is chosen
(a) Linked lines and junctions (b) Lines and oriented normals

Figure

7: Linked lines detected using the new approach (a) and oriented normals (b). Lines are
drawn in white while junctions are displayed as black crosses and normals as black lines.
as the correct orientation. This procedure ensures that the normal always points to the right of
the line as it is traversed from start to end.
With a slight modification the algorithm is able to deal with multiple responses if it is assumed
that no more than three parallel responses are generated. For the facet model, for exam-
ple, no such case has been encountered for mask sizes of up to 13 \Theta 13. Under this assumption,
the algorithm can proceed as above. Additionally, if there are multiple responses to the line in
the direction perpendicular to the line, e.g., the pixels
in the example
above, they are marked as processed if they have roughly the same orientation as
termination criterion for lines has to be modified to stop at processed line points instead of line
points that are contained in another line.
3.2 Example

Figure

7(a) shows the result of linking the line points in Fig. 6 into lines. The results are overlaid
onto the original image. In this case, the upper threshold was set to zero, i.e., all lines, no matter
how faint, were selected. It is apparent that the lines obtained with the proposed approach are
very smooth and the sub-pixel location of the line is quite precise. Figure 7(b) displays the way
the normals to the line were oriented for this example.
3.3 Parameter Selection
The selection of thresholds is very important to make an operator generally useable. Ideally,
semantically meaningful parameters should be used to select salient objects. For the proposed
line detector, these are the line width w and its contrast h. However, as was described above,
salient lines are defined by their second directional derivative along n(t). To convert thresholds
on w and h into thresholds the operator can use, first a oe should be chosen according to (13).
Then, oe, w, and h can be plugged into (12) to yield an upper threshold for the operator.

Figure

8 exemplifies this procedure and shows that the presented line detector can be scaled
(a) Aerial image (b) Detected lines

Figure

8: Lines detected (b) in an aerial image (a) of ground resolution 1m.
arbitrarily. In Fig. 8(a) a larger part of the aerial image in Fig. 7 is displyed, but this time with a
ground resolution of 1 m, i.e., twice the resolution. If 7 pixel wide lines are to be detected, i.e.,
3:5, according to (13), a oe - 2:0207 should be selected. In fact, oe = 2:2 was used for
this image. If lines with a contrast of h - 70 are to be selected, (12) shows that these lines will
have a second derivative of - \Gamma5:17893. Therefore, the upper threshold for the absolute value
of the second derivative was set to 5, while the lower threshold was 0:8. Figure 8(b) displays the
lines that were detected with these parameters. As can be seen, all of the roads were detected.
4 Determination of the Line Width
The width of a line is an important feature in its own right. Many applications, especially in
remote sensing tasks, are interested in obtaining the width of an object, e.g., a road or a river, as
precisely as possible. Furthermore, the width can, for instance, be used in perceptual grouping
processes to avoid the grouping of lines that have incompatible widths. However, the main
reason that width is important in the proposed approach is that it is needed to obtain an estimate
of the true line width such that the bias that is introduced by asymmetrical lines can be removed.
4.1 Extraction of Edge Points
From the discussion in Section 2.2 it follows that a line is bounded by an edge on each side.
Hence, to detect the width of the line, for each line point the closest points in the image, to
the left and to the right of the line point, where the absolute value of the gradient takes on
its maximum value need to be determined. Of course, these points should be searched for
exclusively along a line in the direction n(t) of the current line point. Only a trivial modification

Figure

9: Lines and their corresponding edges in an image of the absolute value of the gradient.
of the Bresenham line drawing algorithm [36] is necessary to yield all pixels that this line will
intersect. The analysis in Section 2.2 shows that it is only reasonable to search for edges in a
restricted neighborhood of the line. Ideally, the line to search would have a length of
3oe. In
order to ensure that almost all of the edge points are detected, the current implementation uses
a slightly larger search line length of 2.5oe.
In an image of the absolute value of the gradient of the image, the desired edges will appear
as bright lines. Figure 9 exemplifies this for the aerial image of Fig. 8(a). In order to extract the
lines from the gradient image
x
y
where
the following coefficients of a local Taylor polynomial need to be computed:
e
e
xx
xy
x
e
e
xy
yy
y
e
This has three main disadvantages. First of all, the computational load increases by almost a
factor of two since four additional partial derivatives with slightly larger mask sizes have to be

Figure

10: Comparison between the locations of edge points extracted using the exact formula
(black crosses) and the 3 \Theta 3 facet model (white crosses).
computed. Furthermore, the third partial derivatives of the image would need to be used. This
is clearly undesirable since they are very susceptible to noise. Finally, the expressions above
are undefined whenever e(x; However, since the only interesting characteristic of the
Taylor polynomial is the zero crossing of its first derivative in one of the principal directions,
the coefficients can be multiplied by e(x; y) to avoid this problem.
It might appear that an approach to solve these problems would be to use the algorithm to
detect line points described in Section 2 on the gradient image in order to detect the edges of
the line with sub-pixel accuracy. However, this would mean that some additional smoothing
would be applied to the gradient image. This is undesireable since it would destroy the correlation
between the location of the line points and the location of the corresponding edge points.
Therefore, the edge points in the gradient image are extracted with a facet model line detector
which uses the same principles as described in Section 2, but uses different convolution masks
to determine the partial derivatives of the image [21, 20, 34]. The smallest possible mask size
(3 \Theta 3) is used since this will result in the most accurate localization of the edge points while
yielding as little of the problems mentioned in Section 1 as possible. It has the additional benefit
that the computational costs are very low. Experiments on a large number of images have
shown that if the coefficients of the Taylor polynomial are computed in this manner, they can,
in some cases, be significantly different than the correct values. However, the positions of the
edge points, especially those of the edges corresponding to salient lines, will only be affected
very slightly. Figure 10 illustrates this on the image of Fig. 6(a). Edge points extracted with
the correct formulas are displayed as black crosses, while those extracted with the 3 \Theta 3 facet
model are displayed as white crosses. It is apparent that because third derivatives are used in
the correct formulas there are many more spurious responses. Furthermore, five edge points
along the salient line in the upper middle part of the image are missed because of this. Finally,
it can be seen that the edge positions corresponding to salient lines differ only minimally, and
therefore the approach presented here seems to be justified.
4.2 Handling of Missing Edge Points
One final important issue is what the algorithm should do when it is unable to locate an edge
point for a given line point. This might happen, for example, if there is a very weak and wide
gradient next to the line, which does not exhibit a well defined maximum. Another case where
this typically happens are the junction areas of lines, where the line width usually grows beyond
the range of 2:5oe. Since the algorithm does not have any other means of locating the edge
points, the only viable solution to this problem is to interpolate or extrapolate the line width
from neighboring line points. It is at this point that the notion of a right and a left side of the
line, i.e., the orientation of the normals of the line, becomes crucial.
The algorithm can be described as follows. First of all, the width of the line is extracted for
each line point. After this, if there is a gap in the extracted widths on one side of the line, i.e.,
if the width of the line is undefined at some line point, but there are some points in front and
behind the current line point that have a defined width, the width for the current line point is
obtained by linear interpolation. This can be formalized as follows. Let i be the index of the
last point and j be the index of the next point with a defined line width, respectively. Let a be
the length of the line from i to the current point k and b be the total line length from i to j. Then
the width of the current point k is given by
a
This scheme can easily be extended to the case where either i or j are undefined, i.e., the line
width is undefined at either end of the line. The algorithm sets w in this case, which
means that if the line width is undefined at the end of a line, it will be extrapolated to the last
defined line width.
4.3 Examples

Figure

11(b) displays the results of the line width extraction algorithm for the example image
of Fig. 8. This image is fairly good-natured in the sense that the lines it contains are rather
symmetrical. From Fig. 11(a) it can be seen that the algorithm is able to locate the edges of the
wider line with very high precision. The only place where the edges do not correspond to the
semantic edges of the road object are in the bottom part of the image, where nearby vegetation
causes a strong gradient and causes the algorithm to estimate the line width too large. Please
note that the width of the narrower line is extracted slightly too large, which is not surprising
when the discussion in Section 2.2 is taken into account. Revisiting Fig. 4 again, it is clear that
an effect like this is to be expected. How to remove this effect is the topic of Section 5. A final
thing to note is that the algorithm extrapolates the line width in the junction area in the middle
of the image, as discussed in Section 4.2. This explains the seemingly unjustified edge points
in this area.

Figure

12(b) exhibits the results of the proposed approach on another aerial image of the
same ground resolution, given in Fig. 12(a). Please note that the line in the upper part of the
image contains a very asymmetrical part in the center part of the line due to shadows of nearby
objects. Therefore, as is predictable from the discussion in Section 2.2, especially Fig. 5, the
line position is shifted towards the edge of the line that posesses the weaker gradient, i.e., the
(a) Aerial image (b) Detected lines and their width

Figure

11: Lines and their width detected (b) in an aerial image (a). Lines are displayed in white
while the corresponding edges are displayed in black.
(a) Aerial image (b) Detected lines and their width

Figure

12: Lines and their width detected (b) in an aerial image (a).
upper edge in this case. Please note also that the line and edge positions are very accurate in the
rest of the image.
5 Removing the Bias from Asymmetric Lines
5.1 Detailed Analysis of Asymmetrical Line Profiles
Recall from the discussion at the end of Section 2.2 that if the algorithm knew the true values of
w and a it could remove the bias in the estimation of the line position and width. Equations (15)-
give an explict scale-space description of the asymmetrical line profile f a
. The position l
of the line can be determined analytically by the zero-crossings of r 0
a (x; oe; w; a) and is given
in (18). The total width of the line, as measured from the left to right edge, is given by the
zero-crossings of r 00
a (x; oe; w; a). Unfortunately, these positions can only be computed by a root
finding algorithm since the equations cannot be solved analytically. Let us call these positions
e l
and e r
. Then the width to the left and right of the line is given by v l
The total width of the line is . The quantities l, e l , and e r have the following useful
property:
Proposition 1 The values of l, e l
, and e r
form a scale-invariant system. This means that if both
oe and w are scaled by the same constant factor c the line and edge locations will be given by
cl, ce l
, and ce r
Proof: Let l 1
be the line location for oe 1
and w 1
for an arbitrary, but fixed a. Let oe
and
. Then l
Hence we have l
cl 1
Now let e 1
be one of the two solution of r 00
a
or e r
, and likewise for
, with oe 1;2
and w 1;2
as above. This expression can be transformed to (a \Gamma 1)(e 2
=oe 2). If we plug in oe 2
and w 2
, we see that this expression can
only be fulfilled for e ce 1
since only then will the factors c cancel everywhere. 2
Of course, this property will also hold for the derived quantities v l
, and v.
The meaning of Proposition 1 is that w and oe are not independent of each other. In fact,
we only need to consider all w for one particular oe, e.g., oe = 1. Therefore, for the following
analysis we only need to discuss values that are normalized with regard to the scale oe, i.e.,
v=oe, and so on. A useful consequence is that the behaviour of f a can be
analyzed for oe = 1. All other values can be obtained by a simple multiplication by the actual
scale oe.
With all this being established, the predicted total line width v oe
can be calculated for all w oe
and a 2 [0; 1].

Figure

13 displays the predicted v oe for w oe 2 [0; 3]. It be seen that v oe can grow
without bounds for w oe
# 0 or a " 1. Furthermore, it can be proved that v oe
Therefore,
in Fig. 13 the contour lines for v oe 2 [2; 6] are also displayed.
Section 4 gave a procedure to extract the quantity v oe
from the image. This is half of the
information required to get to the true values of w and a. However, an additional quantity is
needed to estimate a. Since the true height h of the line profile hf a
is unknown this quantity
needs to be independent of h. One such quantity is the ratio of the gradient magnitude at e r
and
e l , i.e., the weak and strong side. This quantity is given by
a
a
It is obvious that the influence of h cancels out. Furthermore, it is easy to convince oneself that r
also remains constant under simultaneous scalings of oe and w. The quantity r has the advantage
that it is easy to extract from the image. Figure 13 displays the predicted r for w oe
[0; 3]. It is
Predicted line width v s5.65.24.84.443.63.22.82.40 0.51.5 23
line width
Predicted gradient ratio r0.80.60.40.20123
(b) Predicted gradient ratio

Figure

13: Predicted behaviour of the asymmetrical line f a
for w oe 2 [0; 3] and a 2 [0; 1]. (a)
Predicted line width v oe . (b) Predicted gradient ratio r.
obvious that r 2 [0; 1]. Therefore, the contour lines for r in this range are displayed in Figure 13
as well. It can be seen that for large w oe
, r is very close to 1 \Gamma a. For small w oe
it will drop to
near-zero for all a.
5.2 Inversion of the Bias Function
The discussion above can be summarized as follows: The true values of w oe
and a are mapped to
the quantities v oe
and r, which are observable from the image. More formally, there is a function
From the discussion in Section 4 it
follows that it is only useful to consider v oe
However, for very small oe it is possible
that an edge point will be found within a pixel in which the center of the pixel is less than 2:5oe
from the line point, but the edge point is farther away than this. Therefore, v oe
[0; 6] is a
good restriction for v oe
. Since the algorithm needs to determine the true values (w oe ; a) from the
observed (v oe ; r), the inverse f \Gamma1 of the map f has to be determined. Figure 14 illustrates that
f is invertible. It displays the contour lines of v oe 2 [2; 6] and r 2 [0; 1]. The contour lines of v oe
are U-shaped with the tightest U corresponding to v 2:1. The contour line corresponding to
actually only the point (0; 0). The contour lines for r run across with the lowermost
visible contour line corresponding to 0:95. The contour line for lies completely on
the w oe
-axis. It can be seen that, for any pair of contour lines from v oe
and r, there will only be
one intersection point. Hence, f is invertible.
To calculate f \Gamma1 , a multi-dimensional root finding algorithm has to be used [31]. To obtain
maximum precision for w oe and a, this root finding algorithm would have to be called at each line
point. This is undesirable for two reasons. Firstly, it is a computationally expensive operation.
More importantly, however, due to the nature of the function f , very good starting values are
required for the algorithm to converge, especially for small v oe
. Therefore, the inverse f \Gamma1 is
computed for selected values of v oe
and r and the true values are obtained by interpolation. The
step size of v oe
was chosen as 0:1, while r was sampled at 0:05 intervals. Hence, the intersection
points of the contour lines in Fig. 14 are the entries in the table of f \Gamma1 . Figure 15 shows the
true values of w oe and a for any given v oe and r. It can be seen that despite the fact that f is very
a

Figure

14: Contour lines of v oe 2 [2; 6] and r 2 [0; 1].
True w s
r0.51.52.5(a) True w oe
True a2.53.54.55.5v s0.20.610.20.61
(b) True a

Figure

15: True values of the line width w oe (a) and the asymmetry a (b).
ill-behaved for small w oe
f \Gamma1 is quite well-behaved. This behaviour leads to the conclusion that
linear interpolation can be used to obtain good values for w oe
and a.
One final important detail is how the algorithm should handle line points where v oe ! 2, i.e.,
f \Gamma1 is undefined. This can happen, for example, because the facet model sometimes gives
a multiple response for an edge point, or because there are two lines very close to each other. In
this case the edge points cannot move as far outward as the model predicts. If this happens, the
line point will have an undefined width. These cases can be handled by the procedure given in
Section 4.2 that fills such gaps.
5.3 Examples

Figure

shows how the bias removal algorithm is able to succesfully adjust the line widths in
the aerial image of Fig. 11. Please note from Fig. 16(a) that because the lines in this image are
fairly symmetrical, the line positions have been adjusted only minimally. Furthermore, it can
be seen that the line widths correspond much better to the true line widths. Figure 16(b) shows
a) Lines detected with bias removal
(b) Detail of (a) (c) Detail of (a) without bias removal

Figure

Lines and their width detected (a) in an aerial image of resolution 1m with the bias
removed. A four times enlarged detail (b) superimposed onto the original image of resolution
m. (c) Comparison to the line extraction without bias removal.
a four times enlarged part of the results superimposed onto the image in its original ground
resolution of 0.25 m, i.e., four times the resolution in which the line extraction was carried out.
For most of the line the edges are well within one pixel of the edge in the larger resolution.

Figure

16(c) shows the same detail without the removal of the bias. In this case, the extracted
edges are about 2-4 pixels from their true locations. The bottom part of Fig. 16(a) shows that
sometimes the bias removal can make the location of one edge worse in favor of improving the
location of the other edge. However, the position of the line is affected only slightly.
a) Lines detected with bias removal
(b) Detail of (a) (c) Detail of (a) without bias removal

Figure

17: Lines and their width detected (a) in an aerial image of resolution 1m with the bias
removed. A four times enlarged detail (b) superimposed onto the original image of resolution
m. (c) Comparison to the line extraction without bias removal.

Figure

17 shows the results of removing the bias from the test image of Fig. 12. Please note
that in the areas of the image where the line is highly asymmetrical the line and edge locations
are much improved. In fact, for a very large part of the road the line position is within one
pixel of the road markings in the center of the road in the high resolution image. Again, a four
times enlarged detail is shown in Fig. 17(b). If this is compared to the detail in Fig. 17(c) the
significant improvement in the line and edge locations becomes apparent.
The final example in the domain of aerial images is a much more difficult image since it
contains much structure. Figure 18(a) shows an aerial image, again of ground resolution 1 m.
This image is very tough to process correctly because it contains a large area where the model
of the line does not hold. There is a very narrow line on the left side of the image that has a very
strong asymmetry in its lower part in addition to another edge being very close. Furthermore,
in its upper part the house roof acts as a nearby line. In such cases, the edge of a line can only
move outward much less than predicted by the model. Unfortunately, due to space limitations
this property cannot be elaborated here. Figure 18(b) shows the result of the line extraction
algorithm with bias removal. Since in the upper part the line edges cannot move as far outward
as the model predicts, the width of the line is estimated as almost zero. The same holds for the
lower part of the line. The reason that the bias removal corrects the line width to near zero is
that small errors in the width extraction lead to a large correction for very narrow lines, i.e., if v oe
is close to 2, as can be seen from Fig. 13(a). Please note, however, that the algorithm is still able
to move the line position to within the true line in its asymmetrical part. This is displayed in

Figures

18(c) and (d). The extraction results are enlarged by a factor of two and superimposed
onto the original image of ground resolution 0.25 m. Please note also that despite the fact that
the width is estimated incorrectly the line positions are not affected by this, i.e., they correspond
very closely to the true line positions in the whole image.
The next example is taken from the domain of medical imaging. Figure 19(a) shows a
magnetic resonance (MR) image of a human head. The results of extracting bright lines with
bias removal are displayed in Fig. 19(b), while a three times enlarged detail from the left center
of the image is given in Fig. 19(c). The extracted line positions and widths are very good
throughout the image. Whether or not they correspond to "interesting" anatomical features
is application dependent. Note, however, that the skull bone and several other features are
extracted with high precision. Compare this to Fig. 19(d), where the line extraction was done
without bias removal. Note that the line positions are much worse for the gyri of the brain since
they are highly asymmetrical lines in this image.
The final example is again from the domain of medical imaging, but this time the input
is an X-ray image. Figure 20 shows the results of applying the proposed approach to a coronary
angiogram. Since the image in Fig. 20(a) has very low contrast, Fig. 20(b) shows the
same image with higher contrast. Figure 20(c) displays the results of extracting dark lines from
Fig. 20(a), the low contrast image, superimposed onto the high contrast image. A three times
enlarged detail is displayed in Fig. 20(d). In particular, it can be seen that the algorithm is very
succesful in delineating the vascular stenosis in the central part of the image. Note also that
the algorithm was able to extract a large part of the coronary artery tree. The reason that some
arteries were not found is that very restrictive thresholds were set for this example. Therefore, it
seems that the presented approach could be used in a system like the one described in [3] to extract
complete coronary trees. However, since the presented algorithm does not generate many
false hypotheses, and since the extracted lines are already connected into lines and junctions,
no complicated perceptual grouping would be necessary, and the rule base would only need to
eliminate false arteries, and could therefore be much smaller.
a) Input image (b) Lines detected with bias removal
(c) Detail of (b) (d) Detail of (b) without bias removal

Figure

18: Lines and their width detected (b) in an aerial image of resolution 1m (a) with bias
removal. A two times enlarged detail (c) superimposed onto the original image of resolution
m. (d) Comparison to the line extraction without bias removal.
6 Conclusions
This paper has presented an approach to extract lines and their widths with very high precision.
A model for the most common type of lines, the asymmetrical bar-shaped line, was developed
from simpler types of lines, namely the parabolic and symmetrical bar-shaped line. A scale-space
analysis was carried out for each of these model profiles. This analysis shows that there
is a strong interaction between a line and its two corresponding edges which cannot be ignored.
The true line width influences the line width occuring in an image, while asymmetry influ-
a) Input image (b) Lines detected with bias removal
(c) Detail of (b) (d) Detail of (b) without bias removal

Figure

19: Lines and their width detected (b) in a MR image (a) with the bias removed. A three
times enlarged detail (c) superimposed onto the original image. (d) Comparison to the line
extraction without bias removal.
ences both the line width and its position. From this analysis an algorithm to extract the line
position and its width was derived. This algorithm exhibits the bias that is predicted by the
model for the asymmetrical line. Therefore, a method to remove this bias was proposed. The
resulting algorithm works very well for a range of images containing lines of different widths
and asymmetries, as was demonstrated by a number of test images. High resolution versions
of the test images were used to check the validity of the obtained results. They show that the
proposed approach is able to extract lines with very high precision from low resolution images.
The extracted line positions and edges correspond to semantically meaningful entities in the im-
(a) Input image (b) Higher contrast version of (a)
(c) Lines and their widths detected in (a) (d) Detail of (c)

Figure

20: Lines detected in the coronary angiogram (a). Since this image has very low con-
trast, the results (c) extracted from (a) are superimposed onto a version of the image with better
contrast (b). A three times enlarged detail of (c) is displayed in (d).
age, e.g., road center lines and roadsides or blood vessels. Although the test images used were
mainly aerial and medical images, the algorithm can be applied in many other domains as well,
e.g., optical character recognition [23]. The approach only uses the first and second directional
derivatives of an image for the extraction of the line points. No specialized directional filters are
needed. The edge point extraction is done by a localized search around the line points already
found using five very small masks. This makes the approach computationally very efficient. For
example, the time to process the MR image of Fig. 19 of size 256 \Theta 256 is about 1.7 seconds
on a HP 735 workstation.
The presented approach shows two fundamental limitations. First of all, it can only be used
to detect lines with a certain range of widths, i.e., between 0 and 2:5oe. This is a problem if the
width of the important lines varies greatly in the image. However, since the bias is removed by
the algorithm, one can in principle select oe large enough to cover all desired line widths and the
algorithm will still yield valid results. This will work if the narrow lines are relatively salient.
Otherwise they will be smoothed away in scale-space. Of course, once oe is selected so large
that neighboring lines will start to influence each other the line model will fail and the results
will deteriorate. Hence, in reality there is a limited range in which oe can be chosen to yield good
results. In most applications this is not a very significant restriction since one is usually only
interested in lines in a certain range of widths. Furthermore, the algorithm could be iterated
through scale-space to extract lines of very different widths. The second problem is that the
definition of salient lines is done via the second directional derivatives. However, one can plug
semantically meaningful values, i.e., the width and height of the line, as well as oe, into (12) to
obtain the desired thresholds. Again, this is not a severe restriction of the algorithm, but only a
matter of convenience.
Finally, it should be stressed that the lines extracted are not ridges in the topographic sense,
i.e., they do not define the way water runs downhill or accumulates [17, 37]. In fact, they are
much more than a ridge in the sense that a ridge can be regarded in isolation, while a line needs
to model its surroundings. If a ridge detection algorithm is used to extract lines, the asymmetry
of the lines will invariably cause it to return biased results.



--R


Update of roads in GIS from aerial imagery: Verification and multi-resolution extraction
An artificial vision system for X-ray images of human coronary trees

Detection of roads and linear structures in low-resolution aerial imagery using a multisource knowledge integration technique
Tracking roads in satellite images by playing twenty questions.
An active testing model for tracking roads in satellite images.
The perception of linear structure: A generic linker.
Linear delineation.
Shape recognition and twenty questions.
Multiscale detection of curvilinear structures in 2-d and 3-d image data
From step edge to line edge: Combining geometric and photometric information.

In So Kweon and Takeo Kanade.
Ridges for image analy- sis
Curves and singularities: A geometrical introduction to singularity theory.

Thin nets and crest lines: Application to satellite data and medical images.
Geometric differentiation for the intelligence of curves and surfaces.
Curve finding by ridge detection and grouping.
Fast recognition of lines in digital images without user-supplied param- eters
Direct gray-scale extraction of features for character recog- nition
Detection of curved and straight segments from gray scale topography.
Computer and Robot Vision
The topographic primal sketch.
Edge detection and ridge detection with automatic scale selection.
Logical/linear operators for image curves.
A common framework for the extraction of lines and edges.
ter Haar Romney
ter Haar Romney
Numerical Recipes in C: The Art of Scientific Computing.
Discrete derivative approximations with scale-space properties: A basis for low-level feature extraction
Recursively implementing the Gaussian and its derivatives.
Extracting curvilinear structures: A differential geometric approach.
A computational approach to edge detection.
Procedural Elements for Computer Graphics.
Tracing crease curves by solving a system of differential equations.
Automatic Extraction of Man-Made Objects from Aerial and Space Images
--TR

--CTR
Jian Chen , Yoshinobu Sato , Shinichi Tamura, Orientation Space Filtering for Multiple Orientation Line Segmentation, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.5, p.417-429, May 2000
Markus Mller , Wolfgang Krger , Gnter Saur, Robust image registration for fusion, Information Fusion, v.8 n.4, p.347-353, October, 2007
Nassir Navab , Yakup Genc , Mirko Appel, Lines in One Orthographic and Two Perspective Views, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.7, p.912-917, July
Jan-Mark Geusebroek , Arnold W. M. Smeulders , Hugo Geerts, A Minimum Cost Approach for Segmenting Networks of Lines, International Journal of Computer Vision, v.43 n.2, p.99-111, July 1, 2001
Thierry Graud , Jean-Baptiste Mouret, Fast road network extraction in satellite images using mathematical morphology and Markov random fields, EURASIP Journal on Applied Signal Processing, v.2004 n.1, p.2503-2514, 1 January 2004
Jong Kwan Lee , Timothy S. Newman , G. Allen Gary, Oriented connectivity-based method for segmenting solar loops, Pattern Recognition, v.39 n.2, p.246-259, February, 2006
Andrew K. C. Wong , Peiyi Niu , Xiang He, Fast acquisition of dense depth data by a new structured light scheme, Computer Vision and Image Understanding, v.98 n.3, p.398-422, June 2005
G. J. Streekstra , R. Van Den Boomgaard , A. W. M. Smeulders, Scale Dependency of Image Derivatives for Feature Measurement in Curvilinear Structures, International Journal of Computer Vision, v.42 n.3, p.177-189, May-June 2001
Antonio M. Lpez , Felipe Lumbreras , Joan Serrat , Juan J. Villanueva, Evaluation of Methods for Ridge and Valley Detection, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.21 n.4, p.327-335, April 1999
E. Cernadas , M. L. Durn , T. Antequera, Recognizing marbling in dry-cured Iberian ham by multiscale analysis, Pattern Recognition Letters, v.23 n.11, p.1311-1321, September 2002
Derek C. Stanford , Adrian E. Raftery, Finding Curvilinear Features in Spatial Point Patterns: Principal Curve Clustering with Noise, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.6, p.601-609, June 2000
lines;medical images;contour linking;feature extraction;aerial images;low-level processing;scale-space;curvilinear structures