A Gentle Introduction to Matlab for Image Processing
Downloads and Preparation
To work through the
tutorial, you will need to go to
Get cnl_mfiles.zip
and Arrays.zip. If you would like
to make your own *.img files into matlab arrays, then you will need spm
installed. A practice set of spm *.img files is also available for download,
cnl_leftles_data.zip. The m-files should
work under unix, but haven't been tested.
You should have
Matlab installed and know how to start it.
This tutorial is designed for Matlab 6.0, which has many new friendly
point and click (Graphical User Interface) features. If you are using an older version of Matlab,
then references to the gui won’t make sense.
To use
“loadimg.m” (in the cnl_mfiles.zip), you
will need to have spm99 installed. If
you don’t have spm, you can use the imgarray and anatarray files (from the
Arrays.zip download), to do everything in the tutorial. These images help you
keep track of left-right flips because they are from a subject with a large
left frontal lesion. Original 2D and 3D structural img files are at: http://merlin.psych.arizona.edu/cgi-bin/wrap/dpat/Public/Imaging/SPM/TestImages/. The entire functional data set in *.img
format is in the cnl_leftles_data.zip download in the matlab area on merlin
(see above URL).
Getting Started
Matlab is a programming language and interface that is
especially good at dealing with matrices of numbers. Read the CNL matlab web page (http://cnl.web.arizona.edu/matlab.htm)
before continuing with this tutorial. When you start Matlab, you will see a
screen divided into several windows:
The large window on the right is the Command
window. Here you can display variables,
type commands etc.
On the lower left, are two tabbed windows, the Current
Directory window which allows you to see and manipulate files in your
current directory, navigate to other directories etc., the Command History
window which keeps track of all the commands you’ve typed as well as copying
those commands to the current command line or rerunning (evaluating) the lines.
On the upper left is the Workspace window. This window displays all current variables
and allows you to perform various operations on them by right clicking and
selecting from the popup window.
On the menu bar, the File menu allows you to create
new files, open files, import data, and set the path (see below). The Help menu gives you access to
extensive online help files.
Conventions: I will show you what to type at the matlab
prompt “>>” in bold black (don’t type the prompt). Comments will adhere to Matlab
conventions: preceded by a “%” and in green.
Matlab will not try to interpret anything on the line after a percent
sign.
|
Setting the Path
Unzip cnl_mfiles.zip in a directory where you want to
work. Matlab finds m-files by looking
in all the directories in its “path”. To tell Matlab where to look for these
new m-files, click file->set path.
You should see the “Set path” window. Click “Add Folder” or “Add with
subfolders” (depending on how many levels of directories you want Matlab to
search for m-files). Browse to the directory you want to add. Select it and
click “Ok” .
|
 |
Now try each of the following commands:
>>
testabc % Should
create 3 variables, testa, testb and testc (if the path was set correctly)
>>
testa % Look at
a variable by typing its name
>>
who % Tells you
what your variables are (same as looking in the workspace window)
>>
whos % Tells you
more about your variables
>> what % Enumerate
matlab files in current directory (make sure you are a directory where there
are some matlab files).
>> ls % Lists files
in current directory (like looking at the Current Directory window)
>>
save testy % Saves current variables in a *.mat file called testy in the
current directory
>>
load testy % Looks for testy.mat and loads its variables into the
workspace.
>>
clear testb % Removes the specified variable from the workspace
>>
clear % Removes all loaded variables from
the workspace
>>
delete testy.mat %removes the
saved testy.mat from the computer
M-files
Use File->Open (Browse) to find and display m-files. As
you progress, m-files should make more sense.
Within m-files, select and right-click for a menu: For keyword help choose “help selection”. To
run parts of the program and see what they do, choose “evaluate
selection”. When you get to
"Control Structures", use File->New->M-file to create your own
blank m-file.
Vectors and Arrays
Matlab works with
rectangles of numbers: vectors, and matrices (arrays). A vector is a 1-D matrix or array, a matrix
or array can have any number of dimensions (Figure 1 illustrates 1-D, 2-D
and 3-D arrays). Arrays can represent images:
A 2-D matrix is like an image slice. A 3-D matrix is like a set of
slices comprising a brain. A 4-D matrix
is like an fMRI dataset representing a series of 3-D brain images through time.
Figure 1
Creating Vectors and Arrays
You can create the row
vector “a”, of values 1 2 3 4, in several ways:
>> a = [1 2 3 4]
>> a = [1:4]
>> a = (1:4)
>> a = [1:1:4]
>> a = (1:1:4)
Typically [ ] are
used for concatenation, e.g., [1 2 3 4].
The other
statements are not really concatenation statements so either () or [ ] seem to
work. The colon operator is used in several ways: In a=[1:4] and a=(1:4)
it means “1 thru 4”, or “start with
‘1’and count up to ‘4’”.
In a=[1:1:4]
and a=(1:1:4) a value is inserted in the
middle that refers to the increment.
These say “Start with ‘1’, count up by ‘1’ (the default), end at 4. Think of this latter as the basic way to
write the statement, and think of
a=[1:4] and a=(1:4) are a shorthand that works because the default
increment is ‘1’.
Column vectors can be created from scratch:
>> a2 = [1;2;3;4]
Or you can transpose
a row vector to create a column vector:
>> a2 = a’
Either of the above
produces a2:
1
2
3
4
The 2-D matrix in
Fig. 1, can be created by putting a hard return between rows:
>> b = [1
4 7 10
2 5 8 11
3 6 9 12]
or by putting
semicolons between rows:
>> b = [1 4 7 10; 2 5 8 11; 3 6 9 12]
A 3-D matrix
is more difficult:
>> c(2,3,2) = 0
creates a 3-D array
with 2 rows, 2 columns and two pages, but it is filled with zeros. How might you create the 3-D array in Figure
1?:
>> c = [6 3 2 5 1 4
7 6 8 4 9 3]
>> c = reshape(c, 2, 3, 2)
Describing Arrays
>> size(a) % Size as an
ordered list: #rows, #columns, #pages etc.
>> ndims(a) % Number of dimensions
>> numel(a) % The number of elements in the array.
>> length(a) % Equivalent to max(size(a)), returns the largest dimension of
“a”.
Counting and Indexing Elements
For a 2-D matrix, start counting elements at
‘1’ in the upper left corner and count down to the bottom of the first column,
then start counting again at the top of the next column (see the 2-D matrix in
Figure 1). For the 3-D matrix, count the same way, starting on page 1 (a
“page’ in Matlab is like a slice in an image), then going to page 2, etc.
Each array element
has a simple index based on the count (e.g., 1, 2, 3 etc.) AND an index based on array
shape (like a Cartesian Coordinate); e.g., the simple index “1” corresponds to
the 2D index “1,1” or the 3D index “1,1,1”.
“2” corresponds to the 2-D “2,1” etc.
Finally, each
element in an array has a value assigned to it. For our grayscale Mri and fMRI images, this
value represents intensity: black is zero, white is the largest value (256 for
an 8 bit image, 65,536 for a 16 bit image).
To find the value
of an element/voxel, simply enter a valid index for it. Consider the 3-D array,
Figure 1:
>> size(c)
% 2 rows, 2 columns and 2 pages (like 2 slices from an image)
>> c(1) %
Produces the answer (ans) “6”. “6” is the value
>> c(1,1,1) % Same as
above
>> c(2) % Has the value “3”
>> c(2,1,1) % Same as above
>> c(7) % Has the value “7”
>>
c(1,1,2) % Same as
above
Extracting Subarrays (pg 39 of Mastering Matlab)
You can also
extract a subset of elements from an array:
>> a(1:end) % Gets the whole
matrix “a”, beginning with the first element
>> a(3:end) % Gets elements
3 through the last element of “a”
>> a(3:4) % Gets the
third and fourth elements of “a”
>> a(4:-1:1) % Start with 4th
element of “a”: goes backwards in increments of 1:
% End with the first element of “a”.
>> a([3 1]) % Selects the
third and first element from “a”
>> b(:,2) % Gets
elements in all rows [“:” means “all”] for column 2 of “b”
>> b(2,:) % Gets
elements in all columns for row 2 of “b”
The above examples
only show you the output values. To extract the subarray into a variable, you need a variable name
(e.g., “d”) and the assignment operator (= is the assignment operator,
not to be confused with ==, the equal sign), e.g.:
>> d =
a(4:-1:1)
produces “d”: d= 4 3
2 1
Here’s another one,
a bit more complex:
>>d2 = c(1:2,2:3) % Extract elements that are in both the first 2 rows of c
and the second and
% third columns of c to create d:
2 5
1 4
Array Concatenation (pg 42 of “Mastering Matlab”)
>>
e = [a d] % This works if
“a” and “d” have the same number of rows.
>> e = c(:,[2 2 2 2]) % Create e
using all rows “:” of column 2 concatenated 4 times.
The next two
commands accomplish exactly the same thing:
>> e =
c(:,2+zeros(1,4))
>> e = repmat(c(:,2),1,4) %
Replicate c (all rows, column 2; create columns 1 thru 4 of “e”.
e =2 2 2 2
5 5 5 5
|
2 5 2 5
1 4 2 5
>> g = cat(3,e,f) % Concatenate “e“ and “f“ along the 3rd dimension to
create “g“
Array Manipulation
(pg 57 of “Mastering Matlab”)
Clear the workspace
and create “a” as a 2-D matrix:
1 2 3
4 5 6
7 8 9
Think of several
ways to create a fresh copy of “a” for each manipulation.
>> a(:,2) = 9 % For all rows in
a, set elements in column 2 = 9:
1 9 3
4 9 6
7 9 9
>> a(2,2) = 9 % In row 2,
column 2 (their intersection), set elements to “9”:
1 2 3
4 9 6
7 8 9
|
>> b = a(3:-1:1,1:3)
7 8 9
4 5 6
1 2 3
|
>> b = a(1: 1:3,1:3)
1 2 3
4 5 6
7 8 9
|
>> b = a(3: -1:3,1:3)
7 8 9
|
8 9
5 6
2 3
Merging Arrays
a= 1 2
3 b=4 5 4
4 5 6 4 5 4
7 8 9 5 4 5
>> b(3:4,:) = a(1:2,:)
creates
|
4 5 4
1 2 3
|
>> g(1:8) = b(:,2:3)
creates:
5 5 2 5 4 4 3 6
|
>> h(:) = a(:,2:3)
2
5
8
3
6
9
Array Reshaping (pg 63, Mastering Matlab)
For reference, “a”
is: 1
2 3
4
5 6
7
8 9
>>
b = a(:) % Reshape “a” into a column “b”
>>
b = reshape(a,1,9) % Reshape “a”
into a 1x9 vector “b” (a row vector)
>>
b = a’ % Transpose a to create b (transpose is not equivalent to
reshape)
1 4
7
2 5
8
3 6
9
>>
b(2,:) = [ ]
%
Remove row 2 from “b”
1 4
7
3 6
9
>>
a(2,:) = b(2,:) % Replace row 2 of
“a” with row 2 of “b” (with most recent “b”):
1 2
3
3 6
9
7 8
9
>>
a(2,:) = 6 %
Replace all elements of row 2 in “a” with “6’s”.
>>
a(2,:) = [ 6 6 6] %
Same as above, but without relying on scalar expansion of “6”
1 2
3
6 6
6
7 8
9
Substituting Values in Arrays
a= 0 0 0
0 0 0
0 0 0
|
1 0
1
0 0
0
0 0
0
|
2 0
0
0 0
0
2 0
0
Standard Arrays (pg 52, Mastering Matlab)
Matlab can create arrays of all zeros, all ones, or all
random numbers between 0 and 1.
>>
a = zeros(1,4)
%
Creates array of zeros, “a”, with one row and 4 columns.
>>
b = ones(1,4,3) %
Creates array of ones, “b”, with one row, 4 columns, and 3 pages
>>
c = rand(2,3,4,5) %
Creates random number array “c”, 2 rows, 3 columns, 4 pages
% and 5 4th dimensional components
(Values from 0 to 1).
More Array Manipulations
>> d = permute(c,[1,2,4,3]) % Permute the
ordering of dimensions in c (a 4-D matrix), by switching the third and 4th
dimension. This particular permutation would change an array of concatenated
*.img files into a “BRIK” style array,
(i.e., img files are sorted first by physical slice (page), then by time; BRIK
files are sorted by time first and physical slice afterwards.
>> d = find(c==0) % Create a vector
“d” with indices to every 0 valued element in imgarray.
>> d = squeeze(c) %Create “d” which
contains the same values as “c” but gets rid of any extraneous singleton
dimensions (Since “c” doesn’t have any dimensions of “1”, this won’t do
anything.
________________________________________________________________
Sample Mathematical Operations on Arrays
>>
c^2 % c-squared (c can be a matrix).
>>
d = (c.*a) %This is dot
multiplication, it does element by element multiplication of 2 equal
% sized
arrays; in this case, “c” and “a”. Use
dot multiplication to apply a binary
% (1,0) mask
to an image.
________________________________________________________________
Cell Arrays
A cell array is a useful structure to know about if you want
to work in SPM. A cell array can hold
different sized vectors in each cell. In
SPM, you can use the cell array to hold a vector of stimulus onsets for each of
several conditions (e.g., the vector for the first condition is in cell 1. The vector for the second condition is in
cell 2 etc.)
To create a cell array:
>>a{1,1}
= [1 2 4 6 8]
>>a{1,2}
= [ 5 77 89]
>>a{1,3}
= [3 4 5 6 2 1 7 8 9 334]
You now have a cell array, a, that contains 3 row vectors
(this is perfect for SPM stimulus onset times) in 3 cells.
To view a description of the cell array:
>>a
To see the contents of cell 1:
>>a{1,1}
________________________________________________________________
Some Basic Control Structures
For Loop: The for loop is generally a
non-preferred structure in Matlab because usually array operations can
substitute. Note that by default, a for
loop outputs scalar values and overwrites them each time (unless you
concatenate them into an array). Nested
for loops are complicatedMatrix operations produce other matrices by default:
for
x = a %
Where “a” is an array
x
+ 5 % Add 5 to each element. Note that x is overwritten by
each calculation
end
Alternative: a + 5 % Add 5 to each element in “a”. This is the “matrix” solution
___________________________________
If-else-end
If statements come in levels of complexity depending on how
many conditions (1,2 or more):
1) if-end,
2) if-else-end,
3) if-elseif-else-end
1) A simple if statement: if apples > 5;
cost = 20;
end;
2) if-else statement, 1 very simple and one incorporating an
“and” statement:
|
if
apples > 5;
cost = 20;
else;
cost = 25;
end;
|
if
(a == 5) & (c == 2); %note equal sign
b = 2; % assignment here
else;
b = 3;
end;
|
3) A more complex if-elseif-else statement if apples > 15;
cost
= 20;
elseif apples > 12;
cost
= 22;
elseif apples > 10;
cost
= 23;
else;
cost = 25;
end;
Putting it All Together
Here’s a program that takes an input array, calculates a
global mean for that array, and then creates a second array (binmask) that puts
a zero in every element less than the threshold and a one in every other
element. Our real solution in binmask.m is a little more complex, since we want
to fill each time vector with 0’s if it averages less than the threshold and
fill all other time vectors with 1’s.
Nevertheless, this is a good starting point. Look at picker.m (for loop that gets voxels
like this one), rowpicker.m (for loop that gets rows), and historypicker.m (for
loop that gets vectors along the 4th dimension).
Here you see how to request user input from the command
line, use functions like size, mean2 and reshape, employ "for" and
"if-else" control structures, and use a number of matrix
manipulations (matrix multiplication, creating an empty matrix, specifying
subarrays, and concatenating and reshaping arrays)
in_array
= input('Enter the name of loaded array > '); %Request user input
[r,c,nas,ntr] =
size(in_array); %dim1=rows, dim2=columns, dim3=# of anatomical %slices,
dim4=# of temporal repetitions
GlobalMean
= mean2(in_array); %Calculate
mean for 4-D array
Thresh
= .75*GlobalMean; %Create a threshold at 75% of the mean
binmask
= []; %Create empty variable binmask
for
vr = 1:r; %define vr as all rows from 1 to r (for all rows)
for vc = 1:c; %for all columns
for vnas = 1:nas; %for all slices
for vntr = 1:ntr; %for all timepoints
a =
in_array(vr,vc,vnas,vntr); %set a to each element in the array
if a > Thresh; %if a is > Thresh
binmask = [ binmask 1 ]; % set the
mask element to 1 and concatenate it to the %existing binmask
else %otherwise
binmask = [
binmask 0 ]; % mask element=0 and
is concatenated onto binmask
end;
end;
end;
end;
end;
binmask = reshape(binmask,r,c,nas,ntr); %reshape
output to be the same shape
%
as in_array
Viewing your Matrix
Providing you have spm and the program loadimg.m, you can
load a 4D array of data into the main workspace and play with it. You will need a 4D data set saved in *.img
format. In these examples, we use imgarray (from Array.zip). imgarray is a 4-D
data set from a patient with a large left frontal lesion. The array is 64x64x17x80. To view any “slice”
of the data in a 4D matrix, use surface or surf:
>>surface(imgarray(:,:,7,1) %
View the slice that includes all rows and all columns of the seventh physical
slice at time 1.
>>surf(imgarray(:,:,7,1) %
View the slice that includes all rows and all columns of the seventh physical
slice at time 1. The difference is that it comes up at a different angle by
default and that surf is a “high level” function that allows the user more
control over display properties.
Note that for both “surface” and “surf” plots, you can
rotate the plot and see it displayed as a topo map. You can also annotate it, display or change
the color bar, etc. If you imported your image from Afni with AFNItoANALYZE,
and your original BRIK was in radiological orientation (as it generally would
be in Afni), then viewing the array with surf or surface will show you an image
in radiological orientation. However,
that same image, viewed with imshow (see below) will be switched into
neurological orientation. In addition, for both imshow and surface, the z axis
(dim 3) will be flipped relative to the Afni slices (e.g., in a 12 slice image, slice 0 in Afni will be
slice 12 in Matlab, Afni slice 1 will be Matlab slice 11 etc.).
>>plot(imgarray(:,32,4,2)) % Make an
intensity line plot of all columns at the point where they intersect row 32 of
the fourth physical slice at the second timepoint (in effect, plot row 32 of
slice 4, timepoint 2).
>>plot(imgarray(32,:,4,2) % Make an
intensity line plot of all rows at the point where they intersect column 32 of
the fourth physical slice at the second timepoint (in effect, plot column 32 of
slice 4, timepoint 2).
>>a=(imgarray(:,:,:,1)
%
Extracts all values (a 3D brain volume) at timepoint 1.
>>slice(a,32,32,10) %Creates
figure of 3 planes from "a" intersecting at x(row)=32, y(column)=32
and z=10.
Right click any variable in the Matlab workspace to find 2D
and 3D graphing options. Keep in mind
that once you have displayed a figure window,
it may simply minimize and you won’t see the results of your graphing
unless you go click the figure window in the start bar or request a new figure
window:
>>figure
Basic Statistics on Arrays
>>mean2(imgarray) %returns the
mean value for an array of any number of dimensions. Requires the image processing toolbox.
At least two other statistics for n-dimensional arrays, corr2
and std2 that come with the image processing toolbox
>>mean(mean(mean(mean(randarray)))) % same result
as above for a 4-D array (4 levels of nesting for the 4D array), but doesn’t
require the image processing toolbox.
You can use the same technique for the standard functions std
(the standard deviation of the array), min (the minimum value in the
array), and max (the maximum value in the array). Normally, these simple
statistics work on 1-D vectors, but if you nest them for the number of
dimensions, as in the above example, they can be made to work on any sized
array. Check out zscore.m in cnl_mfiles.
Fun with Guide
You can create graphical user interfaces with guide. cnlimg
is one such interface.
>>cnlimg % Brings up the figure window with active buttons.
>>guide
cnlimg % Brings
up cnlimg in edit mode.
In edit mode, you can right click any button and choose
“Inspect properties”. Simply set the
Callback property to be the name of the m-file that the button should run when
pressed.
You may wish to set “tooltip string” to be the explanatory
help that comes up when someone mouses over the button.
>>guide % Brings up the
gui development environment so you can make your own interface.
Fun with the Image Processing Toolbox
If you have the optional image processing toolbox, then you
can do some other nice things.
>>ver %This
will tell you the version of matlab you have and which toolboxes you have
The image processing toolbox will not automatically treat
the results of running loadimg (which requires spm) as images, because it
expects images to contain certain structures.
For example, the image type we are primarily interested in, intensity
images, must have values between 0 and 1. imgarray.mat in Arrays.zip does not
contain the correct values. You can use
mat2gray to convert it:
>>imgarray2
= mat2gray(imgarray);
A slice out of imgarray2 can now be directly displayed as a
proper intensity image:
>>imshow
(imgarray2(:,:,7,3)); %A figure window should come up. If it does not, type:
>>figure, imshow (imgarray2(:,:,7,3));
You can get imshow to interpret your array values and
display a slice without converting with mat2gray:
>>figure, imshow (imgarray2(:,:,7,3), []); %
add an empty matrix [] to the statement.
But if you are interested in playing with other image
processing functions, you will run into trouble unless you do the mat2gray
conversion. anatarray is a 256x256x17 mri structural image, already converted
with mat2gray. Note that the simple
mat2gray conversion produced an image
that was way too dark. By default, the
min value and max value in the array are used as the limits (In this case 0 and
1454). However, a histogram of one of
the slices revealed that most values in the image were less than 300 and only a
few ranged between 300 and 750 on a typical slice. Explicitly setting a tighter range, based on
the histogram, brightened up the image considerably:
>>anatarray2
= mat2gray(anatarray, [0, 500]);
Other ranges are possible. 1000 looked better than 1454, but
I liked 500 best. You should be very
reluctant to crop the upper part of the range for a functional image, since
that range provides vital information.
However, cropping the range for a structural should be less crucial,
especially if it improves your ability to resolve features by eye.
Now that you have a figure window up, displaying a slice,
type:
>>roipoly
% This allows you to draw a shape on a
displayed image which it automatically makes into a binary mask, so if you draw
around the brain, then an image will be produced in which everything outside
the brain is black (0) and everything inside is white (1).
>>impixel % This
allows you to select points (a little star will appear at each point) on a
displayed image. When you hit enter, the
intensity values are returned for each star (they are presented in RGB format,
for an intensity image, this means you'll get a row of 3 equivalent values
R=G=B for each star). Similar to picker.m.
>>improfile %
This allows you to draw a line on a displayed image, or series of line segments
on the image and get back an intensity profile.
A complex set of line segments results in a 3-D histogram. Similar to rowpicker.m and historypicker.m.
>>imcontour (imgarray(:,:,7,3)) % This draws an edge contour of the slice, even if the
slice is not from a mat2gray converted array.
>>imhist (imgarray(:,:,7,3)) % Creates an image histogram, seems to be finicky about the
range of values (It may require a true intensity image converted with
mat2gray). "hist" will work on
any slice of an array, whether or not it is in the "image" range.
|
Montage creates a
montage of multiple images. Mov creates a movie. Both montage and mov expect
the multiple images to be arranged on the 4th dimension. montage requires that the image be
converted with mat2gray first. mov
require that the image be converted from an intensity image to an indexed
image. The following steps should work
to create a montage and movie respectively:
>>anatmontage = reshape(anatarray,
[256,256,1,17]); % Create anatmontage by
reshaping anatarray so that it has an extra singleton dimension and the
slices are arranged on the 4th dimension. Remember anatarray has already been
converted to a true intensity image with mat2gray.
>>figure, montage(anatmontage) %This should bring up a figure window and display the
montage.
|
 |
To create a movie:
>>[ind_anat,
gray] = gray2ind(anatmontage, 256); %This converts the
image from intensity to indexed. Indexed images havetwo parts, an image and a
colormap. In this example,
"ind_anat" is the indexed image and "gray" is the colormap.
>>anat_mov=immovie(ind_anat,
gray); % This creates the movie (and should play it once)
>>movie(anat_mov) %
This plays the movie
Summary
This should get you started.
Look on-line for m-files, and look at the m-files you have. The ones you’ve downloaded with this tutorial
are pretty straightforward (load imgarray or anatarray and then type the name
of the m-file you want to run). Look at
simple ones like mkmatrixrand.m, and zscore.m.
See corall9.m to see a switch-case statement. Let me (dpat@u.arizona.edu) know when you have
suggestions or find errors. Thanks,
Dianne.
Functions Covered
(*=Image Processing Toolbox)
cat
clear
corr2 *
delete
figure
find
for
gray2ind *
guide
hist
if-else
imcontour *
imhist *
immovie *
impixel *
improfile *
imshow *
length
load
ls
mat2gray *
max
mean
mean2 *
min
montage *
movie *
ndims
numel
ones
permute
plot
rand
repmat
reshape
roipoly *
save
size
squeeze
std
std2 *
surf
surface
ver
what
who
whos
zeros
No comments:
Post a Comment