2009 Grace Tang & Kelly Hennigan Appendix
Main Code
%%Source: http://sites.google.com/site/mvlombardo/matlab-tutorials
%%comments that don't start with 'Grace:' were written by website author
SubjArray={'F-1', 'F-2', 'F-4', 'I-6', 'H-5', 'C-1', 'C-4', 'C-5', 'C-7', 'J-5'}; %Grace: change subject names here! SubjArray2={'F1', 'F2', 'F4', 'I6', 'H5', 'C1', 'C4', 'C5', 'C7', 'J5'}; %Grace: ...and here, if your segmentation file's subject names were different Datamat=zeros(157, 189, 136, length(SubjArray)); % Grace: creates matrix where last number is subject index number
for i=1:length(SubjArray)
datadir = ['/Users/grace/Desktop/PSYCH204A/' SubjArray{i} '/']; % datadir = ['/Users/grace/Desktop/PSYCH204A/']; fname = ['w' SubjArray2{i} '_NAcc_seg_clean.nii']; Data2Read = fullfile(datadir,fname);
% What I've done in Step 1 is to construct a filename such as /data/fMRI/12001.nii from its filepath (the variable called datadir) and filename (the variable called fname). Both these variables are string variables and I've used the function fullfile to essentially concatenate these two strings together. The end result is a string variable called Data2Read, which is essential a string that looks like this: /data/fMRI/12001.nii
% Step 2. Now pass Data2Read into the function spm_vol in order to read in the header information for the file. Note that the spm_vol function comes with the SPM toolbox in Matlab that is a standard software package for neuroimagers. Having SPM is essential for this part.
HeaderInfo = spm_vol(Data2Read);
% In Step 2, the header information for the data you want to read in is now in your workspace. This is required for Step 3, where you will use the header information to actually read in the entire data.
% Step 3. Now use spm_read_vols to read in the data Data = spm_read_vols(HeaderInfo);
%% Grace: converting non zero matrix values to 1
indx=find(Data); Data(indx)=1; Datamat(:,:,:,i)=Data; %Grace: Datamat now contains all your matrix data
% %% saving all values in file to '1'
% HeaderInfo.fname = ['w' SubjArray2{i} '_NAcc_seg_allten.nii']; % This is where you fill in the new filename
% HeaderInfo.private.dat.fname = HeaderInfo.fname; % This just replaces the old filename in another location within the header.
% %
% % Step 2. Now use spm_write_vol to write out the new data. You need to give spm_write_vol the new header information and corresponding data matrix
% spm_write_vol(HeaderInfo,Data);
end
combinemat_g=sum(Datamat,4); %Grace: sum across subjects... now values indicate overlap across your subjects
%%%
%% histogram
index2=find(combinemat_g); %G: find non zero elements, so histogram doesn't have massive bar at 0; toplot=combinemat_g(index2); hist(toplot,9); % G: hardcoded to 9 bins, because maximum overlap was 9 ylabel('# voxels with x subjects overlapping') xlabel('x - number of subjects per voxel')
figure(2)
rhist(toplot,9);
ylabel('relative frequency of voxels with x subjects overlapping') xlabel('x - number of subjects per voxel')
%% figure 3: cumulative hist
figure(3) n_elements=histc(toplot, 1:9); c_elements=cumsum(n_elements); bar(1:9, c_elements); ylabel('cumulative number of voxels with overlap above x number of subjects') xlabel('x number of subjects')
%%% % % %% writing data out to nifti % % %%%% % % Step 1. Take the header information from a previous file with similar dimensions and voxel sizes and change the filename in the header. % HeaderInfo.fname = 'Combinemat_N10_allten_smooth.nii'; % This is where you fill in the new filename % HeaderInfo.private.dat.fname = HeaderInfo.fname; % This just replaces the old filename in another location within the header. % % % Step 2. Now use spm_write_vol to write out the new data. You need to give spm_write_vol the new header information and corresponding data matrix % spm_write_vol(HeaderInfo,combinemat_g); % where HeaderInfo is your header information for the new file, and Data is the image matrix corresponding to the image you'll be writing out.
More Code
%% Analysis
% rgb color code for labels in ITK-Snap (match bar colors with labels)
% 0 0 0 0 0 0 0 "Clear Label"
% 1 85 36 255 1 1 1 "Label 1"
% 2 7 81 255 1 1 1 "Label 2"
% 3 35 252 255 1 1 1 "Label 3"
% 4 0 255 145 1 1 1 "Label 4"
% 5 87 255 21 1 1 1 "Label 5"
% 6 224 255 10 1 1 1 "Label 6"
% 7 231 178 0 1 1 1 "Label 7"
% 8 228 96 0 1 1 1 "Label 8"
% 9 255 2 0 1 1 1 "Label 9"
% 10 252 255 249 1 1 1 "Label 10"
%
% Specify data sets with regular voxel size and double voxel size
dat = reg_vox_10; % 1x1x1 voxel size, matrix for each subject dat_sum = sum(reg_vox_10,4); % 1x1x1 voxel size, summed matrix across subjects
dat2 = MNIeach_ten; % 2x2x2 voxel size, matrix for each subject dat2_sum = sum(MNIeach_ten,4); % 2x2x2 voxel size, summed matrix across subjects dat_rep = MNI_NAcc_rep; % SPM's representative subject NAc mask
%% 1x1x1 voxels
% Size - finds the # of voxels in each subjects mask then calculates the average and the s.e.
for k = 1:size(dat,4); % 2x2x2
indx = find(dat(:,:,:,k)); subj_vox_count(k) = length(indx);
end
mean_nvox = mean(subj_vox_count); se_nvox = std(subj_vox_count)/(size(dat,4)^.5);
%% PLOTS
forhist = [];
for i=1:9
indx = find(dat_sum==i); forhist(i) = length(indx);
end forhist = forhist';
% ITK Snap Label colors mycolors = [255 2 0; 228 96 0; 231 178 0; 224 255 10; 87 255 21; 0 255 145; 35 252 255; 7 81 255; 85 36 255]./255;
% Histogram figure(1), h = bar(forhist) colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('number of NAc voxels') title('Subject overlap in NAc voxels (1 x 1 x 1 mm)') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(forhist,1); for j = 1:length(forhist)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
clear h ch fvd fvcd zs izs j
% Cumulative Histogram c_elements = cumsum(forhist); figure(2), h = bar(c_elements) colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('number of NAc voxels') title('Cumulative Histogram of subject overlap in NAc voxels (1 x 1 x 1 mm)') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(c_elements,1); for j = 1:length(c_elements)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
clear h ch fvd fvcd zs izs j
% proportional histogram (rhist) rhist_bars = forhist./sum(forhist); figure(3), h = bar(rhist_bars); colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('proportion of total NAc voxels') title('Relative histogram of subject overlap in NAc voxels (1 x 1 x 1 mm)') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(rhist_bars,1); for j = 1:length(rhist_bars)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
clear h ch fvd fvcd zs izs j c_elements forhist i indx indx2 k temp
%% 2x2x2 voxels
% Size - finds the # of voxels in each subjects mask then calculates the average and the s.e.
for k = 1:size(dat2,4); % 2x2x2
indx = find(dat2(:,:,:,k)); subj_vox_count2(k) = length(indx);
end
mean_nvox2 = mean(subj_vox_count2); se_nvox2 = std(subj_vox_count2)/(size(dat2,4)^.5);
%% PLOTS
forhist = [];
for i=1:9
indx = find(dat_sum==i); forhist(i) = length(indx);
end forhist = forhist';
% ITK Snap Label colors mycolors = [255 2 0; 228 96 0; 231 178 0; 224 255 10; 87 255 21; 0 255 145; 35 252 255; 7 81 255; 85 36 255]./255;
% Histogram figure(4), h = bar(forhist) colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('number of NAc voxels') title('Subject overlap in NAc voxels (2 x 2 x 2 mm)') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(forhist,1); for j = 1:length(forhist)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
clear h ch fvd fvcd zs izs j
% Cumulative Histogram c_elements = cumsum(forhist); figure(5), h = bar(c_elements) colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('number of NAc voxels') title('Cumulative Histogram of subject overlap in NAc voxels (2 x 2 x 2 mm)') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(c_elements,1); for j = 1:length(c_elements)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
clear h ch fvd fvcd zs izs j
% proportional histogram (rhist) rhist_bars2 = forhist./sum(forhist); figure(6), h = bar(rhist_bars2); colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('proportion of total NAc voxels') title('Relative histogram of subject overlap in NAc voxels (2 x 2 x 2 mm)') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(rhist_bars2,1); for j = 1:length(rhist_bars2)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
%% Compare to single subject rep NAc
in_indx = find(dat_rep); % # of NAc voxels in Colin mask in_rep_mask = dat_sum(in_indx); out_indx = find(dat_rep==0); out_temp = dat_sum(out_indx); out_indx2 = find(out_temp); % # of NAc voxels out of Colin mask out_rep_mask = out_temp(out_indx2);
mean_in_mask = mean(in_rep_mask); % mean subject overlap/NAc voxel within Colin mask mean_out_mask = mean(out_rep_mask); % mean subject overlap/NAc voxel outside of mask
clear h ch fvd fvcd zs izs j c_elements forhist i indx indx2 k temp
% plot relative hist for reg and double size voxels next to each other
% figure(7) % ax(1) = subplot(1,2,1); % rgb = imread('ngc6543a.jpg'); % image(rgb); title('RGB image') % ax(2) = subplot(122); % im = mean(rgb,3); % image(im); title('Intensity Heat Map') % colormap(hot(256)) % linkaxes(ax,'xy') % axis(ax,'image')
% proportional histogram (rhist)
figure(7) subplot(121), h = bar(rhist_bars); colormap(mycolors); xlabel('number of subjects overlap/ NAc voxel') ylabel('proportion of total NAc voxels') title('1 x 1 x 1 mm voxels ') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(rhist_bars,1); for j = 1:length(rhist_bars)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd) subplot(122),h = bar(rhist_bars2); colormap(mycolors); title('2 x 2 x 2 mm voxels') ch = get(h,'Children'); %getting a handle on the kids fvd = get(ch,'Faces'); fvcd = get(ch,'FaceVertexCData'); [zs, izs] = sortrows(rhist_bars2,1); for j = 1:length(rhist_bars2)
row = izs(j); fvcd(fvd(row,:)) = j;
end set(ch,'FaceVertexCData',fvcd)
Smoothing Code
%%Source: %%http://sites.google.com/site/mvlombardo/matlab-tutorials %%comments that don't start with 'Grace:' were written by website author
SubjArray={'F-1', 'F-2', 'F-4', 'I-6', 'H-5', 'C-1', 'C-4', 'C-5', 'C-7', 'J-5'}; %Grace: change subject names here! SubjArray2={'F1', 'F2', 'F4', 'I6', 'H5', 'C1', 'C4', 'C5', 'C7', 'J5'}; %Grace: ...and here, if your segmentation file's subject names were different Datamat=zeros(157, 189, 136, length(SubjArray)); % Grace: creates matrix where last number is subject index number
smoothmat_=zeros(157,189,136,length(SubjArray));
for i=1:length(SubjArray)
%datadir = ['/Users/grace/Desktop/PSYCH204A/' SubjArray{i} '/']; datadir = '/Users/grace/Desktop/PSYCH204A/'; fname = ['w' SubjArray2{i} '_NAcc_seg_allone.nii']; Data2Read = fullfile(datadir,fname);
% What I've done in Step 1 is to construct a filename such as /data/fMRI/12001.nii from its filepath (the variable called datadir) and filename (the variable called fname). Both these variables are string variables and I've used the function fullfile to essentially concatenate these two strings together. The end result is a string variable called Data2Read, which is essential a string that looks like this: /data/fMRI/12001.nii
% Step 2. Now pass Data2Read into the function spm_vol in order to read in the header information for the file. Note that the spm_vol function comes with the SPM toolbox in Matlab that is a standard software package for neuroimagers. Having SPM is essential for this part.
HeaderInfo = spm_vol(Data2Read);
% In Step 2, the header information for the data you want to read in is now in your workspace. This is required for Step 3, where you will use the header information to actually read in the entire data.
% Step 3. Now use spm_read_vols to read in the data Data = spm_read_vols(HeaderInfo);
%% Grace: converting non zero matrix values to 1
indx=find(Data); Data(indx)=1; Datamat(:,:,:,i)=Data; %Grace: Datamat now contains all your matrix data
%% smoothing with smooth3
%smoothmat(:,:,:,i)=smooth3(Datamat(:,:,:,i), 'gaussian'); %smoothmat2(:,:,:,i)=smooth3(Datamat(:,:,:,i), 'gaussian', [1 1 1]); % smoothing kernel of 1, since voxel size is 1mm cube smoothmat3(:,:,:,i)=smooth3(Datamat(:,:,:,i), 'gaussian', [5 5 5]); %
for x=1:157
for y=1:189
for z=1:136
if smoothmat3(x,y,z,i)>=0.3 % vary threshold here
smoothmat_(x,y,z,i)=1;
else
smoothmat_(x,y,z,i)=0;
end
end
end
end
%%
%
% % %% saving all values in file to '1'
% % HeaderInfo.fname = ['w' SubjArray2{i} '_NAcc_seg_allten.nii']; % This is where you fill in the new filename
% % HeaderInfo.private.dat.fname = HeaderInfo.fname; % This just replaces the old filename in another location within the header.
% % %
% % % Step 2. Now use spm_write_vol to write out the new data. You need to give spm_write_vol the new header information and corresponding data matrix
% % spm_write_vol(HeaderInfo,Data);
%
end
combine_smooth_rethresh2=sum(smoothmat_,4); % % combinemat_g=sum(Datamat,4); %Grace: sum across subjects... now values indicate overlap across your subjects
%%
%%%
%% histogram
index2=find(combine_smooth_rethresh2); %G: find non zero elements, so histogram doesn't have massive bar at 0; toplot=combine_smooth_rethresh2(index2); hist(toplot,9); % G: hardcoded to 9 bins, because maximum overlap was 9 ylabel('# voxels with x subjects overlapping') xlabel('x - number of subjects per voxel')
figure(2)
rhist(toplot,9);
ylabel('relative frequency of voxels with x subjects overlapping') xlabel('x - number of subjects per voxel')
% %%%
% %
% % %% writing data out to nifti
% %
% % %%%%
% % Step 1. Take the header information from a previous file with similar dimensions and voxel sizes and change the filename in the header.
% HeaderInfo.fname = 'Combinesmooth_rethresh111_03.nii'; % This is where you fill in the new filename
% HeaderInfo.private.dat.fname = HeaderInfo.fname; % This just replaces the old filename in another location within the header.
%
% % Step 2. Now use spm_write_vol to write out the new data. You need to give spm_write_vol the new header information and corresponding data matrix
% spm_write_vol(HeaderInfo,combine_smooth_rethresh2); % where HeaderInfo is your header information for the new file, and Data is the image matrix corresponding to the image you'll be writing out.