avw_center_mass - find center of mass for an image volume [center] = avw_center_mass(avw) avw - an Analyze 7.5 data struct, see avw_read center is the Cartesian coordinates for the volume center of mass. It is a struct with coordinates in both voxels and mm (1x3). see also avw_center
0001 function center = avw_center_mass(avw) 0002 0003 % avw_center_mass - find center of mass for an image volume 0004 % 0005 % [center] = avw_center_mass(avw) 0006 % 0007 % avw - an Analyze 7.5 data struct, see avw_read 0008 % 0009 % center is the Cartesian coordinates for 0010 % the volume center of mass. It is a struct with 0011 % coordinates in both voxels and mm (1x3). 0012 % 0013 % see also avw_center 0014 % 0015 0016 % $Revision: 1.1 $ $Date: 2004/11/12 01:30:25 $ 0017 0018 % Licence: GNU GPL, no implied or express warranties 0019 % History: 09/2004, Darren.Weber_at_radiology.ucsf.edu 0020 % Tania Boniske 0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0022 0023 version = '[$Revision: 1.1 $]'; 0024 fprintf('\nAVW_CENTER_MASS [v%s]\n',version(12:16)); tic; 0025 0026 % Extract info from avw.hdr 0027 xdim = double(avw.hdr.dime.dim(2)); 0028 ydim = double(avw.hdr.dime.dim(3)); 0029 zdim = double(avw.hdr.dime.dim(4)); 0030 0031 xpix = double(avw.hdr.dime.pixdim(2)); 0032 ypix = double(avw.hdr.dime.pixdim(3)); 0033 zpix = double(avw.hdr.dime.pixdim(4)); 0034 0035 % For 1 dimensional space, the center of gravity is: 0036 % 0037 % x_pos = ( sum(i = 1:n) ( mi * xi ) ) / sum(i = 1:n) mi; 0038 % y_pos = ( sum(i = 1:n) ( mi * yi ) ) / sum(i = 1:n) mi; 0039 % z_pos = ( sum(i = 1:n) ( mi * zi ) ) / sum(i = 1:n) mi; 0040 % 0041 % If mi = 1, then the above is simply the mean(xi) = 1/n * sum(xi). 0042 % 0043 % note how the sum of mi is in both the numberator/denominator, so the 0044 % units of the weights cancel out, leaving only the spatial units of xi. 0045 % 0046 % For 2 dimensions, the mi becomes 2D, so mi = sum(1:i,1:j) m(i,j); 0047 % For 3 dimensions, the mi becomes 3D, so mi = sum(1:i,1:j,1:k) m(i,j,k); 0048 % 0049 % mx(i) is intended to caluculate the scalar yz mass at each x(i). 0050 % xpos is then the division of the sum of the product of these scalar masses by the 0051 % position by the sum of the scalar masses. sum(mass*position)/sum(mass). 0052 % if the mass = 1 then the COG is the average position. 0053 0054 xi = 1:xdim; 0055 yj = 1:ydim; 0056 zk = 1:zdim; 0057 0058 mx = zeros(xdim,1); 0059 my = zeros(ydim,1); 0060 mz = zeros(zdim,1); 0061 0062 for i = xi, 0063 mx(i,1) = sum(sum(avw.img(i,:,:))); 0064 end 0065 xpos = (xi*mx)/sum(mx); 0066 0067 for j = yj, 0068 my(j,1) = sum(sum(avw.img(:,j,:))); 0069 end 0070 ypos = (yj*my)/sum(my); 0071 0072 for k = zk, 0073 mz(k,1) = sum(sum(avw.img(:,:,k))); 0074 end 0075 zpos = (zk*mz)/sum(mz); 0076 0077 % Find center voxel of volume 0078 center.voxels = [xpos, ypos, zpos]; 0079 0080 % Find mm coordinates of that voxel 0081 center.mm = center.voxels .* [xpix, ypix, zpix]; 0082 0083 t =toc; fprintf('...done (%6.2f sec)\n\n',t); 0084 0085 return