RE: [Freesurfer] neighborhood of a vertex

Hi Ron,

"fs_fread3" is actually just what I renamed the fread3.m found in $FREESURFER_HOME/matlab. I actually made no changes; just wanted my own self-contained set of matlab scripts. Anyway, here it is:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [retval] = fs_fread3(fid) % [retval] = fd3(fid) % read a 3 byte integer out of a file b1 = fread(fid, 1, 'uchar') ; b2 = fread(fid, 1, 'uchar') ; b3 = fread(fid, 1, 'uchar') ; retval = bitshift(b1, 16) + bitshift(b2,8) + b3 ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

From: "R. Chen" eemath2002@yahoo.com To: Don Hagler dhaglerjr@hotmail.com CC: freesurfer@nmr.mgh.harvard.edu Subject: RE: [Freesurfer] neighborhood of a vertex Date: Tue, 13 Jun 2006 07:50:04 -0700 (PDT)

Thank you, Don. However, I try your code and matlab tells me 'Undefined function or variable 'fs_fread3''. Could you post other matlab function which are used by your 'fs_read_surf.m' and ' fs_find_neighbors.m'?

Ron

Don Hagler dhaglerjr@hotmail.com wrote: >From: "R. Chen"

...

To: freesurfer@nmr.mgh.harvard.edu Subject: [Freesurfer] neighborhood of a vertex Date: Mon, 12 Jun 2006 11:02:59 -0700 (PDT)

I want to do some customrized smoothing and need to know the neighborhood of a vertex. I think FreeSurfer also uses this neighborhood information in random field modelling. Could you provide some information about how to get the neighborhood?

Take a look at $FREESURFER_HOME/matlab/read_surf.m. Darren Webber at UCSF has something similar.

Moo Chung at UWisc has a function called mni_getmesh.m that find the neighbors to enable smoothing and he has a function for heat kernel smoothing.

I've used those sources and made my own modifications that some might find useful.

Here is my modification of read_surf: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [surf] = fs_read_surf(fname) % fs_read_surf - read a freesurfer surface file % % [surf] = fs_read_surf(fname) % % Reads the vertex coordinates (mm) and face lists from a surface file % Then finds neighbors for each vertex % % surf is a structure containg: % nverts: number of vertices % nfaces: number of faces (triangles) % faces: vertex numbers for each face (3 corners) % coords: x,y,z coords for each vertex % nbrs: vertex numbers of neighbors for each vertex % % created: 03/02/06 Don Hagler % last modified: 05/09/06 Don Hagler % % code for reading surfaces taken from Darren Weber's freesurfer_read_surf % % see also: fs_read_trisurf, fs_find_neighbors, fs_calc_triarea %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

funcname = 'fs_read_surf';

%QUAD_FILE_MAGIC_NUMBER = (-1 & 0x00ffffff) ; %NEW_QUAD_FILE_MAGIC_NUMBER = (-3 & 0x00ffffff) ;

TRIANGLE_FILE_MAGIC_NUMBER = 16777214 ; QUAD_FILE_MAGIC_NUMBER = 16777215 ;

% open it as a big-endian file fid = fopen(fname, 'rb', 'b') ; if (fid < 0), str = sprintf('%s: could not open surface file %s.',funcname,fname) ; error(str) ; end

magic = fs_fread3(fid) ;

if (magic == QUAD_FILE_MAGIC_NUMBER), surf.nverts = fs_fread3(fid) ; surf.nfaces = fs_fread3(fid) ; fprintf('%s: reading %d quad file vertices...',funcname,surf.nverts); tic; surf.coords = fread(fid, surf.nverts*3, 'int16') ./ 100 ; t=toc; fprintf('done (%0.2f sec)\n',t); fprintf('%s: reading %d quad file faces (please wait)...\n',... funcname,surf.nfaces); tic; surf.faces = zeros(surf.nfaces,4); for iface = 1:surf.nfaces, for n=1:4, surf.faces(iface,n) = fs_fread3(fid) ; end if(~rem(iface, 10000)), fprintf(' %7.0f',iface); end if(~rem(iface,100000)), fprintf('\n'); end end t=toc; fprintf('\ndone (%0.2f sec)\n',t); elseif (magic == TRIANGLE_FILE_MAGIC_NUMBER), fprintf('%s: reading triangle file...',funcname); tic; tline = fgets(fid); % read creation date text line tline = fgets(fid); % read info text line

