Computational Vision
U. Minn. Psy 5036

Assignment #3
Linear systems, neural systems & edge detection

Initialization

Off[General::spell1];
SetOptions[ArrayPlot,Mesh->False,AspectRatio->Automatic,PlotRange->All,ColorFunction->"GrayTones",DataReversed->True,ImageSize->Small];
SetOptions[ListPlot3D,Mesh->False,AspectRatio->Automatic,PlotRange->All];

Goals: Learn: 
	1) Review vectors and matrices
	2) practice using convolution for spatial filtering
	3) practice using the Fast Fourier Transform for spatial filtering
	4) learn one example for how to extract edges from the outputs of model V1 simple cells

Linear algebra review

There is a lot of content in this section, but no programming exercises. To complete the Linear algebra review “exercise” here,  you could just select all the cells and execute them. But I encourage you to go through the material cell by cell as either a review, or as a basic introduction to linear algebra through vectors and matrices.

I will check to see if you’ve executed the input cells.  And if you’ve explored by changing values or functions, that is fine--in fact, great!

Remember that if you change variables or function definitions, this could have unexpected effects later if the same symbols are used. You can use Clear[ ] to reset variables and functions.

Vector operations

Dimension of a vector.

You can get the dimensionality of a vector using  Dimensions[], or Length[].

v = {2.1, 3, -0.45, 4.9};

Dimensions[v]

Dimensions[], will give you the dimensions of a matrix, while  Length[] tells you the number of elements in the list. For example,

M = {{2,4,2}, {1,6,4}};

Length[M]

Try comparing  Length[M] with Dimensions[M].

Transpose of a vector.

The transpose of a column vector is just the same vector arranged in a row. However, because of the way Mathematica uses lists to represent vectors you don't have to distinguish between row and column vectors. In standard math notation, transpose of a vector x, is often written xT. You can see a vector in column form by typing v//MatrixForm, or:

MatrixForm[v]

Vector addition is accomplished by simply adding the components of each vector to make a new vector. Note that the vectors all have the same dimension.

a = {3,1,2};
b = {2,4,8};
c = a + b

Vectors can be multiplied by a constant. We saw an example of this earlier.

2 a

Metric length of a vector.

It is unfortunate terminology, but Length[] does NOT give you the metrical length of the vector. In order to get the length of a vector, you calculate the Euclidean distance from the origin to the end-point of the vector. We get this by squaring each component, adding up the squares, and taking the square root. First, we will do this using the Apply[] function, where the Plus operation is applied to all the elements of the list. Note that the operation of exponentiation is "listable", that is it is applied to each element of the vector:

a^2

What is a a ?

N[Sqrt[Apply[Plus, a^2]]]

If you wish, you can define your own function to apply to the list. What we have just calculated is the square root of the dot product or inner product of a with itself. The length of a vector a is often written as |a| in standard math notation. In the next section, we use the inner or dot product to calculate the metric length of a vector.

Inner product (or “dot” product). To calculate the inner product of two vectors, you multiply the corresponding components and add them up:

u = {u1,u2,u3,u4};
v = {v1,v2,v3,v4};
u.v

The inner product is also called the dot product. Later we will see what is meant by outer product. The inner product between two vectors a and b is written either as:
                                                   a.b or [a,b], or aTb
Mathematica  uses the dot notation.

One use of the inner product is to calculate the length of a vector. a.a  is just the sum of the squares of the elements of a, so gives us another way of calculating the length of a vector.

N[Sqrt[a.a]]

Let's define a function that will return the length of a vector, x:

vectorlength[x_] := N[Sqrt[x.x]]

Compare this function with the built-in function, Norm[]

Clear[x1,x2,x3]
Norm[{x1,x2,x3}]
Sqrt[{x1,x2,x3}.{x1,x2,x3}]

Projection. The dot product, a.b, is equal to:
                                                   |a| |b| cos(angle between a and b)

The output of a linear neuron model is the dot product between an input vector and a weight vector. Both the weight and input lists can be thought of as vectors in an n-dimensional space. Suppose the weight vector has unit length.  Recall that you can normalize any vector to unit length by dividing by its length:

