=================[Matlab Script: StartUp]======================

% freesurfer start

 

fshome = getenv('FREESURFER_HOME');

fsmatlab = sprintf('%s/matlab',fshome);

if (exist(fsmatlab) == 7)

    path(path,fsmatlab);

end

clear fshome fsmatlab;

=================[Matlab Script: make_outer_surface]================

function make_outer_surface (filled_volume, se_diameter, output_surface)

 

% Original Author: Marie Schaer

% Date: 2007/11/14

% 

% This function takes as an input the binary volume resulting of the

% filling of the any surface (usually the pial one) using mris_fill, and

% will close the sulci using morphological operation, with a sphere as the

% structural element.

% 

% Parameters: 

% se_diameter is the diameter of the sphere (in mm), use 15mm by default 

% to close the sulci.  

% 

% Utilities to write the surface to the freesurfer format are modified from

% "freesurfer_write_surf" (from Darren Weber's bioelectromagnetism toolbox), 

% according to a suggestion by Don Hagler in FreeSurfer's mailing list on 

% August 3, 2007.

%

% Example: make_outer_surface('lh.pial.mgz',15,'lh.outer-pial')

 

 

    fprintf('reading filled volume...\n');

    vol=MRIread(filled_volume);

    volume=vol.vol;

    volume(volume==1)=255;

    fprintf('closing volume...\n');

    

% first apply a very soft gaussian filter, with sigma = 1mm, in order to

% facilitate the closing

Gaussian = fspecial('gaussian',[2 2],1);

image_f=zeros(256,256,256);

for slice=1:256

    temp = double(reshape(volume(:,:,slice),256,256));

    image_f(:,:,slice) = conv2(temp,Gaussian,'same');

end

image2=zeros(size(image_f));

image2(image_f<=25)=0;

image2(image_f>25)=255;

       

    se=strel('ball',se_diameter,se_diameter);

    BW2=imclose(image2,se);

    thresh = max(BW2(:))/2;

    i=find(BW2<=thresh);

    BW2(i)=0;

    i=find(BW2>thresh);

    BW2(i)=255;

    

    [f,v] = isosurface(BW2,100);

 

    v2=[129-v(:,1) v(:,3)-129 129-v(:,2)]; % in order to cope with the different orientation 

    v=v2;

 

    fprintf('morphological closing done.\n');

    fprintf('writing outer surface...\n');

    

    fname=output_surface;

    vert = v;

    face = f - 1;

    vnum = size(vert,1);  

    fnum = size(face,1);  

    

    % open it as a big-endian file

    fid = fopen(fname, 'wb', 'b');

    TRIANGLE_FILE_MAGIC_NUMBER = 16777214;

    fwrite3(fid, TRIANGLE_FILE_MAGIC_NUMBER);

    

    % Ouput a couple of text lines with creation date

    str = sprintf('created from matlab on %s\n',datestr(now));

    fwrite(fid, str,'char');

    fwrite(fid, vnum,'int32');

    fwrite(fid, fnum,'int32');

 

    % reshape vert into column array and write

    vert = reshape(vert',size(vert,1)*size(vert,2),1);

    fwrite(fid, vert,'float32');

 

    % reshape face into column array and write

    face = reshape(face',size(face,1)*size(face,2),1);

    fwrite(fid, face,'int32');

    fclose(fid) ;

return

end

=================[Matlab Script: make_roi_paths]===============

function make_roi_paths(pial, outersmoothed, radius, stepsize, outdir, flagf)

 

% Original Author: Marie Schaer

% Date: 2007/11/15

%

% This function creates multiple circular regions of interest on the outer

% surface, measures their area (saved in matlab format) and save path file

% containing a set of

% corresponding points on the pial surface. The path files are required as

% an input to the mri_path2label function

%

% Parameters:

% radius: radius of the circular region of interest on the outer surface.

% Most appropriate values between 20 and 25 mm, further documentation

% available in the validation publication: "A surface-based approach to

% quantify local cortical gyrification" by Schaer M, Bach Cuadra M,

% Tamarit L, Eliez E, Thiran JP, IEEE Transactions on Medical Imaging, 2007

%

% Utilities to transfer points on the pial surface calls "mesh_vertex_nearest",

% a function written by Darren Weber. Utilities to compute matrix adjacency

% for the pial surface calls "mesh_adjacency", a program first written by

% Darren Weber, then modified by G. Peyre (source: toolbox_graph on

% mathworks)

%

% Example: make_roi_paths('lh.pial','lh.outer-pial-smoothed',25,100,/tmp)

 

t0 = cputime;

 

disp('loading datas ...')

[mesh_total.vertices, mesh_total.faces]=freesurfer_read_surf(pial);

[mesh_outer.vertices, mesh_outer.faces]=freesurfer_read_surf(outersmoothed);

 

disp('preparing outer mesh structure ...')

mesh_outer = createMeshFacesOfVertex (mesh_outer.vertices, mesh_outer.faces)

 

disp('preparing pial mesh structure ...')

A = mesh_adjacency(mesh_total);

 

% Keep track of the area for each individual region of interest on the

% outer surface:

nbVertices = size(mesh_outer.vertices,1);

areasOuterROI = zeros(nbVertices,1);

areasOuterROI = sparse (areasOuterROI);

clear nbVertices

 

 

% Circular regions of interest are defined each 100 vertex on the outer

% surface. Due to the high resolution of the outer mesh (average face area

% of ~0.3mm2), calculations are computed each 1 on 100 vertex, to avoid

% high redundancies and optimize calculation time.

 

for iV = 1:stepsize: length(mesh_outer.vertices)

 

    disp(['... creating path file for vertex ',num2str(iV),' / ',num2str(length(mesh_outer.vertices))])

 

    % ------ Part 1: find the ROI on the enveloppe --------------

 

    [verticesInROI, facesInROI] = getVerticesAndFacesInSphere(mesh_outer, iV, radius);

 

    % In rare case, the intersection of the enveloppe with the sphere result

    % in two regions of interest: the one that we expect, circular and

    % centered in iV; and an aberrant one at the bottom of the sphere

    % (e.g. near the superior sagittal vault). The next step is used to

    % control for that and keep only the radial region centered at iV.

    facesListOuterGeo = MakeGeodesicOuterROI (mesh_outer, iV, facesInROI);

 

    % Find the perimeter of the outer ROI

    verticesOuterGeo=unique(facesListOuterGeo(:));

    perim=setdiff(verticesOuterGeo,verticesInROI);

 

    % Measure its area

    areaOuterROI = getFacesArea(mesh_outer,facesListOuterGeo);

 

    areasOuterROI(iV) = areaOuterROI;

 

    if (mod(iV,50*stepsize) == 1)

        fname = [outdir '/' pial '.area_intermed.mat'];

        save(fname, 'areasOuterROI') ;

        disp(['area file for outer ROIs saved at ',num2str(iV)])

    end

 

    clear verticesInROI facesInROI facesListOuterGeo verticesOuterGeo areaOuterROI

 

 

    % ----- Part2: Define the corresponding set of points on the pial surface:

    % Transfer a few points of the perimeter of the hulls ROI from

    % the envelope to the pial surface, and save them as path file.

    step = 7;

    [verticeslist] = SearchProjectionOnPial(mesh_total, mesh_outer, perim, step);

 

    % reorganize the set of points in the right order to input them to

    % mri_path2label

    reorglist = reorganize_verticeslist (mesh_total, A, mesh_outer, perim, verticeslist, step);

    lindex = reorglist';

    lxyz = mesh_total.vertices(lindex,:); lindex= lindex-1;

 

 

    p = sprintf ('%06d', iV);

    write_path(lxyz, lindex,[outdir '/' pial '.' p '.path']);

 

    clear verticeslist lindex lxyz

 

end

 

disp(['saving area file in matlab format... '])

fname = [outdir '/' pial '.outer_ROIs_area.mat'];

save(fname, 'areasOuterROI')

 

disp('DONE.')

deltaT = cputime - t0

 

% indicate successful completion by deleting flagfile

delete(flagf);

=======================[find corresponding_center_FSformat]=========================

t0 = cputime;    

 

[mesh_pial.vertices, mesh_pial.faces] = freesurfer_read_surf(pial);

[mesh_outer.vertices, mesh_outer.faces] = freesurfer_read_surf(outersmoothed);

 

fidv = fopen([outdir '/' pial '.center.vertices'], 'w') ;

 

for iV = 1 :stepsize: length(mesh_outer.vertices) 

        

    centerSeed=mesh_vertex_nearest(mesh_pial.vertices,mesh_outer.vertices(iV,:));

    centerSeed = centerSeed - 1; %FreeSurfer-s vertex are 0-based

    p = sprintf ('%06d', iV);

      

    fid = fopen([outdir '/' pial '.' p '.center'], 'w') ;

    fprintf(fid,'%d\n',centerSeed);

    fclose(fid) ;

 

    fprintf(fidv,'%s\n',p);

end

    

fclose(fidv) ;

 

deltaT = cputime - t0

 

% indicate successful completion by deleting flagfile

delete(flagfile);

===================================================================================