{"id":5043,"date":"2014-01-03T09:00:02","date_gmt":"2014-01-03T14:00:02","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=5043"},"modified":"2014-01-03T15:28:48","modified_gmt":"2014-01-03T20:28:48","slug":"fit-a-sphere","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2014\/01\/03\/fit-a-sphere\/","title":{"rendered":"Fit a Sphere!"},"content":{"rendered":"
\r\n \r\n

Sean<\/a>'s pick this week is Sphere Fit<\/a> by Alan Jennings<\/a>.\r\n <\/p>\r\n <\/introduction>\r\n

Fitting a Sphere to Data<\/a><\/h3>\r\n

I recently had some data and wanted to fit a sphere through it so that I could find the radius of this sphere. As I started\r\n writing out an objective function for one of the Optimization Toolbox<\/a> optimizers (yes I was taking the way too big a hammer\r\n approach), a quick query on the File Exchange brought up sphereFit<\/tt>, a pleasant discovery.\r\n <\/p>\r\n

After downloading the files, I looked at Alan's published example<\/a>, it looked like it would work!\r\n <\/p>\r\n

Let's see it in action:<\/p>

% Load mri of a human head<\/span>\r\nS = load('mri'<\/span>);\r\nD = squeeze(S.D);\r\n\r\n% Generate data points on the edge of the head<\/span>\r\nfor<\/span> ii = size(D,3):-1:1;\r\n    % Perimeter of convex hull of each slice<\/span>\r\n    M(:,:,ii) = bwperim(bwconvhull(D(:,:,ii)>50));\r\nend<\/span><\/pre>
 <\/p>
% Find sub indices in our mask<\/span>\r\nidx = find(M); % find points in M<\/span>\r\n[rr,cc,pp] = ind2sub(size(M),idx); % sub indices<\/span>\r\npp = pp.*floor(size(D,1).\/size(D,3)); % rescale third dimension<\/span>\r\n\r\n% View data<\/span>\r\nfigure;\r\nscatter3(rr,cc,pp,10,pp);\r\ndaspect([1,1,1]);\r\nview(-121,36);\r\naxis tight<\/span>;<\/pre> 
% Fit the Sphere:<\/span>\r\n[cent,radius] = sphereFit([rr, cc, pp]);\r\nfprintf(1,'\\nRadius of sphere is %3.1f\\nIt is centered at [%3.1f %3.1f %3.1f]\\n'<\/span>,radius,cent);<\/pre>
\r\nRadius of sphere is 54.3\r\nIt is centered at [72.9 63.4 51.7]\r\n<\/pre>
Use Alan's example code to show the sphere through the points<\/p>
scatter3(cc,rr,pp,25,pp,'*'<\/span>); %points<\/span>\r\nhold on<\/span>; % hold<\/span>\r\ndaspect([1,1,1]); % equal axis so a sphere looks like a sphere<\/span>\r\n[Base_rr,Base_cc,Base_pp] = sphere(20);\r\nsurf(radius*Base_rr+cent(2),...<\/span>\r\n     radius*Base_cc+cent(1),...<\/span>\r\n     radius*Base_pp+cent(3),'faceAlpha'<\/span>,0.3,'Facecolor'<\/span>,'m'<\/span>)\r\naxis tight<\/span>;\r\nview(-121,36);<\/pre> 
Comments<\/a><\/h3>\r\n   
Give it a try and let us know what you think here<\/a> or leave a comment<\/a> for Alan.\r\n   <\/p>