% mrisegment() - Segment anatomical MRI images using Matlab segmentation % routine isosurface() % % Usage: % [mesh norm vol vol2] = mrisegment(mri, 'key', 'val', ...); % % Inputs: % mri - MRI volume or anatomical MRI filename. % If a filename, the function first tries to read it as an AFNI file % (do not enter the file extension) if the AFNI Matlab librairy % is present (http://afni.nimh.nih.gov/afni/matlab/). % Then it tries to read it using the read_fcdc_mri() function of the % Fieldtrip toolbox (this also requires SPM2 (not SPM) to be % instaled). % % Optional inputs: % 'subsampl' - subsampling factor (default is 2) % 'subsampl2' - subsampling after smoothing (default is none) % 'smooth' - 3-D smoothing in pixel (default is 3) % 'thresh' - value for extracting the surface (default is 30) % 'plot' - ['on'|'off'] plot 3-D mesh. default is 'off'. % % Outputs: % mesh - matlab mesh with fields 'vertices' and 'faces' % norm - normals to the mesh (for plotting) % vol - original MRI volume % vol2 - subsampled and smoothed MRI volume % % Author: Arnaud Delorme, Salk, SCCN, UCSD, CA, March 22, 2004 %123456789012345678901234567890123456789012345678901234567890123456789012 % Copyright (C) 2004 Arnaud Delorme % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA function [fv, in, V, D] = mrisegment(filename, varargin); coordinc = 2; % 2mm if nargin < 1 [filename filepath] = uigetfile('*.HEAD;*.head', 'Pick an AFNI header file'); if isnumeric(filename), return; end; filename = [filepath filename(1:end-5)]; end; %load('/home/arno/matlab/MNI/VolumeMNI.mat', '-mat'); result = 1; if isstr(filename) try, try, addpath('/home/arno/matlab/afni/'); catch, end; disp('Loading Brikfile...'); [err, V, Info, ErrMessage] = BrikLoad(filename); planes = Info.ORIENT_SPECIFIC; %set to 0 3 4 on write. The three numbers refer to x,y,z and have the following values: %{RL=0, LR, PA, AP, IS, SI} where %R=right, P=posterior, I=inferior, etc. rl = find(planes == 0 | planes == 1); ap = find(planes == 2 | planes == 3); is = find(planes == 4 | planes == 5); V = permute(V, [rl ap is]); if planes(rl) == 1, V = flipdim(V, 1); end; % lr if planes(ap) == 3, V = flipdim(V, 2); end; % ap if planes(is) == 5, V = flipdim(V, 3); end; % si catch, try, addpath('/home/duann/matlab/spm2'); catch, end; disp('Cannot load BRIK file'); disp(lasterr); V = read_fcdc_mri(filename); V = V.anatomy; end; else V = filename; end; if nargin < 2 uilist = { { 'style' 'text' 'string' 'Subsample (integer)' } ... { 'style' 'edit' 'string' '2' } { } ... { 'style' 'text' 'string' 'Smoothing (pixel)' } ... { 'style' 'edit' 'string' '3' } { } ... { 'style' 'text' 'string' 'Threshold color' }, ... { 'style' 'edit' 'string' '30' } ... { 'style' 'pushbutton' 'string' 'Histogram' ... 'callback' 'figure; tmp = get(gcbf, ''userdata''); hist(tmp(:), 100); clear tmp;' } ... { 'style' 'text' 'string' 'Plot now' } ... { 'style' 'checkbox' 'string' '' 'value' 0 } {} }; uigeom = { [ 2 1 1 ] [ 2 1 1 ] [ 2 1 1 ] [ 2 0.4 1.6 ] }; result = inputgui( uigeom, uilist, '', 'Head mesh generation parameters', V); options = { 'subsampl' str2num( result{1} ) 'smooth' str2num( result{2} ) 'thresh' str2num( result{3} ) }; if result{4}, options = { options{:} 'plot' 'on' }; end; else options = varargin; end; g = finputcheck(options, { 'subsampl' 'integer' [1 Inf] 2; 'subsampl2' 'integer' [1 Inf] 1; 'smooth' 'real' [0 Inf] 3; 'thresh' 'real' [0 Inf] 30; 'plot' 'string' { 'on' 'off' } 'off' }); if isstr(g), error(g); end; disp('Generating head mesh'); disp('Sub-sampling MRI data...'); %coordx = [coordinc/2:coordinc:size(V,1)*coordinc]-size(V,1)/2*coordinc; %coordy = [coordinc/2:coordinc:size(V,2)*coordinc]-size(V,2)/2*coordinc; %coordz = [coordinc/2:coordinc:size(V,3)*coordinc]-size(V,3)/2*coordinc; coordx = [0.5:1:size(V,1)]-size(V,1)/2; coordy = [0.5:1:size(V,2)]-size(V,2)/2; coordz = [0.5:1:size(V,3)]-size(V,3)/2; D = reducevolume(V, [g.subsampl g.subsampl g.subsampl]); coordx = coordx(1:g.subsampl:end); coordy = coordy(1:g.subsampl:end); coordz = coordz(1:g.subsampl:end); %coordx = repmat(reshape(coordx, [size(D,1) 1 1]), [1 size(D,2) size(D,3)]); %coordy = repmat(reshape(coordy, [1 size(D,2) 1]), [size(D,1) 1 size(D,3)]); %coordz = repmat(reshape(coordz, [1 1 size(D,3)]), [size(D,1) size(D,2) 1]); disp('Smoothing MRI data...'); D = smooth3(D, 'gaussian', g.smooth); D = reducevolume(D, [g.subsampl2 g.subsampl2 g.subsampl2]); coordx = coordx(1:g.subsampl2:end); coordy = coordy(1:g.subsampl2:end); coordz = coordz(1:g.subsampl2:end); [coordx coordy coordz] = meshgrid(coordy, coordx, coordz); %D = smooth3(D, 'gaussian', g.smooth); disp('Computing Head surface...'); fv = isosurface(coordx, coordy, coordz, D, g.thresh); disp('Compute normals to mesh'); in = isonormals(coordx, coordy, coordz, D, fv.vertices); disp('Use plot head mesh to plot head mesh with electrodes'); if strcmpi(g.plot, 'on') figure; hiso = patch(fv, ... 'facecolor', [1,.75,.65], ... 'edgecolor', 'none', 'vertexnormals', in); view([45 30]); axis tight; axis equal %daspet([1 1 .4]); lightangle(45,30); set(gcf, 'renderer', 'zbuffer'); %set(hcap, 'ambientstrength', .6); set(hiso, 'specularcolorreflectance', 0, 'specularexponent',100); %, 'specularstrength', 0); end;