text:
function val = fgGet(fg,param,varargin)
%Get values from a fiber group structure
%
%  val = fgGet(fg,param,varargin)
%
% Parameters
% General
%    'name'
%    'type'
%    'colorrgb'
%    'thickness'
%    'visible'
%
% Fiber related
%    'nfibers'- Number of fibers in this group
%    'nodes per fiber'  - Number of nodes per fiber.
%    'fibers' - Fiber coordinates
%    'fibernames'
%    'fiberindex'
%
% ROI and image coord related
%    'unique image coords'
%    'nodes to imagecoords' -
%    'voxel2fiber node pairs' - For each roi coord, an Nx2 matrix of
%        (fiber number,node number)
%    'nodes in voxels' - Nodes inside the voxels of roi coords
%    'voxels in fg'  - Cell array of the roiCoords touched by each fiber
%    'voxels2fibermatrix' - Binary matrix (voxels by fibers).  1s when a
%        fiber is in a voxel of the roiCoords (which are, sadly, implicit).
%
% Tensor and tractography related
%      'tensors'     - Tensors for each node
%
%
% See also: dwiGet/Set, fgCreate; fgSet
%
% (c) Stanford VISTA Team

% NOTES:
% Programming TODO:
%   We should store the transforms needed to shift the fg coordinates
%   between acpc and image space.

% I eliminated these checks because this function is now called many many
% times and this slows the computations (Franco).
%
%if notDefined('fg'),error('fiber group required.'); end
% if notDefined('param'), error('param required.'); end

val = [];

