vb_calc_roi_current

 Calculate ROI-current according to area file

 - Input
   current_file : Current file (.curr.mat)
   area_file : Area file (.area.mat) defining ROIs
   brain_file : Brain file (.brain.mat)
   mode : 1 -> Average (default)
          2 -> First principal component
          3 -> Maximum current

 - Output
   Z : ROI-current (roi x time x trial)
   ix_roi : Index of center vertex in each ROI

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [Z, ix_roi] = vb_calc_roi_current(current_file, area_file, brain_file, mode)
0002 % Calculate ROI-current according to area file
0003 %
0004 % - Input
0005 %   current_file : Current file (.curr.mat)
0006 %   area_file : Area file (.area.mat) defining ROIs
0007 %   brain_file : Brain file (.brain.mat)
0008 %   mode : 1 -> Average (default)
0009 %          2 -> First principal component
0010 %          3 -> Maximum current
0011 %
0012 % - Output
0013 %   Z : ROI-current (roi x time x trial)
0014 %   ix_roi : Index of center vertex in each ROI
0015 %
0016 % Copyright (C) 2011, ATR All Rights Reserved.
0017 % License : New BSD License(see VBMEG_LICENSE.txt)
0018 
0019 % Load data
0020 [~, current] = vb_load_current(current_file);
0021 load(area_file, 'Area')
0022 [V, F, xx, inf_C] = vb_load_cortex(brain_file);
0023 
0024 if ~exist('mode', 'var')
0025     mode = 1;
0026 end
0027 
0028 % Calculate ROI-current
0029 Nroi = length(Area);
0030 [~, Nt, Ntr] = size(current);
0031 Z = zeros(Nroi, Nt, Ntr);
0032 switch mode
0033     case 1% Average
0034         for roi = 1:Nroi
0035             ix = Area{roi}.Iextract;
0036             Z(roi, :, :) = mean(current(ix, :, :), 1);
0037         end
0038     case 2% Principal component
0039         for roi = 1:Nroi
0040             ix = Area{roi}.Iextract;
0041             Nix = length(ix);
0042             c = reshape(current(ix, :, :), Nix, Nt*Ntr)';
0043             [u, s] = svd(c, 0);
0044             us = u(:, 1)*s(1, 1);
0045             mc = mean(current(ix, :, :), 1);
0046             mc = mc(:)';
0047             us = us*sign(mc*us);
0048             Z(roi, :, :) = reshape(us', 1, Nt, Ntr);
0049         end
0050     case 3% Maximum current
0051         for roi = 1:Nroi
0052             ix = Area{roi}.Iextract;
0053             Nix = length(ix);
0054             c = reshape(current(ix, :, :), Nix, Nt*Ntr);
0055             pc = sum(c.^2, 2);
0056             [~, ix2] = max(pc);
0057             Z(roi, :, :) = current(ix(ix2), :, :);
0058         end
0059 end
0060 
0061 % Obtain index of center vertex in each ROI
0062 ix_roi = zeros(Nroi, 1);
0063 for roi = 1:Nroi
0064     ix = Area{roi}.Iextract;
0065     Nix = length(ix);
0066     center = mean(V(ix,:), 1);
0067     distance = sum((V(ix, :)-repmat(center, Nix, 1)).^2, 2);
0068     [~, b] = min(distance);
0069     ix_roi(roi) = ix(b);
0070 end
0071 
0072 % Check ix_roi
0073 if 0
0074     figure
0075     color = hsv(Nroi);
0076     [~, ix] = sortrows(rand(Nroi, 1));
0077     color = color(ix, :);
0078     plot_parm = vb_set_plot_parm;
0079     vb_plot_cortex(plot_parm, V, F, inf_C);
0080     for roi = 1:Nroi
0081         ix = Area{roi}.Iextract;
0082         Nix = length(ix);
0083         for i = 1: Nix
0084             plot3(V(ix(i), 1), V(ix(i), 2), V(ix(i), 3), 'o', 'markersize', 5,...
0085                 'markerfacecolor', color(roi, :), 'markeredgecolor', color(roi, :))
0086         end
0087         plot3(V(ix_roi(roi), 1), V(ix_roi(roi), 2), V(ix_roi(roi), 3), 'o', 'markersize', 7,...
0088             'markerfacecolor', color(roi,:), 'markeredgecolor', [1 1 1])
0089     end
0090     axis off
0091 end
0092 
0093 
0094