Cornell University
BioNB 441
Using Images

Introduction.

Images and image analysis has a large role in neurobiology including:

Matlab makes it easy to produce and handle images algorithmically. It fills a different niche from a program like Adobe Photoshop which allows human artistic interaction with an image. You might use Matlab image abilities to process a few hundred microscope images in order to filter, convert to grayscale, and store as a movie. Such a process would be very tedious using an program designed for artistic manipulations.

The examples will use Matlab and the Matlab Image toolbox.

Examples

  1. Moving sinewave grating images are often used as simple visual stimuli. The following program generates a vertical grating, stores it as a movie, then plays it back several times. 
    clf
clear all
set(gcf,'doublebuffer','on');

%Build the grid which determines the image resolution
[x,y] = meshgrid(1:256,1:128);

%initialize the movie
n=40;
mov=moviein(n);

%build the frames
speed=.05;freq=.2;
for i=1:n
   %the frame is just a sinewave grating
   %intensity modulated to--to-botton
   img=128 * (sin(x*freq+speed*i)+1) .* y.^2/128^2;
   image(img);
   %add a frame number
   text(10,10,num2str(i),'color','white')
   %Force real pixels
   truesize(gcf)
   %transfer the image to a movie frame
   mov(i)=getframe; 
end

%clear the image
image(0*x);axis image;
pause(1);
%playback the movie 
%--five times, with reversal
%--at 30 frames/sec
movie(mov,-5,30);
image(0*x);axis image;
    One frame from the animation is shown below.

  2. We can use signal processing techniques to produce an image of a fractal coastline. A random 2D array of complex frequencies is constructed, then filtered by multiplying by a power law of 1/fc. The flitered array of complex frequencies is then transformed back to a space-array, and added to a bias function to make an island or a lake. The array of altitudes is then turned into an rgb image and displayed. 
    % generate a fractal coastline via FFT
clf;
clear all;

%grid size
n = 256;
[x,y] = meshgrid(1:n,1:n);

%Generate a random frequency domain grid
s = sum(100*clock);
rand('seed',s);
freq = (rand(n,n)-.5) + i* (rand(n,n)-.5);  
% Set the DC level to zero
freq(1,1) = 0;	 

% Slope of cutoff filter in freq domain
%less negative => rougher outline 
%-1.5 to -2.5 is reasonable for real coestlines
c = -2.1; 
%
% make a matrix with entries proportional to the freq
dist = sqrt( x.^2 + y.^2 ) ;

% Reduce high freq components by a power law f**c
%to make a fractal surface
filtered = freq .* dist .^c ;
surface = ( real(ifft2(filtered)) ) ; 

%add a bias to make a lake or a mountain
%just change the sign when adding the bias to the surface
bias = 6e-8*sqrt( (x-n/2).^2 + (y-n/2).^2 );
bias = bias-mean(mean(bias));
surface = surface - bias ;

%Convert to a rgb (truecolor) image
surface = cat(3,...
   .8*(surface>0), ...
   .6*(surface>0),...
   surface<=0);

image(surface);
truesize(gcf);
    An example is shown below.

  3. Another signal processing application is generating band-limited, random, moving visual noise to use as a background for visual stimuli. An example animation shows an image bandlimited in time to about 5 frames and bandlimited in space to about 2-10 cycles. The program generates random complex numbers as a basis for moving Fourier synthesised patterns. 
    figure(1)
clf;
set(gcf,'doublebuffer','on')
clear all;

%grid size
n = 128;

%Generate a random frequency domain grid
s = sum(100*clock);
rand('seed',s);
freq = (rand(n,n)-.5) + i* (rand(n,n)-.5);  
%Set the DC level to zero
newfreq(n/2,n/2) = 0;

%low pass parameter--bigger means faster pattern
alpha = .1;
%lowest spatial frequency -- units are cycles/frame width
lowcut = 2;
%highest spatial frequency
hicut = 10;
 
% make a matrix with entries proportional to the freq
%for use in filtering in freq domain
[x,y] = meshgrid(1:n,1:n);
dist = sqrt( (x-n/2).^2 + (y-n/2).^2 ) ;
    
for j=1:120
    %get the next random input
    newfreq = (rand(n,n)-.5) + i* (rand(n,n)-.5);  
    % Set the DC level to zero
    newfreq(n/2,n/2) = 0;	 
    %lowpass filter combines old and new state
    freq = freq *(1-alpha) + newfreq*alpha;
    
    %band-pass filter the spatial frequency
    filtered = freq .* ( (dist>=lowcut) & (dist<=hicut) );
    
    %make the intensity distribution
    surface = abs( (ifft2(fftshift(filtered))) ) ; 
    %and scale it to 0-1 range
    surface = (surface - min(min(surface)))/(max(max(surface)) - min(min(surface)));
    
    %Convert to a rgb (truecolor) image
    dis = cat(3,surface, surface, surface);
    %display the grayscale image
    imshow(dis);
    truesize(gcf);
    colormap(gray);
    mov(j) = getframe;
    %drawnow
end

figure(2)
movie(mov,1,30);
  4. This example uses binary images to compute the overlap of two regions. The regions represent the ranges of two different animal species computed as the convex hull of the positions in which animals are actually seen. The following program computes two animal distrbutions, then steps through the area analysis visually. 
    %make two animal distibutions
dist1x=randn(100,1);
dist1y=randn(100,1);
dist2x=randn(100,1)+2;
dist2y=randn(100,1);

%compute the convex hull of each
c1=convhull(dist1x,dist1y);
c2=convhull(dist2x,dist2y);

%plot the two distributions and 
%their respective convex hulls
figure(1)
clf
hold on
plot(dist1x,dist1y,'r.');
plot(dist1x(c1),dist1y(c1),'r');

plot(dist2x,dist2y,'b.');
plot(dist2x(c2),dist2y(c2),'b');

%will need the image size below
xdist=get(gca,'xlim');
ydist=get(gca,'ylim');

%=============
%fill the convex hull polygons
%convert them to binary images
%then compute the intersection area
%and compute the area ratios

figure(2)
clf
hold off
axis image
fill(dist1x(c1),dist1y(c1),'r');
set(gca,'xlim',xdist);
set(gca,'ylim',ydist);
pause(1)

img=getframe;
d1bw=~im2bw(img.cdata,.9);
imshow (d1bw);
pause(1)

fill(dist2x(c2),dist2y(c2),'r');
set(gca,'xlim',xdist);
set(gca,'ylim',ydist);
img=getframe;

d2bw=~im2bw(img.cdata,.9);
imshow (d2bw);
pause(1)

%trim the image edges because of 
%bogus pixels at the edges
[a b]=size(d1bw)
d1bw=(d1bw(5:a-10,5:b-10));
d2bw=(d2bw(5:a-10,5:b-10));

imshow(d2bw & d1bw)

disp('areas')
bwarea(d2bw & d1bw)
bwarea(d2bw & d1bw)/bwarea(d2bw)
bwarea(d2bw & d1bw)/bwarea(d1bw)
    An example of two disributions and their computed overlap region is shown below.