%SHOWMIP Show a SPM figure

% File name of the image to display
fname = 'D:\bizzell\SPM2test\Analysis\06143\snraV0001b.hdr';
% Voxels above this threshhold are considered "brain"
threshhold = 250;
% Percentage of brain voxels to plot on "glass brain"
P = 50;

vol = readmr(fname);

% Set up the transformation matrix and dimensions
M = zeros(4,4);
DIM = [];
M(4,4) = 1;
for i = 1:3
  DIM(i) = vol.info.dimensions(i).size;
  M(i,i) = vol.info.dimensions(i).spacing;
  if isnan(M(i,i)), error('Voxel dimension spacing undefined.'); end
  M(i,4) = -M(i,i) * vol.info.dimensions(i).origin;
  if isnan(M(i,4)) | M(i,4)==0
    M(i,4) = -M(i,i) * round(DIM(i)/2);
  end
end

% Calculate XYZ and Z matrices
idx = find(vol.data > threshhold);
idx = idx(1:round(100/P):length(idx));

XYZ = zeros(3,length(idx));
for i=1:length(idx)
  XYZ(3,i) = floor(idx(i)/(DIM(1)*DIM(2)));
  XYZ(2,i) = floor((idx(i) - (XYZ(3,i))*(DIM(1)*DIM(2)))/DIM(1));
  XYZ(1,i) = (idx(i) - XYZ(3,i)*(DIM(1)*DIM(2))) - XYZ(2,i)*DIM(1);
end
XYZ(3,:) = XYZ(3,:) + 1;
XYZ(2,:) = XYZ(2,:) + 1;

XYZ = M(1:3,:)*[XYZ; ones(1,size(XYZ,2))];
%Z = ones(1,size(XYZ,2));
Z = rand(1,size(XYZ,2))*50;

F = figure;
colormap(gray);
spm_mip(Z,XYZ,M,DIM)