From: "R. Chen" <[EMAIL PROTECTED]>
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');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
