## Tagged Region Analysis with MATLAB

MATLAB is very well suited for analysis of data collected through the
use of the `display` program. This demo page gives a short
overview of some of the analysis that can be performed on this data.
This is a demonstration of analyzing paint (tag) data with MATLAB.
First, we must select the study to analyze. This can be done with the
EMMA `openimage` function, which returns a handle to a study:

h=openimage('/data/disk2/wolforth/minc/farduharson_leslie_tal.mnc');

We can also load the tags from the tag file (generated by
`display`).
tags = loadtagfile('/data/disk2/wolforth/minc/farduharson_right44.tag');

The point coordinates read from the tag file are expressed in the same
space in which the volume they are derived from is expressed.
Therefore, if the volume is stored in Talairach space (as it is in
this case), the coordinates in the tag file will be in Talairach
space. In order to identify which voxels are tagged, we must convert
from the world space of the volume (Talairach space in this case) to
voxel space. This can be done with the `world2voxel` function:
[x,y,z] = world2voxel(h,tags(:,1),tags(:,2),tags(:,3));

Since the volume we are using is stored as transverse slices, we can
find out which slices are interesting to us simply by examining the
`z` coordinate:
>> min(z)
ans =
44
>> max(z)
ans =
85

By examining the `z` coordinate, we see that slice 65 falls
among the slices that have points painted on them. We can load this
image and display it using the EMMA `getimages` and
`viewimage` functions:
MRI = getimages(h,65);
viewimage(MRI);

Now, let's make a mask where the pixels that were painted are set to
some non-zero value, and the non-painted pixels are set to zero. We
start by initializing the mask to all zero:
mask = zeros(256,256);

Next, find all of the points that lie within slice 65 (the slice that
we are interested in):
index=find(z==65);
xi=x(index);
yi=y(index);

Finally, set the voxels in the mask that are painted to some
convenient value:
for i=1:length(xi);
mask(xi(i),yi(i)) = max(MRI)+500;
end

Now that we have the mask, we can view the painted image superimposed
on the underlying structure. Note that we must reshape the mask
before adding it to the variable `MRI` since `MRI` is a
column vector image, and the mask is a square matrix.
viewimage (reshape(mask,length(MRI),1) + MRI);

We can also get a histogram of the intensity of the tagged pixels:

[no,xo] = gettaggedhist(h,tags);

This histogram can be viewed with the MATLAB bar function. We exclude
the first bin since it contains all of the zeros from the mask:
figure;
bar (xo(2:length(xo)), no(2:length(xo)));

We can also generate a histogram for the pixel intensity over the
whole image (this takes some time). Here we specify that we want a
histogram with 50 bins:

[n1,x1] = getvolumehist (h,50);

This histogram can also be viewed with the MATLAB `bar`
function:
figure;
bar (x1, n1);

This page was created by Mark Wolforth
(wolforth@pet.mni.mcgill.ca)