6.3. Subpixel Localization of Local Structure
Image data is measured only in the sample points. You may think that because of this we can also localize the structure in images only with pixel distance accuracy. This is not the case fortunately. We can determine the location of local structure with sub pixel accuracy in many cases.
The ability to localize local structure at subpixel accuracy is important in several practical applications:
Finding characteristic points in an image based on a local structure description is dependent on the accurate localization of those points. E.g. in the SIFT procedure the extrema in \(\ell(\v x,s)=s^2( f\ast \nabla^2 G^s)(\v x)\) (the scale normalized Laplacian) are calculated with subpixel accuracy much the same way as done in this section.
In calibrating camera’s and when reconstructing the 3D structure from several 2D views of the scene, it is extremely important to localize points in images with great accuracy. In this section we look at an example to localize X-corners on the block patterns used to calibrate camera’s.
In medical image processing the same part of the human anatomy (say the head) is often ‘seen’ using several imaging devices (e.g. an X-ray and an MRI scanner). For diagnosis it is important to match these two views such that for one point in the head we are able to determine where in the 3D images this point can be found. This registration of the images needs to be done with subpixel accuracy to make accurate measurements and thus a reliable diagnosis possible.
These are just a few examples where the need for subpixel accurate localization arises. In this section we look at sub pixel accurate localization methods for contour lines (isophotes), edges, lines and X-corners (saddle points).
6.3.1. Isophotes — Contour Lines
Consider the isophote \(f(\v x)=c\) where \(c\) is some constant in the range of the scalar function \(f\). Let \(\partial(\set X)\) be the contour of a set \(\set X\subset\setR^d\) (i.e. the set minus its interior), then the (set of) isophote(s) is given by \(\partial(f\geq c)=\partial(f\leq c)\).
For a discrete image \(F\) (let’s assume that \(F\) is the sampled version of \(f\)) the inner contour \(\partial(F\geq c)\) and the outer contour \(\partial(F\leq c)\) are not equal.
Assume that \(\v a\) is a grid point on either the inner or outer contour. A first order Taylor series approximation along the line from \(\v a\) in the gradient direction \(\v e_w\) then is:
where \(w\) is the distance from \(\v a\) in the gradient direction. We are looking for that value of \(w\) such that \(f(\v a + w\v e_w)=c\):
Solving for \(w\) we find:
The subpixel localization of the contour point then is
Expressing the derivatives with respect to the grid aligned coordinates \(x\) and \(y\) we get
Show code for figure
1#
2# Subpixel Accurate Localization of Local Structure
3#
4import numpy as np
5import matplotlib.pyplot as plt
6from scipy.ndimage import binary_erosion, binary_dilation
7from scipy.ndimage import gaussian_filter as gD
8
9plt.close('all')
10
11
12# Isophotes
13
14def innercontour_points(f, threshold):
15 X = f>=threshold
16 Xinner = binary_erosion(X)
17 C = np.logical_xor(X, Xinner)
18 return np.where(C)
19
20def outercontour_points(f, threshold):
21 X = f>=threshold
22 Xouter = binary_dilation(X)
23 C = np.logical_xor(X, Xouter)
24 return np.where(C)
25
26
27N = 32
28X, Y = np.meshgrid(np.arange(N), np.arange(N))
29f = 255 - (X - N // 2)**2 - (Y - N // 2)**2
30
31fig_isophotes, ax = plt.subplots()
32ax.imshow(f, cmap=plt.cm.gray)
33
34for c in [142, 197, 229]:
35 print(np.sqrt(255 - c))
36 circle = plt.Circle((N // 2, N // 2), np.sqrt(255 - c),
37 fill=False, color='green')
38 ax.add_patch(circle)
39 xi, yi = innercontour_points(f, c)
40 ax.plot(xi, yi, 's',
41 markersize=5, markerfacecolor='green', markeredgecolor='k')
42 xo, yo = outercontour_points(f, c)
43 ax.plot(xo, yo, 's',
44 markersize=5, markerfacecolor='red', markeredgecolor='k')
45
46 fx = gD(f, 1, (0, 1))
47 fy = gD(f, 1, (1, 0))
48
49 fi = f[yi, xi]
50 fxi = fx[yi, xi]
51 fyi = fy[yi, xi]
52 sxi = xi + (c - fi) / (fxi**2 + fyi**2) * fxi
53 syi = yi + (c - fi) / (fxi**2 + fyi**2) * fyi
54 ax.plot(sxi, syi, '.', markersize=5, markerfacecolor='green')
55
56 fo = f[yo, xo]
57 fxo = fx[yo, xo]
58 fyo = fy[yo, xo]
59 sxo = xo + (c - fo) / (fxo**2 + fyo**2) * fxo
60 syo = yo + (c - fo) / (fxo**2 + fyo**2) * fyo
61 ax.plot(sxo, syo, '.', markersize=5, color='red')
62
63 plt.axis([6.5, 16.8, 14.5, 3.9])
64 plt.axis('off')
65
66
67# Canny Edges
68
69fig_canny, ax = plt.subplots(figsize=(10,6))
70c = 120
71# g = gD((f >= c).astype(np.float64), 1)
72# g = 1.0 * (f>=c)
73
74c = 10.5
75f = np.abs((X - N // 2)) + np.abs((Y - N // 2))
76g = gD(1.0 * (f <= c), 0.7)
77
78# hidden code below (would spoil the fun of implementing the Canny edge detector)
79
80
81
82if __name__=="__main__":
83 fig_isophotes.show()
84 fig_canny.show()
85fig_isophotes.savefig('source/images/subpixelisophotes.png')
In the example the scale for the Gaussian differentiation was set at 1. Experiments did indicate that accuracy of the subpixel localization in this case was not very dependent on the scale however.
For the discrete contour (either the inner or the outer contour) it is easy to derive the ordering of the points in such a way that the contour is traversed as a curve, in other words it is easy to ‘connect the dots’. The same ordering of points may be used to order the subpixel accurate localized points on the curve.
The classical algorithm to to find a subpixel accurate contour of sampled data is the marching squares algorithm. The marching squares algorithm is somewhat related to the presented one as it is also based on a first order interpolation scheme. In the marching squares algorithm a list of straight line segments is found representing the contour, circumventing the need to find an explicit ordering of the contour points.
6.3.2. Edges
The Canny edge detector is defined with the two conditions:
i.e. we are looking for points with high gradient norm and being zerocrossings in the second order derivative in the (gradient) \(w\)-direction. Let \(\v a\) be a sample grid point where \(f_{w}\gg0\) and where \(f_{ww}\) changes sign in the local neighborhood (a test to check whether the value of \(f_{ww}\) is close to zero is unwise as its value may wildly fluctuate between sample points).
In a local neighborhood of point \(\v a\) we can approximate \(f\) along a line in the gradient direction with a third order Taylor polynomial:
where \(w\) is the coordinate along the line with greyvalues \(g(w)=f(\v a+w \v e_w)\). The zero crossing in the second order derivative \(g_{ww}\) thus is \(g''(w)=f_{ww}(\v a)+w f_{www}(\v a)\) leading to a zerocrossing at:
and thus the zero crossing relative to \(\v a\) is found at:
In previous sections we have allready encountered the derivatives up to order 3 in the \(vw\)-coordinate system expressed as derivatives in \(xy\) (sampling grid axis aligned) coordinate system:
Substitution into the expression for \(\v a_{zc}\) leads to:
6.3.3. Lines
The line detector defined in the 2nd order curvature gauge system is:
The above conditions define bright lines on a dark background. Along the line the curvature is small (the \(\v p\)-direction) and across the line (the \(\v q\)-direction) the curvature is high. To reverse the contrast the role of \(p\) and \(q\) should be interchanged. To locate the midpoint of the line at subpixel accuracy we need to find an analytical description of the luminance function \(f\) in the \(\v e_q\) direction (i.e.\ perpendicular to the line). The Taylor series expansion is:
The line midpoint is characterized as the local maximum in the \(\v e_q\)-direction, i.e.\ \(f_{q}=0\). From the above:
Setting this equal to zero and solving for \(q\) results in the midpoint:
We leave it as an exercise to derive the subpixel line localization in cartesian coordinates.
6.3.4. X-Corners
An X-Corner is what you see on checkerboards. Two bright squares meeting in a corner on a dark background. When we observe the local structure of an X-corner using Gaussian derivatives at scale \(s\), we `see’ a saddlepoint at the corner position.
A saddlepoint is where the principal curvatures are of opposite signs. That is, we are looking for the points where:
Locally around point \(\v a\) we can write the luminance function as:
Let \(\v a\) be a candidate (sampling grid) point where \(f_{pp}(\v a)f_{qq}(\v a)<0\) then we may search for the saddle point location in the neighborhood of \(\v a\) by looking at the point \(\v x\) where the gradient of the above Taylor polynomial vanashish:
Solving for \(\v x\) we find the sub pixel accurate estimate of the X-corner:
6.3.5. Concluding (computational) remarks
In the preceding subsections we have been looking at the location of an \(n\)-th order extremum (zero order for the isophote, second order for the edge, first order for the line and again first order for the X-corner) and in all cases we needed an \((n+1)\)-th order Taylor approximation in a local neighborhood of the extremum.
In all examples it was possible to select candidate sample points (say the point \(\v a\)). For each of these candidate points we calculated the \((n+1)\)-th Gaussian jet and from that the subpixel deviation \(\Delta \v a\) leading to the subpixel accurate location \(\v a + \Delta \v a\). In all examples we have used a very simple test to see whether the resulting estimate was correct. Assuming a unit square sampling grid we accept the subpixel estimate in case \(\|\Delta \v a\|_\infty<\half\), i.e.\ the estimate is accepted in case the subpixel location is close to \(\v a\) (“lies within the pixel corresponding with sample point \(\v a\)”). If this is not the case we expect that there will be another grid point \(\v a'\) that is closer to the subpixel estimate.