Computational Vision
U. Minn. Psy 5036W

Problem set 1: Mathematica

Expand, read and/or execute each cell of this notebook. Answer all the questions, and then upload the completed notebook to the class Moodle site.

The most recent result of a calculation can be accessed by %, the one before by %%, and so forth.  Try it.

There’s a lot to explore stylistically in Notebooks, but one of the most common things you will do is to select and manipulate the brackets on the far right. These brackets serve to organize text and calculations into a Notebook with outlining features. You can group or ungroup cells for text, graphs, and expressions in various ways to present your calculations. Explore these options under Cell/Grouping in the menu.  You can see the possible cell types under the Style menu.

The arithmetic example seems pretty dumb so far...but you can see that Mathematica's default symbolic handling of arithmetic can make a difference in the next examples:

a)  calculate the square roots of: 12345678987654321 and 12345678987654321.0.

b) Compare (2^.000000000001)^1000000000000 with (2^(1/1000000000000))^1000000000000

If you don't explicitly tell Mathematica that your expressions involve floating point numbers, it will treat the numbers as true integers. And in general, Mathematica will do symbolic processing as a default, and assume you want as general an answer as possible. Because symbolic processing demands considerably more computer resources than numerical processing, this results in annoying consequences, like slow or no responsiveness. So if you are interested in numerical answers, remember to include a decimal point somewhere, or use N[ ] to force a floating point representation.

Front-end stuff

By ending an expression with ; you can suppress the output:

(3/4)/6;
(3 4)/6

Before Mathematica 6, you had to be particularly careful ... without a semi-colon before evaluating, you could end up with a Notebook full of 50,000 "hi"s. It behaves nicely with later versions. Execute the next cell:

a=Table["hi",{50000}]

Built-in functions

Mathematica has a very large library of built-in functions. They all begin with an uppercase letter and the arguments are enclosed by square brackets. Knowing that, you can often guess the form of a function. Or under the Help menu, go to Find Selected Function.

If you know the function, you can enter: ?Log to get more info. And then click on >> to take you to the documentation page.  Try the next cell:

?Plot

If you type two question marks before a function, ??Plot, you'll get more information. Try it.

Defining your own functions

Let's illustrate function definition by building a simple model of a neuron’s input/output function--its visual receptive field. We'll see the justification later.  Suppose a neuron takes a set of inputs, say x1, x2, x3,  (e.g. from the outputs of three other neurons in the network, or they may represent image intensities at 3 points in an image) and produces an output signal, call it y (which may represent the spiking rate of a ganglion cell).

For a  linear model, the output is the weighted sum of the inputs. We'll assume the weights are fixed and given by w1, w2, w3.

w1=-1; w2 = 2; w3=-1;

Clear[y]; 
y[x1_,x2_,x3_]:= w1*x1 + w2*x2 +w3*x3 + 30;

The underscore in x1_ is important because it tells Mathematica that x1 represents a slot for a variable, not an expression. (Clear[ ] is not needed unless you have earlier defined y[] to be something else, intentionally or not. The reason is that Mathematica allows one to build up function definitions incrementally, but more on that later.)

:= vs.  =

Note that when defining a function, we used a colon followed by equals ( := ) instead of just an equals sign (=). When you use an equals sign, the value is calculated and assigned immediately. When there is a colon in front of the equals, the value is calculated only when called on later. So we use := for function definition because we need to define the function for later use and evaluation,  when we may have new values for its arguments.

(A double equals (==) has yet another meaning and is used to represent a symbolic equation which evaluates to True if left and right hand sides are identical, and False if not.)

xv:=RandomInteger[10];
xf=RandomInteger[10];

a) Now evaluate xv and xf three times each. What is the difference between the two definitions?

{xv, xf}
{xv, xf}
{xv, xf}

You can do the same thing by making a list of lists: {{xv,xf},{xv,xf},{xv,xf}} (lists are covered below).

b) Compare repeated evaluations of the above with repeated evaluations of the two expressions below.

y[xv,xv,xv]
y[xf,xf,xf]

Why does our neuron model always output a 30 when the inputs  all have the same value? In general, what family of inputs will all produce a constant output of 30?

Hint: The weights we defined can be thought of as a discrete approximation to a 2nd derivative operator in differential calculus.

Graphics & more function definitions

Visual neurons, and neurons in general can’t fire faster than an upper limit of several hundred spikes per second. And often it takes a sizeable input strength to get a neuron to fire at all--i.e. to get past its threshold. We can model the suppression of small outputs and the clamping of large outputs to a maximum level using a squashing function. Here is an example:

