Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# vsoch/flattenSphere.m

Last active Dec 23, 2015
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
to join this conversation on GitHub. Already have an account? Sign in to comment