mesh_center_mass - find center of mass for a group of vertices [center] = mesh_center_mass(xyz) xyz - Nx3 Cartesian coordinates center is the Cartesian coordinates for the center of mass.
0001 function center = mesh_center_mass(xyz) 0002 0003 % mesh_center_mass - find center of mass for a group of vertices 0004 % 0005 % [center] = mesh_center_mass(xyz) 0006 % 0007 % xyz - Nx3 Cartesian coordinates 0008 % 0009 % center is the Cartesian coordinates for 0010 % the center of mass. 0011 % 0012 0013 % $Revision: 1.1 $ $Date: 2004/11/12 01:30:25 $ 0014 0015 % Licence: GNU GPL, no implied or express warranties 0016 % History: 09/2004, Darren.Weber_at_radiology.ucsf.edu 0017 % Tania Boniske 0018 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0019 0020 version = '[$Revision: 1.1 $]'; 0021 fprintf('\nMESH_CENTER_MASS [v%s]\n',version(12:16)); tic; 0022 0023 % For 1 dimensional space, the center of gravity is: 0024 % 0025 % x_pos = ( sum(i = 1:n) ( mi * xi ) ) / sum(i = 1:n) mi; 0026 % y_pos = ( sum(i = 1:n) ( mi * yi ) ) / sum(i = 1:n) mi; 0027 % z_pos = ( sum(i = 1:n) ( mi * zi ) ) / sum(i = 1:n) mi; 0028 % 0029 % If mi = 1, then the above is simply the mean(xi) = 1/n * sum(xi). 0030 % 0031 % note how the sum of mi is in both the numberator/denominator, so the 0032 % units of the weights cancel out, leaving only the spatial units of xi. 0033 % 0034 % For 2 dimensions, the mi becomes 2D, so mi = sum(1:i,1:j) m(i,j); 0035 % For 3 dimensions, the mi becomes 3D, so mi = sum(1:i,1:j,1:k) m(i,j,k); 0036 % 0037 % mx(i) is intended to caluculate the scalar yz mass at each x(i). 0038 % xpos is then the division of the sum of the product of these scalar masses by the 0039 % position by the sum of the scalar masses. sum(mass*position)/sum(mass). 0040 % if the mass = 1 then the COG is the average position. 0041 0042 0043 printf('in development, not functioning correctly yet\n'); 0044 center = NaN; 0045 return 0046 0047 % This function is only useful if we have potential values at each vertex, 0048 % in which case we can return a center of "activity". 0049 0050 Nvertices = length(xyz); 0051 0052 mx = zeros(Nvertices,1); 0053 my = zeros(Nvertices,1); 0054 mz = zeros(Nvertices,1); 0055 0056 for i = 1:Nvertices, 0057 mx(i,1) = sum(sum(xyz(i,:,:))); 0058 my(i,1) = sum(sum(xyz(:,i,:))); 0059 mz(i,1) = sum(sum(xyz(:,:,i))); 0060 end 0061 xpos = (xi*mx)/sum(mx); 0062 ypos = (yj*my)/sum(my); 0063 zpos = (zk*mz)/sum(mz); 0064 0065 % Find center voxel of volume 0066 center = [xpos, ypos, zpos]; 0067 0068 t=toc; fprintf('...done (%6.2f sec)\n\n',t); 0069 0070 return