squash[x_] := N[1/(1 + Exp[-x+4])];

Also note that our squashing function was defined with N[ ]. Again, remember that Mathematica tries to keep everything exact as long as possible and thus will try to do symbol manipulation if we don't explicitly tell it that we want numerical representations and calculations.

Let's plot a graph of the squash function for -5<x<10, using the syntax we learned about above:

Plot[squash[x],{x,-5,10}]

We'll use this and similar functions to model the small-signal compression and large signal saturation characteristics of neural outputs. This kind of non-linearity is sometimes called a point-wise non-linearity.

(Since Mathematica version 10, this sigmoidal function is built in: MouseOver)

Even though Exp grows exponentially fast with increasing x (by definition!), squash “keeps a lid” on it.

Plot[]

It can be important to check your definitions like this. One reason is that, as we will see later, Mathematica definitions can be built up line-by-line, incorporating multiple criteria/constraints. And sometimes you might add to a function unwittingly and it appears to misbehave. You can check your definition by asking Mathematica for it.

Mathematica has useful and easy-to-program GUI tools. One of the most useful is  Manipulate[ ]

Using the sliders for λ and b, find a) values corresponding to a linear regime, and sharp threshold regime, 
as shown in the figures.

Check your definition of dsquash using ?

We’ll often want to visualize the receptive field of a visual neuron, or the effective synaptic weights of a model neuron. Some neurons in the primary visual cortex have 2D spatial pattern of weights that are well-modeled by a gaussian windowed sinusoid (called a gabor function):

Turn the above into a Mathematica expression and visualize the function by replacing zzz with it in the routine below:

Manipulate[Plot3D[zzz,{x,-2,2},{y,-2,2},PlotRange->{-2,2}, ImageSize->Tiny],{a,1,6},{b,1,6}]

Your output should look something like this ():

Try using your mouse and mouse button to rotate the figure in your plot. What does the Option or Alt key do?

We will uses lists a lot. Mathematica doesn’t enforce type casting. Their elements can be numbers, symbols, other lists, and a list can have mixed types.

{1.1, "hi", {Pi, 33/2}}

Mostly we will use lists to represent vectors and matrices. You can quickly make a list  using the Table function:

avector = Table[i^2, {i,1,6}]

The syntax avector[[i]]  outputs the (i^th) value of avector. Access the fourth value of avector.

amatrix = Table[i+3 j,{i,1,3},{j,1,3}]

You can view the list in matrix form using:

Table[i+3 j,{i,1,3},{j,1,3}]//MatrixForm

You can access the element in the  (i^th) row and  (j^th) column either by: amatrix[[i,j]] or  amatrix[[i]][[j]]. Use these syntaxes to output the value of the element in the 2nd row and 3rd column.

amatrix[[2,3]]

We will use matrices (a list of lists of pixels) to represent images. Visualize what the circles 2D list below looks like using Image[ ]

circles= Table[N[Cos[Sqrt[i^2 + j^2]]],{i,-30,30,.1},{j,-30,30,.1}];

Image[circles]

For fun, try executing the following lines which adds a phase term to the circles, and then animates the modified function, rings.

rings[phase_]:=Table[N[Cos[phase+Sqrt[(i)^2 + (j)^2]]],{i,-10,10,.1},{j,-10,10,.1}];
Animate[Image[rings[phase]],{phase,0,2 Pi}]

This stimulus produces a classic motion adaptation effect.

Find an image of your choice (online or on your computer), and drag the image into the argument slot of the function Image[ ] . Here’s an example with a picture of Richard Feynman:

Image[];

You can select, copy, and the paste the image into other functions that take an image as an argument (see MouseOver).

Images are often really big. Drag an image into the argument of  ImageDimensions[ ] to check the size of image without extracting the pixel intensities.

Extract the data in list form to get a list with the intensities:

idata = ImageData[];

Check the format and size of the image data using Dimensions[ ]

So it isn’t a simple matrix. The first two entries tell us the number of rows and columns, respectively. Note that although an image might look black and white, it’s pixels have still three values, one for each color channel. This is because it is in color format. Convert it to black and white using ColorConvert[]

grayimagedata = ColorConvert[,"Grayscale"]

Dimensions[ImageData[ColorConvert[,"Grayscale"]]]

