diff --git a/toolbox/anatomy/CalcVertexNormals.m b/toolbox/anatomy/CalcVertexNormals.m
new file mode 100644
index 0000000000..f795301551
--- /dev/null
+++ b/toolbox/anatomy/CalcVertexNormals.m
@@ -0,0 +1,138 @@
+function [VertexNormals,VertexArea,FaceNormals,FaceArea]=CalcVertexNormals(FV,FaceNormals) % ,up,vp
+ % Very slightly modified for efficiency and convenience, and added normr. Marc L.
+%% Summary
+%Author: Itzik Ben Shabat
+%Last Update: July 2014
+
+%summary: CalcVertexNormals calculates the normals and voronoi areas at each vertex
+%INPUT:
+%FV - triangle mesh in face vertex structure
+%N - face normals
+%OUTPUT -
+%VertexNormals - [Nv X 3] matrix of normals at each vertex
+%Avertex - [NvX1] voronoi area at each vertex
+%Acorner - [NfX3] slice of the voronoi area at each face corner
+
+%% Code
+
+if nargin < 2 || isempty(FaceNormals)
+ [FaceNormals, FaceArea] = CalcFaceNormals(FV);
+end
+
+% disp('Calculating vertex normals... Please wait');
+% Get all edge vectors
+e0=FV.Vertices(FV.Faces(:,3),:)-FV.Vertices(FV.Faces(:,2),:);
+e1=FV.Vertices(FV.Faces(:,1),:)-FV.Vertices(FV.Faces(:,3),:);
+e2=FV.Vertices(FV.Faces(:,2),:)-FV.Vertices(FV.Faces(:,1),:);
+% Normalize edge vectors
+% e0_norm=normr(e0);
+% e1_norm=normr(e1);
+% e2_norm=normr(e2);
+
+%normalization procedure
+%calculate face Area
+%edge lengths
+de0=sqrt(e0(:,1).^2+e0(:,2).^2+e0(:,3).^2);
+de1=sqrt(e1(:,1).^2+e1(:,2).^2+e1(:,3).^2);
+de2=sqrt(e2(:,1).^2+e2(:,2).^2+e2(:,3).^2);
+l2=[de0.^2 de1.^2 de2.^2];
+
+%using ew to calulate the cot of the angles for the voronoi area
+%calculation. ew is the triangle barycenter, I later check if its inside or
+%outide the triangle
+ew=[l2(:,1).*(l2(:,2)+l2(:,3)-l2(:,1)) l2(:,2).*(l2(:,3)+l2(:,1)-l2(:,2)) l2(:,3).*(l2(:,1)+l2(:,2)-l2(:,3))];
+
+s=(de0+de1+de2)/2;
+%Af - face area vector
+FaceArea=sqrt(max(0, s.*(s-de0).*(s-de1).*(s-de2)));%herons formula for triangle area, could have also used 0.5*norm(cross(e0,e1))
+% if any(~Af) || any(~FaceArea)
+% error('Degenerate faces.');
+% end
+
+%calculate weights
+Acorner=zeros(size(FV.Faces,1),3);
+VertexArea=zeros(size(FV.Vertices,1),1);
+
+% Calculate Vertice Normals
+VertexNormals=zeros([size(FV.Vertices,1) 3]);
+
+% up=zeros([size(FV.Vertices,1) 3]);
+% vp=zeros([size(FV.Vertices,1) 3]);
+for i=1:size(FV.Faces,1)
+ %Calculate weights according to N.Max [1999]
+
+ wfv1=FaceArea(i)/(de1(i)^2*de2(i)^2);
+ wfv2=FaceArea(i)/(de0(i)^2*de2(i)^2);
+ wfv3=FaceArea(i)/(de1(i)^2*de0(i)^2);
+
+ VertexNormals(FV.Faces(i,1),:)=VertexNormals(FV.Faces(i,1),:)+wfv1*FaceNormals(i,:);
+ VertexNormals(FV.Faces(i,2),:)=VertexNormals(FV.Faces(i,2),:)+wfv2*FaceNormals(i,:);
+ VertexNormals(FV.Faces(i,3),:)=VertexNormals(FV.Faces(i,3),:)+wfv3*FaceNormals(i,:);
+ %Calculate areas for weights according to Meyer et al. [2002]
+ %check if the tringle is obtuse, right or acute
+
+ if ew(i,1)<=0
+ Acorner(i,2)=-0.25*l2(i,3)*FaceArea(i)/(e0(i,:)*e2(i,:)');
+ Acorner(i,3)=-0.25*l2(i,2)*FaceArea(i)/(e0(i,:)*e1(i,:)');
+ Acorner(i,1)=FaceArea(i)-Acorner(i,2)-Acorner(i,3);
+ elseif ew(i,2)<=0
+ Acorner(i,3)=-0.25*l2(i,1)*FaceArea(i)/(e1(i,:)*e0(i,:)');
+ Acorner(i,1)=-0.25*l2(i,3)*FaceArea(i)/(e1(i,:)*e2(i,:)');
+ Acorner(i,2)=FaceArea(i)-Acorner(i,1)-Acorner(i,3);
+ elseif ew(i,3)<=0
+ Acorner(i,1)=-0.25*l2(i,2)*FaceArea(i)/(e2(i,:)*e1(i,:)');
+ Acorner(i,2)=-0.25*l2(i,1)*FaceArea(i)/(e2(i,:)*e0(i,:)');
+ Acorner(i,3)=FaceArea(i)-Acorner(i,1)-Acorner(i,2);
+ else
+ ewscale=0.5*FaceArea(i)/(ew(i,1)+ew(i,2)+ew(i,3));
+ Acorner(i,1)=ewscale*(ew(i,2)+ew(i,3));
+ Acorner(i,2)=ewscale*(ew(i,1)+ew(i,3));
+ Acorner(i,3)=ewscale*(ew(i,2)+ew(i,1));
+ end
+ VertexArea(FV.Faces(i,1))=VertexArea(FV.Faces(i,1))+Acorner(i,1);
+ VertexArea(FV.Faces(i,2))=VertexArea(FV.Faces(i,2))+Acorner(i,2);
+ VertexArea(FV.Faces(i,3))=VertexArea(FV.Faces(i,3))+Acorner(i,3);
+
+% %Calculate initial coordinate system
+% up(FV.Faces(i,1),:)=e2_norm(i,:);
+% up(FV.Faces(i,2),:)=e0_norm(i,:);
+% up(FV.Faces(i,3),:)=e1_norm(i,:);
+end
+VertexNormals=normr(VertexNormals);
+
+% %Calculate initial vertex coordinate system
+% for i=1:size(FV.Vertices,1)
+% up(i,:)=cross(up(i,:),VertexNormals(i,:));
+% up(i,:)=up(i,:)/norm(up(i,:));
+% vp(i,:)=cross(VertexNormals(i,:),up(i,:));
+% end
+
+% disp('Finished Calculating vertex normals');
+end
+
+
+function [FaceNormals, FaceArea]=CalcFaceNormals(FV)
+%% Summary
+%Author: Itzik Ben Shabat
+%Last Update: July 2014
+
+%CalcFaceNormals recives a list of vrtexes and Faces in FV structure
+% and calculates the normal at each face and returns it as FaceNormals
+%INPUT:
+%FV - face-vertex data structure containing a list of Vertices and a list of Faces
+%OUTPUT:
+%FaceNormals - an nX3 matrix (n = number of Faces) containng the norml at each face
+%% Code
+% Get all edge vectors
+e0=FV.Vertices(FV.Faces(:,3),:)-FV.Vertices(FV.Faces(:,2),:);
+e1=FV.Vertices(FV.Faces(:,1),:)-FV.Vertices(FV.Faces(:,3),:);
+% Calculate normal of face
+FaceNormalsA=cross(e0,e1);
+FaceArea = sqrt(sum(FaceNormalsA.^2, 2)) / 2;
+FaceNormals=normr(FaceNormalsA);
+end
+
+function V = normr(V)
+ V = bsxfun(@rdivide, V, sqrt(sum(V.^2, 2)));
+end
+
diff --git a/toolbox/anatomy/SurfaceSmooth.m b/toolbox/anatomy/SurfaceSmooth.m
new file mode 100755
index 0000000000..90be3eb749
--- /dev/null
+++ b/toolbox/anatomy/SurfaceSmooth.m
@@ -0,0 +1,544 @@
+function V = SurfaceSmooth(Surf, Faces, VoxSize, DisplTol, IterTol, Freedom, Verbose)
+ % Smooth closed triangulated surface to remove "voxel" artefacts.
+ %
+ % V = SurfaceSmooth(Surf, [], VoxSize, DisplTol, IterTol, Freedom, Verbose)
+ % V = SurfaceSmooth(Vertices, Faces, VoxSize, DisplTol, IterTol, Freedom, Verbose)
+ %
+ % Smooth a triangulated surface, trying to optimize getting rid of blocky
+ % voxel segmentation artefacts but still respect initial segmentation.
+ % This is achieved by restricting vertex displacement along the surface
+ % normal to half the voxel size, and by compensating normal displacement
+ % of a voxel by an opposite distributed shift in its neighbors. That
+ % way, the total normal displacement is approximately zero (before
+ % potential restrictions are applied).
+ %
+ % Tangential motion is currently left unrestricted, which means the mesh
+ % will readjust over many iterations, much more than necessary to obtain
+ % a smoothed surface. On the other hand, this produces a more uniform
+ % triangulation, which may be desirable in some cases, e.g. after a
+ % reducepatch operation. This tangential motion may also make the normal
+ % restriction a bit less accurate. This all depends on how irregular the
+ % mesh was to start with. To avoid long running times and some of the
+ % tangential deformation, DisplTol and IterTol can be used to limit the
+ % number of iterations.
+ %
+ % To separate these two effects (normal smoothing and tangential mesh
+ % uniformization), the function can be run to achieve each type of motion
+ % separately, by setting Freedom to the appropriate value (see below).
+ % Even if both effects are desired, but precision in the normal
+ % displacement restriction is preferred over running time, I would
+ % suggest running twice (norm., tang.) or three times (tang., norm.,
+ % tang.), but only once in the normal direction. Note that tangential
+ % motion is not perfect and may cause a small amount of smoothing as
+ % well.
+ %
+ % Input variables:
+ % Surf: Instead of Vertices and Faces, a single structure can be given
+ % with fields 'Vertices' and 'Faces' (lower-case v, f also work). In this
+ % case, leave Faces empty [].
+ % Vertices [nV, 3]: Point 3d coordinates.
+ % Faces [nF, 3]: Triangles, i.e. 3 point indices.
+ % VoxSize (default inf): Length of voxels, this determines the amount of
+ % smoothing. For a voxel size of 1, vertices are allowed to move only
+ % 0.5 units in the surface normal direction. This is somewhat optimal
+ % for getting rid of artefacts: it allows steps to become flat even at
+ % shallow angles and a single voxel cube would be transformed to a
+ % sphere of identical voxel volume.
+ % DisplTol (default 0.01*VoxSize): Once the maximum displacement of
+ % vertices is less than this distance, the algorithm stops. If two
+ % values are given, e.g. [0.01, 0.01], the second value is compared to
+ % normal displacement only. The first limit encountered stops
+ % iterating. This allows stopping earlier if only smoothing is
+ % desired and not mesh uniformity.
+ % IterTol (default 100): If the algorithm did not converge, it will stop
+ % after this many iterations.
+ % Freedom (default 2): Indicate which motion is allowed by the
+ % algorithm with an integer value: 0 for (restricted) normal
+ % smoothing, 1 for (unrestricted) tangential motion to get a more
+ % uniform triangulation, or 2 for both at the same time.
+ % Verbose (default 0): If 1, writes initial and final volumes and
+ % areas on the command line. Also gives the number of iterations and
+ % final displacement when the algorithm converged. (A warning is
+ % always given if convergence was not obtained in IterTol iterations.) If
+ % > 1, details are given at each iteration.
+ %
+ % Output: Modified voxel coordinates [nV, 3].
+ %
+ % Written by Marc Lalancette, Toronto, Canada, 2014-02-04
+ % Volume calculation from divergence theorem idea:
+ % http://www.mathworks.com/matlabcentral/fileexchange/26982-volume-of-a-surface-triangulation
+
+ % Note: Although this seems to work relatively well, it is still very new
+ % and not fully tested. Despite the description above which is what was
+ % intended, the algorithm had the tendency to drive growing oscillations
+ % (from iteration to iteration) on the surface. Thus a basic damping
+ % mechanism was added: I simply multiply each movement by a fraction that
+ % seems to avoid oscillations and still converge rapidly enough.
+
+ % Attempt at damping oscillations. Reduce any movement by a certain
+ % fraction. (Multiply movements by this factor.)
+ DampingFactor = 0.91;
+
+ % Add visualizations to debug.
+ isDebugFigures = true;
+
+ if ~isstruct(Surf)
+ if nargin < 2 || isempty(Faces)
+ error('Faces required as second input or "faces" field of first input.');
+ else
+ SurfV = Surf;
+ clear 'Surf';
+ Surf.Vertices = SurfV;
+ Surf.Faces = Faces;
+ clear SurfV Faces;
+ end
+ else
+ if isfield(Surf, 'faces')
+ Surf.Faces = Surf.faces;
+ Surf = rmfield(Surf, 'faces');
+ end
+ if isfield(Surf, 'vertices')
+ Surf.Vertices = Surf.vertices;
+ Surf = rmfield(Surf, 'vertices');
+ elseif ~isfield(Surf, 'Vertices')
+ error('Surf.Vertices field required when second input is empty.');
+ end
+ end
+ if nargin < 3 || isempty(VoxSize)
+ VoxSize = inf;
+ if nargin < 5 || isempty(IterTol)
+ error(['Unrestricted smoothing (no VoxSize) would lead to a sphere of similar volume, ', ...
+ 'unless limited by the number of iterations.']);
+ end
+ end
+ if nargin < 4 || isempty(DisplTol)
+ DisplTol = 0.01 * VoxSize;
+ end
+ if numel(DisplTol) == 1
+ % Only stop when total displacement reaches the limit, normal displacement alone
+ % isn't checked for stopping.
+ DisplTol = [DisplTol, 0];
+ end
+ if nargin < 5 || isempty(IterTol)
+ IterTol = 100;
+ end
+ if nargin < 6 || isempty(Freedom)
+ Freedom = 2; % 0=norm, 1=tang, 2=both.
+ end
+ if nargin < 7 || isempty(Verbose)
+ Verbose = false;
+ end
+
+ % Verify surface is a triangulation.
+ if size(Surf.Faces, 2) > 3
+ error('SurfaceSmooth only works with a triangulated surface.');
+ end
+
+ % Optimal allowed normal displacement, in units of voxel side length.
+ % Based on turning a single voxel into a sphere of same volume: max
+ % needed displacement is in corner:
+ % sqrt(3)/2 - 1/(4/3*pi)^(1/3) = 0.2457
+ % In middle of face it is rather:
+ % 1/(4/3*pi)^(1/3) - 1/2 = 0.1204
+ % Based on very gentle sloped staircase, it would be 0.5, but for 45
+ % degree steps, we only need cos(pi/4)/2 = 0.3536. So something along
+ % those lines seems like a good compromize. For now try to make steps
+ % completely disappear.
+ MaxNormDispl = 0.5 * VoxSize;
+ % MaxDispl = 2 * VoxSize; % To avoid large scale slow flows tangentially, which could distort.
+
+ nV = size(Surf.Vertices, 1);
+ %nF = size(Surf.Faces, 1);
+
+ % Remove duplicate faces. Not necessary considering we have to use
+ % unique on the edges later anyway.
+ % Surf.Faces = unique(Surf.Faces, 'rows');
+
+ if Verbose
+ if isDebugFigures
+ [FdA, VdA, FN, VN] = CalcVertexNormals(Surf);
+ ViewSurfWithNormals(Surf.Vertices, Surf.Faces, VN, FN, VdA, FdA)
+ else
+ [FdA, VdA, FN] = CalcAreas(Surf);
+ end
+ FaceCentroidZ = ( Surf.Vertices(Surf.Faces(:, 1), 3) + ...
+ Surf.Vertices(Surf.Faces(:, 2), 3) + Surf.Vertices(Surf.Faces(:, 3), 3) ) /3;
+ Pre.Volume = FaceCentroidZ' * (FN(:, 3) .* FdA);
+ Pre.Area = sum(FdA);
+ fprintf('Total enclosed volume before smoothing: %g\n', Pre.Volume);
+ fprintf('Total area before smoothing: %g\n', Pre.Area);
+ end
+
+ % Calculate connectivity matrix.
+
+ % Euler characteristic (2-2handles) = V - E + F
+ % For simple closed surface, E = 3*F/2
+ % disp(nV)
+ % disp(2 + nF/2)
+
+ % This expression works when each edge is found once in each direction,
+ % i.e. as long as all normals are consistently pointing in (or out).
+ % C = sparse(Faces(:), [Faces(:, 2); Faces(:, 3); Faces(:, 1)], true);
+ % Seems users had patches that didn't satisfy this restriction, or had
+ % duplicate faces or had possibly intersecting surfaces with 3 faces
+ % sharing an edge.
+ [Edges, ~, iE] = unique(sort([Surf.Faces(:), ...
+ [Surf.Faces(:, 2); Surf.Faces(:, 3); Surf.Faces(:, 1)]], 2), 'rows'); % [Surf.Faces...] = Edges(iE,:)
+ % Look for boundaries of open surface.
+ isBoundE = false(size(Edges, 1), 1);
+ isBoundV = false(nV, 1);
+ iE = sort(iE);
+ n = 1;
+ for i = 2:numel(iE)
+ if iE(i) ~= iE(i-1)
+ if n == 1
+ % Only one copy, boundary edge.
+ isBoundE(iE(i-1)) = true;
+ else
+ n = 1;
+ end
+ else
+ if n == 2
+ % This makes a 3rd copy of the same edge. Strange surface.
+ isBoundE(iE(i)) = true;
+ end
+ n = n + 1;
+ end
+ end
+ % This was very slow for many edges.
+ % for i = 1:size(Edges, 1)
+ % isBoundE(i) = sum(iE == i) < 2;
+ % end
+ if any(isBoundE)
+ if Verbose
+ warning('Open surface detected. Results may be unexpected.');
+ end
+ isBoundV(Edges(isBoundE, :)) = true;
+ end
+ iBoundV = find(isBoundV);
+ iBulkV = setdiff(1:nV, iBoundV);
+ % Vertex connectivity matrix, does not include diagonal (self)
+ C = sparse(Edges, Edges(:, [2,1]), true); % Fills symetrically 1>2, 2>1
+ %C = C | C';
+ % Logical matrix would be huge, so use sparse. However tests in R2011b
+ % indicate that using logical sparse indexing is sometimes slightly
+ % faster (possibly when using linear indexing) but sometimes noticeably
+ % slower. Seems here using a cell array is better.
+ CCell = cell(nV, 1);
+ CCellBulk = cell(nV, 1);
+ for v = 1:nV
+ CCell{v} = find(C(:, v));
+ CCellBulk{v} = setdiff(CCell{v}, iBoundV);
+ end
+ clear C
+ % Number of connected neighbors at each vertex.
+ % nC = full(sum(C, 1));
+
+ V = Surf.Vertices;
+ LastMaxDispl = [inf, inf]; % total, normal only
+ Iter = 0;
+ NormDispl = zeros(nV, 1);
+ while LastMaxDispl(1) > DisplTol(1) && LastMaxDispl(2) > DisplTol(2) && ...
+ Iter < IterTol
+ Iter = Iter + 1;
+ [~, VdA, ~, N] = CalcVertexNormals(Surf);
+ % Double boundary vertex areas to balance their "pull". But not very precise, depends on boundary shape.
+ VdA(isBoundV) = 2 * VdA(isBoundV);
+ VWeighted = bsxfun(@times, VdA, Surf.Vertices);
+
+ % Moving step. (This is slow.)
+ switch Freedom
+ case 2 % Both.
+ for v = iBulkV
+ % Neighborhood average. Improved to weigh by area element to avoid
+ % tangential deformation based on number of neighbors (e.g. shrinking
+ % towards vertices with fewer neighbors).
+ NeighdA = sum(VdA(CCell{v}));
+ NeighdABulk = sum(VdA(CCellBulk{v}));
+ NeighborAverage = sum(VWeighted(CCell{v}, :) / NeighdA, 1); % / nC(v);
+ % Neighborhood correction displacement along normal. Volume
+ % corresponding to this point's normal movement, distributed (divided)
+ % over neighborhood area + itself, which will be shifted inversely.
+ NormalDisplCorr = (NeighborAverage - Surf.Vertices(v, :)) * N(v, :)' / (NeighdABulk/VdA(v) + 1); % / (nC(v) + 1);
+ % Central point is moved to average of neighbors, but shifted back a
+ % bit as they all will be.
+ V(v, :) = V(v, :) + DampingFactor * ( NeighborAverage - Surf.Vertices(v, :) - NormalDisplCorr * N(v, :) );
+ % Neighbors are shifted a bit too along their own normals, such that
+ % the total change in volume (normal displacement times surface area)
+ % is close to zero.
+ V(CCellBulk{v}, :) = V(CCellBulk{v}, :) - DampingFactor * NormalDisplCorr * N(CCellBulk{v}, :);
+ end
+ case 0 % Normal motion only.
+ for v = iBulkV
+ NeighdA = sum(VdA(CCell{v}));
+ NeighdABulk = sum(VdA(CCellBulk{v}));
+ NeighborAverage = sum(VWeighted(CCell{v}, :) / NeighdA, 1); % / nC(v);
+ % d * a / (b+a) = d / (b/a + 1)
+ NormalDisplCorr = (NeighborAverage - Surf.Vertices(v, :)) * N(v, :)' / (NeighdABulk/VdA(v) + 1); % / (nC(v) + 1);
+ % The vertex should move the projected distance minus the correction.
+ % 1 - a/(b+a) = b/(b+a) = 1/(1+a/b) = (b/a) * 1/(b/a+1)
+ V(v, :) = V(v, :) + DampingFactor * NeighdABulk/VdA(v) * NormalDisplCorr * N(v, :);
+ V(CCellBulk{v}, :) = V(CCellBulk{v}, :) - DampingFactor * NormalDisplCorr * N(CCellBulk{v}, :);
+ end
+ case 1 % Tangential motion only. Unrestricted.
+ for v = iBulkV
+ NeighdA = sum(VdA(CCell{v}));
+ NeighborAverage = sum(VWeighted(CCell{v}, :) / NeighdA, 1); % / nC(v);
+ NormalDisplacement = (NeighborAverage - Surf.Vertices(v, :)) * N(v, :)'; % / (nC(v) + 1);
+ V(v, :) = V(v, :) + DampingFactor * ( (NeighborAverage - Surf.Vertices(v, :)) - NormalDisplacement * N(v, :) );
+ % No compensation among neighbors.
+ end
+ otherwise
+ error('Unrecognized Freedom parameter. Should be 0, 1 or 2.');
+ end
+ % Restricting step.
+ % Displacements along normals (N at last positions, Surf.Vertices), added to
+ % previous normal displacement since we want to restrict total normal
+ % displacement.
+ D = NormDispl + dot((V - Surf.Vertices), N, 2);
+ % New restricted total normal displacement.
+ NormDispl = sign(D) .* min(abs(D), MaxNormDispl);
+ % Amounts to move back if greater than allowed.
+ D = D - NormDispl;
+ Where = abs(D) > DisplTol(1) * 1e-6; % > 0, but ignore precision errors.
+ % Fix.
+ if any(Where)
+ V(Where, :) = V(Where, :) - [D(Where), D(Where), D(Where)] .* N(Where, :);
+ end
+ % New restriction on tangential displacement. [Not implemented.]
+ % MaxDispl
+
+ if Verbose > 1
+ [LastMaxDispl(1), iMax(1)] = max(sqrt( sum((V - Surf.Vertices).^2, 2)) );
+ [LastMaxDispl(2), iMax(2)] = max(abs(dot(V - Surf.Vertices, N, 2)));
+ TangDisplVec = CrossProduct(V - Surf.Vertices, N);
+ [LastMaxDispl(3), iMax(3)] = max(sqrt(TangDisplVec(:,1).^2 + TangDisplVec(:,2).^2 + TangDisplVec(:,3).^2));
+ fprintf('Iter %d: max displ %1.4g at vox %d; norm %1.4g (vox %d); tang %1.4g (vox %d)\n', ...
+ Iter, LastMaxDispl(1), iMax(1), ...
+ sign((V(iMax(2),:) - Surf.Vertices(iMax(2),:)) * N(iMax(2),:)')*LastMaxDispl(2), iMax(2), ...
+ sign(TangDisplVec(iMax(3), 1))*LastMaxDispl(3), iMax(3));
+ % Signs are to see if these are oscillations or translations.
+ else
+ LastMaxDispl(1) = sqrt( max(sum((V - Surf.Vertices).^2, 2)) );
+ LastMaxDispl(2) = max(dot(V - Surf.Vertices, N, 2));
+ end
+ Surf.Vertices = V;
+ end
+
+ if Iter >= IterTol && Verbose
+ warning('SurfaceSmooth did not converge within %d iterations. \nLast max point displacement = %f', ...
+ IterTol, LastMaxDispl(1));
+ elseif Verbose
+ fprintf('SurfaceSmooth converged in %d iterations. \nLast max point displacement = %f\n', ...
+ Iter, LastMaxDispl(1));
+ end
+ if Verbose && IterTol > 0
+ [FdA, ~, FN] = CalcVertexNormals(Surf);
+ FaceCentroidZ = ( Surf.Vertices(Surf.Faces(:, 1), 3) + ...
+ Surf.Vertices(Surf.Faces(:, 2), 3) + Surf.Vertices(Surf.Faces(:, 3), 3) ) /3;
+ Post.Volume = FaceCentroidZ' * (FN(:, 3) .* FdA);
+ Post.Area = sum(FdA);
+ fprintf('Total enclosed volume after smoothing: %g\n', Post.Volume);
+ fprintf('Relative volume change: %g %%\n', ...
+ 100 * (Post.Volume - Pre.Volume)/Pre.Volume);
+ fprintf('Total area after smoothing: %g\n', Post.Area);
+ fprintf('Relative area change: %g %%\n', ...
+ 100 * (Post.Area - Pre.Area)/Pre.Area);
+ end
+
+
+
+end
+
+% Much faster than using the Matlab version.
+function c = CrossProduct(a, b)
+ c = [a(:,2).*b(:,3)-a(:,3).*b(:,2), ...
+ a(:,3).*b(:,1)-a(:,1).*b(:,3), ...
+ a(:,1).*b(:,2)-a(:,2).*b(:,1)];
+end
+
+% ----------------------------------------------------------------------
+function [FaceArea, VertexArea, FaceNormals, VertexNormals] = CalcVertexNormals(Surf)
+ % Get face and vertex normals and areas
+
+ % First, see if Brainstorm function is present.
+ isBst = exist('tess_normals', 'file') == 2;
+
+ % Get areas first.
+ if isBst
+ [FaceArea, VertexArea] = CalcAreas(Surf);
+ % Might need to use connectivity matrix if we get bad normals.
+ [VertexNormals, FaceNormals] = tess_normals(Surf.Vertices, Surf.Faces); % , VertConn
+ else
+ [FaceArea, VertexArea, FaceNormals, VertexNormals] = CalcAreas(Surf);
+ end
+end
+
+% % Calculate dA normal vectors to each vertex.
+% function [N, VdA, FN, FdA] = CalcVertexNormals(S)
+% N = zeros(nV, 3);
+% % Get face normal vectors with length the size of the face area.
+% FNdA = CrossProduct( (S.Vertices(S.Faces(:, 2), :) - S.Vertices(S.Faces(:, 1), :)), ...
+% (S.Vertices(S.Faces(:, 3), :) - S.Vertices(S.Faces(:, 2), :)) ) / 2;
+% % For vertex normals, add adjacent face normals, then normalize. Also
+% % add 1/3 of each adjacent area element for vertex area.
+% FdA = sqrt(FNdA(:,1).^2 + FNdA(:,2).^2 + FNdA(:,3).^2);
+% VdA = zeros(nV, 1);
+% for ff = 1:size(S.Faces, 1) % (This is slow.)
+% N(S.Faces(ff, :), :) = N(S.Faces(ff, :), :) + FNdA([ff, ff, ff], :);
+% VdA(S.Faces(ff, :), :) = VdA(S.Faces(ff, :), :) + FdA(ff)/3;
+% end
+% N = bsxfun(@rdivide, N, sqrt(N(:,1).^2 + N(:,2).^2 + N(:,3).^2));
+% FN = bsxfun(@rdivide, FNdA, FdA);
+% end
+
+% % Calculate areas of faces and vertices.
+% function [FdA, VdA, FN] = CalcAreas(S)
+% % Get face normal vectors with length the size of the face area.
+% FN = CrossProduct( (S.Vertices(S.Faces(:, 2), :) - S.Vertices(S.Faces(:, 1), :)), ...
+% (S.Vertices(S.Faces(:, 3), :) - S.Vertices(S.Faces(:, 2), :)) ) / 2;
+% FdA = sqrt(FN(:,1).^2 + FN(:,2).^2 + FN(:,3).^2); % no sum for speed
+% if nargout > 2
+% FN = bsxfun(@rdivide, FN, FdA);
+% end
+% % For vertex areas, add 1/3 of each adjacent area element.
+% VdA = zeros(size(S.Vertices, 1), 1);
+% for iV = 1:3
+% VdA(s.Faces(:, iV)) = VdA(s.Faces(:, iV)) + FdA / 3;
+% end
+% end
+
+% % Find boundary vertices.
+% function isBound = FindBoundary(Faces)
+% nF = size(Faces, 1);
+% Found = logical(nF);
+% Inside = logical(nF);
+% for f = 1:nF
+% for e = 1:3
+% if Found(Faces(f, e), Faces(f, mod(e, 3)+1))
+% Inside(Faces(f, e), Faces(f, mod(e, 3)+1)) = true;
+% else
+% Found(Faces(f, e), Faces(f, mod(e, 3)+1)) = true;
+% end
+% end
+% end
+% isBound = Found & ~Inside;
+% end
+
+
+function ViewSurfWithNormals(Vertices, Faces, VNorm, FNorm, VArea, FArea)
+ figure;
+ hold on;
+ axis equal;
+
+ % Create surface patch
+ patch('Faces', Faces, 'Vertices', Vertices, ...
+ 'FaceColor', 'grey', 'EdgeColor', 'k', 'FaceAlpha', 0.6);
+
+ % Compute face centers
+ face_centers = mean(reshape(Vertices(Faces', :), [], 3, size(Faces, 2)), 2);
+ face_centers = squeeze(face_centers);
+
+ % Scale normals for visualization
+ normal_length = 0.05 * mean(range(Vertices)); % Adjust scaling as needed
+ FNorm = normal_length * bsxfun(@times, FNorm, FArea) / max(FArea);
+ VNorm = normal_length * bsxfun(@times, VNorm, VArea) / max(VArea);
+
+ % Plot face normals (Red)
+ quiver3(face_centers(:,1), face_centers(:,2), face_centers(:,3), ...
+ FNorm(:,1), FNorm(:,2), FNorm(:,3), ...
+ 'r', 'LineWidth', 1.5, 'MaxHeadSize', 0.5);
+
+ % Plot vertex normals (Blue)
+ quiver3(Vertices(:,1), Vertices(:,2), Vertices(:,3), ...
+ VNorm(:,1), VNorm(:,2), VNorm(:,3), ...
+ 'b', 'LineWidth', 1, 'MaxHeadSize', 0.5);
+
+ % Lighting and view settings
+ camlight;
+ lighting gouraud;
+ view(3);
+end
+
+
+function [FaceArea, VertexArea, FaceNormals, VertexNormals] = CalcAreas(S)
+ % Compute face areas, and divides them to assign vertex areas.
+ % Face areas are split into three parts such that each point is assigned to the
+ % vertex it is closest to, similar to Voronoi diagrams, using triangle
+ % circumcenters.
+
+ nV = size(S.Vertices, 1);
+ nF = size(S.Faces, 1);
+
+ % Extract triangle vertex positions
+ Vertices = reshape(S.Vertices(S.Faces(:), :), [nF, 3, 3]); % nF, 3 (x,y,z), 3 (iV)
+ % Edge vectors, ordered such that each is opposite to the correspondingly indexed
+ % vertex (edge 1 is from v2 to v3, opposite vertex 1).
+ % circshift +1 moves elements to the next (larger) index, so the last element
+ % gets to position 1. So here we're doing, e.g. v3 - v2 in first position.
+ Edges = circshift(Vertices, 1, 3) - circshift(Vertices, -1, 3); % nF, 3 (x,y,z), 3 (iE)
+
+ % Get face normal vectors with length the size of the face area.
+ FaceNormals = CrossProduct(Edges(:,:,3), Edges(:,:,1)) / 2;
+ FaceArea = sqrt(FaceNormals(:,1).^2 + FaceNormals(:,2).^2 + FaceNormals(:,3).^2); % no sum for speed
+ if nargout > 2
+ FaceNormals = bsxfun(@rdivide, FaceNormals, FaceArea);
+ end
+
+ % Edge lengths
+ EdgeSq = squeeze(Edges(:,1,:).^2 + Edges(:,2,:).^2 + Edges(:,3,:).^2); % nF, 3 (iE)
+ % Circumcenter barycentric weights
+ BaryWeights = EdgeSq .* (circshift(EdgeSq, 1, 2) + circshift(EdgeSq, -1, 2) - EdgeSq); % nF, 3 (iE)
+
+ FaceVertAreas = zeros(nF, 3);
+
+ % Process in 4 batches
+ % Logical index for faces that have not been processed
+ isRemain = true(nF, 1);
+ % First three cases: circumcenter outside triangular face, past one of 3 edges
+ % This divides the area into two triangles, and one pentagon (or rectangle if
+ % original face has a right angle).
+ for iLongest = 1:3
+ isDo = BaryWeights(:,iLongest) <= 0;
+ isRemain = isRemain & ~isDo;
+ iSh = circshift(1:3, 1-iLongest); % place iLongest in first position
+ % Two vertices of long edge; areas are triangles with right angle
+ % -1/4 face area * ratio of adjoining short edge to projection of long edge on that short edge.
+ % (show with similar right triangle and edge length ratios) Minus sign because of
+ % "dot product" with vectors always making obtuse angle.
+ FaceVertAreas(isDo,iSh(2)) = 0.25 * EdgeSq(isDo,iSh(3)) .* FaceArea(isDo,:) ./ -sum(Edges(isDo,:,iSh(1)) .* Edges(isDo,:,iSh(3)), 2);
+ FaceVertAreas(isDo,iSh(3)) = 0.25 * EdgeSq(isDo,iSh(2)) .* FaceArea(isDo,:) ./ -sum(Edges(isDo,:,iSh(1)) .* Edges(isDo,:,iSh(2)), 2);
+ % Vertex opposite long edge; area has 4 or 5 sides, just assign remaining area.
+ FaceVertAreas(isDo,iSh(1)) = FaceArea(isDo,:) - FaceVertAreas(isDo,iSh(2)) - FaceVertAreas(isDo,iSh(3));
+ end
+ % Last case: circumcenter is inside face, each region is a quadrilateral
+ % Normalize weights and scale with half face area
+ BaryWeights = BaryWeights ./ (BaryWeights(:,1) + BaryWeights(:,2) + BaryWeights(:,3));
+ FaceVertAreas(isRemain,:) = 1/2 * FaceArea(isRemain,:) .* (circshift(BaryWeights(isRemain,:), 1, 2) + circshift(BaryWeights(isRemain,:), -1, 2));
+
+ % Now sum back to single vertex area list
+ VertexArea = accumarray(S.Faces(:), FaceVertAreas(:), [nV, 1]);
+
+ % Use areas as weights to average face normals, if requested
+ if nargout > 3
+ VertexNormals = zeros(nV, 3);
+ % Have to do each coordinate (x,y,z) sequentially with accumarray
+ for i = 1:3
+ % Weight face normals by face vertex area, which we will normalize at the
+ % end by dividing by total vertex area.
+ WeightedFaceVertN = bsxfun(@times, FaceVertAreas, FaceNormals(:,i));
+ VertexNormals(:,i) = accumarray(S.Faces(:), WeightedFaceVertN(:), [nV, 1]);
+ end
+ % Normalize for the area weights
+ VertexNormals = bsxfun(@rdivide, VertexNormals, VertexArea);
+
+ % Final check no NaN
+ if any(isnan(VertexNormals))
+ error('NaN values in VertexNormals.');
+ end
+ end
+ if any(isnan(VertexArea))
+ error('NaN values in VertexArea.');
+ end
+
+end
diff --git a/toolbox/anatomy/cs_convert.m b/toolbox/anatomy/cs_convert.m
index 250fa170b6..706e4cb158 100644
--- a/toolbox/anatomy/cs_convert.m
+++ b/toolbox/anatomy/cs_convert.m
@@ -16,7 +16,7 @@
% DESCRIPTION: https://neuroimage.usc.edu/brainstorm/CoordinateSystems
% - voxel : X=left>right, Y=posterior>anterior, Z=bottom>top
% Coordinate of the center of the first voxel at the bottom-left-posterior of the MRI volume: (1,1,1)
-% - mri : Same as 'voxel' but in millimeters instead of voxels: mriXYZ = voxelXYZ * Voxsize
+% - mri : Same as 'voxel' but in meters (not mm here) instead of voxels: mriXYZ = voxelXYZ * Voxsize
% - scs : Based on: Nasion, left pre-auricular point (LPA), and right pre-auricular point (RPA).
% Origin: Midway on the line joining LPA and RPA
% Axis X: From the origin towards the Nasion (exactly through)
@@ -122,7 +122,7 @@
scs2captrak = [tCapTrak.R, tCapTrak.T; 0 0 0 1];
end
-% ===== CONVERT SRC => MRI =====
+% ===== CONVERT SRC => MRI (m) =====
% Evaluate the transformation to apply
switch lower(src)
case 'voxel'
@@ -174,7 +174,7 @@
error(['Invalid coordinate system: ' src]);
end
-% ===== CONVERT MRI => DEST =====
+% ===== CONVERT MRI (m) => DEST =====
% Evaluate the transformation to apply
switch lower(dest)
case 'voxel'
diff --git a/toolbox/anatomy/mri_histogram.m b/toolbox/anatomy/mri_histogram.m
index 6d985006fd..579008bfa5 100644
--- a/toolbox/anatomy/mri_histogram.m
+++ b/toolbox/anatomy/mri_histogram.m
@@ -23,11 +23,11 @@
% |- max[4] : array of the 4 most important maxima (structure)
% | |- x : intensity of the given maximum
% | |- y : amplitude of this maximum (number of MRI voxels with this value)
-% | |- power : difference of this maximum and the adjacent minima
+% | |- power : difference of this maximum and the adjacent minimum
% |- min[4] : array of the 3 most important minima (structure)
% | |- x : intensity of the given minimum
% | |- y : amplitude of this minimum (number of MRI voxels with this value)
-% | |- power : difference of this minimum and the adjacent maxima
+% | |- power : difference of this minimum and the adjacent maximum
% |- bgLevel : intensity value that separates the background and the objects (estimation)
% |- whiteLevel : white matter threshold
% |- intensityMax : maximum value in the volume
@@ -102,11 +102,11 @@
end
% Construct a regular Histogram function
-% Suppress all indices that has zero-values (to avoid previous normalizations)
-% NOTA : Do not consider the values at the intensity value 0, it may
+% Suppress all indices that have zero values (to avoid previous normalizations)
+% NOTA : Do not consider the first bin (intensity 0), it may
% not correspond to the real image Histogram...
index = find(Histogram.fncY > 10); % PREVIOUSLY: 100 instead of 10
-index = index(2:length(index));
+index = index(2:length(index)); % discard first bin
histoX = [0 Histogram.fncX(index)];
histoY = [0 Histogram.fncY(index)];
@@ -136,7 +136,7 @@
minIndex(diff(minIndex) == 1) = [];
end
-% Detect and deleting all "wrong" extrema (that are too close to each other)
+% Detect and delete all "wrong" extrema (that are too close to each other)
epsilon = max(histoX)*.02;
i = 1;
while(i <= length(maxIndex))
@@ -179,7 +179,7 @@
Histogram.max(i).y = histoY(maxIndex(i));
if(length(minIndex)>=1)
% If there is at least a minimum, power = distance between
- % maximum and adjacent minima
+ % maximum and adjacent minimum
Histogram.max(i).power = histoY(maxIndex(i)) - (histoY(minIndex(max(1, i-1))) + histoY(minIndex(min(length(minIndex), i))))./2;
else
% Else power = maximum value
@@ -223,26 +223,28 @@
defaultWhite = round(interp1(unikCumulFncY, unikFncX, .8));
Histogram.bgLevel = defaultBg;
Histogram.whiteLevel = defaultWhite;
- % Detect if the background has already been removed :
- % ie. if there is a unique 0 valued interval a the beginning of the Histogram
- % Practically : - nzero = length of the first 0-valued interval
- % - nnonzero = length of the first non-0-valued interval
- % - bg removed if : (nzero > 1) and (nnonzero > nzero)
- nzero = find(Histogram.fncY(2:length(Histogram.fncY)) ~= 0);
- nnonzero = find(Histogram.fncY((nzero(1)+1):length(Histogram.fncY)) == 0);
- if ((nzero(1)>2) && ~isempty(nnonzero) && (nnonzero(1) > nzero(1)))
- Histogram.bgLevel = nzero(1);
+ % Detect if the background has already been removed, i.e. if a range of low intensity values
+ % were replaced by zeros in the volume, producing an interval of empty bins at the beginning
+ % of the Histogram, after the first (zero-intensity) bin. Require also that the zero bin be
+ % the largest, as it seems denoising can produce this background removal effect but with a
+ % very low threshold which is not appropriate.
+ [~, iMaxBin] = max(Histogram.fncY);
+ % First non-empty bin after first bin.
+ iNonZero = find(Histogram.fncY(2:end) ~= 0, 1) + 1;
+ if iMaxBin == 1 && ~isempty(iNonZero) && iNonZero > 2
+ % There was an interval of empty bins. Set threshold at last empty bin intensity.
+ Histogram.bgLevel = Histogram.fncX(iNonZero - 1);
% Else, background has not been removed yet
- % If there is less than two maxima : use the default background threshold
+ % If there are fewer than two maxima : use the default background threshold
elseif (length(cat(1,Histogram.max.x)) < 2)
Histogram.bgLevel = defaultBg;
Histogram.whiteLevel = defaultWhite;
- % Else if there is more than one maxima :
+ % Else if there is more than one maximum :
else
- % If the highest maxima is > (3*second highest maxima) :
- % it is a background maxima : use the first minima after the
- % background maxima as background threshold
- % (and if this minima exist)
+ % If the highest maximum is > (3*second highest maximum) :
+ % it is a background maximum : use the first minimum after the
+ % background maximum as background threshold
+ % (and if this minimum exists)
[orderedMaxVal, orderedMaxInd] = sort(cat(1,Histogram.max.y), 'descend');
if ((orderedMaxVal(1) > 3*orderedMaxVal(2)) && (length(Histogram.min) >= orderedMaxInd(1)))
Histogram.bgLevel = Histogram.min(orderedMaxInd(1)).x;
@@ -251,7 +253,20 @@
Histogram.bgLevel = defaultBg;
end
end
-
+
+ case 'headgrad'
+ dX = mean(diff(Histogram.smoothFncX));
+ %Deriv = gradient(Histogram.smoothFncY, dX);
+ %SecondDeriv = gradient(Deriv, dX);
+ RemainderCumul = Histogram.smoothFncY ./ cumsum(Histogram.smoothFncY, 2, 'reverse');
+ DerivRC = gradient(RemainderCumul, dX);
+ %DerivRC2 = Deriv ./ cumsum(Histogram.smoothFncY, 2, 'reverse');
+ %figure; plot(Histogram.smoothFncX', [RemainderCumul', DerivRC']); legend({'hist/remaining', 'derivative'});
+ % Pick point where things flatten.
+ Histogram.bgLevel = Histogram.smoothFncX(find(DerivRC > -0.005 & DerivRC < DerivRC([2:end, end]), 1) + 2);
+ % Can't get white matter with gradient.
+ Histogram.whiteLevel = 0;
+
case 'brain'
% Determine an intensity value for the background/gray matter limit
% and the gray matter/white matter level
@@ -328,5 +343,4 @@
end
-
\ No newline at end of file
diff --git a/toolbox/anatomy/tess_check.m b/toolbox/anatomy/tess_check.m
new file mode 100644
index 0000000000..a3bc191a12
--- /dev/null
+++ b/toolbox/anatomy/tess_check.m
@@ -0,0 +1,201 @@
+function [isOk, Info] = tess_check(Vertices, Faces, isVerbose, isOpenOk, isShow)
+% TESS_CHECK: Check the integrity of a tesselation.
+%
+% USAGE: [isOk, Details] = tess_check(Vertices, Faces, Verbose)
+%
+% DESCRIPTION:
+% Check if a surface mesh is simple, closed, non self-intersecting, well oriented, duplicate
+% vertices or faces, etc. There are some custom checks, and if available, use the Matlab Lidar
+% toolbox. Could add meshcheckrepair from iso2mesh toolbox, and possibly others.
+%
+% INPUTS:
+% - Vertices : Mx3 double matrix
+% - Faces : Nx3 double matrix
+% - isOpenOk : An open surface is considered ok, otherwise flag non-closed as an issue.
+% - isVerbose : Write details to command window if any unexpected features.
+% - isShow : Display surface in new figure
+% OUTPUTS:
+% - isOk : All checks look good
+% - Details : Structure with all check statuses
+
+% @=============================================================================
+% This function is part of the Brainstorm software:
+% https://neuroimage.usc.edu/brainstorm
+%
+% Copyright (c) University of Southern California & McGill University
+% This software is distributed under the terms of the GNU General Public License
+% as published by the Free Software Foundation. Further details on the GPLv3
+% license can be found at http://www.gnu.org/copyleft/gpl.html.
+%
+% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
+% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
+% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
+% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
+% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
+%
+% For more information type "brainstorm license" at command prompt.
+% =============================================================================@
+%
+% Authors: Marc Lalancette, 2025
+
+% Parse inputs
+if nargin < 4 || isempty(isOpenOk)
+ isOpenOk = false;
+end
+if nargin < 3 || isempty(isVerbose)
+ isVerbose = true;
+end
+% Check matrices orientation
+if (size(Vertices, 2) ~= 3) || (size(Faces, 2) ~= 3)
+ error('Faces and Vertices must have 3 columns (X,Y,Z).');
+end
+
+isOk = true;
+Info = [];
+
+% Some custom checks first based on face edge connectivity.
+% First check duplicate faces (based on vertex indices, not coordinates)
+[~, FaceConn3] = tess_faceconn(Faces, 3); % 3 vertices in common = duplicate
+nAdjFaces = sum(FaceConn3, 2) - 1;
+Info.nDuplicate = sum(nAdjFaces > 0) / 2;
+% Then check for other issues
+[~, FaceConn] = tess_faceconn(Faces, 2);
+nAdjFaces = sum(FaceConn, 2) - 1;
+Info.nEdgeDisconnectedFaces = sum(nAdjFaces == 0);
+Info.nEdgeOverConnectedFaces = sum(nAdjFaces > 3); % 5, 7, 9 found after reducepatch
+% If there is an even number of adjacent faces > 3, there's something strange. An odd number could
+% be two lobes of the same surface just touching on an edge or face, not necessarily intersecting.
+Info.nEvenEdgeOverConnectedFaces = sum(nAdjFaces > 3 & ~mod(nAdjFaces,2));
+Info.nBoundaryFaces = sum(nAdjFaces == 1) + sum(nAdjFaces == 2);
+
+if Info.nDuplicate > 0
+ isOk = false;
+ if isVerbose
+ fprintf('BST>Surface has %d duplicate faces.\n', Info.nDuplicate);
+ end
+end
+if Info.nEdgeDisconnectedFaces > 0 % yes many after reducepatch
+ isOk = false;
+ if isVerbose
+ fprintf('BST>Surface has %d disconnected faces (no shared edges, but possibly one shared vertex).\n', Info.nEdgeDisconnectedFaces);
+ end
+end
+if Info.nEdgeOverConnectedFaces > 0
+ isOk = false;
+ if isVerbose
+ fprintf('BST>Surface intersects or has touching/duplicate edges or faces (%d).\n', Info.nEdgeOverConnectedFaces);
+ if Info.nEvenEdgeOverConnectedFaces > 0
+ fprintf('BST>Surface has edges shared by 3,5,7,... faces, indicating strange topology (%d).\n', Info.nEvenEdgeOverConnectedFaces);
+ end
+ end
+end
+% Note openness if ok so far
+if isOk
+ if Info.nBoundaryFaces > 0
+ Info.isOpen = true;
+ else
+ Info.isOpen = false;
+ end
+else % don't bother defining if surface is weird
+ Info.isOpen = [];
+end
+% Check for boundary faces either way if we want closed
+if Info.nBoundaryFaces > 0 && ~isOpenOk
+ isOk = false;
+ if isVerbose
+ fprintf('BST>Surface has %d boundary edges.\n', Info.nBoundaryFaces);
+ end
+end
+
+% Check orientation if ok so far; so we should have each edge shared by 2 faces, or 1 if open.
+Info.isOriented = [];
+if isOk
+ Edges = [Faces(:), [Faces(:, 2); Faces(:, 3); Faces(:, 1)]];
+ isEdgeFlip = Edges(:, 1) > Edges(:, 2);
+ [~, ~, iE] = unique(sort(Edges, 2), 'rows');
+ % Look for boundaries of open surface.
+ isBoundE = false(size(Edges, 1), 1);
+ [iE, iSort] = sort(iE);
+ % Add one more row for the loop to also evaluate the last real edge.
+ iE(end+1,:) = 0;
+ n = 1;
+ for i = 2:numel(iE)
+ if iE(i) ~= iE(i-1)
+ % Evaluate previous edge
+ if n == 1
+ % Only one copy, boundary edge.
+ isBoundE(iE(i-1)) = true;
+ elseif n == 2
+ % Two faces, were the orientations different?
+ if sum(isEdgeFlip(iSort([i-1,i-2]))) ~= 1 % should be sum([0,1])
+ isOk = false;
+ if isVerbose
+ fprintf('BST>Surface not well oriented (face normals are mixed pointing in and out).\n');
+ end
+ break;
+ end
+ % Reset for new edge
+ n = 1;
+ else
+ % Previously undetected issue with edge shared among more than 2 faces.
+ isOk = false;
+ if isVerbose
+ fprintf('BST>Surface has edge shared among more than 2 faces.\n');
+ end
+ break;
+ end
+ else
+ n = n + 1;
+ end
+ end
+end
+
+if isShow
+ FigureId = db_template('FigureId');
+ FigureId.Type = '3DViz';
+ hFig = figure_3d('CreateFigure', FigureId);
+ figure_3d('PlotSurface', hFig, Faces, Vertices, [1,1,1], 0); % color, transparency required
+ figure_3d('ViewAxis', hFig, true); % isVisible
+ hFig.Visible = "on";
+end
+
+
+% -------------------------------------------
+% Check if Lidar Toolbox is installed (requires image processing + computer vision)
+isLidarToolbox = exist('surfaceMesh', 'file') == 2;
+if ~isLidarToolbox
+ fprintf('BST>tess_downsize method "simplify" requires Matlab''s Lidar Toolbox, which was not found.\n');
+ return;
+end
+% Create mesh object
+oMesh = surfaceMesh(Vertices, Faces);
+
+% Check all mesh features Lidar Toolbox offers
+Info.isVertexManifold = isVertexManifold(oMesh); % Check if surface mesh is vertex-manifold
+Info.isEdgeManifold = isEdgeManifold(oMesh, true); % allow boundary edges
+Info.isClosedManifold = isEdgeManifold(oMesh, false); % allow boundary edges
+Info.isOrientable = isOrientable(oMesh); % Check if surface mesh is orientable
+% Self-intersecting test is slow (not sure how long, didn't wait more than 15 minutes, uses one core only)
+%Info.isSelfIntersecting = isSelfIntersecting(oMesh); % Check if surface mesh is self-intersecting
+Info.isWatertight = isWatertight(oMesh); % Check if surface mesh is watertight
+% removeDefects could be used in tess_clean
+
+if isVerbose
+ if (isOpenOk && ~Info.isEdgeManifold) || (~isOpenOk && ~Info.isClosedManifold)
+ fprintf('BST>Surface not "edge manifold" (each edge has at most one face on each side).\n');
+ end
+ if ~Info.isVertexManifold
+ fprintf('BST>Surface not "vertex manifold" (like a "fan" at each vertex).\n');
+ end
+ if ~Info.isOrientable
+ fprintf('BST>Surface not well oriented (face normals are mixed pointing in and out).\n');
+ end
+ % if Info.isSelfIntersecting
+ % fprintf('BST>Surface self-intersects.\n');
+ % end
+end
+
+end
+
+
+
diff --git a/toolbox/anatomy/tess_downsize.m b/toolbox/anatomy/tess_downsize.m
index a97f30bca2..1a62bcd880 100644
--- a/toolbox/anatomy/tess_downsize.m
+++ b/toolbox/anatomy/tess_downsize.m
@@ -2,9 +2,12 @@
% TESS_DOWNSIZE: Reduces the number of vertices in a surface file.
%
% USAGE: [NewTessFile, iSurface, I, J] = tess_downsize(TessFile, newNbVertices=[ask], Method=[ask]);
+% [NewTessMat, iSurface, I, J] = tess_downsize(TessMat, newNbVertices=[ask], Method=[ask]);
%
% INPUT:
% - TessFile : Full path to surface file to decimate
+% - TessMat : Already loaded surface structure, a structure is returned in that case and no
+% file is saved.
% - newNbVertices : Desired number of vertices
% - Method : {'reducepatch', 'reducepatch_subdiv', 'iso2mesh', 'iso2mesh_project'}
% OUTPUT:
@@ -56,11 +59,17 @@
I = [];
J = [];
+% Surface structure now accepted as input
+isTessInput = isstruct(TessFile);
%% ===== ASK FOR MISSING OPTIONS =====
% Get the number of vertices
-VarInfo = whos('-file',file_fullpath(TessFile),'Vertices');
-oldNbVertices = VarInfo.size(1);
+if isTessInput
+ oldNbVertices = size(TessFile.Vertices, 1);
+else
+ VarInfo = whos('-file',file_fullpath(TessFile),'Vertices');
+ oldNbVertices = VarInfo.size(1);
+end
% If new number of vertices was not provided: ask user
if isempty(newNbVertices)
% Ask user the new number of vertices
@@ -82,6 +91,9 @@
return;
end
+% Check if Lidar Toolbox is installed (requires image processing + computer vision)
+isLidarToolbox = exist('surfaceMesh', 'file') == 2;
+
% Ask for resampling method
if isempty(Method)
% Downsize methods strings
@@ -99,6 +111,12 @@
% ' | - Homogeneous mesh but possible topological problems
' ...
% ' | - Damages the atlases and the subject co-registration']},
+ if isLidarToolbox
+ methods_str{end+1} = ...
+ ['Matlab''s simplify, from Lidar Toolbox:
' ...
+ ' | - Possibly better at preserving shape topology than reducepatch
' ...
+ ' | - Deletes the atlases and the subjects co-registration'];
+ end
% Identify textured surfaces (color info is present) and show available methods for them
VarInfo = whos('-file',file_fullpath(TessFile), 'Color');
if ~isempty(VarInfo) && all(VarInfo.size ~= 0)
@@ -114,7 +132,8 @@
case 1, Method = 'reducepatch';
case 2, Method = 'reducepatch_subdiv';
case 3, Method = 'iso2mesh';
- case 4, Method = 'iso2mesh_project';
+ case 4, Method = 'simplify';
+ % case 4, Method = 'iso2mesh_project';
end
end
@@ -128,17 +147,27 @@
end
%% ===== LOAD FILE =====
-% Progress bar
-bst_progress('start', 'Resample surface', 'Loading file...');
-% Load file
-TessMat = in_tess_bst(TessFile);
-% Prepare variables
+if isTessInput
+ TessMat = TessFile;
+ if ~isfield(TessMat, 'Comment')
+ TessMat.Comment = 'iso head';
+ end
+else
+ % Progress bar
+ bst_progress('start', 'Resample surface', 'Loading file...');
+ % Load file
+ TessMat = in_tess_bst(TessFile);
+ % Prepare variables
+end
TessMat.Faces = double(TessMat.Faces);
TessMat.Vertices = double(TessMat.Vertices);
-TessMat.Color = double(TessMat.Color);
+if isfield(TessMat, 'Color')
+ TessMat.Color = double(TessMat.Color);
+else
+ TessMat.Color = [];
+end
dsFactor = newNbVertices / size(TessMat.Vertices, 1);
-
%% ===== RESAMPLE =====
bst_progress('start', 'Resample surface', ['Resampling surface: ' TessMat.Comment '...']);
% Resampling methods
@@ -422,8 +451,29 @@
NewTessMat.Faces = Faces;
NewTessMat.Vertices = Vertices;
MethodTag = '_iso2mesh_proj';
+
+ % ===== SIMPLIFY =====
+ % Matlab's Lidar Toolbox simplify (since R2022b)
+ case 'simplify'
+ if ~isLidarToolbox
+ fprintf('BST>tess_downsize method "simplify" requires Matlab''s Lidar Toolbox, which was not found.\n');
+ end
+ % Create mesh object
+ oMesh = surfaceMesh(TessMat.Vertices, TessMat.Faces);
+
+ % Reduce number of vertices
+ simplify(oMesh, 'TargetNumFaces', (newNbVertices - 2) * 2); % no output variable for this toolbox!
+ NewTessMat.Faces = oMesh.Faces;
+ NewTessMat.Vertices = oMesh.Vertices;
+
end
+if isTessInput
+ % This should go after "remove folded faces" if that step was desired... skipping for now for
+ % tess_isohead as it does its own checks and corrections.
+ NewTessFile = NewTessMat;
+ return;
+end
%% ===== REMOVE FOLDED FACES =====
% Find equal faces
@@ -540,15 +590,18 @@
%% ===== UPDATE DATABASE =====
-% Save downsized surface file
-bst_save(NewTessFile, NewTessMat, 'v7');
-% Make output filename relative
-NewTessFile = file_short(NewTessFile);
-% Get subject
-[sSubject, iSubject] = bst_get('SurfaceFile', TessFile);
-% Register this file in Brainstorm database
-iSurface = db_add_surface(iSubject, NewTessFile, NewComment);
-
+% Save downsized surface file, or return structure
+if isTessInput
+ NewTessFile = NewTessMat;
+else
+ bst_save(NewTessFile, NewTessMat, 'v7');
+ % Make output filename relative
+ NewTessFile = file_short(NewTessFile);
+ % Get subject
+ [~, iSubject] = bst_get('SurfaceFile', TessFile);
+ % Register this file in Brainstorm database
+ iSurface = db_add_surface(iSubject, NewTessFile, NewComment);
+end
% Close progress bar
bst_progress('stop');
diff --git a/toolbox/anatomy/tess_faceconn.m b/toolbox/anatomy/tess_faceconn.m
index 98b999ad5a..d22f6d26b9 100644
--- a/toolbox/anatomy/tess_faceconn.m
+++ b/toolbox/anatomy/tess_faceconn.m
@@ -1,12 +1,15 @@
-function [VertFacesConn, FaceConn] = tess_faceconn(Faces)
+function [VertFacesConn, FaceConn] = tess_faceconn(Faces, nVert)
% TESS_FACECONN: Computes faces connectivity.
%
-% USAGE: [VertFacesConn, FaceConn] = tess_faceconn(Faces);
+% USAGE: [VertFacesConn, FaceConn] = tess_faceconn(Faces, nVert=1);
%
% INPUT:
% - Faces : Nx3 double matrix
+% - nVert : Number of vertices that a pair of faces need to share to be considered connected.
+% nVert=1 by default, but for finding only edge-adjacent faces, use 2.
% OUTPUT:
-% - FacesConn : sparse matrix [nVertices x nFaces]
+% - VertFacesConn : sparse matrix [nVertices x nFaces]
+% - FacesConn : sparse matrix [nFaces x nFaces]
% @=============================================================================
% This function is part of the Brainstorm software:
@@ -28,6 +31,11 @@
%
% Authors: Anand Joshi, Dimitrios Pantazis, November 2007
% Francois Tadel, 2008-2010
+% Marc Lalancette, 2025
+
+if nargin < 2 || isempty(nVert)
+ nVert = 1;
+end
% Check matrices orientation
if (size(Faces, 2) ~= 3)
@@ -43,7 +51,7 @@
% Build FacesConn
if (nargout > 1)
- FaceConn = (VertFacesConn' * VertFacesConn) > 0;
+ FaceConn = (VertFacesConn' * VertFacesConn) >= nVert;
end
diff --git a/toolbox/anatomy/tess_isohead.m b/toolbox/anatomy/tess_isohead.m
index dc81e29e1d..c0f4944347 100644
--- a/toolbox/anatomy/tess_isohead.m
+++ b/toolbox/anatomy/tess_isohead.m
@@ -1,220 +1,710 @@
-function [HeadFile, iSurface] = tess_isohead(iSubject, nVertices, erodeFactor, fillFactor, bgLevel, Comment)
-% TESS_GENERATE: Reconstruct a head surface based on the MRI, based on an isosurface
-%
-% USAGE: [HeadFile, iSurface] = tess_isohead(iSubject, nVertices=10000, erodeFactor=0, fillFactor=2, bgLevel=GuessFromHistorgram, Comment)
-% [HeadFile, iSurface] = tess_isohead(MriFile, nVertices=10000, erodeFactor=0, fillFactor=2, bgLevel=GuessFromHistorgram, Comment)
-% [Vertices, Faces] = tess_isohead(sMri, nVertices=10000, erodeFactor=0, fillFactor=2)
-%
-% If input is loaded MRI structure, no surface file is created and the surface vertices and faces are returned instead.
-
-% @=============================================================================
-% This function is part of the Brainstorm software:
-% https://neuroimage.usc.edu/brainstorm
-%
-% Copyright (c) University of Southern California & McGill University
-% This software is distributed under the terms of the GNU General Public License
-% as published by the Free Software Foundation. Further details on the GPLv3
-% license can be found at http://www.gnu.org/copyleft/gpl.html.
-%
-% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
-% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
-% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
-% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
-% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
-%
-% For more information type "brainstorm license" at command prompt.
-% =============================================================================@
-%
-% Authors: Francois Tadel, 2012-2022
-
-%% ===== PARSE INPUTS =====
-% Initialize returned variables
-HeadFile = [];
-iSurface = [];
-isSave = true;
-% Parse inputs
-if (nargin < 6)
- if nargin == 5
- % Handle legacy call: tess_isohead(iSubject, nVertices, erodeFactor, fillFactor, Comment)
- if ischar(bgLevel)
- Comment = bgLevel;
- bgLevel = [];
- % Parameter 'bgLevel' is provided: tess_isohead(iSubject, nVertices, erodeFactor, fillFactor, bgLevel)
+function [HeadFile, iSurface] = tess_isohead(iSubject, nVertices, erodeFactor, fillFactor, bgLevel, Comment, isGradient, Method)
+ % TESS_GENERATE: Reconstruct a head surface based on the MRI, based on an isosurface
+ %
+ % USAGE: [HeadFile, iSurface] = tess_isohead(iSubject, nVertices=10000, erodeFactor=0, fillFactor=2, bgLevel=GuessFromHistorgram, Comment)
+ % [HeadFile, iSurface] = tess_isohead(MriFile, nVertices=10000, erodeFactor=0, fillFactor=2, bgLevel=GuessFromHistorgram, Comment)
+ % [Vertices, Faces] = tess_isohead(sMri, nVertices=10000, erodeFactor=0, fillFactor=2)
+ %
+ % If input is loaded MRI structure, no surface file is created and the surface vertices and faces are returned instead.
+
+ % @=============================================================================
+ % This function is part of the Brainstorm software:
+ % https://neuroimage.usc.edu/brainstorm
+ %
+ % Copyright (c) University of Southern California & McGill University
+ % This software is distributed under the terms of the GNU General Public License
+ % as published by the Free Software Foundation. Further details on the GPLv3
+ % license can be found at http://www.gnu.org/copyleft/gpl.html.
+ %
+ % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
+ % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
+ % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
+ % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
+ % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
+ %
+ % For more information type "brainstorm license" at command prompt.
+ % =============================================================================@
+ %
+ % Authors: Francois Tadel, 2012-2022
+ % Marc Lalancette, 2022-2025
+
+ % To visualize steps for debugging.
+ isDebugVis = false;
+ nDebugVisSlices = 9; %#ok
+
+ %% ===== PARSE INPUTS =====
+ % Initialize returned variables
+ HeadFile = [];
+ iSurface = [];
+ isSave = true;
+ % Parse inputs
+ if (nargin < 8) || isempty(Method)
+ Method = 'simplify';
+ end
+ if strcmpi(Method, 'simplify')
+ % Check if Lidar Toolbox is installed (requires image processing + computer vision)
+ isLidarToolbox = exist('surfaceMesh', 'file') == 2;
+ if ~isLidarToolbox
+ bst_error('Lidar toolbox required for method ''simplify''.');
+ end
+ end
+ if (nargin < 7) || isempty(isGradient)
+ isGradient = false;
+ end
+ if (nargin < 6)
+ if nargin == 5
+ % Handle legacy call: tess_isohead(iSubject, nVertices, erodeFactor, fillFactor, Comment)
+ if ischar(bgLevel)
+ Comment = bgLevel;
+ bgLevel = [];
+ % Parameter 'bgLevel' is provided: tess_isohead(iSubject, nVertices, erodeFactor, fillFactor, bgLevel)
+ else
+ Comment = [];
+ end
+ % Call tess_isohead(iSubject, nVertices, erodeFactor, fillFactor)
else
+ bgLevel = [];
Comment = [];
end
- % Call tess_isohead(iSubject, nVertices, erodeFactor, fillFactor)
+ end
+ % MriFile instead of subject index
+ sMri = [];
+ if ischar(iSubject)
+ MriFile = file_short(iSubject);
+ [sSubject, iSubject] = bst_get('MriFile', MriFile);
+ elseif isnumeric(iSubject)
+ % Get subject
+ sSubject = bst_get('Subject', iSubject);
+ MriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
+ elseif isstruct(iSubject)
+ sMri = iSubject;
+ MriFile = sMri.FileName;
+ [sSubject, iSubject] = bst_get('MriFile', MriFile);
+ % Don't save a surface file, instead return surface directly.
+ isSave = false;
else
- bgLevel = [];
- Comment = [];
+ error('Wrong input type.');
end
-end
-% MriFile instead of subject index
-sMri = [];
-if ischar(iSubject)
- MriFile = file_short(iSubject);
- [sSubject, iSubject] = bst_get('MriFile', MriFile);
-elseif isnumeric(iSubject)
- % Get subject
- sSubject = bst_get('Subject', iSubject);
- MriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
-elseif isstruct(iSubject)
- sMri = iSubject;
- MriFile = sMri.FileName;
- [sSubject, iSubject] = bst_get('MriFile', MriFile);
- % Don't save a surface file, instead return surface directly.
- isSave = false;
-else
- error('Wrong input type.');
-end
-%% ===== LOAD MRI =====
-isProgress = ~bst_progress('isVisible');
-if isempty(sMri)
- % Load MRI
- bst_progress('start', 'Generate head surface', 'Loading MRI...');
- sMri = bst_memory('LoadMri', MriFile);
+ %% ===== LOAD MRI =====
+ isProgress = ~bst_progress('isVisible');
+ if isempty(sMri)
+ % Load MRI
+ bst_progress('start', 'Generate head surface', 'Loading MRI...');
+ sMri = bst_memory('LoadMri', MriFile);
+ if isProgress
+ bst_progress('stop');
+ end
+ end
+ % Save current scouts modifications
+ panel_scout('SaveModifications');
+ % If subject is using the default anatomy: use the default subject instead
+ if sSubject.UseDefaultAnat
+ iSubject = 0;
+ end
+ % Check layers
+ if isempty(sSubject.iAnatomy) || isempty(sSubject.Anatomy)
+ bst_error('The generate of the head surface requires at least the MRI of the subject.', 'Head surface', 0);
+ return
+ end
+ % Check that everything is there
+ if ~isfield(sMri, 'Histogram') || isempty(sMri.Histogram) || isempty(sMri.SCS) || isempty(sMri.SCS.NAS) || isempty(sMri.SCS.LPA) || isempty(sMri.SCS.RPA)
+ bst_error('You need to set the fiducial points in the MRI first.', 'Head surface', 0);
+ return
+ end
+ % Guess background level
+ if isempty(bgLevel) && ~isGradient
+ bgLevel = sMri.Histogram.bgLevel;
+ end
+
+ %% ===== ASK PARAMETERS =====
+ % Ask user to set the parameters if they are not set
+ if (nargin < 4) || isempty(erodeFactor) || isempty(nVertices)
+ res = java_dialog('input', {'Number of vertices [integer]:', 'Erode factor [0,1,2,3]:', 'Fill holes factor [0,1,2,3]:', 'Background threshold:
(guessed from MRI histogram)'}, 'Generate head surface', [], {'15000', '0', '0', num2str(bgLevel)});
+ % If user cancelled: return
+ if isempty(res)
+ return
+ end
+ % Get new values
+ nVertices = str2double(res{1});
+ erodeFactor = str2double(res{2});
+ fillFactor = str2double(res{3});
+ bgLevel = str2double(res{4});
+ if isempty(bgLevel)
+ bgLevel = sMri.Histogram.bgLevel;
+ end
+ end
+ % Check parameters values
+ if isempty(nVertices) || (nVertices < 50) || (nVertices ~= round(nVertices)) || isempty(erodeFactor) || ~ismember(erodeFactor,[0,1,2,3]) || isempty(fillFactor) || ~ismember(fillFactor,[0,1,2,3])
+ bst_error('Invalid parameters.', 'Head surface', 0);
+ return
+ end
+
+
+ %% ===== CREATE HEAD MASK =====
+ % Progress bar
+ bst_progress('start', 'Generate head surface', 'Creating head mask...', 0, 100);
+ % Threshold mri to the level estimated in the histogram
+ if isGradient
+ isGradLocalMax = false;
+ % Compute gradient
+ % Find appropriate threshold from gradient histogram
+ % TODO: need to find a robust way to do this. Only verified on one
+ % relatively bad MRI sequence with preprocessing (debias, denoise).
+ [Grad, VectGrad] = NormGradient(sMri.Cube(:,:,:,1)); %#ok
+ if isempty(bgLevel)
+ Hist = mri_histogram(Grad, [], 'headgrad');
+ bgLevel = Hist.bgLevel;
+ end
+ if isGradLocalMax
+ % Index gymnastics... Is there a simpler way to do this (other than looping)?
+ nVol = [1, cumprod(size(Grad))]';
+ [unused, UpDir] = max(abs(reshape(VectGrad, nVol(4), [])), [], 2); % (nVol, 1)
+ UpDirSign = sign(VectGrad((1:nVol(4))' + (UpDir-1) * nVol(4)));
+ % Get neighboring value of the gradient in the increasing gradiant direction.
+ % Using linear indices shaped as 3d array, which will give back a 3d array.
+ iUpGrad = zeros(size(Grad));
+ iUpGrad(:) = UpDirSign .* nVol(UpDir); % change in index: +-1 along appropriate dimension for each voxel, in linear indices
+ % Removing problematic indices at edges.
+ iUpGrad([1, end], :, :) = 0;
+ iUpGrad(:, [1, end], :) = 0;
+ iUpGrad(:, :, [1, end]) = 0;
+ iUpGrad(:) = iUpGrad(:) + (1:nVol(4))'; % adding change to each element index
+ UpGrad = Grad(iUpGrad);
+ headmask = Grad > bgLevel & Grad >= UpGrad;
+ else
+ headmask = Grad > bgLevel;
+ end
+ else
+ headmask = sMri.Cube(:,:,:,1) > bgLevel;
+ end
+ % Closing all the faces of the cube
+ headmask(1,:,:) = 0;
+ headmask(end,:,:) = 0;
+ headmask(:,1,:) = 0;
+ headmask(:,end,:) = 0;
+ headmask(:,:,1) = 0;
+ headmask(:,:,end) = 0;
+ if isDebugVis
+ view_mri_slices(headmask, 'x', nDebugVisSlices); %#ok<*UNRCH>
+ end
+
+ % Fill neck holes (bones, etc.) where it is cut at edge of volume.
+ bst_progress('text', 'Filling holes and removing disconnected parts...');
+ % Brainstorm reorients MRI so voxels are in RAS. But do all faces in case the bounding box was too
+ % small and another part is cut (e.g. nose).
+
+ % Number of slices to average to smooth out noise in low SNR regions (e.g. around neck and chin).
+ % 4,3 worked ok in noisy scan, but probably best to denoise entire scan first.
+ nSlices = 1; % 1 = no averaging.
+ FillThresh = 1; %min(nSlices, floor(nSlices/2)+1);
+ if FillThresh > nSlices || FillThresh < floor(nSlices/2)
+ error('Bad hard-coded FillThresh.');
+ end
+ for iDim = 1:3
+ % Swap slice dimension into first position. For a single swap, the permutation is it's own inverse.
+ switch iDim
+ case 1
+ Perm = 1:3;
+ case 2
+ Perm = [2, 1, 3];
+ case 3
+ Perm = [3, 2, 1];
+ end
+ TempMask = permute(headmask, Perm);
+ % Edit second and second-to-last slices. Flip the array to reuse code with same indices.
+ for isFlip = [false, true]
+ if isFlip
+ TempMask = flip(TempMask, 1);
+ end
+ % Skip if just background (e.g. above or behind head)
+ if ~any(any(squeeze(TempMask(2, :, :))))
+ % Flip back and move on
+ if isFlip
+ TempMask = flip(TempMask, 1);
+ end
+ continue;
+ end
+ if isDebugVis
+ figure; imagesc(squeeze(TempMask(2,:,:))); colormap('gray'); axis equal; title(sprintf('Dim %d, Flip %d', iDim, isFlip));
+ end
+ Slice = sum(TempMask(2:2+nSlices-1, :, :), 1);
+ if isDebugVis && nSlices > 1
+ figure; imagesc(squeeze(Slice)); colormap('gray'); axis equal; title(sprintf('Dim %d, Flip %d, Avg %d slices', iDim, isFlip, nSlices));
+ end
+ Slice = Slice >= FillThresh;
+ % Skip if just background (previous check had just some noise)
+ if ~any(any(squeeze(Slice)))
+ if isFlip
+ TempMask = flip(TempMask, 1);
+ end
+ continue;
+ end
+ Slice = FillConcaveVolume(Slice, true); % isClean
+ % Slice = Slice | (Fill(Slice, 2) & Fill(Slice, 3));
+ if isDebugVis
+ figure; imagesc(squeeze(Slice)); colormap('gray'); axis equal; title(sprintf('Dim %d, Flip %d, Filled', iDim, isFlip));
+ end
+ [Slice, isFail] = CenterSpread(Slice);
+ % Avoid warnings for slices other than neck.
+ if isFail && iDim == 3 && isFlip == false % inferior slice
+ warning('CenterSpread failed for filling "neck" slice. Resulting head surface may be problematic.');
+ % Keep original.
+ else
+ % Keep filled in slice.
+ TempMask(2, :, :) = Slice;
+ end
+ if isDebugVis
+ figure; imagesc(squeeze(Slice)); colormap('gray'); axis equal; title(sprintf('Dim %d, Flip %d, Center spread', iDim, isFlip));
+ end
+ % Flip back this dimension
+ if isFlip
+ TempMask = flip(TempMask, 1);
+ end
+ end
+ % Permute back dimensions to original order.
+ headmask = permute(TempMask, Perm);
+ end
+ % Fill holes
+ headmask = FillConcaveVolume(headmask, true); % clean, which may remove a few more original 1-voxel-wide or noise bits
+ if isDebugVis
+ view_mri_slices(headmask, 'x', nDebugVisSlices); title('Filled');
+ end
+ % Keep only central connected volume (trim "beard" or bubbles)
+ headmask = CenterSpread(headmask);
+ bst_progress('inc', 15);
+
+ if isDebugVis
+ view_mri_slices(headmask, 'x', nDebugVisSlices); title('Center spread');
+ end
+
+
+ %% ===== CREATE SURFACE =====
+ % Compute isosurface
+ bst_progress('text', 'Creating isosurface...');
+
+ switch Method
+ case 'iso2mesh'
+ method = 'cgalsurf';
+ opt.radbound = 4; % max radius of the Delaunay sphere - adjust to get desired vertex numbers
+ opt.distbound = 1; % max distance from isosurface
+ dofix = 1; % don't know if needed
+
+ [sHead.Vertices, sHead.Faces, regions, holes] = vol2surf(headmask, ...
+ 1:size(headmask,1), 1:size(headmask,2), 1:size(headmask,3), opt, dofix, method); % ,isovalues
+ if size(regions, 1) > 1
+ bst_error('Multiple regions returned.\n');
+ return;
+ elseif ~isempty(holes)
+ bst_error('Holes present.\n');
+ return;
+ end
+ % Remove region label
+ sHead.Faces(:,4) = [];
+ if isDebugVis
+ fprintf('iso2mesh surface\n');
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); % verbose, not open, show
+ title('iso2mesh surface', 'color', 'white');
+ end
+ bst_progress('inc', 45);
+
+ case {'reducepatch', 'simplify'}
+ % Could have avoided x-y flip by specifying XYZ in isosurface...
+ [sHead.Faces, sHead.Vertices] = mri_isosurface(headmask, 0.5);
+ % Flip x-y back to our voxel coordinates
+ sHead.Vertices = sHead.Vertices(:, [2, 1, 3]);
+ % Flip to have desired face orientations (seems inconsistent if needed or not).
+ % sHead.Faces = sHead.Faces(:, [2, 1, 3]);
+ if isDebugVis
+ fprintf('mri_isohead surface\n');
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); % verbose, not open, show
+ title('isosurface', 'color', 'white');
+ end
+ bst_progress('inc', 10);
+
+ % Remove small objects
+ bst_progress('text', 'Removing small patches...');
+ nVertTemp = size(sHead.Vertices, 1);
+ [sHead.Vertices, sHead.Faces] = tess_remove_small(sHead.Vertices, sHead.Faces);
+ if isDebugVis && nVertTemp > size(sHead.Vertices, 1) % only if something was removed in previous step
+ fprintf('BST>Some disconnected small patches removed (%d vertices).\n', nVertTemp - size(sHead.Vertices, 1));
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); %#ok % verbose, not open, show
+ title('isosurface & small removed', 'color', 'white');
+ end
+ bst_progress('inc', 15);
+
+ % TODO: No existing functions, including from iso2mesh, correctly fix topology issues, which
+ % are present after downsampling with all methods tested (including again iso2mesh).
+ % tess_clean is very strange, it doesn't look at face locations, only the normals. And after
+ % isosurface, many faces are parallel.
+
+ % Smooth voxel artefacts, but preserve shape and volume.
+ bst_progress('text', 'Smoothing voxel artefacts...');
+ % Should normally use 1 as voxel size, but using a larger value smooths.
+ % Restrict iterations to make it faster, smooth a bit more (normal to surface
+ % only) after downsampling.
+ sHead.Vertices = SurfaceSmooth(sHead.Vertices, sHead.Faces, 2, [], 5, [], false); % voxel/smoothing size, iterations, verbose
+ if isDebugVis
+ fprintf('mildly smoothed isosurface\n');
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); % verbose, not open, show
+ title('mildly smoother isosurface', 'color', 'white');
+ end
+ bst_progress('inc', 20);
+ end
+
+ % Downsampling surface
+ if (length(sHead.Vertices) > 1.5* nVertices)
+ bst_progress('text', 'Downsampling surface...');
+ % Modified tess_downsize to accept sHead
+ sHead = tess_downsize(sHead, nVertices, Method);
+ if isDebugVis
+ fprintf('reduced surface (%s)\n', Method);
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); % verbose, not open, show
+ title('reduced', 'color', 'white')
+ end
+ % Fix this patch
+ if ~strcmpi(Method, 'iso2mesh') % I don't think iso2mesh returns multiple disconnected regions.
+ nVertTemp = size(sHead.Vertices, 1);
+ [sHead.Vertices, sHead.Faces] = tess_remove_small(sHead.Vertices, sHead.Faces);
+ if isDebugVis && nVertTemp > size(sHead.Vertices, 1) % only if something was removed in previous step
+ fprintf('BST>Some disconnected small patches removed (%d vertices).\n', nVertTemp - size(sHead.Vertices, 1));
+ fprintf('reduced surface (small disconnected parts removed)\n');
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); %#ok % verbose, not open, show
+ title('reduced & small removed', 'color', 'white')
+ end
+ end
+ end
+ bst_progress('inc', 15);
+
+ bst_progress('text', 'Smoothing...');
+ sHead.Vertices = SurfaceSmooth(sHead.Vertices, sHead.Faces, 1.5, [], 45, 0, false); % voxel/smoothing size, iterations, freedom (normal), verbose
+ if isDebugVis
+ fprintf('final smoothed surface\n');
+ [isOk, Info] = tess_check(sHead.Vertices, sHead.Faces, true, false, true); % verbose, not open, show
+ end
+ bst_progress('inc', 10);
+
+ % Convert to SCS
+ sHead.Vertices = cs_convert(sMri, 'voxel', 'scs', sHead.Vertices);
+ % Flip face order to Brainstorm convention
+ sHead.Faces = sHead.Faces(:,[2,1,3]);
+
+
+ %% ===== SAVE FILES =====
+ if isSave
+ bst_progress('text', 'Saving new file...');
+ % Create output filenames
+ ProtocolInfo = bst_get('ProtocolInfo');
+ SurfaceDir = bst_fullfile(ProtocolInfo.SUBJECTS, bst_fileparts(MriFile));
+ HeadFile = file_unique(bst_fullfile(SurfaceDir, 'tess_head_mask.mat'));
+ % Save head
+ if ~isempty(Comment)
+ sHead.Comment = Comment;
+ else
+ sHead.Comment = sprintf('head mask (%d,%d,%d,%d)', nVertices, erodeFactor, fillFactor, round(bgLevel));
+ end
+ sHead = bst_history('add', sHead, 'bem', 'Head surface generated with Brainstorm');
+ bst_save(HeadFile, sHead, 'v7');
+ iSurface = db_add_surface( iSubject, HeadFile, sHead.Comment);
+ else
+ % Return surface
+ HeadFile = sHead.Vertices;
+ iSurface = sHead.Faces;
+ end
+
+ % Close, success
if isProgress
bst_progress('stop');
end
end
-% Save current scouts modifications
-panel_scout('SaveModifications');
-% If subject is using the default anatomy: use the default subject instead
-if sSubject.UseDefaultAnat
- iSubject = 0;
-end
-% Check layers
-if isempty(sSubject.iAnatomy) || isempty(sSubject.Anatomy)
- bst_error('The generate of the head surface requires at least the MRI of the subject.', 'Head surface', 0);
- return
-end
-% Check that everything is there
-if ~isfield(sMri, 'Histogram') || isempty(sMri.Histogram) || isempty(sMri.SCS) || isempty(sMri.SCS.NAS) || isempty(sMri.SCS.LPA) || isempty(sMri.SCS.RPA)
- bst_error('You need to set the fiducial points in the MRI first.', 'Head surface', 0);
- return
-end
-% Guess background level
-if isempty(bgLevel)
- bgLevel = sMri.Histogram.bgLevel;
-end
-%% ===== ASK PARAMETERS =====
-% Ask user to set the parameters if they are not set
-if (nargin < 4) || isempty(erodeFactor) || isempty(nVertices)
- res = java_dialog('input', {'Number of vertices [integer]:', 'Erode factor [0,1,2,3]:', 'Fill holes factor [0,1,2,3]:', 'Background threshold:
(guessed from MRI histogram)'}, 'Generate head surface', [], {'10000', '0', '2', num2str(bgLevel)});
- % If user cancelled: return
- if isempty(res)
- return
+%% ===== Subfunctions =====
+function mask = FillConcaveVolume(mask, isClean)
+ % Try to fill the interior of a concave volume. For a head, expects the neck "cut" to be
+ % previously filled. This method still depends on the orientation of the object vs the volume
+ % dimensions (x,y,z axes). But for a head shape, not too important. Mostly noticeable behind
+ % ears or in nostrils for example, depending on their angles.
+
+ if nargin < 2 || isempty(isClean)
+ % Default to not remove any of the original mask. But called with true in this file.
+ isClean = false;
end
- % Get new values
- nVertices = str2num(res{1});
- erodeFactor = str2num(res{2});
- fillFactor = str2num(res{3});
- bgLevel = str2num(res{4});
- if isempty(bgLevel)
- bgLevel = sMri.Histogram.bgLevel;
+
+ % First, fill thin holes with boolean kernel convolution
+ % This fills voxels surrounded by 1s with specific patterns.
+ mask = KernelClean(mask, true);
+
+ % Main filling step
+ % Fill from first to last 1-voxels along each direction and keep intersection across 3 dimensions.
+ % This can result in very narrow deep "tunnels", which we fix with the "surround" fill again next.
+ mask = (Fill(mask, 1, true) & Fill(mask, 2, true) & Fill(mask, 3, true));
+
+ % "Surround" and "sandwich" fills again
+ mask = KernelClean(mask, true);
+ % Repeat for filling intersections. Twice ok for 2d or 3d.
+ mask = KernelClean(mask, true);
+
+ if isClean
+ % Apply inverse "surround" to remove noise and small protrusions.
+ % Erase voxel if surrounded by 0 in a plane. Could do before "first to last" above to clean
+ % noise first, but could also erase more parts in low SNR areas. We later keep only
+ % connected central part so no need to iterate this step.
+ mask = KernelClean(mask, false);
end
end
-% Check parameters values
-if isempty(nVertices) || (nVertices < 50) || (nVertices ~= round(nVertices)) || isempty(erodeFactor) || ~ismember(erodeFactor,[0,1,2,3]) || isempty(fillFactor) || ~ismember(fillFactor,[0,1,2,3])
- bst_error('Invalid parameters.', 'Head surface', 0);
- return
-end
+function mask = KernelClean(mask, Value)
+ % Boolean convolution, with predetermined kernels, to fill in thin strands or cracks.
+ % If Value = false, inverts the mask before and after filling, essentially eroding away thin
+ % structures instead of filling thin holes. Works for 3d volume or 2d slice, with a different
+ % set of kernel patterns for each.
+ if nargin < 2 || isempty(Value)
+ Value = true;
+ end
+ % Flip the mask if we're looking for false.
+ if ~Value
+ mask = ~mask;
+ end
-%% ===== CREATE HEAD MASK =====
-% Progress bar
-bst_progress('start', 'Generate head surface', 'Creating head mask...', 0, 100);
-% Threshold mri to the level estimated in the histogram
-headmask = (sMri.Cube(:,:,:,1) > bgLevel);
-% Closing all the faces of the cube
-headmask(1,:,:) = 0*headmask(1,:,:);
-headmask(end,:,:) = 0*headmask(1,:,:);
-headmask(:,1,:) = 0*headmask(:,1,:);
-headmask(:,end,:) = 0*headmask(:,1,:);
-headmask(:,:,1) = 0*headmask(:,:,1);
-headmask(:,:,end) = 0*headmask(:,:,1);
-% Erode + dilate, to remove small components
-if (erodeFactor > 0)
- headmask = headmask & ~mri_dilate(~headmask, erodeFactor);
- headmask = mri_dilate(headmask, erodeFactor);
-end
-bst_progress('inc', 10);
-% Fill holes
-bst_progress('text', 'Filling holes...');
-headmask = (mri_fillholes(headmask, 1) & mri_fillholes(headmask, 2) & mri_fillholes(headmask, 3));
-bst_progress('inc', 10);
-
-% view_mri_slices(headmask, 'x', 20)
-
-
-%% ===== CREATE SURFACE =====
-% Compute isosurface
-bst_progress('text', 'Creating isosurface...');
-[sHead.Faces, sHead.Vertices] = mri_isosurface(headmask, 0.5);
-bst_progress('inc', 10);
-% Downsample to a maximum number of vertices
-maxIsoVert = 60000;
-if (length(sHead.Vertices) > maxIsoVert)
- bst_progress('text', 'Downsampling isosurface...');
- [sHead.Faces, sHead.Vertices] = reducepatch(sHead.Faces, sHead.Vertices, maxIsoVert./length(sHead.Vertices));
- bst_progress('inc', 10);
-end
-% Remove small objects
-bst_progress('text', 'Removing small patches...');
-[sHead.Vertices, sHead.Faces] = tess_remove_small(sHead.Vertices, sHead.Faces);
-bst_progress('inc', 10);
-
-% Downsampling isosurface
-if (length(sHead.Vertices) > nVertices)
- bst_progress('text', 'Downsampling surface...');
- [sHead.Faces, sHead.Vertices] = reducepatch(sHead.Faces, sHead.Vertices, nVertices./length(sHead.Vertices));
- bst_progress('inc', 10);
+ % Dimensions that have thickness, enough for convolution with 3-wide kernel.
+ isThk = size(mask, [1,2,3]) > 2;
+ nD = sum(isThk);
+ if nD == 1
+ error('Unexpected 1d "volume".');
+ end
+
+ % All kernels should have mirror symmetry on one plane through the middle, since they're only
+ % applied in one orientation for each dimension. For now only 3x3x3 kernels, but should work
+ % with larger "square" ones as well.
+ switch nD
+ case 3
+ % Kernels for 3d
+ % "Surround": to get rid of thin tunnels, "hairs"
+ % 4 adjacent voxels in a plane (not diagonals) are 1s.
+ Kernels{1} = zeros(3,3,3); Kernels{1}(:,:,2) = [0,1,0;1,0,1;0,1,0];
+ % "Sandwich"/"Tie fighter": to remove thin cracks, wedges
+ % both side planes (3x3) are 1s
+ Kernels{2} = ones(3,3,3); Kernels{2}(:,:,2) = zeros(3,3);
+ case 2
+ % Kernels for 2d
+ % "Surround"/"sandwich" (same for 2d)
+ Kernels{1} = zeros(3,3); Kernels{1}(2,:) = [1,0,1];
+ otherwise
+ error('Unexpected multi-dim (>3) mask.');
+ end
+ nK = numel(Kernels);
+
+ switch nD
+ case 3
+ for iK = 1:nK
+ % To avoid cumulative effects that would depend on the order in which we orient the
+ % kernel, we apply all dimensions on the original mask before "adding" them (with "or").
+ FilledMask = mask;
+ for iD = 1:3
+ % Re-orient kernel along each dimension
+ K = permute(Kernels{iK}, circshift(1:3, iD-1));
+ N = sum(K(:));
+ FilledMask = FilledMask | convn(mask, K, 'same') == N;
+ end
+ % Because of the shapes of our kernels, we don't have to enforce keeping "zero" on
+ % our mask volume boundary slices.
+ % mask(2:end-1,2:end-1,2:end-1) = FilledMask(2:end-1,2:end-1,2:end-1);
+ mask = FilledMask;
+ end
+ case 2
+ % Rotate mask to have flat dimension last
+ iD = 1:3;
+ Perm = [iD(isThk), iD(~isThk)]; % This is not always a single swap, so we need the inverse.
+ [~, PermInv] = sort(Perm);
+ mask = permute(mask, Perm);
+ for iK = 1:nK
+ FilledMask = mask;
+ for iD = 1:2
+ % Re-orient kernel along each dimension
+ if iD == 1
+ K = Kernels{iK};
+ else % iD == 2 permute
+ K = Kernels{iK}';
+ end
+ N = sum(K(:));
+ FilledMask = FilledMask | convn(mask, K, 'same') == N;
+ end
+ mask = FilledMask;
+ end
+ mask = permute(mask, PermInv);
+ end
+ % Inverse mask back if needed.
+ if ~Value
+ mask = ~mask;
+ end
end
-% Convert to millimeters
-sHead.Vertices = sHead.Vertices(:,[2,1,3]);
-sHead.Faces = sHead.Faces(:,[2,1,3]);
-sHead.Vertices = bst_bsxfun(@times, sHead.Vertices, sMri.Voxsize);
-% Convert to SCS
-sHead.Vertices = cs_convert(sMri, 'mri', 'scs', sHead.Vertices ./ 1000);
-
-% Reduce the final size of the meshed volume
-erodeFinal = 3;
-% Fill holes in surface
-%if (fillFactor > 0)
- bst_progress('text', 'Filling holes...');
- [sHead.Vertices, sHead.Faces] = tess_fillholes(sMri, sHead.Vertices, sHead.Faces, fillFactor, erodeFinal);
- bst_progress('inc', 30);
+
+
+% function mask = Surrounded(mask, Value)
+% % Find voxels that are surrounded by a value and add them (if Value=true) or remove them (false).
+% % 4 adjacent voxels in a plane (not diagonals) for 3d, 2 adjacent voxels in a line for 2d.
+%
+% if nargin < 2 || isempty(Value)
+% Value = true;
+% end
+% % Flip the mask if we're looking for false.
+% if ~Value
+% mask = ~mask;
+% end
+%
+% % Indices for dimensions excluding ends (2:end-1), except if thin dimension, then it's just 1 or [1 2].
+% nVox = size(mask, [1,2,3]);
+% iVox = {min(2, nVox(1)):max(nVox(1)-1, 1), min(2, nVox(2)):max(nVox(2)-1, 1), min(2, nVox(3)):max(nVox(3)-1, 1)};
+% S = zeros(nVox - (nVox > 2)*2);
+%
+% % Loop so it works on 2d or 3d
+% nDim = 0;
+% for iDim = 1:3
+% if size(mask, iDim) > 2 % Skip singleton dimensions
+% nDim = nDim + 1;
+% switch iDim
+% case 1
+% S = S + (mask(1:end-2,iVox{2},iVox{3}) & mask(3:end,iVox{2},iVox{3}));
+% case 2
+% S = S + (mask(iVox{1},1:end-2,iVox{3}) & mask(iVox{1},3:end,iVox{3}));
+% case 3
+% S = S + (mask(iVox{1},iVox{2},1:end-2) & mask(iVox{1},iVox{2},3:end));
+% end
+% end
+% end
+%
+% % Modify original mask by adding (true) or removing (false)
+% mask(iVox{1},iVox{2},iVox{3}) = mask(iVox{1},iVox{2},iVox{3}) | S >= nDim - 1;
+% if ~Value
+% mask = ~mask;
+% end
% end
-%% ===== SAVE FILES =====
-if isSave
- bst_progress('text', 'Saving new file...');
- % Create output filenames
- ProtocolInfo = bst_get('ProtocolInfo');
- SurfaceDir = bst_fullfile(ProtocolInfo.SUBJECTS, bst_fileparts(MriFile));
- HeadFile = file_unique(bst_fullfile(SurfaceDir, 'tess_head_mask.mat'));
- % Save head
- if ~isempty(Comment)
- sHead.Comment = Comment;
- else
- sHead.Comment = sprintf('head mask (%d,%d,%d,%d)', nVertices, erodeFactor, fillFactor, round(bgLevel));
- end
- sHead = bst_history('add', sHead, 'bem', 'Head surface generated with Brainstorm');
- bst_save(HeadFile, sHead, 'v7');
- iSurface = db_add_surface( iSubject, HeadFile, sHead.Comment);
-else
- % Return surface
- HeadFile = sHead.Vertices;
- iSurface = sHead.Faces;
+function mask = Fill(mask, dim, isFullSingl)
+ % Modified to exclude boundaries, so we can get rid of external junk as well as
+ % internal holes easily.
+ if nargin < 3 || isempty(isFullSingl)
+ % Return original mask for singleton dim by default.
+ isFullSingl = false;
+ end
+
+ % Initialize two accumulators, for the two directions
+ acc1 = mask;
+ acc2 = mask;
+ n = size(mask,dim);
+ % Skip singleton dimensions
+ if n == 1
+ if isFullSingl
+ % Return all true, e.g. to combine with Fill in other directions.
+ mask = true(size(mask));
+ end
+ return;
+ end
+
+ % Process in required direction
+ switch dim
+ case 1
+ for i = 2:n
+ acc1(i,:,:) = acc1(i,:,:) | acc1(i-1,:,:);
+ end
+ for i = n-1:-1:1
+ acc2(i,:,:) = acc2(i,:,:) | acc2(i+1,:,:);
+ end
+ case 2
+ for i = 2:n
+ acc1(:,i,:) = acc1(:,i,:) | acc1(:,i-1,:);
+ end
+ for i = n-1:-1:1
+ acc2(:,i,:) = acc2(:,i,:) | acc2(:,i+1,:);
+ end
+ case 3
+ for i = 2:n
+ acc1(:,:,i) = acc1(:,:,i) | acc1(:,:,i-1);
+ end
+ for i = n-1:-1:1
+ acc2(:,:,i) = acc2(:,:,i) | acc2(:,:,i+1);
+ end
+ end
+ % Combine two accumulators
+ mask = acc1 & acc2;
end
-% Close, success
-if isProgress
- bst_progress('stop');
+
+% function mask = Dilate(mask)
+% % Dilate by 1 voxel in 6 directions, except at volume edges
+% % Indices for dimensions excluding ends (2:end-1), except if thin dimension, then it's just 1 or [1 2].
+% nVox = size(mask, [1,2,3]);
+% iVox = {min(2, nVox(1)):max(nVox(1)-1, 1), min(2, nVox(2)):max(nVox(2)-1, 1), min(2, nVox(3)):max(nVox(3)-1, 1)};
+%
+% % Loop so it works on 2d or 3d
+% for iDim = 1:3
+% if nVox(iDim) > 2 % Skip thin dimensions (size 1 or 2)
+% switch iDim
+% case 1
+% DilateMask = mask(1:end-2,iVox{2},iVox{3}) | mask(3:end,iVox{2},iVox{3});
+% case 2
+% DilateMask = mask(iVox{1},1:end-2,iVox{3}) | mask(iVox{1},3:end,iVox{3});
+% case 3
+% DilateMask = mask(iVox{1},iVox{2},1:end-2) | mask(iVox{1},iVox{2},3:end);
+% end
+% mask(iVox{1},iVox{2},iVox{3}) = mask(iVox{1},iVox{2},iVox{3}) | DilateMask;
+% end
+% end
+% end
+
+
+function [OutMask, isFail] = CenterSpread(InMask)
+ % Similar to Fill, but from a central starting point and intersecting with the input "reference" mask.
+ % This should work on slices as well as volumes.
+ isFail = false;
+ OutMask = false(size(InMask));
+ iStart = max(1,round(size(OutMask)/2));
+ nVox = size(OutMask);
+ % Force starting center point to be 1, and spread from there. But this will still fail if it's fully
+ % surrounded by 0s.
+ OutMask(iStart(1), iStart(2), iStart(3)) = true;
+ nPrev = 0;
+ nOut = 1;
+ while nOut > nPrev
+ % Dilation loop was very slow.
+ % OutMask = OutMask | (Dilate(OutMask) & InMask);
+ % Instead, propagate as far as possible in each direction (3 dim, forward & back) at each step
+ % of the main loop.
+ for x = 2:nVox(1)
+ OutMask(x,:,:) = OutMask(x,:,:) | (OutMask(x-1,:,:) & InMask(x,:,:));
+ end
+ for x = nVox(1)-1:-1:1
+ OutMask(x,:,:) = OutMask(x,:,:) | (OutMask(x+1,:,:) & InMask(x,:,:));
+ end
+ for y = 2:nVox(2)
+ OutMask(:,y,:) = OutMask(:,y,:) | (OutMask(:,y-1,:) & InMask(:,y,:));
+ end
+ for y = nVox(2)-1:-1:1
+ OutMask(:,y,:) = OutMask(:,y,:) | (OutMask(:,y+1,:) & InMask(:,y,:));
+ end
+ for z = 2:nVox(3)
+ OutMask(:,:,z) = OutMask(:,:,z) | (OutMask(:,:,z-1) & InMask(:,:,z));
+ end
+ for z = nVox(3)-1:-1:1
+ OutMask(:,:,z) = OutMask(:,:,z) | (OutMask(:,:,z+1) & InMask(:,:,z));
+ end
+ nPrev = nOut;
+ nOut = sum(OutMask(:));
+ end
+ if nOut == 1
+ % Remove "forced" initial vertex, everything else is now gone.
+ OutMask(iStart(1), iStart(2), iStart(3)) = false;
+ isFail = true;
+ % warning('CenterSpread failed: starting center point is not part of the mask.');
+ end
end
+function [Vol, Vect] = NormGradient(Vol)
+ % Norm of the spatial gradient vector field in a regular 3D volume.
+ [x,y,z] = gradient(Vol);
+ Vect = cat(4,x,y,z);
+ Vol = sqrt(sum(Vect.^2, 4));
+end
diff --git a/toolbox/core/bst_colormaps.m b/toolbox/core/bst_colormaps.m
index ad1078a205..8d0d448663 100644
--- a/toolbox/core/bst_colormaps.m
+++ b/toolbox/core/bst_colormaps.m
@@ -376,25 +376,30 @@ function SetMaxCustom(ColormapType, DisplayUnits, newMin, newMax)
case {'3DViz', 'MriViewer'}
% Get surfaces defined in this figure
TessInfo = getappdata(sFigure.hFigure, 'Surface');
- % Find 1st surface that match this ColormapType
- iTess = find(strcmpi({TessInfo.ColormapType}, ColormapType), 1);
+ % Find surfaces that match this ColormapType
+ iSurfaces = find(strcmpi({TessInfo.ColormapType}, ColormapType));
DataFig = [];
- if ~isempty(iTess) && ~isempty(TessInfo(iTess).DataSource.Type)
- DataFig = TessInfo(iTess).DataMinMax;
- DataType = TessInfo(iTess).DataSource.Type;
- % For Data: use the modality instead
- if strcmpi(DataType, 'Data') && ~isempty(ColormapInfo.Type) && ismember(ColormapInfo.Type, {'eeg', 'meg', 'nirs'})
- DataType = upper(ColormapInfo.Type);
- % sLORETA: Do not use regular source scaling (pAm)
- elseif strcmpi(DataType, 'Source')
- if ~isempty(strfind(lower(TessInfo(iTess).DataSource.FileName), 'sloreta'))
- DataType = 'sLORETA';
- elseif ~isempty(strfind(lower(TessInfo(iTess).DataSource.FileName), 'headmodel'))
- DataType = 'headmodel';
- else
- [~, iResult] = bst_memory('LoadResultsFile', TessInfo(iTess).DataSource.FileName, 0);
- if ~isempty(strfind(lower(GlobalData.DataSet(iDS).Results(iResult).Function), 'sloreta'));
+ for i = 1:length(iSurfaces)
+ iTess = iSurfaces(i);
+ if ~isempty(TessInfo(iTess).DataSource.Type)
+ DataFig = [min([DataFig(:); TessInfo(iTess).DataMinMax(:)]), ...
+ max([DataFig(:); TessInfo(iTess).DataMinMax(:)])];
+ % We'll keep the last non-empty DataType
+ DataType = TessInfo(iTess).DataSource.Type;
+ % For Data: use the modality instead
+ if strcmpi(DataType, 'Data') && ~isempty(ColormapInfo.Type) && ismember(ColormapInfo.Type, {'eeg', 'meg', 'nirs'})
+ DataType = upper(ColormapInfo.Type);
+ % sLORETA: Do not use regular source scaling (pAm)
+ elseif strcmpi(DataType, 'Source')
+ if ~isempty(strfind(lower(TessInfo(iTess).DataSource.FileName), 'sloreta'))
DataType = 'sLORETA';
+ elseif ~isempty(strfind(lower(TessInfo(iTess).DataSource.FileName), 'headmodel'))
+ DataType = 'headmodel';
+ else
+ [~, iResult] = bst_memory('LoadResultsFile', TessInfo(iTess).DataSource.FileName, 0);
+ if ~isempty(strfind(lower(GlobalData.DataSet(iDS).Results(iResult).Function), 'sloreta'));
+ DataType = 'sLORETA';
+ end
end
end
end
@@ -1338,7 +1343,7 @@ function SetColormapRealMin(ColormapType, status)
% Fire change notificiation to all figures (3DViz and Topography)
FireColormapChanged(ColormapType);
end
-function SetMaxMode(ColormapType, maxmode, DisplayUnits)
+function SetMaxMode(ColormapType, maxmode, DisplayUnits, varargin)
% Parse inputs
if (nargin < 3) || isempty(DisplayUnits)
DisplayUnits = [];
@@ -1349,7 +1354,7 @@ function SetMaxMode(ColormapType, maxmode, DisplayUnits)
end
% Custom: ask for custom values
if strcmpi(maxmode, 'custom')
- SetMaxCustom(ColormapType, DisplayUnits);
+ SetMaxCustom(ColormapType, DisplayUnits, varargin{:});
else
% Update colormap
sColormap = GetColormap(ColormapType);
diff --git a/toolbox/db/db_set_channel.m b/toolbox/db/db_set_channel.m
index d5d0764bc9..d866171e49 100644
--- a/toolbox/db/db_set_channel.m
+++ b/toolbox/db/db_set_channel.m
@@ -183,12 +183,10 @@
[ChannelMat, R, T, isSkip, isUserCancel, strReport, Tolerance] = channel_align_auto(OutputFile, [], 0, isConfirm, Tolerance);
% User validated: keep this answer for the next round (force alignment for next call)
if ~isSkip
- if isUserCancel
+ if isUserCancel || isempty(ChannelMat)
ChannelAlign = 0;
- elseif ~isempty(ChannelMat)
+ elseif ChannelAlign < 2
ChannelAlign = 2;
- else
- ChannelAlign = 0;
end
end
end
diff --git a/toolbox/gui/figure_3d.m b/toolbox/gui/figure_3d.m
index 384b261374..658e97899a 100644
--- a/toolbox/gui/figure_3d.m
+++ b/toolbox/gui/figure_3d.m
@@ -4094,18 +4094,22 @@ function ViewAxis(hFig, isVisible)
isVisible = isempty(findobj(hAxes, 'Tag', 'AxisXYZ'));
end
if isVisible
+ % Works but now default zoom is too tight >:(
% Get dimensions of current axes
+ set(hAxes, 'XLimitMethod', 'tight', 'YLimitMethod', 'tight', 'ZLimitMethod', 'tight');
XLim = get(hAxes, 'XLim');
- YLim = get(hAxes, 'XLim');
- ZLim = get(hAxes, 'XLim');
- d = max(abs([XLim(:); YLim(:); ZLim(:)]));
+ YLim = get(hAxes, 'YLim');
+ ZLim = get(hAxes, 'ZLim');
% Draw axis lines
+ d = XLim(2) * 1.04;
line([0 d], [0 0], [0 0], 'Color', [1 0 0], 'Marker', '>', 'Parent', hAxes, 'Tag', 'AxisXYZ');
+ text(d * 1.04, 0, 0, 'X', 'Color', [1 0 0], 'Parent', hAxes, 'Tag', 'AxisXYZ');
+ d = YLim(2) * 1.04;
line([0 0], [0 d], [0 0], 'Color', [0 1 0], 'Marker', '>', 'Parent', hAxes, 'Tag', 'AxisXYZ');
+ text(0, d * 1.04, 0, 'Y', 'Color', [0 1 0], 'Parent', hAxes, 'Tag', 'AxisXYZ');
+ d = ZLim(2) * 1.04;
line([0 0], [0 0], [0 d], 'Color', [0 0 1], 'Marker', '>', 'Parent', hAxes, 'Tag', 'AxisXYZ');
- text(d+0.002, 0, 0, 'X', 'Color', [1 0 0], 'Parent', hAxes, 'Tag', 'AxisXYZ');
- text(0, d+0.002, 0, 'Y', 'Color', [0 1 0], 'Parent', hAxes, 'Tag', 'AxisXYZ');
- text(0, 0, d+0.002, 'Z', 'Color', [0 0 1], 'Parent', hAxes, 'Tag', 'AxisXYZ');
+ text(0, 0, d * 1.04, 'Z', 'Color', [0 0 1], 'Parent', hAxes, 'Tag', 'AxisXYZ');
% Enforce camera target at (0,0,0)
% camtarget(hAxes, [0,0,0]);
else
diff --git a/toolbox/gui/figure_mri.m b/toolbox/gui/figure_mri.m
index 9a9d66f7d8..e8d3bd764b 100644
--- a/toolbox/gui/figure_mri.m
+++ b/toolbox/gui/figure_mri.m
@@ -2537,10 +2537,22 @@ function ButtonSave_Callback(hFig, varargin)
%% ===== SAVE MRI =====
+% Update MRI fiducials and adjust surfaces.
+% Input can be figure handle or sMri.
function [isCloseAccepted, MriFile] = SaveMri(hFig)
ProtocolInfo = bst_get('ProtocolInfo');
% Get MRI
- sMri = panel_surface('GetSurfaceMri', hFig);
+ if ishandle(hFig)
+ sMri = panel_surface('GetSurfaceMri', hFig);
+ isUser = true;
+ elseif isstruct(hFig)
+ sMri = hFig;
+ isUser = false;
+ else
+ bst_error('SaveMri: Unexpected input: %s.', class(hFig));
+ isCloseAccepted = 0;
+ return;
+ end
MriFile = sMri.FileName;
MriFileFull = bst_fullfile(ProtocolInfo.SUBJECTS, MriFile);
% Do not accept "Save" if user did not select all the fiducials
@@ -2558,21 +2570,26 @@ function ButtonSave_Callback(hFig, varargin)
warning('off', 'MATLAB:load:variableNotFound');
sMriOld = load(MriFileFull, 'SCS');
warning('on', 'MATLAB:load:variableNotFound');
- % If the fiducials were modified
+ % Check if the fiducials were modified (> 1um), to realign surfaces below
if isfield(sMriOld, 'SCS') && all(isfield(sMriOld.SCS,{'NAS','LPA','RPA'})) ...
&& ~isempty(sMriOld.SCS.NAS) && ~isempty(sMriOld.SCS.LPA) && ~isempty(sMriOld.SCS.RPA) ...
- && ((max(sMri.SCS.NAS - sMriOld.SCS.NAS) > 1e-3) || ...
- (max(sMri.SCS.LPA - sMriOld.SCS.LPA) > 1e-3) || ...
- (max(sMri.SCS.RPA - sMriOld.SCS.RPA) > 1e-3))
+ && ((max(abs(sMri.SCS.NAS - sMriOld.SCS.NAS)) > 1e-3) || ...
+ (max(abs(sMri.SCS.LPA - sMriOld.SCS.LPA)) > 1e-3) || ...
+ (max(abs(sMri.SCS.RPA - sMriOld.SCS.RPA)) > 1e-3))
% Nothing to do...
else
+ % sMri.SCS.R, T and Origin are updated before calling this function.
sMriOld = [];
end
% === HISTORY ===
% History: Edited the fiducials
if ~isfield(sMriOld, 'SCS') || ~isequal(sMriOld.SCS, sMri.SCS) || ~isfield(sMriOld, 'NCS') || ~isequal(sMriOld.NCS, sMri.NCS)
- sMri = bst_history('add', sMri, 'edit', 'User edited the fiducials');
+ if isUser
+ sMri = bst_history('add', sMri, 'edit', 'User edited the fiducials');
+ else
+ sMri = bst_history('add', sMri, 'edit', 'Applied digitized anatomical fiducials'); % string used to verify elsewhere
+ end
end
% ==== SAVE MRI ====
diff --git a/toolbox/gui/view_channels.m b/toolbox/gui/view_channels.m
index 074b8fdf8b..d22e2bfdbf 100644
--- a/toolbox/gui/view_channels.m
+++ b/toolbox/gui/view_channels.m
@@ -189,7 +189,7 @@
selChannels = setdiff(selChannels, iNoLoc);
end
if isempty(selChannels) && ~isempty(iNoLoc)
- error('None of the channels have positions defined.');
+ error('None of the channels have positions defined for modality %s.', Modality);
end
end
% Set figure selected sensors
diff --git a/toolbox/gui/view_mri_histogram.m b/toolbox/gui/view_mri_histogram.m
index b2a13ff586..cd2043c4e5 100644
--- a/toolbox/gui/view_mri_histogram.m
+++ b/toolbox/gui/view_mri_histogram.m
@@ -1,10 +1,12 @@
-function hFig = view_mri_histogram( MriFile )
+function hFig = view_mri_histogram(MriFile, isInteractive)
% VIEW_MRI_HISTOGRAM: Compute and view the histogram of a brainstorm MRI.
%
-% USAGE: hFig = view_mri_histogram( MriFile );
+% USAGE: hFig = view_mri_histogram(MriFile, isInteractive=false);
%
% INPUT:
% - MriFile : Full path to a brainstorm MRI file
+% - isInteractive : If true, clicking on the figure will update the background threshold value,
+% and offer saving when closing the figure.
% OUTPUT:
% - hFig : Matlab handle to the figure where the histogram is displayed
% @=============================================================================
@@ -25,19 +27,30 @@
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
-% Authors: Francois Tadel, 2006-2020
+% Authors: Francois Tadel, 2006-2020, Marc Lalancette 2025
-%% ===== COMPUTE HISTOGRAM =====
+if nargin < 2 || isempty(isInteractive)
+ isInteractive = false;
+end
+
+%% ===== LOAD OR COMPUTE HISTOGRAM =====
% Display progress bar
bst_progress('start', 'View MRI historgram', 'Computing histogram...');
-% Load full MRI
-MRI = load(MriFile);
-% Compute histogram
-Histogram = mri_histogram(MRI.Cube(:,:,:,1));
-% Save histogram
-s.Histogram = Histogram;
-bst_save(MriFile, s, 'v7', 1);
-
+if isstruct(MriFile) && isfield(MriFile, 'intensityMax')
+ Histogram = MriFile;
+else
+ % Load full MRI
+ MRI = load(MriFile);
+ % Compute histogram if missing
+ if ~isfield(MRI, 'Histogram') || isempty(MRI.Histogram) || ~isfield(MRI.Histogram, 'intensityMax')
+ Histogram = mri_histogram(MRI.Cube(:,:,:,1));
+ % Save histogram
+ s.Histogram = Histogram;
+ bst_save(MriFile, s, 'v7', 1); % isAppend
+ else
+ Histogram = MRI.Histogram;
+ end
+end
%% ===== DISPLAY HISTOGRAM =====
% Create figure
@@ -81,17 +94,56 @@
yLimits = [0 maxVal*1.3];
ylim(yLimits);
% Display background and white matter thresholds
-line([Histogram.bgLevel, Histogram.bgLevel], yLimits, 'Color','b');
+% Keep original level to check if it changes.
+bgLevel = Histogram.bgLevel;
+hBg = line([Histogram.bgLevel, Histogram.bgLevel], yLimits, 'Color','b');
line([Histogram.whiteLevel, Histogram.whiteLevel], yLimits, 'Color','y');
h = legend('MRI hist.','Smoothed hist.','Cumulative hist.','Maxima','Minima',...
'Scalp or grey thresh.','White m thresh.');
set(h, 'FontSize', bst_get('FigFont'), ...
'FontUnits', 'points');
+% Set interactive callbacks
+if isInteractive
+ set(hAxes, 'ButtonDownFcn', @clickCallback);
+ set(hFig, 'CloseRequestFcn', @(src, event) closeFigureCallback());
+end
+
% Hide progress bar
bst_progress('stop');
+function clickCallback(~, event)
+ % Extract the x-coordinate of the click
+ bgLevel = event.IntersectionPoint(1);
+ % fprintf('Clicked at x = %.4f\n', x);
+ if bgLevel ~= Histogram.bgLevel
+ set(hBg, 'xdata', [bgLevel, bgLevel]);
+ % drawnow % needed?
+ end
+end
+
+function closeFigureCallback()
+ if bgLevel ~= Histogram.bgLevel
+ % Request save confirmation.
+ [Proceed, isCancel] = java_dialog('confirm', sprintf(...
+ 'MRI background intensity threshold changed (%d > %d). Save?', round(Histogram.bgLevel), round(bgLevel)), ...
+ 'MRI background threshold');
+ if isCancel
+ return;
+ end
+ if Proceed
+ % Save histogram
+ Histogram.bgLevel = bgLevel;
+ s.Histogram = Histogram;
+ bst_save(MriFile, s, 'v7', 1); % isAppend
+ end
+ end
+ % Close figure
+ delete(hFig);
+end
+
+end
diff --git a/toolbox/io/bst_save_coregistration.m b/toolbox/io/bst_save_coregistration.m
new file mode 100644
index 0000000000..733d1dfa90
--- /dev/null
+++ b/toolbox/io/bst_save_coregistration.m
@@ -0,0 +1,604 @@
+function [isSuccess, OutFilesMri, OutFilesMeg] = bst_save_coregistration(iStudies, isBids, ...
+ RecreateMegCoordJson, isOverwrite, isDryRun)
+ % Save MRI-MEG coregistration info in imported raw BIDS dataset, or MRI fiducials only if not BIDS.
+ % IMPORTANT: isSuccess currently not working, can often be true despite skipping studies or subjects.
+ %
+ % [isSuccess, OutFilesMri, OutFilesMeg] = bst_save_coregistration(iStudies, isBids=)
+ %
+ % Save MRI-MEG coregistration by adding AnatomicalLandmarkCoordinates to the
+ % _T1w.json MRI metadata, in 0-indexed voxel coordinates, and to the
+ % _coordsystem.json files for functional data, in native coordinates (e.g. CTF).
+ % The points used are the anatomical fiducials marked in Brainstorm on the MRI
+ % that define the Brainstorm subject coordinate system (SCS).
+ %
+ % iStudies is a list of Brainstorm study indices. Shown e.g. in pop-up when hovering over a
+ % study folder (or data file within it) in the tree.
+ %
+ % RecreateMegCoordJson: if true, this file will be deleted and recreated with
+ % BidsBuildRecordingFiles before adding coregistration info to it (was added because fiducial
+ % descriptions were updated).
+ %
+ % If the raw data is not BIDS, the anatomical fiducials are saved in a
+ % fiducials.m file next to the raw MRI file, in Brainstorm MRI coordinates.
+ %
+ % Discussion about saving MRI-MEG coregistration in BIDS:
+ % https://groups.google.com/g/bids-discussion/c/BeyUeuNGl7I
+
+ % TODO: some dependencies from BIDS code
+ % TODO: apply isOverwrite to other things?, for now, only MRI json
+
+ % This could become an option, depending on where this process is called form.
+ isInteractive = true;
+ % Another potential option, to back up json files under derivatives
+ isBackup = true;
+
+ if nargin < 5 || isempty(isDryRun)
+ isDryRun = false;
+ end
+ if nargin < 4 || isempty(isOverwrite)
+ isOverwrite = true;
+ end
+ if nargin < 3 || isempty(RecreateMegCoordJson)
+ RecreateMegCoordJson = false;
+ end
+ if nargin < 2 || isempty(isBids)
+ isBids = [];
+ end
+ sSubjects = bst_get('ProtocolSubjects');
+ if nargin < 1 || isempty(iStudies)
+ % Try to get all subjects from currently loaded protocol.
+ nSub = numel(sSubjects.Subject);
+ iSubjects = 1:nSub;
+ sStudies = [];
+ else
+ sStudies = bst_get('Study', iStudies);
+ for iiStudy = 1:numel(iStudies)
+ [~, iSubForStudies(iiStudy)] = bst_get('Subject', sStudies(iiStudy).BrainStormSubject);
+ end
+ iSubjects = unique(iSubForStudies);
+ nSub = numel(iSubjects);
+ end
+
+ bst_progress('start', 'Save co-registration', ' ', 0, nSub);
+
+ % Avoid useless warnings. I don't understand why they warn about this, a cell of char vectors is what I want.
+ warning('off', 'MATLAB:table:PreallocateCharWarning');
+
+ OutFilesMri = cell(nSub, 1);
+ OutFilesMeg = cell(nSub, 1);
+ isSuccess = false(nSub, 1);
+ BidsRoot = '';
+ for iOutSub = 1:nSub
+ iSub = iSubjects(iOutSub);
+ fprintf('%3d %s\n', iSub, sSubjects.Subject(iSub).Name);
+ % Get anatomical file.
+ if ~contains(sSubjects.Subject(iSub).Anatomy(sSubjects.Subject(iSub).iAnatomy).Comment, 'MRI', 'ignorecase', true) && ...
+ ~contains(sSubjects.Subject(iSub).Anatomy(sSubjects.Subject(iSub).iAnatomy).Comment, 't1w', 'ignorecase', true)
+ warning('Selected anatomy is not ''MRI''. Skipping subject %s.', sSubjects.Subject(iSub).Name);
+ continue;
+ end
+ sMri = load(file_fullpath(sSubjects.Subject(iSub).Anatomy(sSubjects.Subject(iSub).iAnatomy).FileName));
+ ImportedFile = strrep(sMri.History{1,3}, 'Import from: ', '');
+ if ~exist(ImportedFile, 'file')
+ warning('Imported anatomy file not found. Skipping subject %s.', sSubjects.Subject(iSub).Name);
+ continue;
+ end
+ % Get studies for this subject, only those provided (unless none, then we'll check all)
+ if ~isempty(sStudies)
+ % sStudies here contains only the requested ones.
+ sStudiesForSub = sStudies(iSubForStudies == iSub);
+ else
+ % Get all linked raw data files.
+ sStudiesForSub = bst_get('StudyWithSubject', sSubjects.Subject(iSub).FileName);
+ end
+ if isempty(isBids)
+ % Try to find root BIDS folder.
+ BidsRoot = bst_fileparts(bst_fileparts(ImportedFile)); % go back through "anat" and subject folders at least (session not mandatory).
+ isBids = true; % changed if not found below
+ while ~exist(fullfile(BidsRoot, 'dataset_description.json'), 'file')
+ if isempty(BidsRoot) || ~exist(BidsRoot, 'dir')
+ isBids = false;
+ fprintf('BST> bst_save_coregistration detected that raw imported data is NOT structured as BIDS.');
+ break;
+ end
+ BidsRoot = bst_fileparts(BidsRoot);
+ end
+ end
+ if isBids
+ if isempty(BidsRoot)
+ BidsRoot = bst_fileparts(bst_fileparts(ImportedFile)); % go back through "anat" and subject folders at least (session not mandatory).
+ while ~exist(fullfile(BidsRoot, 'dataset_description.json'), 'file')
+ if isempty(BidsRoot) || ~exist(BidsRoot, 'dir')
+ error('Cannot find BIDS root folder and dataset_description.json file; subject %s.', sSubjects.Subject(iSub).Name);
+ end
+ BidsRoot = bst_fileparts(BidsRoot);
+ end
+ end
+
+ % MRI _t1w.json
+ % Save anatomical landmarks in Nifti voxel coordinates
+ [MriPath, MriName, MriExt] = bst_fileparts(ImportedFile);
+ if strcmpi(MriExt, '.gz')
+ [~, MriName, MriExt2] = fileparts(MriName);
+ MriExt = [MriExt2, MriExt]; %#ok
+ end
+ if ~strncmpi(MriExt, '.nii', 4)
+ warning('Imported anatomy not BIDS. Skipping subject %s.', sSubjects.Subject(iSub).Name);
+ continue;
+ end
+ MriJsonFile = fullfile(MriPath, [MriName, '.json']);
+ if ~exist(MriJsonFile, 'file')
+ warning('Imported anatomy BIDS json file not found for subject %s. Creating new file.', sSubjects.Subject(iSub).Name);
+ sMriJson = struct();
+ else
+ sMriJson = bst_jsondecode(MriJsonFile, false);
+ % Make backup in derivatives folder.
+ if isBackup
+ BakMriJsonFile = replace(MriJsonFile, BidsRoot, fullfile(BidsRoot, 'derivatives'));
+ if ~exist(BakMriJsonFile, 'file')
+ BakFolder = fileparts(BakMriJsonFile);
+ if ~exist(BakFolder, 'dir')
+ [isOk, Msg] = mkdir(BakFolder);
+ if ~isOk, warning(Msg); end
+ if ~exist(BakFolder, 'dir')
+ warning('Unable to create backup folder %s. Skipping subject %s.', BakFolder, sSubjects.Subject(iSub).Name);
+ continue;
+ end
+ end
+ [isOk, Msg] = copyfile(MriJsonFile, BakMriJsonFile);
+ if ~isOk, warning(Msg); end
+ if ~exist(BakMriJsonFile, 'file')
+ warning('Unable to back up anatomy BIDS json file. Skipping subject %s.', sSubjects.Subject(iSub).Name);
+ continue;
+ end
+ end
+ end
+ end
+ % We ignore other fiducials after the 3 we use for coregistration.
+ % These were likely not placed well, only automatically by initial
+ % linear template alignment.
+ BstFids = {'NAS', 'LPA', 'RPA'}; %, 'AC', 'PC', 'IH'};
+ % We need to go to original Nifti voxel coordinates, but Brainstorm may have
+ % flipped/permuted dimensions to bring voxels to RAS orientation. If it did, it modified
+ % all sMRI fields accordingly, including under .Header, and it saved the transformation
+ % under .InitTransf 'reorient'.
+ iTransf = find(strcmpi(sMri.InitTransf(:,1), 'reorient'));
+ if ~isempty(iTransf)
+ tReorient = sMri.InitTransf{iTransf(1),2}; % Voxel 0-based transformation, from original to Brainstorm
+ tReorientInv = inv(tReorient);
+ tReorientInv(4,:) = [];
+ end
+
+ if isfield(sMriJson, 'AnatomicalLandmarkCoordinates')
+ isPrevJsonLandmarks = true;
+ % Keep copy to check if anything changed and we should save.
+ PrevJsonLandmarks = sMriJson.AnatomicalLandmarkCoordinates;
+ % Remove previous landmarks regardless. The existing backup in derivatives may or
+ % may not have those - only one original backup is kept - but no practical use in
+ % keeping multiple previous aligments.
+ % Set empty instead of removing field, so that it keeps its order. Could use
+ % orderfields, but more complicated when list of fields differ (set diff first but need all fields).
+ % sMriJson = rmfield(sMriJson, 'AnatomicalLandmarkCoordinates');
+ sMriJson.AnatomicalLandmarkCoordinates = [];
+ else
+ isPrevJsonLandmarks = false;
+ end
+ isLandmarksFound = true;
+ for iFid = 1:numel(BstFids)
+ if iFid < 4
+ CS = 'SCS';
+ else
+ CS = 'NCS';
+ end
+ Fid = BstFids{iFid};
+ % Voxel coordinates (Nifti: 0-indexed, but orientation not standardized, world coords are RAS)
+ % Bst MRI coordinates are in mm and voxels are 1-indexed, so subtract 1 voxel after going from mm to voxels.
+ if isfield(sMri, CS) && isfield(sMri.(CS), Fid) && ~isempty(sMri.(CS).(Fid)) && any(sMri.(CS).(Fid))
+ % Round to 0.001 voxel.
+ % Voxsize has 3 elements, ok for non-isotropic voxel size
+ FidCoord = round(1000 * (sMri.(CS).(Fid)./sMri.Voxsize - 1)) / 1000;
+ if ~isempty(iTransf)
+ % Go from Brainstorm RAS-oriented voxels, back to original Nifti voxel orientation.
+ % Both are 0-indexed in this transform.
+ FidCoord = [FidCoord, 1] * tReorientInv';
+ end
+ sMriJson.AnatomicalLandmarkCoordinates.(Fid) = FidCoord;
+ else
+ isLandmarksFound = false;
+ break;
+ end
+ end
+ % Only save if something changed.
+ isSaveMri = false;
+ if ~isLandmarksFound
+ if isPrevJsonLandmarks
+ if isOverwrite
+ warning('MRI landmark coordinates not found, but previously saved in T1w.json file. Removing field and skipping subject %s.', sSubjects.Subject(iSub).Name);
+ isSaveMri = true;
+ else
+ warning('MRI landmark coordinates not found, but previously saved in T1w.json file. Skipping subject %s.', sSubjects.Subject(iSub).Name);
+ end
+ else
+ warning('MRI landmark coordinates not found. Skipping subject %s.', sSubjects.Subject(iSub).Name);
+ end
+ % In case some were found but not all, just remove them all again.
+ sMriJson.AnatomicalLandmarkCoordinates = [];
+ elseif isPrevJsonLandmarks
+ % Verify if different and replacing is needed, otherwise no warning and continue to MEG.
+ % Vector orientation may differ because of json formatting, check values only.
+ if ~isequal(PrevJsonLandmarks.NAS(:), sMriJson.AnatomicalLandmarkCoordinates.NAS(:)) || ...
+ ~isequal(PrevJsonLandmarks.LPA(:), sMriJson.AnatomicalLandmarkCoordinates.LPA(:)) || ...
+ ~isequal(PrevJsonLandmarks.RPA(:), sMriJson.AnatomicalLandmarkCoordinates.RPA(:))
+ if isOverwrite
+ fprintf('Replacing previous MRI landmark coordinates for subject %s.\n', sSubjects.Subject(iSub).Name);
+ isSaveMri = true;
+ else
+ warning('Previous MRI landmark coordinates do not match current ones, but asked to not overwrite, so skipping subject %s.\n', sSubjects.Subject(iSub).Name);
+ isLandmarksFound = false;
+ end
+ end
+ end
+ if isSaveMri
+ % Remove field if empty. There are no other anat landmark fields in MRI json (as opposed to MEG, see below).
+ if isempty(sMriJson.AnatomicalLandmarkCoordinates)
+ sMriJson = rmfield(sMriJson, 'AnatomicalLandmarkCoordinates');
+ end
+ if isDryRun
+ JsonText = bst_jsonencode(sMriJson, true); % indent -> force bst
+ disp(MriJsonFile);
+ disp(JsonText);
+ else
+ WriteJson(MriJsonFile, sMriJson);
+ end
+ end
+ OutFilesMri{iOutSub} = MriJsonFile;
+ if ~isLandmarksFound
+ continue;
+ end
+
+
+ % --------------------------------------------------------------------------------------
+ % MEG _coordsystem.json
+ % Save MRI anatomical landmarks in SCS coordinates and link to MRI.
+ % This includes coregistration refinement using head points, if used.
+
+ % Convert from mri to scs.
+ for iFid = 1:3
+ Fid = BstFids{iFid};
+ % cs_convert mri is in meters
+ sMriScs.SCS.(Fid) = cs_convert(sMri, 'mri', 'scs', sMri.SCS.(Fid) ./ 1000);
+ end
+ sMriNative = sMriScs; % transformed below
+
+ % MEG _coordsystem.json are shared between studies (recordings) within a session.
+ % Get (session, study, recording) for each study requested, in the form of
+ % coordsystem.json file, channel file, and first linked raw data file.
+ MegList = table('Size', [0, 3], 'VariableNames', {'CoordJson', 'Channel', 'Recording'}, 'VariableTypes', {'char', 'char', 'char'});
+ for iStudy = 1:numel(sStudiesForSub)
+ % Try to find the first link to raw file in this study.
+ isLinkToRaw = false;
+ for iData = 1:numel(sStudiesForSub(iStudy).Data)
+ if strcmpi(sStudiesForSub(iStudy).Data(iData).DataType, 'raw')
+ isLinkToRaw = true;
+ break;
+ end
+ end
+ % Skip study if no raw file found.
+ if ~isLinkToRaw
+ continue;
+ end
+ Recording = load(file_fullpath(sStudiesForSub(iStudy).Data(iData).FileName));
+ Recording = Recording.F.filename;
+ % Skip empty-room noise recordings
+ if contains(Recording, 'sub-emptyroom') || contains(Recording, 'task-noise')
+ % No warning, just skip.
+ continue;
+ end
+ % Skip if original MEG file not found.
+ if ~exist(Recording, 'file')
+ warning('Missing original raw MEG file. Skipping study %s.', Recording);
+ continue;
+ end
+
+ % Find MEG _coordsystem.json, should be 1 per session.
+ if strcmpi(Recording(end-4:end), '.meg4')
+ Recording = fileparts(Recording);
+ end
+ MegPath = fileparts(Recording);
+ MegCoordJsonFile = dir(fullfile(MegPath, '*_coordsystem.json'));
+ % Skip if MEG json file not found, or more than one.
+ if isempty(MegCoordJsonFile)
+ warning('MEG BIDS _coordsystem.json file not found. Skipping study %s.', Recording);
+ continue;
+ elseif numel(MegCoordJsonFile) > 1
+ warning('MEG BIDS issue: found multiple _coordsystem.json files. Skipping study %s.', Recording);
+ continue;
+ end
+ % Add to be processed.
+ MegList(end+1, :) = {fullfile(MegCoordJsonFile.folder, MegCoordJsonFile.name), ...
+ sStudiesForSub(iStudy).Channel.FileName, Recording}; %#ok
+ end
+
+ % Keep only one row per session. Unique by default returns the index of the first
+ % occurrence prior to any sorting, i.e. as requested.
+ % TODO: would bst_process sort inputs before getting to the process function?
+ [~, iKeepFirsts] = unique(MegList.CoordJson);
+ % Warn if there were duplicates in the requested studies, but not if we're processing
+ % the entire protocol (no requested studies).
+ if ~isempty(iStudies) && numel(iKeepFirsts) < size(MegList, 1)
+ if isInteractive
+ % Request confirmation.
+ [Proceed, isCancel] = java_dialog('confirm', [...
+ 'Coregistration in BIDS applies globally in a recording session, while in ' 10 ...
+ 'Brainstorm there is no session-level organization and different studies ' 10 ...
+ 'can be coregistered independently. Multiple studies were provided ' 10 ...
+ 'corresponding to the same BIDS session (_coordsystem.json file), but only ' 10 ...
+ 'the alignment of the first (in the order provided, excluding emtpy-room ' 10 ...
+ 'noise) will be saved. ' 10 10 ...
+ 'Proceed with the first study for each session?' 10], 'Save coregistration');
+ if ~Proceed || isCancel
+ % isCancel = true;
+ return;
+ end
+ else
+ % Do not proceed.
+ Message = [...
+ 'Coregistration in BIDS applies globally in a recording session, while in ' ...
+ 'Brainstorm there is no session-level organization and different studies ' ...
+ 'can be coregistered independently. Multiple studies were provided ' ...
+ 'corresponding to the same BIDS session (_coordsystem.json file). Either ' ...
+ 'run the save coregistration process in interactive mode to confirm only ' ...
+ 'the first study''s alignment per session should be saved, or only provide ' ...
+ 'a single study per session.']; %#ok
+ %isCancel = true;
+ warning(Message);
+ return;
+ end
+ end
+ % Sort and unique table rows
+ MegList = MegList(iKeepFirsts, :);
+ nChan = size(MegList, 1);
+ %ChanNativeTransf = zeros(3, 4);
+ iOutMeg = 0;
+ for iChan = 1:nChan
+ % If same json as previous, move on.
+ % Same if same channel file, but not expected to happen for CTF.
+ if iChan > 1 && ( strcmp(MegList.CoordJson{iChan}, MegList.CoordJson{iChan-1}) || ...
+ strcmp(MegList.Channel{iChan}, MegList.Channel{iChan-1}) )
+ continue;
+ end
+ ChannelMat = in_bst_channel(MegList.Channel{iChan});
+ % ChannelMat.SCS are *digitized* anatomical landmarks (if present, otherwise might be
+ % digitized head coils) in Brainstorm/SCS coordinates (defined as CTF but with
+ % anatomical landmarks). They are NOT updated after refining with head points, so we
+ % don't rely on them but use those saved in sMri, and update them now with
+ % UpdateChannelMatScs.
+ %
+ % We applied MRI=>SCS (from sMri) to the MRI anat landmarks above, and now need to apply
+ % SCS=>Native (from ChannelMat). We ignore head motion related adjustments, which are
+ % dataset specific. We need original raw Native coordinates. UpdateChannelMatScs also
+ % adds a "Native" copy of .SCS, which represents the digitized anatomical fiducials in
+ % Native coordinates. Here we just use the transformation, as we want the MRI anat
+ % fids, not the digitized ones (which can be different if another session).
+ ChannelMat = process_adjust_coordinates('UpdateChannelMatScs', ChannelMat);
+ % If no digitization, native transformation will be missing. Could warn and use
+ % identity, assuming MRI fids actually match the head coils, not anatomical points.
+ % But we would still need to modify the descriptions to indicate the points are
+ % coils, which can only be true if there was no previous session with points - and
+ % we should therefore check that. For now, just warn and move on.
+ % TODO: button to display initial MEG coils as fid points on coreg edit figure.
+ if ~isfield(ChannelMat, 'Native')
+ %warning('Could not get native transformation, which possibly indicates missing digitized head points. Assuming that MRI fiducials are then actually head coils and not anatomical points, as is customary without digitized points. %s', MegList.Channel{iChan});
+ %ChanNativeTransf = [eye(3) zeros(3,1)];
+ warning('Could not get native transformation, which possibly indicates missing digitized head points. For single session, coregistration could be saved with head coils on MRI, but this is not yet implemented here. %s', MegList.Channel{iChan});
+ continue;
+ end
+
+ % TODO, TEMPORARY HACK: we've only coregistered the first session for each subject so far.
+ % Skip other sessions.
+ % if iChan > 1 % implies not first session here
+ % % Skipping other sessions.
+ % break;
+ % end
+
+ % New json, store SCS>Native transformation to compare with next channel files in this session.
+ ChanNativeTransf = [ChannelMat.Native.R, ChannelMat.Native.T];
+
+ % Convert MRI fids from (possibly adjusted) SCS to Native, and m to cm.
+ for iFid = 1:3
+ Fid = BstFids{iFid};
+ sMriNative.SCS.(Fid)(:) = 100 * ChanNativeTransf * [sMriScs.SCS.(Fid)'; 1];
+ % Round to um.
+ sMriNative.SCS.(Fid) = round(sMriNative.SCS.(Fid) * 10000) / 10000;
+ end
+
+ % Check MRI-digitized anat fids match, and prepare values.
+ % The description here will be appended to what already describes the coils and anat
+ % fids, which should explain cases when one or the other are missing. The wording
+ % added here is valid for these cases, though this could be slightly confusing (yet
+ % should be inconsequential).
+ [~, isMriUpdated, isMriMatch, isSessionMatch] = process_adjust_coordinates('CheckPrevAdjustments', ChannelMat, sMri);
+ if ~isMriUpdated
+ % For now we skip if the mri was not updated, as this is our workflow for coreg.
+ warning('MRI landmarks have not been updated. This is unexpected in our current workflow and should be verified. Study %s.', Recording);
+ %continue; % Don't skip, we had one participant where the initial fids were good!
+ end
+ if ~isMriMatch && isSessionMatch
+ % The sMri and digitized anat fids match in terms of shape, but they are not
+ % aligned. Probably some other alignment was performed after updating the sMri.
+ % This is unexpected and should be checked.
+ warning('MRI and digitized anat fids are from the same session, but not aligned. This is unexpected and should be verified. Skipping study %s.', Recording);
+ continue;
+ end
+ IntendedForMri = strrep(ImportedFile, [BidsRoot filesep], 'bids::');
+
+ % Make backup in derivatives folder. May not be very useful since it can be
+ % recreated easily except for coreg data which is not yet saved.
+ if isBackup
+ BakMegJsonFile = replace(MegList.CoordJson{iChan}, BidsRoot, fullfile(BidsRoot, 'derivatives'));
+ if ~exist(BakMegJsonFile, 'file')
+ BakFolder = fileparts(BakMegJsonFile);
+ if ~exist(BakFolder, 'dir')
+ [isOk, Msg] = mkdir(BakFolder);
+ if ~isOk, warning(Msg); end
+ if ~exist(BakFolder, 'dir')
+ warning('Unable to create backup folder %s. Skipping session.', BakFolder);
+ continue;
+ end
+ end
+ [isOk, Msg] = copyfile(MegList.CoordJson{iChan}, BakMegJsonFile);
+ if ~isOk, warning(Msg); end
+ if ~exist(BakMegJsonFile, 'file')
+ warning('Unable to back up MEG coordinates BIDS json file. Skipping session %s.', MegList.CoordJson{iChan});
+ continue;
+ end
+ end
+ end
+ % Load existing json
+ sMegJson = bst_jsondecode(MegList.CoordJson{iChan});
+ % Update MEG json file.
+ % Possibly fully recreate the content, e.g. if descriptions were updated.
+ if RecreateMegCoordJson
+ % EEG may not be present in all runs, so check if it was there
+ if isfield(sMegJson, 'EEGCoordinateSystem')
+ isEeg = true;
+ else
+ isEeg = false;
+ end
+ [~, BidsInfo] = BidsParseRecordingName(MegList.Recording{iChan}, [], false);
+ sMegJson = BidsBuildRecordingFiles(MegList.Recording{iChan}, BidsInfo, [], false); % Don't save, will just output structure.
+ % Keep only content of _coordsystem.json
+ sMegJson = sMegJson.CoordSystem;
+ if ~isfield(sMegJson, 'EEGCoordinateSystem') && isEeg
+ % Double check digitized points are present, but should be the case.
+ if ~isfield(sMegJson, 'DigitizedHeadPointsCoordinateSystem')
+ warning('EEG were present, but no head points when trying to recreate %s', MegList.CoordJson{iChan});
+ else
+ % Add back
+ sMegJson.EEGCoordinateSystem = sMegJson.DigitizedHeadPointsCoordinateSystem;
+ sMegJson.EEGCoordinateUnits = sMegJson.DigitizedHeadPointsCoordinateUnits;
+ sMegJson.EEGCoordinateSystemDescription = sMegJson.DigitizedHeadPointsCoordinateSystemDescription;
+ end
+ end
+ end
+ % Check if the coreg was done with anat or coils, before adding or modifying anat.
+ isAnatFids = isfield(sMegJson, 'AnatomicalLandmarkCoordinates');
+ if isMriMatch % implies session matches
+ if isAnatFids
+ LandmarkDescrip = 'They correspond to digitized landmarks from this session. ';
+ else
+ LandmarkDescrip = 'They correspond to digitized head coils from this session, as no anatomical landmarks were digitized. They are still named NAS/LPA/RPA for consistency throughout this dataset, therefore simplifying importing coregistration info. ';
+ end
+ elseif ~isSessionMatch
+ % We still use the sMri fids, whether they were updated (different session) or not.
+ LandmarkDescrip = 'They correspond to a set of digitized landmarks from another session, used for coregistration for this subject. They are usually anatomical points, but sometimes head coils in older data. ';
+ end
+ AddFidDescrip = [' The anatomical landmarks saved here match those in the associated T1w image (see IntendedFor field). ', ...
+ LandmarkDescrip, ...
+ 'Coregistration with the T1w image was performed with Brainstorm before defacing, ', ...
+ 'initially with an automatic procedure fitting head points to the scalp surface, but ', ...
+ 'often adjusted manually, and validated with pictures of MEG head coils on the participant when available. ', ...
+ 'As such, these landmarks and the corresponding alignment should be preferred.'];
+ % Remove previous coordinates, though for MEG it may never get saved unless it's filled below.
+ % Still good practice in case it was done wrong previously (e.g. bad field names)
+ % Set empty instead of removing field, so that it keeps its order, however we use
+ % orderfields later.
+ sMegJson.AnatomicalLandmarkCoordinates = [];
+ % Here we want to only point to the aligned MRI, even if there are multiple MRIs in
+ % this BIDS subject and they were all listed previously. But inform about any change.
+ if isfield(sMegJson, 'IntendedFor') && ~isempty(sMegJson.IntendedFor)
+ if iscell(sMegJson.IntendedFor)
+ if numel(sMegJson.IntendedFor)
+ fprintf('Replaced "IntendedFor" MEG json field (had multiple). %s\n', MegList.CoordJson{iChan});
+ sMegJson.IntendedFor = IntendedForMri; % to simplify following checks, but we do it again below for every case.
+ else % single cell; we don't expect this case
+ sMegJson.IntendedFor = sMegJson.IntendedFor{1};
+ end
+ end
+ % Check if it's different and not just the new BIDS URI convention.
+ if ~strcmpi(sMegJson.IntendedFor, IntendedForMri) && ...
+ ~contains(sMegJson.IntendedFor, strrep(IntendedForMri, 'bids::', ''))
+ fprintf('Replaced "IntendedFor" MEG json field: %s > %s in %s\n', sMegJson.IntendedFor, IntendedForMri, MegList.CoordJson{iChan});
+ end
+ end
+ % Save the single aligned MRI.
+ sMegJson.IntendedFor = IntendedForMri;
+ % Save native coordinates (rounded to um above).
+ for iFid = 1:3
+ Fid = BstFids{iFid};
+ sMegJson.AnatomicalLandmarkCoordinates.(Fid) = sMriNative.SCS.(Fid);
+ end
+ if ~isfield(sMegJson, 'HeadCoilCoordinateSystem')
+ warning('No head coils when trying to copy for anat landmark system fields %s', MegList.CoordJson{iChan});
+ else
+ sMegJson.AnatomicalLandmarkCoordinateSystem = sMegJson.HeadCoilCoordinateSystem;
+ sMegJson.AnatomicalLandmarkCoordinateUnits = sMegJson.HeadCoilCoordinateUnits;
+ sMegJson.AnatomicalLandmarkCoordinateSystemDescription = sMegJson.HeadCoilCoordinateSystemDescription;
+ end
+ % System descriptions are saved in BidsBuildRecordingFiles.
+
+ if ~isfield(sMegJson, 'FiducialsDescription')
+ sMegJson.FiducialsDescription = '';
+ end
+ % Append description, if not already there.
+ if isempty(strfind(sMegJson.FiducialsDescription, AddFidDescrip))
+ % Remove possibly inaccurate description of anat fids from other session.
+ if ~isSessionMatch
+ % This regexp matches the given words, followed by max number of non-period characters, then a period (escaped) '\.', and 0 or 1 space ' ?'
+ sMegJson.FiducialsDescription = regexprep(sMegJson.FiducialsDescription, 'The anatomical landmarks[^.]*\. ?', '');
+ end
+ sMegJson.FiducialsDescription = strtrim([sMegJson.FiducialsDescription, AddFidDescrip]);
+ end
+ % Remove fields if anat coords empty, though should not really happen.
+ if isempty(sMegJson.AnatomicalLandmarkCoordinates)
+ AnatLandmarkFields = {'AnatomicalLandmarkCoordinates', 'AnatomicalLandmarkCoordinateSystem', 'AnatomicalLandmarkCoordinateUnits', 'AnatomicalLandmarkCoordinateSystemDescription'};
+ isFieldFound = ismember(AnatLandmarkFields, fieldnames(sMriJson));
+ sMriJson = rmfield(sMriJson, AnatLandmarkFields(isFieldFound));
+ warning('No anat landmark coordinates in %s\n', MegList.CoordJson{iChan});
+ end
+ % Reorder fields according to this full list of all possible fields.
+ CoordJsonFields = {'MEGCoordinateSystem', 'MEGCoordinateUnits', 'MEGCoordinateSystemDescription', ...
+ 'EEGCoordinateSystem', 'EEGCoordinateUnits', 'EEGCoordinateSystemDescription', ...
+ 'DigitizedHeadPoints', 'DigitizedHeadPointsCoordinateSystem', 'DigitizedHeadPointsCoordinateUnits', 'DigitizedHeadPointsCoordinateSystemDescription', ...
+ 'HeadCoilCoordinates', 'HeadCoilCoordinateSystem', 'HeadCoilCoordinateUnits', 'HeadCoilCoordinateSystemDescription', ...
+ 'AnatomicalLandmarkCoordinates', 'AnatomicalLandmarkCoordinateSystem', 'AnatomicalLandmarkCoordinateUnits', 'AnatomicalLandmarkCoordinateSystemDescription', ...
+ 'FiducialsDescription', 'IntendedFor'};
+ % But we must give a matching list, no extras.
+ % Warn if there are additional fields not accounted for here, and don't sort.
+ % (Should not happen in our current workflow.)
+ if any(~ismember(fieldnames(sMegJson), CoordJsonFields))
+ warning('Unexpected fields in coordsystem.json stucture %s', MegList.CoordJson{iChan});
+ else
+ sMegJson = orderfields(sMegJson, CoordJsonFields(ismember(CoordJsonFields, fieldnames(sMegJson))));
+ end
+ % Save
+ if isDryRun
+ JsonText = bst_jsonencode(sMegJson, true); % indent -> force bst
+ disp(MegList.CoordJson{iChan});
+ disp(JsonText);
+ else
+ WriteJson(MegList.CoordJson{iChan}, sMegJson);
+ end
+ iOutMeg = iOutMeg + 1;
+ OutFilesMeg{iOutSub}{iOutMeg,1} = MegList.CoordJson{iChan};
+ end % channel file (studies/sessions) loop
+ else
+ % Not BIDS, save in fiducials.m file.
+ FidsFile = fullfile(bst_fileparts(ImportedFile), 'fiducials.m');
+ FidsFile = figure_mri('SaveFiducialsFile', sMri, FidsFile);
+ if ~exist(FidsFile, 'file')
+ warning('Fiducials.m file not written for subject %s.', sSubjects.Subject(iSub).Name);
+ continue;
+ end
+ OutFilesMri{iOutSub} = FidsFile;
+ end
+
+ isSuccess(iOutSub) = true;
+ bst_progress('inc', 1);
+ end % subject loop
+
+ bst_progress('stop');
+
+end
+
+
diff --git a/toolbox/math/bst_meshfit.m b/toolbox/math/bst_meshfit.m
index 5ad40a2562..e4ccc2baa9 100644
--- a/toolbox/math/bst_meshfit.m
+++ b/toolbox/math/bst_meshfit.m
@@ -1,18 +1,19 @@
-function [R, T, newP, distFinal] = bst_meshfit(Vertices, Faces, P)
+function [R, T, newP, distFinal] = bst_meshfit(Vertices, Faces, P, Outliers)
% BST_MESHFIT: Find the best possible rotation-translation to fit a point cloud on a mesh.
%
% USAGE: [R, T, newP, distFinal] = bst_meshfit(Vertices, Faces, P)
%
-% DESCRIPTION:
+% DESCRIPTION:
% A Gauss-Newton method is used for the optimization of the distance points/mesh.
-% The Gauss-Newton algorithm used here was initially implemented by
+% The Gauss-Newton algorithm used here was initially implemented by
% Qianqian Fang (fangq at nmr.mgh.harvard.edu) and distributed under a GPL license
-% as part of the Metch toolbox (http://iso2mesh.sf.net, regpt2surf.m).
+% as part of the Metch toolbox (http://iso2mesh.sf.net, regpt2m).
%
% INPUTS:
% - Vertices : [Mx3] double matrix
% - Faces : [Nx3] double matrix
% - P : [Qx3] double matrix, points to fit on the mesh defined by Vertices/Faces
+% - Outliers : proportion of outlier points to ignore (between 0 and 1)
%
% OUTPUTS:
% - R : [3x3] rotation matrix from the original P to the fitted positions.
@@ -23,12 +24,12 @@
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
-%
+%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
-%
+%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
@@ -40,20 +41,65 @@
%
% Authors: Qianqian Fang, 2008
% Francois Tadel, 2013-2021
+% Marc Lalancette, 2022
+
+% Coordinates are in m.
+PenalizeInside = true;
+SquareDistCost = false;
+if nargin < 4 || isempty(Outliers)
+ Outliers = 0;
+end
+
+% nV = size(Vertices, 1);
+nF = size(Faces, 1);
+nP = size(P, 1);
+Outliers = ceil(Outliers * nP);
+
+% Edges as indices
+Edges = unique(sort([Faces(:,[1,2]); Faces(:,[2,3]); Faces(:,[3,1])], 2), 'rows');
+% Edge direction "doubly normalized" so that later projection should be between 0 and 1.
+EdgeDir = Vertices(Edges(:,2),:) - Vertices(Edges(:,1),:);
+EdgeL = sqrt(sum(EdgeDir.^2, 2));
+EdgeDir = bsxfun(@rdivide, EdgeDir, EdgeL);
+% Edges as vectors
+EdgesV = zeros(nF, 3, 3);
+EdgesV(:,:,1) = Vertices(Faces(:,2),:) - Vertices(Faces(:,1),:);
+EdgesV(:,:,2) = Vertices(Faces(:,3),:) - Vertices(Faces(:,2),:);
+EdgesV(:,:,3) = Vertices(Faces(:,1),:) - Vertices(Faces(:,3),:);
+% First edge to second edge: counter clockwise = up
+FaceNormals = cross(EdgesV(:,:,1), EdgesV(:,:,2));
+FaceNormals = bsxfun(@rdivide, FaceNormals, sqrt(sum(FaceNormals.^2, 2))); % normr
+%FaceArea = sqrt(sum(FaceNormals.^2, 2));
+% Perpendicular vectors to edges, pointing inside triangular face.
+for e = 3:-1:1
+ EdgeTriNormals(:,:,e) = cross(FaceNormals, EdgesV(:,:,e));
+end
+FaceVertices = zeros(nF, 3, 3);
+FaceVertices(:,:,1) = Vertices(Faces(:,1),:);
+FaceVertices(:,:,2) = Vertices(Faces(:,2),:);
+FaceVertices(:,:,3) = Vertices(Faces(:,3),:);
-% Calculate norms
-VertNorm = tess_normals(Vertices, Faces);
% Calculate the initial error
-[distInit, dt] = get_distance(Vertices, VertNorm, P, []);
-errInit = sum(abs(distInit));
+InitParams = zeros(6,1);
+errInit = CostFunction(InitParams);
% Fit points
-[R,T,newP] = fit_points(Vertices, VertNorm, P, dt);
+% [R,T,newP] = fit_points(Vertices, VertNorm, P, dt);
+% Do optimization
+% Stop at 0.1 mm total distance, or 0.02 mm displacement.
+OptimOptions = optimoptions(@fminunc, 'MaxFunctionEvaluations', 1000, 'MaxIterations', 200, ...
+ 'FiniteDifferenceStepSize', 1e-3, ...
+ 'FunctionTolerance', 1e-4, 'StepTolerance', 2e-2, 'Display', 'none'); % 'OptimalityTolerance', 1e-15, 'final-detailed'
+
+BestParams = fminunc(@CostFunction, InitParams, OptimOptions);
+[R,T,newP] = Transform(BestParams, P);
+T = T';
+distFinal = PointSurfDistance(newP);
% Calculate the final error
-distFinal = get_distance(Vertices, VertNorm, newP, dt);
-errFinal = sum(abs(distFinal));
+errFinal = CostFunction(BestParams);
% If the error is larger than at the beginning: cancel the modifications
+% Should no longer occur.
if (errFinal > errInit)
disp('BST> The optimization failed finding a better fit');
R = [];
@@ -61,88 +107,111 @@
newP = P;
end
-end
-
+% Better cost function for points fitting: higher cost for points inside the
+% head > 1mm, (better distance calculation).
+ function [Cost, Dist] = CostFunction(Params)
+ [~,~,Points] = Transform(Params, P);
+ Dist = PointSurfDistance(Points);
+ if PenalizeInside
+ isInside = inpolyhedron(Faces, Vertices, Points, 'FaceNormals', FaceNormals, 'FlipNormals', true);
+ %patch('Faces',Faces,'Vertices',Vertices); hold on; light; axis equal;
+ %quiver3(FaceVertices(:,1,1), FaceVertices(:,2,1), FaceVertices(:,3,1), FaceNormals(:,1,1), FaceNormals(:,2,1), FaceNormals(:,3,1));
+ %scatter3(Points(1,1),Points(1,2),Points(1,3));
+ iSquare = isInside & Dist > 0.001;
+ Dist(iSquare) = Dist(iSquare).^2 *1e3; % factor for "squaring in mm units"
+ iOutside = find(~isInside);
+ end
+ if SquareDistCost
+ Cost = sum(Dist.^2);
+ else
+ Cost = sum(Dist);
+ end
+ for iP = 1:Outliers
+ if PenalizeInside
+ % Only remove outside points.
+ [MaxD, iMaxD] = max(Dist(~isInside));
+ Dist(iOutside(iMaxD)) = 0;
+ else
+ [MaxD, iMaxD] = max(Dist);
+ Dist(iMaxD) = 0;
+ end
+ if SquareDistCost
+ Cost = Cost - MaxD.^2;
+ else
+ Cost = Cost - MaxD;
+ end
+ end
+ end
+%% TODO Slow, look for alternatives.
+% This seems similar: https://www.mathworks.com/matlabcentral/fileexchange/52882-point2trimesh-distance-between-point-and-triangulated-surface
% ===== COMPUTE POINTS/MESH DISTANCE =====
-% Approximates the distance to the mesh by the projection on the norm vector of the nearest neighbor
-function [dist,dt] = get_distance(Vertices, VertNorm, P, dt)
- % Find the nearest neighbor
- [iNearest, dist_pt, dt] = bst_nearest(Vertices, P, 1, 0, dt);
- % Distance = projection of the distance between the point and its nearest
- % neighbor in the surface on the vertex normal
- % As the head surface is supposed to be very smooth, it should be a good approximation
- % of the distance from the point to the surface.
- dist = abs(sum(VertNorm(iNearest,:) .* (P - Vertices(iNearest,:)),2));
-end
-
-% ===== COMPUTE TRANSFORMATION =====
-% Gauss-Newton optimization algorithm
-% Based on work from Qianqian Fang, 2008
-function [R,T,newP] = fit_points(Vertices, VertNorm, P, dt)
- % Initial parameters: no rotation, no translation
- C = zeros(6,1);
- newP = P;
- % Sensitivity
- delta = 1e-4;
- % Maximum number of iterations
- maxiter = 20;
- % Initialize error at the previous iteration
- errPrev = Inf;
- % Start Gauss-Newton iterations
- for iter = 1:maxiter
- % Calculate the current residual: the sum of distances to the surface
- dist0 = get_distance(Vertices, VertNorm, newP, dt);
- err = sum(abs(dist0));
- % If global error is going up: stop
- if (err > errPrev)
- break;
+% For exact distance computation, we need to check all 3: vertices, edges and faces.
+ function Dist = PointSurfDistance(Points)
+ Epsilon = 1e-9; % nanometer
+ % Check distance to vertices
+ if license('test','statistics_toolbox')
+ DistVert = pdist2(Vertices, Points, 'euclidean', 'Smallest', 1)';
+ else
+ DistVert = zeros(nP, 1);
+ for iP = 1:nP
+ % Find closest surface vertex.
+ DistVert(iP) = sqrt(min(sum(bsxfun(@minus, Points(iP, :), Vertices).^2, 2)));
+ end
end
- errPrev = err;
- % fprintf('iter=%d error=%f\n', iter, err);
- % Build the Jacobian (sensitivity) matrix
- J = zeros(length(dist0),length(C));
- for i = 1:length(C)
- dC = C;
- if (C(i))
- dC(i) = C(i) * (1+delta);
- else
- dC(i) = C(i) + delta;
+ % Check distance to faces
+ DistFace = inf(nP, 1);
+ for iP = 1:nP
+ % Considered Möller and Trumbore 1997, Ray-Triangle Intersection (https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d), but this is simpler still.
+ % Vectors from triangle vertices to point.
+ Pyramid = bsxfun(@minus, Points(iP, :), FaceVertices);
+ % Does the point project inside each face?
+ InFace = all(sum(Pyramid .* EdgeTriNormals, 2) > -Epsilon, 3);
+ if any(InFace)
+ DistFace(iP) = min(abs(sum(Pyramid(InFace,:,1) .* FaceNormals(InFace,:), 2)));
end
- % Apply this new transformation to the points
- [tmpR,tmpT,tmpP] = get_transform(dC, P);
- % Calculate the distance for this new transformation
- dist = get_distance(Vertices, VertNorm, tmpP, dt);
- % J=dL/dC
- J(:,i) = (dist-dist0) / (dC(i)-C(i));
end
- % Weight the matrix (normalization)
- wj = sqrt(sum(J.*J));
- J = J ./ repmat(wj,length(dist0),1);
- % Calculate the update: J*dC=dL
- dC = (J\dist0) ./ wj';
- C = C - 0.5*dC;
- % Get the updated positions with the calculated A and b
- [R,T,newP] = get_transform(C, P);
+ % Check distance to edges
+ DistEdge = inf(nP, 1);
+ for iP = 1:nP
+ % Vector from first edge vertex to point.
+ Pyramid = bsxfun(@minus, Points(iP, :), Vertices(Edges(:, 1), :));
+ Projection = sum(Pyramid .* EdgeDir, 2);
+ InEdge = Projection > -Epsilon & Projection < (EdgeL + Epsilon);
+ if any(InEdge)
+ DistEdge(iP) = sqrt(min(sum((Pyramid(InEdge,:) - bsxfun(@times, Projection(InEdge), EdgeDir(InEdge,:))).^2, 2)));
+ end
+ end
+
+ Dist = min([DistVert, DistEdge, DistFace], [], 2);
end
-end
+
+
+% % Approximates the distance to the mesh by the projection on the norm vector of the nearest neighbor
+% function [dist,dt] = get_distance(Vertices, VertNorm, P, dt)
+% % Find the nearest neighbor
+% [iNearest, dist_pt, dt] = bst_nearest(Vertices, P, 1, 0, dt);
+% % Distance = projection of the distance between the point and its nearest
+% % neighbor in the surface on the vertex normal
+% % As the head surface is supposed to be very smooth, it should be a good approximation
+% % of the distance from the point to the surface.
+% dist = abs(sum(VertNorm(iNearest,:) .* (P - Vertices(iNearest,:)),2));
+% end
% ===== GET TRANSFORMATION =====
-function [R,T,P] = get_transform(params, P)
- % Get values
- mx = params(1); my = params(2); mz = params(3); % Translation parameters
- x = params(4); y = params(5); z = params(6); % Rotation parameters
- % Rotation
- Rx = [1 0 0 ; 0 cos(x) sin(x) ; 0 -sin(x) cos(x)]; % Rotation over x
- Ry = [cos(y) 0 -sin(y); 0 1 0; sin(y) 0 cos(y)]; % Rotation over y
- Rz = [cos(z) sin(z) 0; -sin(z) cos(z) 0; 0 0 1]; % Rotation over z
- R = Rx*Ry*Rz;
- % Translation
- T = [mx; my; mz];
- % Apply to points
- if (nargin >= 2)
- P = (R * P' + T * ones(1,size(P,1)))';
+ function [R,T,Points] = Transform(Params, Points)
+ % Translation in mm (to use default TypicalX of 1)
+ T = Params(1:3)'/1e3;
+ % Rotation in degrees (again for expected order of magnitude of 1)
+ x = Params(4)*pi/180; y = Params(5)*pi/180; z = Params(6)*pi/180; % Rotation parameters
+ Rx = [1 0 0 ; 0 cos(x) sin(x) ; 0 -sin(x) cos(x)]; % Rotation over x
+ Ry = [cos(y) 0 -sin(y); 0 1 0; sin(y) 0 cos(y)]; % Rotation over y
+ Rz = [cos(z) sin(z) 0; -sin(z) cos(z) 0; 0 0 1]; % Rotation over z
+ R = Rx*Ry*Rz;
+ % Apply to points
+ Points = bsxfun(@plus, Points * R', T);
end
+
end
diff --git a/toolbox/process/functions/process_adjust_coordinates.m b/toolbox/process/functions/process_adjust_coordinates.m
index bf861a2476..b145c9e44d 100644
--- a/toolbox/process/functions/process_adjust_coordinates.m
+++ b/toolbox/process/functions/process_adjust_coordinates.m
@@ -1,9 +1,10 @@
function varargout = process_adjust_coordinates(varargin)
-% PROCESS_ADJUST_COORDINATES: Adjust, recompute, or remove various coordinate transformations.
+% PROCESS_ADJUST_COORDINATES: Adjust, recompute, or remove various coordinate transformations, primarily for CTF MEG datasets.
%
-% Native coordinates are based on system fiducials (e.g. MEG head coils),
-% whereas Brainstorm's SCS coordinates are based on the anatomical fiducial
-% points from the .pos file.
+% Native coordinates are based on system fiducials (e.g. MEG head coils), whereas Brainstorm's SCS
+% coordinates are based on the anatomical fiducial points. After alignment between MRI and
+% headpoints, the anatomical fiducials on the MRI side define the SCS and the ones in the channel
+% files (ChannelMat.SCS) are ignored.
% @=============================================================================
% This function is part of the Brainstorm software:
@@ -23,14 +24,13 @@
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
-% Authors: Marc Lalancette, 2018-2020
+% Authors: Marc Lalancette, 2018-2022
eval(macro_method);
end
-
-function sProcess = GetDescription() %#ok
+function sProcess = GetDescription()
% Description of the process
sProcess.Comment = 'Adjust coordinate system';
sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/HeadMotion#Adjust_the_reference_head_position';
@@ -42,7 +42,7 @@
sProcess.OutputTypes = {'raw', 'data'};
sProcess.nInputs = 1;
sProcess.nMinFiles = 1;
- % Option [to do: ignore bad segments]
+ % Option
sProcess.options.reset.Type = 'checkbox';
sProcess.options.reset.Comment = 'Reset coordinates using original channel file (removes all adjustments: head, points, manual).';
sProcess.options.reset.Value = 0;
@@ -52,7 +52,7 @@
sProcess.options.format.Type = 'combobox';
sProcess.options.format.Comment = 'For reset option, specify the channel file format:';
sProcess.options.format.Value = {1, FileFormatsChan(:, 2)'};
- sProcess.options.format.Class = 'Reset';
+ % sProcess.options.format.Class = 'Reset';
sProcess.options.head.Type = 'checkbox';
sProcess.options.head.Comment = 'Adjust head position to median location - CTF only.';
sProcess.options.head.Value = 0;
@@ -60,10 +60,18 @@
sProcess.options.bad.Type = 'checkbox';
sProcess.options.bad.Comment = 'For adjust option, exclude bad segments.';
sProcess.options.bad.Value = 1;
- sProcess.options.bad.Class = 'Adjust';
+ sProcess.options.bad.Class = 'Adjust';
sProcess.options.points.Type = 'checkbox';
sProcess.options.points.Comment = 'Refine MRI coregistration using digitized head points.';
sProcess.options.points.Value = 0;
+ sProcess.options.points.Controller = 'Refine';
+ sProcess.options.tolerance.Comment = 'Tolerance (outlier points to ignore):';
+ sProcess.options.tolerance.Type = 'value';
+ sProcess.options.tolerance.Value = {0, '%', 0};
+ sProcess.options.tolerance.Class = 'Refine';
+ sProcess.options.scs.Type = 'checkbox';
+ sProcess.options.scs.Comment = 'Replace MRI nasion and ear points with digitized landmarks (cannot undo).
Requires selecting channel file format above.';
+ sProcess.options.scs.Value = 0;
sProcess.options.remove.Type = 'checkbox';
sProcess.options.remove.Comment = 'Remove selected adjustments (if present) instead of adding them.';
sProcess.options.remove.Value = 0;
@@ -82,7 +90,7 @@
function OutputFiles = Run(sProcess, sInputs)
-
+ OutputFiles = {};
isDisplay = sProcess.options.display.Value;
nInFiles = length(sInputs);
@@ -93,25 +101,58 @@
bst_memory('UnloadAll', 'Forced'); % Close all the existing figures.
end
- [UniqueChan, iUniqFiles, iUniqInputs] = unique({sInputs.ChannelFile});
+ [~, iUniqFiles, iUniqInputs] = unique({sInputs.ChannelFile});
nFiles = numel(iUniqFiles);
- if ~sProcess.options.remove.Value && sProcess.options.head.Value && ...
- nFiles < nInFiles
+ % Special cases when replacing MRI fids: don't allow resetting or refining with head points,
+ % probably a user mistake. Still possible if done in two separate calls. Adjusting head position
+ % ignored with a warning only (would be lost anyway).
+ % Also only a single channel file per subject is allowed.
+ if ~sProcess.options.remove.Value && sProcess.options.scs.Value
+ % TODO: how to go back to process panel, like for missing TF options in some processes?
+ if sProcess.options.reset.Value
+ bst_report('Error', sProcess, sInputs, ...
+ ['Incompatible options: "Reset coordinates" would be applied first and remove coregistration adjustments.' 10, ...
+ '"Replace MRI landmarks" automatically resets all channel files for selected subject(s) after MRI update.']);
+ return;
+ end
+ if sProcess.options.points.Value
+ bst_report('Error', sProcess, sInputs, ...
+ ['Incompatible options: "Refine with head points" not currently allowed in same call as' 10, ...
+ '"Replace MRI landmarks" to avoid user error and possibly loosing manual coregistration adjustments.' 10, ...
+ 'If you really want to do this, run this process separately for each operation.']);
+ return;
+ end
+ if sProcess.options.head.Value
+ bst_report('Warning', sProcess, sInputs, ...
+ ['Incompatible options: "Adjust head position" ignored since it would be lost after MRI update.' 10, ...
+ '"Replace MRI landmarks" automatically resets all channel files for selected subject(s) after MRI update.']);
+ sProcess.options.head.Value = 0;
+ end
+
+ % Check for multiple channel files per subject.
+ UniqueSubs = unique({sInputs(iUniqFiles).SubjectFile});
+ if numel(UniqueSubs) < nFiles
+ bst_report('Error', sProcess, sInputs, ...
+ '"Replace MRI landmarks" can only be run with a single channel file per subject.');
+ return;
+ end
+ end
+
+ if ~sProcess.options.remove.Value && sProcess.options.head.Value && nFiles < nInFiles
bst_report('Info', sProcess, sInputs, ...
'Multiple inputs were found for a single channel file. They will be concatenated for adjusting the head position.');
end
- bst_progress('start', 'Adjust coordinate system', ...
- ' ', 0, nFiles);
- % If resetting, in case the original data moved, and because the same
- % channel file may appear in many places for processed data, keep track
- % of user file selections.
+
+ bst_progress('start', 'Adjust coordinate system', ' ', 0, nFiles);
+ % If resetting, in case the original data moved, and because the same channel file may appear in
+ % many places for processed data, keep track of user file selections.
NewChannelFiles = cell(0, 2);
for iFile = iUniqFiles(:)' % no need to repeat on same channel file.
ChannelMat = in_bst_channel(sInputs(iFile).ChannelFile);
% Get the leading modality
- [tmp, DispMod] = channel_get_modalities(ChannelMat.Channel);
+ [~, DispMod] = channel_get_modalities(ChannelMat.Channel);
% 'MEG' is added when 'MEG GRAD' or 'MEG MAG'.
ModPriority = {'MEG', 'EEG', 'SEEG', 'ECOG', 'NIRS'};
iMod = find(ismember(ModPriority, DispMod), 1, 'first');
@@ -132,52 +173,69 @@
% ----------------------------------------------------------------
if sProcess.options.reset.Value
- % The main goal of this option is to fix a bug in a previous
- % version: when importing a channel file, when going to SCS
- % coordinates based on digitized coils and anatomical
- % fiducials, the channel orientation was wrong. We wish to fix
- % this but keep as much pre-processing that was previously
- % done. Thus we will re-import the channel file, and copy the
- % projectors (and history) from the old one.
+ % The original goal of this option was to fix data affected by a previous bug while
+ % keeping as much pre-processing that was previously done. We re-import the channel
+ % file, and copy the projectors (and history) from the old one.
- [ChannelMat, NewChannelFiles, Failed] = ...
- ResetChannelFile(ChannelMat, NewChannelFiles, sInputs(iFile), sProcess);
- if Failed
+ [ChannelMat, NewChannelFiles, isError] = ResetChannelFile(ChannelMat, ...
+ NewChannelFiles, sInputs(iFile), sProcess);
+ if isError
continue;
end
% ----------------------------------------------------------------
elseif sProcess.options.remove.Value
- % Because channel_align_manual does not consistently apply the
- % manual transformation to all sensors or save it in both
- % TransfMeg and TransfEeg, it could lead to confusion and
- % errors when playing with transforms. Therefore, if we detect
- % a difference between the MEG and EEG transforms when trying
- % to remove one that applies to both (currently only refine
- % with head points), we don't proceed and recommend resetting
- % with the original channel file instead.
+ % Because channel_align_manual does not consistently apply the manual transformation to
+ % all sensors or save it in both TransfMeg and TransfEeg, it could lead to confusion and
+ % errors when playing with transforms. Therefore, if we detect a difference between the
+ % MEG and EEG transforms when trying to remove one that applies to both (currently only
+ % refine with head points), we don't proceed and recommend resetting with the original
+ % channel file instead.
Which = {};
if sProcess.options.head.Value
- Which{end+1} = 'AdjustedNative';
+ Which{end+1} = 'AdjustedNative'; %#ok
end
if sProcess.options.points.Value
- Which{end+1} = 'refine registration: head points';
+ Which{end+1} = 'refine registration: head points'; %#ok
end
for TransfLabel = Which
- TransfLabel = TransfLabel{1};
+ TransfLabel = TransfLabel{1}; %#ok
ChannelMat = RemoveTransformation(ChannelMat, TransfLabel, sInputs(iFile), sProcess);
end % TransfLabel loop
+
+ % We cannot change back the MRI fiducials, but in order to be able to update it again
+ % from digitized fids without warnings, edit the MRI history.
+ if sProcess.options.scs.Value
+ % Get subject in database, with subject directory
+ sSubject = bst_get('Subject', sInputs(iFile).SubjectFile);
+ MriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName;
+ sMri = in_mri_bst(sSubject.Anatomy(sSubject.iAnatomy).FileName);
+ % Slightly change the string we use to verify if it was done: append " (hidden)".
+ iHist = find(strcmpi(sMri.History(:,3), 'Applied digitized anatomical fiducials'));
+ if ~isempty(iHist)
+ for iH = 1:numel(iHist)
+ sMri.History{iHist(iH),3} = [sMri.History{iHist(iH),3}, ' (hidden)'];
+ end
+ try
+ bst_save(file_fullpath(MriFile), sMri, 'v7');
+ catch
+ bst_report('Error', sProcess, sInputs(iFile), ...
+ sprintf('Unable to save MRI file %s.', MriFile));
+ continue;
+ end
+ end
+ end
end % reset channel file or remove transformations
% ----------------------------------------------------------------
if ~sProcess.options.remove.Value && sProcess.options.head.Value
% Complex indexing to get all inputs for this same channel file.
- [ChannelMat, Failed] = AdjustHeadPosition(ChannelMat, ...
+ [ChannelMat, isError] = AdjustHeadPosition(ChannelMat, ...
sInputs(iUniqInputs == iUniqInputs(iFile)), sProcess);
- if Failed
+ if isError
continue;
end
end % adjust head position
@@ -187,23 +245,57 @@
% Redundant, but makes sense to have it here also.
bst_progress('text', 'Fitting head surface to points...');
- [ChannelMat, R, T, isSkip] = ...
- channel_align_auto(sInputs(iFile).ChannelFile, ChannelMat, 0, 0); % No warning or confirmation
- % ChannelFile needed to find subject and scalp surface, but not
- % used otherwise when ChannelMat is provided.
- if isSkip
- bst_report('Error', sProcess, sInputs(iFile), ...
- 'Error trying to refine registration using head points.');
+ % If called externally without a tolerance value, set isWarning true so it asks.
+ if isempty(sProcess.options.tolerance.Value)
+ isWarning = true;
+ Tolerance = 0;
+ else
+ isWarning = false;
+ Tolerance = sProcess.options.tolerance.Value{1} / 100;
+ end
+ [ChannelMat, ~, ~, isSkip, isUserCancel, strReport] = channel_align_auto(sInputs(iFile).ChannelFile, ...
+ ChannelMat, isWarning, 0, Tolerance); % No confirmation
+ % ChannelFile needed to find subject and scalp surface, but not used otherwise when
+ % ChannelMat is provided.
+ if ~isempty(strReport)
+ bst_report('Info', sProcess, sInputs(iFile), strReport);
+ elseif isSkip
+ bst_report('Warning', sProcess, sInputs(iFile), ...
+ 'Refine registration using head points, failed finding a better fit.');
continue;
+ elseif isUserCancel
+ bst_report('Info', sProcess, sInputs(iFile), 'User cancelled registration with head points.');
+ continue
end
end % refine registration with head points
% ----------------------------------------------------------------
- % Save channel file.
+ % Save channel file.
+ % Before potiential MRI update since that function takes ChannelFile, not ChannelMat.
bst_save(file_fullpath(sInputs(iFile).ChannelFile), ChannelMat, 'v7');
- isFileOk(iFile) = true;
+ % ----------------------------------------------------------------
+ if ~sProcess.options.remove.Value && sProcess.options.scs.Value
+ % This is now a subfunction in this process
+ [~, isCancel, isError] = channel_align_scs(sInputs(iFile).ChannelFile, eye(4), ...
+ false, false, sInputs(iFile), sProcess); % interactive warnings but no confirmation
+ % Head surface was unloaded from memory, in case we want to display "after" figure.
+ % TODO Maybe modify like channel_align_auto to take in and return ChannelMat.
+ if isError
+ return;
+ elseif isCancel
+ continue;
+ end
+ % If something happened and channel files were not reset, there was a pop-up error and
+ % we'll just let the user deal with it for now (i.e. try again to reset the channel file).
+
+ % TODO Verify
+ end
+
+ % ----------------------------------------------------------------
+ isFileOk(iFile) = true;
+
if isDisplay && ~isempty(Modality)
% Display "after" results, besides the "before" figure.
hFigAfter = channel_align_manual(sInputs(iFile).ChannelFile, Modality, 0);
@@ -212,10 +304,9 @@
end % file loop
bst_progress('stop');
- % Return the input files that were processed properly. Include those
- % that were removed due to sharing a channel file, where appropriate.
- % The complicated indexing picks the first input of those with the same
- % channel file, i.e. the one that was marked ok.
+ % Return the input files that were processed properly. Include those that were removed due to
+ % sharing a channel file, where appropriate. The complicated indexing picks the first input of
+ % those with the same channel file, i.e. the one that was marked ok.
OutputFiles = {sInputs(isFileOk(iUniqInputs(iUniqFiles))).FileName};
end
@@ -266,21 +357,32 @@
% end
-function [ChannelMat, NewChannelFiles, Failed] = ...
- ResetChannelFile(ChannelMat, NewChannelFiles, sInput, sProcess)
- if nargin < 4
+function [ChannelMat, NewChannelFiles, isError] = ResetChannelFile(ChannelMat, NewChannelFiles, sInput, sProcess)
+ % Reload a channel file, but keep projectors and history. First look for original file from
+ % history, and if it's no longer there, user will be prompted. User selections are noted as
+ % pairs {old, new} in NewChannelFiles for potential reuse (e.g. same original data at multiple
+ % pre-processing steps).
+ % This function does not save the file, but only returns the updated structure.
+ if nargin < 4 || isempty(sProcess)
sProcess = [];
+ isReport = false;
+ else
+ isReport = true;
end
- Failed = false;
+ if nargin < 3
+ sInput = [];
+ end
+ if nargin < 2 || isempty(NewChannelFiles)
+ NewChannelFiles = cell(0,2);
+ end
+ isError = false;
bst_progress('text', 'Importing channel file...');
% Extract original data file from channel file history.
- if any(size(ChannelMat.History) < [1, 3]) || ...
- ~strcmp(ChannelMat.History{1, 2}, 'import')
+ if any(size(ChannelMat.History) < [1, 3]) || ~strcmp(ChannelMat.History{1, 2}, 'import')
NotFound = true;
ChannelFile = '';
else
- ChannelFile = regexp(ChannelMat.History{1, 3}, ...
- '(?<=: )(.*)(?= \()', 'match');
+ ChannelFile = regexp(ChannelMat.History{1, 3}, '(?<=: )(.*)(?= \()', 'match');
if isempty(ChannelFile)
NotFound = true;
else
@@ -298,19 +400,26 @@
if NewFound
ChannelFile = NewChannelFiles{iNew, 2};
NotFound = false;
- bst_report('Info', sProcess, sInput, ...
- sprintf('Using channel file in new location: %s.', ChannelFile));
+ if isReport
+ bst_report('Info', sProcess, sInput, ...
+ sprintf('Using channel file in new location: %s.', ChannelFile));
+ end
end
end
- FileFormatsChan = bst_get('FileFilters', 'channel');
- FileFormat = FileFormatsChan{sProcess.options.format.Value{1}, 3};
+ if isfield(sProcess, 'options') && isfield(sProcess.options, 'format')
+ FileFormatsChan = bst_get('FileFilters', 'channel');
+ FileFormat = FileFormatsChan{sProcess.options.format.Value{1}, 3};
+ else
+ FileFormat = [];
+ end
if NotFound
- bst_report('Info', sProcess, sInput, ...
- sprintf('Could not find original channel file: %s.', ChannelFile));
- % import_channel will prompt the user, but they will not
- % know which file to pick! And prompt is modal for Matlab,
- % so likely can't look at command window (e.g. if
- % Brainstorm is in front).
+ if isReport
+ bst_report('Info', sProcess, sInput, ...
+ sprintf('Could not find original channel file: %s.', ChannelFile));
+ end
+ % import_channel will prompt the user, but they would not know which file to pick! And
+ % prompt is modal for Matlab, so likely can't look at command window (e.g. if Brainstorm is
+ % in front). So add another pop-up with the needed info.
[ChanPath, ChanName, ChanExt] = fileparts(ChannelFile);
MsgFig = msgbox(sprintf('Select the new location of channel file %s %s to reset %s.', ...
ChanPath, [ChanName, ChanExt], sInput.ChannelFile), ...
@@ -318,50 +427,63 @@
movegui(MsgFig, 'north');
figure(MsgFig); % bring it to front.
% Adjust default format to the one selected.
- DefaultFormats = bst_get('DefaultFormats');
- DefaultFormats.ChannelIn = FileFormat;
- bst_set('DefaultFormats', DefaultFormats);
+ if ~isempty(FileFormat)
+ DefaultFormats = bst_get('DefaultFormats');
+ DefaultFormats.ChannelIn = FileFormat;
+ bst_set('DefaultFormats', DefaultFormats);
+ end
- [NewChannelMat, NewChannelFile] = import_channel(...
- sInput.iStudy, '', FileFormat, 0, 0, 0, [], []);
+ [NewChannelMat, NewChannelFile] = import_channel(sInput.iStudy, '', FileFormat, 0, 0, 0, [], []);
else
% Import from original file.
- [NewChannelMat, NewChannelFile] = import_channel(...
- sInput.iStudy, ChannelFile, FileFormat, 0, 0, 0, [], []);
+ [NewChannelMat, NewChannelFile] = import_channel(sInput.iStudy, ChannelFile, FileFormat, 0, 0, 0, [], 0);
% iStudies, ChannelFile, FileFormat, ChannelReplace, ChannelAlign, isSave, isFixUnits, isApplyVox2ras)
% iStudy index is needed to avoid error for noise recordings with missing SCS transform.
% ChannelReplace is for replacing the file, only if isSave.
% ChannelAlign is for headpoints, but also ONLY if isSave. We do it later if user selected.
+ % isApplyVox2ras should be false to avoid the pop-up asking if we want to use the MRI world
+ % transformation even before checking for digitized fids (it would then apply MRI world > MRI vox > MRI SCS).
end
% See if it worked.
if isempty(NewChannelFile)
- bst_report('Error', sProcess, sInput, ...
- 'No file channel file selected.');
- Failed = true;
+ if isReport
+ bst_report('Error', sProcess, sInput, 'No channel file selected.');
+ else
+ bst_error('No channel file selected.');
+ end
+ isError = true;
return;
elseif isempty(NewChannelMat)
- bst_report('Error', sProcess, sInput, ...
- sprintf('Unable to import channel file: %s', NewChannelFile));
- Failed = true;
+ if isReport
+ bst_report('Error', sProcess, sInput, sprintf('Unable to import channel file: %s', NewChannelFile));
+ else
+ bst_error(sprintf('Unable to import channel file: %s', NewChannelFile));
+ end
+ isError = true;
return;
elseif numel(NewChannelMat.Channel) ~= numel(ChannelMat.Channel)
- bst_report('Error', sProcess, sInput, ...
- 'Original channel file has different channels than current one, aborting.');
- Failed = true;
+ if isReport
+ bst_report('Error', sProcess, sInput, ...
+ 'Original channel file has different channels than current one, aborting.');
+ else
+ bst_error('Original channel file has different channels than current one, aborting.');
+ end
+ isError = true;
return;
elseif NotFound && ~isempty(ChannelFile)
% Save the selected new location.
NewChannelFiles(end+1, :) = {ChannelFile, NewChannelFile};
end
- % Copy the new old projectors and history to the new structure.
+ % Copy the old projectors and history to the new structure.
NewChannelMat.Projector = ChannelMat.Projector;
NewChannelMat.History = ChannelMat.History;
ChannelMat = NewChannelMat;
- % clear NewChannelMat
+
% Add number of channels to comment, like in db_set_channel.
ChannelMat.Comment = [ChannelMat.Comment, sprintf(' (%d)', length(ChannelMat.Channel))];
+ % Add history
ChannelMat = bst_history('add', ChannelMat, 'import', ...
['Reset from: ' NewChannelFile ' (Format: ' FileFormat ')']);
end % ResetChannelFile
@@ -380,9 +502,8 @@
% Need to check for empty, otherwise applies to all channels!
else
iChan = []; % All channels.
- % Note: NIRS doesn't have a separate set of
- % transformations, but "refine" and "SCS" are applied
- % to NIRS as well.
+ % Note: NIRS doesn't have a separate set of transformations, but "refine" and "SCS" are
+ % applied to NIRS as well.
end
while ~isempty(iUndoMeg)
if isMegOnly && isempty(iChan)
@@ -441,8 +562,8 @@
end % RemoveTransformation
-function [ChannelMat, Failed] = AdjustHeadPosition(ChannelMat, sInputs, sProcess)
- Failed = false;
+function [ChannelMat, isError] = AdjustHeadPosition(ChannelMat, sInputs, sProcess)
+ isError = false;
% Check the input is CTF.
isRaw = (length(sInputs(1).FileName) > 9) && ~isempty(strfind(sInputs(1).FileName, 'data_0raw'));
if isRaw
@@ -453,13 +574,12 @@
if ~strcmp(DataMat.Device, 'CTF')
bst_report('Error', sProcess, sInputs, ...
'Adjust head position is currently only available for CTF data.');
- Failed = true;
+ isError = true;
return;
end
- % The data could be changed such that the head position could be
- % readjusted (e.g. by deleting segments). This is allowed and the
- % previous adjustment will be replaced.
+ % The data could be changed such that the head position could be readjusted (e.g. by deleting
+ % segments). This is allowed and the previous adjustment will be replaced.
if isfield(ChannelMat, 'TransfMegLabels') && iscell(ChannelMat.TransfMegLabels) && ...
ismember('AdjustedNative', ChannelMat.TransfMegLabels)
bst_report('Info', sProcess, sInputs, ...
@@ -477,7 +597,7 @@
[Locs, HeadSamplePeriod] = process_evt_head_motion('LoadHLU', sInputs(iIn), [], false);
if isempty(Locs)
% No HLU channels. Error already reported. Skip this file.
- Failed = true;
+ isError = true;
return;
end
% Exclude all bad segments.
@@ -513,7 +633,7 @@
Locs(:, ismember(iHeadSamples, iBad)) = [];
end
end
- Locations = [Locations, Locs];
+ Locations = [Locations, Locs]; %#ok
end
% If a collection was aborted, the channels will be filled with zeros. Remove these.
@@ -521,57 +641,52 @@
bst_progress('text', 'Correcting head position...');
MedianLoc = MedianLocation(Locations);
+ if isempty(MedianLoc)
+ isError = true;
+ return;
+ end
% disp(MedianLoc);
- % Also get the initial reference position. We only use it to see
- % how much the adjustment moves.
+ % Also get the initial reference position. We only use it to estimate how much the adjustment moves.
InitRefLoc = ReferenceHeadLocation(ChannelMat, sInputs);
if isempty(InitRefLoc)
% There was an error, already reported. Skip this file.
- Failed = true;
+ isError = true;
return;
end
- % Extract transformations that are applied before and after the
- % head position adjustment. Any previous adjustment will be
- % ignored here and replaced later.
+ % Extract transformations that are applied before and after the head position adjustment. Any
+ % previous adjustment will be ignored here and replaced later.
[TransfBefore, TransfAdjust, TransfAfter, iAdjust, iDewToNat] = ...
GetTransforms(ChannelMat, sInputs);
if isempty(TransfBefore)
% There was an error, already reported. Skip this file.
- Failed = true;
+ isError = true;
return;
end
% Compute transformation corresponding to coil position.
[TransfMat, TransfAdjust] = LocationTransform(MedianLoc, ...
TransfBefore, TransfAdjust, TransfAfter);
- % This TransfMat would automatically give an identity
- % transformation if the process is run multiple times, and
- % TransfAdjust would not change.
-
- % Apply this transformation to the current head position.
- % This is a correction to the 'Dewar=>Native'
- % transformation so it applies to MEG channels only and not
- % to EEG or head points, which start in Native.
+ % This TransfMat would automatically give an identity transformation if the process is run
+ % multiple times, and TransfAdjust would not change.
+
+ % Apply this transformation to the current head position. This is a correction to the
+ % 'Dewar=>Native' transformation so it applies to MEG channels only and not to EEG or head
+ % points, which start in Native.
iMeg = sort([good_channel(ChannelMat.Channel, [], 'MEG'), ...
good_channel(ChannelMat.Channel, [], 'MEG REF')]);
ChannelMat = channel_apply_transf(ChannelMat, TransfMat, iMeg, false); % Don't apply to head points.
ChannelMat = ChannelMat{1};
- % After much thought, it was decided to save this
- % adjustment transformation separately and at its logical
- % place: between 'Dewar=>Native' and
- % 'Native=>Brainstorm/CTF'. In particular, this allows us
- % to use it directly when displaying head motion distance.
- % This however means we must correctly move the
- % transformation from the end where it was just applied to
- % its logical place. This "moved" transformation is also
- % computed in LocationTransform above.
+ % After much thought, it was decided to save this adjustment transformation separately and at
+ % its logical place: between 'Dewar=>Native' and 'Native=>Brainstorm/CTF'. In particular, this
+ % allows us to use it directly when displaying head motion distance. This however means we must
+ % correctly move the transformation from the end where it was just applied to its logical place.
+ % This "moved" transformation is also computed in LocationTransform above.
if isempty(iAdjust)
iAdjust = iDewToNat + 1;
- % Shift transformations to make room for the new
- % adjustment, and reject the last one, that we just
- % applied.
+ % Shift transformations to make room for the new adjustment, and reject the last one, that
+ % we just applied.
ChannelMat.TransfMegLabels(iDewToNat+2:end) = ...
ChannelMat.TransfMegLabels(iDewToNat+1:end-1); % reject last one
ChannelMat.TransfMeg(iDewToNat+2:end) = ChannelMat.TransfMeg(iDewToNat+1:end-1);
@@ -591,24 +706,22 @@
AfterRefLoc = ReferenceHeadLocation(ChannelMat, sInputs);
if isempty(AfterRefLoc)
% There was an error, already reported. Skip this file.
- Failed = true;
+ isError = true;
return;
end
DistanceAdjusted = process_evt_head_motion('RigidDistances', AfterRefLoc, InitRefLoc);
fprintf('Head position adjusted by %1.1f mm.\n', DistanceAdjusted * 1e3);
bst_report('Info', sProcess, sInputs, ...
- sprintf('Head position adjusted by %1.1f mm.\n', DistanceAdjusted * 1e3));
+ sprintf('Head position adjusted by %1.1f mm.', DistanceAdjusted * 1e3));
end % AdjustHeadPosition
-
function [InitLoc, Message] = ReferenceHeadLocation(ChannelMat, sInput)
% Compute initial head location in Dewar coordinates.
- % Here we want to recreate the correct triangle shape from the relative
- % head coil locations and in the position saved as the reference
- % (initial) head position according to Brainstorm coordinate
- % transformation matrices.
+ % Here we want to recreate the correct triangle shape from the relative head coil locations and
+ % in the position saved as the reference (initial) head position according to Brainstorm
+ % coordinate transformation matrices.
if nargin < 2
sInput = [];
@@ -616,23 +729,24 @@
sInput = sInput(1);
end
Message = '';
-
- % These aren't exactly the coil positions in the .hc file, which are not saved
- % anywhere in Brainstorm, but was verified to give the same transformation.
- % The SCS coil coordinates are from the digitized coil positions.
- if isfield(ChannelMat, 'SCS') && all(isfield(ChannelMat.SCS, {'NAS','LPA','RPA'})) && ...
- (length(ChannelMat.SCS.NAS) == 3) && (length(ChannelMat.SCS.LPA) == 3) && (length(ChannelMat.SCS.RPA) == 3)
- % % Use the SCS distances from origin, with left and right PA points symmetrical.
- % LeftRightDist = sqrt(sum((ChannelMat.SCS.LPA - ChannelMat.SCS.RPA).^2));
- % NasDist = ChannelMat.SCS.NAS(1);
- InitLoc = [ChannelMat.SCS.NAS(:), ChannelMat.SCS.LPA(:), ChannelMat.SCS.RPA(:); ones(1, 3)];
- elseif ~isempty(sInput) && isfield(sInput, 'header') && isfield(sInput.header, 'hc') && isfield(sInput.header.hc, 'SCS') && ...
+
+ % From recent investigations, digitized locations are probably not as robust/accurate as those
+ % measured by the MEG. So use the .hc positions if available.
+ if ~isempty(sInput) && isfield(sInput, 'header') && isfield(sInput.header, 'hc') && isfield(sInput.header.hc, 'SCS') && ...
all(isfield(sInput.header.hc.SCS, {'NAS','LPA','RPA'})) && length(sInput.header.hc.SCS.NAS) == 3
% Initial head coil locations from the CTF .hc file, but in dewar coordinates, NOT in SCS coordinates!
InitLoc = [sInput.header.hc.SCS.NAS(:), sInput.header.hc.SCS.LPA(:), sInput.header.hc.SCS.RPA(:)]; % 3x3 by columns
InitLoc = InitLoc(:);
return;
- %InitLoc = TransfAdjust * TransfBefore * [InitLoc; ones(1, 3)];
+ % ChannelMat.SCS are not the coil positions in the .hc file, which are not saved in Brainstorm,
+ % but the digitized coil positions, if present. However, both are saved in "Native" coordinates
+ % and thus give the same transformation.
+ elseif isfield(ChannelMat, 'SCS') && all(isfield(ChannelMat.SCS, {'NAS','LPA','RPA'})) && ...
+ (length(ChannelMat.SCS.NAS) == 3) && (length(ChannelMat.SCS.LPA) == 3) && (length(ChannelMat.SCS.RPA) == 3)
+ % % Use the SCS distances from origin, with left and right PA points symmetrical.
+ % LeftRightDist = sqrt(sum((ChannelMat.SCS.LPA - ChannelMat.SCS.RPA).^2));
+ % NasDist = ChannelMat.SCS.NAS(1);
+ InitLoc = [ChannelMat.SCS.NAS(:), ChannelMat.SCS.LPA(:), ChannelMat.SCS.RPA(:); ones(1, 3)];
else
% Just use some reasonable distances, with a warning.
Message = 'Exact reference head coil locations not available. Using reasonable (adult) locations according to head position.';
@@ -643,11 +757,9 @@
% InitLoc above is in Native coordiates (if pre head loc didn't fail).
% Bring it back to Dewar coordinates to compare with HLU channels.
%
- % Take into account if the initial/reference head position was
- % "adjusted", i.e. replaced by the median position throughout the
- % recording. If so, use all transformations from 'Dewar=>Native' to
- % this adjustment transformation. (In practice there shouldn't be any
- % between them.)
+ % Take into account if the initial/reference head position was "adjusted", i.e. replaced by the
+ % median position throughout the recording. If so, use all transformations from 'Dewar=>Native'
+ % to this adjustment transformation. (In practice there shouldn't be any between them.)
[TransfBefore, TransfAdjust] = GetTransforms(ChannelMat, sInput);
InitLoc = TransfBefore \ (TransfAdjust \ InitLoc);
InitLoc(4, :) = [];
@@ -655,25 +767,25 @@
end % ReferenceHeadLocation
-function [TransfBefore, TransfAdjust, TransfAfter, iAdjust, iDewToNat] = ...
- GetTransforms(ChannelMat, sInputs)
- % Extract transformations that are applied before and after the head
- % position adjustment we are creating now. We keep the 'Dewar=>Native'
- % transformation intact and separate from the adjustment for no deep
- % reason, but it is the only remaining trace of the initial head coil
+function [TransfBefore, TransfAdjust, TransfAfter, iAdjust, iDewToNat] = GetTransforms(ChannelMat, sInputs)
+ % Extract transformations that are applied before and after the head position adjustment we are
+ % creating now. We keep the 'Dewar=>Native' transformation intact and separate from the
+ % adjustment for no deep reason, but it is the only remaining trace of the initial head coil
% positions in Brainstorm.
- % The reason this function was split from LocationTransform is that it
- % can be called only once outside the head sample loop in process_sss,
- % whereas LocationTransform is called many times within the loop.
+ % The reason this function was split from LocationTransform is that it can be called only once
+ % outside the head sample loop in process_sss, whereas LocationTransform is called many times
+ % within the loop.
- % When this is called from process_sss, we are possibly working on a
- % second head adjustment, this time based on the instantaneous head
- % position, so we need to keep the global adjustment based on the entire
- % recording if it is there.
+ % When this is called from process_sss, we are possibly working on a second head adjustment,
+ % this time based on the instantaneous head position, so we need to keep the global adjustment
+ % based on the entire recording if it is there.
- % Check order of transformations. These situations should not happen
- % unless there was some manual editing.
+ if nargin < 2
+ sInputs = [];
+ end
+ % Check order of transformations. These situations should not happen unless there was some
+ % manual editing.
iDewToNat = find(strcmpi(ChannelMat.TransfMegLabels, 'Dewar=>Native'));
iAdjust = find(strcmpi(ChannelMat.TransfMegLabels, 'AdjustedNative'));
TransfBefore = [];
@@ -724,11 +836,9 @@
end % GetTransforms
-function [TransfMat, TransfAdjust] = LocationTransform(Loc, ...
- TransfBefore, TransfAdjust, TransfAfter)
- % Compute transformation corresponding to head coil positions.
- % We want this to be as efficient as possible, since used many times by
- % process_sss.
+function [TransfMat, TransfAdjust] = LocationTransform(Loc, TransfBefore, TransfAdjust, TransfAfter)
+ % Compute transformation corresponding to head coil positions. We want this to be as efficient
+ % as possible, since used many times by process_sss.
% Check for previous version.
if nargin < 4
@@ -736,10 +846,10 @@
end
% Transformation matrices are in m, as are HLU channels.
- % The HLU channels (here Loc) are in dewar coordinates. Bring them to
- % the current system by applying all saved transformations, starting with
- % 'Dewar=>Native'. This will save us from having to use inverse
- % transformations later.
+ %
+ % The HLU channels (here Loc) are in dewar coordinates. Bring them to the current system by
+ % applying all saved transformations, starting with 'Dewar=>Native'. This will save us from
+ % having to use inverse transformations later.
Loc = TransfAfter(1:3, :) * TransfAdjust * TransfBefore * [reshape(Loc, 3, 3); 1, 1, 1];
% [[Loc(1:3), Loc(4:6), Loc(5:9)]; 1, 1, 1]; % test if efficiency difference.
@@ -759,14 +869,12 @@
TransfMat(1:3,1:3) = [X, Y, Z]';
TransfMat(1:3,4) = - [X, Y, Z]' * Origin;
- % TransfMat at this stage is a transformation from the current system
- % back to the now adjusted Native system. We thus need to reapply the
- % following tranformations.
+ % TransfMat at this stage is a transformation from the current system back to the now adjusted
+ % Native system. We thus need to reapply the following tranformations.
if nargout > 1
- % Transform from non-adjusted native coordinates to newly adjusted native
- % coordinates. To be saved in channel file between "Dewar=>Native" and
- % "Native=>Brainstorm/CTF".
+ % Transform from non-adjusted native coordinates to newly adjusted native coordinates. To
+ % be saved in channel file between "Dewar=>Native" and "Native=>Brainstorm/CTF".
TransfAdjust = TransfMat * TransfAfter * TransfAdjust;
end
@@ -781,6 +889,8 @@
if size(Locations, 1) ~= 9
bst_error('Expecting 9 HLU channels in first dimension.');
+ MedianLoc = [];
+ return;
end
nSxnT = size(Locations, 2) * size(Locations, 3);
@@ -796,34 +906,19 @@
%
% M = GeoMedian(X, Precision)
%
- % Calculate the geometric median: the point that minimizes sum of
- % Euclidean distances to all points. size(X) = [n, d, ...], where n is
- % the number of data points, d is the number of components for each point
- % and any additional array dimension is treated as independent sets of
- % data and a median is calculated for each element along those dimensions
- % sequentially; size(M) = [1, d, ...]. This is an approximate iterative
- % procedure that stops once the desired precision is achieved. If
- % Precision is not provided, 1e-4 of the max distance from the centroid
- % is used.
- %
- % Weiszfeld's algorithm is used, which is a subgradient algorithm; with
- % (Verdi & Zhang 2001)'s modification to avoid non-optimal fixed points
- % (if at any iteration the approximation of M equals a data point).
- %
- %
- % (c) Copyright 2018 Marc Lalancette
- % The Hospital for Sick Children, Toronto, Canada
+ % Calculate the geometric median: the point that minimizes sum of Euclidean distances to all
+ % points. size(X) = [n, d, ...], where n is the number of data points, d is the number of
+ % components for each point and any additional array dimension is treated as independent sets of
+ % data and a median is calculated for each element along those dimensions sequentially; size(M)
+ % = [1, d, ...]. This is an approximate iterative procedure that stops once the desired
+ % precision is achieved. If Precision is not provided, 1e-4 of the max distance from the
+ % centroid is used.
%
- % This file is part of a free repository of Matlab tools for MEG
- % data processing and analysis .
- % 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 3 of the License, or (at your option) a later version.
+ % Weiszfeld's algorithm is used, which is a subgradient algorithm; with (Verdi & Zhang 2001)'s
+ % modification to avoid non-optimal fixed points (if at any iteration the approximation of M
+ % equals a data point).
%
- % This program is distributed WITHOUT ANY WARRANTY.
- % See the LICENSE file, or for details.
- %
- % 2012-05
+ % Marc Lalancette 2012-05
nDims = ndims(X);
XSize = size(X);
@@ -850,17 +945,15 @@
Precision = bsxfun(@rdivide, Precision, Scale); % Precision ./ Scale; % [1, 1, nSets]
end
- % Initial estimate: median in each dimension separately. Though this
- % gives a chance of picking one of the data points, which requires
- % special treatment.
+ % Initial estimate: median in each dimension separately. Though this gives a chance of picking
+ % one of the data points, which requires special treatment.
M2 = median(X, 1);
- % It might be better to calculate separately each independent set,
- % otherwise, they are all iterated until the worst case converges.
+ % It might be better to calculate separately each independent set, otherwise, they are all
+ % iterated until the worst case converges.
for s = 1:nSets
- % For convenience, pick another point far enough so the loop will always
- % start.
+ % For convenience, pick another point far enough so the loop will always start.
M = bsxfun(@plus, M2(:, :, s), Precision(:, :, s));
% Iterate.
while sum((M - M2(:, :, s)).^2 , 2) > Precision(s)^2 % any()scalar
@@ -868,8 +961,7 @@
% Distances from M.
% R = sqrt(sum( (M(ones(n, 1), :) - X(:, :, s)).^2 , 2 )); % [n, 1]
R = sqrt(sum( bsxfun(@minus, M, X(:, :, s)).^2 , 2 )); % [n, 1]
- % Find data points not equal to M, that we use in the computation
- % below.
+ % Find data points not equal to M, that we use in the computation below.
Good = logical(R);
nG = sum(Good);
if nG % > 0
@@ -883,10 +975,10 @@
end
% New estimate.
- % Note the possibility of D = 0 and (n - nG) = 0, in which case 0/0
- % should be 0, but here gives NaN, which the max function ignores,
- % returning 0 instead of 1. This is fine however since this
- % multiplies D (=0 in that case).
+ %
+ % Note the possibility of D = 0 and (n - nG) = 0, in which case 0/0 should be 0, but
+ % here gives NaN, which the max function ignores, returning 0 instead of 1. This is fine
+ % however since this multiplies D (=0 in that case).
M2(:, :, s) = M - max(0, 1 - (n - nG)/sqrt(sum( D.^2 , 2 ))) * ...
D / sum(1 ./ R, 1);
end
@@ -903,3 +995,560 @@
end % GeoMedian
+function [AlignType, isMriUpdated, isMriMatch, isSessionMatch, ChannelMat] = CheckPrevAdjustments(ChannelMat, sMri)
+ % Flag if auto or manual registration performed, and if MRI fids updated. Print to command
+ % window for now, if no output arguments. Also make sure to update ChannelMat.SCS if outputting
+ % ChannelMat.
+ AlignType = [];
+ isMriUpdated = [];
+ isMriMatch = [];
+ isSessionMatch = [];
+ isPrint = nargout == 0;
+ if any(~isfield(ChannelMat, {'History', 'HeadPoints'}))
+ % Nothing to check.
+ return;
+ end
+ if nargin < 2 || isempty(sMri) || ~isfield(sMri, 'History') || isempty(sMri.History)
+ iMriHist = [];
+ else
+ % History string is set in figure_mri SaveMri.
+ iMriHist = find(strcmpi(sMri.History(:,3), 'Applied digitized anatomical fiducials'), 1, 'last');
+ end
+ % Can also be reset, so check for 'import' action and ignore previous alignments.
+ iImport = find(strcmpi(ChannelMat.History(:,2), 'import'));
+ iAlign = find(strcmpi(ChannelMat.History(:,2), 'align'));
+ iAlign(iAlign < iImport(end)) = [];
+ if numel(iImport) > 1
+ AlignType = 'none/reset';
+ else
+ AlignType = 'none';
+ end
+ while ~isempty(iAlign)
+ % Check which adjustment was done last.
+ switch lower(ChannelMat.History{iAlign(end),3}(1:5))
+ case 'remov' % ['Removed transform: ' TransfLabel]
+ % Removed a previous step. Ignore corresponding adjustment and look again.
+ if strncmpi(ChannelMat.History{iAlign(end),3}(20:24), 'AdjustedNative', 5)
+ iAlignRemoved = find(cellfun(@(c)strcmpi(c(1:5), 'added'), ChannelMat.History(iAlign,3)), 1, 'last');
+ elseif strncmpi(ChannelMat.History{iAlign(end),3}(20:24), 'refine registration: head points', 5)
+ iAlignRemoved = find(cellfun(@(c)strcmpi(c(1:5), 'refin'), ChannelMat.History(iAlign,3)), 1, 'last');
+ elseif strncmpi(ChannelMat.History{iAlign(end),3}(20:24), 'manual correction', 5)
+ iAlignRemoved = find(cellfun(@(c)strcmpi(c(1:5), 'align'), ChannelMat.History(iAlign,3)), 1, 'last');
+ else
+ bst_error('Unrecognized removed transformation in history.');
+ return;
+ end
+ iAlign(end) = [];
+ if isempty(iAlignRemoved)
+ bst_error('Missing removed transformation in history.');
+ return;
+ else
+ iAlign(iAlignRemoved) = [];
+ end
+ case 'added' % 'Added adjustment to Native coordinates based on median head position'
+ % This alignment is between points and functional dataset, ignore here.
+ iAlign(end) = [];
+ case 'refin' % 'Refining the registration using the head points:'
+ % Automatic MRI-points alignment
+ AlignType = 'auto';
+ break;
+ case 'align' % 'Align channels manually:'
+ % Manual MRI-points alignment
+ AlignType = 'manual';
+ break;
+ case 'non-l' % 'Non-linear transformation'
+ AlignType = 'non-linear';
+ break;
+ otherwise
+ AlignType = 'unrecognized';
+ break;
+ end
+ end
+ if isPrint
+ disp(['BST> Previous registration adjustment: ' AlignType]);
+ end
+ if ~isempty(iMriHist)
+ isMriUpdated = true;
+ % Compare digitized fids to MRI fids (in MRI coordinates, mm). ChannelMat.SCS fids are NOT
+ % kept up to date when adjusting registration (manual or auto), so get them from head points
+ % again.
+ % Get the three fiducials in the head points
+ ChannelMat = UpdateChannelMatScs(ChannelMat);
+ % Check if coordinates differ by more than 1 um.
+ if isempty(ChannelMat.SCS.NAS) || isempty(ChannelMat.SCS.LPA) || isempty(ChannelMat.SCS.RPA)
+ isMriMatch = false;
+ isSessionMatch = false;
+ if isPrint
+ disp('BST> MRI fiducials previously updated, but different session than current (missing) digitized fiducials.');
+ end
+ elseif any(abs(sMri.SCS.NAS - cs_convert(sMri, 'scs', 'mri', ChannelMat.SCS.NAS) .* 1000) > 1e-3) || ...
+ any(abs(sMri.SCS.LPA - cs_convert(sMri, 'scs', 'mri', ChannelMat.SCS.LPA) .* 1000) > 1e-3) || ...
+ any(abs(sMri.SCS.RPA - cs_convert(sMri, 'scs', 'mri', ChannelMat.SCS.RPA) .* 1000) > 1e-3)
+ isMriMatch = false;
+ % Check if just different alignment, or if different set of fiducials (different
+ % session), using inter-fid distances.
+ DiffMri = [sMri.SCS.NAS - sMri.SCS.LPA, sMri.SCS.LPA - sMri.SCS.RPA, sMri.SCS.RPA - sMri.SCS.NAS];
+ DiffChannel = [ChannelMat.SCS.NAS - ChannelMat.SCS.LPA, ChannelMat.SCS.LPA - ChannelMat.SCS.RPA, ChannelMat.SCS.RPA - ChannelMat.SCS.NAS];
+ if any(abs(DiffMri - DiffChannel) > 1e-3)
+ isSessionMatch = false;
+ if isPrint
+ disp('BST> MRI fiducials previously updated, but different session than current digitized fiducials.');
+ end
+ else
+ isSessionMatch = true;
+ if isPrint
+ disp('BST> MRI fiducials previously updated, same session but not aligned with current digitized fiducials.');
+ end
+ end
+ else
+ isMriMatch = true;
+ isSessionMatch = true;
+ if isPrint
+ disp('BST> MRI fiducials previously updated, and match current digitized fiducials.');
+ end
+ end
+ else
+ if nargout > 4
+ % Update SCS for consistency.
+ % Get the three fiducials in the head points
+ ChannelMat = UpdateChannelMatScs(ChannelMat);
+ end
+ isMriUpdated = false;
+ isMriMatch = false;
+ isSessionMatch = false;
+ end
+end
+
+
+function [DistHead, DistSens, Message] = CheckCurrentAdjustments(ChannelMat, ChannelMatRef)
+ % Display max displacement from registration adjustments, in command window.
+ % If second ChannelMat is provided as reference, get displacement between the two.
+ isPrint = nargout == 0;
+ if nargin < 2 || isempty(ChannelMatRef)
+ ChannelMatRef = [];
+ end
+
+ % Update SCS from head points if present.
+ ChannelMat = UpdateChannelMatScs(ChannelMat);
+
+ if ~isempty(ChannelMatRef)
+ ChannelMatRef = UpdateChannelMatScs(ChannelMatRef);
+ % For head displacement, we use the "rigid distance" from the head motion code, basically
+ % the max distance of any point on a simplified spherical head.
+ DistHead = process_evt_head_motion('RigidDistances', ...
+ [ChannelMat.SCS.NAS(:); ChannelMat.SCS.LPA(:); ChannelMat.SCS.RPA(:)], ...
+ [ChannelMatRef.SCS.NAS(:); ChannelMatRef.SCS.LPA(:); ChannelMatRef.SCS.RPA(:)]);
+ DistSens = max(sqrt(sum(([ChannelMat.Channel.Loc] - [ChannelMatRef.Channel.Loc]).^2)));
+ else
+ % Implicitly using actual (MRI) SCS as reference, this includes all adjustments.
+ DistHead = process_evt_head_motion('RigidDistances', ...
+ [ChannelMat.SCS.NAS(:); ChannelMat.SCS.LPA(:); ChannelMat.SCS.RPA(:)]);
+ % Get equivalent transform for all adjustments to "undo" on sensors for comparison. The
+ % adjustments we want come after 'Native=>Brainstorm/CTF'
+ iNatToScs = find(strcmpi(ChannelMat.TransfMegLabels, 'Native=>Brainstorm/CTF'));
+ if iNatToScs < numel(ChannelMat.TransfMeg)
+ Transf = eye(4);
+ for t = iNatToScs+1:numel(ChannelMat.TransfMeg)
+ Transf = ChannelMat.TransfMeg{t} * Transf;
+ end
+ Loc = [ChannelMat.Channel.Loc];
+ % Inverse transf: subtract translation first, then rotate the "other way" (transpose).
+ LocRef = Transf(1:3,1:3)' * bsxfun(@minus, Loc, Transf(1:3,4));
+ DistSens = max(sqrt(sum((Loc - LocRef).^2)));
+ else
+ DistSens = 0;
+ end
+ end
+
+ Message = sprintf('BST> Max displacement for registration adjustment:\n head: %1.1f mm\n sensors: %1.1f cm\n', ...
+ DistHead*1000, DistSens*100);
+ if isPrint
+ fprintf(Message);
+ end
+
+end
+
+
+function ChannelMat = UpdateChannelMatScs(ChannelMat)
+ % Update the coordinates of the digitized anatomical fiducials in the channel SCS field, after
+ % potentially having edited the coregistration such that these points no longer define the SCS
+ % for the functional data - it's still defined by the MRI anatomical fiducials.
+ if ~isfield(ChannelMat, 'HeadPoints')
+ return;
+ end
+ % Get the three anatomical fiducials in the head points
+ iNas = find(strcmpi(ChannelMat.HeadPoints.Label, 'Nasion') | strcmpi(ChannelMat.HeadPoints.Label, 'NAS'));
+ iLpa = find(strcmpi(ChannelMat.HeadPoints.Label, 'Left') | strcmpi(ChannelMat.HeadPoints.Label, 'LPA'));
+ iRpa = find(strcmpi(ChannelMat.HeadPoints.Label, 'Right') | strcmpi(ChannelMat.HeadPoints.Label, 'RPA'));
+ if ~isempty(iNas) && ~isempty(iLpa) && ~isempty(iRpa)
+ ChannelMat.SCS.NAS = mean(ChannelMat.HeadPoints.Loc(:,iNas)', 1); %#ok<*UDIM>
+ ChannelMat.SCS.LPA = mean(ChannelMat.HeadPoints.Loc(:,iLpa)', 1);
+ ChannelMat.SCS.RPA = mean(ChannelMat.HeadPoints.Loc(:,iRpa)', 1);
+ % The SCS.R, T and Origin fields no longer have any use, except for missing digitized fids
+ % (see below), but keep them updated always for consistency.
+ [~, ChannelMat] = cs_compute(ChannelMat, 'scs');
+ end
+ % Do the same with head coils, used when exporting coregistration to BIDS
+ % Also, if only the head coils were digitized and no anatomical points, use these as SCS, which
+ % is what Brainstorm will implicitly use as SCS coordinates to align with the MRI anyway.
+ iHpiN = find(strcmpi(ChannelMat.HeadPoints.Label, 'HPI-N'));
+ iHpiL = find(strcmpi(ChannelMat.HeadPoints.Label, 'HPI-L'));
+ iHpiR = find(strcmpi(ChannelMat.HeadPoints.Label, 'HPI-R'));
+ if ~isempty(iHpiN) && ~isempty(iHpiL) && ~isempty(iHpiR)
+ % Temporarily put the head coils there to calculate transform.
+ ChannelMat.Native.NAS = mean(ChannelMat.HeadPoints.Loc(:,iHpiN)', 1);
+ ChannelMat.Native.LPA = mean(ChannelMat.HeadPoints.Loc(:,iHpiL)', 1);
+ ChannelMat.Native.RPA = mean(ChannelMat.HeadPoints.Loc(:,iHpiR)', 1);
+ % Get "current" SCS to Native transformation.
+ TmpChanMat = ChannelMat;
+ TmpChanMat.SCS = ChannelMat.Native;
+ % cs_compute doesn't change coordinates, only adds the R,T,Origin fields
+ % Digitized points are normally saved in Native coordinates, and converted to SCS if anat
+ % fids present. This transformation would then go back from SCS to Native.
+ [~, TmpChanMat] = cs_compute(TmpChanMat, 'scs');
+ ChannelMat.Native = TmpChanMat.SCS;
+ % If SCS missing (no anat points), Native matches SCS. Explicitly save SCS, which is missing
+ % from initial import.
+ if ~isfield(ChannelMat, 'SCS') || ~isfield(ChannelMat.SCS, 'NAS') || isempty(ChannelMat.SCS.NAS)
+ ChannelMat.SCS = ChannelMat.Native;
+ else
+ % Now apply the transform to the digitized anat fiducials. These are not used anywhere yet,
+ % only the transform, but might as well be consistent and save the same points as in .SCS
+ % (digitized anat fids). Still in meters, not cm, despite actual native CTF being in cm.
+ ChannelMat.Native.NAS(:) = [ChannelMat.Native.R, ChannelMat.Native.T] * [ChannelMat.SCS.NAS'; 1];
+ ChannelMat.Native.LPA(:) = [ChannelMat.Native.R, ChannelMat.Native.T] * [ChannelMat.SCS.LPA'; 1];
+ ChannelMat.Native.RPA(:) = [ChannelMat.Native.R, ChannelMat.Native.T] * [ChannelMat.SCS.RPA'; 1];
+ end
+ else
+ if ~isempty(iNas) && ~isempty(iLpa) && ~isempty(iRpa)
+ % Missing digitized MEG head coils, probably the anatomical points are actually coils.
+ % No study, subject or file name available to print here.
+ disp('BST> Missing digitized MEG head coils (in channel file), assuming NAS/LPA/RPA are actually head coils, but they should be renamed.');
+ ChannelMat.Native = ChannelMat.SCS;
+ else
+ ChannelMat.Native.R = [];
+ ChannelMat.Native.T = [];
+ disp('BST> No digitized fiducials, neither anatomical nor MEG head coils.');
+ end
+ end
+end
+
+
+% Decided to bring this back as subfunction of this process, as it is the only place to run it from for now.
+function [Transform, isCancel, isError] = channel_align_scs(ChannelFile, Transform, isInteractive, isConfirm, sInput, sProcess)
+% CHANNEL_ALIGN_SCS: Saves new MRI anatomical points after manual or auto registration adjustment.
+%
+% USAGE: Transform = channel_align_scs(ChannelFile, Transform=eye(4), isInteractive=1, isConfirm=1)
+%
+% DESCRIPTION:
+% After modifying registration between digitized head points and MRI (with "refine with head
+% points" or manually), this function allows saving the change in the MRI fiducials so that
+% they exactly match the digitized anatomical points (nasion and ears). This would replace
+% having to save a registration adjustment transformation for each functional dataset sharing
+% this set of digitized points. This affects all files registered to the MRI and should
+% therefore be done as one of the first steps after importing, and with only one set of
+% digitized points (one session). Surfaces are adjusted to maintain alignment with the MRI.
+% Additional sessions for the same Brainstorm subject, with separate digitized points, will
+% still need the usual "per dataset" registration adjustment to align with the same MRI.
+%
+% This function will not modify an MRI that it changed previously without user confirmation
+% (if both isInteractive and isConfirm are false). In that case, Transform is returned unaltered.
+%
+% INPUTS:
+% - ChannelFile : Channel file to align with its anatomy
+% - Transform : Transformation matrix from digitized SCS coordinates to MRI SCS coordinates,
+% after some alignment is made (auto or manual) and the two no longer match.
+% This transform should not already be saved in the ChannelFile, though the
+% file may already contain similar adjustments, in which case Transform would be
+% an additional adjustment to add. (This will typically be empty or identity, it
+% was intended for calling from manual alignment panel, but now done after.)
+% - isInteractive : If true, display dialog in case of errors, or if this was already done
+% previously for this MRI.
+% - isConfirm : If true, ask the user for confirmation before proceeding.
+%
+% OUTPUTS:
+% - Transform : If the MRI fiducial points and coordinate system are updated, the transform
+% becomes the identity. If the MRI was not updated, the input Transform is
+% returned. The idea is that the returned Transform applied to the *reset*
+% channels would maintain the registration. If channel files were not reset
+% (error or cancellation in this function), this will no longer be true and the
+% user should verify the registration of all affected studies.
+% - isCancel : If true, nothing was changed nor saved.
+% - isError : An error occurred that can affect registration of MRI and functional studies,
+% e.g. the MRI was updated, but some channel files of that subject were not
+% reset.
+
+% The Transform output is currently unused. It was changed (below is the previous behavior) since it
+% depended on CTF MEG specific transform labels.
+% - Transform : If the MRI fiducial points and coordinate system are updated, and the channel
+% file is reset, the transform becomes the identity. If the channel file is not
+% reset, Transform will be the inverse of all previous manual or automatic
+% adjustments. If the MRI was not updated, the input Transform is returned. The
+% idea is that the returned Transform applied to the channels would maintain the
+% registration.
+
+% Authors: Marc Lalancette 2022-2025
+
+if nargin < 6 || isempty(sProcess)
+ isReport = false;
+else
+ isReport = true;
+end
+if nargin < 4 || isempty(isConfirm)
+ isConfirm = true;
+end
+if nargin < 3 || isempty(isInteractive)
+ isInteractive = true;
+end
+if nargin < 2 || isempty(Transform)
+ Transform = eye(4);
+end
+if nargin < 1 || isempty(ChannelFile)
+ bst_error('ChannelFile argument required.');
+end
+
+isError = false;
+isCancel = false;
+% Get study
+sStudy = bst_get('ChannelFile', ChannelFile);
+% Get subject
+sSubject = bst_get('Subject', sStudy.BrainStormSubject);
+% Check if default anatomy.
+if sSubject.UseDefaultAnat
+ Message = 'Digitized nasion and ear points cannot be applied to default anatomy.';
+ if isReport
+ bst_report('Error', sProcess, sInput, Message);
+ elseif isInteractive
+ bst_error(Message, 'Apply digitized anatomical fiducials to MRI', 0);
+ else
+ disp(['BST> ' Message]);
+ end
+ isCancel = true;
+ return;
+end
+% Get Channels
+ChannelMat = in_bst_channel(ChannelFile);
+
+% Check if digitized anat points present, saved in ChannelMat.SCS.
+% Note that these coordinates are NOT currently updated when doing refine with head points (below).
+% They are in "initial SCS" coordinates, updated in channel_detect_type.
+if ~all(isfield(ChannelMat.SCS, {'NAS','LPA','RPA'})) || ~(length(ChannelMat.SCS.NAS) == 3) || ~(length(ChannelMat.SCS.LPA) == 3) || ~(length(ChannelMat.SCS.RPA) == 3)
+ Message = 'Digitized nasion and ear points not found.';
+ if isReport
+ bst_report('Error', sProcess, sInput, Message);
+ elseif isInteractive
+ bst_error(Message, 'Apply digitized anatomical fiducials to MRI', 0);
+ else
+ disp(['BST> ' Message]);
+ end
+ isCancel = true;
+ return;
+end
+
+% Check if already adjusted
+sMriOld = in_mri_bst(sSubject.Anatomy(sSubject.iAnatomy).FileName);
+% This Check function also updates ChannelMat.SCS with the saved (possibly previously adjusted) head
+% points IF isMriUpdated. (We don't consider isMriMatch here because we still have to apply the
+% provided Transformation.)
+[~, isMriUpdated, ~, ~, ChannelMat] = CheckPrevAdjustments(ChannelMat, sMriOld);
+% Get user confirmation
+if isMriUpdated
+ % Already done previously.
+ if isInteractive || isConfirm
+ % Request confirmation.
+ [Proceed, isCancel] = java_dialog('confirm', ['The MRI fiducial points NAS/LPA/RPA were previously updated from a set of' 10 ...
+ 'aligned digitized points. Updating them again will break any previous alignment' 10 ...
+ 'with other sets of digitized points and associated functional datasets.' 10 10 ...
+ 'Proceed and overwrite previous alignment?' 10], 'Head points/anatomy registration');
+ if ~Proceed || isCancel
+ isCancel = true;
+ return;
+ end
+ else
+ % Do not proceed.
+ Message = 'Digitized nasion and ear points previously applied to this MRI. Not applying again.';
+ if isReport
+ bst_report('Warning', sProcess, sInput, Message);
+ else
+ disp(['BST> ' Message]);
+ end
+ isCancel = true;
+ return;
+ end
+elseif isConfirm
+ % Request confirmation.
+ [Proceed, isCancel] = java_dialog('confirm', ['Updating the MRI fiducial points NAS/LPA/RPA to match a set of' 10 ...
+ 'aligned digitized points is mainly used for exporting registration to a BIDS dataset.' 10 ...
+ 'It will break any previous alignment of this subject with all other functional datasets!' 10 10 ...
+ 'Proceed and update MRI now?' 10], 'Head points/anatomy registration');
+ if ~Proceed || isCancel
+ isCancel = true;
+ return;
+ end
+end
+%% TEMPORARY BYPASS OF EEG POP-UP
+if isConfirm
+% If EEG, warn that only linear transformation would be saved this way.
+if ~isempty([good_channel(ChannelMat.Channel, [], 'EEG'), good_channel(ChannelMat.Channel, [], 'SEEG'), good_channel(ChannelMat.Channel, [], 'ECOG')])
+ [Proceed, isCancel] = java_dialog('confirm', ['Updating the MRI fiducial points NAS/LPA/RPA will only save' 10 ...
+ 'global rotations and translations. Any other changes to EEG channels will be lost.' 10 10 ...
+ 'Proceed and update MRI now?' 10], 'Head points/anatomy registration');
+ if ~Proceed || isCancel
+ isCancel = true;
+ return;
+ end
+end
+end
+
+% Convert digitized fids to MRI SCS coordinates.
+% Here, ChannelMat.SCS already may contain some auto/manual adjustment, and we're adding a new one (possibly identity).
+% Apply the transformation provided.
+sMri = sMriOld;
+% Intermediate step, these are not valid coordinates for sMri.
+sMri.SCS.NAS = (Transform(1:3,:) * [ChannelMat.SCS.NAS'; 1])';
+sMri.SCS.LPA = (Transform(1:3,:) * [ChannelMat.SCS.LPA'; 1])';
+sMri.SCS.RPA = (Transform(1:3,:) * [ChannelMat.SCS.RPA'; 1])';
+% Then convert to MRI coordinates (mm), this is how sMri.SCS is saved.
+% cs_convert mri is in meters
+sMri.SCS.NAS = cs_convert(sMriOld, 'scs', 'mri', sMri.SCS.NAS) .* 1000;
+sMri.SCS.LPA = cs_convert(sMriOld, 'scs', 'mri', sMri.SCS.LPA) .* 1000;
+sMri.SCS.RPA = cs_convert(sMriOld, 'scs', 'mri', sMri.SCS.RPA) .* 1000;
+% Re-compute transformation in this struct, which goes from MRI to SCS (but the fids stay in MRI coords in this struct).
+[~, sMri] = cs_compute(sMri, 'scs');
+
+% Compare with existing MRI fids, replace if changed (> 1um), and update surfaces.
+sMri.FileName = sSubject.Anatomy(sSubject.iAnatomy).FileName;
+figure_mri('SaveMri', sMri);
+% At minimum, we must unload surfaces that have been modified, but we want to avoid closing figures
+% for when we show "before" and "after" figures.
+bst_memory('UnloadAll'); % not 'forced' so won't close figures, but won't unload what we want most likely.
+bst_memory('UnloadSurface'); % this unloads the head surface as we want.
+% Now that MRI is saved, update Transform to identity.
+Transform = eye(4);
+
+% MRI SCS now matches digitized-points-defined SCS (defined from same points), but registration is
+% now broken with all channel files that were adjusted! Reset channel file, and all others for this
+% anatomy.
+isError = ResetChannelFiles(ChannelMat, sSubject, isConfirm, sInput, sProcess);
+
+% Removed this output Transform for now as GetTransform only works on CTF MEG data and the rest of
+% this function can work more generally.
+% if isError
+% % Get the equivalent overall registration adjustment transformation previously saved.
+% [~, ~, TransfAfter] = GetTransforms(ChannelMat);
+% % Return its inverse as it's now part of the MRI and should be removed from the channel file.
+% Transform = inverse(TransfAfter);
+% else
+% Transform = eye(4);
+% end
+
+end % main function
+
+% (This function was based on channel_align_manual CopyToOtherFolders).
+function [isError, Message] = ResetChannelFiles(ChannelMatSrc, sSubject, isConfirm, sInput, sProcess)
+ if nargin < 5 || isempty(sProcess)
+ sProcess = [];
+ isReport = false;
+ else
+ isReport = true;
+ end
+ % Confirmation: ask the first time
+ if nargin < 3 || isempty(isConfirm)
+ isConfirm = true;
+ end
+
+ NewChannelFiles = cell(0,2);
+ % First, always reset the "source" channel file.
+ [ChannelMatSrc, NewChannelFiles, isError] = ResetChannelFile(ChannelMatSrc, NewChannelFiles, sInput, sProcess);
+ if isError
+ Message = sprintf(['Unable to reset channel file for subject: %s\n' ...
+ 'MRI registration for all their functional studies should be verified!'], sSubject.Name);
+ if isReport
+ bst_report('Error', sProcess, sInput, Message);
+ end
+ % This is very important so always show it interactively.
+ java_dialog('msgbox', Message);
+ return;
+ end
+ bst_save(file_fullpath(sInput.ChannelFile), ChannelMatSrc, 'v7');
+
+ % If the subject is configured to share its channel files, nothing to do
+ if (sSubject.UseDefaultChannel >= 1)
+ return;
+ end
+ % Get all the dependent studies
+ sStudies = bst_get('StudyWithSubject', sSubject.FileName);
+ % List of channel files to update
+ ChannelFiles = {};
+ % Loop on the other folders
+ for i = 1:length(sStudies)
+ % Skip studies without channel files
+ if isempty(sStudies(i).Channel) || isempty(sStudies(i).Channel(1).FileName)
+ continue;
+ end
+ % Add channel file to list of files to process
+ ChannelFiles{end+1} = sStudies(i).Channel(1).FileName; %#ok
+ end
+ % Unique files and skip "source".
+ ChannelFiles = setdiff(unique(ChannelFiles), sInput.ChannelFile);
+ if ~isempty(ChannelFiles)
+ % Ask confirmation to the user
+ if isConfirm
+ Proceed = java_dialog('confirm', ...
+ sprintf('Reset all %d other channel files for this subject (typically recommended)?', numel(ChannelFiles)), 'Reset channel files');
+ if ~Proceed
+ Message = sprintf(['User cancelled resetting %d other channel files for subject: %s\n' ...
+ 'MRI registration for all functional studies should be verified!'], numel(ChannelFiles), sSubject.Name);
+ if isReport
+ bst_report('Warning', sProcess, sInput, Message);
+ else
+ java_dialog('msgbox', Message);
+ end
+ return;
+ end
+ end
+ % Progress bar
+ bst_progress('start', 'Reset channel files', 'Updating other studies...');
+ strMsg = [sInput.ChannelFile 10];
+ strErr = '';
+ for iChan = 1:numel(ChannelFiles)
+ ChannelFile = ChannelFiles{iChan};
+ % Load channel file
+ ChannelMat = in_bst_channel(ChannelFile);
+ % Reset & save
+ % Need correct sInput for each study here, not the "source" study. Only 3 fields used.
+ [sStudyTmp, iStudy] = bst_get('ChannelFile', ChannelFile);
+ sInputTmp.FileName = sStudyTmp(1).Data.FileName; % for report only
+ sInputTmp.iStudy = iStudy; % this is the study getting reset
+ sInputTmp.ChannelFile = ChannelFile; % this is the file being read
+ [ChannelMat, NewChannelFiles, isError] = ResetChannelFile(ChannelMat, NewChannelFiles, sInputTmp, sProcess);
+ if isError
+ strErr = [strErr ChannelFile 10]; %#ok
+ else
+ strMsg = [strMsg ChannelFile 10]; %#ok
+ bst_save(file_fullpath(ChannelFile), ChannelMat, 'v7');
+ end
+ end
+ bst_progress('stop');
+ % Give report to the user
+ if ~isempty(strErr)
+ Message = sprintf(['Unable to reset channel file(s) for subject %s:\n%s\n' ...
+ 'MRI registration should be verified for these studies!'], sSubject.Name, strErr);
+ if isReport
+ bst_report('Error', sProcess, sInput, Message);
+ end
+ % This is very important so always show it interactively.
+ java_dialog('msgbox', Message);
+ return;
+ end
+ Message = sprintf('%d channel files reset for subject %s:\n%s', numel(ChannelFiles)+1, sSubject.Name, strMsg);
+ if isReport
+ bst_report('Info', sProcess, sInput, Message);
+ elseif isConfirm
+ java_dialog('msgbox', Message);
+ else
+ disp(Message);
+ end
+ end
+end
+
diff --git a/toolbox/process/functions/process_evt_head_motion.m b/toolbox/process/functions/process_evt_head_motion.m
index f319e7eb71..6fbd2f54be 100644
--- a/toolbox/process/functions/process_evt_head_motion.m
+++ b/toolbox/process/functions/process_evt_head_motion.m
@@ -291,6 +291,12 @@
% Load and downsample continuous head localization channels.
% HeadSamplePeriod is in (MEG) samples per (head) sample, not seconds.
% Locations are in meters, [nChannels, nSamples, nEpochs] possibly converted to continuous.
+ % Data is returned for each epoch (concatenated if converted to continuous), for the provided
+ % sample bounds.
+ % The coordinate system is "CTF dewar", with X & Y axes at 45 degrees, and the origin inside the
+ % dewar.
+ % However, with CTF software version 3847, the extra coil is mislocalized near the coordinate
+ % system origin, which doesn't make sense for any coordinate system.
% For now removing bad segments is done in process_adjust_coordinates only.
% , RemoveBadSegments
@@ -467,8 +473,7 @@
if nargin < 3 || isempty(StopThreshold)
StopThreshold = false;
end
-
- if size(Locations, 1) ~= 9 || size(Reference, 1) ~= 9
+ if size(Locations, 1) ~= 9 %|| size(Reference, 1) ~= 9
bst_error('Expecting 9 HLU channels in first dimension.');
end
nS = size(Locations, 2);
@@ -476,15 +481,24 @@
% Calculate distances.
- Reference = reshape(Reference, [3, 3]);
- % Reference "head origin" and inverse "orientation matrix".
- [YO, YR] = RigidCoordinates(Reference);
- % Sphere radius.
- r = max( sqrt(sum((Reference - YO(:, [1, 1, 1])).^2, 1)) );
- if any(YR(:)) % any ignores NaN and returns false for empty.
- YI = inv(YR); % Faster to calculate inverse once here than "/" in loop.
+ if nargin < 2 || isempty(Reference)
+ % Assume reference defines coordinate system.
+ YO = zeros(3, 1);
+ YI = eye(3);
+ % Use first location to estimate sphere radius.
+ r = max(sqrt(sum(bsxfun(@minus, reshape(Locations(1:9), [3,3]), ...
+ (Locations(4:6) + Locations(7:9))/2).^2, 1)));
else
- YI = YR;
+ Reference = reshape(Reference, [3, 3]);
+ % Reference "head origin" and inverse "orientation matrix".
+ [YO, YR] = RigidCoordinates(Reference);
+ % Sphere radius, estimate from reference.
+ r = max( sqrt(sum((Reference - YO(:, [1, 1, 1])).^2, 1)) );
+ if any(YR(:)) % any ignores NaN and returns false for empty.
+ YI = inv(YR); % Faster to calculate inverse once here than "/" in loop.
+ else
+ YI = YR;
+ end
end
% SinHalf = zeros([nS, 1, nT]);
@@ -499,21 +513,8 @@
% it is a rotation around an axis through the real origin).
R = XR * YI; % %#ok
- % Sine of half the rotation angle.
- % SinHalf = sqrt(3 - trace(R)) / 2;
- % For very small angles, this formula is not accurate compared to
- % w, since diagonal elements are around 1, and eps(1) = 2.2e-16.
- % This will be the order of magnitude of non-diag. elements due to
- % errors. So we should get SinHalf from w.
- % Rotation axis with amplitude = SinHalf (like in rotation quaternions).
- w = [R(3, 2) - R(2, 3); R(1, 3) - R(3, 1); R(2, 1) - R(1, 2)] / ...
- (2 * sqrt(1 + R(1, 1) + R(2, 2) + R(3, 3)));
- SinHalf = sqrt(sum(w.^2));
- TNormSq = sum(T.^2);
- % Maximum sphere distance for translation + rotation, as described
- % above.
- D(s, t) = sqrt( TNormSq + (2 * r * SinHalf)^2 + ...
- 4 * r * sqrt(TNormSq * SinHalf^2 - (T' * w)^2) );
+ % Maximum sphere distance for translation + rotation, as described above.
+ D(s, t) = RigidDistTransform(R, T, r);
% CHECK should be comparable AND >= to max coil movement.
% Option to interrupt when past a distance threshold.
@@ -527,6 +528,28 @@
+function D = RigidDistTransform(R, T, rad)
+ % Maximum sphere distance for translation + rotation, as described above.
+ if isempty(T) && size(R,1) == 4
+ T = R(1:3, 4);
+ R = R(1:3, 1:3);
+ end
+ % Sine of half the rotation angle.
+ % SinHalf = sqrt(3 - trace(R)) / 2;
+ % For very small angles, this formula is not accurate compared to w, since diagonal
+ % elements are around 1, and eps(1) = 2.2e-16. This will be the order of magnitude of
+ % non-diag. elements due to errors. So we should get SinHalf from w.
+ % Rotation axis with amplitude = SinHalf (like in rotation quaternions).
+ w = [R(3, 2) - R(2, 3); R(1, 3) - R(3, 1); R(2, 1) - R(1, 2)] / ...
+ (2 * sqrt(1 + R(1, 1) + R(2, 2) + R(3, 3)));
+ SinHalf = sqrt(sum(w.^2));
+ TNormSq = sum(T.^2);
+
+ D = sqrt( TNormSq + (2 * rad * SinHalf)^2 + 4 * rad * sqrt(TNormSq * SinHalf^2 - (T' * w)^2) );
+end % RigidDistTransform
+
+
+
function [O, R] = RigidCoordinates(FidsColumns)
% Convert head coil locations to origin position and rotation matrix.
% Works with 9x1 or 3x3 (columns) input.
diff --git a/toolbox/process/functions/process_headpoints_refine.m b/toolbox/process/functions/process_headpoints_refine.m
index 9cf995a6f7..f6e4f6db9c 100644
--- a/toolbox/process/functions/process_headpoints_refine.m
+++ b/toolbox/process/functions/process_headpoints_refine.m
@@ -41,7 +41,7 @@
sProcess.options.title.Comment = [...
'Refine the MEG/MRI registration using digitized head points.
' ...
'If (tolerance > 0): fit the head points, remove the digitized points the most
' ...
- 'distant to the scalp surface, and fit again the the head points on the scalp.
'];
+ 'distant to the scalp surface, and fit again the head points on the scalp.
'];
sProcess.options.title.Type = 'label';
% Tolerance
sProcess.options.tolerance.Comment = 'Tolerance (outlier points to ignore):';
@@ -61,13 +61,13 @@
% Get options
tolerance = sProcess.options.tolerance.Value{1} / 100;
% Get all the channel files
- uniqueChan = unique({sInputs.ChannelFile});
+ [uniqueChan, iUniqFiles] = unique({sInputs.ChannelFile});
% Loop on all the channel files
for i = 1:length(uniqueChan)
% Refine registration
[ChannelMat, R, T, isSkip, isUserCancel, strReport] = channel_align_auto(uniqueChan{i}, [], 0, 0, tolerance);
if ~isempty(strReport)
- bst_report('Info', sProcess, sInputs, strReport);
+ bst_report('Info', sProcess, sInputs(iUniqFiles(i)), strReport);
end
end
% Return all the files in input
diff --git a/toolbox/process/functions/process_import_bids.m b/toolbox/process/functions/process_import_bids.m
index c4ae818a1c..3edf4c00dc 100644
--- a/toolbox/process/functions/process_import_bids.m
+++ b/toolbox/process/functions/process_import_bids.m
@@ -533,7 +533,7 @@
errorMsg = [errorMsg, 10, errMsg];
end
% Generate head surface
- tess_isohead(iSubject, 10000, 0, 2);
+ tess_isohead(iSubject, 15000, 0, 0);
else
MrisToRegister{end+1} = BstMriFile;
end
diff --git a/toolbox/sensors/channel_align_auto.m b/toolbox/sensors/channel_align_auto.m
index 86671ea897..71d8fcd48e 100644
--- a/toolbox/sensors/channel_align_auto.m
+++ b/toolbox/sensors/channel_align_auto.m
@@ -14,7 +14,7 @@
%
% INPUTS:
% - ChannelFile : Channel file to align on its anatomy
-% - ChannelMat : If specified, do not read or write any information from/to ChannelFile
+% - ChannelMat : If specified, do not read or write any information from/to ChannelFile (except to get scalp surface).
% - isWarning : If 1, display warning in case of errors (default = 1)
% - isConfirm : If 1, ask the user for confirmation before proceeding
% - tolerance : Percentage of outliers head points, ignored in the final fit
@@ -101,7 +101,7 @@
sSubject = bst_get('Subject', sStudy.BrainStormSubject);
if isempty(sSubject) || isempty(sSubject.iScalp)
if isWarning
- bst_error('No scalp surface available for this subject', 'Align EEG sensors', 0);
+ bst_error('No scalp surface available for this subject', 'Automatic EEG-MEG/MRI registration', 0);
else
disp('BST> No scalp surface available for this subject.');
end
@@ -162,19 +162,19 @@
%% ===== FIND OPTIMAL FIT =====
% Find best possible rigid transformation (rotation+translation)
-[R,T,tmp,dist] = bst_meshfit(SurfaceMat.Vertices, SurfaceMat.Faces, HP);
+[R,T,tmp,dist] = bst_meshfit(SurfaceMat.Vertices, SurfaceMat.Faces, HP, tolerance);
% Remove outliers and fit again
-if ~isempty(dist) && ~isempty(tolerance) && (tolerance > 0)
- % Sort points by distance to scalp
- [tmp__, iSort] = sort(dist, 1, 'descend');
- iRemove = iSort(1:nRemove);
- % Remove from list of destination points
- HP(iRemove,:) = [];
- % Fit again
- [R,T,tmp,dist] = bst_meshfit(SurfaceMat.Vertices, SurfaceMat.Faces, HP);
-else
- nRemove = 0;
-end
+% if ~isempty(dist) && ~isempty(tolerance) && (tolerance > 0)
+% % Sort points by distance to scalp
+% [tmp__, iSort] = sort(dist, 1, 'descend');
+% iRemove = iSort(1:nRemove);
+% % Remove from list of destination points
+% HP(iRemove,:) = [];
+% % Fit again
+% [R,T,tmp,dist] = bst_meshfit(SurfaceMat.Vertices, SurfaceMat.Faces, HP);
+% else
+% nRemove = 0;
+% end
% Current position cannot be optimized
if isempty(R)
bst_progress('stop');
@@ -190,24 +190,32 @@
' | Number of outlier points removed: ' sprintf('%d (%d%%)', nRemove, round(tolerance*100)), 10 ...
' | Initial number of head points: ' num2str(size(HeadPoints.Loc,2))];
+% Create [4,4] transform matrix from digitized SCS to MRI SCS according to this fit.
+DigToMriTransf = eye(4);
+DigToMriTransf(1:3,1:3) = R;
+DigToMriTransf(1:3,4) = T;
+
%% ===== ROTATE SENSORS AND HEADPOINTS =====
-for i = 1:length(ChannelMat.Channel)
- % Rotate and translate location of channel
- if ~isempty(ChannelMat.Channel(i).Loc) && ~all(ChannelMat.Channel(i).Loc(:) == 0)
- ChannelMat.Channel(i).Loc = R * ChannelMat.Channel(i).Loc + T * ones(1,size(ChannelMat.Channel(i).Loc, 2));
+if ~isequal(DigToMriTransf, eye(4))
+ for i = 1:length(ChannelMat.Channel)
+ % Rotate and translate location of channel
+ if ~isempty(ChannelMat.Channel(i).Loc) && ~all(ChannelMat.Channel(i).Loc(:) == 0)
+ ChannelMat.Channel(i).Loc = R * ChannelMat.Channel(i).Loc + T * ones(1,size(ChannelMat.Channel(i).Loc, 2));
+ end
+ % Only rotate normal vector to channel
+ if ~isempty(ChannelMat.Channel(i).Orient) && ~all(ChannelMat.Channel(i).Orient(:) == 0)
+ ChannelMat.Channel(i).Orient = R * ChannelMat.Channel(i).Orient;
+ end
end
- % Only rotate normal vector to channel
- if ~isempty(ChannelMat.Channel(i).Orient) && ~all(ChannelMat.Channel(i).Orient(:) == 0)
- ChannelMat.Channel(i).Orient = R * ChannelMat.Channel(i).Orient;
+ % Rotate and translate head points
+ if isfield(ChannelMat, 'HeadPoints') && ~isempty(ChannelMat.HeadPoints) && ~isempty(ChannelMat.HeadPoints.Loc)
+ ChannelMat.HeadPoints.Loc = R * ChannelMat.HeadPoints.Loc + ...
+ T * ones(1, size(ChannelMat.HeadPoints.Loc, 2));
end
end
-% Rotate and translate head points
-if isfield(ChannelMat, 'HeadPoints') && ~isempty(ChannelMat.HeadPoints) && ~isempty(ChannelMat.HeadPoints.Loc)
- ChannelMat.HeadPoints.Loc = R * ChannelMat.HeadPoints.Loc + ...
- T * ones(1, size(ChannelMat.HeadPoints.Loc, 2));
-end
%% ===== SAVE TRANSFORMATION =====
+% We could decide to skip this if the transformation is identity.
% Initialize fields
if ~isfield(ChannelMat, 'TransfEeg') || ~iscell(ChannelMat.TransfEeg)
ChannelMat.TransfEeg = {};
@@ -221,13 +229,9 @@
if ~isfield(ChannelMat, 'TransfEegLabels') || ~iscell(ChannelMat.TransfEegLabels) || (length(ChannelMat.TransfEeg) ~= length(ChannelMat.TransfEegLabels))
ChannelMat.TransfEegLabels = cell(size(ChannelMat.TransfEeg));
end
-% Create [4,4] transform matrix
-newtransf = eye(4);
-newtransf(1:3,1:3) = R;
-newtransf(1:3,4) = T;
% Add a rotation/translation to the lists
-ChannelMat.TransfMeg{end+1} = newtransf;
-ChannelMat.TransfEeg{end+1} = newtransf;
+ChannelMat.TransfMeg{end+1} = DigToMriTransf;
+ChannelMat.TransfEeg{end+1} = DigToMriTransf;
% Add the comments
ChannelMat.TransfMegLabels{end+1} = 'refine registration: head points';
ChannelMat.TransfEegLabels{end+1} = 'refine registration: head points';
diff --git a/toolbox/sensors/channel_align_manual.m b/toolbox/sensors/channel_align_manual.m
index 6f4de79bc1..dfb9186ccd 100644
--- a/toolbox/sensors/channel_align_manual.m
+++ b/toolbox/sensors/channel_align_manual.m
@@ -203,7 +203,7 @@
% ===== DISPLAY HEAD POINTS =====
% Display head points
-figure_3d('ViewHeadPoints', hFig, 1);
+figure_3d('ViewHeadPoints', hFig, 1, 1); % visible and color-coded
% Get patch and vertices
hHeadPointsMarkers = findobj(hFig, 'Tag', 'HeadPointsMarkers');
hHeadPointsLabels = findobj(hFig, 'Tag', 'HeadPointsLabels');
@@ -216,7 +216,9 @@
HeadPointsHpiLoc = [];
if isHeadPoints
% More transparency to view points inside.
- panel_surface('SetSurfaceTransparency', hFig, 1, 0.5);
+ panel_surface('SetSurfaceTransparency', hFig, 1, 0.4);
+ % Hide MEG helmet
+ set(hHelmetPatch, 'visible', 'off');
% Get markers positions
HeadPointsMarkersLoc = get(hHeadPointsMarkers, 'Vertices');
% Hide HeadPoints when looking at EEG and number of EEG channels is the same as headpoints
@@ -244,7 +246,7 @@
% ===== DISPLAY MRI FIDUCIALS =====
% Get the fiducials positions defined in the MRI volume
-sMri = load(file_fullpath(sSubject.Anatomy(sSubject.iAnatomy).FileName), 'SCS');
+sMri = load(file_fullpath(sSubject.Anatomy(sSubject.iAnatomy).FileName), 'SCS', 'History');
if ~isempty(sMri.SCS.NAS) && ~isempty(sMri.SCS.LPA) && ~isempty(sMri.SCS.RPA)
% Convert coordinates MRI => SCS
MriFidLoc = [cs_convert(sMri, 'mri', 'scs', sMri.SCS.NAS ./ 1000); ...
@@ -425,6 +427,8 @@
bst_progress('stop');
end
+% Check and print to command window if previously auto/manual registration, and if MRI fids updated.
+process_adjust_coordinates('CheckPrevAdjustments', in_bst_channel(ChannelFile), sMri);
end
%% ===== MOUSE CALLBACKS =====
@@ -776,6 +780,7 @@ function AlignKeyPress_Callback(hFig, keyEvent)
end
% Ask if needed to update also the other modalities
if isempty(isAll)
+ % TODO We might have < 10 but still want to update electrodes. Verify instead if there are real EEG (not just ECG EOG)
if (gChanAlign.isMeg || gChanAlign.isNirs) && (length(iEeg) > 10)
isAll = java_dialog('confirm', 'Do you want to apply the same transformation to the EEG electrodes ?', 'Align sensors');
elseif ~gChanAlign.isMeg && ~isempty(iMeg)
@@ -890,49 +895,48 @@ function AlignKeyPress_Callback(hFig, keyEvent)
function AlignClose_Callback(varargin)
global gChanAlign;
if gChanAlign.isChanged
- % Sensors or headpoints changed
+ isCancel = false;
strTitle = 'Align sensors';
strType = 'sensor';
- if gChanAlign.isHp
- strTitle = 'Head points';
- strType = 'head points';
- end
-
% Ask user to save changes (only if called as a callback)
if (nargin == 3)
- SaveChanged = 1;
+ SaveChanges = 1;
else
- SaveChanged = java_dialog('confirm', ['The ' strType ' locations changed.' 10 10 ...
- 'Would you like to save changes? ' 10 10], 'Align sensors');
+ [SaveChanges, isCancel] = java_dialog('confirm', ['The ' strType ' locations changed.' 10 10 ...
+ 'Would you like to save changes? ' 10 10], strTitle);
+ end
+ % Don't close figure if cancelled.
+ if isCancel
+ return;
+ end
+ % Get new positions
+ [ChannelMat, Transf, iChannels] = GetCurrentChannelMat();
+ % Load original channel file
+ ChannelMatOrig = in_bst_channel(gChanAlign.ChannelFile);
+ % Report (in command window) max head and sensor displacements from changes.
+ if SaveChanges || gChanAlign.isHeadPoints
+ process_adjust_coordinates('CheckCurrentAdjustments', ChannelMat, ChannelMatOrig);
end
% Save changes to channel file and close figure
- if SaveChanged
+ if SaveChanges
% Progress bar
bst_progress('start', strTitle, 'Updating channel file...');
% Restore standard close callback for 3DViz figures
set(gChanAlign.hFig, 'CloseRequestFcn', gChanAlign.Figure3DCloseRequest_Bak);
drawnow;
- % Get new positions
- [ChannelMat, Transf, iChannels] = GetCurrentChannelMat();
- % Load original channel file
- ChannelMatOrig = in_bst_channel(gChanAlign.ChannelFile);
% Save new electrodes positions in ChannelFile
bst_save(gChanAlign.ChannelFile, ChannelMat, 'v7');
% Get study associated with channel file
[sStudy, iStudy] = bst_get('ChannelFile', gChanAlign.ChannelFile);
% Reload study file
db_reload_studies(iStudy);
+ % Apply to other recordings with same sensor locations in the same subject
+ CopyToOtherFolders(ChannelMatOrig, iStudy, Transf, iChannels);
bst_progress('stop');
end
- else
- SaveChanged = 0;
end
% Only close figure
gChanAlign.Figure3DCloseRequest_Bak(varargin{1:2});
- % Apply to other recordings with same sensor locations in the same subject
- if SaveChanged && ~gChanAlign.isHp
- CopyToOtherFolders(ChannelMatOrig, iStudy, Transf, iChannels);
- end
end
@@ -1439,7 +1443,7 @@ function SelectHeadpoint()
%% ===== DELETE HEADPOINTS =====
function RemoveHeadpoints(varargin)
- global GlobalData gChanAlign;
+ global gChanAlign;
if isempty(gChanAlign.HeadPointsMarkersSel)
% Nothing to do
return
diff --git a/toolbox/sensors/channel_apply_transf.m b/toolbox/sensors/channel_apply_transf.m
index 75f312b628..a7871c45b5 100644
--- a/toolbox/sensors/channel_apply_transf.m
+++ b/toolbox/sensors/channel_apply_transf.m
@@ -47,7 +47,7 @@
if isnumeric(Transf)
R = Transf(1:3,1:3);
T = Transf(1:3,4);
- Transf = @(Loc)(R * Loc + T * ones(1, size(Loc,2)));
+ TransfFunc = @(Loc)(R * Loc + T * ones(1, size(Loc,2)));
else
R = [];
end
@@ -82,7 +82,7 @@
Orient = ChannelMat.Channel(iChan(i)).Orient;
% Update location
if ~isempty(Loc) && ~isequal(Loc, [0;0;0])
- ChannelMat.Channel(iChan(i)).Loc = Transf(Loc);
+ ChannelMat.Channel(iChan(i)).Loc = TransfFunc(Loc);
end
% Update orientation
if ~isempty(Orient) && ~isequal(Orient, [0;0;0])
@@ -95,7 +95,7 @@
end
% If needed: transform the digitized head points
if isHeadPoints && ~isempty(ChannelMat.HeadPoints.Loc)
- ChannelMat.HeadPoints.Loc = Transf(ChannelMat.HeadPoints.Loc);
+ ChannelMat.HeadPoints.Loc = TransfFunc(ChannelMat.HeadPoints.Loc);
end
% If a TransfMeg field with translations/rotations available