Xu Cui

Here are some parameters in matlab lighting. I first use patch command to plot the surface, then use camlight to shed some light on the surface.

Lighting in MatLab

In your AIR application, if you overwrite the application menu, you might find that the keyboard shortcuts (e.g. Cmd-C, V) doesn’t work under Mac. It seems the “Edit” menu in the default application is essential to these functions. So you should at least keep the “Edit” menu if you need to use your own menu.

SVM is mostly commonly used for binary classifications. But one branch of SVM, SVM regression or SVR, is able to fit a continuous function to data. This is particularly useful when the predicted variable is continuous. Here I tried some very simple cases using libsvm matlab package:

1. Feature 1D, use 1st half to train, 2nd half to test. The fitting is pretty good.

linear 1D

2. Still 1D, but apparently the data is nonlinear. So I use nonlinear SVR (radial basis). The fitting is good.

nonlinear 1D

3. What if we have a lot of dimensions? Here I tried feature space with up to 100 dimensions and calculated the correlation between predicted values and the actual values. For linear SVR (blue), the number of dimension doesn’t affect the correlation much. (red: nonlinear, blue:linear, same data for both cases)

curse of dimension

One property of SVR I like is that, when two features are similar (i.e. highly correlated), their weights are similar. This is in contrast with “winner take all” property of general linear model (GLM). This property is desired in brain imaging analysis: neighbor voxels have highly correlated signals and you want them to have similar weights.

About performance: If different features have different scales, then normalization of data will improve the speed of libsvm. Also, the cost parameter c also affects the speed. The larger c is, the slower libsvm is. For the simulated data I used, the parameters don’t affect the accuracy.

MatLab code: test_svr.m

libsvm: http://www.csie.ntu.edu.tw/~cjlin/libsvm/#matlab

options:

-s svm_type : set type of SVM (default 0)

	0 -- C-SVC

	1 -- nu-SVC

	2 -- one-class SVM

	3 -- epsilon-SVR

	4 -- nu-SVR

-t kernel_type : set type of kernel function (default 2)

	0 -- linear: u'*v

	1 -- polynomial: (gamma*u'*v + coef0)^degree

	2 -- radial basis function: exp(-gamma*|u-v|^2)

	3 -- sigmoid: tanh(gamma*u'*v + coef0)

-d degree : set degree in kernel function (default 3)

-g gamma : set gamma in kernel function (default 1/num_features)

-r coef0 : set coef0 in kernel function (default 0)

-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)

-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)

-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)

-m cachesize : set cache memory size in MB (default 100)

-e epsilon : set tolerance of termination criterion (default 0.001)

-h shrinking: whether to use the shrinking heuristics, 0 or 1 (default 1)

-b probability_estimates: whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)

-wi weight: set the parameter C of class i to weight*C, for C-SVC (default 1)

The k in the -g option means the number of attributes in the input data.

>

This is an example of brain activation plotted on surface. In many circumstances surface view is much more straightforward than a slice view. Here is how I created such a plot using MatLab and SPM.

Environment and Tools:

1. Windows XP
2. MatLab (v7.6)
3. SPM (v8)

Step 1. Generate “faces” and “vertices” of the surface using SPM

If you prefer to use SPM’s interface, launch SPM, then click “Render …” box and click “Xtract surface”. You are asked to input an image file. Select an image file and press Done. Select “Save Extracted Surface” and a few seconds later you will find a file called “surf_???.mat” saved in the current directory. Load this mat file and you will find two variables, faces and vertices. These two variables will be used in the future steps.

If you prefer to use command line, you can simple call function

spm_surf(’xxx.img’, threshold)

About the input image: If you want to plot surface of brain, you should use a segmentation tool (such as FSL) to get the gray matter image, and use this image as input to spm_surf. If you want to plot the surface of skull, you can input the original structural image but with careful choice of threshold.

Step 2. Plot the surface

Plot using MatLab’s patch function. Here is a sample code

p=patch(’faces’,faces,’vertices’,vertices, ‘facecolor’, ‘flat’,  ‘edgecolor’, ‘none’, ‘facealpha’, 1);
facecolor = repmat([1 1 1], length(faces), 1);
set(p,’FaceVertexCData’,facecolor);

daspect([1 1 1])
view(3); axis tight
view([50 -40 100])
camlight
lighting gouraud

Step 3. Overlay a functional image
The strategy here is that you first get the coordinates of the activated voxels and the corresponding colors you want to use for each voxel (presumablly dependent on the intensity of each voxel). Then find the vertices and faces that are “close” to the coordinates, and change the face color.

a. get coordinates

v = spm_vol(imagefile); % functional image file
m = spm_read_vols(v);
pos = find(m>=threshold);

[xyzcor(:,1) xyzcor(:,2) xyzcor(:,3)] = ind2sub(size(m), pos);
xyz = cor2mni(xyzcor, v.mat);
intensity = m(pos);

b. convert intensity to color

tmp = hot; % using ‘hot’ color scheme
mx = max(intensity);
mn = min(intensity);
intensity = round((intensity - mn)/(mx - mn) * 63 + 1);
colors = tmp(intensity,:);

c. plot