switch mrvParamFormat(param)

  % Basic fiber parameters
  case 'name'
    val = fg.name;
  case 'type'  % Should always be fibergroup
    val = fg.type;

    % Fiber visualization settings.
  case 'colorrgb'
    val = fg.colorRgb;
  case 'thickness'
    val = fg.thickness;
  case 'visible'
    val = fg.visible;

    % Simple fiber properties --
  case {'fibers'}
    % val = fgGet(fg,'fibers',fList);
    %
    % Returns a 3xN matrix of fiber coordinates corresponding to the
    % fibers specified in the integer vector, fList.  This differs from
    % the dtiH (mrDiffusion) representation, where fiber coordinates
    % are stored as a set of cell arrays for each fiber.
    if ~isempty(varargin)
      list = varargin{1};
      val = cell(length(list),1);
      for ii=1:length(list)
        val{ii} = fg.fibers{ii};
      end
    else
      val = fg.fibers;
    end
  case 'fibernames'
    val = fg.fiberNames;
  case 'fiberindex'
    val = fg.fiberIndex;
  case 'nfibers'
    val = length(fg.fibers);
  case {'nodesperfiber','nsamplesperfiber','nfibersamples'}
    % fgGet(fg,'n samples per fiber ')
    % How many samples per fiber.  This is about equal to
    % their length in mm, though we need to write the fiber lengths
    % routine to actually calculate this.
    nFibers = fgGet(fg,'n fibers');
    val = zeros(1,nFibers);
    for ii=1:nFibers
      val(ii) = length(fg.fibers{ii});
    end

    % Fiber group (subgroup) properties.
    % These are used when we classify fibers into subgroups.  We should
    % probably clean up this organization which is currently
    %
    %   subgroup - length of fibers, an index of group identity
    %   subgroupNames()
    %    .subgroupIndex - Probably should go away and the index should
    %          just be
    %    .subgroupName  - Probably should be moved up.
    %
  case {'ngroups','nsubgroups'}
    val = length(fg.subgroupNames);
  case {'groupnames'}
    val = cell(1,fgGet(fg,'n groups'));
    for ii=1:nGroups
      val{ii} = fg.subgroupNames(ii).subgroupName;
    end

    % DTI properties
  case 'tensors'
    val = fg.tensors;

    % Fiber to coord calculations
  case {'imagecoords'}
    % c = fgGet(fgAcpc,'image coords',fgList,xForm);
    % c = fgGet(fgAcpc,'image coords',fgList,xForm);
    %
    % Return the image coordinates of a specified list of fibers
    % Returns a matrix that is fgList by 3 of the image coordinates for
    % each node of each fiber.
    %
    % Fiber coords are represented at fine resolution in ACPC space.
    % These coordinates are rounded and in image space
    if ~isempty(varargin)
      fList = varargin{1};
      if length(varargin) > 1
        xForm = varargin{2};
        % Put the fiber coordinates into image space
        fg = dtiXformFiberCoords(fg,xForm);
      end
    else
      % In this case, the fiber coords should already be in image
      % space.
      nFibers = fgGet(fg,'n fibers');
      fList = 1:nFibers;
    end

    % Pull out the coordinates and floor them.  These are in image
    % space.
    nFibers = length(fList);
    val = cell(1,nFibers);
    if nFibers == 1
      %val = round(fg.fibers{fList(1)}');
      val = floor(fg.fibers{fList(1)}');

    else
      for ii=1:nFibers
        %val{ii} = round(fg.fibers{fList(ii)}');
        val{ii} = floor(fg.fibers{fList(ii)}');

      end
    end

  case {'uniqueimagecoords'}
    %   coords = fgGet(fgIMG,'unique image coords');
    %
    % The fg input must be in IMG space.
    %
    % Returns the unique image coordinates of all the fibers as an Nx3
    % matrix of integers.
    % val = round(horzcat(fg.fibers{:})');
    val = floor(horzcat(fg.fibers{:})');
    val = unique(val,'rows');

  case {'nodes2voxels'}
    %   nodes2voxels = fgGet(fgImg,'nodes2voxels',roiCoords)
    %
    % The roiCoords are a matrix of Nx3 coordinates.  They describe a
    % region of interest, typically in image space or possibly in acpc
    % space.
    %
    % We return a cell array that is a mapping of fiber nodes to voxels in
    % the roi.  The roi is specified as an Nx3 matrix of coordinates.
    % The returned cell array, nodes2voxels, has the same number of
    % cells as there are fibers.
    %
    % Unlike the fiber group cells, which have a 3D coordinate of each
    % node, this cell array has an integer that indexes the row of
    % roiCoords that contains the node. If a node is not in any of the
    % roiCoords, the entry in node2voxels{ii} for that node is zero.
    % This means that node is outside the 'roiCoords'.
    %
    % Once again: The cell nodes2voxels{ii} specifies whether each
    % node in the iith fiber is inside a voxel in the roiCoords.  The
    % value specifies the row in roiCoords that contains the node.
    %
    if isempty(varargin), error('roiCoords required');
    else
      roiCoords = varargin{1};
    end

    % Find the roiCoord for each node in each fiber.
    nFiber = fgGet(fg,'n fibers');
    val    = cell(nFiber,1);
    for ii=1:nFiber
      % if ~mod(ii,200), fprintf('%d ',ii); end
      % Node coordinates in image space
      nodeCoords = fgGet(fg,'image coords',ii);

      % The values in loc are the row of the coords matrix that contains
      % that sample point in a fiber.  For example, if the number 100 is
      % in the 10th position of loc, then the 10th sample point in the
      % fiber passes through the voxel in row 100 of coords.
      [~, val{ii}] = ismember(nodeCoords, roiCoords, 'rows');
    end

  case {'voxel2fibernodepairs','v2fn'}
    % voxel2FNpairs = fgGet(fgImg,'voxel 2 fibernode pairs',roiCoords);
    % voxel2FNpairs = fgGet(fgImg,'voxel 2 fibernode pairs',roiCoords,nodes2voxels);
    %
    % The return is a cell array whose size is the number of voxels.
    % The cell is a Nx2 matrix of the (fiber, node) pairs that pass
    % through it.
    %
    % The value N is the number of nodes in the voxel.  The first
    % column is the fiber number.  The second column reports the indexes
    % of the nodes for each fiber in each voxel.
    tic
    fprintf('\n[fgGet] Computing fibers/nodes pairing in each voxel...')
    if length(varargin) < 1, error('Requires the roiCoords.');
    else
      roiCoords = varargin{1};
      nCoords = size(roiCoords,1);
    end
    if length(varargin) < 2
      % We assume the fg and the ROI coordinates are in the same
      % coordinate frame.
      nodes2voxels    = fgGet(fg,'nodes 2 voxels',roiCoords);
    else nodes2voxels = varargin{2};
    end

    nFibers      = fgGet(fg,'nFibers');
    voxelsInFG   = fgGet(fg,'voxels in fg',nodes2voxels);
    roiNodesInFG = fgGet(fg,'nodes in voxels',nodes2voxels);
    val = cell(1,nCoords);
    for thisFiber=1:nFibers
      voxelsInFiber = voxelsInFG{thisFiber};   % A few voxels, in a list
      nodesInFiber  = roiNodesInFG{thisFiber}; % The corresponding nodes

      % Then add a row for each (fiber,node) pairs that pass through
      % the voxels for this fiber.
      for jj=1:length(voxelsInFiber)
        thisVoxel = voxelsInFiber(jj);
        % Print out roi coord and fiber coord to verify match
        % roiCoords(thisVoxel,:)
        % fg.fibers{thisFiber}(:,nodesInFiber(jj))
        % Would horzcat be faster?
        val{thisVoxel} = cat(1,val{thisVoxel},[thisFiber,nodesInFiber(jj)]);
      end
    end
    fprintf('process completed in: %2.3fs.\n',toc)

  case {'nodesinvoxels'}
    % nodesInVoxels = fgGet(fg,'nodes in voxels',nodes2voxels);
    %
    % This cell array is a modified form of nodes2voxels (see above).
    % In that cell array every node in every fiber has a number
    % referring to its row in roiCoords, or a 0 when the node is not in
    % any roiCoord voxel.
    %
    % This cell array differs only in that the 0s removed.  This
    % is used to simplify certain calculations.
    %
    if length(varargin) <1
      error('Requires nodes2voxels cell array.');
    end

    nodes2voxels = varargin{1};
    nFibers = fgGet(fg,'nFibers');
    val = cell(1,nFibers);

    % For each fiber, this is a list of the nodes that pass through
    % a voxel in the roiCoords
    for ii = 1:nFibers
      % For each fiber, this is a list of the nodes that pass through
      % a voxel in the roiCoords
      lst = (nodes2voxels{ii} ~= 0);
      val{ii} = find(lst);
    end

  case 'voxelsinfg'
    % voxelsInFG = fgGet(fgImg,'voxels in fg',nodes2voxels);
    %
    % A cell array length n-fibers. Each cell has a list of the voxels
    % (rows of roiCoords) for a fiber.
    %
    % This routine eliminates the 0's in the nodes2voxels lists.
    %
    if length(varargin) < 1, error('Requires nodes2voxels cell array.'); end

    nodes2voxels = varargin{1};
    nFibers = fgGet(fg,'nFibers');
    val = cell(1,nFibers);
    for ii = 1:nFibers
      % These are the nodes that pass through a voxel in the
      % roiCoords
      lst = (nodes2voxels{ii} ~= 0);
      val{ii} = nodes2voxels{ii}(lst);
    end

  case {'voxels2fibermatrix','v2fm'}
    %   v2fm = fgGet(fgImg,'voxels 2 fiber matrix',roiCoords);
    % Or,
    %   v2fnPairs = fgGet(fgImg,'v2fn',roiCoords);
    %   v2fm = fgGet(fgImg,'voxels 2 fiber matrix',roiCoords, v2fnPairs);
    %
    % mrvNewGraphWin; imagesc(v2fm)
    %
    % Returns a binary matrix of size Voxels by Fibers.
    % When voxel ii has at least one node from fiber jj, there is a one
    % in v2fm(ii,jj).  Otherwise, the entry is zero.
    %

    % Check that the fg is in the image coordspace:
    if isfield(fg, 'coordspace') && ~strcmp(fg.coordspace, 'img')
      error('Fiber group is not in the image coordspace, please xform');
    end

    if isempty(varargin), error('roiCoords required');
    else
      roiCoords = varargin{1};
      nCoords   = size(roiCoords,1);
      if length(varargin) < 2
        v2fnPairs = fgGet(fg,'v2fn',roiCoords);
      else
        v2fnPairs = varargin{2};
      end
    end

    % Allocate matrix of voxels by fibers
    val = zeros(nCoords,fgGet(fg,'n fibers'));

    % For each coordinate, find the fibers.  Set those entries to 1.
    for ii=1:nCoords
      if ~isempty(v2fnPairs{ii})
        f = unique(v2fnPairs{ii}(:,1));
      end
      val(ii,f) = 1;
    end

  case {'fibersinroi','fginvoxels','fibersinvoxels'}
    % fList = fgGet(fgImg,'fibersinroi',roiCoords);
    %
    % v2fn = fgGet(fgImg,'v2fn',roiCoords);
    % fList = fgGet(fgImg,'fibersinroi',roiCoords,v2fn);
    %
    % Returns an integer vector of the fibers with at least
    % one node in a region of interest.
    %
    % The fg and roiCoords should be in the same coordinate frame.
    %
    if isempty(varargin), error('roiCoords required');
    elseif length(varargin) == 1
      roiCoords = varargin{1};
      v2fnPairs = fgGet(fg,'v2fn',roiCoords);
    elseif length(varargin) > 1
      roiCoords = varargin{1};
      v2fnPairs = varargin{2};
    end

    val = []; nCoords = size(roiCoords,1);
    for ii=1:nCoords
      if ~isempty(v2fnPairs{ii})
        val = cat(1,val,v2fnPairs{ii}(:,1));
      end
    end
    val = sort(unique(val),'ascend');

  case {'coordspace','fibercoordinatespace','fcspace'}
    % In some cases, the fg might contain information telling us in which
    % coordinate space its coordinates are set. This information is set
    % as a struct. Each entry in the struct can be either a 4x4 xform
    % matrix from the fiber coordinates to that space (with eye(4) for
    % the space in which the coordinates are defined), or (if the xform
    % is not know) an empty matrix.

    cspace_fields = fields(fg.coordspace);
    val = [];
    for f=1:length(cspace_fields)
      this_field = cspace_fields{f};
      if isequal(getfield(fg.coordspace, this_field), eye(4))
        val = this_field;
      end
    end

  otherwise
    error('Unknown fg parameter: "%s"\n',param);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [I,O] = tone2interval(T,key)
% MUSIC.TONE2INTERVAL Returns the interval and octave of notes in a key.
%    I = MUSIC.TONE2INTERVAL(T) returns the interval at which the tones in T are
%    found in the key of 'C'. T is a vector of semitones defined relative to C4,
%    where C4 corresponds to tone 0.
%
%    I = MUSIC.TONE2INTERVAL(T,KEY) uses the key of KEY. KEY may be an interval
%    offset from 'C', or a character note (e.g., 'A', 'F#').
%
%    [I,O] = MUSIC.TONE2INTERVAL(...) also returns the octave number of each
%    tone.
%
%    Examples
%       I     = music.tone2interval([7 19])      % returns [7 7]
%       [I,O] = music.tone2interval([7 19])      % returns [7 7], [4 5]
%       [I,O] = music.tone2interval([7 19],'G')  % returns [0 0], [4 5]
%
%    See also music.tone2freq, music.tone2note,  music.interval2tone.

%    Author: E. Johnson
%    Copyright 2010 The MathWorks, Inc.

if nargin < 2
    key = 0;
end
if ischar(key)
    key = music.note2interval(key);
end

I = mod(T,12) - key;
I = mod(I,12);

O = 4 + floor(T / 12);

author: Sable

repo: mcbench-benchmarks


text:
function D = EuDist2(fea_a,fea_b,bSqrt)
%EUDIST2 Efficiently Compute the Euclidean Distance Matrix by Exploring the
%Matlab matrix operations.
%
%   D = EuDist(fea_a,fea_b)
%   fea_a:    nSample_a * nFeature
%   fea_b:    nSample_b * nFeature
%   D:      nSample_a * nSample_a
%       or  nSample_a * nSample_b
%
%    Examples:
%
%       a = rand(500,10);
%       b = rand(1000,10);
%
%       A = EuDist2(a); % A: 500*500
%       D = EuDist2(a,b); % D: 500*1000
%
%   version 2.1 --November/2011
%   version 2.0 --May/2009
%   version 1.0 --November/2005
%
%   Written by Deng Cai (dengcai AT gmail.com)

if ~exist('bSqrt','var')
    bSqrt = 1;
end

if (~exist('fea_b','var')) || isempty(fea_b)
    aa = sum(fea_a.*fea_a,2);
    ab = fea_a*fea_a';

    if issparse(aa)
        aa = full(aa);
    end

    D = bsxfun(@plus,aa,aa') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
    D = max(D,D');
else
    aa = sum(fea_a.*fea_a,2);
    bb = sum(fea_b.*fea_b,2);
    ab = fea_a*fea_b';

    if issparse(aa)
        aa = full(aa);
        bb = full(bb);
    end

    D = bsxfun(@plus,aa,bb') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [orderedPoints,biggest,goodPerimEdges]=orderMeshPerimeterPointsAll_new(mesh,perimeterEdges)
%
% orderedPoints=orderMeshPerimeterPointsAll(perimeterEdges)
%
% AUTHOR:  Wade
% DATE: 062100
% PURPOSE:
% Takes a list of edges defined as point pairs..
% Returns the points in the list in their connection order...
% If more than one (unlinked) perimeter exists, return a structure containing all the separate ordered perimeters
%
%  N.B.  Each vertex should appear in exactly two rows.
%  Starting with an initial position on the perimeter, startNode, the routine finds
%  the other node in the row. By selectively deleting entries in the connection matrix
%  as we progress, we can use a simple find operation each time to give us the next node.
%  ARW 031401 : This routine is critical to the success of the opration. It will break when
%  a perimeter point is connected to more than two other perimeter points. How can this happen? (you ask).
%  Well it can happen when two perimeters share a common point - like interlocking rings. In theory, I think you could have two perimeters sharing many
%  points. We deal with this in some sense in findGroupPerimeter when we eliminate nodes that are not part of any face but the 'just touching' condition
%  is a headache. Here's one solution...
%  The perimeter connection matrix should have 2 entries per row. If it doesn't then iteratively generate new ones with fewer connections until you have a
% (large?) set of perimeter connection matrices with 2 entries per row. E.G. if row t has entries at x,y,z then generate three copies of the conmat,
% one with entries at x,y one with y,z and one with x,z.
% If it has four entries, then you generate  nchoosek(1:4,2) new matrices.
% (You can see how this might mount up if we have many bad points -  if you have two bad points, each generating
% n conMats then you end up with n^2 conmats - fortunately they seem to be rare
% and anyway, connection matrices are sparse and perimeter connection matrices are small and sparse.
% Then just run the regular perimeter detection algorithm on all the conmats and pick the biggest final perimeter.

nVerts=length(mesh.connectionMatrix); % Total number of vertices in the full connection matrix

edgePoints=unique(perimeterEdges(:)); % List of nodes that are on the perimeter
[y x]=size(mesh.connectionMatrix); % = nVerts

% Generate a connection matrix using only the edges - make sure it's symmetric
ec=sparse(perimeterEdges(:,1),perimeterEdges(:,2),1,nVerts,nVerts);
ec=ec+ec';
edgeConMat=sparse((ec~=0));
nnz(edgeConMat)
[eY,eX]=size(edgeConMat);

% Now do something interesting: go through the edge connecton matrix and find all the rows that have more than 2 entries
% Split these up and do some record keeping - we generate new nodes
% Used for checking which perimeter is biggest
nConnects=sum(edgeConMat,2); % How many points on each row?

badPoints=find(nConnects>2);
nBadPoints=length(badPoints);
fprintf('There are %d bad points',nBadPoints);
conMatStack{1}=edgeConMat;

for thisBadPoint=1:nBadPoints
    fprintf('\nLegalizing bad point number %d',thisBadPoint);

	conMatStack=legaliseConMatRow(conMatStack,badPoints(thisBadPoint));
end

nConMats=length(conMatStack);
fprintf('\n%d different edge connection matrices',nConMats);

biggest=-Inf;
biggestSize=-Inf;
goodEdgeConMat=-1;

thisPerim=1; % Index of the current perimeter
nextRow=-99999; % Dummy

foundAllPerimeters=0; % Flag to say we've finished

% Now loop over all the connection matrices
for thisConMat=1:nConMats
	edgeConMat=conMatStack{thisConMat};
	nnz(edgeConMat)
	while (~foundAllPerimeters)
   counter=1;
   [startRow startCol]=find(edgeConMat);
   if (~isempty(startRow)) % If there is one...

      startRow=startRow(1); % ...pick the first non-zero row in the edge connection matrix
      startCol=startCol(1); % Pick the first non-zero column.
      orderedPoints{thisPerim}.points=startRow;

      thisRow=startRow;
      deleteCol=startCol;
      nextRow=-Inf;
      foundEnd=0;
      edgeConMat(thisRow,deleteCol)=0; % Zero one of the entries in this row

      while ((~foundEnd)&(counter<10000))

         % Do the deletions
         edgeConMat(thisRow,deleteCol)=0; % Zero one of the entries in this row
         edgeConMat(deleteCol,thisRow)=0; % Zero one of the entries in this row

         [thisCol]=find(edgeConMat(thisRow,:));

         thisCol=unique(thisCol(:));
         if (length(thisCol)>1)
            disp ('Warning!! - thiscol has more than 1 entry:');
            disp (thisCol);
            disp ('Perim node is connected to many other perim nodes. Carrying on but expect trouble....');
            % Can we just ignore superfluous ones? We'll press ahead regardless...
         	thisCol=thisCol(1);
         end

         if (~isempty(thisCol)) % If there >is< another point in this perimeter
            orderedPoints{thisPerim}.points=[orderedPoints{thisPerim}.points;thisCol];
             deleteCol=thisRow;
             thisRow=thisCol;

         else % We've found the end of the perimeter.
           % fprintf ('\n.');

            foundEnd=1;
            orderedPoints{thisPerim}.size=counter;

         end % endif
         counter=counter+1;

      end % end while not startRow

      if (orderedPoints{thisPerim}.size>biggestSize) % Have to return the index of the larges tperim
         biggestSize=orderedPoints{thisPerim}.size;
         biggest=thisPerim;
		 goodEdgeConMat=conMatStack{thisConMat};
      end % end size check

      fprintf('\n%d points');

   end % End if anything found to start with

     foundAllPerimeters=(nnz(edgeConMat)==0);
     thisPerim=thisPerim+1;
	end % End while not all perims found
end % Next con mat

fprintf('\nFound %d perimeters...',thisPerim-1);
[goodY,goodX]=find(triu(goodEdgeConMat));
goodPerimEdges=[goodY(:),goodX(:)];


text:
%
Nfilt = size(rez.W,2);
ik = ceil(rand * Nfilt);
% ik = 1;
%
sp = find(rez.st3(:,2)==ik);
sp = sort(sp);

st = rez.st3(sp, 1);

clp = rez.cProjPC(sp, :, :);
clp = clp - mean(clp,1);
% clp = clp - my_conv2(clp, 50, 1);

nspikes = size(clp,1);
[u s v] = svdecon(clp(:,:)');

figure(1)

plotmatrix(v(:,1:4), '.')

figure(2)
vplot = clp(:,:) * u(:,1:2);
plot(st, my_conv2(vplot(:,1:2), 1, 1))

figure(3)
clp = rez.cProjPC(sp, :, :);
clp = clp - mean(clp,1);

nspikes = size(clp,1);
[u s v] = svdecon(clp(:,:)');

cls = (clp(:,:) * u(:,1:4)) * u(:,1:4)';

clp = rez.cProjPC(sp, :, :);
clp = clp - mean(clp,1);

for i = 1:16
    subplot(4,4,i)
    l = randperm(size(clp,1), 2);
    w = cls(l(1), :) - cls(l(2), :);
    w = zscore(w);

    Y = clp(:,:) * w';

    hist(Y, 100)
end
%%

hist(Y, 100)

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function y = inverseGammaPdf(x,alpha,beta)
% The inverse Gamma PDF.
% INPUT
% x  (>0)
% alpha - shape parameter
% beta  - scale parameter
% OUTPUT
% y = the pdf [nrX nrBeta]
% BK - 2018

[nrX,nrBetaInX] = size(x);
[nrBeta,nrXInBeta] = size(beta);
if nrXInBeta<nrBetaInX
    beta =repmat(beta,[1 nrBetaInX]);
end

if nrX==1 && nrBeta>1
    x = repmat(x,[1 nrBeta]);
end

z = x<0;
x(z) =NaN;
y = (beta.^alpha)./gamma(alpha).*(1./x).^(alpha+1).*exp(-beta./x);
y(z) = 0;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [occGrad,occlusionErrors] = occlusionGradOptimal(S,state,p2d)

%mask = double(state.mask);
[Y,X]=size(state.mask);
%mkp = mean(state.kps,2);% we had shifted everything by this when computing projected points so shift back now
%occlusionErrors = zeros(numpoints,1);

%% Finding points projected outside the image
%points2d = state.cameraScale*state.cameraRot*S'; %% Verify this later
%points2d = points2d' + repmat(mkp',size(S,1),1);
%points2d = round(points2d(:,1:2)); %% Verify correctness later
if(nargin<3)
    p2d = transform2d(S,state.cameraRot,state.cameraScale,state.translation);
    p2d = round(p2d); %% Verify correctness later
end

badX = (p2d(:,1)>X);p2d(badX,1)=X;
badX = (p2d(:,1)<1);p2d(badX,1)=1;
badY = (p2d(:,2)>Y);p2d(badY,2)=Y;
badY = (p2d(:,2)<1);p2d(badY,2)=1;

%% Computing nearest silhoutte point and gradient
%tic
diff3d = zeros(size(S));
IDX_2 = state.occlusionIDX;
sil_nbridx = IDX_2(sub2ind([Y X],p2d(:,2),p2d(:,1)));
[Ys,Xs] = ind2sub([Y X],sil_nbridx);
diff2d = double([Xs Ys]) - p2d;
diff3d(:,1:2)=diff2d;

occlusionErrors = double(sum(abs(diff3d),2)>=2);
occGrad = 1/state.cameraScale*diff3d*state.cameraRot;

end

repo: CategoryShapes


text:
function f = flops_solve_tri(T,b,c)
% FLOPS_SOLVE_TRI   Flops for triangular left division.
% FLOPS_SOLVE_TRI(T,b) returns the number of flops for solve_tri(T,b).
% FLOPS_SOLVE_TRI(n,n,m) returns the number of flops for
% solve_tri(eye(n),ones(n,m)).

if nargin == 2
  f = flops_solve_tri(rows(T),cols(T),cols(b));
  return;
end
% number of multiplies+adds is
% sum(i=1..n) sum(k=i-1..1) 2 = sum(i=1..n) 2*(i-1) = n^2-n
% number of divides is n
f = (T*b + b*(flops_div-1))*c;

author: Cloud-CV

repo: object-proposals


text:
function spdemovarout
% SPDEMOVAROUT   Sparse interpolation demo for multiple outputs
%    Sparse grid interpolation demo for a function with multiple
%    output parameters and non-vectorized processing. In this case,
%    we use the Maximum-norm-based grid.

% Author : Andreas Klimke, Universitaet Stuttgart
% Date   : January 27, 2004
% Version: 1.3

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

% Define problem dimension
d = 2;

% Define number of output parameters of fvarout
nout = 3;

% Get options structure for sparse interpolation
options = spset('NumberOfOutputs', nout, ...
								'GridType', 'Maximum');

% Create full grid for plotting
gs = 33;
[X,Y] = meshgrid(linspace(0,2,gs),linspace(-1,1,gs));

% Compute sparse grid weights over domain [0,2]x[-1,1]
z = spvals(@fvarout, d, [0 2; -1 1], options, nout);

% Compute inpterpolated values at full grid
for k = 1:nout
	z.selectOutput = k;
	ip{k} = spinterp(z, X, Y);
end

% Plot interpolated results
for k = 1:nout
	subplot(1,nout,k);
	mesh(X, Y, ip{k});
	title(['interpolated out',num2str(k)]);
end

%------------------------------------------------------------------
function [varargout] = fvarout(x,y,nout)
% FVARARGOUT    function with multiple output arguments, the number
% of output arguments is determined by the parameter nout.

out = zeros(1,nout);
out(1:nout) = 1./((x-(1:nout)/nout).^2+(y)^2+(1:nout)./nout);
varargout = num2cell(out);

author: johannesgerer

repo: jburkardt-m


text:
% POP_CHANCENTER - recenter cartesian X,Y,Z channel coordinates
%
% Usage:
%    >> chanlocs = pop_chancenter(chanlocs); % pop up interactive window
%    >> [chanlocs centerloc] = pop_chancenter(chanlocs, center, omitchan);
%
% Inputs:
%    chanlocs  = eeglab channel location structure (see READLOCS)
%    center    = [X Y Z] known center different from [0 0 0]
%                [] will optimize the center location according
%                to the best sphere. Default is [0 0 0].
%    omitchan  = indices of channel to omit when computing center
%
% Outputs:
%    chanlocs  = updated channel location structure
%    centerloc = 3-D location of the new center (in the old coordinate
%                frame).
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, Feb 2004
%
% See also: CHANCENTER, SPHERROR, CART2TOPO

% Copyright (C) 2004, Arnaud Delorme, SCCN/INC/UCSD, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [ chanlocs, newcenter, com] = pop_chancenter( chanlocs, center, omitchans)

optim = 0;

if nargin<1
    help pop_chancenter
    return;
end

com = '';
newcenter = [];
if nargin < 3
    omitchans = [];
end
if nargin < 2
    cb_browse = [ 'tmpchans = get(gcbf, ''userdata'');' ...
                  'set(findobj(gcbf, ''tag'', ''chans''), ''string'', ' ...
                  'int2str(pop_chansel( { tmpchans.labels } )));' ];
    cb_opt    = [ 'if get(gcbo, ''value''), ' ...
                  '    set(findobj(gcbf, ''tag'', ''center''), ''enable'', ''off'');' ...
                  'else,' ...
                  '    set(findobj(gcbf, ''tag'', ''center''), ''enable'', ''on'');' ...
                  'end;' ];
    geometry = { [1.3 0.28 1 1] [1] [1] [2 1] };
    uilist = { { 'Style', 'text',       'string', 'Optimize center location', 'fontweight', 'bold'   } ...
			   { 'Style', 'checkbox',   'value',  1  'callback' cb_opt } ...
               { 'Style', 'text',       'string', 'or specify center', 'fontweight', 'bold'  } ...
               { 'Style', 'edit',       'string', '0 0 0', 'tag' 'center' 'enable' 'off' } ...
               { } ...
			   { 'Style', 'text',       'string', 'Channel indices to ignore for best-sphere matching'  } ...
			   { 'Style', 'edit',       'string', '', 'tag', 'chans' } ...
			   { 'Style', 'pushbutton', 'string', 'Browse', 'callback', cb_browse } };

    results = inputgui( geometry, uilist, 'pophelp(''pop_chancenter'');', ...
                        'Convert channel locations -- pop_chancenter()', chanlocs );
	if isempty(results), return; end
	if results{1}
        center = [];
    else
        center  = eval( [ '[' results{2} ']' ] );
    end
    if ~isempty(results{3})
        omitchans =  eval( [ '[' results{3} ']' ] );
    end
end

% remove channels
% ---------------
c = setdiff_bc([1:length(chanlocs)], union(omitchans, find(cellfun('isempty', { chanlocs.theta }))));

% optimize center
% ---------------
[X, Y, Z, newcenter]= chancenter( [ chanlocs(c).X ]', [ chanlocs(c).Y ]', [ chanlocs(c).Z ]', center);
for index = 1:length(c)
    chanlocs(c(index)).X  = X(index);
    chanlocs(c(index)).Y  = Y(index);
    chanlocs(c(index)).Z  = Z(index);
end
disp('Note: automatically convert XYZ coordinates to spherical and polar');
chanlocs = convertlocs(chanlocs, 'cart2all');
if ~isempty(omitchans)
    disp('Important warning: the location of omitted channels has not been modified');
end
if nargout > 2
    com = sprintf('%s = pop_chancenter( %s, %s);', inputname(1), inputname(1), vararg2str({ center omitchans }));
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [meshr,Ir] = mshMidpoint(mesh,I)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshMidpoint.m                                 |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 25.11.2018                                    |
%| ( === ) |   SYNOPSIS   : Refine triangular mesh with midpoint algorithm|
%|  `---'  |                                                              |
%+========================================================================+

% Check dimenion
if (size(mesh,2) ~= 3)
    error('mshMidpoint : unavailable case 1')
end

% Save color and replace by hierarchy
col      = mesh.col;
mesh.col = (1:length(mesh))';

% Submeshing (triangle)
meshs = mesh.sub(I);

% Edge meshes
edgs       = meshs.edg;
[edg,Itri] = mesh.edg;

% Interface with edge multiplicity for triangle
[int,Iedg] = intersect(edg,edgs);
ind        = ismember(Itri,Iedg);
mlt        = sum(ind,2);

% Security
tmp = setdiff(int,edgs);
if (size(tmp.elt,1) ~= 0)
    error('mshMidpoint.m : unavailable case 2');
end

% Initialize refined mesh for element without refinement
meshr = mesh.sub(mlt==0);

% Subdivision with 1 common edge
tmp   = mesh.sub(mlt==1);
tmp   = mshMidpoint1(tmp,int);
meshr = union(meshr,tmp);

% Subdivision with 2 common edges
tmp   = mesh.sub(mlt==2);
tmp   = mshMidpoint2(tmp,int);
meshr = union(meshr,tmp);

% Subdivision with 3 common edges
tmp   = mesh.sub(mlt==3);
tmp   = mshMidpoint3(tmp);
meshr = union(meshr,tmp);

% Parent indices and replace colours
Ir        = meshr.col;
meshr.col = col(Ir);

% Security
if (sum(mesh.ndv)-sum(meshr.ndv))/sum(mesh.ndv) > 1e-15*length(meshr)
    error('mshMidpoint.m : unavailable case 3');
end
end

function mesh = mshMidpoint1(mesh,int)
% Mesh nodes and edges center
[nds,ctr] = data(mesh);

% Interface center
Xctr = int.ctr;

% Refined mesh initialization
Nvtx = size(mesh.vtx,1);
Nelt = length(mesh);
col  = mesh.col;
mesh = mesh.sub([]);

% Loop over nodes
for i = 1:3
    % Neighbours
    ip1 = mod(i,3)+1;
    ip2 = mod(ip1,3)+1;

    % Selected center are inside subdivided mesh
    I = find(ismember(single(ctr{i}),single(Xctr),'rows'));

    % First elements
    vtx  = [nds{i}(I,:) ; nds{ip1}(I,:) ; ctr{i}(I,:)];
    elt  = reshape((1:3*length(I))',length(I),3);
    tmp  = msh(vtx,elt,col(I));
    mesh = union(mesh,tmp);

    % Second elements
    vtx  = [nds{i}(I,:) ; ctr{i}(I,:) ; nds{ip2}(I,:)];
    elt  = reshape((1:3*length(I))',length(I),3);
    tmp  = msh(vtx,elt,col(I));
    mesh = union(mesh,tmp);

    % Mesh fusion with previous submeshes
    mesh = union(mesh,tmp);
end

% Security
if size(mesh.elt,1) ~= 2*Nelt
    error('mshMidpoint1.m : unavailable case 1')
end
if size(mesh.vtx,1) ~= Nvtx+Nelt
    error('mshMidpoint1.m : unavailable case 2')
end
end

function mesh = mshMidpoint2(mesh,int)
% Mesh nodes and edges center
[nds,ctr] = data(mesh);

% Interface center
Xctr = int.ctr;

% Refined mesh initialization
Nvtx = size(mesh.vtx,1);
Nelt = length(mesh);
col  = mesh.col;
mesh = mesh.sub([]);

% Loop over nodes
for i = 1:3
    % Neighbours
    ip1 = mod(i,3)+1;
    ip2 = mod(ip1,3)+1;

    % Selected nodes center not inside subdivided mesh
    I = find(~ismember(single(ctr{i}),single(Xctr),'rows'));

    % First elements
    vtx  = [nds{i}(I,:) ; ctr{ip2}(I,:) ; ctr{ip1}(I,:)];
    elt  = reshape((1:3*length(I))',length(I),3);
    tmp  = msh(vtx,elt,col(I));
    mesh = union(mesh,tmp);

    % Second elements
    vtx  = [nds{ip1}(I,:) ; nds{ip2}(I,:) ; ctr{ip1}(I,:)];
    elt  = reshape((1:3*length(I))',length(I),3);
    tmp  = msh(vtx,elt,col(I));
    mesh = union(mesh,tmp);

    % Third elements
    vtx  = [nds{ip1}(I,:) ; ctr{ip1}(I,:) ; ctr{ip2}(I,:) ; ];
    elt  = reshape((1:3*length(I))',length(I),3);
    tmp  = msh(vtx,elt,col(I));
    mesh = union(mesh,tmp);
end

% Security
if size(mesh.elt,1) ~= 3*Nelt
    error('mshMidpoint2.m : unavailable case 1')
end
if size(mesh.vtx,1) ~= Nvtx+2*Nelt
    error('mshMidpoint2.m : unavailable case 2')
end
end

function mesh = mshMidpoint3(mesh)
% Mesh nodes and edges center
[nds,ctr] = data(mesh);

% Refined mesh initialization with centered triangle
Nelt = length(mesh);
vtx  = [ctr{1} ; ctr{2} ; ctr{3}];
elt  = reshape((1:3*Nelt)',Nelt,3);
col  = mesh.col;
mesh = msh(vtx,elt,col);

% For each node
for i = 1:3
    % Neighbours
    ip1 = mod(i,3)+1;
    ip2 = mod(ip1,3)+1;

    % New submesh with nodes triangles
    vtx = [nds{i} ; ctr{ip2} ; ctr{ip1}];
    tmp = msh(vtx,elt,col);

    % Mesh fusion with previous submeshes
    mesh = union(mesh,tmp);
end

% Security
if size(mesh.elt,1) ~= 4*Nelt
    error('mshMidpoint3.m : unavailable case 1')
end
end

function [nds,ctr] = data(mesh)
% Triangle nodes
nds{1} = mesh.vtx(mesh.elt(:,1),:);
nds{2} = mesh.vtx(mesh.elt(:,2),:);
nds{3} = mesh.vtx(mesh.elt(:,3),:);

% Edges center
ctr{1} = 0.5 * (nds{2} + nds{3});
ctr{2} = 0.5 * (nds{3} + nds{1});
ctr{3} = 0.5 * (nds{1} + nds{2});
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [outStr] = paddStr(str, desLength)
%paddStr Summary of this function goes here
%   Detailed explanation goes here
    outStr = str;
    while(length(outStr) < desLength)
        outStr = [outStr, ' '];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef PTKAdmin < MimScript
    % PTKAdmin. Script for running functions from the PTK object
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        InterfaceVersion = '1'
        Version = '1'
        Category= 'Admin'
    end

    methods (Static)
        function output = RunScript(ptk_obj, reporting, varargin)
            ptk_obj.(varargin{1})(varargin{2:end});
            output = [];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function value = octahedron_unit_volume_nd ( n )

%*****************************************************************************80
%
%% OCTAHEDRON_UNIT_VOLUME_ND returns the volume of the unit octahedron in ND.
%
%  Integration region:
%
%    Points X(1:N) such that:
%
%      Sum ( Abs ( X(1:N) ) ) <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the space.
%
%    Output, real OCTAHEDRON_UNIT_VOLUME_ND, the volume of
%    the unit octahedron.
%
  value = 2^n / prod ( 1 : n );

  return
end


text:
function [mem]=tt_mem(tt)
% Memory required to store the tensor TT
%   [MEM]=TT_MEM(TT) Computes the total number of elements in the TT-tensor
%   in the TT1.0 format. Please avoid its usage: it will be removed in
%   future releases. Use mem() from the object-oriented version.
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
mem=0;
d=size(tt,1);
for k=1:d
    mem=mem+numel(tt{k});
end
end

author: oseledets

repo: TT-Toolbox


text:
function pdf = cardioid_pdf ( x, a, b )

%*****************************************************************************80
%
%% CARDIOID_PDF evaluates the Cardioid PDF.
%
%  Discussion:
%
%    The cardioid PDF can be thought of as being applied to points on
%    a circle.  Compare this distribution with the "Cosine PDF".
%
%    PDF(A,B;X) = ( 1 / ( 2 * PI ) ) * ( 1 + 2 * B * COS ( X - A ) )
%    for 0 <= B <= 1/2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    N I Fisher,
%    Statistical Analysis of Circular Data,
%    Cambridge, 1993.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 <= B <= 0.5.
%
%    Output, real PDF, the value of the PDF.
%
  pdf = ( 1.0 + 2.0 * b * cos ( x - a ) ) / ( 2.0 * pi );

  return
end


text:
%THIS FUMCTION ILLUSTRATES HOW TO USE THE USER-DEFINED ALGORITHM
function [AH,XH] = user_alg1(Y,r,X)
%
% The example of implementing the user-defined algorithm (this is the regularized Fixed-Point algorithm based)
%
% INPUTS:
% Y - mixed signals (matrix of size [m by T])
% r - number of estimated signals
% X - true source signals
%
% OUTPUTS
% AH - estimated mixing matrix (matrix of size [m by r])
% XH - estimated source signals (matrix of size [r by T])
%
% #########################################################################
% Initialization
[m,T]=size(Y);
Y(Y <=0) = eps; % this enforces the positive value in the data
AH=rand(m,r);
XH=rand(r,T);

IterNo = 1000; % number of alternating steps

% Iterations
for k = 1:IterNo

    alpha_reg = 20*exp(-k/10); % regularization parameter
    XH = max(1E6*eps,pinv(AH'*AH +  alpha_reg)*AH'*Y);
    AH = max(1E6*eps, Y*XH'*pinv(XH*XH' + alpha_reg));
    AH = AH*diag(1./(sum(AH,1) + eps));

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% WAVELET_LAYER_2D Compute the wavelet transform of a scattering layer
%
% Usage
%    [U_phi, U_psi] = WAVELET_LAYER_2D(U, filters, options)
%
% Input
%    U (struct): input scattering layer
%    filters (struct): filter bank
%    options (struct): same as wavelet_2d
%
% Output
%    A (struct): Averaged wavelet coefficients
%    V (struct): Wavelet coefficients of the next layer
%
% Description
%    This function has a pivotal role between WAVELET_2D (which computes a
%    single wavelet transform), and WAVELET_FACTORY_2D (which creates the
%    whole cascade). Given inputs modulus wavelet coefficients
%    corresponding to a layer, WAVELET_LAYER_2D computes the wavelet
%    transform coefficients of the next layer using WAVELET_2D.
%
% See also
%   WAVELET_2D, WAVELET_FACTORY_2D, WAVELET_LAYER_1D

function [U_phi, U_psi] = wavelet_layer_2d(U, filters, options)

    calculate_psi = (nargout>=2); % do not compute any convolution
    % with psi if the user does get U_psi

    if ~isfield(U.meta,'theta')
        U.meta.theta = zeros(0,size(U.meta.j,2));
    end

    if ~isfield(U.meta, 'resolution'),
        U.meta.resolution = 0;
    end

    p2 = 1;
    for p = 1:numel(U.signal)
        x = U.signal{p};
        if (numel(U.meta.j)>0)
            j = U.meta.j(end,p);
        else
            j = -1E20;
        end

        % compute mask for progressive paths
        options.psi_mask = calculate_psi & ...
            (filters.psi.meta.j >= j + filters.meta.Q);

        % set resolution of signal
        options.x_resolution = U.meta.resolution(p);

        % compute wavelet transform
        [x_phi, x_psi, meta_phi, meta_psi] = wavelet_2d(x, filters, options);

        % copy signal and meta for phi
        U_phi.signal{p} = x_phi;
        U_phi.meta.j(:,p) = [U.meta.j(:,p); filters.phi.meta.J];
        U_phi.meta.theta(:,p) = U.meta.theta(:,p);
        U_phi.meta.resolution(1,p) = meta_phi.resolution;

        % copy signal and meta for psi
        for p_psi = find(options.psi_mask)
            U_psi.signal{p2} = x_psi{p_psi};
            U_psi.meta.j(:,p2) = [U.meta.j(:,p);...
                filters.psi.meta.j(p_psi)];
            U_psi.meta.theta(:,p2) = [U.meta.theta(:,p);...
                filters.psi.meta.theta(p_psi)];
            U_psi.meta.resolution(1,p2) = meta_psi.resolution(p_psi);
            p2 = p2 +1;
        end

    end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
 function [sn,psdx,ff] = get_noise_fft(Y,options)

        defoptions = CNMFSetParms;
        if nargin < 2 || isempty(options); options = defoptions; end

        if ~isfield(options,'noise_range'); options.noise_range = defoptions.noise_range; end
        range_ff = options.noise_range;
        if ~isfield(options,'noise_method'); options.noise_method = defoptions.noise_method; end
        method = options.noise_method;
        if ~isfield(options,'block_size'); options.block_size = defoptions.block_size; end
        block_size = options.block_size;
        if ~isfield(options,'split_data'); options.split_data = defoptions.split_data; end
        split_data = options.split_data;
        if ~isfield(options,'max_timesteps') || isempty(options.max_timesteps);
            options.max_timesteps = defoptions.max_timesteps;
        end

        dims = ndims(Y);
        sizY = size(Y);
        N = min(sizY(end),options.max_timesteps);
        if N < sizY(end)
           %Y = reshape(Y,prod(sizY(1:end-1)),[]);
           switch ndims(Y),
               case 2,
                    Y(:,N+1:end) = [];
               case 3,
                    Y(:,:,N+1:end) = [];
               case 4,
                    Y(:,:,:,N+1:end) = [];
           end
        end

        Fs = 1;
        ff = 0:Fs/N:Fs/2;
        indf=ff>range_ff(1);
        indf(ff>range_ff(2))=0;
        if dims > 1
            d = prod(sizY(1:dims-1));
            Y = reshape(Y,d,N);
            Nb = prod(block_size);
            SN = cell(ceil(d/Nb),1);
            PSDX = cell(ceil(d/Nb),1);
            if ~split_data
                for ind = 1:ceil(d/Nb);
                    xdft = fft(Y((ind-1)*Nb+1:min(ind*Nb,d),:),[],2);
                    xdft = xdft(:,1: floor(N/2)+1); % FN: floor added.
                    psdx = (1/(Fs*N)) * abs(xdft).^2;
                    psdx(:,2:end-1) = 2*psdx(:,2:end-1) + eps;
                    %SN{ind} = mean_psd(psdx(:,indf),method);
                    switch lower(method)
                        case 'mean'
                            SN{ind}=sqrt(mean(psdx(:,indf)/2,2));
                        case 'median'
                            SN{ind}=sqrt(median(psdx(:,indf)/2,2));
                        case 'logmexp'
                            SN{ind} = sqrt(exp(mean(log(psdx(:,indf)/2),2)));
                        otherwise
                            error('unknown method for averaging noise..')
                    end
                    PSDX{ind} = psdx;
                end
            else
                nc = ceil(d/Nb);
                Yc = mat2cell(Y,[Nb*ones(nc-1,1);d-(nc-1)*Nb],N);
                parfor ind = 1:ceil(d/Nb);
                    xdft = fft(Yc{ind},[],2);
                    xdft = xdft(:,1:floor(N/2)+1);
                    psdx = (1/(Fs*N)) * abs(xdft).^2;
                    psdx(:,2:end-1) = 2*psdx(:,2:end-1) + eps;
                    Yc{ind} = [];
                    switch lower(method)
                        case 'mean'
                            SN{ind}=sqrt(mean(psdx(:,indf)/2,2));
                        case 'median'
                            SN{ind}=sqrt(median(psdx(:,indf)/2,2));
                        case 'logmexp'
                            SN{ind} = sqrt(exp(mean(log(psdx(:,indf)/2),2)));
                        otherwise
                            error('unknown method for averaging noise..')
                    end

                end
            end
            sn = cell2mat(SN);
        else
            xdft = fft(Y);
            xdft = xdft(:,1:floor(N/2)+1);
            psdx = (1/(Fs*N)) * abs(xdft).^2;
            psdx(:,2:end-1) = 2*psdx(:,2:end-1) + eps;
            switch lower(method)
                case 'mean'
                    sn = sqrt(mean(psdx(:,indf)/2,2));
                case 'median'
                    sn = sqrt(median(psdx(:,indf)/2,2));
                case 'logmexp'
                    sn = sqrt(exp(mean(log(psdx(:,indf)/2),2)));
                otherwise
                    error('unknown method for averaging noise..')
            end
        end
        psdx = cell2mat(PSDX);
        if dims > 2
            sn = reshape(sn,sizY(1:dims-1));
        end
 end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function ReinhardTMOv(hdrv, filenameOutput, tmo_alpha, tmo_white, tmo_gamma, tmo_quality, tmo_video_profile)
%
%
%      ReinhardTMOv(hdrv, filenameOutput, tmo_alpha, tmo_white, tmo_gamma, tmo_quality, tmo_video_profile)
%
%
%       Input:
%           -hdrv: a HDR video structure; use hdrvread to create a hdrv
%           structure
%           -filenameOutput: output filename (if it has an image extension,
%           single files will be generated)
%           -tmo_alpha:
%           -tmo_white:
%           -tmo_gamma: gamma for encoding the frame
%           -tmo_quality: the output quality in [1,100]. 100 is the best quality
%           1 is the lowest quality.%
%           -tmo_video_profile: the compression profile (encoder) for compressing the stream.
%           Please have a look to the profile of VideoWriter from the MATLAB
%           help. Depending on the version of MATLAB some profiles may be not
%           be present.
%
%     Copyright (C) 2013-15 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/>.
%
%     The paper describing this technique is:
%     "Real-time Automated Tone Mapping System for HDR Video"
% 	  by Chris Kiser, Erik Reinhard, Mike Tocci and Nora Tocci
%     in IEEE International Conference on Image Processing, 2012
%
%

if(~exist('tmo_alpha', 'var'))
    tmo_alpha = 0.5;
end

if(~exist('tmo_white', 'var'))
    tmo_white = 1e6;
end

if(~exist('tmo_gamma', 'var'))
    tmo_gamma = -2.2;
end

if(tmo_gamma < 0)
    bsRGB = 1;
else
    bsRGB = 0;
end

if(~exist('tmo_quality', 'var'))
    tmo_quality = 95;
end

if(~exist('tmo_video_profile', 'var'))
    tmo_video_profile = 'MPEG-4';
end

name = RemoveExt(filenameOutput);
ext  = fileExtension(filenameOutput);

bVideo = 0;
writerObj = 0;

if(strcmp(ext, 'avi') == 1 | strcmp(ext, 'mp4') == 1)
    bVideo = 1;
    writerObj = VideoWriter(filenameOutput, tmo_video_profile);
    writerObj.FrameRate = hdrv.FrameRate;
    writerObj.Quality = tmo_quality;
    open(writerObj);
end

hdrv = hdrvopen(hdrv, 'r');

disp('Tone Mapping...');

beta = 0.5;

for i=1:hdrv.totalFrames
    disp(['Processing frame ', num2str(i)]);
    [frame, hdrv] = hdrvGetFrame(hdrv, i);

    %Only physical values
    frame = RemoveSpecials(frame);
    frame(frame < 0) = 0;

    %computing statistics for the current frame
    L = lum(frame);
    Lav_cur = logMean(L);

    if(i == 1)
        Lav_prev = Lav_cur;
    end

    Lav = Lav_prev * (1 - beta) + beta * Lav_cur;

    %tone mapping
    [frameOut, ~, ~] = ReinhardTMO(frame, tmo_alpha, tmo_white, 'global', -1, Lav);

    %Gamma/sRGB encoding
    if(bsRGB)
        frameOut = ConvertRGBtosRGB(frameOut, 0);
    else
        frameOut = GammaTMO(frameOut, tmo_gamma, 0, 0);
    end

    frameOut = ClampImg(frameOut, 0.0, 1.0);

    if(bVideo)
        writeVideo(writerObj, frameOut);
    else
        imwrite(frameOut, [name, sprintf('%.10d',i), '.', ext]);
    end

    %updating for the next frame
    Lav_prev = Lav;
end
disp('OK');

if(bVideo)
    close(writerObj);
end

hdrvclose(hdrv);

end

author: banterle


text:
%% Set bounds for optimization problem

% Set Bounds
model_bounds = rabbit.getLimits();
model_bounds.states.x.lb = [-10, 0.5, -1, 2, 0.5, 2, 0.5];
model_bounds.states.x.ub = [10, 2, 1, 5, 2, 5, 2];
model_bounds.states.dx.lb = -20*ones(1,7);
model_bounds.states.dx.ub = 20*ones(1,7);
model_bounds.states.ddx.lb = -100*ones(1,7);
model_bounds.states.ddx.ub = 100*ones(1,7);
bounds = struct();

% Right Stance
bounds.RightStance = model_bounds;

bounds.RightStance.time.t0.lb = 0;
bounds.RightStance.time.t0.ub = 0;
bounds.RightStance.time.t0.x0 = 0;

bounds.RightStance.time.tf.lb = 0.25;
bounds.RightStance.time.tf.ub = 0.75;
bounds.RightStance.time.tf.x0 = 0.75;

bounds.RightStance.time.duration.lb = 0.25;
bounds.RightStance.time.duration.ub = 0.75;
bounds.RightStance.time.duration.x0 = 0.75;

bounds.RightStance.inputs.ConstraintWrench.fRightToe.lb = -1000;
bounds.RightStance.inputs.ConstraintWrench.fRightToe.ub = 1000;
bounds.RightStance.inputs.ConstraintWrench.fRightToe.x0 = 100;

bounds.RightStance.inputs.Control.u.lb = -100*ones(4,1);
bounds.RightStance.inputs.Control.u.ub = 100*ones(4,1);
bounds.RightStance.inputs.Control.u.x0 = zeros(4,1);

bounds.RightStance.params.pRightToe.lb = -0*ones(3,1);
bounds.RightStance.params.pRightToe.ub = 0*ones(3,1);
bounds.RightStance.params.pRightToe.x0 = zeros(3,1);

bounds.RightStance.params.atime.lb = -10*ones(6*4,1);
bounds.RightStance.params.atime.ub = 10*ones(6*4,1);
bounds.RightStance.params.atime.x0 = zeros(6*4,1);

bounds.RightStance.params.ptime.lb = [bounds.RightStance.time.tf.lb, bounds.RightStance.time.t0.lb];
bounds.RightStance.params.ptime.ub = [bounds.RightStance.time.tf.ub, bounds.RightStance.time.t0.ub];
bounds.RightStance.params.ptime.x0 = [bounds.RightStance.time.t0.x0, bounds.RightStance.time.tf.x0];

bounds.RightStance.time.kp = 100;
bounds.RightStance.time.kd = 20;

% Right Impact
bounds.RightImpact = model_bounds;

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function allen_ccf_npx_4shank_spherical(tv,av,st)
% allen_ccf_npx(tv,av,st)
% Andy Peters (peters.andrew.j@gmail.com)
% bunch of 4-shank modifications: Flora

% GUI for planning Neuropixels trajectories with the Allen CCF
% Part of repository: https://github.com/cortex-lab/allenCCF
% - directions for installing atlas in repository readme
% - some dependent functions from that repository
%
% (optional inputs - if CCF path written in Line 22, loaded automatically)
% tv, av, st = CCF template volume, annotated volume, structure tree

% Initialize gui_data structure
gui_data = struct;

% Allen CCF-bregma transform (estimated from eyeballing Paxinos->CCF)
% [AP,DV,ML]
bregma = [540,0,570];

% If not already loaded in, load in atlas
if nargin < 3
    allen_atlas_path = 'C:\Users\Andrew\OneDrive for Business\Documents\Atlases\AllenCCF';
    if isempty(allen_atlas_path)
        error('Enter path where Allen CCF is stored at Line 23');
    end
    tv = readNPY([allen_atlas_path filesep 'template_volume_10um.npy']); % grey-scale "background signal intensity"
    av = readNPY([allen_atlas_path filesep 'annotation_volume_10um_by_index.npy']); % the number at each pixel labels the area, see note below
    st = loadStructureTree([allen_atlas_path filesep 'structure_tree_safe_2017.csv']); % a table of what all the labels mean
end

% Load the colormap (located in the repository, find by associated fcn)
allenCCF_path = fileparts(which('allenCCFbregma'));
cmap_filename = [allenCCF_path filesep 'allen_ccf_colormap_2017.mat'];
load(cmap_filename);

% Set up the gui
probe_atlas_gui = figure('Toolbar','none','Menubar','none','color','w', ...
    'Name','Atlas-probe viewer','Units','normalized','Position',[0.2,0.2,0.7,0.7]);

% Set up the atlas axes
axes_atlas = subplot(1,5,1);
[~, brain_outline] = plotBrainGrid([],axes_atlas);
hold(axes_atlas,'on');
axis vis3d equal off manual
view([-30,25]);
caxis([0 300]);
[ap_max,dv_max,ml_max] = size(tv);
xlim([-10,ap_max+10])
ylim([-10,ml_max+10])
zlim([-10,dv_max+10])

% 4 shank parameters
shank_spacing=0:20:60;
shank_no=4;

for i=1:shank_no
% Set up the probe area axes                    % add 4 of these
    axes_probe_areas(i) = subplot(1,5,i+1);
    axes_probe_areas(i).ActivePositionProperty = 'position';
    set(axes_probe_areas(i),'FontSize',11);
    yyaxis(axes_probe_areas(i),'left');
    probe_areas_plot(i) = image(0);
    set(axes_probe_areas(i),'XTick','','YLim',[0,2880],'YColor','k','YDir','reverse');
    ylabel(axes_probe_areas(i),'');

    colormap(axes_probe_areas(i),cmap);
    caxis([1,size(cmap,1)])
    yyaxis(axes_probe_areas(i),'right');
    set(axes_probe_areas(i),'XTick','','YLim',[0,2880],'YColor','k','YDir','reverse');
    set(axes_probe_areas(i),'YTicklabels',[])
    title(axes_probe_areas(i),sprintf('shank %d',i));
end
yyaxis(axes_probe_areas(1),'left');
ylabel(axes_probe_areas(1),'Depth (\mum)');

% Position the axes              % change positioning so that there are 4
set(axes_atlas,'Position',[-0.15,-0.1,.8,1.2]);
set(axes_probe_areas(1),'Position',[.55,0.05,0.01,0.8]);
set(axes_probe_areas(2),'Position',[.67,0.05,0.01,0.8]);
set(axes_probe_areas(3),'Position',[.79,0.05,0.01,0.8]);
set(axes_probe_areas(4),'Position',[.91,0.05,0.01,0.8]);

% Set the current axes to the atlas (dirty, but some gca requirements)
axes(axes_atlas);

% Set up the probe reference/actual % again set up the reference so that 4
% shans are plotted
probe_ref_top = [bregma(1),bregma(3),0];
probe_ref_bottom = [bregma(1),bregma(3),size(tv,2)];
probe_ref_vector = [probe_ref_top',probe_ref_bottom'];
shank_spacing_vector=repmat(shank_spacing,2,1)';
allshanks_ref_vector=repmat(probe_ref_vector(2,:),shank_no,1)+shank_spacing_vector;

probe_ref_line = line(probe_ref_vector(1,:),allshanks_ref_vector,probe_ref_vector(3,:), ...
    'linewidth',1.5,'color','r','linestyle','--');

% probe_ref_line = line(probe_ref_vector(1,:),probe_ref_vector(2,:),probe_ref_vector(3,:), ...
%     'linewidth',1.5,'color','r','linestyle','--');

probe_length = 288.8; % IMEC phase 2.0 (in 10 ums)
probe_vector = [probe_ref_vector(:,1),diff(probe_ref_vector,[],2)./ ...
    norm(diff(probe_ref_vector,[],2))*probe_length + probe_ref_vector(:,1)];

allshanks_vector=repmat(probe_vector(2,:),shank_no,1)+shank_spacing_vector;
probe_line = line(probe_vector(1,:),allshanks_vector,probe_vector(3,:),'linewidth',3,'color','b','linestyle','-');
set(probe_line(1),'color','g'); % set shank 0 to be a different color

%probe_line = line(probe_vector(1,:),probe_vector(2,:),probe_vector(3,:),'linewidth',3,'color','b','linestyle','-');

% Set up the text to display coordinates
probe_coordinates_text = uicontrol('Style','text','String','', ...
    'Units','normalized','Position',[0,0.95,1,0.05], ...
    'BackgroundColor','w','HorizontalAlignment','left','FontSize',12);

% Store data
gui_data.tv = tv; % Intensity atlas
gui_data.av = av; % Annotated atlas
gui_data.st = st; % Labels table
gui_data.cmap = cmap; % Atlas colormap
gui_data.bregma = bregma; % Bregma for external referencing
gui_data.probe_length = probe_length; % Length of probe
gui_data.structure_plot_idx = []; % Plotted structures
gui_data.probe_angle = [0;0]; % Probe angles in azimuth/elevation
gui_data.shank_no=shank_no;
gui_data.spin=0;
gui_data.shank_spacing=shank_spacing;

%Store handles
gui_data.handles.cortex_outline = brain_outline;
gui_data.handles.structure_patch = []; % Plotted structures
gui_data.handles.axes_atlas = axes_atlas; % Axes with 3D atlas
gui_data.handles.axes_probe_areas = axes_probe_areas; % Axes with probe areas
gui_data.handles.slice_plot = surface('EdgeColor','none'); % Slice on 3D atlas
gui_data.handles.slice_volume = 'tv'; % The volume shown in the slice
gui_data.handles.probe_ref_line = probe_ref_line; % Probe reference line on 3D atlas
gui_data.handles.probe_line = probe_line; % Probe reference line on 3D atlas
gui_data.handles.probe_areas_plot = probe_areas_plot; % Color-coded probe regions
gui_data.probe_coordinates_text = probe_coordinates_text; % Probe coordinates text

% Make 3D rotation the default state (toggle on/off with 'r')
h = rotate3d(axes_atlas);
h.Enable = 'on';
% Update the slice whenever a rotation is completed
h.ActionPostCallback = @update_slice;

% Set functions for key presses
hManager = uigetmodemanager(probe_atlas_gui);
[hManager.WindowListenerHandles.Enabled] = deal(false);
set(probe_atlas_gui,'KeyPressFcn',@key_press);
set(probe_atlas_gui,'KeyReleaseFcn',@key_release);

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

% Display the first slice and update the probe position
update_slice(probe_atlas_gui);
update_probe_coordinates(probe_atlas_gui);

% Display controls
display_controls;

end

function key_press(probe_atlas_gui,eventdata)

% Get guidata
gui_data = guidata(probe_atlas_gui);

switch eventdata.Key

    case 'uparrow'
        if isempty(eventdata.Modifier)
            % Up: move probe anterior
            ap_offset = -10;

            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'XData',get(gui_data.handles.probe_ref_line(myshank,1),'XData') + ap_offset);
            set(gui_data.handles.probe_line(myshank,1),'XData',get(gui_data.handles.probe_line(myshank,1),'XData') + ap_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-up: increase DV angle
            angle_change=-1;
            update_rotation(probe_atlas_gui,angle_change,'Elevation');
            % update gui value
            new_angle = gui_data.probe_angle(2) + angle_change;
            gui_data.probe_angle(2) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % Alt-up: raise probe
            probe_offset = -10;

            for myshank=1:numel(gui_data.handles.probe_line)
            old_probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

            move_probe_vector = diff(old_probe_vector,[],2)./ ...
                norm(diff(old_probe_vector,[],2))*probe_offset;

            new_probe_vector = bsxfun(@plus,old_probe_vector,move_probe_vector);

            set(gui_data.handles.probe_line(myshank,1),'XData',new_probe_vector(1,:), ...
                'YData',new_probe_vector(2,:),'ZData',new_probe_vector(3,:));
            end
        end

    case 'downarrow'
        if isempty(eventdata.Modifier)
            % Down: move probe posterior
            ap_offset = 10;
            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'XData',get(gui_data.handles.probe_ref_line(myshank,1),'XData') + ap_offset);
            set(gui_data.handles.probe_line(myshank,1),'XData',get(gui_data.handles.probe_line(myshank,1),'XData') + ap_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-down: decrease DV angle

            angle_change=1;
            update_rotation(probe_atlas_gui,angle_change,'Elevation');
            new_angle = gui_data.probe_angle(2) + angle_change;
            gui_data.probe_angle(2) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % Alt-down: lower probe
            probe_offset = 10;
            for myshank=1:numel(gui_data.handles.probe_line)  % this could be yoked together
            old_probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

            move_probe_vector = diff(old_probe_vector,[],2)./ ...
                norm(diff(old_probe_vector,[],2))*probe_offset;

            new_probe_vector = bsxfun(@plus,old_probe_vector,move_probe_vector);

            set(gui_data.handles.probe_line(myshank,1),'XData',new_probe_vector(1,:), ...
                'YData',new_probe_vector(2,:),'ZData',new_probe_vector(3,:));
            end
        end

    case 'rightarrow'
        if isempty(eventdata.Modifier)
            % Right: move probe right
            ml_offset = 10;
            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'YData',get(gui_data.handles.probe_ref_line(myshank,1),'YData') + ml_offset);
            set(gui_data.handles.probe_line(myshank,1),'YData',get(gui_data.handles.probe_line(myshank,1),'YData') + ml_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-right: increase vertical angle
%             angle_change = [1;0];
%             new_angle = gui_data.probe_angle + angle_change;
%             gui_data.probe_angle = new_angle;
%             update_probe_angle(probe_atlas_gui);
            angle_change=-1;
            update_rotation(probe_atlas_gui,angle_change,'Azimuth');
            new_angle = gui_data.probe_angle(1) + angle_change;
            gui_data.probe_angle(1) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % alt-right : rotate clockwise
            angle_change=-1;
            update_rotation(probe_atlas_gui,angle_change,'spin');
            new_angle=gui_data.spin+angle_change;
            gui_data.spin=new_angle;
        end

    case 'leftarrow'
        if isempty(eventdata.Modifier)
            % Left: move probe left
            ml_offset = -10;
            for myshank=1:numel(gui_data.handles.probe_ref_line)
            set(gui_data.handles.probe_ref_line(myshank,1),'YData',get(gui_data.handles.probe_ref_line(myshank,1),'YData') + ml_offset);
            set(gui_data.handles.probe_line(myshank,1),'YData',get(gui_data.handles.probe_line(myshank,1),'YData') + ml_offset);
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-right: increase vertical angle
%             angle_change = [-1;0];
%             new_angle = gui_data.probe_angle + angle_change;
%             gui_data.probe_angle = new_angle;
%             update_probe_angle(probe_atlas_gui);
            angle_change=1;
            update_rotation(probe_atlas_gui,angle_change,'Azimuth');
            new_angle = gui_data.probe_angle(1) + angle_change;
            gui_data.probe_angle(1) = new_angle;

        elseif any(strcmp(eventdata.Modifier,'alt'))
            % alt-left : rotate anticlockwise
            angle_change=1;
            update_rotation(probe_atlas_gui,angle_change,'spin');
            new_angle=gui_data.spin+angle_change;
            gui_data.spin=new_angle;
        end

    case 'c'
        % Bring up controls again
        display_controls;

    case 'b'
        % Toggle brain outline visibility
        current_visibility = gui_data.handles.cortex_outline.Visible;
        switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
        set(gui_data.handles.cortex_outline,'Visible',new_visibility);

    case 'a'
        % Toggle plotted structure visibility
        if ~isempty(gui_data.structure_plot_idx)
            current_visibility = get(gui_data.handles.structure_patch(1),'Visible');
            switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
            set(gui_data.handles.structure_patch,'Visible',new_visibility);
        end

    case 's'
        % Toggle slice volume/visibility
        slice_volumes = {'tv','av','none'};
        new_slice_volume = slice_volumes{circshift( ...
            strcmp(gui_data.handles.slice_volume,slice_volumes),[0,1])};

        if strcmp(new_slice_volume,'none')
            set(gui_data.handles.slice_plot,'Visible','off');
        else
            set(gui_data.handles.slice_plot,'Visible','on');
        end

        gui_data.handles.slice_volume = new_slice_volume;
        guidata(probe_atlas_gui, gui_data);

        update_slice(probe_atlas_gui);

    case 'p'
        % Toggle probe visibility
        current_visibility = gui_data.handles.probe_ref_line.Visible;
        switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
        set(gui_data.handles.probe_ref_line,'Visible',new_visibility);
        set(gui_data.handles.probe_line,'Visible',new_visibility);

    case 'r'
        % Toggle 3D rotation
        h = rotate3d(gui_data.handles.axes_atlas);
        switch h.Enable
            case 'off'
                h.Enable = 'on';
                % Update the slice whenever a rotation is completed
                h.ActionPostCallback = @update_slice;
                %(need to restore key-press functionality with rotation)
                hManager = uigetmodemanager(probe_atlas_gui);
                [hManager.WindowListenerHandles.Enabled] = deal(false);
                set(probe_atlas_gui,'KeyPressFcn',@key_press);
            case 'on'
                h.Enable = 'off';
        end

    case 'm'
        % Set probe angle
        set_probe_position(probe_atlas_gui);
        % Get updated guidata
        gui_data = guidata(probe_atlas_gui);

    case {'equal','add'}
        % Add structure(s) to display
        slice_spacing = 10;

        % Prompt for which structures to show (only structures which are
        % labelled in the slice-spacing downsampled annotated volume)

        if ~any(strcmp(eventdata.Modifier,'shift'))
            % (no shift: list in native CCF order)
            parsed_structures = unique(reshape(gui_data.av(1:slice_spacing:end, ...
                1:slice_spacing:end,1:slice_spacing:end),[],1));
            plot_structures_parsed = listdlg('PromptString','Select a structure to plot:', ...
                'ListString',gui_data.st.safe_name(parsed_structures),'ListSize',[520,500]);
            plot_structures = parsed_structures(plot_structures_parsed);

            if ~isempty(plot_structures)
                for curr_plot_structure = reshape(plot_structures,1,[])
                    % If this label isn't used, don't plot
                    if ~any(reshape(gui_data.av( ...
                            1:slice_spacing:end,1:slice_spacing:end,1:slice_spacing:end),[],1) == curr_plot_structure)
                        disp(['"' gui_data.st.safe_name{curr_plot_structure} '" is not parsed in the atlas'])
                        continue
                    end

                    gui_data.structure_plot_idx(end+1) = curr_plot_structure;

                    plot_structure_color = hex2dec(reshape(gui_data.st.color_hex_triplet{curr_plot_structure},2,[])')./255;
                    structure_3d = isosurface(permute(gui_data.av(1:slice_spacing:end, ...
                        1:slice_spacing:end,1:slice_spacing:end) == curr_plot_structure,[3,1,2]),0);

                    structure_alpha = 0.2;
                    gui_data.handles.structure_patch(end+1) = patch('Vertices',structure_3d.vertices*slice_spacing, ...
                        'Faces',structure_3d.faces, ...
                        'FaceColor',plot_structure_color,'EdgeColor','none','FaceAlpha',structure_alpha);
                end
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % (shift: use hierarchy search)
            plot_structures = hierarchicalSelect(gui_data.st);

            if ~isempty(plot_structures) % will be empty if dialog was cancelled
                % get all children of this one
                thisID = gui_data.st.id(plot_structures);
                idStr = sprintf('/%d/', thisID);
                theseCh = find(cellfun(@(x)contains(x,idStr), gui_data.st.structure_id_path));

                % plot the structure
                slice_spacing = 5;
                plot_structure_color = hex2dec(reshape(gui_data.st.color_hex_triplet{plot_structures},3,[]))./255;
                structure_3d = isosurface(permute(ismember(gui_data.av(1:slice_spacing:end, ...
                    1:slice_spacing:end,1:slice_spacing:end),theseCh),[3,1,2]),0);

                structure_alpha = 0.2;
                gui_data.structure_plot_idx(end+1) = plot_structures;
                gui_data.handles.structure_patch(end+1) = patch('Vertices',structure_3d.vertices*slice_spacing, ...
                    'Faces',structure_3d.faces, ...
                    'FaceColor',plot_structure_color,'EdgeColor','none','FaceAlpha',structure_alpha);

            end

        end

    case {'hyphen','subtract'}
        % Remove structure(s) already plotted
        if ~isempty(gui_data.structure_plot_idx)
            remove_structures = listdlg('PromptString','Select a structure to remove:', ...
                'ListString',gui_data.st.safe_name(gui_data.structure_plot_idx));
            delete(gui_data.handles.structure_patch(remove_structures))
            gui_data.structure_plot_idx(remove_structures) = [];
            gui_data.handles.structure_patch(remove_structures) = [];
        end

    case 'x'
        % Export the probe coordinates in Allen CCF to the workspace
        probe_vector = cell2mat(get(gui_data.handles.probe_line,{'XData','YData','ZData'})');
        probe_vector_ccf = round(probe_vector([1,3,2],:))';
        assignin('base','probe_vector_ccf',probe_vector_ccf)
        disp('Copied probe vector coordinates to workspace');

    case 'h'
        % Load probe histology points, plot line of best fit
        probe_histology(probe_atlas_gui);
end

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function key_release(probe_atlas_gui,eventdata)

% Get guidata
gui_data = guidata(probe_atlas_gui);

switch eventdata.Key
    case {'rightarrow','leftarrow','uparrow','downarrow'}
        % Update the probe info/slice on arrow release
        update_probe_coordinates(probe_atlas_gui);
        update_slice(probe_atlas_gui);
end

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function update_slice(probe_atlas_gui,varargin)

% Get guidata
gui_data = guidata(probe_atlas_gui);

% Only update the slice if it's visible
if strcmp(gui_data.handles.slice_plot(1).Visible,'on')

    % Get current position of camera
    curr_campos = campos;

    % Get probe vector
    probe_ref_top = [gui_data.handles.probe_ref_line(1,1).XData(1), ...
        gui_data.handles.probe_ref_line(1,1).YData(1),gui_data.handles.probe_ref_line(1,1).ZData(1)];
    probe_ref_bottom = [gui_data.handles.probe_ref_line(1).XData(2), ...
        gui_data.handles.probe_ref_line(1,1).YData(2),gui_data.handles.probe_ref_line(1,1).ZData(2)];
    probe_vector = probe_ref_top - probe_ref_bottom;

    % Get probe-camera vector
    probe_camera_vector = probe_ref_top - curr_campos;

    % Get the vector to plot the plane in (along with probe vector)
    plot_vector = cross(probe_camera_vector,probe_vector);

    % Get the normal vector of the plane
    normal_vector = cross(plot_vector,probe_vector);

    % Get the plane offset through the probe
    plane_offset = -(normal_vector*probe_ref_top');

    % Define a plane of points to index
    % (the plane grid is defined based on the which cardinal plan is most
    % orthogonal to the plotted plane. this is janky but it works)
    slice_px_space = 3;
    %[~,cam_plane] = max(abs((campos - camtarget)./norm(campos - camtarget)));

    [~,cam_plane] = max(abs(normal_vector./norm(normal_vector)));

    switch cam_plane

        case 1
            [plane_y,plane_z] = meshgrid(1:slice_px_space:size(gui_data.tv,3),1:slice_px_space:size(gui_data.tv,2));
            plane_x = ...
                (normal_vector(2)*plane_y+normal_vector(3)*plane_z + plane_offset)/ ...
                -normal_vector(1);

        case 2
            [plane_x,plane_z] = meshgrid(1:slice_px_space:size(gui_data.tv,1),1:slice_px_space:size(gui_data.tv,2));
            plane_y = ...
                (normal_vector(1)*plane_x+normal_vector(3)*plane_z + plane_offset)/ ...
                -normal_vector(2);

        case 3
            [plane_x,plane_y] = meshgrid(1:slice_px_space:size(gui_data.tv,1),1:slice_px_space:size(gui_data.tv,3));
            plane_z = ...
                (normal_vector(1)*plane_x+normal_vector(2)*plane_y + plane_offset)/ ...
                -normal_vector(3);

    end

    % Get the coordiates on the plane
    x_idx = round(plane_x);
    y_idx = round(plane_y);
    z_idx = round(plane_z);

    % Find plane coordinates in bounds with the volume
    use_xd = x_idx > 0 & x_idx < size(gui_data.tv,1);
    use_yd = y_idx > 0 & y_idx < size(gui_data.tv,3);
    use_zd = z_idx > 0 & z_idx < size(gui_data.tv,2);
    use_idx = use_xd & use_yd & use_zd;

    curr_slice_idx = sub2ind(size(gui_data.tv),x_idx(use_idx),z_idx(use_idx),y_idx(use_idx));

    % Find plane coordinates that contain brain
    curr_slice_isbrain = false(size(use_idx));
    curr_slice_isbrain(use_idx) = gui_data.av(curr_slice_idx) > 1;

    % Index coordinates in bounds + with brain
    grab_pix_idx = sub2ind(size(gui_data.tv),x_idx(curr_slice_isbrain),z_idx(curr_slice_isbrain),y_idx(curr_slice_isbrain));

    % Grab pixels from (selected) volume
    curr_slice = nan(size(use_idx));
    switch gui_data.handles.slice_volume
        case 'tv'
            curr_slice(curr_slice_isbrain) = gui_data.tv(grab_pix_idx);
            colormap(gui_data.handles.axes_atlas,'gray');
            caxis([0,255]);
        case 'av'
            curr_slice(curr_slice_isbrain) = gui_data.av(grab_pix_idx);
            colormap(gui_data.handles.axes_atlas,gui_data.cmap);
            caxis([1,size(gui_data.cmap,1)]);
    end

    % Update the slice display
    set(gui_data.handles.slice_plot,'XData',plane_x,'YData',plane_y,'ZData',plane_z,'CData',curr_slice);

    % Upload gui_data
    guidata(probe_atlas_gui, gui_data);

end

end

function set_probe_position(probe_atlas_gui,varargin)

% Get guidata
gui_data = guidata(probe_atlas_gui);

% Prompt for angles
prompt_text = { ...
    'AP position (mm from bregma)', ...
    'ML position (mm from bregma)', ...
    'ML angle (relative to midline)', ....
    'DV angle (relative to coronal cut)',...
    'probe spin'};

new_probe_position = cellfun(@str2num,inputdlg(prompt_text,'Set probe position',1));

% Convert probe position: mm->CCF and degrees->radians
probe_ccf_coordinates = round(gui_data.bregma([1,3])' - new_probe_position(1:2)*100);
%probe_angle_rad = (new_probe_position(3:4)/360)*2*pi;

% Update the probe and trajectory reference
[ap_max,dv_max,ml_max] = size(gui_data.tv);

max_ref_length = sqrt(sum(([ap_max,dv_max,ml_max].^2)));

%[~,y,z] = sph2cart(probe_angle_rad(1),probe_angle_rad(2),max_ref_length);

%[~,x,~] = sph2cart(probe_angle_rad(2),probe_angle_rad(1),max_ref_length);

probe_ref_top = [probe_ccf_coordinates(1),probe_ccf_coordinates(2),0];
probe_ref_bottom = [probe_ccf_coordinates(1),probe_ccf_coordinates(2),max_ref_length];
probe_ref_vector = [probe_ref_top;probe_ref_bottom]';
probe_vector = [probe_ref_vector(:,1),diff(probe_ref_vector,[],2)./ ...
    norm(diff(probe_ref_vector,[],2))*gui_data.probe_length + probe_ref_vector(:,1)];

shank_spacing_vector=reshape([gui_data.shank_spacing;gui_data.shank_spacing],[],1)';

probe_ref_vector=repmat(probe_ref_vector,1,gui_data.shank_no);
probe_ref_vector(2,:)=probe_ref_vector(2,:)+shank_spacing_vector;

probe_vector=repmat(probe_vector,1,gui_data.shank_no);
probe_vector(2,:)=probe_vector(2,:)+shank_spacing_vector;

gui_data.probe_angle(1) =new_probe_position(3);
gui_data.probe_angle(2) =new_probe_position(4);
gui_data.spin=new_probe_position(5);
% % spin
% rotation_matrix=get_rotation_matrix(gui_data.spin,'spin');
% [probe_ref_vector,probe_vector]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix);
%
% % rotate ML
% rotation_matrix=get_rotation_matrix(gui_data.probe_angle(1),'ML');
% [probe_ref_vector,probe_vector]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix);
%
% % rotate AP
% rotation_matrix=get_rotation_matrix(gui_data.probe_angle(2),'AP');
% [probe_ref_vector,probe_vector]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix);

for myshank=1:gui_data.shank_no
    st_ix=2*myshank-1;
    set(gui_data.handles.probe_ref_line(myshank),'XData',probe_ref_vector(1,st_ix:(st_ix+1)), ...
            'YData',probe_ref_vector(2,st_ix:(st_ix+1)), ...
            'ZData',probe_ref_vector(3,st_ix:(st_ix+1)));

    set(gui_data.handles.probe_line(myshank),'XData',probe_vector(1,st_ix:(st_ix+1)), ...
        'YData',probe_vector(2,st_ix:(st_ix+1)), ...
        'ZData',probe_vector(3,st_ix:(st_ix+1)));
end

% Upload gui_data
guidata(probe_atlas_gui, gui_data);

% Update the slice and probe coordinates
update_slice(probe_atlas_gui);
update_probe_coordinates(probe_atlas_gui);

end

function update_probe_angle(probe_atlas_gui,varargin)
% not used right now
% Get guidata
gui_data = guidata(probe_atlas_gui);

% Get the positions of the probe and trajectory reference
for myshank=1:numel(gui_data.handles.probe_ref_line)
    probe_ref_vector = cell2mat(get(gui_data.handles.probe_ref_line(myshank,1),{'XData','YData','ZData'})');
    probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

    % Update the probe trajectory reference (rotate about trajectory top)
    [ap_max,dv_max,ml_max] = size(gui_data.tv);

    max_ref_length = sqrt(sum(([ap_max,dv_max,ml_max].^2)));

    probe_angle_rad = (gui_data.probe_angle./360)*2*pi;
    [x,y,z] = sph2cart(pi-probe_angle_rad(1),probe_angle_rad(2),max_ref_length);

    new_probe_ref_top = [probe_ref_vector(1,1),probe_ref_vector(2,1),0];
    new_probe_ref_bottom = new_probe_ref_top + [x,y,z];
    new_probe_ref_vector = [new_probe_ref_top;new_probe_ref_bottom]';

    set(gui_data.handles.probe_ref_line(myshank,1),'XData',new_probe_ref_vector(1,:), ...
        'YData',new_probe_ref_vector(2,:), ...
        'ZData',new_probe_ref_vector(3,:));

    % Update probe (retain depth)
    new_probe_vector = [new_probe_ref_vector(:,1),diff(new_probe_ref_vector,[],2)./ ...
        norm(diff(new_probe_ref_vector,[],2))*gui_data.probe_length + new_probe_ref_vector(:,1)];

    probe_depth = sqrt(sum((probe_ref_vector(:,1) - probe_vector(:,1)).^2));
    new_probe_vector_depth = (diff(new_probe_vector,[],2)./ ...
        norm(diff(new_probe_vector,[],2))*probe_depth) + new_probe_vector;

    set(gui_data.handles.probe_line(myshank,1),'XData',new_probe_vector_depth(1,:), ...
        'YData',new_probe_vector_depth(2,:),'ZData',new_probe_vector_depth(3,:));
end
% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function update_probe_coordinates(probe_atlas_gui,varargin)

% Get guidata
gui_data = guidata(probe_atlas_gui);

for myshank=1:numel(gui_data.handles.probe_ref_line)
    % Get the positions of the probe and trajectory reference
    probe_ref_vector = cell2mat(get(gui_data.handles.probe_ref_line(myshank,1),{'XData','YData','ZData'})');
    probe_vector = cell2mat(get(gui_data.handles.probe_line(myshank,1),{'XData','YData','ZData'})');

    trajectory_n_coords = max(abs(diff(probe_ref_vector,[],2)));
    [trajectory_xcoords,trajectory_ycoords,trajectory_zcoords] = deal( ...
        linspace(probe_ref_vector(1,1),probe_ref_vector(1,2),trajectory_n_coords), ...
        linspace(probe_ref_vector(2,1),probe_ref_vector(2,2),trajectory_n_coords), ...
        linspace(probe_ref_vector(3,1),probe_ref_vector(3,2),trajectory_n_coords));

    probe_n_coords = sqrt(sum(diff(probe_vector,[],2).^2));
    [probe_xcoords,probe_ycoords,probe_zcoords] = deal( ...
        linspace(probe_vector(1,1),probe_vector(1,2),probe_n_coords), ...
        linspace(probe_vector(2,1),probe_vector(2,2),probe_n_coords), ...
        linspace(probe_vector(3,1),probe_vector(3,2),probe_n_coords));

    % Get brain labels across the probe and trajectory, and intersection with brain
    pixel_space = 5;
    trajectory_areas = interp3(single(gui_data.av(1:pixel_space:end,1:pixel_space:end,1:pixel_space:end)), ...
        round(trajectory_zcoords/pixel_space),round(trajectory_xcoords/pixel_space),round(trajectory_ycoords/pixel_space),'nearest');
    trajectory_brain_idx = find(trajectory_areas > 1,1);
    trajectory_brain_intersect = ...
        [trajectory_xcoords(trajectory_brain_idx),trajectory_ycoords(trajectory_brain_idx),trajectory_zcoords(trajectory_brain_idx)]';

    probe_areas = interp3(single(gui_data.av(1:pixel_space:end,1:pixel_space:end,1:pixel_space:end)), ...
        round(probe_zcoords/pixel_space),round(probe_xcoords/pixel_space),round(probe_ycoords/pixel_space),'nearest')';
    probe_area_boundaries = intersect(unique([find(~isnan(probe_areas),1,'first'); ...
        find(diff(probe_areas) ~= 0);find(~isnan(probe_areas),1,'last')]),find(~isnan(probe_areas)));
    probe_area_centers = probe_area_boundaries(1:end-1) + diff(probe_area_boundaries)/2;
    probe_area_labels = gui_data.st.safe_name(probe_areas(round(probe_area_centers)));

    % Get position of brain intersect relative to bregma
    probe_bregma_coordinate = round((gui_data.bregma([1,3])' - trajectory_brain_intersect(1:2))*10);

    % Get the depth of the bottom of the probe (sign: hack by z offset)
    probe_depth = round(sqrt(sum((trajectory_brain_intersect - probe_vector(:,2)).^2))*10)* ...
        sign(probe_vector(3,2)-trajectory_brain_intersect(3));

    % Update the text
    if myshank==1
        probe_text = ['Shank 1 insertion: ' ....
            num2str(probe_bregma_coordinate(1)) ' AP, ', ...
            num2str(-probe_bregma_coordinate(2)) ' ML, ', ...
            num2str(probe_depth) ' Probe-axis, ' ...
            num2str(gui_data.probe_angle(1)+90) char(176) ' from midline, ' ...
            num2str(-(gui_data.probe_angle(2)+90)) char(176) ' from horizontal,' ...
            num2str(gui_data.spin) char(176) ' spin '];
        set(gui_data.probe_coordinates_text,'String',probe_text);
    end

    % Update the probe areas
    yyaxis(gui_data.handles.axes_probe_areas(myshank),'right');
    set(gui_data.handles.probe_areas_plot(myshank),'YData',[1:length(probe_areas)]*10,'CData',probe_areas);
    set(gui_data.handles.axes_probe_areas(myshank),'YTick',probe_area_centers*10,'YTickLabels',probe_area_labels','YTickLabelRotation',45);

end
% Upload gui_data
guidata(probe_atlas_gui, gui_data);

end

function probe_histology(probe_atlas_gui)
% Load histology points
% UNDER CONSTRUCTION
% (used to use SHARP-Track, now using mine)

% Get guidata
gui_data = guidata(probe_atlas_gui);

[probe_file,probe_path] = uigetfile('*.mat','Choose probe coordinate file');
load([probe_path,probe_file]);

if exist('pointList','var')
    histology_points = pointList.pointList{1};
elseif exist('probe_ccf','var')
    histology_points = probe_ccf(1).points; % only use first probe
end

r0 = mean(histology_points,1);
xyz = bsxfun(@minus,histology_points,r0);
[~,~,V] = svd(xyz,0);
histology_probe_direction = V(:,1);

probe_eval_points = [-1000,1000];
probe_line_endpoints = bsxfun(@plus,bsxfun(@times,probe_eval_points',histology_probe_direction'),r0);

% Philip's GUI: not saved in native CCF order?
% plot3(histology_points(:,3),histology_points(:,1),histology_points(:,2),'.b','MarkerSize',20);
% line(P(:,3),P(:,1),P(:,2),'color','k','linewidth',2)

% % Mine: saved in native CCF order [AP,DV,ML]
plot3(histology_points(:,1),histology_points(:,3),histology_points(:,2),'.b','MarkerSize',20);
line(probe_line_endpoints(:,1),probe_line_endpoints(:,3),probe_line_endpoints(:,2),'color','k','linewidth',2)

% Place the probe on the histology best-fit axis
[ap_max,dv_max,ml_max] = size(gui_data.tv);

probe_ref_top = probe_line_endpoints(1,[1,3,2]);
probe_ref_bottom = probe_line_endpoints(2,[1,3,2]);
probe_ref_vector = [probe_ref_top;probe_ref_bottom]';

set(gui_data.handles.probe_ref_line,'XData',probe_ref_vector(1,:), ...
    'YData',probe_ref_vector(2,:), ...
    'ZData',probe_ref_vector(3,:));

probe_vector = [probe_ref_vector(:,1),diff(probe_ref_vector,[],2)./ ...
    norm(diff(probe_ref_vector,[],2))*gui_data.probe_length + probe_ref_vector(:,1)];
set(gui_data.handles.probe_line,'XData',probe_vector(1,:), ...
    'YData',probe_vector(2,:),'ZData',probe_vector(3,:));

% Upload gui_data
[theta,phi] = cart2sph(diff(probe_ref_vector(1,:)),diff(probe_ref_vector(2,:)),diff(probe_ref_vector(3,:)));
gui_data.probe_angle = ([theta,phi]/(2*pi))*360;
guidata(probe_atlas_gui, gui_data);

% Update the slice and probe coordinates
update_slice(probe_atlas_gui);
update_probe_coordinates(probe_atlas_gui);
end

function display_controls

% Print controls
CreateStruct.Interpreter = 'tex';
CreateStruct.WindowStyle = 'non-modal';
msgbox( ...
    {'\fontsize{12}' ...
    '\bf Probe: \rm' ...
    'Arrow keys : translate probe' ...
    'Alt/Option up/down : raise/lower probe' ...
    'Alt/Option left/right : spin probe' ...
    'Shift arrow keys : rotate probe (around top)' ...
    'm : set probe location manually', ...
    '\bf 3D brain areas: \rm' ...
    ' =/+ : add (list selector)' ...
    ' Shift =/+ : add (hierarchy selector)' ...
    ' - : remove', ...
    '\bf Visibility: \rm' ...
    's : atlas slice (toggle tv/av/off)' ...
    'b : brain outline' ...
    'p : probe' ...
    'a : 3D brain areas' ...
    '\bf Other: \rm' ...
    'r : toggle clickable rotation' ...
    'x : export probe coordinates to workspace' ...
    'h : load and plot histology-defined trajectory', ...
    'c : bring up controls box'}, ...
    'Controls',CreateStruct);

end

function [rotation_matrix]=get_rotation_matrix(probe_vector,change_angle,current_angle,direction)
% Zhiwen modified
% Convert Axis-Angle Rotation to Rotation Matrix
angle_rad=deg2rad(change_angle);
current_angle_rad = deg2rad(current_angle);
probeV = probe_vector(:,2)-probe_vector(:,1);
probehV = probe_vector(:,3)-probe_vector(:,1);
    switch direction
        case 'Azimuth'
            % Azimuth
            % Convert Axis-Angle Rotation to Rotation Matrix
            % DV line as rotation axis
            axang = [0 0 288.8 angle_rad];
        case 'Elevation'
            % Elevation
            % Convert Axis-Angle Rotation to Rotation Matrix
            % plane of 4-shank probe is the plane to rotate on;
            % which means it rotates around vector orthogonal to this
            % plane;
            axang1 = [0 0 1 current_angle_rad];
            rotation_matrix1 = axang2rotm1(axang1);
            rotated_azimuth = rotation_matrix1*[1; 0; 0];
            %%
            othorg = cross(probeV,rotated_azimuth);
            axang = [othorg(1) othorg(2) othorg(3) angle_rad];
        case 'spin'
            % spin
            % Convert Axis-Angle Rotation to Rotation Matrix
            % spin around reference vector;
            axang = [probeV(1,1) probeV(2,1) probeV(3,1) angle_rad];
    end
    rotation_matrix = axang2rotm1(axang);
end

function [rotated_ref_bregma,rotated_probe_bregma]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix)
    zeropoint=probe_vector(:,1);
    zerocentered_ref=probe_ref_vector-zeropoint;

    zerocentered_probe=probe_vector-zeropoint;

    rotated_ref=rotation_matrix*zerocentered_ref;
    rotated_ref_bregma=rotated_ref+zeropoint;

    rotated_probe=rotation_matrix*zerocentered_probe;
    rotated_probe_bregma=rotated_probe+zeropoint;
end

function update_rotation(probe_atlas_gui,angle_change,mydir)
% this is the function when we update the rotation of the probe either in
% ML or AP
    gui_data = guidata(probe_atlas_gui);
    probe_ref_vector = cell2mat(get(gui_data.handles.probe_ref_line,{'XData','YData','ZData'})');
    probe_vector = cell2mat(get(gui_data.handles.probe_line,{'XData','YData','ZData'})');
    current_angle = gui_data.probe_angle(1)+90;

    rotation_matrix=get_rotation_matrix(probe_vector,angle_change,current_angle,mydir);
    [rotated_ref_bregma,rotated_probe_bregma]=please_rotate(probe_ref_vector,probe_vector,rotation_matrix);

    for myshank=1:gui_data.shank_no
        st_ix=2*myshank-1;
        set(gui_data.handles.probe_ref_line(myshank),'XData',rotated_ref_bregma(1,st_ix:(st_ix+1)), ...
                'YData',rotated_ref_bregma(2,st_ix:(st_ix+1)), ...
                'ZData',rotated_ref_bregma(3,st_ix:(st_ix+1)));

        set(gui_data.handles.probe_line(myshank),'XData',rotated_probe_bregma(1,st_ix:(st_ix+1)), ...
            'YData',rotated_probe_bregma(2,st_ix:(st_ix+1)), ...
            'ZData',rotated_probe_bregma(3,st_ix:(st_ix+1)));
    end
end

function R = axang2rotm1( axang )
%AXANG2ROTM Convert axis-angle rotation representation to rotation matrix
%   R = AXANG2ROTM(AXANG) converts a 3D rotation given in axis-angle form,
%   AXANG, to an orthonormal rotation matrix, R. AXANG is an N-by-4
%   matrix of N axis-angle rotations. The first three elements of every
%   row specify the rotation axis and the last element defines the rotation
%   angle (in radians).
%   The output, R, is an 3-by-3-by-N matrix containing N rotation matrices.
%   Each rotation matrix has a size of 3-by-3 and is orthonormal.
%
%   Example:
%      % Convert a rotation from axis-angle to rotation matrix
%      axang = [0 1 0 pi/2];
%      R = axang2rotm(axang)

% For a single axis-angle vector [ax ay az theta] the output rotation
% matrix R can be computed as follows:
% R =  [t*x*x + c	  t*x*y - z*s	   t*x*z + y*s
%       t*x*y + z*s	  t*y*y + c	       t*y*z - x*s
%       t*x*z - y*s	  t*y*z + x*s	   t*z*z + c]
% where,
% c = cos(theta)
% s = sin(theta)
% t = 1 - c
% x = normalized axis ax coordinate
% y = normalized axis ay coordinate
% z = normalized axis az coordinate

% Normalize the axis
v = bsxfun(@times, axang(:,1:3), 1./sqrt(sum(axang(:,1:3).^2,2)));
% Extract the rotation angles and shape them in depth dimension
numInputs = size(axang,1);
theta = zeros(1,1,numInputs);
theta(1,1,:) = axang(:,4);

% Compute rotation matrices
cth = cos(theta);
sth = sin(theta);
vth = (1 - cth);

% Preallocate input vectors
vx = zeros(1,1,numInputs,'like',axang);
vy = vx;
vz = vx;

% Shape input vectors in depth dimension
vx(1,1,:) = v(:,1);
vy(1,1,:) = v(:,2);
vz(1,1,:) = v(:,3);

% Explicitly specify concatenation dimension
tempR = cat(1, vx.*vx.*vth+cth,     vy.*vx.*vth-vz.*sth, vz.*vx.*vth+vy.*sth, ...
               vx.*vy.*vth+vz.*sth, vy.*vy.*vth+cth,     vz.*vy.*vth-vx.*sth, ...
               vx.*vz.*vth-vy.*sth, vy.*vz.*vth+vx.*sth, vz.*vz.*vth+cth);

R = reshape(tempR, [3, 3, length(vx)]);
R = permute(R, [2 1 3]);

end

author: cortex-lab

repo: allenCCF

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