From 247e2e8ce1f43ee390ceb32e7e02460bd9463412 Mon Sep 17 00:00:00 2001 From: Robert Seymour Date: Mon, 8 Oct 2018 13:57:38 +1100 Subject: [PATCH] Fix MEMES3 blurb --- MEMES3.m | 232 +++++++++++++++++++++++++++---------------------------- 1 file changed, 113 insertions(+), 119 deletions(-) diff --git a/MEMES3.m b/MEMES3.m index 2619f48..81faf26 100644 --- a/MEMES3.m +++ b/MEMES3.m @@ -48,7 +48,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % Example function call: % MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,... -% mesh_library,initial_mri_realign,'','best',[0.98:0.01:1.02],8,'no') +% '','best',[0.98:0.01:1.02],8,'no') % I have introduced a variable scaling parameter for the MRIs to % help with coregistration. For example to apply -2% to +2% scaling to @@ -56,9 +56,6 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % % However NOTE: the more scaling factors you apply the longer it will take % -% This script is specific for adult MEG-HCP data (95 MRIs, headmodels and -% sourcemodels), but could easily be adapted for other datasets. -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fprintf('\nThis is MEMES v0.3\n\nMake sure you have asked Robert for an MRI library\n\n'); @@ -110,17 +107,17 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad dirFlags = [files.isdir]; % Extract only those that are directories. subFolders = files(dirFlags); - + % Now these names to a variable called subject subject = []; - + for sub = 1 : length(subFolders) subject{sub} = subFolders(sub).name; end - + fprintf('%d subjects found in the MRI library: from %s to %s\n',... length(subject),subject{1}, subject{end}); - + catch warning('Something is wrong with your MRI library... Check the path!\n'); end @@ -134,7 +131,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad load([path_to_MRI_library subject{1} '/sourcemodel3d_8mm.mat']); clear mesh headmodel mri_realigned sourcemodel3d fprintf('...Subject %s is organised correctly!\n',subject{1}); - + catch warning('Your MRI library is not organised correctly'); disp('Each folder should contain: mesh.mat, headmodel.mat, mri_realigned.mat, sourcemodel3d_8mm.mat'); @@ -166,34 +163,34 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad markers = mrk.fid.pos([2 3 1 4 5],:);%reorder mrk to match order in shape [R,T,Yf,Err] = rot3dfit(markers,shape.fid.pnt(4:end,:));%calc rotation transform meg2head_transm = [[R;T]'; 0 0 0 1];%reorganise and make 4*4 transformation matrix - + disp('Performing re-alignment'); grad_trans = ft_transform_geometry_PFS_hacked(meg2head_transm,grad_con); %Use my hacked version of the ft function - accuracy checking removed not sure if this is good or not grad_trans.fid = shape; %add in the head information save grad_trans grad_trans - + % Else if there is a bad marker else fprintf(''); disp('TAKING OUT BAD MARKER(S)'); - + badcoilpos = []; - + % Identify the bad coil for num_bad_coil = 1:length(bad_coil) pos_of_bad_coil = find(ismember(shape.fid.label,bad_coil{num_bad_coil}))-3; badcoilpos(num_bad_coil) = pos_of_bad_coil; end - + % Re-order mrk file to match elp file markers = mrk.fid.pos([2 3 1 4 5],:);%reorder mrk to match order in shape % Now take out the bad marker(s) when you realign markers(badcoilpos,:) = []; - + % Get marker positions from elp file - fids_2_use = shape.fid.pnt(4:end,:); + fids_2_use = shape.fid.pnt(4:end,:); % Now take out the bad marker(s) when you realign fids_2_use(badcoilpos,:) = []; - + % If there are two bad coils use the ICP method, if only one use % rot3dfit as usual disp('Performing re-alignment'); @@ -242,14 +239,14 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % For each subject... for m = 1:length(subject) - + % Load the mesh load([path_to_MRI_library subject{m} '/mesh.mat']) - + numiter = 30; count2 = 1; - + trans_matrix_temp = []; error_2 = []; - + % Perform ICP fit with different scaling factors for scale = scaling fprintf('Completed iteration %d of %d ; %d of %d MRIs\n',count2,length(scaling),m,length(subject)); @@ -260,19 +257,19 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad trans_matrix_temp{count2} = inv([real(R) real(t);0 0 0 1]); count2 = count2+1; end - + % Find scaling factor with smallest error min_error = min(error_2); % Add error to error_term error_term(m) = min_error; - + % Add transformation matrix to trans_matrix_library trans_matrix_library{m} = trans_matrix_temp{find(error_2==min_error)}; % Add scaling factor scaling_factor_all(m) = scaling(find(error_2==min_error)); - + fprintf('Best scaling factor is %.2f\n',scaling(find(error_2==min_error))); - + % Clear mesh for next loop clear mesh end @@ -299,7 +296,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad 0 scaling_factor_all(concat(i)) 0 0; ... 0 0 scaling_factor_all(concat(i)) 0; 0 0 0 1],mesh_spare.pos); mesh_spare.pos = ft_warp_apply(trans_matrix_library{(concat(i))}, mesh_spare.pos); - + subplot(3,3,i) ft_plot_mesh(mesh_spare,'facecolor',[238,206,179]./255,'EdgeColor','none','facealpha',0.8); hold on; camlight; hold on; view([-270,-10]); @@ -310,11 +307,11 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad elseif ismember(i,7:9) title(sprintf('WORST: %d', error_term((concat(i))))); end - + ft_plot_headshape(headshape_downsampled); - + clear mesh - + if i == 9 print('best_middle_worst_examples','-dpng','-r100'); end @@ -327,22 +324,22 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad figure;hist(scaling_factor_all,length(scaling)); ylabel('Count'); xlabel('Scaling Parameter'); - + % Get information about the same % histogram by returning arguments [n,x] = hist(scaling_factor_all,5); % Create strings for each bar count barstrings = num2str(n'); - + barstrings2 = num2str(scaling'); - + % Create text objects at each location ylim([0 max(n)+5]); text(x,n,barstrings,'horizontalalignment','center','verticalalignment','bottom'); - + xticks(scaling); xTick = get(gca,'xtick'); - + h = findobj(gca,'Type','patch'); h.FaceColor = [0 0.5 0.5]; h.EdgeColor = 'w'; @@ -498,18 +495,18 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % % % fprintf('\nCOMPLETED - check the output for quality control\n'); - - + + case 'best' - + % Find the MRI with the lowest ICP error between Polhemus points % and 3D scalp mesh winner = find(error_term == min(min(error_term))); fprintf('\nThe winning MRI is number %d of %d : %s\n',winner,length(subject),subject{winner}); - + % Get the transformation matrix of the winner trans_matrix = trans_matrix_library{winner}; - + % Get facial mesh of winner load([path_to_MRI_library subject{winner} '/mesh.mat']) mesh_spare = mesh; @@ -517,50 +514,50 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad scaling_factor_all(winner) 0 0; 0 0 scaling_factor_all(winner) 0;... 0 0 0 1],mesh_spare.pos); mesh_spare.pos = ft_warp_apply(trans_matrix, mesh_spare.pos); - + % Get MRI of winning subject mri_realigned = load([path_to_MRI_library... subject{winner} '/mri_realigned.mat']); - + %% Create Headmodel (in mm) fprintf(' Creating Headmodel in mm\n'); - + load([path_to_MRI_library subject{winner} '/headmodel.mat']); - + % Scale headmodel.bnd.pos = ft_warp_apply([scaling_factor_all(winner) 0 0 0;0 ... scaling_factor_all(winner) 0 0; 0 0 scaling_factor_all(winner) 0; 0 0 0 1],... headmodel.bnd.pos); - + % Transform (MESH --> coreg via ICP adjustment) headmodel.bnd.pos = ft_warp_apply(trans_matrix,headmodel.bnd.pos); - + figure; ft_plot_vol(headmodel); ft_plot_headshape(headshape_downsampled); - + %% Create Sourcemodel (in mm) fprintf('Creating an %dmm Sourcemodel in mm\n',sourcemodel_size); - - % Load specified sized sourcemodel + + % Load specified sized sourcemodel load([path_to_MRI_library ... subject{winner} '/sourcemodel3d_' num2str(sourcemodel_size) 'mm.mat']); - + % Scale sourcemodel3d.pos = ft_warp_apply([scaling_factor_all(winner) 0 0 0;0 scaling_factor_all(winner) 0 0; 0 0 scaling_factor_all(winner) 0; 0 0 0 1],sourcemodel3d.pos); - + % Transform (MESH --> coreg via ICP adjustment) sourcemodel3d.pos = ft_warp_apply(trans_matrix,sourcemodel3d.pos); - + % Create figure to check headodel and sourcemodel match figure; ft_plot_vol(headmodel, 'facecolor', 'cortex', 'edgecolor', 'none'); alpha 0.4; camlight; ft_plot_mesh(sourcemodel3d.pos(sourcemodel3d.inside,:),'vertexsize',5); view([0 0]); - + view_angle = [0 90 180 270]; - + % Create figure to show final coregiration figure; hold on; for rep = 1:4 @@ -573,9 +570,9 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad ft_plot_mesh(mesh_spare,'facecolor',[238,206,179]./255,'EdgeColor','none','facealpha',0.5); camlight; lighting phong; material dull; end - + print('coregistration_volumetric_quality_check','-dpng','-r100'); - + % %% Create coregistered 3D cortical mesh % mesh = ft_read_headshape({[path_to_MRI_library ... % subject{winner} '/MEG/anatomy/' subject{winner} '.L.midthickness.4k_fs_LR.surf.gii'],... @@ -599,19 +596,19 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % ft_plot_headshape(headshape_downsampled) %plot headshape % ft_plot_mesh(mesh,'facealpha',0.8); camlight; hold on; view([100 4]); % print('headmodel_3D_cortical_mesh_quality','-dpng'); - + %% SAVE fprintf('\nSaving the necessary data\n'); - + save headmodel headmodel save trans_matrix trans_matrix save grad_trans grad_trans save sourcemodel3d sourcemodel3d % save mesh mesh - - + + fprintf('\nCOMPLETED - check the output for quality control\n'); - + otherwise fprintf('Something went wrong - did you specify *average* or *best*') end @@ -628,66 +625,66 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [shape] = parsePolhemus(elpfile,hspfile) - + fid1 = fopen(elpfile); C = fscanf(fid1,'%c'); fclose(fid1); - + E = regexprep(C,'\r','xx'); E = regexprep(E,'\t','yy'); - + returnsi = strfind(E,'xx'); tabsi = strfind(E,'yy'); sensornamesi = strfind(E,'%N'); fiducialsstarti = strfind(E,'%F'); lastfidendi = strfind(E(fiducialsstarti(3):fiducialsstarti(length(fiducialsstarti))+100),'xx'); fiducialsendi = fiducialsstarti(1)+strfind(E(fiducialsstarti(1):fiducialsstarti(length(fiducialsstarti))+lastfidendi(1)),'xx'); - + NASION = E(fiducialsstarti(1)+4:fiducialsendi(1)-2); NASION = regexprep(NASION,'yy','\t'); NASION = str2num(NASION); - + LPA = E(fiducialsstarti(2)+4:fiducialsendi(2)-2); LPA = regexprep(LPA,'yy','\t'); LPA = str2num(LPA); - + RPA = E(fiducialsstarti(3)+4:fiducialsendi(3)-2); RPA = regexprep(RPA,'yy','\t'); RPA = str2num(RPA); - + LPAredstarti = strfind(E,'LPAred'); LPAredendi = strfind(E(LPAredstarti(1):LPAredstarti(length(LPAredstarti))+45),'xx'); LPAred = E(LPAredstarti(1)+11:LPAredstarti(1)+LPAredendi(2)-2); LPAred = regexprep(LPAred,'yy','\t'); LPAred = str2num(LPAred); - + RPAyelstarti = strfind(E,'RPAyel'); RPAyelendi = strfind(E(RPAyelstarti(1):RPAyelstarti(length(RPAyelstarti))+45),'xx'); RPAyel = E(RPAyelstarti(1)+11:RPAyelstarti(1)+RPAyelendi(2)-2); RPAyel = regexprep(RPAyel,'yy','\t'); RPAyel = str2num(RPAyel); - + PFbluestarti = strfind(E,'PFblue'); PFblueendi = strfind(E(PFbluestarti(1):PFbluestarti(length(PFbluestarti))+45),'xx'); PFblue = E(PFbluestarti(1)+11:PFbluestarti(1)+PFblueendi(2)-2); PFblue = regexprep(PFblue,'yy','\t'); PFblue = str2num(PFblue); - + LPFwhstarti = strfind(E,'LPFwh'); LPFwhendi = strfind(E(LPFwhstarti(1):LPFwhstarti(length(LPFwhstarti))+45),'xx'); LPFwh = E(LPFwhstarti(1)+11:LPFwhstarti(1)+LPFwhendi(2)-2); LPFwh = regexprep(LPFwh,'yy','\t'); LPFwh = str2num(LPFwh); - + RPFblackstarti = strfind(E,'RPFblack'); RPFblackendi = strfind(E(RPFblackstarti(1):end),'xx'); RPFblack = E(RPFblackstarti(1)+11:RPFblackstarti(1)+RPFblackendi(2)-2); RPFblack = regexprep(RPFblack,'yy','\t'); RPFblack = str2num(RPFblack); - + allfids = [NASION;LPA;RPA;LPAred;RPAyel;PFblue;LPFwh;RPFblack]; fidslabels = {'NASION';'LPA';'RPA';'LPAred';'RPAyel';'PFblue';'LPFwh';'RPFblack'}; - + fid2 = fopen(hspfile); C = fscanf(fid2,'%c'); fclose(fid2); @@ -695,18 +692,18 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad E = regexprep(E,'\t','yy'); %replace tabs with "yy" returnsi = strfind(E,'xx'); tabsi = strfind(E,'yy'); - + headshapestarti = strfind(E,'position of digitized points'); headshapestartii = strfind(E(headshapestarti(1):end),'xx'); headshape = E(headshapestarti(1)+headshapestartii(2)+2:end); headshape = regexprep(headshape,'yy','\t'); headshape = regexprep(headshape,'xx',''); headshape = str2num(headshape); - + shape.pnt = headshape; shape.fid.pnt = allfids; shape.fid.label = fidslabels; - + %convert to BESA style coordinates so can use the .pos file or sensor %config from .con % shape.pnt = cat(2,fliplr(shape.pnt(:,1:2)),shape.pnt(:,3)).*1000; @@ -721,7 +718,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % shape.fid.pnt(:,2) = neg2; shape.unit='m'; % shape = ft_convert_units(shape,'cm'); - + new_name2 = ['shape.mat']; save (new_name2,'shape'); end @@ -741,7 +738,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % % rot3dfit: Frank Evans, NHLBI/NIH, 30 November 2001 % - + % ROT3DFIT uses the method described by K. S. Arun, T. S. Huang,and % S. D. Blostein, "Least-Squares Fitting of Two 3-D Point Sets", % IEEE Transactions on Pattern Analysis and Machine Intelligence, @@ -754,42 +751,42 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % % Special cases, e.g. colinear and coplanar points, are not % implemented. - + %error(nargchk(2,2,nargin)); narginchk(2,2); %PFS Change to update if size(X,2) ~= 3, error('X must be N x 3'); end; if size(Y,2) ~= 3, error('Y must be N x 3'); end; if size(X,1) ~= size(Y,1), error('X and Y must be the same size'); end; - + % mean correct - + Xm = mean(X,1); X1 = X - ones(size(X,1),1)*Xm; Ym = mean(Y,1); Y1 = Y - ones(size(Y,1),1)*Ym; - + % calculate best rotation using algorithm 12.4.1 from % G. H. Golub and C. F. van Loan, "Matrix Computations" % 2nd Edition, Baltimore: Johns Hopkins, 1989, p. 582. - + XtY = (X1')*Y1; [U,S,V] = svd(XtY); R = U*(V'); - + % solve for the translation vector - + T = Ym - Xm*R; - + % calculate fit points - + Yf = X*R + ones(size(X,1),1)*T; - + % calculate the error - + dY = Y - Yf; Err = norm(dY,'fro'); % must use Frobenius norm end function [output] = ft_transform_geometry_PFS_hacked(transform, input) - + % FT_TRANSFORM_GEOMETRY applies a homogeneous coordinate transformation to % a structure with geometric information, for example a volume conduction model % for the head, gradiometer of electrode structure containing EEG or MEG @@ -806,7 +803,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % output = ft_transform_geometry(transform, input) % % See also FT_WARP_APPLY, FT_HEADCOORDINATES - + % Copyright (C) 2011, Jan-Mathijs Schoffelen % % This file is part of FieldTrip, see http://www.fieldtriptoolbox.org @@ -826,18 +823,18 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % along with FieldTrip. If not, see . % % $Id: ft_transform_geometry.m$ - + % flg rescaling check allowscaling = ~ft_senstype(input, 'meg'); - + % determine the rotation matrix rotation = eye(4); rotation(1:3,1:3) = transform(1:3,1:3); - + if any(abs(transform(4,:)-[0 0 0 1])>100*eps) error('invalid transformation matrix'); end - + %%### get rid of this accuracy checking below as some of the transformation %%matricies will be a bit hairy### if ~allowscaling @@ -847,19 +844,19 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad %error('only a rigid body transformation without rescaling is allowed'); %end end - + if allowscaling % FIXME build in a check for uniform rescaling probably do svd or so % FIXME insert check for nonuniform scaling, should give an error end - + tfields = {'pos' 'pnt' 'o' 'coilpos' 'chanpos' 'chanposold' 'chanposorg' 'elecpos', 'nas', 'lpa', 'rpa', 'zpoint'}; % apply rotation plus translation rfields = {'ori' 'nrm' 'coilori' 'chanori' 'chanoriold' 'chanoriorg'}; % only apply rotation mfields = {'transform'}; % plain matrix multiplication recfields = {'fid' 'bnd' 'orig'}; % recurse into these fields % the field 'r' is not included here, because it applies to a volume % conductor model, and scaling is not allowed, so r will not change. - + fnames = fieldnames(input); for k = 1:numel(fnames) if ~isempty(input.(fnames{k})) @@ -905,7 +902,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function rmatx=rotate_about_z(deg) - + deg = deg2rad(deg); rmatx = [cos(deg) sin(deg) 0 0;-sin(deg) cos(deg) 0 0;0 0 1 0;0 0 0 1]; end @@ -983,16 +980,16 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad function [headshape_downsampled] = downsample_headshape(path_to_headshape,... numvertices,varargin) - + % If not specified include the facial points if isempty(varargin) include_facial_points = 'yes'; - + else include_facial_points = varargin{1}; end - - + + % Get headshape headshape = ft_read_headshape(path_to_headshape); % Convert to cm @@ -1000,14 +997,14 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad % Convert to BESA co-ordinates % headshape.pos = cat(2,fliplr(headshape.pos(:,1:2)),headshape.pos(:,3)); % headshape.pos(:,2) = headshape.pos(:,2).*-1; - + % Get indices of facial points (up to 3cm above nasion) - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Is 3cm the correct distance? % Possibly different for child system? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + count_facialpoints = find(headshape.pos(:,3)<3); if isempty(count_facialpoints) disp('CANNOT FIND ANY FACIAL POINTS - COREG BY ICP MAY BE INACCURATE'); @@ -1016,20 +1013,20 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad rrr = 1:4:length(facialpoints); facialpoints = facialpoints(rrr,:); clear rrr; end - + % Remove facial points for now headshape.pos(count_facialpoints,:) = []; - + % Create mesh out of headshape downsampled to x points specified in the % function call cfg.numvertices = numvertices; cfg.method = 'headshape'; cfg.headshape = headshape.pos; mesh = ft_prepare_mesh(cfg, headshape); - + % Replace the headshape info with the mesh points headshape.pos = mesh.pos; - + % Create figure for quality checking figure; subplot(2,2,1);ft_plot_mesh(mesh); hold on; title('Downsampled Mesh'); @@ -1044,7 +1041,7 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad title('Downsampled Headshape View 3'); view(180,0); print('headshape_quality','-dpdf'); - + % Add the facial points back in (default) or leave out if user specified % 'no' in function call if strcmp(include_facial_points,'yes') @@ -1058,26 +1055,23 @@ function MEMES3(dir_name,elpfile,hspfile,confile,mrkfile,path_to_MRI_library,bad headshape.pos = headshape.pos; disp('Not adding facial points back into headshape'); end - + %Add in names of the fiducials from the sensor headshape.fid.label = {'NASION','LPA','RPA'}; - + % Convert fiducial points to BESA % headshape.fid.pos = cat(2,fliplr(headshape.fid.pos(:,1:2)),headshape.fid.pos(:,3)); % headshape.fid.pos(:,2) = headshape.fid.pos(:,2).*-1; - + % Plot for quality checking figure;%ft_plot_sens(sensors) %plot channel position : between the 1st and 2nd coils ft_plot_headshape(headshape) %plot headshape view(0,0); print('headshape_quality2','-dpdf'); - + % Export filename headshape_downsampled = headshape; - + end end - - -