Description
This assignment introduces several new concepts, so be patient and careful as you read. You will hopefully learn something new.
You will need the file data.mat and a MATLAB function, all available in a zip file on the homework page.
Image classification and regression problems take as inputs descriptors of the image or of parts of it. These descriptors are called
features. As not all parts of an image are interesting, it is useful to develop an operator that selects image regions that are worth
describing. This selection is difficult in general, because what is “interesting” is a semantic notion that depends on the application. In
this assignment, we will develop a syntactic “interest operator” instead. This operator is given a window (a small square of contiguous
pixels) out of a gray-level image and returns a measure of the local distinctiveness of the image within the window.
In this context, a window W is “locally distinctive” if nearby windows are sufficiently different from W, in a sense to be made more
precise later. To formalize this notion, we first need to define what it means for two image windows to be similar or dissimilar form each
other. Let
D(x, y) = X
z∈Z
2
w(z) [I(z + y) − I(z + x)]2
for x, y ∈ R
2
(1)
be a measure of the dissimilarity between two windows centered at x and y in image I. In this expression, Z is the set of integers, the
summation is over two scalar variables z = (z1, z2)
T
(where the subscript T denotes transposition), and w(z) is a discrete, truncated,
normalized Gaussian function with parameter σ:
w(z) = (
ce
−
z
2
1+z
2
2
2σ2 for − h ≤ z1 ≤ h and − h ≤ z2 ≤ h
0 elsewhere
.
Here,
h = d2.5σe
is the smallest integer that is no smaller than 2.5σ, and c is a real-valued constant such that
X
z∈Z
2
w(z) = 1 .
The “window” is then an n × n square with n = 2h + 1. When x or y are not integer-valued, image values between pixels can be
computed by bilinear interpolation. However, we do not worry about the values of D(x, y) for real-valued arguments in this exercise,
other than in a conceptual sense.
The meaning of D(x, y) is straightforward. The truncated Gaussian w(z) selects small areas of the image around x and y. Discrepancies between corresponding pixel values in those areas are measured by the square of their difference, and weighted (through
multiplication by w(z)) according to how close the pixels are from the center of each window. The weighted discrepancies are then
added to each other to produce a single number D(x, y). If the image values in the two small areas were identical to each other, this
number would be zero, and the number increases as the contents between the two windows grow different from each other.
1. In this problem we get some intuition as to what the function D(x, y) looks like when x is kept fixed and y is varied around x. To
do this, we set
y = x + d
so that d is the displacement between a window centered at x and another one centered at y. We then pick three qualitatively different
windows centered at x
(1)
, x
(2)
, x
(3) in a test image, and plot D(x
(i)
, x
(i) + d) for the 81 integer displacement vectors
d ∈ R = {d | d ∈ Z
2
and d = (d1, d2)
T with − 4 ≤ d1 ≤ 4 and − 4 ≤ d2 ≤ 4}
(the letter ‘R’ stands for the range of displacements). The set R is a square in “displacement space” and D(x
(i)
, x
(i) + d) for fixed x
(i)
is a function from the integer grid in this space into the non-negative reals.
COMPSCI 527
(a) Write a MATLAB function with header
function [w1, u] = gauss(sigma)
that computes a discrete, truncated, one-dimensional n × 1 Gaussian kernel w1 with parameter σ, such that
w = w1 * w1’;
is the two-dimensional Gaussian kernel w(z) defined above. Note that if w1 is normalized to sum to 1, then w sums to 1 as well.
[You may want to check this.]
The second output argument u is the vector -h:h where h is the half-width of the window (same as h in the problem statement
above).
Hand in your code for gauss and show a plot of w1 versus u for σ = 3. Label the axes appropriately.
(b) Write a MATLAB function with header
function D = dissimilarity(I, x, y, sigma)
that computes the dissimilarity D(x, y). The input parameter I is a gray-level image with pixels of type uint8. The vectors x
and y are 2 × 1 vectors with integer coordinates (it is fine to use variables of type double for x and y, but the values must be
integer). The positive scalar sigma is the parameter σ used in the definition of D. The output D is a scalar of type double. Use
your function gauss to generate the weight mask w(z).
The first coordinate of x, y, d is a row index, and the second coordinate is a column index.
Your code should abort with an error message if either x or y is too close to the image boundary and parts of the corresponding
window would fall outside the image. It is OK to use for loops to scan the range R. Make sure you cast either the whole image
or the relevant parts of it to double before doing any calculations with pixel values.
Hand in your code for the function dissimilarity, as well as three mesh plots of D(x
(i)
, x
(i)+d)for the image shadow.jpg
provided with this assignment and for i = 1, 2, 3 and for d in the range R defined earlier. For x
(i) use the three columns of the
matrix X in the data.mat file provided with this assignment (use load data to load X into your workspace). Set σ = 3 for
all plots.
To display a mesh plot use the MATLAB function mesh (use the MATLAB help to figure out how mesh works). Before printing
the mesh plot, set its viewing direction (by hand or with the view function) so that the plot is well visible. What MATLAB calls
the y axis should be marked with label d1 and the x axis should be marked with label d2. The tick marks on the axes should be
labeled from −4 to 4 on each axis. Do not hand in the code used to make the plots.
(c) The three windows centered at x
(i)
for i = 1, 2, 3 in the previous question contain respectively a rather flat part of the image
(i = 1), an edge (i = 2), and a corner (i = 3). You may want to visually explore these windows using MATLAB’s display
facilities.
In what way do the three functions D(x
(i)
, x
(i) +d) differ from each other qualitatively? More specifically, the function D(x, x+
d) is nonnegative and is zero for d = 0, so in some sense it is supposed to look like a bowl (a convex paraboloid) in the vicinity
of d = 0. How well-defined is the bottom of the bowl in each of the three different situations? In answering this question, pay
attention also to the actual values of D (the values on the vertical axis of your plots), and keep in mind that the image is noisy.
Because of this, small fluctuations in the value of D are not significant.
2. We now use the intuition developed with the previous exercise to come up with a definition of “interestingness.” An image window
is “interesting” if the bottom of the bowl formed by D(x, x+d) as a function of d is well defined. When this is the case, there is enough
going on in the window around x that sliding the window in any direction makes the window look different from its original version. In
particular, a flat part of the image is not interesting, nor is a window along a straight edge. A corner is interesting, and so is an image
detail with rich textural content (for instance, a polka-dot pattern, gravel, or a fine checkerboard).
To measure this notion of “interestingness,” one could in principle build the bowl for each window in the image and then somehow
measure how well defined the bottom of the bowl is. A more efficient way to achieve the same result is to note that D(x, x+d) must look
like a paraboloid in an infinitesimally small region around d = 0, because in such a small region one can approximate the function D by
its Taylor series truncated after the second-order term. The derivation that follows computes the smallest curvature of this second-order
term, and calls this curvature the “distinctiveness” (or “interestingness”) of the window centered at x.
Imagine fixing x and seeing D(x, x + d) as a function of d only. If you take d = 0, the two windows centered at x and x + d are
then the same window, and
D(x, x) = 0 .
COMPSCI 527
In addition, since D cannot get any smaller than zero, d = 0 is also a minimum of D(x, x + d) with respect to d, and therefore the
gradient there is zero:
dD(x, x + d)
dd
d=0
= 0 .
If you move away from d = 0 in any direction, the function D cannot decrease. It either stays constant, which means that the window
at x + d looks exactly the same as the window at x, or it increases. In other words, the function D(x, x + d) with x fixed looks like the
bottom of a bowl in a small neighborhood of d = 0: The gradient at the bottom is zero, and the curvature of the bowl may be different
as you move in different directions away from d = 0, but is always nonnegative.
The distinctiveness of the window at x is defined as the curvature of the bowl at d = 0 in the direction in which it curves the least.
Specifically, we say that a window centered at x has distinctiveness λmin(x) if
min
{u∈R
2
| kuk=1}
1
2
d
2D(x, x + au)
da2
a=0
= λmin(x) . (2)
Let us parse this definition. Imagine starting at the bottom of the bowl at d = 0. Now pick an arbitrary direction away from that point,
and keep track of that direction with a unit vector u on the plane. As the positive, real number a increases away from 0, the point
y = x + au
moves away from x in the direction of u. The function D has zero gradient, so to first order it does not change. However, the
second derivative may be positive, which means that the window centered at y is different from the starting window centered at x, and
D(x, x+au) increases, at least initially. Measure the rate of increase with one half1 of the second derivative in the definition (2). Along
some directions u, the second derivative at a = 0 is smaller than along others. Pick the direction u that gives the smallest increase. One
half of the second derivative in that direction is the distinctiveness of the window at x. At this point, you may want to go back to your
plots of D(x, x + d) and check your intuition.
You will now prove that λmin(x) is the smaller of the two eigenvalues of the 2 × 2 matrix
A(x) = X
z∈Z
2
w(z) ∇I(z + x)[∇I(z + x)]T
(3)
where
∇I(x) = ”
∂I
∂x1
∂I
∂x2
#
is the gradient of the image I at x = (x1, x2)
T
. I will take you through the proof, so do not worry: All you need to do is to fill the
blanks.
The expression (1) for D contains the argument z + y, which for y = x + d becomes z + x + d. We start by approximating the
image at z + x + d by its Taylor series around z + x, truncated to the linear term in d:
I(z + x + d) ≈ I(z + x) + [∇I(z + x)]T d .
In this way, D(x, x + d) is approximately a quadratic function of d for small values of d:
D(x, x + d) ≈ Q(x, x + d) = X
z∈Z
2
w(z)
I(z + x) + [∇I(z + x)]T d − I(z + x)
2
=
X
z∈Z
2
w(z)
[∇I(z + x)]T d
2
.
(a) Compute one half of the partial derivatives of Q(x, x + d) with respect to d. Since d has two components (call them d1 and d2),
your result will be a vector with two derivatives:
q(x, d) = 1
2
∂Q(x, x + d)
∂d
=
1
2
” ∂Q(x,x+d)
∂d1
∂Q(x,x+d)
∂d2
#
[Hint: If you are familiar with differential calculus applied to vectors, this is a one-line derivation. If not, just compute the two
derivatives separately, and then package the result into a 2 × 1 vector.]
(b) Find a matrix A(x) such that your previous answer has the form
q(x, d) = A(x) d .
Show all your steps clearly. [Hints: (1) Depending on how you wrote the expression for q(x, d), it may be useful to use the fact
that if s is a scalar and f is a vector, then sf = fs. (2) The matrix in this question is not called A(x) by coincidence.]
1The factor 1/2 is introduced merely for convenience in later calculations.
COMPSCI 527
(c) The Hessian of Q(x, x + d) is the 2 × 2 matrix of second derivatives of Q with respect to d:
H =
∂
2Q(x, x + d)
∂d∂dT
=
∂
2Q(x,x+d)
∂d2
1
∂
2Q(x,x+d)
∂d1∂d2
∂
2Q(x,x+d)
∂d2∂d1
∂
2Q(x,x+d)
∂d2
1
and can be found by differentiating the two entries of q(x, d) each with respect to d1 and d2. That is, if the two entries of q are q1
and q2, we have
1
2
H =
”
∂q1
∂d1
∂q1
∂d2
∂q2
∂d1
∂q2
∂d2
#
.
Use your answer to the previous question to show that
1
2
H = A(x) .
[Hint: Again, if you know how to differentiate with respect to vectors, this is an instant derivation. If not, just spell out the product
A(x)d in terms of the entries of A(x) and d and compute the four derivatives separately.]
Since D ≈ Q, the Hessian of D is approximately A(x). It is not too hard to show that one half of the second derivative in equation
(2) is then
1
2
d
2D(x, x + au)
da2
a=0
≈ u
T A(x)u
(the equality would be exact if the Hessian of D were exactly A(x)). Also, the expression u
T A(x)u is minimized over all unit vectors
u when
u = vmin ,
the eigenvector of A(x) that corresponds to the smaller (λmin(x)) of its two real eigenvalues, and
v
T
minA(x)vmin = λmin(x) .
You need not prove any of these statements. Together, they show that the distinctiveness λmin(x) of a window at x is the smaller
eigenvalue of the matrix A(x) defined in equation (3). This matrix, which is only approximately the Hessian of the dissimilarity D, is
called the Gramian of D weighted by w(z). You may recognize this expression as the image structure tensor defined in the textbook in
section 13.2.2. The selection criterion in equation (13.15) of the book is slightly different from λmin(x), but the idea is very similar.
We can now compute λmin(x) everywhere in an image, and say that a window is interesting (or distinctive) if λmin(x) is greater than
some predefined threshold τ .
3. The function grad provided with this assignment computes the gradient of an image. If you say
g = grad(I);
then g is a cell array with two entries: g{1} is an image that contains the component of the gradient in the vertical direction at every
pixel and g{2} contains the component in the horizontal direction. The function grad also takes an optional argument for the parameter
σ of the kernel, but we will use the default value in this exercise.
(a) Use the function grad and the MATLAB builtin function conv2 with the ’same’ option to write a function with header
function lambdaMin = smallEigenvalue(I, sigma)
that computes an image with the smallest eigenvalue λmin(x) of the Gramian A(x) defined in equation (3) at every pixel in image
I. The second argument sigma is the parameter σ used in the definition of w(z). There is no need to cast I to double, because
grad already does that for its own computations. Use your function gauss to compute the weighting mask w1 such that w =
w1 * w1’.
If you are surprised that conv2 shows up here, note that expression (3) for A(x) is a convolution. Specifically, the product
P(x) = ∇I(x)[∇I(x)]T =
a(x) d(x)
d(x) b(x)
is a 2 × 2 symmetric matrix. So your code will first call grad to compute the two components of ∇I everywhere in the image,
then use matrix operations to compute the three images a(x), b(x), d(x) in P(x). The three distinct entries of
A(x) =
p(x) r(x)
r(x) q(x)
COMPSCI 527
are the convolutions of the three distinct entries of P(x) with the Gaussian kernel w. Make sure you use the separability of the
Gaussian to make these convolutions efficient.
Finally, you can compute the image lambdaMin by recalling that the smaller eigenvalue of A(x) is
λmin(x) = p(x) + q(x) −
p
[p(x) − q(x)]2 + 4r
2(x) .
Use matrix operations to do this computation at once over the entire image, without explicit for loops.
Hand in your code and the image lambdaMin obtained on the image shadow.jpg provided with this assignment. To save ink,
invert the image so that large values are mapped to black and small values to white. You can do this as follows:
clf
imagesc(lambdaMin)
colormap gray
cmap = 1 – colormap;
colormap(cmap);
axis image
axis off
(b) Comment on your image lambdaMin. In particular, what do interesting windows in shadow.jpg look like?
COMPSCI 527