v = v/Sqrt[v.v] ;

Geometrically, we can think of the output of a neuron as the projection of the activity of the neuron input activity vector onto the weight vector direction. Suppose the input vector is already perpendicular to the weight vector, then the output of the neuron is zero, because the cosine of 90 degrees is zero. As you found or will find with the cross-correlator of Problem Set 1, the further the input pattern is away from the weight vector, as measured by the cosine between them, the poorer the "match" between input and weight vectors, and the lower the response.

Here are three lines of code that calculate the two-dimensional vector z in the direction of w,  with a length determined by "how much of x projects in the w direction":

x = {1,2};
w = N[{2/Sqrt[5],1/Sqrt[5]}];
z = (x.w) w ;

Show[Graphics[{Line[{{0,0},x}],Line[{{0,0},z}],{Dashing[{0.03`,0.03`}],Line[{x,z}]},Text[Style["w",{"Helvetica-Bold",18}],w,{0,-1}],Text[Style["x",{"Helvetica-Bold",18}],x,{-2,0}],Text[Style["z",{"Helvetica-Bold",18}],z,{0,-1}],{AbsoluteThickness[3],Line[{{0,0},w}]}}],Axes->True,AspectRatio->1]

Note the built-in function:

?Projection

Angle between two vectors and orthogonality: Similarity measure between patterns

Often we will want some measure of the similarity between two patterns of neural firings. As we have just seen, one measure of comparison is the degree to which the two state vectors point in the same direction. The cosine of the angle between two vectors is one possible measure:

cosine[x_,y_] := x.y/(vectorlength[x] vectorlength[y])
cosine[a,b]

Note that if two vectors point in the same direction, the cosine of the angle between them is 1:

a = {2,1,3,6};
b = {6, 3, 9, 18};
cosine[a,b]

Try verifying that w and z from the previous section point in the same direction.
If two vectors point in the opposite directions, the cosine of the angle between them is -1:

a = {-2,-1,-3,-6};
b = {6, 3, 9, 18};
cosine[a,b]

Two vectors may point in the same direction, but could be quite different because they have different lengths. Another measure of similarity is the length of the difference between two vectors:

vectorlength[a - b]
Norm[a-b]//N

Orthogonality. The case where vectors are at right angles to each other is an important special case that is worth spending a little time on. Consider an 8-dimensional space. One  very familiar set of orthogonal vectors is the following:

u1 = {1,0,0,0,0,0,0,0};
u2 = {0,1,0,0,0,0,0,0};
u3 = {0,0,1,0,0,0,0,0};
u4 = {0,0,0,1,0,0,0,0};
u5 = {0,0,0,0,1,0,0,0};
u6 = {0,0,0,0,0,1,0,0};
u7 = {0,0,0,0,0,0,1,0};
u8 = {0,0,0,0,0,0,0,1};

Each vector has unit length, and it is easy to see just by inspection that the inner product between any two is zero. On the other hand, here is another set of 8 vectors in 8-space for which it is not immediately obvious that they are all orthogonal. These vectors are called Walsh functions:

v1 = {1, 1, 1, 1, 1, 1, 1, 1};
v2 = {1,-1,-1, 1, 1,-1,-1, 1};
v3 = {1, 1,-1,-1,-1,-1, 1, 1};
v4 = {1,-1, 1,-1,-1, 1,-1, 1};
v5 = {1, 1, 1, 1,-1,-1,-1,-1};
v6 = {1,-1,-1, 1,-1, 1, 1,-1};
v7 = {1, 1,-1,-1, 1, 1,-1,-1};
v8 = {1,-1, 1,-1, 1,-1, 1,-1};

You can calculate the inner products between any two, and you will find out that they are all zero. Note that with the first set of vectors, {ui}, you can tell which vector it is just by looking for where the 1 is. For the second set, {vi}, you can't tell by looking at just one component. For example, the first component of all of the Walsh functions has a 1. You have to look at the pattern to tell which Walsh function you are looking at.

u3.u7
v3.v7

