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

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