Description of dmri_cortex

0001 function dmri_cortex_parcel(brain_file, fs_info_file, num_vertex, parcel_file)
0002 % Get indices of parcellated brain
0003 % [Usage]
0004 %    dmri_cortex_parcel(brain_file, fs_info_file, num_vertex, parcel_file);
0005 %
0006 % [Input]
0007 %      brain_file : VBMEG BRAIN-MAT file(.brain.mat)
0008 %    fs_info_file : Freesurfer information file.
0009 %      num_vertex : The number of verticies.
0010 %     parcel_file : Parcel file.
0011 %
0012 % Original code is vb_original_normal_statics.m
0013 % 2012/12/02 M.Fukushima
0014 %
0015 % Copyright (C) 2011, ATR All Rights Reserved.
0016 % License : New BSD License(see VBMEG_LICENSE.txt)
0017 
0018 %
0019 % --- Previous check
0020 %
0021 if nargin ~= 4
0022     error('Please check input arguments.');
0023 end
0024 if exist(brain_file, 'file') ~= 2
0025     error('brain_file not found.');
0026 end
0027 if exist(fs_info_file, 'file') ~= 2
0028     error('fs_info_file not found.');
0029 end
0030 
0031 %
0032 % --- Main Procedure
0033 %
0034 disp('Cortex parcellation..');
0035 start = tic;
0036 dmri_vbmeg_vertex_select(brain_file, fs_info_file, num_vertex, parcel_file);
0037 
0038 % Max radius of neighbor search
0039 Rmax = 0.004;    % 4mm
0040 
0041 %
0042 % --- Load original brain
0043 %
0044 load(fs_info_file);
0045 
0046 brain_parm.FS_left_file  = FS_wm.lh_smoothwm_asc_file;
0047 brain_parm.FS_right_file = FS_wm.rh_smoothwm_asc_file;
0048 
0049 vb_disp('Load original brain')
0050 [V0L,F0L,n0L,V0R,F0R,n0R] = vb_load_orig_brain_surf(brain_parm);
0051 %
0052 % --- SPM cordinate in [m] , Normal vector is outward.
0053 %
0054 vb_disp('Change coordinate to [m] ');
0055 
0056 V0L = V0L/1000;
0057 V0R = V0R/1000;
0058 
0059 %
0060 % --- Set normal direction outward
0061 %     normal vector for disconnected vertex is set to zero
0062 tic
0063 vb_disp_nonl('Set normal direction outward: ')
0064 [F0L, V0L, xxL] = vb_out_normal_surf(F0L,V0L);
0065 
0066 xxsig = sign( sum( n0L .* xxL , 2) );
0067 n0L = n0L .* repmat(xxsig, [1 3]);
0068 
0069 [F0R, V0R, xxR] = vb_out_normal_surf(F0R,V0R);
0070 
0071 xxsig = sign( sum( n0R .* xxR , 2) );
0072 n0R = n0R .* repmat(xxsig, [1 3]);
0073 vb_disp(sprintf('%f[sec]',toc));
0074 %vb_ptime(toc)
0075 
0076 % Original vertex index for cortex
0077 load(brain_file, 'subj');
0078 [tmp1, tmp2, BV_index] = vb_load_cortex_info(brain_file, 'subj');
0079 
0080 NV0L  = size(V0L,1);
0081 indxL = BV_index.Left;
0082 indxR = BV_index.Right;
0083 
0084 if min(indxR) > NV0L,
0085     indxR = indxR - NV0L;
0086 end
0087 
0088 %
0089 % Search neighbor points along cortex sheet
0090 % and calculate normal direction statistics
0091 %
0092 
0093 %
0094 % --- Make neighbor member list of original brain (Left)
0095 %
0096 tic
0097 vb_disp_nonl('Find neighbor index (Left): ')
0098 [xxDL,xxFL]   = vb_next_distance( F0L, V0L ); 
0099 vb_disp(sprintf('%f[sec]',toc));
0100 %vb_ptime(toc)
0101 
0102 vb_disp_nonl('Make neighbor member list of original brain (Left): ')
0103 tic
0104 Llist  = vb_find_near_member(xxFL, xxDL, indxL, Rmax);
0105 vb_disp(sprintf('%f[sec]',toc));
0106 %vb_ptime(toc)
0107 
0108 %
0109 % --- Make neighbor member list of original brain (Right)
0110 %
0111 tic
0112 vb_disp_nonl('Find neighbor index (Right): ')
0113 [xxDR,xxFR]   = vb_next_distance( F0R, V0R ); 
0114 vb_disp(sprintf('%f[sec]',toc));
0115 %vb_ptime(toc)
0116 
0117 vb_disp_nonl('Make neighbor member list of original brain (Right): ')
0118 tic
0119 Rlist  = vb_find_near_member(xxFR, xxDR, indxR, Rmax);
0120 vb_disp(sprintf('%f[sec]',toc));
0121 %vb_ptime(toc)
0122 
0123 %% Whole Brain Parcellation
0124 % Search radious for parcellation (as large as possible)
0125 Rmax = 0.040; % 4cm
0126 
0127 % Load ix_area_mod
0128 load(parcel_file);
0129 
0130 Nroot_l = length(indxL);
0131 tmp = find(vbmeg_area_ix<=Nroot_l);
0132 half = tmp(end);
0133 
0134 % Parcellation
0135 [lh_id] = vb_find_near_member...
0136   (xxFL, xxDL, indxL(vbmeg_area_ix(1:half)), Rmax);
0137 [rh_id] = vb_find_near_member...
0138   (xxFR, xxDR, indxR(vbmeg_area_ix(half+1:end)-Nroot_l), Rmax);
0139 
0140 
0141 %%
0142 lh_cortex_id = cell(length(lh_id), 1);
0143 rh_cortex_id = cell(length(rh_id), 1);
0144 
0145 
0146 %% 1st element contains sub cortex indicies
0147 %lh_cortex_id{1} = FS_wm.lh_subcortex_index;
0148 %rh_cortex_id{1} = FS_wm.rh_subcortex_index;
0149 
0150 % (lh/rh)_cortex_id contains only cortex index(remove subcortex elements)
0151 for n=1:length(lh_id)
0152     lh_cortex_id{n} = lh_id{n};
0153     [tmp, ind] = intersect(lh_id{n}, FS_wm.lh_subcortex_index);
0154     lh_cortex_id{n}(ind) = [];
0155 end
0156 for n=1:length(rh_id)
0157     rh_cortex_id{n} = rh_id{n};
0158     [tmp, ind] = intersect(rh_id{n}, FS_wm.rh_subcortex_index);
0159     rh_cortex_id{n}(ind) = [];
0160 end
0161 
0162 % Save index
0163 save(parcel_file, '-append', 'lh_id', 'rh_id', 'lh_cortex_id', 'rh_cortex_id', 'brain_file');
0164 
0165 % Create membership matrix(parcels to brain model(V))
0166 % ix = find(parcel(1, :) ~= 0); means
0167 [membershipmat] = make_membershipmat(brain_file, parcel_file);
0168 save(parcel_file, '-append', 'membershipmat');
0169 
0170 disp(sprintf('\nSaved:%s', parcel_file));
0171 
0172 % Check parcellation convers all verticies or not.
0173 is_cover_all = dmri_parcels_cover_all_vertices(parcel_file, ...
0174                                             brain_file, ...
0175                                             fs_info_file);
0176                                         
0177 if is_cover_all == false
0178     warning(['Parcellation result may not be correct. ', ...
0179              'There are vertices that are not affiliated ', ...
0180              'with either parcels. Please try to increase ', ...
0181              'parcellation Number(current:' num2str(num_vertex) ').']);
0182 end
0183 toc(start);
dmri_cortex_parcel

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