Suppose for the moment that we want to assign meaning to each of the patterns--each pattern is a code for some thing, like "grandma Tompkins", "grandma Wilke", and so forth. If we use the u's, then we could look for the one neuron that lights up to find out which grandma it is representing--then neuron activity represented, for example, by the third element of the pattern could mean "grandma Wilke". This strategy wouldn't work if we encoded a collection of grandmothers using the v's. The v's give us a simple example of what is sometimes referred to as a distributed code. The w's are examples of a grandmother cell code. The reason for this obscure terminology can be traced to earlier debates on whether there may be single cells in the brain whose firing uniquely determines the recognition of one's grandmother.

Orthonormality. The Walsh set is orthogonal, but they are not of unit length. We have already seen some of the advantages of working with unit length vectors. The general issue of normalization comes up all the time in neural networks both in terms of limiting overall neural activity, and limiting synaptic weights. So it is sometimes convenient to normalize an orthogonal set, producing what is known as an orthonormal set of vectors:

w1 = v1/vectorlength[v1];
w2 = v2/vectorlength[v2];
w3 = v3/vectorlength[v3];
w4 = v4/vectorlength[v4];
w5 = v5/vectorlength[v5];
w6 = v6/vectorlength[v6];
w7 = v7/vectorlength[v7];
w8 = v8/vectorlength[v8];

Executing the next cell illustrates the above points:

{u4.u4,Norm[u4]}
{v4.v4,Norm[v4],Sqrt[v4.v4]}
{w4.w4,Norm[w4]}

Vector representations, linear algebra

The issue of how information is to be represented is fundamental in the information sciences generally, as well as for neural network theory. A pattern of activity over a set of neurons is presumed to mean something, and there are different ways of coding the same meaning. But different codes have different properties. A code may not be sufficient to uniquely code all the possible things we need to represent. A code could be redundant and have more than one way of representing the same thing. This section continues with our review of the basics of vector and linear algebra by going a little more deeply into the subject. The pay-off will be some mathematics that provides intuition about issues of neural representation. You can think of this as a first lesson in the "psychology of linear algebra".

Basis sets

It is pretty clear that given any vector whatsoever in 8-space, you can specify how much of it gets projected in each of the eight directions specified by the unit vectors v1, v2, ...v8. But you can also build back up an arbitrary vector by adding up all the contributions from each of the component vectors. This is a consequence of vector addition and can be easily seen to be true in 2 dimensions. We can verify it ourselves. Pick an arbitrary vector g, project it onto each of the basis vectors, and then add them back up again:

g = {2,6,1,7,11,4,13, 29} ;

(g.u1) u1  + (g.u2) u2  +(g.u3) u3  +(g.u4) u4  +
(g.u5) u5  +(g.u6) u6  +(g.u7) u7  +(g.u8) u8

What happens if you project g onto the normalized Walsh basis set defined by {w1,w2,...}  above, and then add up all 8 components?

(g.w1) w1  + (g.w2) w2  +(g.w3) w3  +(g.w4) w4  +
(g.w5) w5  +(g.w6) w6  +(g.w7) w7  +(g.w8) w8

The projections, g.ui are sometimes called the spectrum of g. This terminology comes from the Fourier basis set used in Fourier analysis. A discrete version of a Fourier basis set is similar to the Walsh set, except that the elements fit a sine wave pattern, and so are not binary-valued.

The orthonormal set of vectors we've defined above  is said to be complete, because any vector in 8-space can be expressed as a linear weighted sum of these basis vectors. The weights are just the projections. If we had only 7 vectors in our set, then we would not be able to express any 8-dimensional vector in terms of this basis set. The seven vector set would be said to be incomplete. A basis set which is orthonormal and complete is very nice from a mathematical point of view. Another bit of terminology is that these seven vectors would not span the 8-dimensional space. But they would span some sub-space, that is of smaller dimension, of the 8-space.
There has been much interest in describing the effective weighting properties of visual neurons in primary visual cortex of higher level mammals (cats, monkeys) in terms of basis vectors. One issue is if the input (e.g. an image) is projected (via a collection of  receptive fields) onto a set of neurons, is information lost? If the set of weights representing the receptive fields of the collection of neurons is complete, then no information is lost.

Linear dependence