When doing calculations on images, it is useful to be able to easily go back and forth between matrix and vector representations of images. Images are usually represented represented as an m x n matrix, where m and n are the number of rows and columns respectively. But in some calculations it is convenient to treat the image as a vector. How do you get back and forth?  Use Flatten[ ] and Partition[ ] to see how they work

m = {{1,2,3},{4,5,6},{7,8,9},{10,11,12}};
%//MatrixForm
v=Flatten[m]
Partition[v,3];
%//MatrixForm

Now turn grayimagedata into a one dimensional list--i.e. a vector using Flatten[ ];

Now turn the flattened one dimensional list back into a matrix using Partition[ ], and then to check to see if you did it right, display the image data using Image[ ]; You’ll have to get the right 2nd argument for Partition.

Michelson contrast is a common measure of contrast for patterns such as sine wave gratings. It is defined as: (((Subscript[L, max])-(Subscript[L, min]))/((Subscript[L, max])+(Subscript[L, min]))), where (Subscript[L, max])and (Subscript[L, min]) are the maximum and minimum intensity values of an image, respectively. Given a sinewave grating defined by:

bmatrix = Table[0.54*N[Cos[Sqrt[i^2 + j^2]]],{i,-6,6},{j,-6,6}]+17,
 
use Mathematica functions to calculate its Michelson contrast. Hint: use Max[ ], Min[ ], and Flatten[ ];

Animal visual systems have distinct channels for light signals that deviate positively or negatively from an average value over some region of space, referred to as “on” and “off” channels. So sometimes it is useful to represent an image in terms of its deviations from a global average. Further, because the visual system is relatively insensitive to mean levels, we can subtract the mean value out. And by normalizing with respect to the mean, we can capture a basic property of vision--that a scene doesn’t look very different with sunglasses on.

For example, given an intensity pattern L(x,y), we can define contrast at a point (x,y) as: contrast(x,y) = ((L(x,y)-(Subscript[L, 0]))/(Subscript[L, 0])), where (Subscript[L, 0]) is the mean value of L(x,y) over the whole image.

grayimagedata[[x,y]] corresponds to Let L(x,y).  Write a function to turn grayimagedata into a matrix called contrast  (= contrast(x,y)). Calculate the minimum, maximum, and standard deviation of contrast. Now do it again for contrast/10 = ((contrast(x,y))/10).

Compare Image[contrast] with ImageAdjust@Image[contrast]

What happens to the contrast matrix if we view grayimagedata through sunglasses that attenuate the luminance by a factor of 10?

Photon emission and absorption are discrete events that can be modeled as a Poisson process. Use Mathematica to define a Poisson distribution p(k) = (((e^(-μ)) (μ^k))/(k!)), where k is the number of photons that might be measured given a mean rate of μ per second. Hint: Mathematica has a built-in function for many distributions, including Poisson distributions (MouseOver).

What is the probability of measuring 10 photons if the mean is 20?

Use RandomVariate[PoissonDistribution[μ], #samples ] to draw 1000 samples from p(k) = (((e^(-20)) (20^k))/(k!)) and plot the results up as a Histogram (using Histogram[<<samples>>,Automatic,”PDF”]. Superimpose the theoretical (population) prediction on top of the sampled distribution (Hint: DiscretePlot[ ] is a good function to use for this). Your result should look something like this:

Check out the Demonstrations center at the Wolfram Mathematica site for some cool examples, such as:

Manipulate[Module[{bluecol,colfunc},
bluecol[n_]:=Blend[{Black,Blue,White},n];
If[usecol, colfunc=bluecol, colfunc=GrayLevel];

Graphics[{colfunc[left],Rectangle[{0,0},{3,3}],colfunc[right],Rectangle[{3,0},{6,3}],colfunc[mid],Rectangle[{1,1},{2,2}],Rectangle[{4,1},{5,2}]}]],{{left,.35,"left square"},.2,.8},{{right,.75,"right square"},.2,.8},{{mid,.5,"middle squares"},.2,.8},{{usecol, True,""},{True->"blues", False->"grays"},ControlPlacement->Bottom}]

"The Simultaneous Contrast Effect" from The Wolfram Demonstrations Project http://demonstrations.wolfram.com/TheSimultaneousContrastEffect/

You don’t need to write anything for this section. Just make sure you understand and can define the following terms: Weber and Michelson contrast (see wikipedia entry for Contrast (vision)), scotopic/photopic, image intensity patterns vs. scenes, generative model, Marr’s levels of analysis.

Make sure you can draw a schematic of the mammalian eye labeling its main anatomical features.