Hi, I convert a thickness file to ascii format and the file looks like
(vertex id, x, y, z, thickness) 000 -36.19709 -19.27060 66.10529 2.61633 001 -16.37191 -69.16402 59.60456 2.18290 002 -10.58665 -10.65137 50.24026 2.06269 003 -22.09574 43.85081 28.52032 1.54455 004 -57.95597 -0.45362 12.97830 2.41985 005 -48.57283 -51.02512 46.88505 2.08371
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? Thank you,
Ron
__________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com
Hi Ron,
convert one of the surfaces to ascii as well. It contains both a list of vertices (like the thickness) and a list of triangles, which has the nbhd info you need.
cheers, Bruce On Mon, 12 Jun 2006, R. Chen wrote:
Hi, I convert a thickness file to ascii format and the file looks like
(vertex id, x, y, z, thickness) 000 -36.19709 -19.27060 66.10529 2.61633 001 -16.37191 -69.16402 59.60456 2.18290 002 -10.58665 -10.65137 50.24026 2.06269 003 -22.09574 43.85081 28.52032 1.54455 004 -57.95597 -0.45362 12.97830 2.41985 005 -48.57283 -51.02512 46.88505 2.08371
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? Thank you,
Ron
Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com
From: "R. Chen" eemath2002@yahoo.com 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'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thanks Don,
one warning: Doug implemented Moo Chung's smoothing and we didn't find quite the same results he reported. There seemed to be some scale factor that was off (that is, the relationship between fwhm and iterations of diffusion smoothing wasn't quite right).
Thanks for posting this Don. Maybe now's also a good time to update people on the open source status. Basically we are working hard to remove code that we're not allowed to distribute from our source environment. This is a slow and painful process of identifying a replacement, writing a unit test for it, then trying to figure out if we can live with the 1e-5 difference in the results (and thanks to Dennis Jen for slogging through it). We're still hoping to post a fully downloadable and buildable open source version in the not-too-distant future. If there was enough interest I guess we could post a version that won't build, but at least will allow people to look at the source code themselves to answer questions.
cheers, Bruce
On Mon, 12 Jun 2006, Don Hagler wrote:
From: "R. Chen" eemath2002@yahoo.com 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
Yes, I did implement Moo's smoother. The problem was that it was impossible to comute the final fwhm based on the number of iterations, so it was not very useful.
doug
thanks Don,
one warning: Doug implemented Moo Chung's smoothing and we didn't find quite the same results he reported. There seemed to be some scale factor that was off (that is, the relationship between fwhm and iterations of diffusion smoothing wasn't quite right).
Thanks for posting this Don. Maybe now's also a good time to update people on the open source status. Basically we are working hard to remove code that we're not allowed to distribute from our source environment. This is a slow and painful process of identifying a replacement, writing a unit test for it, then trying to figure out if we can live with the 1e-5 difference in the results (and thanks to Dennis Jen for slogging through it). We're still hoping to post a fully downloadable and buildable open source version in the not-too-distant future. If there was enough interest I guess we could post a version that won't build, but at least will allow people to look at the source code themselves to answer questions.
cheers, Bruce
On Mon, 12 Jun 2006, Don Hagler wrote:
From: "R. Chen" eemath2002@yahoo.com 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
Freesurfer mailing list Freesurfer@nmr.mgh.harvard.edu https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
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
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
freesurfer@nmr.mgh.harvard.edu