What if we had 9 vectors in our basis set used to represent vectors in 8-space? For the u's, it is easy to see that in a sense we have too many, because we could express the 9th in terms of a sum of the others. This set of nine vectors would be said to be linearly dependent. A set of vectors is linearly dependent if one or more of them can be expressed as a linear combination of some of the others.  Sometimes there is an advantage to having an "over-complete" basis set (e.g. more than 8 vectors for 8-space; cf. Simoncelli et al., 1992).

	Theorem: A set of mutually orthogonal vectors is linearly independent. 

However, the implication doesn’t go the other way. It is quite possible to have a linearly independent set of vectors which are not orthogonal to each other.  Imagine 3-space and 3 vectors which do not jointly lie on a plane. This set is linearly independent.

If we have a linearly independent set, say of 8 vectors for our 8-space, then no member can be dropped without a loss in the dimensionality of the space spanned.

It is useful to think about the meaning of linear independence in terms of geometry. A set of three linearly independent vectors can completely span 3-space. So any vector in 3-space can be represented as a weighted sum of these 3. If one of the members in our set of three can be expressed in terms of the other two, the set is not linearly independent and the set only spans a 2-dimensional subspace. That is, the set can only represent vectors which lay on a plane in 3-space. This can be easily seen to be true for the set of u's, but is also true for the set of v's.

Thought exercise

Suppose there are three inputs feeding into three neurons in a simple linear network. If the weight vectors of the three neurons are not linearly independent, do we lose information?

Basic matrix arithmetic

Definition of a matrix: a list of scalar lists

An m × n matrix has m rows, and n columns. Here is a 3x4 matrix:

Clear[W,w];
MatrixForm[
	Table[w[i,j],{i,1,3},{j,1,4}]]

And as we have seen, Mathematica represents a matrix as a list of lists:

%//StandardForm

Remember, the symbol % in Mathematica  stands for the most recent output, %% stands for the second to last output, and so forth.

The matrix can be written in traditional mathematical form using subscripts as:

MatrixForm[
	Table[(Subscript[w, i,j]),{i,1,3},{j,1,4}]]

Adding, subtracting and multiplying by a scalar

As with vectors, matrices are added, subtracted, and multiplied by a scalar component by component:

Clear[A,B,a,b,c,d,e,f,x,y,u,v,w,p,q,r,s,t]
A = {{a,b},{c,d}};
B = {{x,y}, {u,v}};

Let's add A to B:

A+B

Subtracting B from A:

A-B

And if we multiply A by 3:

3 A

Multiplying two matrices

We have already seen how to multiply a vector by a matrix: we replace the ith row of the output vector by the inner product of the ith row of the matrix with the vector.
In order to multiply a matrix A, by another matrix B to get C = AB, we calculate the ijth component of the output matrix by taking the inner product of the ith row of A with the jth column of B:

A.B

Note that AB is not necessarily equal to BA:

B.A

Laws of commutation, association and distribution

Look at the element in the upper left of the matrix BA above--there is no reason, in general, for ax+bu to equal ax + cy. That is, matrix multiplication does not commute.

Apart from commutation for matrix multiplication, the usual laws of commutation, association, and distribution that hold for scalars hold for matrices. Matrix addition and subtraction do commute. Matrix multiplication is associative, so (AB)C = A(BC). The distributive law works too: 
A(B+C) = AB + AC

Non-square matrices

It is not necessary for A and B to be square matrices (i.e. have the same number of rows as columns) to multiply them. But if A is an mxn matrix, then B has to be an nxp matrix in order for AB to make sense. For example, here F is a 3x2 matrix, and G is a 2x4 matrix.

F = {{a,b},{c,d},{e,f}};
G = {{p,q,r,s},{t,u,v,w}};

Dimensions[F]

Dimensions[G]

Because F has 2 columns, and G has 2 rows, it makes sense to multiply G by F:

F.G

However, because the number of columns of G (4)  do not match the number of rows of F (3), G.F is not well-defined, so executing this next cell should give you an error message:

G.F

Inverse of a Matrix

Dividing a matrix by a matrix: the identity matrix & matrix inverses

