text:
function [flowField, int_dF, errorData, errorReg, poincare, interval] = ...
  bst_opticalflow(F, FV, Time, tStart, tEnd, hornSchunck)
% BST_OPTICALFLOW: Computes the optical flow of MEG/EEG activities on the cortical surface.
%
% USAGE:  [flow, dEnergy, int_dI, errorData, errorReg, index] = ...
%   bst_opticalflow(dataFile, channelFile, tStart, tEnd, hornSchunck)
%
% @=============================================================================
% 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: Julien Lef�vre, 2006-2010
%          Syed Ashrafulla, 2010

% INPUTS
%   F                 - Reconstructed sources
%   FV                - Tesselation for calculating spatial derivatives
%   Time              - Time points of sampling of recordings
%   tStart            - First time point for optical flow analysis
%   tEnd              - Last time point for optical flow analysis
%   hornSchunck       - Parameter to tune optical flow calculation
%   segment           - Flag to split data into stable and transient states
% OUTPUTS
%   flowField         - Optical flow field
%                       dimension (# of vertices) X length(tStart:tEnd)
%   int_dF            - Constant term in variational formulation
%   errorData         - Error in fit to data
%   errorReg          - Energy in regularization
%   poincare          - Poincar� index

%/---Script Authors---------------------\
%|                                      |
%|   *** J.Lef�vre, PhD                 |
%|   julien.lefevre@chups.jussieu.fr    |
%|                                      |
%\--------------------------------------/

dimension = 3; % 2 for projected maps
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
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
if isempty(tEndIndex)
    [tmp, tEndIndex] = max(Time);
    tEndIndex = tEndIndex + 1;
end
if tEndIndex > size(Time, 2)
    tEndIndex = size(Time, 2);
end
interval = Time(tStartIndex:tEndIndex); % Interval of flow calculations
intervalLength = tEndIndex-tStartIndex+1; % Length of time interval for calculations
M = max(max(abs(F(:,tStartIndex-1:tEndIndex)))); F = F/M; % Scale values to avoid precision error
flowField = zeros(nVertices, dimension, intervalLength);
dEnergy = zeros(1, intervalLength);
errorData = zeros(1, intervalLength);
errorReg = zeros(1, intervalLength);
int_dF = zeros(1, intervalLength);

[regularizerOld, gradientBasis, tangentPlaneBasisCell, tangentPlaneBasis, ...
  triangleAreas, FaceNormals, sparseIndex1, sparseIndex2]= ...
  regularizing_matrix(Faces, Vertices, VertNormals, dimension); % 2 for projected maps
regularizer = hornSchunck * regularizerOld;

% Projection of flow on triangle (for Poincar� index)
Pn = zeros(3,3,nFaces);
for facesIdx = 1:nFaces
    Pn(:,:,facesIdx) = eye(3) - ...
      (FaceNormals(facesIdx,:)'*FaceNormals(facesIdx,:));
end
poincare = zeros(nFaces, intervalLength);

% Optical flow calculations
bst_progress('start', 'Optical Flow', ...
  'Computing optical flow ...', tStartIndex-1, tEndIndex);
for timePoint = tStartIndex:tEndIndex
  timeIdx = timePoint-tStartIndex+1;

  % Solve for flow
  dF = F(:,timePoint)-F(:,timePoint-1);
  [dataFit, B] = data_fit_matrix(Faces, nVertices, ...
    gradientBasis, triangleAreas, tangentPlaneBasisCell, ...
    F(:,timePoint), dF, dimension, ...
    sparseIndex1, sparseIndex2);
  X = (dataFit + regularizer) \ B;

  % Save flows correctly (for 3D surface, have to project back to 3D)
  if dimension == 3 % X coordinates are in tangent space
      flowField(:,:,timeIdx) = ...
        repmat(X(1:end/2), [1,3]) .* tangentPlaneBasis(:,:,1) ...
        + repmat(X(end/2+1:end), [1,3]) .* tangentPlaneBasis(:,:,2);
  else % X coordinates are in R^2
      flowField(:,:,timeIdx) = [X zeros(nVertices,1)];
  end

  errorData(timeIdx)= X'*dataFit*X - 2*B'*X; % Data fit error
  errorReg(timeIdx)= X'*regularizer*X; % Regularization error

  % Variational formulation constant
  dF12=(dF(Faces(:,1),:)+dF(Faces(:,2))).^2;
  dF23=(dF(Faces(:,2),:)+dF(Faces(:,3))).^2;
  dF13=(dF(Faces(:,1),:)+dF(Faces(:,3))).^2;
  int_dF(timeIdx) = sum(triangleAreas.*(dF12+dF23+dF13)) / 24;

  % Precompute flowfield with faces to save time in the loop.
  FacesFlowField = reshape(flowField(Faces', :, timeIdx)', [3,3,nFaces]);

  % Poincare Index of each triangle
  for facesIdx=1:nFaces
      poincare(facesIdx,timeIdx) = ... % projection of flowField(f,:,t) on triangle f
          poincare_index(Pn(:,:,facesIdx) * FacesFlowField(:,:,facesIdx));
  end

  % Displacement energy
  v12 = sum((flowField(Faces(:,1),:,timeIdx)+flowField(Faces(:,2),:,timeIdx)).^2,2) / 4;
  v23 = sum((flowField(Faces(:,2),:,timeIdx)+flowField(Faces(:,3),:,timeIdx)).^2,2) / 4;
  v13 = sum((flowField(Faces(:,1),:,timeIdx)+flowField(Faces(:,3),:,timeIdx)).^2,2) / 4;
  dEnergy(timeIdx) = sum(triangleAreas.*(v12+v23+v13));

  bst_progress('inc', 1); % Update progress bar
end

bst_progress('stop');

end

% =========================================================================
% =====  EXTERNAL FUNCTIONS ===============================================
% =========================================================================

%% ===== TESSELATION NORMALS =====
function [gradientBasis, triangleAreas, FaceNormals] = ...
  geometry_tesselation(Faces, Vertices, dimension)
% GEOMETRY_TESSELATION    Computes some geometric quantities from a surface
%
% INPUTS:
%   Faces           - triangles of tesselation
%   Vertices        - coordinates of nodes
%   dimension       - 3 for scalp or cortical surface (default)
%                     2 for plane (channel cap, etc)
% OUTPUTS:
%   gradientBasis   - gradient of basis function (FEM) on each triangle
%   triangleAreas 	- area of each triangle
%   FaceNormals    - normal of each triangle

% Edges of each triangles
u = Vertices(Faces(:,2),:)-Vertices(Faces(:,1),:);
v = Vertices(Faces(:,3),:)-Vertices(Faces(:,2),:);
w = Vertices(Faces(:,1),:)-Vertices(Faces(:,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,dimension)).*v;
h2 = u-((wu./ww)*ones(1,dimension)).*w;
h3 = v-((uv./uu)*ones(1,dimension)).*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
gradientBasis = cell(1,dimension);
gradientBasis{1} = h1./(hh1*ones(1,dimension));
gradientBasis{2} = h2./(hh2*ones(1,dimension));
gradientBasis{3} = h3./(hh3*ones(1,dimension));

% Remove pathological gradients
indices1 = find(sum(gradientBasis{1}.^2,2)==0|isnan(sum(gradientBasis{1}.^2,2)));
indices2 = find(sum(gradientBasis{2}.^2,2)==0|isnan(sum(gradientBasis{2}.^2,2)));
indices3 = find(sum(gradientBasis{3}.^2,2)==0|isnan(sum(gradientBasis{3}.^2,2)));

min_norm_grad = min([ ...
  sum(gradientBasis{1}(sum(gradientBasis{1}.^2,2) > 0,:).^2,2); ...
  sum(gradientBasis{2}(sum(gradientBasis{2}.^2,2) > 0,:).^2,2); ...
  sum(gradientBasis{3}(sum(gradientBasis{3}.^2,2) > 0,:).^2,2) ...
  ]);

gradientBasis{1}(indices1,:) = repmat([1 1 1]/min_norm_grad, length(indices1), 1);
gradientBasis{2}(indices2,:) = repmat([1 1 1]/min_norm_grad, length(indices2), 1);
gradientBasis{3}(indices3,:) = repmat([1 1 1]/min_norm_grad, length(indices3), 1);

% Area of each face
triangleAreas = sqrt(hh1.*vv)/2;
triangleAreas(isnan(triangleAreas)) = 0;

% Calculate normals to surface at each face
if dimension == 3
    FaceNormals = cross(w,u);
    FaceNormals = FaceNormals./repmat(sqrt(sum(FaceNormals.^2,2)),1,3);
else
    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

%% ===== TESSELATION TANGENT BUNDLE =====
function tangentPlaneBasis = basis_vertices(VertNormals)
% BASIS_VERTICES  Gives an orthonormal basis orthogonal to several vectors
%
% INPUTS:
%   VertNormals       - list of 3D vectors
%   type              - 'uniform' for normal structure of orthonormal
%                       'polar' for R/theta structure of the orthonormal
% OUTPUTS:
%   tangentPlaneBasis - for each vertex, a pair of vectors defining
%                       the orthonormal basis to that vertex using
%                       the normal to the surface at that vertex
% The cheat: if [x y z] is the normal-to-surface, then the tangent plane
% includes the vector [z-y x-z y-x]

nVertices = size(VertNormals,1); VertNormals = -VertNormals;
tangentPlaneBasis = zeros(nVertices, 3, 2); % Initialize

% First vector in basis: [3-2 1-3 2-1]
tangentPlaneBasis(:,:,1) = diff(VertNormals(:, [2 3 1 2]).').';

% Correct for [1 1 1]-ish vertices: use [y -x 0]
bad = abs(dot(VertNormals, ones(nVertices,3)/sqrt(3), 2)) > 0.97;
tangentPlaneBasis(bad,:,1) = ...
    [VertNormals(bad,2) -VertNormals(bad,1) zeros(sum(bad),1)];

% Second vector in basis by cross product
tangentPlaneBasis(:,:,2) = cross(VertNormals, tangentPlaneBasis(:,:,1));

% Normalize to get orthonormal basis
tangentPlaneBasisNorm = sqrt(sum(tangentPlaneBasis.^2,2));
tangentPlaneBasis = tangentPlaneBasis ./ tangentPlaneBasisNorm(:,[1 1 1],:);

end

%% ===== HORN-SCHUNCK REGULARIZATION MATRIX (FOR MANIFOLD) =====
function [regularizer, gradientBasis, tangentPlaneBasisCell, ...
  tangentPlaneBasis, triangleAreas, FaceNormals, ...
  sparseIndex1,sparseIndex2] = ...
  regularizing_matrix(Faces, Vertices, VertNormals, dimension)
% REGULARIZING_MATRIX: Computation of the regularizing part in the
%                      variationnal approach (SS grad(v_k)grad(v_k')) and
%                      other geometrical quantities.
%
% USAGE: [regularizer, gradientBasis, tangentPlaneBasisCell, ...
%        tangentPlaneBasis, triangleAreas, FaceNormals, ...
%        VertNormals, sparseIndex1, sparseIndex2] = ...
%        regularizing_matrix(tri, coord, dim)
%
% INPUTS
%   Faces                   - triangles of the tesselation
%   Vertices                - vertices of the tesselation
%   VertNormals             - normals to surface at each vertex
%   dimension               - 3 for scalp or cortical surface
%                             2 for flat surface
% OUTPUTS :
%   regularizer             - regularizing matrix
%   gradientBasis           - gradient of the basis functions in FEM
%   tangentPlaneBasisCell   - basis of the tangent plane at a surface node
%                             --> nodes listed according to tesselation
%   tangentPlaneBasis       - basis of the tangent plane at a node
%   triangleAreas           - area of the triangles
%   FaceNormals             - normal of each triangle
%   sparseIndex1            - 1st index of sparse tangent basis magnitudes
%   sparseIndex2            - 2nd index of sparse tangent basis magnitudes
%
%/---Script Authors---------------------\
%|                                      |
%|   *** J.Lef�vre, PhD                 |
%|   julien.lefevre@chups.jussieu.fr    |
%|                                      |
%\--------------------------------------/

nVertices = size(Vertices,1); % Number of nodes
[gradientBasis, triangleAreas, FaceNormals] = ...
  geometry_tesselation(Faces, Vertices, dimension);

% Basis of the tangent plane at each vertex
tangentPlaneBasis = basis_vertices(VertNormals);
tangentPlaneBasisCell = cell(2,3); % similar structure to gradientBasis
% 2 = # of basis vectors, 3 = # of nodes in each element (triangle)

tangentPlaneBasisCell{1,1} = tangentPlaneBasis(Faces(:,1),:,1);
tangentPlaneBasisCell{1,2} = tangentPlaneBasis(Faces(:,2),:,1);
tangentPlaneBasisCell{1,3} = tangentPlaneBasis(Faces(:,3),:,1);

tangentPlaneBasisCell{2,1} = tangentPlaneBasis(Faces(:,1),:,2);
tangentPlaneBasisCell{2,2} = tangentPlaneBasis(Faces(:,2),:,2);
tangentPlaneBasisCell{2,3} = tangentPlaneBasis(Faces(:,3),:,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Regularizing matrix SS grad(v_k)grad(v_k') %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sparseIndex1 = [Faces(:,1) Faces(:,1) Faces(:,2)];
sparseIndex2 = [Faces(:,2) Faces(:,3) Faces(:,3)];
gradientBasisSum = [sum(gradientBasis{1}.*gradientBasis{2},2) ...
  sum(gradientBasis{1}.*gradientBasis{3},2) ...
  sum(gradientBasis{2}.*gradientBasis{3},2)];

tang_scal_11 = [ ...
  sum(tangentPlaneBasisCell{1,1}.*tangentPlaneBasisCell{1,2},2) ...
  sum(tangentPlaneBasisCell{1,1}.*tangentPlaneBasisCell{1,3},2) ...
  sum(tangentPlaneBasisCell{1,2}.*tangentPlaneBasisCell{1,3},2) ...
  ] .* gradientBasisSum .* repmat(triangleAreas, [1 3]);
tang_scal_12 = [ ...
  sum(tangentPlaneBasisCell{1,1}.*tangentPlaneBasisCell{2,2},2) ...
  sum(tangentPlaneBasisCell{1,1}.*tangentPlaneBasisCell{2,3},2) ...
  sum(tangentPlaneBasisCell{1,2}.*tangentPlaneBasisCell{2,3},2) ...
  ] .* gradientBasisSum .* repmat(triangleAreas, [1 3]);
tang_scal_21 = [ ...
  sum(tangentPlaneBasisCell{2,1}.*tangentPlaneBasisCell{1,2},2) ...
  sum(tangentPlaneBasisCell{2,1}.*tangentPlaneBasisCell{1,3},2) ...
  sum(tangentPlaneBasisCell{2,2}.*tangentPlaneBasisCell{1,3},2) ...
  ] .* gradientBasisSum .* repmat(triangleAreas, [1 3]);
tang_scal_22 = [ ...
  sum(tangentPlaneBasisCell{2,1}.*tangentPlaneBasisCell{2,2},2) ...
  sum(tangentPlaneBasisCell{2,1}.*tangentPlaneBasisCell{2,3},2) ...
  sum(tangentPlaneBasisCell{2,2}.*tangentPlaneBasisCell{2,3},2) ...
  ] .* gradientBasisSum .* repmat(triangleAreas, [1 3]);

termes_diag = repmat(triangleAreas, [1 3]) .* [ ...
  sum(gradientBasis{1}.^2,2) ...
  sum(gradientBasis{2}.^2,2) ...
  sum(gradientBasis{3}.^2,2) ];

D = sparse(Faces, Faces, termes_diag, nVertices, nVertices);
E11=sparse(sparseIndex1, sparseIndex2, tang_scal_11, nVertices, nVertices);
E11=E11+E11'+D;
E22=sparse(sparseIndex1,sparseIndex2,tang_scal_22,nVertices,nVertices);
E22=E22+E22'+D;
E12=sparse(sparseIndex1,sparseIndex2,tang_scal_12,nVertices,nVertices);
E21=sparse(sparseIndex1,sparseIndex2,tang_scal_21,nVertices,nVertices);

regularizer = [E11 E12+E21'; ...
               E12'+E21 E22];

end

%% ===== HORN-SCHUNCK DATA FIT MATRIX (FOR MANIFOLD) =====
function [dataFit, B] = data_fit_matrix(Faces, nVertices, ...
  gradientBasis, triangleAreas, tangentPlaneBasisCell, F, dF, ...
  dimension, sparseIndex1, sparseIndex2)
% DATA_FIT_MATRIX   Computation of data-fit matrices of the variational
%                   formulation
% INPUTS:
%   Faces                   - triangles
%   nVertices               - number of nodes
%   gradientBasis           - gradient of the basis functions
%   triangleAreas           - area of each triangle
%   tangentPlaneBasisCell 	- basis of each tangent plane
%   F                       - activity at time step t
%   dF                      - change in activity at time step t
%   dimension               - 3 for scalp or cortical surface
%                             2 for channel surface or cap
%   sparseIndex1            - 1st index of sparse tangent basis magnitudes
%   sparseIndex2            - 2nd index of sparse tangent basis magnitudes
% OUTPUTS:
%   dataFit                 - fit to data matrix:
%                             SS (grad(F).w_k)(grad(F).w_k')
%   B                       - fit to data vector:
%                             -2SS(dF/dt)(grad(F).v_k)

% Gradient of intensity obtained through interpolation
grad_F = repmat(F(Faces(:,1)), 1, dimension) .* gradientBasis{1} ...
  + repmat(F(Faces(:,2)), 1, dimension) .* gradientBasis{2} ...
  + repmat(F(Faces(:,3)), 1, dimension) .* gradientBasis{3};

% Projection of the gradient of F on the tangent space
P_grad_F=cell(1,3); % same structure as gradientBasis : size = nFaces,3 ;
for s=1:3
    for k=1:2
       P_grad_F{s}(:,k)=sum(grad_F.*tangentPlaneBasisCell{k,s},2);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Construction of B %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% m�thode Guillaume Obosinski

B = zeros(2*nVertices, 1);
for k = 1:2
   for s = 1:3
     B(Faces(:,s)+(k-1)*nVertices) = ...
       B(Faces(:,s)+(k-1)*nVertices) + ...
       (-1/12 * triangleAreas .* (P_grad_F{s}(:,k)) .* (dF(Faces(:,s))+sum(dF(Faces),2)));
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   Construction of dataFit   %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% scal_F_11=[];
% scal_F_22=[];
% scal_F_12=[];
% scal_F_21=[];
% scal_F_diag_11=[];
% scal_F_diag_22=[];
% scal_F_diag_12=[];

scal_F_11 = [ ...
  P_grad_F{1}(:,1).*P_grad_F{2}(:,1) ...
  P_grad_F{1}(:,1).*P_grad_F{3}(:,1) ...
  P_grad_F{2}(:,1).*P_grad_F{3}(:,1) ...
  ] .* repmat(triangleAreas, [1 3])/12;
scal_F_12 = [ ...
  P_grad_F{1}(:,1).*P_grad_F{2}(:,2) ...
  P_grad_F{1}(:,1).*P_grad_F{3}(:,2) ...
  P_grad_F{2}(:,1).*P_grad_F{3}(:,2) ...
  ] .* repmat(triangleAreas, [1 3])/12;
scal_F_21 = [ ...
  P_grad_F{1}(:,2).*P_grad_F{2}(:,1) ...
  P_grad_F{1}(:,2).*P_grad_F{3}(:,1) ...
  P_grad_F{2}(:,2).*P_grad_F{3}(:,1) ...
  ] .* repmat(triangleAreas, [1 3])/12;
scal_F_22 = [ ...
  P_grad_F{1}(:,2).*P_grad_F{2}(:,2) ...
  P_grad_F{1}(:,2).*P_grad_F{3}(:,2) ...
  P_grad_F{2}(:,2).*P_grad_F{3}(:,2) ...
  ] .* repmat(triangleAreas, [1 3])/12;

scal_F_diag_11 = [ ...
  P_grad_F{1}(:,1).*P_grad_F{1}(:,1) ...
  P_grad_F{2}(:,1).*P_grad_F{2}(:,1) ...
  P_grad_F{3}(:,1).*P_grad_F{3}(:,1) ...
  ] .* repmat(triangleAreas, [1 3])/6;
scal_F_diag_22 = [ ...
  P_grad_F{1}(:,2).*P_grad_F{1}(:,2) ...
  P_grad_F{2}(:,2).*P_grad_F{2}(:,2) ...
  P_grad_F{3}(:,2).*P_grad_F{3}(:,2) ...
  ] .* repmat(triangleAreas, [1 3])/6;
scal_F_diag_12 = [ ...
  P_grad_F{1}(:,1).*P_grad_F{1}(:,2) ...
  P_grad_F{2}(:,1).*P_grad_F{2}(:,2) ...
  P_grad_F{3}(:,1).*P_grad_F{3}(:,2) ...
  ] .* repmat(triangleAreas, [1 3])/6;

S11=sparse(sparseIndex1, sparseIndex2, scal_F_11, nVertices, nVertices);
S22=sparse(sparseIndex1, sparseIndex2, scal_F_22, nVertices, nVertices);
S12=sparse(sparseIndex1, sparseIndex2, scal_F_12, nVertices, nVertices);
S21=sparse(sparseIndex1, sparseIndex2, scal_F_21, nVertices, nVertices);

D11=sparse(Faces, Faces, scal_F_diag_11, nVertices, nVertices);
D22=sparse(Faces, Faces, scal_F_diag_22, nVertices, nVertices);
D12=sparse(Faces, Faces, scal_F_diag_12, nVertices, nVertices);

dataFit = [S11+S11'+D11 S12+S21'+D12; ...
  S12'+S21+D12 S22+S22'+D22];

end

%% ===== POINCARE INDEX OF FLOW FIELD =====
function index = poincare_index(flowField)
% Compute the index of a vector field VF along a curve (whose it is not
% necessary to give the coordinates !!)
% VF has dimension 2*nbr vectors
  theta = myangle(flowField);
  difftheta = diffangle([theta(2),theta(3),theta(1)], theta);
  index = sum(difftheta)/(2*pi);
end

% gives the good angle of a vector (between 0 and 2pi)
function theta = myangle(flowField)
  normv=sqrt(sum(flowField.^2,1));
  c=flowField(1,:)./normv;
  s=flowField(2,:)./normv;
  theta = acos(c);

  for ii=1:size(flowField,2)
    if s(ii) < 0
      theta(ii)= -theta(ii) + 2*pi;
    end
  end
end

% Difference of two angles (result between -pi and pi)
function theta = diffangle(theta2, theta1)
  theta = theta2 - theta1;
  for ii=1:length(theta)
    if theta(ii) < -pi
      theta(ii) = theta(ii) + 2*pi;
    elseif theta(ii) > pi
        theta(ii) = theta(ii) - 2*pi;
    end
  end
end

author: brainstorm-tools


text:
function CSPFilter = spatial_filtering( x_flt, oldSMC, numPatterns, verbose )

if verbose == 1
%     fprintf( '\tSpatial Filtering...\n\t\t' );
end

CSPFilter = cell( 1, oldSMC.numBands );
for i=1:oldSMC.numBands
    if verbose == 1
        if mod(i, 5) == 0
%             fprintf( '%d', i );
        else
%             fprintf( '.' );
        end
        if mod(i, 100) == 0
%             fprintf( '\n' );
        end
    end

    D1 = x_flt{1, i};
    D2 = x_flt{2, i};
    [W, D] = myTrainCSP( D1, D2 );
    CSPFilter{i}.W = W( :, [1:numPatterns, end-numPatterns+1:end] );
    Dd = diag(D);
    CSPFilter{i}.D = Dd([1:numPatterns, end-numPatterns+1:end]);
end
% fprintf( '\n' );

author: PatternRecognition

repo: OpenBMI


text:
function klu_demo
% KLU demo
%
% Example:
%   klu_demo
%
% See also klu, btf

% Copyright 2004-2009 Timothy A. Davis, Univ. of Florida

load west0479
A = west0479 ;

n = size (A,1) ;
b = rand (n,1) ;

clf
subplot (2,2,1) ;
spy (A)
title ('west0479') ;

subplot (2,2,2) ;
[p, q, r] = btf (A) ;
drawbtf (A, p, q, r) ;
title ('BTF form') ;

[x,info,c] = klu (A, '\', b) ;
matlab_condest = condest (A) ;
matlab_cond = cond (full (A)) ;
fprintf ('MATLAB condest: %g KLU condest: %g cond: %g\n', ...
    matlab_condest, c, matlab_cond) ;

fprintf ('\nKLU with scaling, AMD ordering and condition number estimate:\n') ;
[LU,info] = klu (A, struct ('ordering',0, 'scale', 1)) ;
x = klu (LU, '\', b) ;
resid = norm (A*x-b,1) / norm (A,1) ;
rgrowth = full (min (max (abs ((LU.R \ A (LU.p,LU.q)) - LU.F)) ./ ...
    max (abs (LU.U)))) ;
fprintf ('resid: %g KLU condest: %g rgrowth: %g\n', resid, c, rgrowth) ;
disp (info) ;

subplot (2,2,3) ;
spy (LU.L + LU.U + LU.F) ;
title ('KLU+AMD factors') ;

fprintf ('\nKLU with COLAMD ordering\n') ;
[LU,info] = klu (A, struct ('ordering',1)) ;
x = klu (LU, '\', b) ;
resid = norm (A*x-b,1) / norm (A,1) ;
fprintf ('resid: %g\n', resid) ;
disp (info) ;

subplot (2,2,4) ;
spy (LU.L + LU.U + LU.F) ;
title ('KLU+COLAMD factors') ;

fprintf ('\nKLU with natural ordering (lots of fillin)\n') ;
[x,info] = klu (A, '\', b, struct ('ordering',2)) ;
resid = norm (A*x-b,1) / norm (A,1) ;
fprintf ('resid: %g\n', resid) ;
disp (info) ;

try

    fprintf ('\nKLU with CHOLMOD(A''*A) ordering\n') ;
    [x,info] = klu (A, '\', b, struct ('ordering',3)) ;
    resid = norm (A*x-b,1) / norm (A,1) ;
    fprintf ('resid: %g\n', resid) ;
    disp (info) ;

    fprintf ('\nKLU with CHOLMOD(A+A'') ordering\n') ;
    [x,info] = klu (A, '\', b, struct ('ordering',4)) ;
    resid = norm (A*x-b,1) / norm (A,1) ;
    fprintf ('resid: %g\n', resid) ;
    disp (info) ;

catch
    fprintf ('KLU test with CHOLMOD skipped (CHOLMOD not installed)\n') ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s = ymdhms_to_s_common ( y, m, d, h, n, second )

%*****************************************************************************80
%
%% YMDHMS_TO_S_COMMON "prints" a Common YMDHMS date into a string.
%
%  Format:
%
%    CE YYYY/MM/DD HH:MM:SS
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, H, N, SECOND, the YMDHMS date.
%
%    Output, string S, contains a representation of the date.
%
  if ( y < 0 )
    s = sprintf ( 'BCE %d/%02d/%02d %02d:%02d:%02d', -y, m, d, h, n, second );
  else
    s = sprintf ( 'CE %d/%02d/%02d %02d:%02d:%02d', y, m, d, h, n, second );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function NsVals=reduceStdRefrac2Spher(NMeas,height,expConst,multConst,varargin)
%%REDUCESTDREFRAC2SPHER Given a measurement of the atmospheric refractivity
%                       at a particular height above sea level, determine
%                       the equivalent refractivity at sea level using  a
%                       standard exponential atmospheric model. The
%                       algorithm can handle sea surface refractivities
%                       from 200 to 450.
%
%INPUTS: NMeas The measured atmospheric refractivity. Note that the
%              refractivity is (n-1)*1e6, where n is the index of
%              refraction.
%       height The height in meters at which the measurement of NMeas was
%              taken. One would expect this to be an orthometic height
%              (with respect to mean sea level). However, the model is low
%              fidelity, so using a height with respect to the Earth's
%              reference ellipsoid would presumably not introduce a
%              significant amount of error.
%  expConst, multConst The two optional, positive  parameters
%              parameterizing the decay constant in the model. Assume that
%              at a point, the refractivity is N. Increasing the height by
%              1km, the model is that the refractivity changes by
%              deltaN=-multConst*exp(expConst*N). deltaN cannot be
%              negative. If these parameters are omitted or empty matrices
%              are passed, then the values fitting the data in [1] are
%              used: expConst=0.005577 and multConst=7.32. Note that the
%              value ce=log(Ns/(Ns+DeltaN))/1000 is the decay constant in
%              inverse meters.
%     varargin Parameters to pass to the fminbnd function. These are
%              standard comma-separated keys and values, such as
%              'TolX',1e-8.
%
%OUTPUTS: Ns The indices of refraction reduced to sea level under a
%            standard exponential atmospheric model. There can be 1-2
%            solutions.
%
%The standard exponential atmospheric model is derived in [1]. This
%function determines the refractivity at the given altitudes for sea
%surface refractivities from 200 to 450 and subtracts the observed
%refractivity. Any zero crossings indicate potential solutions. The
%function fminbnd is then used to find each of the zeros given bounded
%regions for the zero crossings.
%
%The use of this type of very basic refraction model is discussed in [2].
%
%REFERENCES:
%[1] B. R. Bean and G. D. Thayer, CRPL Exponential Reference Atmosphere.
%     Washington, D.C.: U. S. Department of Commerce, National Bureau of
%     Standards, Oct. 1959. [Online]. Available:
%     http://digicoll.manoa.hawaii.edu/techreports/PDF/NBS4.pdf
%[2] D. F. Crouse, "Basic tracking using 3D monostatic and bistatic
%    measurements in refractive environments," IEEE Aerospace and
%    Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 54-75, Aug.
%    2014.
%
%June 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(expConst))
    expConst=0.005577;
end

if(nargin<4||isempty(multConst))
	multConst=7.32;
end

%Evaluate the refractivity on a grid.
numPoints=100;
Ns=linspace(200,450,numPoints);
theVals=Ns.*(Ns./(Ns-multConst*exp(expConst*Ns))).^(-height/1000)-NMeas;

%Find where the zero crossings are
diffIdx=find(diff(theVals>0));

%Each crossing gives a region to search to find the zero. Because it is
%bounded, we will use fminbnd to find the minimum of the squared value in
%the bounded region.
numNs=length(diffIdx);
NsVals=zeros(numNs,1);

costFun=@(Ns)(Ns.*(Ns./(Ns-multConst*exp(expConst*Ns))).^(-height/1000)-NMeas)^2;
for curNs=1:numNs
    NsCur=fminbnd(costFun,Ns(diffIdx(curNs)),Ns(diffIdx(curNs)+1),varargin{:});
    NsVals(curNs)=NsCur;
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function membership=paretoset(X)

% PARETOSET  To get the Pareto set from a given set of points.
% synopsis:           membership =paretoset (objectiveMatrix)
% where:
%   objectiveMatrix: [number of points X number of objectives] array
%   membership:      [number of points X 1] logical vector to indicate if ith
%                    point belongs to the Pareto set (true) or not (false).
%
% by Yi Cao, Cranfield University, 02 June 2007
% Revised by Yi Cao on 17 October 2007
% Version 3, 21 October 2007, new sorting scheme to improve speed.
% Bugfix, 25 July 2008, devided by zero error is fixed.
%
% Examples: see paretoset_examples
%

m=size(X,1);
Xmin=min(X);
X1=X-Xmin(ones(m,1),:);     %make sure X1>=0;
Xmean=mean(X1);
%sort X1 so that dominated points can be removed quickly
[x,checklist]=sort(max(X1./(Xmean(ones(m,1),:)+1)));
Y=X(checklist,:);
membership=false(m,1);
while numel(checklist)>1
    k=checklist(1);
    [membership(k),checklist,Y]=paretosub(Y,checklist);
end
membership(checklist)=true;

function [ispareto,nondominated,X]=paretosub(X,checklist)

Z=X-X(ones(size(X,1),1),:);
nondominated=any(Z<0,2);                    %retain nondominated points from the check list
ispareto=all(any(Z(nondominated,:)>0,2));   %check if current point belongs to pareto set
X=X(nondominated,:);
nondominated=checklist(nondominated);

author: Sable

repo: mcbench-benchmarks


text:
function pick = non_maximum_supression(boxes, overlap_threshold,type)

    %NMS
	if isempty(boxes)
        pick = [];
        return;
    end

    % Compute the corners of boxes and the area
	x1 = boxes(:,1);
	y1 = boxes(:,2);
	x2 = boxes(:,3);
	y2 = boxes(:,4);
	s = boxes(:,5);
	area = (x2-x1+1) .* (y2-y1+1);

    % Sorting based on confidence scores
    [vals, I] = sort(s);

	pick = zeros(numel(s),1);

	counter = 1;
	while ~isempty(I)
        last = length(I);
        i = I(last);
        pick(counter) = i;
        counter = counter + 1;

        xx1 = max(x1(i), x1(I(1:last-1)));
        yy1 = max(y1(i), y1(I(1:last-1)));
        xx2 = min(x2(i), x2(I(1:last-1)));
        yy2 = min(y2(i), y2(I(1:last-1)));
        w = max(0.0, xx2-xx1+1);
        h = max(0.0, yy2-yy1+1);
        inter = w.*h;

        if strcmp(type,'Min')
            o = inter ./ min(area(i),area(I(1:last-1)));
        else
            o = inter ./ (area(i) + area(I(1:last-1)) - inter);
        end
        I = I(find(o<=overlap_threshold));
    end

	pick = pick(1:(counter-1));
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function lstAct = mlAct_Matrix2IndexList(mlAct, ml)
lstAct = [];
for iCh = 1:size(ml,1)
    if size(ml,2)==4
        k = find(mlAct(:,1)==ml(iCh,1)  &  mlAct(:,2)==ml(iCh,2)  &  mlAct(:,4)==ml(iCh,4));
    else
        k = find(mlAct(:,1)==ml(iCh,1)  &  mlAct(:,2)==ml(iCh,2));
    end
    if mlAct(k,3)==1
        lstAct = [lstAct; iCh];
    end
end

author: BUNPC


text:
function ell = inertiaEllipse(points)
%INERTIAELLIPSE Inertia ellipse of a set of points
%
%   ELL = inertiaEllipse(PTS);
%   where PTS is a N*2 array containing coordinates of N points, computes
%   the inertia ellispe of the set of points.
%
%   The result has the form:
%   ELL = [XC YC A B THETA],
%   with XC and YC being the center of mass of the point set, A and B are
%   the lengths of the inertia ellipse (see below), and THETA is the angle
%   of the main inertia axis with the horizontal (counted in degrees
%   between 0 and 180).
%   A and B are the standard deviations of the point coordinates when
%   ellipse is aligned with the inertia axes.
%
%   Example
%   pts = randn(100, 2);
%   pts = transformPoint(pts, createScaling(5, 2));
%   pts = transformPoint(pts, createRotation(pi/6));
%   pts = transformPoint(pts, createTranslation(3, 4));
%   ell = inertiaEllipse(pts);
%   figure(1); clf; hold on;
%   drawPoint(pts);
%   drawEllipse(ell, 'linewidth', 2, 'color', 'r');
%
%   See also
%   ellipses2d, drawEllipse
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-02-21,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% HISTORY
% 2009-07-29 take into account ellipse orientation
% 2011-03-12 rewrite using inertia moments

% ellipse center
xc = mean(points(:,1));
yc = mean(points(:,2));

% recenter points
x = points(:,1) - xc;
y = points(:,2) - yc;

% number of points
n = size(points, 1);

% inertia parameters
Ixx = sum(x.^2) / n;
Iyy = sum(y.^2) / n;
Ixy = sum(x.*y) / n;

% compute ellipse semi-axis lengths
common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
ra = sqrt(2) * sqrt(Ixx + Iyy + common);
rb = sqrt(2) * sqrt(Ixx + Iyy - common);

% compute ellipse angle in degrees
theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
theta = rad2deg(theta);

% create the resulting inertia ellipse
ell = [xc yc ra rb theta];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function data = fetchDaily(symbol,apikey)
% Copyright 2020 The MathWorks, Inc.

url = ['https://www.alphavantage.co/query?function=TIME_SERIES_DAILY&symbol=',symbol,'&outputsize=full&apikey=',apikey,'&datatype=json'];

s = webread(url);
fnames = fieldnames(s.TimeSeries_Daily_);
C = struct2cell(s.TimeSeries_Daily_);

for i = 1:length(C)
    Close_Date(i,:) = datetime(fnames{i}(2:end),'InputFormat','yyyy_MM_dd');
    Close_Price(i,:) = str2num(C{i}.x4_Close);
end

data = table(Close_Date,Close_Price);
data = table2timetable(data);

author: matlab-deep-learning


text:
% Find all paths from a start to an end node,
% bounded by a constant, using depth-first-search.
%
% Source: Idea from 6.002x, Spring 2014.
% Note: Uses recursion.
%
% INPUTs: graph structure (a nxn adjacency matrix)
%         s - start node
%         e - end node
%         upperBound = 5, some constant bounding the
%                    length of the path; typically size(adj,1)-1
%         allPaths = {}, path = []; are by default empty,
%                                  serving the recursion
%
% OUTPUTs: a list {} of all shortest paths from "s" to "e",
%                                shorter than "upperBound"
%
% Other routines used: kneighbors.m, DFS.m
% GB: last updated Oct 27 2014

function allPaths = DFS(graph, s, e, allPaths = {}, path = [], upperBound = 5)

        path = [path s];

        if s == e
            if length(path)-1<=upperBound   % this can be modified to include edge weights [(1,2,{}),..]
                allPaths{length(allPaths)+1} = path;
            end
        else
            if length(path)-1<=upperBound

                kneigh = kneighbors(graph, s, 1);

                for node=1:length(kneigh)

                    node = kneigh(node);

                    if length(find(path==node))==0   % avoid cycles
                        allPaths = DFS(graph, node, e, allPaths, path, upperBound);
                    end

                end
            end
        end

author: aeolianine

repo: octave-networks-toolbox


text:
function Archive = UpdateArchive(N,combinePopulation)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Huangke Chen

    % Remove the dominated solutions
    Archive = combinePopulation(NDSort(combinePopulation.objs,1)==1);

    % Update the Archive outPopulation, if the number of solutions is larger than the population size
    if length(Archive) > N
        PopObj = Archive.objs;

        Choose = false(1,size(PopObj,1));
        % Select the extreme solutions
        [~,extreme]     = min(pdist2(PopObj,eye(size(PopObj,2)),'cosine'),[],1);
        Choose(extreme) = true;

        %% Lp-norm-distances between each two solutions
        LpNormD = pdist2(PopObj,PopObj,'minkowski',0.5);
        while sum(Choose) < N
            Remain   = find(~Choose);
            [~, rho] = max(min(LpNormD(Remain,Choose),[],2));
            Choose(Remain(rho)) = true;
        end
        Archive = Archive(Choose);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function ssim  =  cal_ssim( im1, im2, b_row, b_col )

[h w ch]  =  size( im1 );
ssim  = 0;
if ch==1
    ssim   = ssim_index( im1(b_row+1:h-b_row, b_col+1:w-b_col), im2( b_row+1:h-b_row, b_col+1:w-b_col) );
else
    for i = 1:ch
        ssim   = ssim + ssim_index( im1(b_row+1:h-b_row, b_col+1:w-b_col, i), im2( b_row+1:h-b_row, b_col+1:w-b_col, i) );
    end
    ssim   =  ssim/3;
end
return;

function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L)

%========================================================================
%SSIM Index, Version 1.0
%Copyright(c) 2003 Zhou Wang
%All Rights Reserved.
%
%The author was with Howard Hughes Medical Institute, and Laboratory
%for Computational Vision at Center for Neural Science and Courant
%Institute of Mathematical Sciences, New York University, USA. He is
%currently with Department of Electrical and Computer Engineering,
%University of Waterloo, Canada.
%
%----------------------------------------------------------------------
%Permission to use, copy, or modify this software and its documentation
%for educational and research purposes only and without fee is hereby
%granted, provided that this copyright notice and the original authors'
%names appear on all copies and supporting documentation. This program
%shall not be used, rewritten, or adapted as the basis of a commercial
%software or hardware product without first obtaining permission of the
%authors. The authors make no representations about the suitability of
%this software for any purpose. It is provided "as is" without express
%or implied warranty.
%----------------------------------------------------------------------
%
%This is an implementation of the algorithm for calculating the
%Structural SIMilarity (SSIM) index between two images. Please refer
%to the following paper:
%
%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
%quality assessment: From error measurement to structural similarity"
%IEEE Transactios on Image Processing, vol. 13, no. 4, Apr. 2004.
%
%Kindly report any suggestions or corrections to zhouwang@ieee.org
%
%----------------------------------------------------------------------
%
%Input : (1) img1: the first image being compared
%        (2) img2: the second image being compared
%        (3) K: constants in the SSIM index formula (see the above
%            reference). defualt value: K = [0.01 0.03]
%        (4) window: local window for statistics (see the above
%            reference). default widnow is Gaussian given by
%            window = fspecial('gaussian', 11, 1.5);
%        (5) L: dynamic range of the images. default: L = 255
%
%Output: (1) mssim: the mean SSIM index value between 2 images.
%            If one of the images being compared is regarded as
%            perfect quality, then mssim can be considered as the
%            quality measure of the other image.
%            If img1 = img2, then mssim = 1.
%        (2) ssim_map: the SSIM index map of the test image. The map
%            has a smaller size than the input images. The actual size:
%            size(img1) - size(window) + 1.
%
%Default Usage:
%   Given 2 test images img1 and img2, whose dynamic range is 0-255
%
%   [mssim ssim_map] = ssim_index(img1, img2);
%
%Advanced Usage:
%   User defined parameters. For example
%
%   K = [0.05 0.05];
%   window = ones(8);
%   L = 100;
%   [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
%
%See the results:
%
%   mssim                        %Gives the mssim value
%   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
%
%========================================================================

if (nargin < 2 | nargin > 5)
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

if (size(img1) ~= size(img2))
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

[M N] = size(img1);

if (nargin == 2)
   if ((M < 11) | (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian', 11, 1.5);	%
   K(1) = 0.01;								      % default settings
   K(2) = 0.03;								      %
   L = 255;                                  %
end

if (nargin == 3)
   if ((M < 11) | (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian', 11, 1.5);
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 4)
   [H W] = size(window);
   if ((H*W) < 4 | (H > M) | (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 5)
   [H W] = size(window);
   if ((H*W) < 4 | (H > M) | (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
window = window/sum(sum(window));
img1 = double(img1);
img2 = double(img2);

mu1   = filter2(window, img1, 'valid');
mu2   = filter2(window, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 & C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;
   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);
end

mssim = mean2(ssim_map);

return

author: csjunxu


text:
% The COBRA Toolbox: testextractConditionSpecificModel.m
%
% Purpose:
%     - test the extractConditionSpecificModel function
%
% Authors:
%     - Loic Marx, December 2018

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which(mfilename));
cd(fileDir);

% define the solver packages to be used to run this test, can't use
% dqq/Minos for the parallel part.
solverPkgs = prepareTest('needsLP', true, 'excludeSolvers', {'dqqMinos', 'quadMinos'});

% load reference data
load('refData_extractConditionSpecificModel.mat')

% load model
model = createToyModelForLifting(true);

% create a parallel pool
try
    minWorkers = 2;
    myCluster = parcluster(parallel.defaultClusterProfile);
    % no parallel pool
    if myCluster.NumWorkers >= minWorkers
        poolobj = gcp('nocreate');  % if no pool, do not create new one.
        if isempty(poolobj)
            parpool(minWorkers);  % launch minWorkers workers
        end
    end
catch
    % No Parallel pool.
end

% Define input
threshold = 1e-6;

% generate data
modelPruned = extractConditionSpecificModel(model, threshold);

% comparison between refData and generated data
assert(isequal(modelPruned_ref, modelPruned))

% change back to the current directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function output = callcbc(interfacedata)

% Standard input interface
model = yalmip2cbc(interfacedata);

options = interfacedata.options;
model.opts = options.cbc;
model.opts.display = options.verbose;

if options.savedebug
    save cbcdebug model
end

showprogress('Calling CBC',options.showprogress);
solvertime = tic;
if nnz(model.H)==0
    [x,fval,exitflag,nodes,xc] = cbc([],model.f, model.A, model.rl, model.ru, model.lb, model.ub, model.xtype,model.sos,model.x0,model.opts);
else
    [x,fval,exitflag,nodes,xc] = cbc(model.H,model.f, model.A, model.rl, model.ru, model.lb, model.ub, model.xtype,model.sos,model.x0,model.opts);
end
solvertime = toc(solvertime);

% No duals
D_struc = [];
if isempty(x)
    x = zeros(length(c),1);
end

% Check, currently not exhaustive...
problem = 0;
switch exitflag
    case {0,2,6}
        problem = 0;
    case {1,8}
        problem = 1;
    case {3,4}
        problem = 3;
    case 5
        problem = 16;
    case 7
        problem = 2;
    otherwise
        problem = -1;
end

% Save all data sent to solver?
if options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fval = fval;
    solveroutput.exitflag = exitflag;
    solveroutput.nodes = nodes;
    solveroutput.xc = xc;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x(:),[],[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% vgg_fit_hplane_to_x  Fitting hyperplane to set of points.
%
% SYNOPSIS
% [A,e] = vgg_fit_hplane_to_x(s), where
%
%   s ... double(N+1,N+1), inv. covariance matrix s = x'*x, the columns of which
%      are (N+1)-vectors of homog. coordinates of points in N-space (s(end,:)==1)
%   A ... double(1,N+1), fitted hyperplane (homog. coordinates).
%     Homogeneous part of A is normalized, norm(A(1,1:N))==1.
%   e ... double(N,1), eignevalues of the fit
%
% Hyperplane A minimizes sum of squared orthogonal distances of points to it.
%
% EXAMPLE
% Let e be (2,?) vector of detected edgels in nonhomog. coords. Then fitting a straight
% line l to this set of edgels is:
%   e(end,:) = 1;
%   l = vgg_fit_hplane_to_x(e*e');

function [A,e] = vgg_fit_hplane_to_x(s)

N = size(s,1) - 1; % space dimension
c = s(1:N,end)/s(end,end); % centroid

[U,e,V] = svd(s(1:N,1:N)-c*s(end,1:N),0);

A = U(:,end)';
A = [A -A*c];
e = diag(e);

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
%kmgident 'Generate identity or unit matrix '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros mgident.pane file
%
% Parameters:
% Integer: r 'Rows', default: 0: 'rows'
% Integer: c 'Cols', default: 0: 'cols'
% OutputFile: o 'Output Matrix', required: 'Output matrix'
%
% Example: o = kmgident( {'r',0;'c',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% mgident - Generate Identity or Unit Matrix
%
%  DESCRIPTION
% .I mgident
% is used to generate a matrix with all diagonal elements equal to one (unity).
% The output is of type DOUBLE.
%
% "mgident\fR is a pane object pointing to \fImgdiag\fR in the MATRIX
% toolbox.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% mgdiag(1)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kmgident(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kmgident(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'r', 0;'c', 0;'o', '__output'};
maxval={1,1,0};
minval={1,1,0};
istoggle=[1,1,0];
was_set=istoggle * 0;
paramtype={'Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'mgdiag"  -cr'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [H,Ip] = imagesAlign( I, Iref, varargin )
% Fast and robust estimation of homography relating two images.
%
% The algorithm for image alignment is a simple but effective variant of
% the inverse compositional algorithm. For a thorough overview, see:
%   "Lucas-kanade 20 years on A unifying framework,"
%   S. Baker and I. Matthews. IJCV 2004.
% The implementation is optimized and can easily run at 20-30 fps.
%
% type may take on the following values:
%  'translation'  - translation only
%  'rigid'        - translation and rotation
%  'similarity'   - translation, rotation and scale
%  'affine'       - 6 parameter affine transform
%  'rotation'     - pure rotation (about x, y and z)
%  'projective'   - full 8 parameter homography
% Alternatively, type may be a vector of ids between 1 and 8, specifying
% exactly the types of transforms allowed. The ids correspond, to: 1:
% translate-x, 2: translate-y, 3: uniform scale, 4: shear, 5: non-uniform
% scale, 6: rotate-z, 7: rotate-x, 8: rotate-y. For example, to specify
% translation use type=[1,2]. If the transforms don't form a group, the
% returned homography may have more degrees of freedom than expected.
%
% Parameters (in rough order of importance): [resample] controls image
% downsampling prior to computing H. Runtime is proportional to area, so
% using resample<1 can dramatically speed up alignment, and in general not
% degrade performance much. [sig] controls image smoothing, sig=2 gives
% good performance, setting sig too low causes loss of information and too
% high will violate the linearity assumption. [epsilon] defines the
% stopping criteria, use to adjust performance versus speed tradeoff.
% [lambda] is a regularization term that causes small transforms to be
% favored, in general any small non-zero setting of lambda works well.
% [outThr] is a threshold beyond which pixels are considered outliers, be
% careful not to set too low. [minArea] determines coarsest scale beyond
% which the image is not downsampled (should not be set too low). [H0] can
% be used to specify an initial alignment. Use [show] to display results.
%
% USAGE
%  [H,Ip] = imagesAlign( I, Iref, varargin )
%
% INPUTS
%  I          - transformed version of I
%  Iref       - reference grayscale double image
%  varargin   - additional params (struct or name/value pairs)
%   .type       - ['projective'] see above for options
%   .resample   - [1] image resampling prior to homography estimation
%   .sig        - [2] amount of Gaussian spatial smoothing to apply
%   .epsilon    - [1e-3] stopping criteria (min change in error)
%   .lambda     - [1e-6] regularization term favoring small transforms
%   .outThr     - [inf] outlier threshold
%   .minArea    - [4096] minimum image area in coarse to fine search
%   .H0         - [eye(3)] optional initial homography estimate
%   .show       - [0] optionally display results in figure show
%
% OUTPUTS
%  H        - estimated homography to transform I into Iref
%  Ip       - tranformed version of I (slow to compute)
%
% EXAMPLE
%  Iref = double(imread('cameraman.tif'))/255;
%  H0 = [eye(2)+randn(2)*.1 randn(2,1)*10; randn(1,2)*1e-3 1];
%  I = imtransform2(Iref,H0^-1,'pad','replicate');
%  o=50; P=ones(o)*1; I(150:149+o,150:149+o)=P;
%  prmAlign={'outThr',.1,'resample',.5,'type',1:8,'show'};
%  [H,Ip]=imagesAlign(I,Iref,prmAlign{:},1);
%  tic, for i=1:30, H=imagesAlign(I,Iref,prmAlign{:},0); end;
%  t=toc; fprintf('average fps: %f\n',30/t)
%
% See also imTransform2
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.61
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get parameters
dfs={'type','projective','resample',1,'sig',2,'epsilon',1e-3,...
  'lambda',1e-6,'outThr',inf,'minArea',4096,'H0',eye(3),'show',0};
[type,resample,sig,epsilon,lambda,outThr,minArea,H0,show] = ...
  getPrmDflt(varargin,dfs,1);
filt = filterGauss(2*ceil(sig*2.5)+1,[],sig^2);

% determine type of transformation to recover
if(isnumeric(type)), assert(length(type)<=8); else
  id=find(strcmpi(type,{'translation','rigid','similarity','affine',...
    'rotation','projective'})); msgId='piotr:imagesAlign';
  if(isempty(id)), error(msgId,'unknown type: %s',type); end
  type={1:2,[1:2 6],[1:3 6],1:6,6:8,1:8}; type=type{id};
end; keep=zeros(1,8); keep(type)=1; keep=keep>0;

% compute image alignment (optionally resample first)
prm={keep,filt,epsilon,H0,minArea,outThr,lambda};
if( resample==1 ), H=imagesAlign1(I,Iref,prm); else
  S=eye(3); S([1 5])=resample; H0=S*H0*S^-1; prm{4}=H0;
  I1=imResample(I,resample); Iref1=imResample(Iref,resample);
  H=imagesAlign1(I1,Iref1,prm); H=S^-1*H*S;
end

% optionally rectify I and display results (can be expensive)
if(nargout==1 && show==0), return; end
Ip = imtransform2(I,H,'pad','replicate');
if(show), figure(show); clf; s=@(i) subplot(2,3,i);
  Is=[I Iref Ip]; ri=[min(Is(:)) max(Is(:))];
  D0=abs(I-Iref); D1=abs(Ip-Iref); Ds=[D0 D1]; di=[min(Ds(:)) max(Ds(:))];
  s(1); im(I,ri,0); s(2); im(Iref,ri,0); s(3); im(D0,di,0);
  s(4); im(Ip,ri,0); s(5); im(Iref,ri,0); s(6); im(D1,di,0);
  s(3); title('|I-Iref|'); s(6); title('|Ip-Iref|');
end

end

function H = imagesAlign1( I, Iref, prm )

% apply recursively if image large
[keep,filt,epsilon,H0,minArea,outThr,lambda]=deal(prm{:});
[h,w]=size(I); hc=mod(h,2); wc=mod(w,2);
if( w*h<minArea ), H=H0; else
  I1=imResample(I(1:(h-hc),1:(w-wc)),.5);
  Iref1=imResample(Iref(1:(h-hc),1:(w-wc)),.5);
  S=eye(3); S([1 5])=2; H0=S^-1*H0*S; prm{4}=H0;
  H=imagesAlign1(I1,Iref1,prm); H=S*H*S^-1;
end

% smooth images (pad first so dimensions unchanged)
O=ones(1,(length(filt)-1)/2); hs=[O 1:h h*O]; ws=[O 1:w w*O];
Iref=conv2(conv2(Iref(hs,ws),filt','valid'),filt,'valid');
I=conv2(conv2(I(hs,ws),filt','valid'),filt,'valid');

% pad images with nan so later can determine valid regions
hs=[1 1 1:h h h]; ws=[1 1 1:w w w]; I=I(hs,ws); Iref=Iref(hs,ws);
hs=[1:2 h+3:h+4]; I(hs,:)=nan; Iref(hs,:)=nan;
ws=[1:2 w+3:w+4]; I(:,ws)=nan; Iref(:,ws)=nan;

% convert weights hardcoded for 128x128 image to given image dims
wts=[1 1 1.0204 .03125 1.0313 0.0204 .00055516 .00055516];
s=sqrt(numel(Iref))/128;
wts=[wts(1:2) wts(3)^(1/s) wts(4)/s wts(5)^(1/s) wts(6)/s wts(7:8)/(s*s)];

% prepare subspace around Iref
[~,Hs]=ds2H(-ones(1,8),wts); Hs=Hs(:,:,keep); K=size(Hs,3);
[h,w]=size(Iref); Ts=zeros(h,w,K); k=0;
if(keep(1)), k=k+1; Ts(:,1:end-1,k)=Iref(:,2:end); end
if(keep(2)), k=k+1; Ts(1:end-1,:,k)=Iref(2:end,:); end
pTransf={'method','bilinear','pad','none','useCache'};
for i=k+1:K, Ts(:,:,i)=imtransform2(Iref,Hs(:,:,i),pTransf{:},1); end
Ds=Ts-Iref(:,:,ones(1,K)); Mref = ~any(isnan(Ds),3);
if(0), figure(10); montage2(Ds); end
Ds = reshape(Ds,[],size(Ds,3));

% iteratively project Ip onto subspace, storing transformation
lambda=lambda*w*h*eye(K); ds=zeros(1,8); err=inf;
for i=1:100
  s=svd(H); if(s(3)<=1e-4*s(1)), H=eye(3); return; end
  Ip=imtransform2(I,H,pTransf{:},0); dI=Ip-Iref; dI0=abs(dI);
  M=Mref & ~isnan(Ip); M0=M; if(outThr<inf), M=M & dI0<outThr; end
  M1=find(M); D=Ds(M1,:); ds1=(D'*D + lambda)^(-1)*(D'*dI(M1));
  if(any(isnan(ds1))), ds1=zeros(K,1); end
  ds(keep)=ds1; H1=ds2H(ds,wts); H=H*H1; H=H/H(9);
  err0=err; err=dI0; err(~M0)=0; err=mean2(err); del=err0-err;
  if(0), fprintf('i=%03i err=%e del=%e\n',i,err,del); end
  if( del<epsilon ), break; end
end

end

function [H,Hs] = ds2H( ds, wts )
% compute homography from offsets ds
Hs=eye(3); Hs=Hs(:,:,ones(1,8));
Hs(2,3,1)=wts(1)*ds(1);                       % 1 x translation
Hs(1,3,2)=wts(2)*ds(2);                       % 2 y translation
Hs(1:2,1:2,3)=eye(2)*wts(3)^ds(3);            % 3 scale
Hs(2,1,4)=wts(4)*ds(4);                       % 4 shear
Hs(1,1,5)=wts(5)^ds(5);                       % 5 scale non-uniform
ct=cos(wts(6)*ds(6)); st=sin(wts(6)*ds(6));
Hs(1:2,1:2,6)=[ct -st; st ct];                % 6 rotation about z
ct=cos(wts(7)*ds(7)); st=sin(wts(7)*ds(7));
Hs([1 3],[1 3],7)=[ct -st; st ct];            % 7 rotation about x
ct=cos(wts(8)*ds(8)); st=sin(wts(8)*ds(8));
Hs(2:3,2:3,8)=[ct -st; st ct];                % 8 rotation about y
H=eye(3); for i=1:8, H=Hs(:,:,i)*H; end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function a = tril(a,k)
%TRIL         Implements  tril(a,k)  for hessians
%
%   c = tril(a,k)
%
% functionality as Matlab function tril for matrices
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin==1
    k = 0;
  end

  a.x = tril(a.x,k);
  index = ( tril( ones(size(a.x)) , k ) == 0 );
  a.dx(:,index) = 0;
  a.hx(:,index) = 0;

repo: OpenMAS


text:
function result = fit_ML_normal( x,hAx )
% fit_ML_normal - Maximum Likelihood fit of the normal distribution of i.i.d. samples!.
%                  Given the samples of a normal distribution, the PDF parameter is found
%
%    fits data to the probability of the form:
%        p(r) = sqrt(1/2/pi/sig^2)*exp(-((r-u)^2)/(2*sig^2))
%    with parameters: u,sig^2
%
% format:   result = fit_ML_normal( x,hAx )
%
% input:    x   - vector, samples with normal distribution to be parameterized
%           hAx - handle of an axis, on which the fitted distribution is plotted
%                 if h is given empty, a figure is created.
%
% output:   result  - structure with the fields
%                      sig^2,u          - fitted parameters
%                      CRB_sig2,CRB_u   - Cram?r-Rao Bound for the estimator value
%                      RMS              - RMS error of the estimation
%                      type             - 'ML'
%
% example:  fit_ML_normal( randn(1,10000)*3 - 1 )
%               or
%           figure;
%           fit_ML_normal( randn(1,10000)*3 - 1,gca );
%

%
% Algorithm
% ===========
%
% We use the ML algorithm to estimate the PDF from the samples.
% The normal destribution is given by:
%
%    p(x;u,sig^2) = sqrt(1/2/pi/sig^2)*exp(-((x-u).^2)/(2*sig^2))
%
%    where x are the samples which distribute by the function p(x;u,sig^2)
%            and are assumed to be i.i.d !!!
%
% The ML estimator is given by:
%
%    a         = parameters vector = [u,sig^2]
%    f(Xn,a)   = sqrt(1/2/pi/sig^2)*exp(-((Xn-u).^2)/(2*sig^2))
%    L(a)      = f(X,a) = product_by_n( f(Xn,a) )
%              = (2*pi*sig^2)^(-N/2)*exp(-sum((Xn-u)^2)/(2*sig^2))
%    log(L(a)) = -N/2*log(2*pi) - N/2*log(sig^2) - sum((Xn-u)^2)/(2*sig^2)
%
%    The maximum likelihood point is found by the derivative of log(L(a)) with respect to "a":
%
%    diff(log(L(a)),u) = -2*sum(Xn-u)/(2*sig^2) = N/sig^2 * ( u - sum(Xn)/N )
%                    = J(u) * (u_estimation - u)
%    diff(log(L(a)),sig^2) = -N/(2*sig^2) + sum((Xn-u)^2)/(2*sig^4)
%                          = N/(2*sig^4) * ( sum((Xn-u)^2)/N - sig^2 )
%                          = J(sig^2) * (sig^2_estimator - sig^2)
%
%    Therefore, the (efficient) estimators are given by:
%
%               u     = sum(Xn)/N
%               sig^2 = sum((Xn-u)^2)/N
%
%    The Cram?r-Rao Bounds for these estimator are:
%
%               VAR( u )     = 1/J(u)     = (sig^2) / N
%               VAR( sig^2 ) = 1/J(sig^2) = (2*sig^4) / N
%
%    NOTE: the ML estimator does not detect a deviation from the model.
%          therefore, check the RMS value !
%

if (nargin<1)
    error( 'fit_ML_normal - insufficient input arguments' );
end

% Estimation
% =============
x       = x(:);                 % should be column vectors !
N       = length(x);
u       = sum(x)/N;
sig2    = (x-u)'*(x-u)/N;
CRB_u   = sig2 / N;
CRB_sig2= (2*sig2^2) / N;
[n,x_c] = hist( x,100 );
n       = n / sum(n*abs(x_c(2)-x_c(1)));
y       = sqrt(1/2/pi/sig2)*exp(-((x_c-u).^2)/(2*sig2));
RMS     = sqrt( (y-n)*((y-n)')/ (x_c(2)-x_c(1))^2 / (length(x_c)-1) );

% finish summarizing results
% ============================
result = struct( 'u',u,'sig2',sig2,'CRB_u',CRB_u,'CRB_sig2',CRB_sig2,'RMS',RMS,'type','ML' );

% plot distribution if asked for
% ===============================
if (nargin>1)
    xspan = linspace(min(x),max(x),100);
    if ishandle( hAx )
        plot_normal( xspan,result,hAx,1 );
    else
        figure;
        plot_normal( xspan,result,gca,1 );
    end
end


text:
function [c,l] = fix_wavedec(x,n)
% Does a haar wavelet decomposition with n scales.
% Avoids using wavelet toolbox

Lo_D = [ 0.7071 0.7071];
Hi_D = [-0.7071 0.7071];

s = size(x); x = x(:)';
c = [];
l = [length(x)];

dwtEXTM = 'sym';
shift = 0;

for k = 1:n
    lf = length(Lo_D);
    lx = length(x);
    lenEXT = lf-1; lenKEPT = lx+lf-1;

    I = getSymIndices(lx,1);
    y  = x(I);

    x = convdown(y,Lo_D,lenKEPT,shift);
    d = convdown(y,Hi_D,lenKEPT,shift);

    c     = [d c];            % store detail
    l     = [length(d) l];    % store length
end

% Last approximation.
c = [x c];
l = [length(x) l];

if s(1)>1, c = c'; l = l'; end

%-----------------------------------------------------%
% Internal Function(s)
%-----------------------------------------------------%
function y = convdown(x,f,lenKEPT,shift)

y = conv2(x(:)',f(:)'); if size(x,1)>1 , y = y'; end

sx = length(y);
begInd = 1;
[first,last,ok] = GetFirstLast(sx,begInd,lenKEPT);
if ok , y = y(first(1):last(1)); end

y = y(2-rem(shift,2):2:end);

%-----------------------------------------------------%
%----------------------------------------------------------------------------%
function I = getSymIndices(lx,lf)

I = [lf:-1:1 , 1:lx , lx:-1:lx-lf+1];
if lx<lf
    K = (I<1);
    I(K) = 1-I(K);
    J = (I>lx);
    while any(J)
        I(J) = 2*lx+1-I(J);
        K = (I<1);
        I(K) = 1-I(K);
        J = (I>lx);
    end
end
%----------------------------------------------------------------------------%
%----------------------------------------------------------------------------%
function [first,last,ok] = GetFirstLast(sx,begInd,varargin)

oneDIM = isequal(begInd,1);
s = varargin{1}(:)';
if ~oneDIM
    K  = find(s>sx);
    s(K) = sx(K);
    m = find((s < 0) | (s ~= fix(s)));
    ok = isempty(m);
else
    ok = (s>=0) & (s<sx) & (s == fix(s));
end
if ok==0 , first = begInd; last = s; return; end

nbarg = length(varargin);
if nbarg<2, o = 'c'; else , o = lower(varargin{2}); end

err = 0;
if ischar(o(1))
    switch o(1)
        case 'c'
            d = (sx-s)/2;
            if nbarg<3
                if length(o)>1 , side = o(2:end); else , side = 'l'; end
            else
                side = varargin{3};
            end
            if oneDIM
                [first,last] = GetFirst1D(side,sx,d);
            else
                if length(side)<2 , side(2) = 0; end
                for k = 1:2
                    [first(k),last(k)] = GetFirst1D(side(k),sx(k),d(k));
                end
            end

        case {'l','u'} , first = begInd; last = s;
        case {'r','d'} , first = sx-s+1; last = sx;
        otherwise      , err = 1;
    end
else
    first = o; last = first+s-1;
    if ~isequal(first,fix(first)) | any(first<1) | any(last>sx)
        err = 1;
    end
end
if err
    errargt(mfilename,'invalid argument','msg');
    error('*');
end
%----------------------------------------------------------------------------%
function [first,last] = GetFirst1D(side,s,d)

switch side
  case {'u','l','0',0} , first = 1+floor(d); last = s-ceil(d);
  case {'d','r','1',1} , first = 1+ceil(d);  last = s-floor(d);
  otherwise    , first = 1+floor(d); last = s-ceil(d);  % Default is left side
end
%----------------------------------------------------------------------------%

author: csn-le


text:
function [sts, ids, xs, Costs, cprojall] = cpuMPmuFEAT(Params,data,fW,WtW, mu, lam1, nu, ops)

nt0 = ops.nt0;

WtW     = permute(WtW, [1 3 2]);

NT      = Params(1);
nFilt   = Params(2);
Th      = Params(3);

fdata   = fft(data, [], 1);
proj    = real(ifft(fdata .* fW(:,:), [], 1));
if ops.Nrank > 1
    proj    = sum(reshape(proj, NT, nFilt, ops.Nrank),3);
end
trange = int32([-(nt0-1):(nt0-1)]);

xs      = zeros(Params(4), 1, 'single');
ids     = zeros(Params(4), 1, 'int32');
sts     = zeros(Params(4), 1, 'int32');
Costs    = zeros(Params(4), 1, 'single');
cprojall = zeros(Params(4), nFilt, 'single');

i0 = 0;
for k = 1:30
    Ci = bsxfun(@plus, proj, (mu.*lam1)');
    Ci = bsxfun(@rdivide, Ci.^2,  1 + lam1');
    Ci = bsxfun(@minus, Ci, (lam1 .* mu.^2)');

    [mX, id] = max(Ci,[], 2);

    maX         = -my_min(-mX, 31, 1);
    id          = int32(id);

    st                   = find((maX < mX + 1e-3) & mX > Th*Th);
    st(st>NT-nt0 | st<nt0) = [];

    if isempty(st)
       break;
    end
    id      = id(st);

    % inds = bsxfun(@plus, st', [1:nt0]');

    x       = zeros(size(id));
    Cost    = zeros(size(id));
    nsp     = zeros(nFilt,1);
    cproj   = zeros(size(id,1), nFilt, 'single');
    for j = 1:numel(id)
        x(j)            = proj(st(j), id(j));
        Cost(j)         = maX(st(j));
        nsp(id(j))      = nsp(id(j)) + 1;

        % subtract off WtW
        cproj(j,:) =  proj(st(j) ,:);
        proj(st(j) + trange,:) = proj(st(j) + trange,:)  - x(j) * WtW(:,:,id(j));
    end

    xs(i0 + [1:numel(st)])          = x;
    sts(i0 + [1:numel(st)])         = st;
    Costs(i0 + [1:numel(st)])       = Cost;
    ids(i0 + [1:numel(st)])         = id;
    cprojall(i0 + [1:numel(st)], :) = cproj;
    i0 = i0 + numel(st);
end

ids     = ids(1:i0);
xs      = xs(1:i0);
Costs   = Costs(1:i0);
sts     = sts(1:i0);
cprojall = cprojall(1:i0, :);
cprojall = cprojall';

ids = ids - 1;

% keyboard

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function [paramsOut] = referenceGUI(hObject, callbackdata, inputData)%#ok<INUSL>
title = 'Reference parameters';
defaultStruct = inputData.userData.reference;
referenceTypeMenu = {'Robust'; 'Average'; 'Specific'; 'None'};
reportingLevelMenu = {'Minimum'; 'Verbose'};
meanEstimateTypeMenu = {'Huber'; 'Median'; 'Mean'; 'None'};
interpolationOrderMenu = {'Post-reference'; 'Pre-reference'};
while(true)
    mainFigure = findobj('Type', 'Figure', '-and', 'Name', inputData.name);
    userdata = get(mainFigure, 'UserData');
    if isempty(userdata) || ~isfield(userdata, 'reference')
        paramsOut = struct();
    else
        paramsOut = userdata.reference;
    end
    [defaultStruct, errors] = checkStructureDefaults(paramsOut, defaultStruct);

    if ~isempty(errors)
        warning('referenceGUI:bad parameters', getMessageString(errors)); %#ok<CTPCT>
    end

    % Creates structure for
    fNamesDefault = fieldnames(defaultStruct);
    for k = 1:length(fNamesDefault)
        textColorStruct.(fNamesDefault{k}) = 'k';
    end

    if defaultStruct.ransacOff.value
        ransacOffValue = 1;
    else
        ransacOffValue = 0;
    end

    referenceTypeValue = typeMenuPosition(...
        defaultStruct.referenceType.value, referenceTypeMenu);
    reportingLevelValue = typeMenuPosition(...
        defaultStruct.reportingLevel.value, reportingLevelMenu);
    meanEstimateTypeValue = typeMenuPosition( ...
        defaultStruct.meanEstimateType.value, meanEstimateTypeMenu);
    closeOpenWindows(title);
    geometry = {[1, 1.75,1], [1, 1.75,1], [1, 1.75,1], 1, [1.35,1,1.75,1], [1.35,1,1.75,1], ...
        [1.35,1,1.75,1], 1, [1.35,1,1.75,1], [1.35,1,1.75,1], [1.35,1,1.75,1], 1, ...
        [1.35,1,1.75,1], [1.35,1,1.75,1], [1.35,1,1.75,1], 1};
    geomvert = [];
    uilist={{'style', 'text', 'string', 'Reference channels', ...
        'TooltipString', defaultStruct.referenceChannels.description}...
        {'style', 'edit', 'string', ...
        vector2str(defaultStruct.referenceChannels.value), 'tag', ...
        'referenceChannels', 'ForegroundColor', ...
        textColorStruct.referenceChannels}...
        {'style', 'text', 'string', ''} ...
        {'style', 'text', 'string', 'Evaluation channels', ...
        'TooltipString', defaultStruct.evaluationChannels.description}...
        {'style', 'edit', 'string', ...
        vector2str(defaultStruct.evaluationChannels.value), ...
        'tag', 'evaluationChannels', 'ForegroundColor', ...
        textColorStruct.evaluationChannels}...
        {'style', 'text', 'string', ''} ...
        {'style', 'text', 'string', 'Re-referenced channels', ...
        'TooltipString', defaultStruct.rereferencedChannels.description}...
        {'style', 'edit', 'string', ...
        vector2str(defaultStruct.rereferencedChannels.value), ...
        'tag', 'rereferencedChannels', 'ForegroundColor', ...
        textColorStruct.rereferencedChannels}...
        {'style', 'text', 'string', ''} ...
        {'style', 'text', 'string', ''}...
        {'style', 'text', 'string', 'Turn ransac off', ...
        'TooltipString', defaultStruct.ransacOff.description}...
        {'style', 'checkbox', 'Value', ransacOffValue, ...
        'tag', 'ransacOff', 'ForegroundColor', textColorStruct.ransacOff}...
        {'style', 'text', 'string', 'Ransac sample size', ...
        'TooltipString', defaultStruct.ransacSampleSize.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.ransacSampleSize.value), ...
        'tag', 'ransacSampleSize', 'ForegroundColor', ...
        textColorStruct.ransacSampleSize}...
        {'style', 'text', 'string', 'Ransac channel fraction', ...
        'TooltipString', defaultStruct.ransacChannelFraction.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.ransacChannelFraction.value), ...
        'tag','ransacChannelFraction', 'ForegroundColor', ...
        textColorStruct.ransacChannelFraction}...
        {'style', 'text', 'string', 'Ransac correlation threshold', ...
        'TooltipString', defaultStruct.ransacCorrelationThreshold.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.ransacCorrelationThreshold.value), ...
        'tag','ransacCorrelationThreshold', 'ForegroundColor', ...
        textColorStruct.ransacCorrelationThreshold}...
        {'style', 'text', 'string','Ransac unbroken time', ...
        'TooltipString', defaultStruct.ransacUnbrokenTime.description}...
        {'style','edit','string', ...
        num2str(defaultStruct.ransacUnbrokenTime.value), ...
        'tag', 'ransacUnbrokenTime', 'ForegroundColor', ...
        textColorStruct.ransacUnbrokenTime}...
        {'style', 'text', 'string', 'Ransac window', ...
        'TooltipString', defaultStruct.ransacWindowSeconds.description}...
        {'style','edit','string', ...
        num2str(defaultStruct.ransacWindowSeconds.value), ...
        'tag', 'ransacWindowSeconds', 'ForegroundColor', ...
        textColorStruct.ransacWindowSeconds}...
        {'style', 'text', 'string', ''}...
        {'style', 'text', 'string', 'Sampling rate', ...
        'TooltipString', defaultStruct.srate.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.srate.value), ...
        'tag', 'srate', 'ForegroundColor', textColorStruct.srate}...
        {'style', 'text', 'string', 'Robust deviation threshold', ...
        'TooltipString', defaultStruct.robustDeviationThreshold.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.robustDeviationThreshold.value), ...
        'tag', 'robustDeviationThreshold', 'ForegroundColor', ...
        textColorStruct.robustDeviationThreshold}...
        {'style', 'text', 'string', 'Correlation window', ...
        'TooltipString', defaultStruct.correlationWindowSeconds.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.correlationWindowSeconds.value), ...
        'tag','correlationWindowSeconds', 'ForegroundColor', ...
        textColorStruct.correlationWindowSeconds}...
        {'style', 'text', 'string', 'High frequency noise threshold', ...
        'TooltipString', defaultStruct.highFrequencyNoiseThreshold.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.highFrequencyNoiseThreshold.value), ...
        'tag', 'highFrequencyNoiseThreshold', 'ForegroundColor', ...
        textColorStruct.highFrequencyNoiseThreshold}...
        {'style', 'text', 'string','Correlation threshold', ...
        'TooltipString',defaultStruct.correlationThreshold.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.correlationThreshold.value), ...
        'tag', 'correlationThreshold', 'ForegroundColor', ...
        textColorStruct.correlationThreshold}...
        {'style', 'text', 'string', 'Bad time threshold', ...
        'TooltipString', defaultStruct.badTimeThreshold.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.badTimeThreshold.value), ...
        'tag', 'badTimeThreshold', 'ForegroundColor', ...
        textColorStruct.badTimeThreshold}...
        {'style', 'text', 'string', ''}...
        {'style', 'text', 'string', 'Max reference iterations', ...
        'TooltipString', defaultStruct.maxReferenceIterations.description}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.maxReferenceIterations.value), ...
        'tag', 'maxReferenceIterations', 'ForegroundColor', ...
        textColorStruct.maxReferenceIterations}...
        {'style', 'text', 'string', 'Reference type', ...
        'TooltipString', defaultStruct.referenceType.description}...
        {'style', 'popupmenu', 'string', 'Robust|Average|Specific', ...
        'value', referenceTypeValue, 'tag', 'referenceType', ...
        'ForegroundColor', textColorStruct.referenceType}...
        {'style', 'text', 'string', 'Reporting level', ...
        'TooltipString', defaultStruct.reportingLevel.description}...
        {'style', 'popupmenu', 'string', 'Minimum|Verbose', ...
        'value', reportingLevelValue, ...
        'tag', 'reportingLevel', 'ForegroundColor', ...
        textColorStruct.reportingLevel}...
        {'style', 'text', 'string', 'Interpolation order', ...
        'TooltipString', defaultStruct.interpolationOrder.description}...
        {'style', 'popupmenu', 'string', 'Post-reference|Pre-reference', ...
        'tag', 'interpolationOrder', 'ForegroundColor', ...
        textColorStruct.interpolationOrder}...
        {'style', 'text', 'string', 'Mean estimate type', ...
        'TooltipString', defaultStruct.meanEstimateType.description}...
        {'style', 'popupmenu', 'string', 'Huber|Median|Mean|None', ...
        'value', meanEstimateTypeValue, 'tag', 'meanEstimateType', ...
        'ForegroundColor', textColorStruct.reportingLevel}...
        {'style', 'edit', 'string', ...
        num2str(defaultStruct.samples.value), ...
        'tag', 'samples', 'visible', 'off'}...
        {'style', 'text', 'string', ''},...
        {'style', 'text', 'string', ''}};

    [~, ~, ~, paramsOut] = ...
        inputgui('geometry', geometry, 'geomvert', geomvert, ...
        'uilist', uilist, 'title', title, 'helpcom', ...
        'pophelp(''pop_prepPipeline'')');
    if(isempty(paramsOut))
        break;
    end
    [paramsOut, typeErrors, fNamesErrors] = ...
        changeType(paramsOut, defaultStruct);
    paramsOut.referenceType = ...
        typeMenuString(paramsOut.referenceType, referenceTypeMenu);
    paramsOut.reportingLevel = ...
        typeMenuString(paramsOut.reportingLevel, reportingLevelMenu);
    paramsOut.interpolationOrder = ...
        typeMenuString(paramsOut.interpolationOrder, interpolationOrderMenu);
    paramsOut.meanEstimateType = ...
        typeMenuString(paramsOut.meanEstimateType, meanEstimateTypeMenu);
    mainFigure = findobj('Type', 'Figure', '-and', 'Name', inputData.name);
    userdata = get(mainFigure, 'UserData');
    userdata.reference = paramsOut;
    set(mainFigure, 'UserData', userdata);
    if isempty(typeErrors)
        break;
    end
    textColorStruct = highlightErrors(fNamesErrors, ...
        fNamesDefault, textColorStruct);
    displayErrors(typeErrors); % Displays the errors and restarts GUI
end

    function position = typeMenuPosition(theString, menuString)
        menuIndex = 1:length(menuString);
        position = strcmpi(theString, menuString);
        position = menuIndex(position);
    end

    function theString = typeMenuString(position, menuString)
        theString = menuString{position};
    end

end

author: VisLab

repo: EEG-Clean-Tools


text:
function install_sedumi( varargin )

%SeDuMi installation script
%
% Heavy rewrite by Michael C. Grant for SeDuMi 1.34
% Now detects if binaries are already
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

need_rebuild = any( strcmp( varargin, '-rebuild' ) );

targets64={...
    'bwblkslv.c sdmauxFill.c sdmauxRdot.c',...
    'choltmpsiz.c',...
    'cholsplit.c',...
    'dpr1fact.c auxfwdpr1.c sdmauxCone.c  sdmauxCmp.c sdmauxFill.c sdmauxScalarmul.c sdmauxRdot.c blkaux.c',...
    'symfctmex.c symfct.c',...
    'ordmmdmex.c ordmmd.c',...
    'quadadd.c',...
    'sqrtinv.c sdmauxCone.c',...
    'givensrot.c auxgivens.c sdmauxCone.c',...
    'urotorder.c auxgivens.c sdmauxCone.c sdmauxTriu.c sdmauxRdot.c',...
    'psdframeit.c reflect.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c',...
    'psdinvjmul.c reflect.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c',...
    'bwdpr1.c sdmauxCone.c sdmauxRdot.c',...
    'fwdpr1.c auxfwdpr1.c sdmauxCone.c sdmauxScalarmul.c',...
    'fwblkslv.c sdmauxScalarmul.c',...
    'qblkmul.c sdmauxScalarmul.c',...
    'blkchol.c blkchol2.c sdmauxFill.c sdmauxScalarmul.c',...
    'vecsym.c sdmauxCone.c',...
    'qrK.c sdmauxCone.c sdmauxRdot.c sdmauxScalarmul.c',...
    'finsymbden.c sdmauxCmp.c',...
    'symbfwblk.c',...
    'ddot.c sdmauxCone.c sdmauxRdot.c sdmauxScalarmul.c',...
    'partitA.c sdmauxCmp.c',...
    'getada1.c sdmauxFill.c',...
    'getada2.c sdmauxCone.c sdmauxRdot.c sdmauxFill.c',...
    'getada3.c spscale.c sdmauxCone.c sdmauxRdot.c sdmauxScalarmul.c sdmauxCmp.c',...
    'adendotd.c sdmauxCone.c',...
    'adenscale.c',...
    'extractA.c',...
    'sortnnz.c sdmauxCmp.c',...
    'iswnbr.c',...
    'incorder.c',...
    'findblks.c sdmauxCone.c sdmauxCmp.c',...
    'invcholfac.c triuaux.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c',...
    };

fs = filesep;
mpath = mfilename('fullpath');
mpath = mpath( 1 : max([1,strfind(mpath,fs)]) - 1 );
ISOCTAVE = exist('OCTAVE_VERSION','builtin');
VERSION  = [1,0.01]*sscanf(version,'%d.%d');
if ISOCTAVE, prog = 'Octave'; else prog = 'Matlab'; end
COMPUTER = computer;
mext = mexext;

%
% We don't want to rebuild the binaries if they're already present, unless
% the user has specifically asked for a rebuild. Frankly, we can't
% guarantee that rebuilding will work.
%

if ISOCTAVE,
    page_output_immediately( true, 'local' );
end

line = '---------------------------------------------------------------------------';
fprintf( '\n%s\nSeDuMi installation script\n   Directory: %s\n   %s %s on %s\n%s\n', ...
    line, mpath, prog, version, COMPUTER, line );

if ~need_rebuild,
    fprintf( 'Looking for existing binaries...' );
    mdir = '';
    if ISOCTAVE && VERSION > 3.08,
        if ispc,
            mdir = 'o_win32';
        elseif ismac,
            mdir = 'o_mac32';
        elseif isunix && any( strfind( COMPUTER, 'linux' ) ),
            mdir = 'o_lin32';
        end
        if ~isempty(mdir) && strncmpi( COMPUTER, 'x86_64', 6 ),
            mdir(end-1:end) = '64';
        end
        if ~exist( [ mpath, fs, mdir ], 'dir' ),
            mdir = '';
        end
    end
    nfound = [ 0, 0 ];
    for k = 1 : length(targets64),
        targ = targets64{k};
        targ = [ targ(1:min(strfind(targ,'.'))), mext ];
        if exist( [ mpath, fs, targ ], 'file' ),
            nfound(1) = nfound(1) + 1;
        elseif ~isempty(mdir) && exist( [ mpath, fs, mdir, fs, targ ], 'file' ),
            nfound(2) = nfound(2) + 1;
        end
    end
    if sum(nfound) == 0,
        fprintf( 'none found; building...\n' );
        need_rebuild = true;
    elseif sum(nfound) < length(targets64),
        fprintf( 'incomplete set found.\n' );
        disp( line );
        disp( 'Some of the binaries for this platform were found, but some' );
        disp( 'were missing as well. This may mean your download was corrupt;' );
        disp( 'consider downloading and unpacking SeDuMi again. Otherwise, to' );
        disp( 'try rebuilding the MEX files yourself, run this command:' );
        disp( '    install_sedumi -rebuild' );
        fprintf( '%s\n\n', line );
        return;
    else
        fprintf( 'found!\n' );
        fprintf( '   If for some reason you need to rebuild the binaries, use this command:\n' );
        fprintf( '      install_sedumi -rebuild\n' );
    end
else
    nfound = [1,0];
end

if need_rebuild,
    disp( 'Attempting to recompile the SeDuMi binaries:' );
    % Note the use of 0.01 here. That's because version 7 had more than 10
    % minor releases, so 7.10-7.14 need to be ordered after 7.01-7.09.
    IS64BIT  = ~ISOCTAVE & strcmp(COMPUTER(end-1:end),'64');
    if ispc,
        flags = {'-DPC'};
    elseif isunix,
        flags = {'-DUNIX'};
    end
    libs = {};
    if ISOCTAVE,
        % Octave has mwSize and mwIndex hardcoded in mex.h as ints.
        % There is no definition for mwSignedIndex so include it here.
        % This means that Octave ignores the -largeArrayDims flag.
        if VERSION < 3.08,
            flags{end+1} = '-DmwSignedIndex=int';
        end
        libs{end+1} = '-lblas';
    else
        flags{end+1} = '-O';
        if IS64BIT && ( VERSION >= 7.03 ),
            flags{end+1} = '-largeArrayDims';
        elseif VERSION < 7.03,
            flags{end+1} = '-DmwIndex=int';
            flags{end+1} = '-DmwSize=int';
            flags{end+1} = '-DmwSignedIndex=int';
        end
        if VERSION >= 7 && ispc,
            if IS64BIT, dirval = 'win64'; else dirval = 'win32'; end
            libdir = [ matlabroot, fs, 'external', fs, 'lib', fs, dirval, fs ];
            if exist( [ libdir, 'microsoft' ], 'dir' ),
                libdir = [ libdir, 'microsoft' ];
                found = true;
            elseif exist( [ libdir, 'msvc60' ], 'dir' ),
                libdir = [ libdir, 'msvc60' ];
                found = true;
            elseif exist( [ libdir, 'lcc' ], 'dir' ),
                libdir = [ libdir, 'lcc' ];
                found = true;
            end
            if found,
                libs{end+1} = [ '-L"', libdir, '"' ];
            end
        end
        if VERSION >= 7.05,
            libs{end+1} = '-lmwblas';
        else
            libs{end+1} = '-lmwlapack';
        end
    end
    libs = sprintf( ' %s', libs{:} );
    flags = sprintf( ' %s', flags{:} );
    olddir = pwd;
    cd( mpath );
    failed = false;
    fprintf( 'Template: mex%s <sources>%s\n', flags, libs );
    for i=1:length(targets64),
        targ = targets64{i};
        mfile = [ targ(1:min(strfind(targ,'.'))), mext ];
        temp = [ 'mex ', flags, ' ', targets64{i}, libs ];
        fprintf( '   %s: %s\n', mfile, targ );
        eval( temp, 'failed=true;' ); %#ok
    end
    cd( olddir );
    if failed,
        fprintf( 'At least one compilation failure occurred.\n' );
        nfound = [0,0];
    else
        fprintf( 'Compilation successful.\n' );
        nfound = [1,0];
    end
end

if any(nfound),
    disp( line );
    fprintf( 'Adding SeDuMi to the %s path:\n', prog );
    ps = pathsep;
    pp = [ ps, path, ps ];
    already = true;
    fprintf( '   Base directory...' );
    if ~any(strfind(pp,[ps,mpath,ps])),
        already = false;
        pp = [ pp, mpath, ps ];
        fprintf( 'added.\n' );
    else
        fprintf( 'already there.\n' );
    end
    if nfound(2),
        fprintf( '   Binaries directory...' );
        if nfound(2) && ~any(strfind(pp,[ps,mpath,fs,mdir,ps])),
            already = false;
            pp = [ pp, mpath, fs, mdir, ps ];
            fprintf( 'added.\n' );
        else
            fprintf( 'already there.\n' );
        end
    end
    fprintf( '   Conversion directory...' );
    if ~any(strfind(pp,[ps,mpath,fs,'conversion',ps])),
        already = false;
        pp = [ pp, mpath, fs, 'conversion', ps ];
        fprintf( 'added.\n' );
    else
        fprintf( 'already there.\n' );
    end
    fprintf( '   Examples directory...'  );
    if ~any(strfind(pp,[ps,mpath,fs,'examples',ps])),
        already = false;
        pp = [ pp, mpath, fs, 'examples', ps ];
        fprintf( 'added.\n' );
    else
        fprintf( 'already there.\n' );
    end
    if ~already,
        path(pp);
        fprintf( 'Please save the %s path if you want to use SeDuMi from any directory.\n', prog );
    end
    disp( line );
    disp('SeDuMi has been succesfully installed.' );
    disp( 'For more information, type "help sedumi" or see the user guide.')
else
    disp( line );
    disp( 'SeDuMi was not successfully installed.' );
    disp( 'Please attempt to correct the errors and try again.' );
end

fprintf('%s\n\n',line);

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function out = ClampImg(img, a, b)
%
%
%       out = ClampImg(img, a, b)
%
%
%       Image is clamped if its values are over [a,b]
%
%       Input:
%           -img: the input img to be clamped
%           -a: the lower bound
%           -b: the upper bound
%
%       Output:
%           -out: output clamped image
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

out = img;

out(img < a) = a;
out(img > b) = b;

end

author: banterle


text:
%=========================================================================
% Sparse Representation-based Classification for Fusion (SRCF) , v1.0
%
% M. Nejati, S. Samavi, S. Shirani, "Multi-focus Image Fusion Using
% Dictionary-Based Sparse Representation", Information Fusion, vol. 25,
% Sept. 2015, pp. 72-84.
%
% Please refer to the above paper if you use this software.
%=========================================================================
clc,clear all,close all;
addpath('Utils','Utils\ompbox10','Utils\Metrics','Utils\GCO\bin','Utils\GCO\');

%====< Test Dataset >====%
dataPath = 'Data\GrayscaleDataset';
% dataPath = 'Data\LytroDataset';
[dataset,Dict] = loadData(dataPath);

%====< Fusion of Mutifocus Image Pairs >====%
res = struct([]);
for i = 1:dataset.numImage

    fprintf('>> start of fusion for image (%d)...\n',i);

    %----< load source image pairs >----%
    imgA0 = imread(fullfile(dataset.dataPath, dataset.imagesA{i}));
    imgB0 = imread(fullfile(dataset.dataPath, dataset.imagesB{i}));
    if strcmp(dataset.imgSet,'Grayscale')
    if (size(imgA0,3)>1), imgA0 = rgb2gray(imgA0); end
    if (size(imgB0,3)>1), imgB0 = rgb2gray(imgB0); end
    end

    %---< Fusion >---%
    imgF0 = SRCF(imgA0, imgB0, Dict);

    if (size(imgA0,3)>1), imgA = rgb2gray(imgA0); else imgA = imgA0; end
    if (size(imgB0,3)>1), imgB = rgb2gray(imgB0); else imgB = imgB0; end
    if (size(imgF0,3)>1), imgF = rgb2gray(imgF0); else imgF = imgF0; end

    %=======================
    % Objective Evaluation
    %=======================
    MI = fusionMI(imgA,imgB,imgF); % no matter double or uint8
    NMI = fusionECC(imgA,imgB,imgF); % no matter double or uint8
    QABF = Qabf_eval(imgA,imgB,imgF); % no matter double or uint8
    VIFF = VIFF_Public(imgA,imgB,imgF); % no matter double or uint8
    fprintf('-------------------------------------------------------\n')
    fprintf(' %s: QABF = %.4f , VIFF = %.4f , NMI = %.4f\n',dataset.imagesA{i}(1:end-4),QABF,VIFF,NMI);
    fprintf('-------------------------------------------------------\n')
    res(i).MI = MI;
    res(i).NMI = NMI;
    res(i).QABF = QABF;
    res(i).VIFF = VIFF;

    %---< Display >---%
%     figure,imshow(imgA0,[]),title('Source Image 1')
%     figure,imshow(imgB0,[]),title('Source Image 2')
    figure,imshow(imgF0,[]),title('Fused Image')
    pause(2);
    close all;

end

author: thfylsty


text:
function varargout = lz4decode(varargin)
%
% output = lz4decode(input)
%    or
% output = lz4decode(input,info)
%
% Decompressing an LZ4-compressed byte-stream to recover the original data
% This function depends on the ZMat toolbox (http://github.com/fangq/zmat)
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      input: a string, int8/uint8 vector or numerical array to store LZ4-compressed data
%      info (optional): a struct produced by the zmat/lz4encode function during
%            compression; if not given, the inputs/outputs will be treated as a
%            1-D vector
%
% output:
%      output: the decompressed byte stream stored in a uint8 vector; if info is
%            given, output will restore the original data's type and dimensions
%
% examples:
%      [bytes, info]=lz4encode(eye(10));
%      orig=lz4decode(bytes,info);
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==0)
    error('you must provide at least 1 input');
end
if(exist('zmat','file')==2 || exist('zmat','file')==3)
    if(nargin>1)
        [varargout{1:nargout}]=zmat(varargin{1},varargin{2:end});
    else
        [varargout{1:nargout}]=zmat(varargin{1},0,'lz4',varargin{2:end});
    end
else
    error('you must install ZMat toolbox to use this feature: http://github.com/fangq/zmat')
end

author: fieldtrip

repo: fieldtrip