diff --git a/toolbox/gui/panel_opticalflow.m b/toolbox/gui/panel_opticalflow.m index 3eb07e7235..5458566363 100644 --- a/toolbox/gui/panel_opticalflow.m +++ b/toolbox/gui/panel_opticalflow.m @@ -58,7 +58,6 @@ % Calculate HHD as well as optical flow jCheckHHD = JCheckBox('Calculate HHD', 0); - jCheckHHD.setEnabled(false); % ======>>>>> RIGHT NOW, HHD IS NOT ENABLED java_setcb(jCheckHHD, 'ActionPerformedCallback', @(h,ev)UpdatePanel()); jPanelSetup.add('br', jCheckHHD); @@ -331,9 +330,15 @@ function Compute(ResultsFile, inputs) %#ok [flowField, int_dF, errorData, errorReg, poincare, timeInterval] = ... bst_call(@bst_opticalflow, F, FV, Time, ... inputs.tStart, inputs.tEnd, inputs.hornSchunck); - if isfield(inputs, 'depthHHD') % ... and, optionally, HHD - % [U A H Vcurl Vdiv index] = ... - % HHD(dataFile, FV, Time, inputs.tStart, inputs.tEnd, inputs.depthHHD); + if inputs.HHDAvailable % ... and, optionally, HHD + [U, A, H, Vcurl, Vdiv, index] = bst_opticalflow_hhd(flowField, FV, inputs.depthHHD); + flowFieldHHD.DivergingPotential = U; + flowFieldHHD.RotatingPotential = A; + flowFieldHHD.HarmonicComponent = H; + flowFieldHHD.RotatingComponent = Vcurl; + flowFieldHHD.DivergingComponent = Vdiv; + else + flowFieldHHD = []; end if inputs.rotate % Rotate flow so that tangent bundle is for circumsphere @@ -343,20 +348,19 @@ function Compute(ResultsFile, inputs) %#ok end % Save results into Results file as latest optical flow calculation - save_flow(ResultsFile, inputs, flowField, flowFieldRotated, ... + save_flow(ResultsFile, inputs, flowField, flowFieldHHD, flowFieldRotated, ... Time, int_dF, errorData, errorReg, poincare); else flowField = ResultsMat.OpticalFlow.flowField; - timeInterval = ResultsMat.OpticalFlow.timeInterval; + timeInterval = linspace(ResultsMat.OpticalFlow.timeInterval(1), ResultsMat.OpticalFlow.timeInterval(2),size(flowField,3)); end % Calculate states if inputs.segment bst_progress('start', 'Optical Flow', 'Segmenting into stable and transition states ...'); - interval = timeInterval(1) : SamplingInterval : timeInterval(2)+2*eps; [stableStates, transientStates, stablePoints, transientPoints, dEnergy] = ... - bst_opticalflow_states(flowField, FV.Faces, FV.Vertices, 3, interval, SamplingInterval, true); + bst_opticalflow_states(flowField, FV.Faces, FV.Vertices, 3, timeInterval, SamplingInterval, true); % Save states Results = in_bst_results(ResultsFile); @@ -400,18 +404,18 @@ function Compute(ResultsFile, inputs) %#ok % flowFieldRotated - Rotate normals-to-vertex to point to sphere, % and then rotate optical flow using same % transformation to get this result - nVertices = size(Vertices,1); - centeredVertices = Vertices - repmat(mean(Vertices), nVertices, 1); + nVertices = size(Vertices,2); + centeredVertices = Vertices - repmat(mean(Vertices,2)', nVertices, 1)'; flowFieldRotated = zeros(size(flowField)); bst_progress('start', 'Optical Flow', ... 'Rotating flows for visualization ... ', 0, nVertices); for m = 1:nVertices - c = VertNormals(m,:); d = -centeredVertices(m,:); + c = VertNormals(:,m); d = -centeredVertices(:,m); current = c/norm(c); desired = d/norm(d); perpCurrent = cross(desired,current)/norm(cross(desired,current)); perpDesired = cross(perpCurrent,desired); - frameChange = [desired' perpDesired' perpCurrent']; + frameChange = [desired perpDesired perpCurrent]; rotation = [dot(current,desired) -dot(current,perpDesired) 0; ... dot(current,perpDesired) dot(current,desired) 0; ... 0 0 1]; @@ -426,7 +430,7 @@ function Compute(ResultsFile, inputs) %#ok end %% ===== SAVE OPTICAL FLOW INTO BRAINSTORM ===== -function save_flow(ResultsFile, inputs, flowField, flowFieldRotated, ... +function save_flow(ResultsFile, inputs, flowField, flowFieldHHD, flowFieldRotated, ... Time, int_dF, errorData, errorReg, poincare) % SAVE_FLOW Save flow results (and possibly publish microstates) % INPUTS: @@ -434,6 +438,12 @@ function save_flow(ResultsFile, inputs, flowField, flowFieldRotated, ... % inputs - inputs to error-check whether results are written % flowField - Optical flow field % dimension (# of vertices) X length(tStart:tEnd) + % flowFieldHHD - flow field Helmholtz Hodge decomposition data + % U - potential field associated to the curl-free component + % A - potential field associated to the div-free component + % H - Harmonic component, suposedly ~0 + % Vcurl - div-free component + % Vdiv - curl-free component % flowFieldRotated - flow field rotated for faux spherical brain % Time - time when activity was reconstructed (including % times for which optical flow is not calculated) @@ -444,31 +454,37 @@ function save_flow(ResultsFile, inputs, flowField, flowFieldRotated, ... opticalFlow.flowField = flowField; % Optical flow results opticalFlow.flowFieldRotatedAvailable = inputs.rotate; + opticalFlow.flowFieldHHDAvailable = inputs.HHDAvailable; + if inputs.HHDAvailable + opticalFlow.depthHHD = inputs.depthHHD; + opticalFlow.flowFieldHHD = flowFieldHHD; + end if inputs.rotate % Results with surface normal rotated towards center of boundary's volume opticalFlow.flowFieldRotated = flowFieldRotated; end - opticalFlow.timeInterval = [inputs.tStart inputs.tEnd]; % Time interval + if isempty(find(Time < inputs.tStart-100*eps, 1, 'last')) % The flow can't be calculated on the first time + opticalFlow.timeInterval = [(inputs.tStart + Time(2)-Time(1)) inputs.tEnd]; + else + opticalFlow.timeInterval = [inputs.tStart inputs.tEnd]; % Time interval + end +% opticalFlow.timeInterval = [inputs.tStart inputs.tEnd]; % Time interval opticalFlow.samplingInterval = Time(2)-Time(1); % Time interval opticalFlow.hornSchunck = inputs.hornSchunck; % Regularization opticalFlow.int_dF = int_dF; opticalFlow.errorData = errorData; % Error in fit to data opticalFlow.errorReg = errorReg; % Error from smooth regularization opticalFlow.poincare = poincare; - opticalFlow.HHDAvailable = inputs.HHDAvailable; - if inputs.HHDAvailable - opticalFlow.depthHHD = inputs.depthHHD; - end % Save optical flow results in original results file Results = in_bst_results(ResultsFile); - if opticalFlow(end).HHDAvailable + if opticalFlow(end).flowFieldHHDAvailable Results = bst_history('add', Results, 'compute', ... - ['Optical flow estimated: [' int2str(inputs.tStart*1000) ... + ['Optical flow & HHD estimated: [' ... + int2str(inputs.tStart*1000) ... ', ' int2str(inputs.tEnd*1000) ']ms']); else Results = bst_history('add', Results, 'compute', ... - ['Optical flow & HHD estimated: [' ... - int2str(inputs.tStart*1000) ... + ['Optical flow estimated: [' int2str(inputs.tStart*1000) ... ', ' int2str(inputs.tEnd*1000) ']ms']); end Results.OpticalFlow = opticalFlow; @@ -495,8 +511,8 @@ function publish_states(ResultsFile, Results, Time, FV) tEndIndex = find(Time < opticalFlow.timeInterval(2)-eps, 1, 'last')+1; % Index of last time point for flow calculation activity = Results.ImageGridAmp(:,tStartIndex:tEndIndex); % Activity in flow-calculated interval extrema = sortrows([ ... - opticalFlow.microstates.stableStates zeros(length(opticalFlow.microstates.stableStates), 1); ... - opticalFlow.microstates.transientStates ones(length(opticalFlow.microstates.transientStates),1) ... + opticalFlow.microstates.stableStates zeros(size(opticalFlow.microstates.stableStates,1), 1); ... + opticalFlow.microstates.transientStates ones(size(opticalFlow.microstates.transientStates, 1),1) ... ]); % Sort extrema for timeline visual % Setup data for visualization @@ -717,7 +733,7 @@ function PlotOpticalFlow(hFig, opticalFlow, currentTime, sSurf) % Process figure (removing old flows if necessary + getting surface axes) nVertices = size(sSurf.Vertices, 1); [ax, currentName] = process_surface(hFig); - + % First check if we need to do anything flagPlotFlow = 0; for n = 1:length(opticalFlow) @@ -768,16 +784,30 @@ function PlotOpticalFlow(hFig, opticalFlow, currentTime, sSurf) % Hold axes to plot on top of surface hold(ax,'on'); if plotRotatedResults - flowField = opticalFlow.flowFieldRotated(:,:,timeIdx); + flowField = opticalFlow.flowFieldRotated; else - flowField = opticalFlow.flowField(:,:,timeIdx); + flowField = opticalFlow.flowField; end - useful = sum(flowField.^2, 2) > max(sum(flowField.^2, 2))*0.1; + + % Leo rescaling parameters + + scalingParameter = mean(sum(flowField.^2, [1 2])); + timeMaxAvg = mean(max(sum(flowField.^2,2),[],1)); + filterParameter = 0.1*timeMaxAvg(1,1); + + flowField = flowField(:,:,timeIdx); + + % scale = sum(flowField.^2,'all')/scalingParameter; + useful = sum(flowField.^2, 2) > filterParameter; % Le filtre s'applique sur chaque image au lieu de s'appliquer sur le tout flowField(~useful,:) = 0; quiver3(ax, ... - Vertices(:,1), Vertices(:,2), Vertices(:,3), ... - flowField(:,1), flowField(:,2), flowField(:,3), ... - 6, 'c', 'LineWidth', 2); % Color is cyan, works well with hot colormap + Vertices(:,1), Vertices(:,2), Vertices(:,3), ... + flowField(:,1), flowField(:,2), flowField(:,3), ... + 6, 'c', 'LineWidth', 2, 'Tag', 'Optical Flow'); % Color is cyan, works well with hot colormap + % quiver3(ax, ... + % Vertices(:,1), Vertices(:,2), Vertices(:,3), ... + % flowField(:,1), flowField(:,2), flowField(:,3), ... + % 6*scale, 'c', 'LineWidth', 2, 'Tag', 'Optical Flow'); % Color is cyan, works well with hot colormap hold(ax,'off'); % Modify figure name if we are in stable/transient state @@ -1105,4 +1135,4 @@ function allow_others(hFig) function CheckboxRotated_Callback(h, ev, hFig, opticalFlow, currentTime, sSurf) PlotOpticalFlow(hFig, opticalFlow, currentTime, sSurf); end -end \ No newline at end of file +end diff --git a/toolbox/math/bst_opticalflow.m b/toolbox/math/bst_opticalflow.m index 22bbbbbcb1..6bad4ce7c6 100644 --- a/toolbox/math/bst_opticalflow.m +++ b/toolbox/math/bst_opticalflow.m @@ -53,12 +53,12 @@ Faces = FV.Faces; Vertices = FV.Vertices; VertNormals = FV.VertNormals; nVertices = size(Vertices,1); % VertNormals = FV.VertNormals'; nFaces = size(Faces,1); -tStartIndex = find(Time < tStart-eps, 1, 'last')+1; % Index of first time point for flow calculation +tStartIndex = find(Time < tStart-100*eps, 1, 'last')+1; % Index of first time point for flow calculation if isempty(tStartIndex) [tmp, tStartIndex] = min(Time); tStartIndex = tStartIndex + 1; end -tEndIndex = find(Time < tEnd-eps, 1, 'last')+1; % Index of last time point for flow calculation +tEndIndex = find(Time < tEnd-100*eps, 1, 'last')+1; % Index of last time point for flow calculation if isempty(tEndIndex) [tmp, tEndIndex] = max(Time); tEndIndex = tEndIndex + 1; diff --git a/toolbox/math/bst_opticalflow_hhd.m b/toolbox/math/bst_opticalflow_hhd.m new file mode 100644 index 0000000000..76efe8d722 --- /dev/null +++ b/toolbox/math/bst_opticalflow_hhd.m @@ -0,0 +1,395 @@ +function [U, A, H, Vcurl, Vdiv, index] = bst_opticalflow_hhd(opticalFlow, FV, depthHHD) + +% BST_OPTICALFLOW_HHD: Computes the Helmhlotz-Hodge Decomposition of the +% optical flow of MEG/EEG activities on the cortical surface. +% +% USAGE: [U, A, H, Vcurl, Vdiv, index] = +% bst_opticalflow_hhd(opticalFlow, FV, depthHHD) +% +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c)2000-2020 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: Léo Nouvelle, 2020 +% Sheraz Khan, 2008 +% Julien Lefevre, 2006-2010 +% +% INPUTS +% opticalFlow - Optical Flow motion field +% FV - Tesselation for calculating spatial derivatives +% depthHHD - Recursive depth of the Helmholtz-Hodge decomposition +% OUTPUTS +% U - Potential Field associated to the curl-free +% component of the HH decomposition +% A - Potential Field associated to the divergence-free +% component of the HH decomposition +% H - Harmonic component of the HH decomposition +% Vcurl - Divergence-free component of the HHD +% Vdiv - Curl-free component of the HHD + +index = 1:size(opticalFlow,3); +nVertices = size(FV.Vertices,1); + +U=zeros(nVertices,length(index)); +A=zeros(nVertices,length(index)); +Vcurl=zeros(nVertices,3,length(index)); +Vdiv=zeros(nVertices,3,length(index)); +H=zeros(nVertices,3,length(index)); + +% Calculating geometric features of the tessalation +[cont_grad_v,cont_grad_vb,grad_v,aires,norm_tri,VertFaceConn,Pn1]=geometry_tesselation(FV.Faces,FV.Vertices,3); + +faces=FV.Faces; +vertices=FV.Vertices; + +bst_progress('start', 'Helmholtz Hodge Decomposition', ... + 'Computing Helmhlotz Hodge Decomposition ...', 0, length(index)); + +%% Désactivé en l'absence de Parralel Computing Toolbox + +% parpool +% try +% poolobj = gcp('nocreate'); +% if isempty(poolobj) +% parpool; +% end +% catch +% disp (' ') +% end +%parfor ii=index(1):index(end) +%Toolbox + +% Ui=zeros(length(vertices),length(index),recur); +% Ai=zeros(length(vertices),length(index),recur); +% Vcurli=zeros(length(vertices),3,length(index),recur); +% Vdivi=zeros(length(vertices),3,length(index),recur); +% Hi=zeros(length(vertices),3,length(index),recur); + +%% + +for ii = index + Ui=zeros(length(vertices),depthHHD); + Ai=zeros(length(vertices),depthHHD); + Vcurli=zeros(length(vertices),3,depthHHD); + Vdivi=zeros(length(vertices),3,depthHHD); + Hi=zeros(length(vertices),3,depthHHD); + Ui=zeros(length(vertices),depthHHD); + for i=1:depthHHD + if i==1 + [Ui(:,i), Ai(:,i), Hi(:,:,i), Vcurli(:,:,i), Vdivi(:,:,i)]=hhdr(opticalFlow(:,:,ii),1,0,cont_grad_v,cont_grad_vb,grad_v,aires,norm_tri,VertFaceConn,Pn1,faces,vertices); + else + [Ui(:,i), Ai(:,i), Hi(:,:,i), Vcurli(:,:,i), Vdivi(:,:,i)]=hhdr(Hi(:,:,i-1),1,0,cont_grad_v,cont_grad_vb,grad_v,aires,norm_tri,VertFaceConn,Pn1,faces,vertices); + end + end %recur + U(:,ii)=sum(Ui,2); + A(:,ii)=sum(Ai,2); + Vcurl(:,:,ii)=sum(Vcurli,3); + Vdiv(:,:,ii)=sum(Vdivi,3); + H(:,:,ii)=Hi(:,:,depthHHD); + + bst_progress('inc', 1); % Update progress bar +end % index +% matlabpool close +bst_progress('stop'); + + +end + +% ========================================================================= +% RECURSIVE FUNCTION +% ========================================================================= + +function [U, A, H, Vcurl, Vdiv, index]=hhdr(V,range,verbose,cont_grad_v,cont_grad_vb,grad_v,aires,norm_tri,VertFaceConn,Pn1,faces,vertices) + +%Routine to calculate U, V and H component of the vector field +% Jflowc contain all the information including optical flow +% range contain Time period on which decomposition required e.g. +% range=500:700 +% U,V,H are the component of HHD of a vector field + + +%% Extracting information from Jflowc and initialization +% Time=Jflowc.t; +FV.Vertices=vertices; +FV.Faces=faces; +clear faces vertices +if nargin<2 +index=1:length(Time); +else +index=range; +end + +U=zeros(length(FV.Vertices),length(index)); +A=zeros(length(FV.Vertices),length(index)); +Vcurl=zeros(length(FV.Vertices),3,length(index)); +Vdiv=zeros(length(FV.Vertices),3,length(index)); +H=zeros(length(FV.Vertices),3,length(index)); + +% calculating U, V and H + +if verbose==1 +h = waitbar(0,'Computing U,V and H'); +end + +for ii=index + if verbose==1 +waitbar((ii-(index(1)-1))/length(index)); + end +FieldV=V(:,:,ii); + +B=zeros(length(FV.Vertices),1); +B2=zeros(length(FV.Vertices),1); + +for i=1:size(FV.Faces,1) + % projection of flows.V(i,:,t) on the triangle i + Pn=eye(3)-(norm_tri(i,:)'*norm_tri(i,:)); + nodes=FV.Faces(i,:); + VV=mean(FieldV(nodes,:),1); + projectV=VV*Pn; + for s=1:3 + B(nodes(s),1)=B(nodes(s),1)+sum(projectV.*grad_v{s}(i,:)*aires(i),2); + B2(nodes(s),1)=B2(nodes(s),1)+sum(projectV.*cross(grad_v{s}(i,:),norm_tri(i,:))*aires(i),2); + end +end + +% Léo 02/2020 : Conditioning phase + +[P,R,C] = equilibrate(cont_grad_v); +[P2,R2,C2] = equilibrate(cont_grad_vb); + +K = P*R*cont_grad_v*C; +Kb = P2*R2*cont_grad_vb*C2; +[L1,U1] = ilu(K,struct('type','ilutp','droptol',1e-2,'thresh',0)); +[L2,U2] = ilu(Kb,struct('type','ilutp','droptol',1e-2,'thresh',0)); +Y = P*R*B; +Y2 = P2*R2*B2; + +% Linear system inversion + +U(:,ii-(index(1)-1))= C*lsqr(K,Y,1e-20,10000, L1, U1); +% % error1 = norm(cont_grad_v*U(:,ii-(index(1)-1))-B)/norm(B); +% disp(append('actual error is ', int2str(error1))) +A(:,ii-(index(1)-1))= C2*lsqr(Kb,Y2,1e-20,10000, L2, U2); +% error2 = norm(cont_grad_vb*A(:,ii-(index(1)-1))-B2)/norm(B2); +% disp(append('actual error is ', int2str(error2))) + +Vd=curl(U(:,ii-(index(1)-1)),FV,grad_v); +Vc=curl(A(:,ii-(index(1)-1)),FV,grad_v,norm_tri); + +Vcurl(:,:,ii-(index(1)-1))=tri2vert(Vc,VertFaceConn,Pn1); +Vdiv(:,:,ii-(index(1)-1))=tri2vert(Vd,VertFaceConn,Pn1); + +H(:,:,ii-(index(1)-1))=V(:,:,ii)-Vcurl(:,:,ii-(index(1)-1))-Vdiv(:,:,ii-(index(1)-1)); + +end % for +if verbose==1 +delete(h) +end + + +end + +% ========================================================================= +% EXTERNAL FUNCTIONS +% ========================================================================= + +function [cont_grad_v,cont_grad_vb,grad_v,aires,norm_tri,VertFaceConn,Pn1]=geometry_tesselation(tri,coord,dim) + +% Computation of the regularizing part in the variationnal approach (SS grad(v_k)grad(v_k')) and +% other geometrical quantities. +% INPUTS : +% tri : triangles of the tesselation +% coord : coordinates of the tesselation +% dim : 3, scalp or cortical surface, 2 flat surface +%%%%% +% OUTPUTS : +% cont_grad_v, cont_grad_vb : regularizing matrices +% grad_v : gradient of the basis functions in Finite Element Methods +% aires : area of the triangles +% norm_tri : normal of each triangle +% Pn1 : Projector matrix on surface's triangular meshes + +nbr_capt=size(coord,1); % Number of nodes + +%% Geometric quantities + +[grad_v,aires,norm_tri]=carac_tri(tri,coord,dim); +norm_coord=carac_coord(tri,coord,norm_tri); + +% Regularizing matrix SS grad(v_k)grad(v_k') + +index1=[]; +index2=[]; +termes_diag=[]; +termes_diag_b=[]; +tang_scal_11=[]; +tang_scal_11_b=[]; + +for k=1:3 + for j=k+1:3 + index1=[index1,tri(:,k)]; + index2=[index2,tri(:,j)]; + tang_scal_11=[tang_scal_11,sum(grad_v{k}.*grad_v{j},2).*aires]; + tang_scal_11_b=[tang_scal_11_b,sum(cross(grad_v{k},norm_tri).*cross(grad_v{j},norm_tri),2).*aires]; + end + termes_diag=[termes_diag,sum(grad_v{k}.^2,2).*aires]; + termes_diag_b=[termes_diag_b,sum(cross(grad_v{k},norm_tri).^2,2).*aires]; +end + +D=sparse(tri,tri,termes_diag,nbr_capt,nbr_capt); +Db=sparse(tri,tri,termes_diag_b,nbr_capt,nbr_capt); + +E11=sparse(index1,index2,tang_scal_11,nbr_capt,nbr_capt); +E11=E11+E11'+D; +cont_grad_v=E11; + +E11b=sparse(index1,index2,tang_scal_11_b,nbr_capt,nbr_capt); +E11b=E11b+E11b'+Db; +cont_grad_vb=E11b; + +% Vertice to triangular mesh connectivity +VertFaceConn=cell(size(coord,1),1); +for tt=1:size(tri,1) + for ind=tri(tt,:) + VertFaceConn{ind,1}=[VertFaceConn{ind,1};tt]; + end +end + +Pn1 = zeros(3,3,size(coord,1)); + +% Projection matrices on triangular mesh +for ii=1:size(coord,1) + Pn1(:,:,ii)=eye(3)-(norm_coord(ii,:)'*norm_coord(ii,:)); +end + +end + +function [grad,aires,vectoriel]=carac_tri(tri,coord,dim) +% Computes some geometric quantities from a surface +% INPUTS +% tri : triangles of tesselation +% coord : coordinates of nodes +% dim : deefault 3, 2 for projection on a plane +% +% OUTPUTS +% grad : gradient of basis function (Finite Elements Method) on each triangle +% aires : area of each triangle +% vectoriel : normal of each triangle + +% Edges of each triangles +u=coord(tri(:,2),:)-coord(tri(:,1),:); +v=coord(tri(:,3),:)-coord(tri(:,2),:); +w=coord(tri(:,1),:)-coord(tri(:,3),:); + +% Length of each edges and angles bewteen edges +uu=sum(u.^2,2); +vv=sum(v.^2,2); +ww=sum(w.^2,2); +uv=sum(u.*v,2); +vw=sum(v.*w,2); +wu=sum(w.*u,2); + +% 3 heights of each triangle and their norm +h1=w-((vw./vv)*ones(1,dim)).*v; +h2=u-((wu./ww)*ones(1,dim)).*w; +h3=v-((uv./uu)*ones(1,dim)).*u; +hh1=sum(h1.^2,2); +hh2=sum(h2.^2,2); +hh3=sum(h3.^2,2); + +% Gradient of the 3 basis functions on a triangle +grad=cell(1,dim); +grad{1}=h1./(hh1*ones(1,dim)); +grad{2}=h2./(hh2*ones(1,dim)); +grad{3}=h3./(hh3*ones(1,dim)); + +% Prevents from pathological gradients + +indices1=find(sum(grad{1}.^2,2)==0|isnan(sum(grad{1}.^2,2))); +indices2=find(sum(grad{2}.^2,2)==0|isnan(sum(grad{2}.^2,2))); +indices3=find(sum(grad{3}.^2,2)==0|isnan(sum(grad{3}.^2,2))); +indices21=find(sum(grad{1}.^2,2)); +indices22=find(sum(grad{2}.^2,2)); +indices23=find(sum(grad{3}.^2,2)); + +min_norm_grad=min([sum(grad{1}(indices21,:).^2,2);sum(grad{2}(indices22,:).^2,2);sum(grad{3}(indices23,:).^2,2)]); + +grad{1}(indices1,:)=repmat([1 1 1]/min_norm_grad,length(indices1),1); +grad{2}(indices2,:)=repmat([1 1 1]/min_norm_grad,length(indices2),1); +grad{3}(indices3,:)=repmat([1 1 1]/min_norm_grad,length(indices3),1); + + +% Area of triangles and normals of each triangle +aires=sqrt(hh1.*vv)/2; +indices=isnan(aires); +aires(indices)=0; + +if dim==3 + vectoriel=cross(w,u); + vectoriel=vectoriel./repmat(sqrt(sum(vectoriel.^2,2)),1,3); +else + vectoriel=[]; +end +end + +function norm_coord=carac_coord(tri,coord,norm_tri) +% Normal at each node of the tesselation +% INPUTS : +% tri : triangles +% coord : coordinates of each node +% norm_tri : normal of each triangle +norm_coord=zeros(size(coord,1),3); +for i=1:size(tri,1) + for k=1:3 + norm_coord(tri(i,k),:)=norm_coord(tri(i,k),:)+norm_tri(i,:); + end +end + +norm_coord=norm_coord./repmat(sqrt(sum(norm_coord.^2,2)),1,3); + +% Not very satisfying solution to the problem of pathological anatomy + +indices=isnan(norm_coord); +norm_coord(indices)=0; +end + +function [V]=curl(X,FV,grad_v,normals) +% Builds curl-free and divergence-free components from the associated +% potential fields + +grad_X=repmat(X(FV.Faces(:,1)),1,3).*grad_v{1}+repmat(X(FV.Faces(:,2)),1,3).*grad_v{2}+repmat(X(FV.Faces(:,3)),1,3).*grad_v{3}; +if nargin<4 %just compute the gradient of the scalar field X + V=grad_X; +else %curl of the scalar field X + V=cross(grad_X,normals); +end +end + +function Vv =tri2vert(Vt,VertFaceConn,Pn) +% Transfers the motion fields calculated on the triangular meshes on the +% vertices + +Vv = zeros(size(VertFaceConn,1),3); +card = zeros(1,size(VertFaceConn,1)); + +for ii=1:size(VertFaceConn,1) + Vv(ii,:)=(sum(Vt(VertFaceConn{ii},:),1))*Pn(:,:,ii); %average + projection + card(ii)=length(VertFaceConn{ii}); +end + +Vv=Vv./repmat(card',1,3); +end \ No newline at end of file diff --git a/toolbox/math/bst_opticalflow_states.m b/toolbox/math/bst_opticalflow_states.m index 2713942737..5dfe88fcb1 100644 --- a/toolbox/math/bst_opticalflow_states.m +++ b/toolbox/math/bst_opticalflow_states.m @@ -44,116 +44,171 @@ % Setup: get displacement energy [syed, triangleAreas] = geometry_tesselation(Faces, Vertices, dimension); dEnergy = zeros(1, size(flowField,3)); -for m = 1:size(flowField,3); - v12 = sum((flowField(Faces(:,1),:,m)+flowField(Faces(:,2),:,m)).^2,2) / 4; - v23 = sum((flowField(Faces(:,2),:,m)+flowField(Faces(:,3),:,m)).^2,2) / 4; - v13 = sum((flowField(Faces(:,1),:,m)+flowField(Faces(:,3),:,m)).^2,2) / 4; - dEnergy(m) = sum(triangleAreas.*(v12+v23+v13)); +% dEnergy2 = zeros(1, size(flowField,3)); +for m = 1:size(flowField,3) +% v12 = sum((flowField(Faces(:,1),:,m)+flowField(Faces(:,2),:,m)).^2,2) / 4; +% v23 = sum((flowField(Faces(:,2),:,m)+flowField(Faces(:,3),:,m)).^2,2) / 4; +% v13 = sum((flowField(Faces(:,1),:,m)+flowField(Faces(:,3),:,m)).^2,2) / 4; +% dEnergy(m) = sum(triangleAreas.*(v12+v23+v13)); + + % Version Leo, 16/03/2020 + dEnergy(m) = sum(triangleAreas.*sum((flowField(Faces(:,1),:,m)+flowField(Faces(:,2),:,m)+flowField(Faces(:,3),:,m)).^2,2)); end dEnergy = sqrt(dEnergy); % Get square root for easier visualization % Find local minima -minima = local_extrema(dEnergy, false, 1, length(dEnergy)); -minima = sortrows([dEnergy(minima); minima]')'; minima = minima(2,:); +[s, minima] = local_extrema(dEnergy, false, 1, length(dEnergy)); +minima = sortrows([s(minima); minima]')'; minima = minima(2,:); % Find local maxima -maxima = local_extrema(dEnergy, true, 1, length(dEnergy)); -maxima = sortrows([-dEnergy(maxima); maxima]')'; maxima = maxima(2,:); +[s2, maxima] = local_extrema(dEnergy, true, 1, length(dEnergy)); +maxima = sortrows([-s2(maxima); maxima]')'; maxima = maxima(2,:); % Display displacement energy and locations of maxima/minima if displayFlag - figure; hEnergy = axes; - plot(hEnergy, interval, dEnergy, 'Color', [0 0 0]) - axis([interval(1) interval(end) 0 max(dEnergy)]); + figure; + hEnergy = axes; + plot(hEnergy, interval, s, 'Color', [0 0 0]) + axis([interval(1) interval(end) 0 max(s)]); hold(hEnergy, 'on'); - plot(hEnergy, interval(minima), dEnergy(minima), 'g*', ... - interval(maxima), dEnergy(maxima), 'r*'); + plot(hEnergy, interval(minima), s(minima), 'g*', ... + interval(maxima), s(maxima), 'r*'); hold(hEnergy, 'off'); xlabel('Time (ms)') - ylabel('\partial Energy') + ylabel('\partial Energy filtered') end % Find transient and stable (defined as not-transient) states -transient = transient_states(dEnergy, maxima, minima); +% Leo 19/03/2020 +transient = transient_states(s, maxima, minima); transient = sortrows(transient); -stable = stable_states(dEnergy, maxima, minima, transient); +stable = stable_states(s, maxima, minima, transient); stable = sortrows(stable); % Get list of transient and stable state intervals +% Il faut revoir clean_flow_states. extrema = clean_flow_states(transient, stable, samplingInterval); -transient = extrema(extrema(:,3) > 1-eps, 1:2); transientPoints = []; -for m = 1:length(transient) +transient = extrema(extrema(:,4) > 1-eps, 1:2); transientPoints = []; +for m = 1:size(transient,1) transientPoints = [transientPoints transient(m,1):transient(m,2)]; end -stable = extrema(extrema(:,3) < eps, 1:2); stablePoints = []; -for m = 1:length(stable) +stable = extrema(extrema(:,4) < eps, 1:2); stablePoints = []; +for m = 1:size(stable,1) stablePoints = [stablePoints stable(m,1):stable(m,2)]; end % Display microstate interval on displaced energy curve if displayFlag for m = 1:size(extrema,1) - if abs(extrema(m,3)) < eps + if abs(extrema(m,4)) < eps type = 'g'; - elseif extrema(m,3) > 1-eps + elseif extrema(m,4) > 1-eps type = 'r'; end hold(hEnergy, 'on'); area(hEnergy, interval(extrema(m,1):extrema(m,2)), ... - dEnergy(extrema(m,1):extrema(m,2)), 'FaceColor', type); + s(extrema(m,1):extrema(m,2)), 'FaceColor', type); hold(hEnergy, 'off'); end end +% Uncomment following lines to display the unfilterd/original figure. +% if displayFlag +% figure; +% h2Energy = axes; +% plot(h2Energy, interval, dEnergy, 'Color', [0 0 0]) +% axis([interval(1) interval(end) 0 max(dEnergy)]); +% hold(h2Energy, 'on'); +% plot(h2Energy, interval(minima), dEnergy(minima), 'g*', ... +% interval(maxima), dEnergy(maxima), 'r*'); +% hold(h2Energy, 'off'); +% xlabel('Time (ms)') +% ylabel('\partial Energy') +% end +% +% if displayFlag +% for m = 1:size(extrema,1) +% if abs(extrema(m,4)) < eps +% type = 'g'; +% elseif extrema(m,4) > 1-eps +% type = 'r'; +% end +% +% hold(h2Energy, 'on'); +% area(h2Energy, interval(extrema(m,1):extrema(m,2)), ... +% dEnergy(extrema(m,1):extrema(m,2)), 'FaceColor', type); +% hold(h2Energy, 'off'); +% end +% end end %% ===== CLEAN SMALL INTERVALS ===== function extrema = clean_flow_states(transient, stable, samplingInterval) -minInterval = floor(0.002/samplingInterval); +% Swallow intervals that are short of length and small relative distance +% with those of longer length + +% Short interval detection size + +minInterval = floor(3*0.0017/samplingInterval); if minInterval < 1 minInterval = 1; end -% Swallow intervals that are of short length into intervals of longer length +% Relative distance threshold + +maximum = max(transient(:,3)); +minimum = min(transient(:,3)); +threshold = 0.075*(maximum-minimum); + +% Swallow intervals + extrema = sortrows([transient ones(size(transient,1), 1); ... stable zeros(size(stable,1), 1)]); + for m = 2:size(extrema,1) - if extrema(m,2) < extrema(m,1) + minInterval - if extrema(m-1,2) > extrema(m,1) - minInterval - if extrema(m,2)-extrema(m,1) > extrema(m-1,2)-extrema(m-1,1) - tag = extrema(m,3); - else - tag = extrema(m-1,3); - end - extrema(m-1,2) = extrema(m,2); % 2nd interval's beginning is 1st interval - extrema(m,1) = extrema(m-1,1); % 1st interval's end is 2nd interval - extrema(m-1,3) = tag; % Both get same label ... - extrema(m,3) = tag; % ... of the bigger interval - elseif m < size(extrema,1) && extrema(m,2) > extrema(m+1,1) - minInterval - if extrema(m,2)-extrema(m,1) > extrema(m+1,2)-extrema(m+1,1) - tag = extrema(m,3); - else - tag = extrema(m+1,3); - end - extrema(m,2) = extrema(m+1,2); % 2nd interval's beginning is 1st interval - extrema(m+1,1) = extrema(m,1); % 1st interval's end is 2nd interval - extrema(m,3) = tag; % Both get same label ... - extrema(m+1,3) = tag; % ... of the bigger interval - else - extrema(m,1) = max(extrema(m,1) - minInterval, 1); % 2nd interval's beginning is 1st interval - extrema(m,2) = min(extrema(m,2) + minInterval, length(dEnergy)); % 1st interval's end is 2nd interval + if extrema(m,2) < extrema(m,1) + minInterval % Small sized interval + % We test the previous interval is close and the relative difference is small + if (extrema(m-1,2) > extrema(m,1) - minInterval) && (abs(extrema(m,3)-extrema(m-1,3)) extrema(m+1,1) - minInterval ... + && (abs(extrema(m,3)-extrema(m-1,3))= extrema(m-1,2)-extrema(m-1,1) + tag = extrema(m,4); + else + tag = extrema(m-1,4); + end + extrema(m-1,2) = extrema(m,2); % 2nd interval's beginning is 1st interval + extrema(m,1) = extrema(m-1,1); % 1st interval's end is 2nd interval + extrema(m-1,4) = tag; % Both get same label ... + extrema(m,4) = tag; % ... of the bigger interval + % We don't actualize the local extremas because we won't be + % using them after this function + end end end end % Merge consecutive states of the same type for m = 2:size(extrema,1) - if abs(extrema(m,3)-extrema(m-1,3)) < eps + if abs(extrema(m,4)-extrema(m-1,4)) < eps n = m+1; - while n <= size(extrema,1) && abs(extrema(n,3)-extrema(n-1,3)) < eps + while n <= size(extrema,1) && abs(extrema(n,4)-extrema(n-1,4)) < eps n = n+1; end extrema((m-1):(n-2),2) = extrema(n-1,2); % 2nd interval's beginning is 1st interval @@ -216,7 +271,7 @@ end % Label microstate - transient = [transient; begin stop]; + transient = [transient; begin stop value]; end end @@ -286,14 +341,14 @@ % Label microstate if stop-begin >= 1 - stable = [stable; begin stop]; + stable = [stable; begin stop -value]; end end end %% ===== LOCAL EXTREMA OF A CURVE ===== -function extrema = local_extrema(signal, maxOrMin, tStart, tEnd) +function [s, extrema] = local_extrema(signal, maxOrMin, tStart, tEnd) % LOCAL_MINIMA Find locations of local extrema in signal % INPUTS: % signal - curve for which we find maxima OR minima @@ -304,25 +359,35 @@ % OUTPUTS: % extrema - Locations of local extrema -% Preprocessing: smoothing of signal for initial extrema +% Preprocessing: smoothing of signal for initial extrema. We will use a +% binomial convolution filter (Gaussian approximation) of size 12 +% (calculated for filtering the signal at 30Hz). We use the mirror +% technique for the signal at the extremities. + +%Leo 19/03/2020 + if maxOrMin % Local maxima == local minima of negative signal = -1*signal; end -valleys = []; smoothingFilter = [1 1 1]/3; s = signal; -s = conv(s, smoothingFilter, 'same'); -s(tStart) = signal(tStart); s(tEnd) = signal(tEnd); -s = conv(s, smoothingFilter, 'same'); -s(tStart) = signal(tStart); s(tEnd) = signal(tEnd); -s = conv(s, smoothingFilter, 'same'); -s(tStart) = signal(tStart); s(tEnd) = signal(tEnd); -s = conv(s, smoothingFilter, 'same'); -s(tStart) = signal(tStart); s(tEnd) = signal(tEnd); +valleys = []; +% Création of the convolution filter +smoothingFilter = (1); +for i = 1:6 + smoother = (1/4)*[1 2 1]; + smoothingFilter = conv(smoothingFilter, smoother); +end +% On créé le signal 'miroir sur lequel on va appliquer le filtre +x = size(signal,2); +s = zeros(1, x+2*6); +s(7:6+x) = signal; +for i = 1:6 + s(i) = signal(7-i); + s(size(s,2)+1-i) = signal(x-6+i); +end s = conv(s, smoothingFilter, 'same'); -s(tStart) = signal(tStart); s(tEnd) = signal(tEnd); - -% Find local minima (of flipped signal if maxima desired) +s = s(7:6+x); -for t = (tStart+2):(tEnd-1) +for t = (tStart+1):(tEnd-1) if s(t) < s(t-1) if s(t) < s(t+1) % Definite dip valleys(end+1) = t; @@ -342,7 +407,7 @@ for m = 1:length(valleys) jitter = valleys(m) + (-2:2); jitter = jitter(1 <= jitter & jitter <= length(signal)); - [syed, idx] = min(signal(jitter)); + [syed, idx] = min(s(jitter)); extrema(m) = idx + jitter(1) - 1; end @@ -350,6 +415,8 @@ extrema(extrema <= tStart | extrema >= tEnd) = []; extrema = unique(extrema); +% plot(hEnergy, extrema, s(extrema), 'g*'); + end %% ===== TESSELATION NORMALS ===== @@ -421,24 +488,4 @@ FaceNormals = []; end -% % Calculate normals to surface at each vertex (from normals at each face) -% VertNormals = zeros(size(Vertices,1),3); -% bst_progress('start', 'Optical Flow', ... -% 'Computing normals to surface at every vertex ...', 1, size(Faces,1)); -% for facesIdx=1:size(Faces,1); -% VertNormals(Faces(facesIdx,:),:) = ... -% VertNormals(Faces(facesIdx,:),:) + ... -% repmat(FaceNormals(facesIdx,:), [3 1]); -% -% if mod(facesIdx,20) == 0 -% bst_progress('inc', 20); % Update progress bar -% end -% end -% bst_progress('stop'); -% -% % Normalize perpendicular-to-surface vectors for each vertex -% VertNormals = VertNormals ./ ... -% repmat(sqrt(sum(VertNormals.^2,2)),1,3); -% VertNormals(isnan(VertNormals)) = 0; % For pathological anatomy - end \ No newline at end of file