Skip to content

Instantly share code, notes, and snippets.

@vsoch

vsoch/flattenSphere.m

Last active Dec 23, 2015
Embed
What would you like to do?
flattenSphere.m takes a 3D image with x y z coordinates defined on the circumference of a circle, and flattens it to a 2D representation (with some obvious distortion). unwrapSphere.m estimates a radius and centroid with least squares, and then writes data to a 3D matrix defined by azimuth, elevation, and radius. Since this sphere is hollow, the…
function flat = flattenSphere(mr)
% First read in sphere files and vox mapping
%vox = spm_read_vols(spm_vol('NDARAK333GZA.lh.vtxvol.nii'));
mr = spm_read_vols(spm_vol(mr));
% First create a random maximum length - we will crop at the end
maxlength = floor(pi*size(mr,1));
% Now go through vox image and save vector of x,y,z coordinates
flat = zeros(size(mr,1),maxlength);
% Create 2D representations of sphere
for i=1:size(mr,3) % All 3 dims should be the same
%map = squeeze(vox(:,:,i));
slice = squeeze(mr(:,:,i));
%slice = slice .* (map~=-1);
% CUT slice image in slices pieces, flatten (sum) up pieces, then concatenate!
piece1 = sum(imrotate(slice(1:128,1:128),-40)); piece1 = piece1(find(piece1,1,'first'):find(piece1,1,'last'));
piece2 = sum(imrotate(slice(129:end,1:128),40)); piece2 = piece2(find(piece2,1,'first'):find(piece2,1,'last'));
piece3 = sum(imrotate(slice(129:end,129:end),-40)); piece3 = piece3(find(piece3,1,'first'):find(piece3,1,'last'));
piece4 = sum(imrotate(slice(1:128,129:end),40)); piece4 = piece4(find(piece4,1,'first'):find(piece4,1,'last'));
vector = [ piece1 piece2 piece3 piece4 ];
diffsize = maxlength - length(vector);
vector = [ zeros(1,ceil(diffsize/2)) vector zeros(1,ceil(diffsize/2))] ;
if ~isempty(vector)
vector = vector(1:maxlength);
end
flat(i,1:length(vector)) = vector;
end
% Crop the image on top and bottom
top = find(sum(flat,2),1,'first');
bot = find(sum(flat,2),1,'last');
up = find(sum(flat,1),1,'first');
dwn = find(sum(flat,1),1,'last');
flat = flat(top:bot,up:dwn);
imagesc(flat);
title('Breaking circumference into pieces, flattening, and concatenating');
end
function [flatlh,flatrh] = projectedSphere(mr)
% This function creates a 2D projection of a sphere, however there are
% overlapping voxels!
% First read in sphere files and vox mapping
%vox = spm_read_vols(spm_vol('NDARAK333GZA.lh.vtxvol.nii'));
mr = spm_read_vols(spm_vol(mr));
% split in half down the middle
mrlh = mr(:,:,1:size(mr,3)/2);
mrrh = mr(:,:,size(mr,3)/2:end);
flatlh = flatten2D(mrlh);
flatrh = flatten2D(mrrh);
function Z = flatten2D(data)
x = [];
y = [];
z = [];
for i=1:size(data,3)
[xtmp ytmp] = find(data(:,:,i) ~= 0);
if size(xtmp,1)~=0
x= [x; xtmp];
y= [y; ytmp];
z= [z; repmat(i,size(xtmp,1),1)];
end
end
R = [x y z]; clear x y z;
% Get min and max of dimensions
xmin = min(R(:,1)); xmax = max(R(:,1));
ymin = min(R(:,2)); ymax = max(R(:,2));
Z = zeros(xmax,ymax);
repeatcount = 0;
for i=1:size(R,1)
if (Z(R(i,1),R(i,2)) ~= 0)
repeatcount=repeatcount+1;
end
Z(R(i,1),R(i,2)) = data(R(i,1),R(i,2),R(i,3));
end
% Find left edge of circle and crop
top = find(sum(Z)==0); top = top(end);
left = find(sum(Z,2)==0); left = left(end);
Z = Z(top:end,left:end);
figure; imagesc(Z);
end
end
function unwrappedSphereFlat = unwrapSphere(mr)
% Unwrap a sphere, and express in terms of azimuth (0 to 360 degrees),
% elevation (-90, 90), and radius. Since matlab indexing cannot be
% negative, we shift the elevation (y) 90 in the positive direction
% REQUIRES: http://www.mathworks.com/matlabcentral/fileexchange/34129-sphere-fit-least-squared
% You may need to tweak the strategy for assigning x,y,z coordinates (eg, round vs ceil vs floor)
% in the case that a coordinate comes out == 0 or greater than the size of the image
% First read in sphere files and vox mapping
%vox = spm_read_vols(spm_vol('NDARAK333GZA.lh.vtxvol.nii'));
mr = spm_read_vols(spm_vol(mr));
%% 1) CALCULATING CENTROID AND RADIUS
% Extract all xyz coordinates
x = [];
y = [];
z = [];
for i=1:size(mr)
[xtmp ytmp] = find(mr(:,:,i) ~= 0);
if size(xtmp,1)~=0
x= [x; xtmp];
y= [y; ytmp];
z= [z; repmat(i,size(xtmp,1),1)];
end
end
R = [x y z]; clear x y z
[center,radius] = sphereFit(R);
%% 2) UNWRAP SPHERE INTO 3D MATRIX
azimuth = 1:360;
elevation = -90:90;
unwrappedSphere = zeros(360,180,ceil(radius));
for a=azimuth
for e=elevation
for r=1:radius
% Calculate x, y, and z coordinates
% This would be for data centered at the origin, so
% we need to add the displacement to index the mr data
x1 = round(r .* cos(e) .* cos(a) + center(1));
y1 = round(r .* cos(e) .* sin(a) + center(2));
z1 = round(r .* sin(e) + center(3));
% This would be to save the coordinates themselves
%x = [ x; x1]; y = [ y; y1]; z = [ z; z1];
unwrappedSphere(a,e+91,r) = mr(x1,y1,z1);
end
end
end
%% 3) VISUALIZE RESULT
% Let's count how many values we have from 90 to radius
imagesc(squeeze(sum(unwrappedSphere(:,:,90:end)~=0,3)))
title('Number of overlapping voxels from 90 to radius')
colorbar
% Try squashing over the radius values (flattening the image)
unwrappedSphereFlat = zeros(size(unwrappedSphere,1),size(unwrappedSphere,2));
for i=1:size(unwrappedSphereFlat,1)
for j=1:size(unwrappedSphereFlat,2)
unwrappedSphereFlat(i,j) = sum(unwrappedSphere(i,j,:));
end
end
imagesc(squeeze(sum(unwrappedSphere(:,:,90:end)~=0,3)))
colorbar; title('Image flattened over radius dimension')
xlabel('elevation + 90'); ylabel('Azimuth')
% The matrix is empty up through radius == ~90, so we should look at
% between 90 and radius to get a sense of our values
figure(1);
for i=90:96
subplot(2,4,i-89)
imagesc(unwrappedSphere(:,:,i));
title([ 'Radius = ' num2str(i) ] )
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.