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…

flattenSphere.m

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

projectedSphere.m

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

unwrapSphere.m

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