Description of dmri_fiber_track_prob_specific

0001 function dmri_fiber_track_prob_specific_pair(...
0002                               parcel_dir, parcel_from_list, parcel_to_list, ...
0003                               mif_file, mask_file, ...
0004                               fa_gz_file, fs_brain_gz_file, ...
0005                               trans_info_dir, dmri_file, ...
0006                               output_dir, process_host, Nworker)
0007 % Stream lines tracking and save connection.
0008 % (parcelN.nii.gz --> parcel_exN.nii.gz)
0009 %
0010 % [Usage]
0011 %     dmri_fiber_track_specific_pair(...
0012 %                             parcel_dir, parcel_from_list, parcel_to_list, ...
0013 %                             mif_file, mask_file, ...
0014 %                             fa_gz_file, trans_info_dir, ...
0015 %                             output_dir, [,process_host][,Nworker]);
0016 %
0017 % [Input]
0018 %           parcel_dir : parcel files in the directory are used.
0019 %                        parcel_dir/parcelN.nii.gz
0020 %     parcel_from_list : calculate path from parcel_from_list(n) to parcel_to_list(n).
0021 %       parcel_to_list : calculate path from parcel_from_list(n) to parcel_to_list(n).
0022 %             mif_file : Fiber orientation density function file(.mif)
0023 %            mask_file : a mask region of interest.
0024 %           fa_gz_file : FA nifti-gz file(.nii.gz) Fractional Anisotropy image.
0025 %     fs_brain_gz_file : Freesurfer nifti-gz file(.nii.gz) Cortex image.
0026 %       trans_info_dir : transform information directory.
0027 %            dmri_file : difusion mri file(.dmri.mat)
0028 %           output_dir : The output directory for fiber tracking result.
0029 %         process_host : [Optional] process hosts name(default: {'localhost'})
0030 %                         stream lines tracking will be executed on specified hosts.
0031 %                        caution: It requires SSH login to the hosts without password.
0032 %                         e.g. {'cbi-node01g', 'cbi-node02g', 'cbi-node03g'};
0033 %              Nworker : [Optional] The number of matlab process which is used for
0034 %                                   tck file processing. (default:8)
0035 %
0036 % [Output]
0037 %    none
0038 %
0039 % [Output files]
0040 %    output_dir/parcelA-B.tck : fiber tracking result from parcel A to parcel B.
0041 %                               A : parcel_from_list(n)
0042 %                               B : parcel_to_list(n)
0043 %
0044 %
0045 %    output_dir/Vmni_connect.mat : Vmni_connection file.
0046 %                                  This is used for visualization brain activity.
0047 %                                  http://www.cns.atr.jp/cbi/?p=596
0048 %
0049 % Copyright (C) 2011, ATR All Rights Reserved.
0050 % License : New BSD License(see VBMEG_LICENSE.txt)
0051 
0052 %
0053 % --- Previous check
0054 %
0055 if nargin < 10
0056     error('Please check input arguments.');
0057 end
0058 if exist(parcel_dir, 'dir') ~= 7
0059     error('Specified parcel directory not found.');
0060 end
0061 if exist(mif_file, 'file') ~= 2
0062     error('Specified mif_file not found.');
0063 end
0064 if exist(mask_file, 'file') ~= 2
0065     error('Specified mask_file not found.');
0066 end
0067 if exist(output_dir, 'dir') ~= 7
0068     mkdir(output_dir);
0069 end
0070 if length(parcel_from_list) ~= length(parcel_to_list)
0071     error('The length of parcel_from_list and parcel_to_list must be same.');
0072 end
0073 if exist(fa_gz_file, 'file') ~= 2
0074     error('Specified fa_gz_file not found.');
0075 end
0076 if exist(fs_brain_gz_file, 'file') ~= 2
0077     error('Specified fs_brain_gz_file not found.');
0078 end
0079 if exist(trans_info_dir, 'dir') ~= 7
0080     error('Specified trans_info_dir not found.');
0081 end
0082 if exist(dmri_file, 'file') ~= 2
0083     error('Specified dmri_file not found.');
0084 end
0085 if ~exist('process_host', 'var') || isempty(process_host)
0086     process_host = {'localhost'};
0087 end
0088 if ischar(process_host)
0089     process_host = {process_host};
0090 end
0091 if ~exist('Nworker', 'var')
0092     Nworker = 8;
0093 end
0094 load(dmri_file, 'brain_file', 'vbmeg_area_ix');
0095 if exist(brain_file, 'file') ~= 2
0096     error('Cannot load brain_file : %s', brain_file);
0097 end
0098 
0099 
0100 dmri_setenv;
0101 check_cmd = which('is_cpu_host_available.sh');
0102 disp('Connecting hosts without password ...');
0103 connection = true;
0104 for k=1:length(process_host)
0105     if strcmpi(process_host{k}, 'localhost')
0106         [a, real_host] = system('hostname');
0107         real_host = real_host(1:end-1); % remove last \n
0108         process_host{k} = real_host;
0109     end
0110     fprintf('%s ...', process_host{k});
0111     ret_msg = '';
0112     [a, ret_msg] = system(['timeout 5 ' check_cmd ' ', process_host{k}]);
0113     
0114     is_ok = true;
0115     if a == 0 || a == 1
0116         cmd_track = ['timeout 5 ssh ' process_host{k} ' "', 'bash -c ''source ' getenv('BASH_PROF') ';streamtrack --help''"'];
0117         [ret_st, ret_msg] = system(cmd_track);
0118         if ret_st ~= 0
0119             is_ok = false;
0120         end
0121     else
0122         is_ok = false;
0123     end
0124     if is_ok
0125         disp('OK.');
0126     else
0127         connection = false;
0128         disp(sprintf('Failed.\n   %s', ret_msg));
0129     end
0130 end
0131 if connection == false
0132     error('Please check the setting of ssh connection.');
0133 end
0134 
0135 %
0136 % --- Main Procedure
0137 %
0138 disp('Now preparing for execution...');
0139 
0140 s = dbstack;
0141 if length(s) > 1 && isequal(s(2).name, mfilename)
0142 else
0143     fprintf('Start Fiber tracking to create Vmni_connect.mat file.\n');
0144 end
0145 
0146 % get absolute output directory path.
0147 pushd = pwd;
0148 cd(output_dir);
0149 output_dir = pwd;
0150 cd(pushd);
0151 
0152 % Working directory
0153 work_dir = tempname(output_dir);
0154 mkdir(work_dir);
0155 process_fiber   = cell(0);
0156 kill_remote_process = true;
0157 
0158 % Define terminate function
0159 function terminate_fcn
0160 
0161     if kill_remote_process
0162         % Kill fiber tracking process on specified hosts
0163         for k=1:length(process_fiber)
0164             process_fiber{k}.destroy();
0165         end
0166     end
0167     % Kill the workers to process fiber tracking result files.
0168     dmri_fiber_tck_file_worker_stop(work_dir);
0169 end
0170 
0171 fprintf('working directory : %s\n', work_dir);
0172 fprintf('(caution: Do not press the ctrl+c when copying the path string.)\n');
0173 
0174 % copy all parcel files into work_dir
0175 parcel_list = unique([parcel_from_list; parcel_to_list]);
0176 for k=1:length(parcel_list)
0177     copyfile(fullfile(parcel_dir, filesep, ['parcel', num2str(parcel_list(k)), '.nii.gz']), work_dir);
0178 end
0179 copyfile(mask_file, work_dir);
0180 
0181 % Unzip mask_file
0182 [p, f, e] = fileparts(mask_file);
0183 mask_nii_gz_file = fullfile(work_dir, [f, e]);
0184 mask_nii_file    = strrep(mask_nii_gz_file, '.nii.gz', '.nii');
0185 cmd = ['gunzip ' mask_nii_gz_file];
0186 system(cmd);
0187 
0188 % Parallel processing
0189 Nhost = length(process_host);
0190 
0191 Npair = length(parcel_from_list);
0192 
0193 script_files    = cell(0, 1);
0194 fids            = cell(0, 1);
0195 Npair_per_host = ceil(Npair / Nhost);
0196 
0197 finished_list_dir = fullfile(work_dir, filesep, 'finished');
0198 mkdir(finished_list_dir);
0199 timeout_sh = which('automatic_shutoff_short.sh');
0200 
0201 
0202 % Create scripts as many as the number of host.
0203 for k=1:Npair
0204     if mod(k, Npair_per_host) == 1 || Npair == 1 || Npair_per_host == 1
0205         Nth_script_file = length(script_files)+1;
0206         script_file     = fullfile(work_dir, filesep, [process_host{Nth_script_file}, '.sh']);
0207 
0208         % create template script
0209         isok = dmri_script_file_create(script_file);
0210         if isok == false, error('Failed to create script.'); end
0211 
0212         % open append mode
0213         fid = fopen(script_file, 'a');
0214         if fid == -1
0215             error('File open failed. processing terminated.');
0216         end
0217 
0218         fids{Nth_script_file}         = fid;
0219         script_files{Nth_script_file} = script_file;
0220     end
0221 
0222 
0223     % Unzip .nii.gz
0224     from = parcel_from_list(k);
0225     to   = parcel_to_list(k);
0226     parcel_a_to_b_str = [num2str(from) '-' num2str(to)];
0227     parcel_n_gz_file1 = fullfile(work_dir, filesep, ['parcel', num2str(from), '.nii.gz']);
0228     parcel_n_gz_file2 = fullfile(work_dir, filesep, ['parcel', num2str(to),   '.nii.gz']);
0229 
0230 
0231     fprintf(fid, 'gunzip %s\n', parcel_n_gz_file1);
0232     fprintf(fid, 'gunzip %s\n', parcel_n_gz_file2);
0233 
0234     % Fiber tracking
0235     tck_file = fullfile(work_dir, filesep, ...
0236                         ['parcel', parcel_a_to_b_str, '.tck']);
0237     parcel_n_nii_file1    = strrep(parcel_n_gz_file1,    '.nii.gz', '.nii');
0238     parcel_n_nii_file2    = strrep(parcel_n_gz_file2,    '.nii.gz', '.nii');
0239     
0240     fprintf(fid, ['echo "Fiber tracking - parcel', ...
0241                   parcel_a_to_b_str, '"\n']);
0242     fprintf(fid, '%s\n', [timeout_sh, ' ', ...
0243                         'streamtrack SD_PROB', ...
0244                         ' '          ,mif_file, ...
0245                         ' -seed '    ,parcel_n_nii_file1, ...
0246                         ' -include ' ,parcel_n_nii_file2, ...
0247                         ' -mask '    ,mask_nii_file, ...
0248                         ' '          ,tck_file, ...
0249                         ' -maxnum 50000000', ...
0250                         ' -num 10', ...
0251                         ' -stop' , ...
0252                         ' -unidirectional', ...
0253                         ' -length 300', ...
0254                         ' -debug']);
0255     fprintf(fid, '%s\n', ['if [ ! -e ' tck_file, ' ]; then']); 
0256     fprintf(fid, '%s\n', ['    echo -n > ' tck_file]);
0257     fprintf(fid, '%s\n', ['fi']);
0258     fprintf(fid, '%s\n', ['mkdir ', fullfile(finished_list_dir, filesep, ['parcel' parcel_a_to_b_str])]);
0259     fprintf(fid, '%s\n', ['if [ ! -d ' work_dir ' ]; then']);
0260     fprintf(fid, 'exit $?;\nfi\n');
0261 end
0262 
0263 
0264 % close script files
0265 for k=1:length(fids), fclose(fids{k}); end
0266 
0267 %
0268 % --- Execution
0269 %
0270 cpu_job_throw = which('cpu_job_throw.sh');
0271 runtime = java.lang.Runtime.getRuntime();
0272 
0273 for k=1:length(script_files)
0274     % Create JOB FILE(one host processes only one job).
0275     JOBS_FILE  = fullfile(work_dir, filesep, [process_host{k}, '_job.txt']);
0276     fid = fopen(JOBS_FILE, 'wt');
0277     if fid == -1, error('Job file write error.'); end
0278     fprintf(fid, '%s\n', script_files{k});
0279     fclose(fid);
0280 
0281     % Create HOSTS FILE
0282     HOSTS_FILE = fullfile(work_dir, filesep, [process_host{k}, '.txt']);
0283     fid = fopen(HOSTS_FILE, 'wt');
0284     if fid == -1, error('hosts file write error.'); end
0285     fprintf(fid, '%s\n', process_host{k});
0286     fclose(fid);
0287 
0288     % Execute
0289     cmd = [cpu_job_throw, ' ', JOBS_FILE, ' ', HOSTS_FILE];
0290     process_fiber{length(process_fiber)+1} = runtime.exec(cmd);
0291     pause(5);
0292     % Wait for creating log directory.
0293     while(1)
0294         Nlog_dir = length(dir(fullfile(work_dir, 'log*')));
0295         if Nlog_dir == k
0296             break;
0297         else
0298             pause(1);
0299         end
0300     end
0301 end
0302 fprintf('Running %d process.\n', Nlog_dir);
0303 
0304 %
0305 % --- Convert TCK files into MAT files.
0306 %
0307 parm_mat_file = fullfile(work_dir, 'parm.mat');
0308 matlab_exe = fullfile(matlabroot, 'bin', 'matlab');
0309 vb_save(parm_mat_file, 'work_dir', 'fa_gz_file', ...
0310                        'fs_brain_gz_file', 'trans_info_dir');
0311 
0312 prog_dir   = fileparts(which('vbmeg.m'));
0313 for n=1:Nworker
0314     worker_script   = fullfile(work_dir, filesep, ['worker', num2str(n), '.m']);
0315     worker_err_file = fullfile(work_dir, filesep, ['worker', num2str(n), '_error.txt']);
0316     fid = fopen(worker_script, 'wt');
0317     if fid == -1, error('Worker file create error.'); end
0318     fprintf(fid, '%s\n', 'try');
0319     fprintf(fid, '%s\n', ['cd(''', prog_dir, ''');']);
0320     fprintf(fid, '%s\n', 'vbmeg;');
0321     fprintf(fid, '%s\n', ['dmri_fiber_tck_file_for_display_worker_start(''', ...
0322                            parm_mat_file, ''', ', num2str(n), ');']);
0323     fprintf(fid, '%s\n', 'catch');
0324     fprintf(fid, '%s\n', ['fid = fopen(''' worker_err_file, ''', ''wt'');']);
0325     fprintf(fid, '%s\n', 'err = lasterror;');
0326     fprintf(fid, '%s\n', 'fprintf(fid, err.message)');
0327     fprintf(fid, '%s\n', 'fclose(fid);');
0328     fprintf(fid, '%s\n', 'end');
0329     fprintf(fid, '%s\n', 'exit;');
0330     
0331     fclose(fid);
0332   
0333     cmd = [matlab_exe, ...
0334            ' -nodisplay -nojvm -nosplash ', ...
0335            ' -singleCompThread ', ...
0336            '-r "cd ', work_dir, '; ' 'worker', num2str(n) ';"'];
0337     [status, result] = system([cmd, '&']);
0338 end
0339 
0340 %
0341 % --- Check the number of created files
0342 %
0343 
0344 % Register clean up function(onCleanup is MATLAB builtin function)
0345 onCleanup(@terminate_fcn);
0346 
0347 mat_files_pre = cell(0);
0348 tic;
0349 err = false;
0350 while(1)
0351     log_info_list = dmri_fiber_log_info_get(work_dir);
0352     
0353     % Check result file and Exit code
0354     for l=1:length(log_info_list)
0355         result_file = log_info_list(l).result_file;
0356         err_file    = log_info_list(l).err_file;
0357         if ~isempty(result_file)
0358             [tmp, msg] = system(['cat ' result_file, '| grep ExitCode=']);
0359             if isempty(msg), continue; end
0360             [parse] = textscan(msg, '%s', 'delimiter', '=,');
0361             host      = parse{1}{4};
0362             exit_code = str2double(parse{1}{6});
0363             if exit_code ~= 0
0364                 fprintf('errorlog(host:%s):\n%s\n', host, err_file);
0365                 err = true;
0366             end
0367         end
0368         if l == length(log_info_list) && err
0369             return;
0370         end
0371     end
0372     d = dir(fullfile(work_dir, filesep, 'fs_parcel*.mat'));
0373     mat_files_now = {d.name}';
0374 
0375     % Display progress
0376     diff = setdiff(mat_files_now, mat_files_pre);
0377     if ~isempty(diff)
0378         for j=1:length(diff)
0379             fprintf('Created (%d/%d): %s\n', length(mat_files_pre)+j, Npair, diff{j});
0380         end
0381         mat_files_pre = mat_files_now;
0382     end
0383 
0384     if length(mat_files_now) == Npair && length(dir(finished_list_dir)) == 2 % 2:'.' and '..'
0385         fprintf('Processing (%d/%d)\nDone.\n', Npair, Npair)
0386         break;
0387     end
0388     pause(5);
0389 end
0390 
0391 % Wait for delete of MATLAB workers.
0392 kill_remote_process = false;
0393 pause(2);
0394 dmri_fiber_tck_file_worker_stop(work_dir);
0395 pause(2);
0396 rmdir(finished_list_dir);
0397 
0398 Vtracks = cell(Npair, 1);
0399 retry_parcel_from_list = [];
0400 retry_parcel_to_list   = [];
0401 empty_path_ix          = [];
0402 for k=1:Npair
0403     from = parcel_from_list(k);
0404     to   = parcel_to_list(k);
0405     path_file = fullfile(work_dir, ['fs_parcel' num2str(from) '-', num2str(to), '.mat']);
0406     load(path_file, 'vtracks');
0407     Vtracks{k} = vtracks;
0408     if isempty(vtracks)
0409         % reverse tracking
0410         retry_parcel_from_list = [retry_parcel_from_list; to];
0411         retry_parcel_to_list   = [retry_parcel_to_list; from];
0412         empty_path_ix = [empty_path_ix; k];
0413     end
0414 end
0415 
0416 if ~isempty(empty_path_ix)
0417     % remove empty path
0418     parcel_from_list(empty_path_ix) = [];
0419     parcel_to_list(empty_path_ix)   = [];
0420     Vtracks(empty_path_ix)          = [];
0421 
0422 %     % retry B-A fiber tracking
0423 %     disp('reverse fiber tracking...');
0424 %     for k=1:length(retry_parcel_from_list)
0425 %         disp(['parcel' num2str(retry_parcel_from_list(k)), '-' ...
0426 %               num2str(retry_parcel_to_list(k))]);
0427 %     end
0428 %     output_dir2 = tempname(output_dir);
0429 %     dmri_fiber_track_prob_specific_pair(...
0430 %                               parcel_dir, retry_parcel_from_list, retry_parcel_to_list, ...
0431 %                               mif_file, mask_file, ...
0432 %                               fa_gz_file, fs_brain_gz_file, ...
0433 %                               trans_info_dir, dmri_file, ...
0434 %                               output_dir2, process_host, Nworker);
0435     retry_parcel_from_list = [];
0436     retry_parcel_to_list   = [];
0437     empty_path_ix = [];
0438 end
0439 
0440 
0441 %
0442 % --- Save Vmni_connect.mat
0443 %
0444 V    = vb_load_cortex(brain_file, 'subj');
0445 Vmni = vb_load_cortex(brain_file, 'MNI');
0446 
0447 Index_from = vbmeg_area_ix(parcel_from_list);
0448 Index_to   = vbmeg_area_ix(parcel_to_list);
0449 V_from    = V(Index_from, :);
0450 V_to      = V(Index_to, :);
0451 Vmni_from = Vmni(Index_from, :);
0452 Vmni_to   = Vmni(Index_to, :);
0453 
0454 if ~isempty(retry_parcel_from_list)
0455   rVmni_connect_file = fullfile(output_dir2, 'Vmni_connect.mat');
0456   r = load(rVmni_connect_file);
0457 
0458   for k=1:length(r.Vtracks)
0459       if ~isempty(r.Vtracks{k})
0460           % save B-A result as A-B
0461           Index_from = [Index_from; r.Index_to(k)];
0462           Index_to   = [Index_to; r.Index_from(k)];
0463           V_from     = [V_from; r.V_to(k)];
0464           V_to       = [V_to;   r.V_from(k)];
0465           Vmni_from  = [Vmni_from; r.Vmni_to(k)];
0466           Vmni_to    = [Vmni_to;   r.Vmni_from(k)];
0467           Vtracks    = [Vtracks, r.Vtracks(k)];
0468       else
0469           disp(['No path : parcel', num2str(r.Index_to(k)) '-' num2str(r.Index_from(k)) '.']);
0470       end
0471   end
0472   rmdir(output_dir2, 's');
0473 end
0474   
0475 
0476 % convert Vtracks to Vtracks_mni
0477 parm.brainfile  = brain_file;
0478 
0479 [Vc, Nc]     = cell_to_mat(Vtracks);
0480 Vtracks_mni  = trans_subj_mni_coord(Vc,parm);
0481 Vtracks_mni  = mat_to_cell(Vtracks_mni, Nc);
0482 
0483 % Save Vmni_connect.mat
0484 Vmni_connect_file = fullfile(output_dir, 'Vmni_connect.mat');
0485 vb_save(Vmni_connect_file, 'Index_from', 'Index_to', ...
0486                            'V_from', 'V_to', ...
0487                            'Vmni_from', 'Vmni_to', ...
0488                            'Vtracks', 'Vtracks_mni');
0489 s = dbstack;
0490 if length(s) > 1 && isequal(s(2).name, mfilename)
0491 else
0492     fprintf('Created Vmni_connect file : %s\n', Vmni_connect_file);
0493 end
0494 
0495 % clean up working directory.
0496 rmdir(work_dir, 's');
0497 
0498 toc;
0499 end
dmri_fiber_track_prob_specific_pair

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
