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 = 85By 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; endNow 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);