distanceThreshold = 2;
p=patch(’faces’,faces,’vertices’,vertices, ‘facecolor’, ‘flat’,  ‘edgecolor’, ‘none’, ‘facealpha’, 1);
facecolor = repmat([1 1 1], length(faces), 1);
pos = [];
for ii=1:size(xyz, 1)
tmp = find(abs(vertices(:,1) - xyz(ii, 1)) <= distanceThreshold & abs(vertices(:,2) - xyz(ii, 2)) <= distanceThreshold & abs(vertices(:,3) - xyz(ii, 3)) <= distanceThreshold  );
pos = [tmp];
posf = find(sum(ismember(faces, pos),2)>0);
facecolor(posf,:) = repmat(colors(ii,:), length(posf), 1);
end

set(p,’FaceVertexCData’,facecolor);

daspect([1 1 1])
view(3); axis tight
view([50 -40 100])
camlight
lighting gouraud

We use TCP/UDP/IP Toolbox 2.0.5 to read and write data from/to a TCPIP port. It’s fast and reliable. The version we use is 2.0.5. Below is a matlab sample script showing how to connect to another computer (called ETG-4000) with TCPIP :

%Connect LAN Port
fid=pnet('tcpconnect','172.17.101.1',51027);

if pnet(fid,'status')==-1
  disp('Connection Failed.');
  return;
end;

%Write to port (i.e. Send Command for ETG-4000)
pnet(fid,'printf','++Hello ETG-4000\r\n');

%Read from port (i.e. Receive Command from ETG-4000)
welcomeMessage = pnet(fid,'read',74); 

disp('Press ETG-4000 START Button!');

%//////////////////////// Get ETG-4000 Data ////////////////////////////////
while(1);
            buff1=pnet(fid,'read',4, 'uint8');
            hsize=bread(buff1,'int32');

            if hsize==12;%Data is comming!
                %disp('hsize is 12');
                %///////////////////////// Data Number ///////////////////////////////
                buff2=pnet(fid,'read',4, 'uint8');
                num=bread(buff2,'int32');%Number of Data
                %disp(['num: ' num]);

                %disp(['num of buff2 ' num2str(num)]);
                %//////////////////////////// Data Size:428 ////////////////////////////
                buff3=pnet(fid,'read',4, 'uint8');
                dsize=bread(buff3,'int32');
                %disp(['dsize ' dsize]);
                %/////////////////////////// Hb Data ////////////////////////////////
                for ch=1:52;%Oxy
                    buff4=pnet(fid,'read',4, 'uint8');
                    oxy(ch,num)=bread(buff4,'single');
                %    disp(oxy(ch,num));
                end;
                for ch=1:52;%Deoxy
                    buff5=pnet(fid,'read',4, 'uint8');
                    deo(ch,num)=bread(buff5,'single');
                end;

                %//////////////////////////// Mark ///////////////////////////////
                buff6=pnet(fid,'read',2, 'uint8');
                mark(num)=bread(buff6,'int16');

                %/////////////////////////// Time ////////////////////////////////
                buff7=pnet(fid,'read',10);
                time(num,:)=char(buff7);
                %disp(time(num,:))

          end;
          %(check status) Push ETG-4000 Stop Button
          stat=pnet(fid,'status');
          if(stat == 0)
            break;
          end
end;

pnet('closeall');
clear fid;

bread is a function to convert binary data to a certain data type:

function y = bread(x, type)
y = typecast(x, type); % work for pnet

One important thing is you need specify the 4th parameter (datatype) when you read data in pnet(’read’). The default value is ‘char’ and if your incoming data is float, you will find any value in a byte which is bigger than 127 is reset to 255. Specifying the 4th parameter to ‘uint8′ solves this problem. I spend a day to find this.

Noise removal methods in NIRS can be divided into 4 categories:

1. reducing noise based on its temporal characteristics: The instrument noise is usually in the high frequency band and thus can be removed by band pass filtering. Band pass filtering can also remove low frequency drift. A real-time version of band pass filtering is exponential moving average (EMA, Utsugi 2007).
2. reducing noise based on its spatial characteristics: motion related noise is assumed to be more spatially spread. The “common” component of the signal across multiple channels (e.g. using PCA) can be treated as noise. (Zhang 2005; Wilcox 2005)
3. reducing noise based on its effect on the correlation between oxy- and deoxy-Hb: motion noise will make the correlation between oxy- and deoxy-Hb, which is typically negative, less negative. (Cui 2010) check out
4. measuring noise independently and subtracting it from the signal. (Zhang 2007, 2009)

Band pass filtering or moving average performs pretty well in reducing non-spike like noise and this method is a standard component in my data analysis. For large spike-like motion artifact, correlation based method works fairly well (even in real-time settings). Of course, for offline analysis, one can also remove these large spikes manually (or semi-automatically).

As I have quite some emails (>13000) in gmail, searching more efficiently becomes a necessity. Fortunately gmail offers some advanced search syntax.

http://mail.google.com/support/bin/answer.py?hl=en&answer=7190

1. Search emails with attachment
   fmri meeting has:attachment
2. Search emails with attachment, and the attachment is a pdf file
   fmri meeting filename:pdf
3. Search emails sent from John
   fmri meeting from:john
4. Search emails sent to John
   fmri meeting to:john
5. Search emails sent before 2009/08/01
   fmri meeting before:2009/08/01
6. Search emails sent after 2009/08/01
   fmri meeting after:2009/08/01