The matrix corresponding to 1 or unity is the identity matrix. Like 1, the identity matrix is fundamental enough, that Mathematica  provides a special function to generate n-dimensional identity matrices. Here is a 2x2:

IdentityMatrix[2]

It is easy to show that the identity matrix plays the role for matrix arithmetic that the scalar 1 plays for scalar arithmetic.

How can one divide one matrix, say B, by another, say A? We can divide numbers, x by y, by multiplying x times the inverse of y, i.e. 1/y. So to do the equivalent of dividing B by A, we need to find a matrix Q such that when A is multiplied by Q, we get the matrix equivalent of unity, i.e. the identity matrix. Then "B/A" can be achieved by calculating the matrix product: B.Q.

A = {{a,b},{c,d}};

Mathematica  provides a built-in function to compute matrix inverses:

Q = Inverse[A]

We can test to see whether the product of a A and Q is the identity matrix, but Mathematica  won't go through the work of simplifying the algebra in this case, unless we specifically ask it to.

Q.A

Simplify[Q.A]

Here is a simple numerical example:

B = {{1,-1},{3,2}};
R = Inverse[B]

B.R

Badly conditioned matrices

What if one row is a scaled version of another row? Then the rows are not linearly independent. In this case, the inverse is not defined.

B1 = {{1.5,1},{3,2.0}}

Inverse[B1]

Sometimes the rows are almost, but not quite, linearly dependent (because the elements are represented as approximate floating point approximations to the actual values). Mathematica may warn you that the matrix is badly conditioned. Mathematica  may try to find a solution to the inverse , but you should be suspicious of the solution. In general, one has to be careful of badly conditioned matrices.

Let's try finding the inverseof the following matrix:

B2 = {{-2,-1},{4.00000000000001,2.0}};
Inverse[B2]

B2.Inverse[B2]

Were we lucky or not?

Determinant of a matrix

There is a scalar function of a matrix called the determinant. If a matrix has an inverse, then its determinant is non-zero. Does B2 have an inverse?

Det[B2]

Why do we get a zero determinant for B2, but not for B1?

Det[B1]

Matrix transpose

We will use the transpose operation quite a bit in this course.  It interchanges the rows and columns of a matrix:

Clear[X,x];
X = Table[(Subscript[x, i,j]),{i,1,3},{j,1,4}]
%//MatrixForm

You may have noticed that in the Cell menu, you can convert to various cell types.
In StandardForm on the input line, the transpose is written:

Transpose[X]

The output default is TraditionalForm. On the input line, in TraditionalForm, the transpose is written:

(X^T)
%//MatrixForm

What do the outputs, X and XT look like in StandardForm, MatrixForm?

Getting parts of a matrix

We can pull out the ith row of a matrix W in Mathematica by simply writing W[[i]]. For example, the 2nd row of X is:

X[[2]]

What about ith column of a matrix?

X[[;;,3]]

This “slicing” notation is similar, but not identical to what you would do in Python or Matlab.

Or, Transpose[W] [[i]] produces the ith column of matrix W. For example, the 3nd column of X is:

Transpose[X][[3]]

Symmetric matrices

If the transpose of a (square) matrix equals itself, HT = H, H is said to be a symmetric matrix. Symmetric matrices occur quite often in physical systems (e.g. the force on particle i by particle j is equal to the force of j on i). This means that the elements of a symmetric matrix H  look like they are reflected about the diagonal.

H = Table[N[Exp[-Abs[i-j]],1],{i,5},{j,5}]

Neural networks and symmetric connections

Do neural systems have symmetric connections? Real neural networks probably do not in general. Although, when the nature of the processing would not be expected to favor a particular asymmetry, we might expect that there should be symmetric connections on average. This assumption is made when setting up neural network weight matrix for lateral inhibition filtering.

Outer product of two vectors

We've already seen that the inner product of two vectors produces a scalar. The outer product of an input and output vector can be used to model synaptic modification. It is can also be used calculate auto-covariances between patterns. Consider two vectors:

Clear[f1,f2,f3,g1,g2,g3,f,g]
f = {f1,f2,f3};
g = {g1,g2,g3};

The outer product is all of the pairwise products of the elements of f and g arranged in a nice order:

