function [variables] = read_lesion_imgs(parameters, variables)

for ni=1:length(variables.SubjectID)
    fname = fullfile(parameters.lesion_img_folder, [variables.SubjectID{ni}, '.nii']);
    if(~exist(fname))
        warning('can not find lesion image file: %s\n', fname);
        continue;
    end
    vo = spm_vol(fname);
    tmp = spm_read_vols(vo);
    
    tmp(find(isnan(tmp(:))==1)) = 0;
    Ldat(:,:,:,ni) = uint8(tmp);%/norm(tmp(:));
end
variables.vo = vo;
variables.vo.name = 'NULL.nii';
variables.vo.dt = [64,0];

%% get a mask based on overlapping map and given mask image
mask_map = sum(Ldat, 4);
variables.l_idx = find(mask_map >= parameters.lesion_thresh);

% get the voxel index within the generated mask
Mdat = reshape(Ldat, length(mask_map(:)), variables.SubNum).';

variables.lesion_dat = double(Mdat(:,variables.l_idx));

%% Read atlas images
vo = spm_vol(parameters.atlas_aal);
variables.atlas_aal_map = uint8(spm_read_vols(vo));
[variables.atlas_aal.id, variables.atlas_aal.name, ba_tmp] = textread(parameters.atlas_aal_name, '%d %s %d');

vo = spm_vol(parameters.atlas_brodmann);
variables.atlas_ba_map = uint8(spm_read_vols(vo));

%% initialize some parameters
variables.neg_idx = [];
variables.pos_idx = [];

end