calc_connection

 Compute connection count and mean fiber length
 Input
  dir_name: directory name where voxel images made by 'fsl2ascii' locate
  trackfile: .tck file made by streamtrack in MRtrix [relative path to dir_name]
  ofile: Output file containing mm coodinate of all the voxels in the mask
  Nlabel: Number of label files
 Output
  connect: Label of the tracks
  count: Connection count for each label
  fiber_length: Fiber length of the tracks
  mfiber_length: mean fiber length for each label
  Noverlap: Number of overlapped labels
  track_parm: copy of 'tracks' other than 'tracks.data'

0001 function [connect, count, fiber_length, mfiber_length, Noverlap, track_parm] = ...
0002   calc_connection(dir_name, trackfile, ofile)
0003 % Compute connection count and mean fiber length
0004 % Input
0005 %  dir_name: directory name where voxel images made by 'fsl2ascii' locate
0006 %  trackfile: .tck file made by streamtrack in MRtrix [relative path to dir_name]
0007 %  ofile: Output file containing mm coodinate of all the voxels in the mask
0008 %  Nlabel: Number of label files
0009 % Output
0010 %  connect: Label of the tracks
0011 %  count: Connection count for each label
0012 %  fiber_length: Fiber length of the tracks
0013 %  mfiber_length: mean fiber length for each label
0014 %  Noverlap: Number of overlapped labels
0015 %  track_parm: copy of 'tracks' other than 'tracks.data'
0016 
0017 % 2012/06/15 M.Fukushima
0018 % 2012/06/26 M.Fukushima
0019 % 2014/05/29 rhayashi     Input argument Nlabel is removed.
0020 %                         ofile contains Nlabel.
0021 %
0022 % Copyright (C) 2011, ATR All Rights Reserved.
0023 % License : New BSD License(see VBMEG_LICENSE.txt)
0024 
0025 % Initialize
0026 Ntrack=0;
0027 load(ofile, 'mcoord', 'label', 'Nlabel');
0028 connect = cell(Ntrack,1);
0029 count = zeros(Nlabel,1);
0030 fiber_length = zeros(Ntrack,1);
0031 mfiber_length = zeros(Nlabel,1);
0032 Noverlap = zeros(Ntrack,1);
0033 track_parm = [];
0034 
0035 % Load files
0036 filename = fullfile(dir_name, filesep, trackfile);
0037 d = dir(filename);
0038 if length(d)==1 && d.bytes > 0
0039     tracks = read_mrtrix_tracks(filename);
0040 else
0041     return;
0042 end
0043 
0044 % Number of tracks
0045 Ntrack = length(tracks.data);
0046 
0047 % Extract coordinates of fiber endpoint
0048 endcoord = zeros(Ntrack,3);
0049 if Ntrack ~= 0
0050   for ii = 1:Ntrack
0051     if ~isempty(tracks.data{ii})
0052       endcoord(ii,:) = tracks.data{ii}(end,:);
0053     end
0054   end
0055 end
0056 
0057 % Shift mask coordinate (from the lowerleft corner to the center)
0058 mcoord = [mcoord(:,1)-1 mcoord(:,2)+1 mcoord(:,3)+1];
0059 
0060 % Step size
0061 step_size = str2num(tracks.step_size);
0062 
0063 % Initialize
0064 connect = cell(Ntrack,1);
0065 count = zeros(Nlabel,1);
0066 fiber_length = zeros(Ntrack,1);
0067 mfiber_length = zeros(Nlabel,1);
0068 Noverlap = zeros(Ntrack,1);
0069 
0070 if Ntrack ~= 0
0071   for ii = 1:Ntrack
0072     % sort distance between mask and fiber (endpoint) coordinate
0073     [inda, indb] = sort(sum((mcoord - repmat(endcoord(ii,:),[size(mcoord,1) 1])).^2,2));
0074     % Number of overlaps
0075     Noverlap(ii) = find(diff(inda)~=0,1);
0076     % Labeling the whole tracks
0077     connect{ii} = label(indb(1:Noverlap(ii)));
0078     % Count for each label
0079     count(connect{ii}) = count(connect{ii}) + 1/length(connect{ii});
0080     % Fiber length
0081     fiber_length(ii) = step_size * size(tracks.data{ii},1);
0082     % Sum of fiber length for each label
0083     mfiber_length(connect{ii}) ...
0084       = mfiber_length(connect{ii}) + fiber_length(ii)/length(connect{ii});
0085   end
0086   tmp = count;
0087   tmp(count==0) = 1;
0088   % Mean fibar length for each label
0089   mfiber_length = mfiber_length./tmp; clear tmp;
0090 end
0091 
0092 % Copy of 'tracks' other than 'tracks.data'
0093 track_parm = tracks;
0094 track_parm.data = [];