h = Outer[Times,f,g]

The outer product is also written in traditional row and column format as: f gT

Compare with the dot product notation, fT g.

Eigenvectors and eigenvalues

Eigenvectors

An eigenvector, x,  of a matrix, A,  is vector that when you multiply it by A, you get an output vector that points in the same direction as x:
Ax = λx
where λ is a scalar that adjusts the length change of x.
There can't be any more than n distinct eigenvectors for an nxn matrix--and there may be less.

A = {{1,2},{3,4}};

The Mathematica function Eigenvectors[A] returns the eigenvectors of matrix A as the rows of a matrix, which we'll call eig:

eig = Eigenvectors[A]

We can verify that eig[[1]] and A.eig[[1]] lie along the same direction by taking the dot product of the unit vectors pointing in the directions of each:

normalize[x_] := x/Sqrt[x.x];

normalize[eig[[1]]].normalize[A.eig[[1]]];
N[%]

Eigenvalues

The eigenvalues are given by:

Eigenvalues[A]

Eigenvalues and eigenvector elements do not have to be real numbers. They can be complex, that is an element can be the sum of a real and imaginary number. In Mathematica, imaginary numbers are represented by multiples (or fractions) of I, the square root of -1:

Sqrt[-1]
Sqrt[-1]//StandardForm

B = {{1,2},{-3,4}};
Eigenvalues[B]

Use convolution to apply a discrete (∇^2)G lateral inhibitory filter to grapes image; 2) locate zero crossings in filtered image

In this exercise, you will build an edge detector based on a “front-end” model of simple cells. You will use the dot product model of a linear visual neuron to spatially filter an input image. You will do this in two ways, and the results should be very similar: convolution, and through the Fast Fourier Transform. 

You will use the filtered images (one from the convolution method and one from the FFT method) to mark edges.

Here is your input image:

ggrapes=Image[];

grapes = ImageData[ggrapes];
size=Dimensions[grapes][[1]]
hsize=size/2

Your final product should mark the edges, and could look something like this, depending on how you decide to mark the edges:

Generate a size x size matrix that represents a 2D low-pass Gaussian filter with standard deviation that is two pixels wide. Put the peak at {hsize+1, hsize+1}. Call the array gaussfilter. You will use this in the FFT problem below. Use Image[] or ArrayPlot[] to show what gaussfilter looks like. It should look like:

Now extract a 13 x13 kernel from the middle (the central "bump").   Call this matrix gausskernel. You will use this in the convolution problem below. Use Image[] or ArrayPlot[] to show what gaussfilter looks like. It should look like:

Generate a size x size array  that represents a discrete Laplacian operator. Center the filter, so that -20 is at {33,33} (See hint below). Call the matrix laplacefilter. 
Big Hint: It is tricky to get a good approximation to a circularly symmetric operator on a discrete rectangular lattice. So I'll give you the specification: The center value should be -20, the left, right, above and below nearest neighbor values should be 4, and the diagonal elements should be 1.  Make sure that you define the matrix to be numerical (i.e. the numbers should be decimals, not integers).
Use Image[] or ArrayPlot[] to display laplacefilter. It should look like:

