batch_create_datamat.m

绝对经典,老外制作的功能强大的matlab实现PLS_TOOBOX

    create_datamat_info.brain_mask_file = brain_mask_file;

    normalize_volume_mean = options.NormalizeVolumeMean;
    create_datamat_info.normalize_volume_mean = normalize_volume_mean;

    curr = pwd;
    if isempty(curr)
       curr = filesep;
    end

    savepwd = curr;
    home = curr;	% pls_data_path;

    %  add path for all subject
    %
    for i = 1:k
        for j = 1:n
            subj_files{i,j} = fullfile(session.dataset_path, subj_files{i,j});
        end
    end

    orient_pattern = [];

    if isempty(orient_pattern)
       [dims, voxel_size, origin] = rri_imginfo(subj_files{1,1});
    else
       dims = getappdata(gcf, 'dims');
       voxel_size = getappdata(gcf, 'voxel_size');
       origin = getappdata(gcf, 'origin');
    end

    dims = [dims(1) dims(2) 1 dims(3)];

    if (use_brain_mask==1) & ~isequal(dims,mask_dims),
       errmsg ='ERROR: Dimensions of the data do not match that of the brain mask!';
       error(errmsg);
    end;

    if (use_brain_mask == 1)			% coords from brain_mask
       coords = find( brain_mask(:) > 0)';
       m = zeros(dims);		
       m(coords) = 1;
       coords = find(m == 1)';
    end;

    num_voxels = prod(dims);

    progress_hdl = rri_progress_status('create', ['Processing "' session_file '"']);

    section1 = n*k/(n*k+10);		% 1st section of progress bar
    factor = 1/(n*k+10);		% factor for the 2nd section

    rri_progress_ui(progress_hdl, '', 0.5*factor);

    %  make tdatamat, which includes non_brain voxels
    %
    tdatamat=[];
    for i = 1:k

        temp=[];

        for j=1:n

            message=['Loading condition ',num2str(i),', subject ',num2str(j),'.'];
            rri_progress_ui(progress_hdl,'',message);

            img = load_nii(subj_files{i,j}, 1);

            v7 = version;
            if str2num(v7(1))<7
               img = reshape(double(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);
            else
               img = reshape(single(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);
            end

            if ~isempty(orient_pattern)
               img = img(orient_pattern);
            end

%            img = reshape(img,[size(img,1)*size(img,2),size(img,4)]);

            temp = [temp; img(:)'];

        end

        tdatamat=[tdatamat;temp];
        rri_progress_ui(progress_hdl, '', ((i-1)*n+j)*factor);

    end

    message = 'Selecting only the brain voxels ...';
    rri_progress_ui(progress_hdl,'',message);

    [dr,dc]=size(tdatamat);
    factor = 0.1/dr;			% factor for the 2nd section
    for i=1:dr
        rri_progress_ui(progress_hdl,'',section1+(i*factor)*(1-section1));
    end

    %  remap data to eliminate non-brain voxels
    %
    datamat = tdatamat(:,coords);

    %  Check zero variance voxels & NaN voxels, and remove them from coords
    %
    check_var = var(datamat);
    bad_coords = find( (check_var==0) | isnan(check_var) );
    coords(bad_coords) = [];
    datamat(:,bad_coords) = [];

    [dr,dc]=size(datamat);		% update dr/dc

    if (normalize_volume_mean == 1)

       % perform whole-brain ratio adjustment
       gmean=mean(datamat,2);		% grand mean for each image

       rri_progress_ui(progress_hdl,'Normalizing datamat',section1+0.2*(1-section1));

       message = 'Normalize datamat with its volume mean ...';
       rri_progress_ui(progress_hdl,'',message);

       factor = 0.8/dc;				% factor for the 2nd section
       checkpoint = floor(dc/10);		% set check point
       check = checkpoint;
       percentage = 10;

       for i=1:dc
           datamat(:,i)=double(datamat(:,i))./gmean;	% normalized on the mean of each img

           if(i==check)
               rri_progress_ui(progress_hdl,'',section1+(0.2+i*factor)*(1-section1));
               message = [num2str(percentage), '% of the volume is normalized.'];
               rri_progress_ui(progress_hdl,'',message);
               check = check + checkpoint;
               percentage = percentage + 10;
           end
       end

    end

    rri_progress_ui(progress_hdl,'',1);
    message = 'Saving to the disk ...';
    rri_progress_ui(progress_hdl,'',message);

%    elapsed_time = toc;
%    disp('Datamat is created ...');

    % save to disk

   datamatfile = fullfile(home, filename);
   datafile = fullfile(home, dataname);

   if(exist(datamatfile,'file')==2)	% datamat file with the same file name exist
      disp(['WARNING: File ',datamatfile,' is overwritten.']);
   end

   if(exist(datafile,'file')==2)  % data file with same filename exist
      disp(['WARNING: File ',datafile,' is overwritten.']);
   end

   create_ver = plsgui_vernum;

    v7 = version;
    if str2num(v7(1))<7
%       datamat = double(datamat);
       singleprecision = 0;
    else
       singleprecision = 1;
    end

   try
          save(datafile,'datamat','create_ver');
   catch
      msg = sprintf('Cannot save STRUCT data to ''%s'' ',datafile);
      error(msg);
   end;

   try
      save(datamatfile,'datafile','coords','behavdata','behavname', ...
		'bad_coords', 'selected_subjects', ...
		'dims','voxel_size','origin','session_file','session_info', ...
		'create_ver','create_datamat_info','singleprecision');
   catch
      msg = sprintf('Cannot save STRUCT datamat to ''%s'' ',datamatfile);
      error(msg);
   end;

   cd(savepwd);
   close(progress_hdl);

   return;					% create_struct_datamat


%---------------------------------------------------------------------------
function create_pet_datamat(session)

   if exist('plslog.m','file')
      plslog('Batch create PET datamat');
   end

   session_info.description = '';
   session_info.pls_data_path = pwd;
   session_info.num_behavior = 0;
   session_info.behavdata = [];
   session_info.behavname = {};
   session_info.datamat_prefix = session.prefix;
   session_info.num_conditions = session.num_cond;
   session_info.condition = session.cond_name;
   session_info.num_subjects = session.num_subj;
   session_info.subject = session.subject;
   session_info.subj_name = session.subj_name;
   session_info.subj_files = session.subj_files;
   session_info.img_ext = session.img_ext;
   session_info.num_subj_init = -1;

   create_ver = plsgui_vernum;
   filename = [session.prefix '_PETsession.mat'];
   session_file = fullfile(session_info.pls_data_path, filename);

   if(exist(session_file, 'file')==2)
      disp(['WARNING: File ',filename,' is overwritten.']);
   end

   try
      save (session_file, 'session_info','create_ver');
   catch
      msg = sprintf('Cannot save session information to ''%s'' ',filename);
      error(msg);
   end;

   %  create datamat
   %  ==============

   if isnumeric(session.brain_region)
      options.UseBrainRegionFile = 0;
      options.BrainRegionFile = [];
      options.Threshold = session.brain_region;
   else
      options.UseBrainRegionFile = 1;
      options.BrainRegionFile = session.brain_region;
      options.Threshold = [];
   end

   options.session_win_hdl = [];
   options.ConsiderAllVoxels = 0;
   options.NormalizeVolumeMean = 1;

   behavdata = [];
   behavname = {};

    pls_data_path = session_info.pls_data_path;
    num_behavior = session_info.num_behavior;
    datamat_prefix = session_info.datamat_prefix;
    num_conditions = session_info.num_conditions;
    condition = session_info.condition;
    num_subjects = session_info.num_subjects;
    subject = session_info.subject;
    subj_name = session_info.subj_name;
    subj_files = session_info.subj_files;

    k = num_conditions;
    n = num_subjects;
    filename = [datamat_prefix, '_PETdatamat.mat'];

    use_brain_mask = options.UseBrainRegionFile;
    brain_mask_file = options.BrainRegionFile;

    brain_mask = [];
    mask_dims = [];

    if use_brain_mask
       brain_mask = load_nii(brain_mask_file, 1);
       brain_mask = reshape(int8(brain_mask.img), [brain_mask.hdr.dime.dim(2:3) 1 brain_mask.hdr.dime.dim(4)]);
       mask_dims = size(brain_mask);

       create_datamat_info.brain_mask_file = brain_mask_file;
       create_datamat_info.brain_coord_thresh = [];
    else
       create_datamat_info.brain_mask_file = '';
       create_datamat_info.brain_coord_thresh = options.Threshold;
    end

    coord_thresh = options.Threshold;
    normalize_volume_mean = options.NormalizeVolumeMean;

    create_datamat_info.normalize_volume_mean = normalize_volume_mean;
    create_datamat_info.consider_all_voxels_as_brain = options.ConsiderAllVoxels;

    curr = pwd;
    if isempty(curr)
       curr = filesep;
    end

    savepwd = curr;
    home = curr;	% pls_data_path;

    %  add path for all subject
    %
    for i = 1:k
        for j = 1:n
            subj_files{i,j} = fullfile(subject{j}, subj_files{i,j});
        end
    end

    orient_pattern = [];

    if isempty(orient_pattern)
       [dims, voxel_size, origin] = rri_imginfo(subj_files{1,1});
    else
       dims = getappdata(gcf, 'dims');
       voxel_size = getappdata(gcf, 'voxel_size');
       origin = getappdata(gcf, 'origin');
    end

    dims = [dims(1) dims(2) 1 dims(3)];

    if (use_brain_mask==1) & ~isequal(dims,mask_dims),
       errmsg ='ERROR: Dimensions of the data do not match that of the brain mask!';
       error(errmsg);
    end;

    if (use_brain_mask == 1)			% coords from brain_mask
       coords = find( brain_mask(:) > 0)';
       m = zeros(dims);		
       m(coords) = 1;
       coords = find(m == 1)';
    end;

    num_voxels = prod(dims);

    if (use_brain_mask == 0)
       coords = zeros(1, num_voxels);		% initial st_coords
    end

    progress_hdl = rri_progress_status('create', ['Processing "' session_file '"']);

    section1 = n*k/(n*k+10);		% 1st section of progress bar
    factor = 1/(n*k+10);		% factor for the 2nd section

    rri_progress_ui(progress_hdl, '', 0.5*factor);

    %  make tdatamat, which includes non_brain voxels
    %
    tdatamat=[];
    for i = 1:k

        temp=[];

        for j=1:n

            message=['Loading condition ',num2str(i),', subject ',num2str(j),'.'];
            rri_progress_ui(progress_hdl,'',message);

            img = load_nii(subj_files{i,j}, 1);

            v7 = version;
            if str2num(v7(1))<7
               img = reshape(double(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);
            else
               img = reshape(single(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);
            end

            if ~isempty(orient_pattern)
               img = img(orient_pattern);
            end

%            img = reshape(img,[size(img,1)*size(img,2),size(img,4)]);

            temp = [temp; img(:)'];

            %  find brain voxel coords, and
            %  accumulated to find common coords for all img
            %
            if (use_brain_mask == 0)
               coords = coords + find_nonbrain_coords(img(:)',coord_thresh,options.ConsiderAllVoxels);
            end

        end

        tdatamat=[tdatamat;temp];
        rri_progress_ui(progress_hdl, '', ((i-1)*n+j)*factor);

    end

    message = 'Selecting only the brain voxels ...';
    rri_progress_ui(progress_hdl,'',message);

    %  determine the coords of the brain region
    %
    if (use_brain_mask == 0)	% coords from thresh by 'find_nonbrain_coords' 
       coords = find(coords == 0);
    end

    [dr,dc]=size(tdatamat);
    factor = 0.1/dr;			% factor for the 2nd section
    for i=1:dr
        rri_progress_ui(progress_hdl,'',section1+(i*factor)*(1-section1));
    end

    %  remap data to eliminate non-brain voxels
    %
    datamat = tdatamat(:,coords);
%    raw_datamat = datamat;

    [dr,dc]=size(datamat);		% update dr/dc

    if (normalize_volume_mean == 1)

       % perform whole-brain ratio adjustment
       gmean=mean(datamat,2);		% grand mean for each image

       rri_progress_ui(progress_hdl,'Normalizing datamat',section1+0.2*(1-section1));

       message = 'Normalize datamat with its volume mean ...';
       rri_progress_ui(progress_hdl,'',message);

       factor = 0.8/dc;				% factor for the 2nd section
       checkpoint = floor(dc/10);		% set check point
       check = checkpoint;
       percentage = 10;

       for i=1:dc
           datamat(:,i)=double(datamat(:,i))./gmean;	% normalized on the mean of each img

           if(i==check)
               rri_progress_ui(progress_hdl,'',section1+(0.2+i*factor)*(1-section1));
               message = [num2str(percentage), '% of the volume is normalized.'];
               rri_progress_ui(progress_hdl,'',message);
               check = check + checkpoint;
               percentage = percentage + 10;
           end
       end

    end

    rri_progress_ui(progress_hdl,'',1);
    message = 'Saving to the disk ...';
    rri_progress_ui(progress_hdl,'',message);

%    elapsed_time = toc;
%    disp('Datamat is created ...');

    % save to disk

   datamatfile = fullfile(home, filename);

   if(exist(datamatfile,'file')==2)	% datamat file with the same file name exist
      disp(['WARNING: File ',datamatfile,' is overwritten.']);
   end

   create_ver = plsgui_vernum;

    v7 = version;
    if str2num(v7(1))<7
%       datamat = double(datamat);
       singleprecision = 0;
    else
       singleprecision = 1;
    end

   try
      save(datamatfile,'datamat','coords','behavdata','behavname', ...
		'dims','voxel_size','origin','session_file','session_info', ...
		'create_ver','create_datamat_info','singleprecision');
   catch
      msg = sprintf('Cannot save PET datamat to ''%s'' ',datamatfile);
      error(msg);
   end;

   cd(savepwd);
   close(progress_hdl);

   return;					% create_pet_datamat