surf.nverts = fread(fid, 1, 'int32') ; % number of vertices
surf.nfaces = fread(fid, 1, 'int32') ; % number of faces

% vertices are read in column format and reshaped below
surf.coords = fread(fid, surf.nverts*3, 'float32');

% faces are read in column format and reshaped
surf.faces = fread(fid, surf.nfaces*3, 'int32') ;
surf.faces = reshape(surf.faces, 3, surf.nfaces)' ;
t=toc; fprintf('done (%0.2f sec)\n',t);

else str = sprintf('%s: unknown magic number in surface file %s.',funcname,fname) ; error(str) ; end

surf.coords = reshape(surf.coords, 3, surf.nverts)' ; fclose(fid) ;

%fprintf('...adding 1 to face indices for matlab compatibility.\n\n'); surf.faces = surf.faces + 1;

return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is my function for getting neighbors: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [surf] = fs_find_neighbors(surf) % fs_read_surf - find neighboring relations between vertices in a surface % % [surf] = fs_find_neighbors(surf) % % Input: % surf is a structure containg: % nverts: number of vertices % nfaces: number of faces (triangles) % faces: vertex numbers for each face (3 corners) % coords: x,y,z coords for each vertex % % Output: % surf is a structure containg: % nverts: number of vertices % nfaces: number of faces (triangles) % faces: vertex numbers for each face (3 corners) % coords: x,y,z coords for each vertex % nbrs: vertex numbers of neighbors for each vertex % % created: 05/09/06 Don Hagler % last modified: 05/09/06 Don Hagler % % code for finding neighbors taken from Moo Chung's mni_getmesh % % see also: fs_read_surf, fs_read_trisurf, fs_calc_triarea %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

funcname = 'fs_find_neighbors';

if ~isfield(surf, 'faces') fprintf('%s: error: input surf must contain faces\n',funcname); return; end;

% compute the maximum degree of node -- number of edges = number of neighbors fprintf('%s: finding number of nearest neighbors...',funcname); tic; num_nbrs=zeros(surf.nverts,1); for i=1:surf.nfaces num_nbrs(surf.faces(i,:))=num_nbrs(surf.faces(i,:))+1; end max_num_nbrs=max(num_nbrs); t=toc; fprintf('done (%0.2f sec)\n',t);

% find nearest neighbors fprintf('%s: finding nearest neighbors...',funcname); tic; surf.nbrs=zeros(surf.nverts,max_num_nbrs); for i=1:surf.nfaces for j=1:3 vcur = surf.faces(i,j); for k=1:3 if (j ~= k) vnbr = surf.faces(i,k); if find(surf.nbrs(vcur,:)==vnbr) ; else n_nbr = min(find(surf.nbrs(vcur,:) == 0)); surf.nbrs(vcur,n_nbr) = vnbr; end; end; end; end; end; t=toc; fprintf('done (%0.2f sec)\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is my function for smoothing: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [smoothedvals] = fs_smooth(surf,vals,niter) %function [smoothedvals] = fs_smooth(surf,vals,niter) % % surf is a structure containing: % coords: x,y,z coords for each surface vertex % nbrs: vertex numbers of neighbors for each vertex % nverts: number of vertices % % vals is vector with nverts members % % niter is number of iterations (smoothing steps) %

funcname = 'fs_smooth';

smoothedvals = [];

if(nargin ~= 3) fprintf('USAGE: [smoothedvals] = %s(surf,vals,niter) \n',funcname); return; end

if size(vals,2)~=1 fprintf('%s: vals must have a single column (it has %d)\n',... funcname,size(vals,2)); end;

if size(vals,1)~=surf.nverts fprintf('%s: number of vals (%d) does not match number of verts (%d)\n',... funcname,size(vals,1),surf.nverts); return; end;

fprintf('%s(%d): smoothing',funcname,niter); vals = [0;vals]; % create dummy vertex with index 1, value 0 maxnbrs = size(surf.nbrs,2); surf.nbrs = [ones(maxnbrs,1)';surf.nbrs + 1]; % nbrs contains zeros, change to 1 num_nbrs = sum(surf.nbrs>1,2)+1; % including vertex and its neighbors for iter=1:niter fprintf('.'); Y=[vals, vals(surf.nbrs(:,:))]; % sum values from vertex and its neighbors vals=sum(Y,2)./num_nbrs; end; smoothedvals = vals(2:end);

fprintf('\n'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

---

Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com