Extract the central 3 x 3 submatrix that represents a discrete Laplacian operator. Call it: laplacekernel. 
 (Consistent with the above hint, laplacekernel should be:{{1,4,1},{4,-20,4},{1,4,1}};
 Use Image[] or ArrayPlot[] to display laplacekernel. It should look like:

Use ListConvolve to make a (∇^2)G operator by convolving laplacekernel with gauskernel. Call the operator delsq. Display the result using Image[] or ArrayPlot[delsq]
Use the no overhang default for ListPlot and your resulting delsq should be an 11x11 matrix.

Use ListConvolve[] to produce a spatially filtered image of grapes. Call the array filteredgrapes. There are various ways of handling the overhang of the kernel with the image, such as by padding the original image. This time handle the overhang using cyclic convolution (see the documentation for ListConvolve). This will preserve the original size of the array. Adjust the third argument of ListConvolve to make sure that the resulting filtered image doesn't shift with respect to the original image. Use Image[filteredgrapes] or ArrayPlot[filteredgrapes] to show your results.

Generate an image where 1 marks the location of a zero-crossing in your filtered image of grapes, and has zeros everywhere else. There are a number of ways of doing this, devise one. 
Then superimpose the zerocrossings image on the original grapes image. For example, a simple way to superimpose them is:
                       ArrayPlot[zerocross + grapes]
Or use Image[]. But you might want to do something fancier. Your answer should look something like the plot shown at the top of this notebook.

Use the Fourier transform to apply a discrete (∇^2)G lateral inhibitory filter to grapes image and locate zero crossings in filtered image

Apply the (∇^2)G operator in the spatial frequency domain to the grapes image (also represented in the spatial frequency domain), and display the filtered result in the space domain.
You can do the filtering operation as follows: Multiply the Fourier transforms of the gaussian filter, the Laplacian operator filter, and the image. (Use Fourier[ ] to take the Fourier transforms of each 64x64 matrix, and multiply them together). Then you need InverseFourier[ ] to produce the filtered image in the space domain.
Use Chop[ ] on the final output of the inverse transform to make sure there aren't any tiny residual imaginary numbers. These can arise due to limits in numerical precision.
One can think of the resulting image as a "neural image" representing the outputs of an array of ganglion off-center cells. The displayed filtered image should look like the filtered image of grapes you produced using convolution.

Apply your earlier zero-crossing detector to generate an image where 1 marks the location of a zero-crossing in your filtered image of grapes, and has zeros everywhere else. Then, as before, superimpose the zerocrossings image on the original grapes image. E.g.
                      Image[zerocross + grapes]
Again, make sure the edge map and the grapes image are in register.

This exercise may look long, but it is easy involving using one function at the right spot (see HINT below).  

An early observation in the recordings of neurons in mammalian visual cortex was that  simple cells corresponding to “bar” and “edge” filtering often seemed to be anatomically close. The bar and edge detectors can be modeled by gabor functions whose underlying sinusoidal waves are shifted 90 degrees relative to each other, as with cosine and sine functions. 

Further, modelers showed that if one squared the outputs of the bar and edge and added them up, one had a response that resembled that of complex cells. 

The exercise provides code for the non-linear combination of two linear filters each representing a simple cell that prefers vertical orientations, with one filter being in sine-phase, and the other in cosine-phase. Such a pair is called a quadrature pair in signal processing. When the outputs of each filter are squared (modeling a simple point-wise non-linearity) and then added, the output emphasizes contrast changes irrespective of polarity. Over small spatial scales, the response of such a “detector” is also approximately invariant to the location of a edge, and can be interpreted as a soft logical “OR” of position.

Modify the code below to show the response to an annulus with the opposite contrast. HINT: Use ColorNegate@

picture=ImageData[ColorConvert[Graphics[Annulus[]],"Grayscale"]];
g9= ImageAdjust@Image[picture];

sgabor[x_,y_,fx_,fy_,sig_]:=N[Exp[(-x^2-y^2)/(2 sig*sig)] Sin[2 Pi (fx x+fy y)]];
cgabor[x_,y_,fx_,fy_,sig_]:=N[Exp[(-x^2-y^2)/(2 sig*sig)] Cos[2 Pi (fx x+fy y)]];

fsize = 32;sfilter=Chop[Table[sgabor[(i-fsize/2),(j-fsize/2),0,1/8,4],{i,0,fsize},{j,0,fsize}]];
g10= ImageAdjust@Image[sfilter];

fsize = 32;cfilter=Chop[Table[cgabor[(i-fsize/2),(j-fsize/2),0,1/8,4],{i,0,fsize},{j,0,fsize}]];
g11=ImageAdjust@Image[cfilter];

fspicture = ListConvolve[sfilter, picture];
g13 = ImageAdjust@Image[fspicture];

fcpicture = ListConvolve[cfilter, picture];
g14 = ImageAdjust@Image[fcpicture];

Look for peaks in local contrast energy

ss = Sqrt[fspicture^2 +fcpicture^2];
g15=ArrayPlot[ss,Mesh->False];
g15=ImageAdjust@Image[ss];
GraphicsRow[{g9,g13,g14,g15}]