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)