text:
function gen_cropped_MPII_test_images()
% Author: Guanghan Ning, Feb 2017
clear all;
clc;

target_dist = 1.45;
boxsize = 368;

% Get [input images paths] and [output images path]
mpii_img_folder_path = '../dataset/MPI/images/';
img_output_folder = 'dataset_mpii/images_test_cropped/qualitative/';

% Load annotations
load('eval_MPII/mpii_Tompson_test', 'annolist');
load('eval_MPII/mpii_Tompson_test', 'keypointsAll');
load('eval_MPII/mpii_Tompson_test', 'RELEASE_img_index');
load('eval_MPII/mpii_Tompson_test', 'RELEASE_person_index');

annolist_test = annolist;
num_test_imgs = size(annolist_test, 2);

load('eval_MPII/mpii_human_pose_v1_u12_1','RELEASE');
annolist_all = RELEASE.annolist;

%% Start from here if everything already read from mat file.
for img_id = 1:num_test_imgs

    img_index = RELEASE_img_index(img_id);

    % (1) Deal with rect ids

        % Get rect id. Note the difference between rect_ct and rect_id
        rect_id = RELEASE_person_index(img_id);

        % Find scale and mid-points in order to convert prediction
        x_mid = annolist_test(img_id).annorect.objpos.x;
        y_mid = annolist_test(img_id).annorect.objpos.y;
        objpos_org = [x_mid, y_mid];
        scale_provided = annolist_test(img_id).annorect.scale;

        % (2) Deal with images
        img_name = annolist_all(img_index).image.name;
        img_path = strcat(mpii_img_folder_path, img_name);
        if ~exist(img_path, 'file')
            continue
        end

        img = imread(img_path);

        scale_in_cpp = target_dist/scale_provided;
        center_s = objpos_org * scale_in_cpp;
        fprintf('scale: %f\n', scale_in_cpp);

        img = imresize(img, scale_in_cpp);
        disp(size(img));
        [img, pad] = padAround(img, boxsize, center_s);
        disp(size(img));

        [pathstr, name, ext] = fileparts(img_path);
        img_output_path = strcat(img_output_folder, name, '_', num2str(rect_id), ext);
        imwrite(img, img_output_path);

        output_scale(img_id) = scale_in_cpp;
end
save('scale_in_cpp_MPII_test.mat', 'output_scale');

function [img_padded, pad] = padAround(img, boxsize, center)
    center = round(center);
    h = size(img, 1);
    w = size(img, 2);
    pad(1) = boxsize/2 - center(2); % up
    pad(3) = boxsize/2 - (h-center(2)); % down
    pad(2) = boxsize/2 - center(1); % left
    pad(4) = boxsize/2 - (w-center(1)); % right

    pad_up = repmat(img(1,:,:)*0, [pad(1) 1 1])+128;
    img_padded = [pad_up; img];
    pad_left = repmat(img_padded(:,1,:)*0, [1 pad(2) 1])+128;
    img_padded = [pad_left img_padded];
    pad_down = repmat(img_padded(end,:,:)*0, [pad(3) 1 1])+128;
    img_padded = [img_padded; pad_down];
    pad_right = repmat(img_padded(:,end,:)*0, [1 pad(4) 1])+128;
    img_padded = [img_padded pad_right];

    center = center + [max(0,pad(2)) max(0,pad(1))];
    img_padded = img_padded(center(2)-(boxsize/2-1):center(2)+boxsize/2, center(1)-(boxsize/2-1):center(1)+boxsize/2, :); %cropping if needed

function [x,y] = findMaximum(map)
    [~,i] = max(map(:));
    [x,y] = ind2sub(size(map), i);

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function mask2conv()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    in = mem.mask_li{curr_layer_idx};
    for m = 1:config.batch_size
        for n = 1:config.chs
            mem.mask_inputs{curr_layer_idx}((n-1)*config.kernel_size(1, 1)*config.kernel_size(1, 2)+1:n*config.kernel_size(1, 1)*config.kernel_size(1, 2), (m-1)*(size(mem.mask_inputs{curr_layer_idx}, 2)/config.batch_size)+1:m*size(mem.mask_inputs{curr_layer_idx}, 2)/config.batch_size) = ...
                    config.IM2COL(in(:,:,n,m), [config.kernel_size(1, 1), config.kernel_size(1, 2)]);
        end
    end
    % only for SR now
    gen_mask_patch_cat_idx_for_super_res(mem.mask_inputs{1});
end

author: jimmy-ren


text:
function gX = wienerKernDiagGradX(kern, X, X2)

% WIENERKERNDIAGGRADX Gradient of WIENER kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the wiener kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : wienerKernParamInit, kernDiagGradX, wienerkernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2009

% KERN

gX = repmat(kern.variance, [1 1 size(X, 1)]);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function DCV=plsldadcv(X,y,A,K,method,OPT,order)
%+++ K-fold double cross validation Cross-validation for PLS-LDA
%+++ Input:  X: m x n  (Sample matrix)
%            y: m x 1  (measured property)
%            A: The max PC for cross-validation
%            K: fold. when K = m, it is leave-one-out CV
%       method: pretreatment method. Contains: autoscaling, center etc.
%          OPT: =1 Print process.
%               =0 No print.
%               pareto,minmax,center or none.
%+++ Order: =1  sorted,default. For CV partition.
%           =0  random.
%+++ Output: Structural data: CV
%+++ Hongdong Li, Oct. 16, 2008.
%+++ Revised in Jan.12, 2009.

if nargin<7;order=1;end;
if nargin<6;OPT=1;end;
if nargin<5;method='autoscaling';end;
if nargin<4;K=10;end;
if nargin<3;A=2;end;

check=0; %+++ status variable:  1: Inf

if order==1
  [y,indexyy]=sort(y);
  X=X(indexyy,:);
else
  indexyy=randperm(length(y));
  X=X(indexyy,:);
  y=y(indexyy);
end

A=min([size(X,1)-ceil(length(y)/K) size(X,2) A]);
yytest=[];YR=[];
[Mx,Nx]=size(X);
groups = 1+rem(0:Mx-1,K);
yytest=[];yp=[];nLV=zeros(K,1);
for group=1:K
    testk = find(groups==group);  calk = find(groups~=group);
    Xcal=X(calk,:);ycal=y(calk);
    Xtest=X(testk,:);ytest=y(testk);

    CV=plsldacv(Xcal,ycal,A,K,method,0,order);
    if CV.check==1;check==1;break;end;
    LDA=plslda(Xcal,ycal,CV.optPC,method);
    ypred=plsldaval(LDA,Xtest,ytest);
    yytest=[yytest;ytest];
    yp=[yp;ypred;];
    nLV(group)=CV.optPC;
    if OPT==1;fprintf('The %dth outer loop finished.\n',group);end;
end

%+++ Find the most frequently chosen nLV.
uniLV=unique(nLV);
for j=1:length(uniLV); freq(j)=length(find(nLV==uniLV(j)));end
[maxf,maxindex]=max(freq);
optPC=uniLV(maxindex(1));

%+++ output
if check==0
  F=roccurve(yytest,yp,0);
  error=sum(sign(yp)~=yytest)/Mx;
  DCV.method=method;
  DCV.check=check;
  DCV.error=error;
  DCV.Sensitivity=F.sensitivity;
  DCV.Specificity=F.specificity;
  DCV.nLV=nLV;
  DCV.optPC=optPC;
elseif check==1
  DCV.method=method;
  DCV.check=check;
end

author: viggin

repo: domain-adaptation-toolbox


text:
 function x = nufft2_adj(X, st)
%function x = nufft2_adj(X, st)
%	SUPERCEDED BY nufft_adj.m
%	Apply adjoint of 2D NUFFT to vector(s) X
%	in:	X	[M,L]
%		st		structure precomputed by nufft2_init()
%	out:	x	[N1,N2,L]
%
%	Copyright 2001-9-17	Jeff Fessler	The University of Michigan

%
%	if no arguments, then run a simple test
%
if nargin < 1
	help(mfilename)

	N1 = 4;
	N2 = 8;
	n_shift = [2.7 3.1];	% random shifts to stress it
	o1 = 2 * pi * [0.0 0.1 0.3 0.4 0.7 0.9]';
	o2 = flipud(o1);
	st = nufft2_init([o1 o2], N1, N2, 4, 6, 2*N1, 2*N2, n_shift, 0, 'best');

	M = length(o1);
	A = zeros(M, N1*N2);
	for n1=1:N1
		for n2=1:N2
			x = zeros(N1,N2); x(n1,n2) = 1;
			jj = n1 + (n2-1) * N1;
			A(:,jj) = nufft2(x, st);
		end
	end

	Aa = zeros(N1*N2, M);
	for ii=1:M
		y = zeros(M,1); y(ii,1) = 1;
		Aa(:,ii) = col(nufft2_adj(y, st));
	end

	printf('max %% = %g', max_percent_diff(A,Aa'))
	clear x
return
end

%	extract attributes from structure
K1 = st.K1;
K2 = st.K2;
if size(X,1) ~= st.M, error size, end

%
%	adjoint of interpolator using precomputed sparse matrix
%
Xk = full(st.p' * X);			% [K1*K2,L]
Xk = reshape(Xk, [K1 K2 size(X,2)]);	% [K1,K2,L]

x = K1 * K2 * ifft2(Xk);		% [K1,K2,L]

% eliminate zero padding from ends
x = x(1:st.N1,1:st.N2,:);		% [N1,N2,L]

if ndims(x) == 2
	x = x .* st.sn;
else
	error '3d not done'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef LevelSetCircle < LevelSetSphereNdim

    methods (Access = public)

        function obj = LevelSetCircle(cParams)
            obj.fracRadius = cParams.fracRadius;
            obj.compute(cParams);
        end
    end

    methods (Access = protected)
        function computeLevelSetValue(obj)
            ls = obj.dist - 1;
            obj.levelSet = ls;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function M_hat = run_mc(params)
  M = params.M;
  Idx = params.Idx;
  [numr,numc] = size(M);
  rank = 10;
  Known = find(Idx);
  data = M(Known);
  [X,Y] = lmafit_mc_adp(numr,numc,rank,Known,data,[]);
  M_hat = X*Y;
end

author: andrewssobral

repo: lrslibrary


text:
% A test script for the Half-vectorization package

% Random matrix
A=ceil(20*rand(5,7))

% Numerotation of diagonals
disp('Numerotation of diagonals');
B=zeros(size(A));
[I J]=itril(size(A),Inf);
B(itril(size(A),Inf))=J-I;
disp(B);

% Get the lower part, shift down by 2
disp('Get the lower part, shift down by 2...');

% Using Matlab
T = tril(ones(size(A)),-2);
HalfV = A(T==1);
disp('Using Matlab TRIL');
disp(HalfV);

% Using c
I = itril(size(A),-2);
HalfV = A(I);
disp('Using indices returned by itril');
disp(HalfV);

% Using Elimination Matrix
HaftV = EliminationM(size(A),-2)*A(:);
disp('Using elimination matrix');
disp(HalfV);

% Replace the diagonal
disp('Replace the diagonal');
Idiag = idiag(size(A));
A(Idiag) = -A(Idiag);
disp(A);

% TOEPLITZ
c = ceil(20*rand(5,1))
r = ceil(20*rand(7,1))
T = mytoeplitz(c, r)

author: Sable

repo: mcbench-benchmarks


text:
function varargout = fastind2sub(siz,ndx)
%"fastind2sub"
%   FAST ind2sub is a faster version of ind2sub, based off the ind2sub that
%   ships with matlab.  Only one line is changed, and the original doc for
%   ind2sub is included below.
%
%   JRA 2/26/04
%
%IND2SUB Multiple subscripts from linear index.
%
%   IND2SUB is used to determine the equivalent subscript values
%   corresponding to a given single index into an array.
%
%   [I,J] = IND2SUB(SIZ,IND) returns the arrays I and J containing the
%   equivalent row and column subscripts corresponding to the index
%   matrix IND for a matrix of size SIZ.
%   For matrices, [I,J] = IND2SUB(SIZE(A),FIND(A>5)) returns the same
%   values as [I,J] = FIND(A>5).
%
%   [I1,I2,I3,...,In] = IND2SUB(SIZ,IND) returns N subscript arrays
%   I1,I2,..,In containing the equivalent N-D array subscripts
%   equivalent to IND for an array of size SIZ.
%
%   See also SUB2IND, FIND.

%   Copyright 1984-2002 The MathWorks, Inc.
%   $Revision: 1.3 $  $Date: 2008/10/08 21:38:59 $
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

nout = max(nargout,1);
if length(siz)<=nout,
  siz = [siz ones(1,nout-length(siz))];
else
  siz = [siz(1:nout-1) prod(siz(nout:end))];
end
n = length(siz);
k = [1 cumprod(siz(1:end-1))];
ndx = ndx - 1;
for i = n:-1:1,
  varargout{i} = floor(ndx/k(i))+1;
  ndx = ndx - (varargout{i}-1)*k(i);

%  This is the old method: slow because it
%  recalculates the division after it has already been found above.
%  ndx = rem(ndx,k(i));
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ cluster_center, seed ] = cluster_initialize_3 ( dim_num, ...
  point_num, cluster_num, point, seed )

%*****************************************************************************80
%
%% CLUSTER_INITIALIZE_3 initializes the cluster centers to random values.
%
%  Discussion:
%
%    In this case, each point is randomly assigned to a cluster, and
%    the cluster centers are then computed as the centroids of the points
%    in the cluster.
%
%    Random integers for this function can be generated by the MATLAB
%    function randi(), or by the source code i4_uniform();
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates
%    of the points.
%
%    Input, integer SEED, a seed for the random
%    number generator.
%
%    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the coordinates of the cluster centers.
%
%    Output, integer SEED, a seed for the random
%    number generator.
%

%
%  Assign one point to each cluster center.
%
  for i = 1 : cluster_num
    cluster_center(1:dim_num,i) = point(1:dim_num,i);
  end

  cluster_population(1:cluster_num) = 1;
%
%  The rest of the points get assigned randomly.
%
  for i = cluster_num + 1 : point_num
    if ( 0 )
      [ j, seed ] = i4_uniform ( 1, cluster_num, seed );
    else
      j = randi ( [ 1, cluster_num ], 1, 1 );
    end
    cluster_center(1:dim_num,j) = cluster_center(1:dim_num,j) + point(1:dim_num,i);
    cluster_population(j) = cluster_population(j) + 1;
  end
%
%  Now average the points to get the centroid.
%
  for i = 1 : cluster_num
    cluster_center(1:dim_num,i) = cluster_center(1:dim_num,i) ...
      / cluster_population(i);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function totalEnergy = getEnergyDeposited(planC,structNum,doseNum)
%function getEnergyDeposited(planC,structNum,doseNum)
%
%This function computes the total Energy deposited within a structure.
%
%APA, 04/21/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS
stateS.optS = opts4Exe([getCERRPath,'CERROptions.json']);
stateS.MLVersion = getMLVersion;

indexS = planC{end};

scanNum = getStructureAssociatedScan(structNum,planC);

%Get Uniformized dose matrix within the structure structNum
dose3M = getUniformDose(doseNum, scanNum, structNum, planC);

numVoxels = length(find(dose3M > 0));

%Get Uniformized scan matrix for scanNum
scan3M = getUniformizedCTScan(1, scanNum, planC);

[xV,yV,zV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));

dx = abs(xV(1)-xV(2));
dy = abs(yV(1)-yV(2));
dz = abs(zV(1)-zV(2));

voxelVolume = dx * dy * dz;

totalEnergy = sum(dose3M(:).*scan3M(:))*voxelVolume*numVoxels;


text:
function p = normpdfln(x, m, S, V)
% NORMPDFLN    log of multivariate normal density.
%   See NORMPDF for argument description.

log2pi = 1.83787706640935;
[d, n] = size(x);
if nargin == 1
  dx = x;
elseif isempty(m)
  dx = x;
else
  % m specified
  sz = size(m);
  if sz(1) ~= d
    error('rows(m) ~= rows(x)')
  end
  nm = sz(2);
  if nm == 1
    dx = x - repmat(m,1,n);
  elseif n == 1
    dx = repmat(x,1,nm) - m;
  elseif nm == n
    dx = x - m;
  else
    error('incompatible number of columns in x and m')
  end
end
if nargin < 3
  % unit variance
  p = -0.5*(d*log2pi + col_sum(dx.*dx));
  return
end
have_inv = 0;
if nargin == 3
  % standard deviation given
  if d == 1
    dx = dx./S;
    p = (-log(S) -0.5*log2pi) - 0.5*(dx.*dx);
    return;
  end
  if S(2,1) ~= 0
    error('S is not upper triangular')
  end
  if any(size(S) ~= [d d])
    error('S is not the right size')
  end
else
  if ischar(V)
    if strcmp(V,'inv')
      % inverse stddev given
      iS = S;
      have_inv = 1;
    else
      error('unknown directive')
    end
  elseif ischar(S)
    if strcmp(S,'inv')
      % inverse variance given
      if d == 1
	iS = sqrt(V);
      else
	iS = chol(V);
      end
      have_inv = 1;
    else
      error('unknown directive')
    end
  else
    % variance given
    if d == 1
      S = sqrt(V);
    else
      S = chol(V);
    end
  end
end
if have_inv
  if d == 1
    dx = iS .* dx;
    logdetiS = log(iS);
  else
    dx = iS*dx;
    logdetiS = sum(log(diag(iS)));
  end
else
  if d == 1
    dx = dx./S;
    logdetiS = -log(S);
  else
    dx = solve_tril(S',dx);
    %dx = S'\dx;
    logdetiS = -sum(log(diag(S)));
  end
end
p = (logdetiS -0.5*d*log2pi) -0.5*col_sum(dx.*dx);

author: Cloud-CV

repo: object-proposals


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

% Useful variables
[m, n] = size(X);

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu = mean(X);
sigma2 = var(X, 1);

mu = mu(:);
sigma2 = sigma2(:);

% =============================================================

end

author: zlotus


text:
function c = opoly(x,y,n)
%
%   opolyfit     OPOLYFIT(x,y,n) finds the n coefficients of a polynomial
%		 formed from the data in row vector x that fits the data
%		 in row vector y in a least-squares sense. Orthogonal
%		 polynomials are used to improve accuracy. See also polyfit.
%
%		 c = opolyfit(x,y,n)
%		 c  ---	coefficient vector, size n
%
%		 References:
%		 G. E. Forsythe, "Generation and use of orthogonal poly-
%		 nomials for data fitting on a digital computer",
%		 J. SIAM, 5(1957), 74-88.
%
%		 P. F. Hultquist, "Numerical Methods for Engineers and
%		 Computer Scientists", Benjamin-Cummings, 1988.
%
%		 Copyright 1992. Gordon	J. Frazer, Signal Processing Research
%		 Laboratory, and Roy J. Hughes, Laboratory for Instrumental and
%		 Developmental Chemistry, QUT, Brisbane, AUSTRALIA, 1992.
%		 Email:	g.frazer@qut.edu.au or	r.hughes@qut.edu.au
n=n+1; % To make it similar to MATLAB polyfit.m

if n<1
 error('polynomial order too small, must be greater than or equal to one');
end
if n+1>length(x)
 error('polynomial order too large for number of data points');
end
% normalise y to improve numerical stability
yscale=0.5/max(abs(y));
yy=yscale*y;
% initialise opoloyfit
gamma_0=length(x);
alpha_0=mean(x);
PC=zeros(n);
PC(1,1)=1;
a=zeros(1,n);
a(1)=mean(yy);
% calculate the	first iteration	manually
PC(2,1:2)=conv([1,-alpha_0],PC(1,1));
gamma_k=sum(polyval(PC(2,1:2),x).*polyval(PC(2,1:2),x));
beta_k=gamma_k/gamma_0;
alpha_k=sum(x.*(polyval(PC(2,1:2),x).*polyval(PC(2,1:2),x)))/gamma_k;
a(2)=sum(yy.*polyval(PC(2,1:2),x))/gamma_k;
% calculate for	the remaining polynomial orders
for k=2:n-1
 PC(k+1,1:k+1)=conv([1,-alpha_k],PC(k,1:k))-beta_k*[0,0,PC(k-1,1:k-1)];
 gamma_km1=gamma_k;
 gamma_k=sum(polyval(PC(k+1,1:k+1),x).*polyval(PC(k+1,1:k+1),x));
 beta_k=gamma_k/gamma_km1;
 alpha_k=sum(x.*(polyval(PC(k+1,1:k+1),x).*polyval(PC(k+1,1:k+1),x)))/gamma_k;
 a(k+1)=sum(yy.*polyval(PC(k+1,1:k+1),x))/gamma_k;
end
% combine the a's to simplify the polynomial
for i=1:n
 PC(i,:)=a(i)*PC(i,:);
end
for i=1:n
 cc(i)=sum(diag(PC,i-n));
end
% unnormalise c	by the y scale factor
c=(1/yscale)*cc;

author: HeLiangHIT


text:
function [bias av_bias] = bias_f(x,y)

% Bias calculator

% Formula: 1 - mean(fused image)/mean(original image)
% (Ideal value = 0)

% 07/03/2010 Version 1.0
% 25/06/2010 Version 1.2    - Excel Output option
% 04/08/2011 Version 1.2F   - Function Version

% Author: Aristidis D. Vaiopoulos

% Find the number of bands
bands = size(x);
if length(bands) == 3
    bands = bands(1,3);
else
    bands = 1;
end
% Preallocation
mx = zeros(1,bands);
my = zeros(1,bands);
% Mean value calculation
for i = 1:bands
    xt = double(x(:,:,i));
    yt = double(y(:,:,i));
    mx(i) = mean(xt(:));
    my(i) = mean(yt(:));
end
% Bias calculation
bias = 1 - (my./mx);
bias = bias';

av_bias = mean(bias);

end


text:
function [ grid_weight, grid_point ] = sparse_grid ( dim_num, level_max, ...
  rule, point_num )

%*****************************************************************************80
%
%% SPARSE_GRID computes a sparse grid.
%
%  Discussion:
%
%    A Smolyak construction is used to create a multidimensional sparse grid.
%
%    The user specifies:
%    * the spatial dimension of the quadrature region,
%    * the level that defines the Smolyak grid.
%    * the 1D quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, controls the size of the final
%    sparse grid.
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Input, integer POINT_NUM, the number of points in the grid,
%    as determined by LEVELS_INDEX_SIZE.
%
%    Output, real GRID_WEIGHT(POINT_NUM), the weights.
%
%    Output, real  GRID_POINT(DIM_NUM,POINT_NUM), the points.
%
  if ( rule == 1 )
    [ grid_weight, grid_point ] = sparse_grid_cfn ( dim_num, level_max, ...
      rule, point_num );
  elseif ( 2 <= rule & rule <= 4 )
    [ grid_weight, grid_point ] = sparse_grid_ofn ( dim_num, level_max, ...
      rule, point_num );
  elseif ( 5 <= rule & rule <= 6 )
    [ grid_weight, grid_point ] = sparse_grid_own ( dim_num, level_max, ...
      rule, point_num );
  elseif ( 7 == rule )
    [ grid_weight, grid_point ] = sparse_grid_onn ( dim_num, level_max, ...
      rule, point_num );
  else
    grid_weight = [];
    grid_point = [];
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPARSE_GRID - Fatal error!\n' );
    fprintf ( 1, '  Illegal input rule index = %d\n', rule );
    error ( 'SPARSE_GRID - Fatal error!' );
  end

  return
end


text:
function errstring = consistent(ds, type)

% CONSISTENT   Check ReBEL data structures for consistentency.
%
%   errstring = consistent(ds, type) checks a data structure DS of supposed TYPE for consistency.
%
%     See also
%     GSSM, GENINFDS
%
%   Copyright (c) Oregon Health & Science University (2006)
%
%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
%   academic use only (see included license file) and can be obtained from
%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
%   software should contact rebel@csee.ogi.edu for commercial licensing information.
%
%   See LICENSE (which should be part of the main toolkit distribution) for more
%   detail.

%===============================================================================================

if (nargin~=2)
  error('Not enough input arguments!');
end

if ~ischar(type)
  error('Type argument should be a string.');
end

% Assume that all is OK as default
errhead = ['  [ consistent(' inputname(1) ',''' type ''') ]   '];

errstring = '';

if ~isstruct(ds)
  errstring = [errhead inputname(1) ' is not a data structure.'];
  return
end

switch (type)

%===============================================================================================================================
%== General State Space Model Data Structure (GSSM)
%===============================================================================================================================
case 'gssm'

  %-- Check existence of required data structure fields
  mf = checkstructfields(ds,'type','tag','statedim','obsdim','paramdim','U1dim','U2dim','Vdim','Ndim','pNoise','oNoise','params','ffun','hfun','setparams');
  if ~isempty(mf), errstring = [errhead ' Data structure does not contain the following fields : ' mf]; return; end

  %-- Check type field
  if ~ischar(ds.type) | ~stringmatch(ds.type,'gssm')
      errstring = [errhead  ' Data structure is not of type ''' type ''' or the type field is not a string'];
      return;
  end

%===============================================================================================================================
%== Inference Data Structure (InferenceDS)
%===============================================================================================================================
case 'InferenceDS'

  %-- Check existence of required data structure fields
  mf = checkstructfields(ds,'type','inftype','statedim','obsdim','U1dim','U2dim','Vdim','Ndim', ...
                            'ffun','hfun','innovation','stateAngleCompIdxVec','obsAngleCompIdxVec');
  if ~isempty(mf), errstring = [errhead ' Data structure does not contain the following fields : ' mf]; return; end

  %-- Check type field
  if (~ischar(ds.type) | ~stringmatch(ds.type,'InferenceDS'))
    errstring = [errhead  ' Data structure is not of type ''' type ''' or the type field is not a string.'];
    return;
  end

  %-- Check inference type field
  if ~stringmatch(ds.inftype,{'state','parameter','joint'})
    errstring = [errhead  ' Inference type ''' type ''' unknown.'];
    return;
  end

%===============================================================================================================================
%== Noise Source Data Structure (NoiseDS)
%===============================================================================================================================
case 'NoiseDS'

  %-- Check existence of required data structure fields
  mf = checkstructfields(ds,'type','ns_type','tag','dim','sample','likelihood');
  if ~isempty(mf), errstring = [errhead ' Data structure does not contain the following fields : ' mf]; return; end

  %-- Check type field
  if (~ischar(ds.type) | ~stringmatch(ds.type,'NoiseDS'))
    errstring = [errhead  ' Data structure is not of type ''' type ''' or the type field is not a string.'];
    return;
  end

%===============================================================================================================================
otherwise
  errstring = [errhead 'Data structure type ' type ' not supported.'];
end

return

repo: ekfmonoslam


text:
classdef ImageRegr < Algorithms.CDAlg
    methods
        function obj = ImageRegr()
            obj@Algorithms.CDAlg('Image Regression', '')
        end
        [DI, coeffs] = detectChange(obj, ~, ~);
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
load ../graphs/bfs_example.mat
[d dt pred] = bfs(A,2);
[ignore order] = sort(dt);
labels(order)
treeplot(pred);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
  function isave = iter_saver(isave, niter)
%|function isave = iter_saver(isave, niter)
%| process 'isave' and 'niter' options of iterative algorithms
%| supporting char options 'all' and 'last' (default if empty)

if nargin < 2, ir_usage, end

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = false;
end

if isempty(isave)
	if ~warned
		warn 'using default isave "last"'
		warned = 1;
	end
	isave = 'last';
end

if streq(isave, 'last')
	isave = niter;
elseif streq(isave, 'all')
	isave = 0:niter;
end

if any(isave < 0) || any(isave > niter)
	pr isave
	pr niter
	fail('bad isave')
end


text:
function mask = voxel2mask(voxels,maskdims)
% :Usage:
% ::
%
%     function mask = voxel2mask(voxels, x y z mask dimensions)
%
% :Voxels:
%   - 3 column vectors
%   - [i j k] = row, column, slice
%   - [x y z] in brain if brain is in analyze format
%
%     (x is rows, y is columns, z is slices)
%
% ..
%    % Tor Wager, 10/17/01
% ..

[n, m] = size(voxels);
wh_bad = false(n, 1);

% Check for illegal voxels
if m ~= 3, warning('voxel2mask: Illegal voxel input list', 'Voxels must be k x 3 matrix'); end

mv = max(voxels);
if any(mv > maskdims)
    warning('voxel2mask: Illegal voxel input list', 'Voxels outside mask, Voxel indices > mask dims');
    wh_bad = any(voxels > mv(ones(n, 1), :), 2);
end

mv = min(voxels);
if any(mv < 1)
    warning('voxel2mask: Illegal voxel input list', 'Voxels outside mask, Voxel indices < 1');
    wh_bad = [wh_bad | any(voxels < ones(n, 3), 2)];
end

if ~isempty(wh_bad)
    wh_bad = find(wh_bad);
    disp('You need to check your images!');
    disp('Offending voxel numbers: '); disp(wh_bad)
    disp('Offending voxel coordinates: '); disp(voxels(wh_bad, :))
    voxels(wh_bad, :) = [];
end

mask = zeros(maskdims);

for i = 1:size(voxels, 1)
    mask(voxels(i,1),voxels(i,2),voxels(i,3)) = 1;
end

mask = double(mask);

return


text:
function representa(robot,qt1,qt2,num,coord,num2)
%REPRESENTA LOS DOS ROBOTS Y LA CAJA SIGUIENDO LA TRAYECTORIA

for i=1:size(qt1,2)

    if num==2

        T= directkinematic(robot.robot1, qt1(:,i));
        Pc = T(1:3,4);
        tangente=atan(T(2,3)/T(1,3));
        robot.robot1.equipment{1}.T0(1:3,4) = Pc;
        robot.robot1.equipment{1}.T0(3,4)=robot.robot1.equipment{1}.T0(3,4)-coord(3,1);
        if num2==1

         robot.robot1.equipment{1}.T0(2,4)=robot.robot1.equipment{1}.T0(2,4)+norm(coord(1:2,1))*sin(tangente);
         robot.robot1.equipment{1}.T0(1,4)=robot.robot1.equipment{1}.T0(1,4)+norm(coord(1:2,1))*cos(tangente);

        robot.robot1.equipment{1}.T0(1:3,1:3)=T(1:3,1:3)*[0 -1  0;
                                                          1  0  0;
                                                          0  0  1];
        else
           robot.robot1.equipment{1}.T0(2,4)=robot.robot1.equipment{1}.T0(2,4)-norm(coord(1:2,1))*sin(tangente);
         robot.robot1.equipment{1}.T0(1,4)=robot.robot1.equipment{1}.T0(1,4)+norm(coord(1:2,1))*cos(tangente);
        end

    end
      draw2robots(robot.robot1,qt1(:,i),robot.robot2,qt2(:,i))

end

end

repo: ARTE


text:
%------------------------------------------------------------------------
% AL-AZHAR UNIVERSITY
% FACULTY OF ENGINEERING
% SYSTEMS & COMPUTERS ENGINEERING Department
%------------------------------------------------------------------------
% Author : Ahmed Samieh Abd El-Wahab
% Date   : 14 December 2006
%------------------------------------------------------------------------
% Shapes Classifier
% Step 1: Read image Read in
% Step 2: Convert image from rgb to gray
% Step 3: Threshold the image
% Step 4: Invert the Binary Image
% Step 5: Find the boundaries Concentrate
% Step 6: Determine Shapes properties
% Step 7: Classify Shapes according to properties
% Square      = 3
% Rectangular = 2
% Circle      = 1
% UNKNOWN     = 0
%------------------------------------------------------------------------
function W = Classify(ImageFile)
% Step 1: Read image Read in
RGB = imread(ImageFile);
figure,
imshow(RGB),
title('Original Image');

% Step 2: Convert image from rgb to gray
GRAY = rgb2gray(RGB);
figure,
imshow(GRAY),
title('Gray Image');

% Step 3: Threshold the image Convert the image to black and white in order
% to prepare for boundary tracing using bwboundaries.
threshold = graythresh(GRAY);
BW = im2bw(GRAY, threshold);
figure,
imshow(BW),
title('Binary Image');

% Step 4: Invert the Binary Image
BW = ~ BW;
figure,
imshow(BW),
title('Inverted Binary Image');

% Step 5: Find the boundaries Concentrate only on the exterior boundaries.
% Option 'noholes' will accelerate the processing by preventing
% bwboundaries from searching for inner contours.
[B,L] = bwboundaries(BW, 'noholes');

% Step 6: Determine objects properties
STATS = regionprops(L, 'all'); % we need 'BoundingBox' and 'Extent'

% Step 7: Classify Shapes according to properties
% Square = 3 = (1 + 2) = (X=Y + Extent = 1)
% Rectangular = 2 = (0 + 2) = (only Extent = 1)
% Circle = 1 = (1 + 0) = (X=Y , Extent < 1)
% UNKNOWN = 0

figure,
imshow(RGB),
title('Results');
hold on
for i = 1 : length(STATS)
  W(i) = uint8(abs(STATS(i).BoundingBox(3)-STATS(i).BoundingBox(4)) < 0.1);
  W(i) = W(i) + 2 * uint8((STATS(i).Extent - 1) == 0 );
  centroid = STATS(i).Centroid;
  switch W(i)
      case 1
          plot(centroid(1),centroid(2),'wO');
      case 2
          plot(centroid(1),centroid(2),'wX');
      case 3
          plot(centroid(1),centroid(2),'wS');
  end
end
return


text:
function test_bug1150

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_sourcestatistics

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1150.mat'));

% the following failed
% stat = ft_sourcestatistics(cfg, temp);
% which could be tracked down to the dimord being incorrect after
% ft_checkdata

% the data seems to be incosistent so the conversion fails.
temp.trial = temp.trial(1:10);
temp.cumtapcnt = temp.cumtapcnt(1:10);
temp.trialinfo = temp.trialinfo(1:10,:);

sourcenew  = ft_checkdata(temp, 'sourcerepresentation', 'new');
if ~isequal(size(sourcenew.pow), [38556 10])
  error('incorrect dimensions');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ w, x ] = gm_rule_set_old ( rule, dim_num, point_num )

%*****************************************************************************80
%
%% GM_RULE_SET_OLD sets a Grundmann-Moeller rule.  (OBSOLETE VERSION)
%
%  Discussion:
%
%    This version of the computation is no longer used.  The direct
%    application of the formula results in overflows and inaccuracies
%    very quickly.
%
%    This rule returns weights and abscissas of a Grundmann-Moeller
%    quadrature rule for the DIM_NUM-dimensional unit simplex.
%
%    The dimension POINT_NUM can be determined by calling GM_RULE_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Axel Grundmann, Michael Moeller,
%    Invariant Integration Formulas for the N-Simplex
%    by Combinatorial Methods,
%    SIAM Journal on Numerical Analysis,
%    Volume 15, Number 2, April 1978, pages 282-290.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%    0 <= RULE.
%
%    Input, integer DIM_NUM, the spatial dimension.
%    1 <= DIM_NUM.
%
%    Input, integer POINT_NUM, the number of points in the rule.
%
%    Output, real W(POINT_NUM), the weights.
%
%    Output, real X(DIM_NUM,POINT_NUM), the abscissas.
%
  s = rule;
  d = 2 * s + 1;
  k = 0;
  n = dim_num;
  one_pm = 1;

  for i = 0 : s

    weight = r8_factorial ( n ) ...
      * ( one_pm *  ( d + n - 2 * i )^d ) ...
      / ( 2^(2*s) * r8_factorial ( i ) * r8_factorial ( d + n - i ) );

    one_pm = - one_pm;

    beta_sum = s - i;
    more = 0;
    beta = [];
    h = 0;
    t = 0;

    while ( 1 )

      [ beta, more, h, t ] = comp_next ( beta_sum, dim_num + 1, ...
        beta, more, h, t );

      k = k + 1;

      w(k) = weight;

      x(1:dim_num,k) = ( 2 * beta(2:dim_num+1)' + 1 ) / ( d + n - 2 * i );

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function B = fevalArrays( A, fHandle, varargin )
% Used to apply the same operation to a stack of array elements.
%
% The only constraint on the function specified in fHandle is that given
% two differrent input arrays a1 and a2, if a1 and a2 have the same
% dimensions then the outputs b1 and b2 must have the same dimensions. For
% long operations shows progress information.
%
% A can have arbitrary dimension.  Suppose A has size d1 x d2 ... x dn.
% Then for the purpose of this function A has dn elements, where
% A(:,:,...,i) is the ith element.  This function then applies the
% operation in fHandle, with paramters given in varargin, to each element
% in A. The results are returned in the array B, of size f1 x f2 x ... x fk
% x dn.  Each of the n element of B of the form B(:,:,...,i) is the the
% result of applying fHandle to A(:,:,...,i).  A may also be a cell array,
% see the last example.
%
% A limitation of fevalArrays is that it does not pass state information
% to fHandle.  For example, fHandle may want to know how many times it's
% been called.  This can be overcome by saving state information inside
% fHandle using 'persistent' variables.  For an example see imwrite2.
%
% USAGE
%  B = fevalArrays( A, fHandle, varargin )
%
% INPUTS
%  A        - input array
%  fHandle  - operation to apply to each 'element' of A
%  params   - [varargin] parameters for each operation specified by fHandle
%
% OUTPUTS
%  B        - output array
%
% EXAMPLE
%  B = fevalArrays( A, @rgb2gray );      % where A is MxNx3xR
%  B = fevalArrays( A, @imresize, .5 );  % where A is MxNxR
%  B = fevalArrays( A, @imNormalize );   % where A has arbitrary dims
%  B = fevalArrays( A, @(x) {imresize(x{1},.5)} ); % A is cell array
%
% See also FEVALIMAGES, FEVALMATS, ARRAYFUN
% IMWRITE2, PERSISTENT, TICSTATUS,
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

nd = ndims(A);  siz = size(A);  n = siz(end);
indsA = {':'}; indsA = indsA(ones(nd-1,1));

ticId = ticStatus('fevalArrays',[],60);
for i=1:n
  % apply fHandle to each element of A
  b = feval( fHandle, A(indsA{:},i), varargin{:} );
  if( i==1 )
    ndb = ndims(b);
    if(ndb==2 && size(b,2)==1); ndb=1; end
    onesNdb = ones(1,ndb);
    B = repmat( b, [onesNdb,n] );
    indsB = {':'}; indsB = indsB(onesNdb);
  else
    B(indsB{:},i) = b;
  end;
  tocStatus( ticId, i/n );
end

author: zouchuhang

repo: LayoutNet


text:
function [gnodes, gedges] = relativeNeighborhoodGraph(points)
%RELATIVENEIGHBORHOODGRAPH Relative Neighborhood Graph of a set of points
%
%   [NODES, EDGES] = relativeNeighborhoodGraph(POINTS)
%   EDGES = relativeNeighborhoodGraph(POINTS)
%
%   The Relative Neighborhood Graph (RNG) is a subgraph of the Delaunay
%   Triangulation computed from the same set of points. The Gabriel graph
%   and the euclidean minimal spanning tree (EMST) are subgraphs of the
%   RNG.
%
%   Example
%     nodes = rand(100, 2) * 100;
%     edges = relativeNeighborhoodGraph(nodes);
%     figure; drawGraph(nodes, edges);
%
%   See also
%     gabrielGraph, euclideanMST
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2016-03-02,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2016 INRA - Cepia Software Platform.

% first compute Delaunay triangulation to reduce further computations
DT = delaunayTriangulation(points);
E = edges(DT);

% compute edge lengths
nEdges = size(E, 1);
edgeLengths = zeros(nEdges, 1);
for i = 1:nEdges
    edgeLengths(i) = distancePoints(points(E(i,1),:), points(E(i,2),:));
end

% identify indices of faces attached to each vertex
vertexFaces = vertexAttachments(DT);

% iterate over edges to check if the should be kept
keepEdge = true(nEdges, 1);
for iEdge = 1:nEdges
    iVertex1 = E(iEdge, 1);
    iVertex2 = E(iEdge, 2);
    vertex1 = points(iVertex1, :);
    vertex2 = points(iVertex2, :);

    % compute indices of faces containing one of the two vertices
    inds = [vertexFaces{iVertex1} vertexFaces{iVertex2}];
    localFaces = DT.ConnectivityList(inds, :);

    % compute indices of vertices is the first neighborhood of the edge
    inds = unique(localFaces);
    inds(ismember(inds, [iVertex1 iVertex2])) = [];

    % compute max of distances to both original vertices
    dists1 = distancePoints(vertex1, points(inds, :));
    dists2 = distancePoints(vertex2, points(inds, :));
    distsMax = max(dists1, dists2);

    % keep edge if all points are outside the "lunule" defined by the edge
    if edgeLengths(iEdge) > min(distsMax)
        keepEdge(iEdge) = false;
    end
end

% filter edges
gedges = E(keepEdge, :);

% format output
gnodes = points;
if nargin == 1
    gnodes = gedges;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function tc = tc_deconvolvedAmps(tc);
% tc = tc_deconvolvedAmps(tc);
% compute the peak-bsl amplitudes of deconvolved
% time courses (and relamps, when I implement that):

%%%%%% get params
frameWindow = unique(round(tc.params.timeWindow./tc.TR));
prestim = -1 * frameWindow(1);
peakFrames = unique(round(tc.params.peakPeriod./tc.TR));
bslFrames = unique(round(tc.params.bslPeriod./tc.TR));
peakFrames = find(ismember(frameWindow,peakFrames));
bslFrames = find(ismember(frameWindow,bslFrames));

%%%%% norm baseline periods
nConds = size(tc.glm.betas, 2);

if tc.params.normBsl==1
    offset = mean(tc.glm.betas(bslFrames,:), 1);
    tc.glm.betas = tc.glm.betas - repmat(offset, [length(frameWindow) 1]);
end

%%%%% calc amplitudes
for i = 1:nConds
    bsl = tc.glm.betas(bslFrames, i);
    peak = tc.glm.betas(peakFrames, i);
    tc.glm.amps(:,i) = (nanmean(peak) - nanmean(bsl))';
    tc.glm.amp_sems(:,i) = nanmean(tc.glm.sems(:,i))';
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function selected = SelectSpikes(spikes,varargin)

%SelectSpikes - Discriminate bursts vs single spikes.
%
%  USAGE
%
%    selected = SelectSpikes(spikes,<options>)
%
%    spikes         spike times
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'mode'        either 'bursts' (default) or 'single'
%     'isi'         max inter-spike interval for bursts (default = 0.006),
%                   or min for single spikes (default = 0.020)
%    =========================================================================
%
%  OUTPUT
%
%    selected       a logical vector indicating for each spike whether it
%                   matches the criterion

% Copyright (C) 2011 by Michaël Zugaro
%
% 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.

% Default values
isiBursts = 0.006;
isiSpikes = 0.020;
mode = 'bursts';

% Check number of parameters
if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help SelectSpikes'' for details).');
end

% Check parameter size
if ~isdvector(spikes),
	error('Incorrect spikes (type ''help SelectSpikes'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+2) ' is not a property (type ''help SelectSpikes'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'mode',
      mode = varargin{i+1};
      if ~isstring_FMAT(mode,'bursts','single'),
        error('Incorrect value for property ''mode'' (type ''help SelectSpikes'' for details).');
      end
    case 'isi',
      isi = varargin{i+1};
      if ~isdscalar(isi,'>0'),
        error('Incorrect value for property ''isi'' (type ''help SelectSpikes'' for details).');
      end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SelectSpikes'' for details).']);
  end
end

t = spikes(:);
dt = diff(t);

if strcmp(mode,'bursts'),
	b = dt<isi;
	selected = logical([0;b])|logical([b;0]); % either next or previous isi < threshold
else
	s = dt>isi;
	selected = logical([0;s])&logical([s;0]); % either next or previous isi < threshold
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function geometry_test203 ( )

%*****************************************************************************80
%
%% TEST203 tests TETRAHEDRON_CENTROID_3D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  tetra = [ ...
     0.000000,  0.942809, -0.333333; ...
    -0.816496, -0.816496, -0.333333; ...
     0.816496, -0.816496, -0.333333; ...
     0.000000,  0.000000,  1.000000 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST203\n' );
  fprintf ( 1, '  For a tetrahedron in 3D,\n' );
  fprintf ( 1, '  TETRAHEDRON_CENTROID_3D computes the centroid;\n' );

  r8mat_transpose_print ( dim_num, 4, tetra, '  Tetrahedron vertices:' );

  centroid = tetrahedron_centroid_3d ( tetra );

  r8vec_print ( dim_num, centroid, '  Centroid:' );

  return
end


text:
function [H s phi T] = RSTLS(X1, X2, normalization)

% [H s phi T] = RSTLS(X1, X2, normalization)
%
% DESC:
% computes the RST transformation between the point pairs X1, X2
%
% VERSION:
% 1.0.1
%
% INPUT:
% X1, X2        = point matches (cartesian coordinates)
% normalization = true (default) or false to enable/disable point
%                 normalzation
%
% OUTPUT:
% H             = homography representing the RST transformation
% s             = scaling
% phi           = rotation angle
% T             = translation vector

% AUTHOR:
% Marco Zuliani, email: marco.zuliani@gmail.com
% Copyright (C) 2011 by Marco Zuliani
%
% LICENSE:
% This toolbox is distributed under the terms of the GNU GPL.
% Please refer to the files COPYING.txt for more information.

% HISTORY
% 1.0.0         08/27/08 - intial version
% 1.0.1         06/09/09 - implemented closed form for the LS estimation
%                          routines

if (nargin < 3)
    normalization = true;
end;

N = size(X1, 2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% checks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (size(X2, 2) ~= N)
    error('RSTLS:inputError', ...
        'The set of input points should have the same cardinality')
end;
if N < 2
    error('RSTLS:inputError', ...
        'At least 2 point correspondences are needed')
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normalize the input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (normalization) && (N > 2)
    % fprintf('\nNormalizing...')
    [X1, T1] = normalize_points(X1);
    [X2, T2] = normalize_points(X2);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (N == 2)

    % fast estimation
    Theta = zeros(4,1);

    % $\mbox{\texttt{MM}} \eqdef M_{:,1} = \vct{y}^{(1)}-\vct{y}^{(2)} = \left[\begin{array}{c} y_1^{(1)}-y_1^{(2)} \\ y_2^{(1)}-y_2^{(2)} \end{array}\right]$
    % 2 additions
    MM = X1(:,1) - X1(:,2);
    % $ \mbox{\texttt{detMM}} \eqdef |M|$
    % 1 additions, 2 multiplication
    detMM = MM(1)*MM(1) + MM(2)*MM(2);
    % $ \mbox{\texttt{MMi}} \eqdef \left[ \begin{array}{c} \left[M^{-1}\right]_{1,1} \\ -\left[M^{-1}\right]_{2,1}\end{array}\right]$
    % 2 multiplications
    MMi = MM / detMM;

    % $ \mbox{\texttt{Delta}} \eqdef \vct{T}_{\vct{\theta}} (\vct{y}^{(1)})-\vct{T}_{\vct{\theta}} (\vct{y}^{(2)})$
    % 2 additions
    Delta = X2(:,1) - X2(:,2);

    % $ \mbox{\texttt{Theta(1:2)}} = M^{-1}\left(\vct{T}_{\vct{\theta}} (\vct{y}^{(1)})-\vct{T}_{\vct{\theta}} (\vct{y}^{(2)})\right)$
    % 1 additions, 2 multiplications
    Theta(1) = MMi(1)*Delta(1) + MMi(2)*Delta(2);
    % 1 additions, 2 multiplications
    Theta(2) = MMi(1)*Delta(2) - MMi(2)*Delta(1);
    % $ \mbox{\texttt{Theta(3:4)}} = -S^{(2)}\vct{\theta}_{1:2}+\vct{T}_{\vct{\theta}} (\vct{y}^{(2)})$
    % 2 additions, 2 multiplications

    Theta(3) = X2(1,2) - Theta(1)*X1(1,2) + Theta(2)*X1(2,2);
    % 2 additions, 2 multiplications
    Theta(4) = X2(2,2) - Theta(1)*X1(2,2) - Theta(2)*X1(1,2);

    % total: 11 additions, 12 multiplications
else

    % Closed form LS solution. Using the tutorial notation.

    % Notation semplification:
    % $\vct{p}^{(i)} = \bar{\vct{y}}^{(i)}$ and $\vct{q}^{(i)} = \overline{\vct{T}_{\vct{\theta}} (\vct{y}^{(i)})}$
    % $a = \sum_{i=1}^N\left( (p_1^{(i)})^2 + (p_2^{(i)})^2 \right)$
    a = sum(X1(:).^2);

    % Explicit LS expansion:
    % $  \theta_1 = \frac{1}{a}\sum_{i=1}^N p_1^{(i)} q_1^{(i)} + p_2^{(i)} q_2^{(i)} $
    % $  \theta_2 = \frac{1}{a}\sum_{i=1}^N -p_2^{(i)} q_1^{(i)} + p_1^{(i)} q_2^{(i)} $
    % $  \theta_3 = 0 $
    % $  \theta_4 = 0 $
    Theta(1) = sum( X1(1, :).*X2(1, :) + X1(2, :).*X2(2, :) ) / a;
    Theta(2) = sum( -X1(2, :).*X2(1, :) + X1(1, :).*X2(2, :) ) / a;
    Theta(3) = 0;
    Theta(4) = 0;

    % Traditional LS
    %
    %     A = zeros(2*N, 4);
    %     b = zeros(2*N, 1);
    %
    %     ind = 1:2;
    %     for n = 1:N
    %
    %         A(ind, 1:2) = [X1(1,n) -X1(2,n); X1(2,n) X1(1,n)];
    %         A(ind, 3:4) = eye(2);
    %
    %         b(ind) = X2(1:2, n);
    %
    %         ind = ind + 2;
    %
    %     end;
    %
    %     % solve the linear system in a least square sense
    %     Theta = A\b;

end;

% compute the corresponding homography
H = [Theta(1) -Theta(2) Theta(3); Theta(2) Theta(1) Theta(4); 0 0 1];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% de-normalize the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (normalization) && (N > 2)
    H = T2\H*T1;
end;
H = H/H(9);

% prepare the output
if nargout > 1

    s       = sqrt(H(1,1)*H(1,1) + H(2,1)*H(2,1));
    phi     = atan2(H(2,1), H(1,1));
    T       = H(1:2, 3);

end;

return

author: RANSAC

repo: RANSAC-Toolbox


text:
function aedat = ImportSecGen2(aedat)
%{
This function works similarly to ImportAedat, but takes instead data from
SEC DVS Gen2. These .bin files have no header info, therefore they are not
compatible with ImportAedat.

Expects the aedat.importParams.filePath parameter to contain a .bin file.
If, however the suffix is not .bin, it assumes that it has been passed a
folder path and takes the most recent file in the folder.
%}

dbstop if error

% Unpack parameters
importParams = aedat.importParams;

%% Open the file

if ~isfield(importParams, 'filePath') || ~strcmp(importParams.filePath(end -3 : end), '.bin')
    % Open the newest file
    % Note: For this to work at the moment, matlab must already be in the
    % correct directory
    directoryListing = dir('*.bin');
    [~, fileIndex] = sort([directoryListing.datenum]);
    filePath = directoryListing(fileIndex(end)).name;
else
    filePath = importParams.filePath;
end

fileHandle = fopen(filePath, 'r');

if fileHandle == -1
    error('file not found')
end

% Go to the EOF to find out how long it is
fseek(fileHandle, 0, 'eof');

% Find the number of packets
numPacketsInFile = ftell(fileHandle) / 4;

%% Check parameters

% Check the startEvent and endEvent parameters - if present, this shall
% actually refer to packets.
if isfield(importParams, 'startPacket')
    startPacket = importParams.startPacket;
elseif isfield(importParams, 'startEvent')
    startPacket = importParams.startEvent;
else
    startPacket = 1;
end
if startPacket > numPacketsInFile
	error([	'The file contains ' num2str(numPacketsInFile) ...
			'; the startEvent/startPacket parameter is ' num2str(startPacket)]);
end
if isfield(importParams, 'endPacket')
    endPacket = importParams.endPacket;
elseif isfield(importParams, 'endEvent')
    endPacket = importParams.endEvent;
else
    endPacket = numPacketsInFile;
end

if endPacket > numPacketsInFile
	disp([	'The file contains ' num2str(numPacketsInFile) ...
			'; the endEvent/endPacket parameter is ' num2str(endPacket) ...
			'; reducing the endPacket parameter accordingly.']);
        endPacket = numPacketsInFile;
end
if startPacket >= endPacket
	error([	'The startEvent/Packet parameter is ' num2str(startEvent) ...
		', but the endEvent/Packet parameter is ' num2str(endEvent) ]);
end

if isfield(importParams, 'startTime')
    startTime = importParams.startTime * 1e6 / 2^10;
else
    startTime = 0;
end

if isfield(importParams, 'endTime')
    endTime = importParams.endTime * 1e6 / 2^10;
else
    endTime = inf;
end

numPacketsToRead = endPacket - startPacket + 1;

%% Read data

disp('Reading data ...')
fseek(fileHandle, 0, 'bof');
allPackets = uint32(fread(fileHandle, numPacketsToRead, 'uint32', 0, 'b'));

%% Prepare to unpack data

% Let's just do this iteratively for now; not sure how to matricise this,
% though see below for a half-baked attempt

%% Prepare data masks

majorTimeStampMask        = bin2dec('0000 0000 0011 1111 1111 1111 1111 1111');
minorTimeStampMask        = bin2dec('0000 0000 0000 1111 1111 1100 0000 0000');
columnAddressMask         = bin2dec('0000 0000 0000 0000 0000 0011 1111 1111');
majorRowAddressMask       = bin2dec('0000 0000 0011 1111 0000 0000 0000 0000');
minorRowAddressAndPolMask = bin2dec('0000 0000 0000 0000 1111 1111 1111 1111');

%This following oneshot mask allows the identification of events to add
%from the one-shot encoding
minorRowAndPolOneShotMask = uint16(2 .^ (0:15)');
oneShotPolarity = [true(8, 1); false(8, 1)];
oneShotMinorRowAddress = uint16(repmat((0:7)', 2, 1));

% Loop variables
numEventsProcessed = 0;
currentLengthOfEventVectors = 1024;
majorTimeStamp = 0;
timeStampOffset = 0;
columnAddress = uint16(0);
minorTimeStamp = uint32(0);

% Data arrays
timeStamp	= zeros(currentLengthOfEventVectors, 1, 'uint32');
x			= zeros(currentLengthOfEventVectors, 1, 'uint16');
y			= zeros(currentLengthOfEventVectors, 1, 'uint16');
polarity    = false(currentLengthOfEventVectors, 1);

%% Unpack events by iterating through packets

for packetIndex = startPacket : endPacket
    if mod(packetIndex, 10000) == 0
		disp(['packet: ' num2str(packetIndex) '; ' ...
              'file position: ' num2str(floor(packetIndex / 2^17)) ' MB; ' ...
              '(' num2str(floor((packetIndex - startPacket + 1) / (endPacket - startPacket + 1) * 100)) '%)'])
	end

    currentPacket = allPackets(packetIndex);
    packetCode = bitshift(currentPacket, -24);
    if packetCode == 102 % 0x66
        % Timestamp packet
        currentTimeStamp = bitand(currentPacket, majorTimeStampMask);
        if majorTimeStamp == 0
            timeStampOffset = currentTimeStamp - 1;
        end
        majorTimeStamp = bitshift(currentTimeStamp - timeStampOffset, 10);
    else
        % apply startTime only to the majorTimeStamp - note that rezero is
        % automatic in this import, so actually comparing to the
        % timeStampOffset

        if startTime <= timeStampOffset
            if endTime < timeStampOffset
                break
            end
            if packetCode == 153 % 0x99
                % Column address packet
                minorTimeStamp = bitshift(bitand(currentPacket, majorTimeStampMask), -10);
                columnAddress = uint16(bitand(currentPacket, columnAddressMask));
            elseif packetCode == 204 % 0xCC
                % Events packet
                majorRowAddress       = uint16(bitshift(bitand(currentPacket, majorRowAddressMask), -13));
                minorRowAndPolOneShot = uint16(bitand(currentPacket, minorRowAddressAndPolMask));
                minorRowAndPolLogical = logical(bitand(minorRowAndPolOneShot, minorRowAndPolOneShotMask));
                numNewEvents          = nnz(minorRowAndPolLogical);
                polarityNew           = oneShotPolarity(minorRowAndPolLogical);
                minorRowAddress       = oneShotMinorRowAddress(minorRowAndPolLogical);
                numEventsProcessed = numEventsProcessed + numNewEvents;
                while numEventsProcessed > currentLengthOfEventVectors
                    timeStamp	= [timeStamp;	zeros(currentLengthOfEventVectors, 1, 'uint32')];
                    x			= [x;			zeros(currentLengthOfEventVectors, 1, 'uint16')];
                    y			= [y;			zeros(currentLengthOfEventVectors, 1, 'uint16')];
                    polarity    = [polarity;	false(currentLengthOfEventVectors, 1)];
                    currentLengthOfEventVectors = currentLengthOfEventVectors * 2;
                end
                rangeStart = numEventsProcessed - numNewEvents + 1;
                timeStamp(rangeStart : numEventsProcessed) = majorTimeStamp + minorTimeStamp;
                x        (rangeStart : numEventsProcessed) = columnAddress;
                y        (rangeStart : numEventsProcessed) = majorRowAddress + minorRowAddress;
                polarity (rangeStart : numEventsProcessed) = polarityNew;

            else
                error('packet formation error')
            end
        end
    end
end

%% Crop the data arrays and pack

aedat.info.deviceAddressSpace = [640 480];

if numEventsProcessed > 0
    polarityArray = polarity;
    polarity = [];
    keepLogical = [true(numEventsProcessed, 1); false(currentLengthOfEventVectors - numEventsProcessed, 1)];
    polarity.timeStamp	= timeStamp(keepLogical);
    polarity.y			= y(keepLogical) - 8; % These subtractions are because of the 1-based encoding in the SEC packet format
    polarity.x			= x(keepLogical) - 1; % These subtractions are because of the 1-based encoding in the SEC packet format
    polarity.polarity	= polarityArray(keepLogical);
    aedat.data.polarity = polarity;

    % Find first and last time stamps
    aedat = FindFirstAndLastTimeStamps(aedat);

    % Add NumEvents field for each data type
    aedat = NumEventsByType(aedat);
end

disp('Import finished')

%% Half-baked approach to matricising computation
%{

msb = bitget(allPackets, 32);
secondMsb = bitget(allPackets, 31);

eventsLogical = msb & secondMsb;
columnLogical = msb & ~secondMsb;
timeStampLogical = ~msb & secondMsb;

timeStamps = zeros(numPackets, 1, 'uint32');
timeStamps(timeStampLogical) = bitshift(allPackets(timeStampLogical), 10);
timeStamps = FloodFillColumnDownwards(timeStamps);

%}

author: panpanfei


text:
%% List al ECGtask availables
% Look for the available ECGtask classes defined in the ecg-kit_rootpath\common
% folder.
%
% Example
%
%   list_all_ECGtask()
%
% Arguments:
%
% Output:
%   All ECGtask classes found and its object handles.
%
% See also ECGwrapper
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Birthdate: 22/10/2014
% Last update: 22/10/2014
% Copyright 2008-2015
%
function [cKnownECGtasks, cKnownECGtasksHdl] = list_all_ECGtask()
common_path = fileparts(mfilename('fullpath'));
common_path = [common_path filesep ];

ECGt_filenames = dir([common_path 'ECGtask*.m' ]);
ECGt_filenames = {ECGt_filenames(:).name};

if( isempty(ECGt_filenames) )
    error('ECGwrapper:NoTasks', 'Could not find any ECGtask in %s', common_path );
else
    % avoid the abstract class
    ECGt_filenames = ECGt_filenames(~strcmpi(ECGt_filenames, 'ECGtask.m'));
    cKnownECGtasksHdl = cellfun(@(a)( eval(a(1:end-2)) ) , ECGt_filenames, 'UniformOutput', false);
    cKnownECGtasks = arrayfun(@(a)( a{1}.name ) , cKnownECGtasksHdl, 'UniformOutput', false);
end

if( nargout < 1 )
    fprintf(1, disp_option_enumeration( 'Task names found:\n', cKnownECGtasks));
    clear cKnownECGtasks
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
new_folder = [];
max_memo_GB = [];
% This code requires the file read_Intan_RHD2000_file.m which is provided
% on the Intan Technologies web site.
% new_folder: output folder with the parsed channels and metadata.
% If it doen't exist or is [], it will be 'wc_data'.
% max_memo_GB is an idea of the number of GB allocated for the data to be
% stored in RAM, so it is used to compute the number of segments in which
% the data should be split for processing
% Based in read_intan_data.m version 1.1, June 26, 2010
% (c) 2010, Intan Technologies, LLC
% For more information, see http://www.intantech.com
% For updates and latest version, see http://www.intantech.com/software.html

with_memory=true;
try
	memory;
catch
	with_memory=false;
end
if with_memory
	[userview,systemview] = memory;
	memo_avaible = floor(systemview.PhysicalMemory.Available*0.80);
	if exist('max_memo_GB','var') && ~isempty(max_memo_GB)
        max_memo = max_memo_GB*(1024)^3;
		if max_memo > memo_avaible
			error('max_memo_GB > 80% of Physical Memory Available')
		end
	else
		max_memo = memo_avaible;
	end
else
	max_memo = max_memo_GB*(1024)^3;
end

read_Intan_RHD2000_file;

sr = frequency_parameters.amplifier_sample_rate;
num_channels = length(amplifier_channels); % info sale de la header file

if ~exist('new_folder','var') || isempty(new_folder)
    new_folder = [pwd filesep 'wc_data'];
end

mkdir(new_folder)

% Read 'info.rhd'

% Read 'amplifier.dat'
fileinfo = dir([path filesep 'amplifier.dat']);
lts = fileinfo.bytes/(num_channels * 2); % int16 = 2 bytes
samples_per_channel = ceil(max_memo/(num_channels * 4)); %single (4 bytes)
num_segments = ceil(lts/samples_per_channel);
channels = 1:num_channels;
save([new_folder filesep 'intan_meta_data'],'sr','lts','channels')

fid = fopen([path filesep 'amplifier.dat'], 'r');

outfile_handles = cell(1,num_channels);

[fname] = 'ch';
for i = 1:num_channels
    outfile_handles{i} = fopen([new_folder filesep  fname '_' num2str(channels(i)) '.intch'],'w','l');
end

for j=1:num_segments
    ini = (j-1)*samples_per_channel;
    fin = min(j*samples_per_channel,lts);
    data2 = single(fread(fid,(fin-ini)*(num_channels),'int16')) * 0.195;
    for ind = 1:num_channels
        fwrite(outfile_handles{ind},data2(ind:num_channels:end),'single');
    end
    fprintf('Segment %d out of %d processed.\n',j,num_segments);
end
fclose('all');
clearvars
%Auxiliar para gui_plot
%time = num_samples / sr; %total time amplifier
%j=-1;
%save('trial','j')

author: csn-le


text:
function val = dualcursor2(state,deltalabelpos,marker_color,datalabelformatfcnh,axh)
%dualcursor2       Add dual vertical cursors to a plot
%
%Syntax
%   dualcursor2
%   dualcursor2('state');
%   dualcursor2([X1 X2]);
%   dualcursor2('state',deltalabelpos);
%   dualcursor2([X1 X2],deltalabelpos);
%   dualcursor2('update', ...)       %Updates existing cursors when data has changed
%   dualcursor2(..., ..., marker_color_spec);
%   dualcursor2(..., ..., ..., fcnhandle);
%   dualcursor2(..., ..., ..., ..., axh);
%   val = dualcursor2        %return values of the 2 cursors on the current axis
%   val = dualcursor2(h)     %return values of the 2 cursors on the axis or line w/handle h
%
%Description
% Easy Mode:
%  dualcursor2 on    %Turns on data cursors
%  dualcursor2 off   %Turns off data cursors
%  dualcursor2       %Toggles the state of data cursors
%  val = dualcursor2 %Return the coordinates of the 2 cursors
%      val = [x1 y1 x2 y2 slope avg stddev]
%
% Interaction:
%  Click on the cursors to drag them around.
%  Click on the cursor label to reposition.
%  If multiple lines are plotted, click on a line to make it active.
%  Right click on a cursor to:
%     - Export the selected region to the workspace
%       (exports as a structure named cursors)
%     - Export the selected region to a new figure
%     - Export cursor data to workspace
%       (exports as a variable named cursordata = [x1 y1 x2 y2 slope avg stddev])
%     - Remove the cursors from the plot
%  Click on the delta calculation labels to reposition them.
%
% Advanced Options:
% 1.  Specify initial x-coordinates for the cursors
%  dualcursor2([x1 x2]);   %Adds cursors at x1, x2
%
% 2.  Specify the location of the text label for displaying x2-x1, y2-y1 results
%  dualcursor2([],deltalabelpos);   %Specifies the location to display DeltaX
%                                         %and DeltaY calculations
%  deltalabelpos = [DeltaX_x, DeltaX_y; DeltaY_x, DeltaY_y] in normalized axis units
%  default = [.65 -.08;.9 -.08]  puts them just below the lower right-hand corner.
%  If this is too hard, you can also reposition the labels with your mouse!
%
% 3.  Specifying the color and marker for the cursors
%  dualcursor2([],[],marker_color_spec);
%      %Turns on data cursor, using specified marker and color
%      %marker_color_spec is a one or two element string, specifying color and/or marker
%      %style.  >>help plot   for valid marker and color specifiers
%  ex 'go'  - red circles
%     's'   - squares, with default color (red)
%
% 4.  Update the existing cursors when the data in the plot has changed.
%   dualcursor2('update');
%     This mode is useful when using the cursors on a live data stream.
%
%  This simply updates the appropriate values to reflect the latest data.  Here's
%  how to do this:
%    - When you initialize the graphics before taking data, make your initial call
%       to dualcursor2 (dualcursor2 on, for instance)
%    - In the code that updates the graphics display with new data, call
%   dualcursor2 update
%
% 5.  Specify custom formatting for the data label showing individual cursor values
%  This is implemented with a function handle (@).  You provide a handle to a
%  specifically formatted function:
%   dualcursor2(..., ..., ..., fcnhandle)
%  %Specifies the text and formatting for the data label
%  %fcnhandle is a handle to the function that defines the formatting.  This function
%  %must have the following argument syntax:
%       function textstring = mytextstring(xv,yv);
%       Input:
%         xv          (scalar)  The x cursor value
%         yv          (scalar)  The y cursor value
%       Output:
%         textstring  (string)  The formatted text string
%                     {'string'} Cell array of strings for multi-line display
%
%       Example
%       Create the following function.  save as frequencystring.m
%          function textstring = frequencystring(xv,yv);
%          textstring = {['Amp: ' num2str(yv,'%2g') ' dB']
%                        ['f: ' num2str(xv,'%2g') ' Hz'];};
%    %Use this code when you are ready to call dualcursor2
%    fcnhandle = @frequencystring;  %Create handle to your function
%    dualcursor2(..., ..., ..., fcnhandle)
%
% 6.  Add cursors to a specific axis (not necessarily current axis).
%   dualcursor2(..., ..., ..., ..., axh);
%    Adds the cursors to the axis with handle axh
%    This mode is useful when axis handles might be hidden
%
%
% Example
% %Set up some interesting data
%   load handel
%   Ns = 2^12;
%   Y = fft(y,Ns);
%   Y = 2/Ns*abs(Y(1:Ns/2));
%   df = Fs/Ns;
%   f = (0:1:Ns/2-1)*df;
%  %Plot it
%   figure;
%   plot(f,100*[Y sqrt(Y)])
%   title('My cursor example');
%   xlabel('Frequency (Hz)');
%   ylabel('Amplitude');
%   axis([0 1200 0 inf]);
%
% %Turn on cursors
%   dualcursor2
%
% %Turn off cursors
%   dualcursor2 off
%
% %Place cursors at x=300 and x=400, place DeltaX/DeltaY display in the upper left
% %  hand corner, and use a green square for the cursor.
%   dualcursor2([300 400],[.05 1.05; .25 1.05],'gs');
%
% %Move the cursors around. Try moving the data label, too
%
% %Now, wasn't that fun?  Finally, get the current cursor positions
%   val = dualcursor2;
%       %val = [x1 y1 x2 y2 slope avg stddev]
%
%
% NOTE: HandleVisibility of axis must be set to 'on' or 'callback'
%
% see also: DATALABEL, LINELABEL

% This function is provided as an example only.  It has not been
% tested, and therefore, it is not officially supported by The
% MathWorks, Inc.

% Written by Scott Hirsch
% shirsch@mathworks.com
% Copyright 2002-2009 The MathWorks, Inc.
% This is a (major) modfication of datalabel, available at MATLAB Central

% Updated with slope, average, and standard devation by Ari Feldman

%Parse input arguments
%If output argument, return the cursor values

%See if the user specified a formatting function for the datalabel

if nargin<4 || isempty(datalabelformatfcnh)    %Nope.  Use the default one included here.
    datalabelformatfcnh = @local_maketextstring;
end;

if nargout
    if nargin==0    %Use current axis
        h = gca;
    else
        h = state;
        if strcmp(get(h,'Type'),'line');
            h = get(h,'Parent');
        end;
    end;
    cursors = findobj(h,'Tag','Cursor');
    if length(cursors)==2       %Should be empty (no cursors), or length=2
        for ii=1:2
            cn = getappdata(cursors(ii),'CursorNumber');
            ind = (cn-1)*2+1:cn*2;      %Index into val
            val(ind) = getappdata(cursors(ii),'Coordinates');
        end;
            avgh = getappdata(gca,'Avgh_Handle');
            slopeh = getappdata(gca,'Slope_Handle');
            stddevh = getappdata(gca,'Stddev_Handle');
            avgstr=get(avgh,'String');
            slopestr=get(slopeh,'String');
            stddevstr=get(stddevh,'String');
            avgstr=regexp(avgstr,'\=','split');
            slopestr=regexp(slopestr,'\=','split');
            stddevstr=regexp(stddevstr,'\=','split');
            val(end+1)=str2num(slopestr{2});
            val(end+1)=str2num(avgstr{2});
            val(end+1)=str2num(stddevstr{2});
    else
        val = [];
        warning('I could not find any cursors');
    end;
    return
elseif nargin<5 || isempty(axh)    %Did the user specify a handle?
    axh = gca;
end;

%If no input arguments, switch state (turn off/on)
if (nargin==0 & nargout==0) | isempty(state)  %Switch state.  Check current state
    dots = findobj(axh,'Type','line','Tag','Cursor');  %See if there are any cursors
    if isempty(dots)    %None found.  Turn cursors on
        state = 'on';
    else
        state = 'off';
    end;
end;

%Check if the first argument is numeric.  The user is specifying
%  the initial x-coordinates of the markers
if nargin>=1 & isnumeric(state) %First input is x coordinates of markers
    x_init = state(:);

    %error check
    if length(x_init)~=2
        error('First input must be 2 element vector of x coordinates');
    end;
    state = 'on';       %Turn on data cursors.
else        %Default position = 1/3, 2/3 x axis limits
    %    xl = xlim;              %X Limits.  this is the letter L, not the number 1
    if strcmp(get(axh,'Type'),'figure') ||  strcmp(get(axh,'Type'),'uipanel') || strcmp(get(axh,'Type'),'root');    %user clicked on the axis itself; do nothing
        return
    end;

    xl = xlim(axh);
    lim = localObjbounds(axh);  % Problem using objbounds with hggroup objects, so I have a simple version of my own
    lim = lim(1:2);         %x values only
    xl(isinf(xl)) = lim(isinf(xl));
    width = diff(xl);     %Axis width
    x_init = xl(1)+[1/3 2/3]*width;

end;

% Get handle to figure just once
hFig = ancestor(axh,'figure');

switch state
    case 'on'
        %Initialization
        % Set the WindowButtonDownFcn
        % Add the cursors.

        %If there are already some data cursors on this plot, delete them!
        dualcursor2('off',[],[],[],axh);

        %Parse user inputs

        %Check for user input of position for delta labels
        if nargin<2 | isempty(deltalabelpos)
            deltalabelpos = [.65 -.08;.9 -.08;.65 -.12;.9 -.12;.65 -.16];            %Use defaults
            %[x1 y1; x2 y2]
        end;

        %Marker and color specification
        if nargin >= 3 ,
            %Parse marker string.  User might specify color, marker, or both
            colors = 'bgrcmyk';
            markers = '+o*.xsdv^><ph';
            for ii=1:length(marker_color)
                col_ind = strfind(colors,marker_color(ii));
                if ~isempty(col_ind)        %It's a color
                    color = marker_color(ii);
                else                        %Try a marker instead
                    mark_ind = strfind(markers,marker_color(ii));
                    if ~isempty(mark_ind)
                        marker = marker_color(ii);
                    end;
                end;
            end;
        end;

        %Handle default marker and color
        if ~exist('color','var'), color = 'r'; end; %set default
        if ~exist('marker','var'), marker = '*'; end; %set default

        %Add the cursors.
        %Ideally, the user specified the 2 x coordinates.
        %If not, just put the cursors somewhere nice.
        lineh = local_findlines(axh);

        %Line data
        if isempty(lineh)
            x1 = 0;
            y1 = 0;

            %             erasemode = 'normal';       %default
        else

            %if multiple lines, define callback to allow for selection

            set(lineh,'ButtonDownFcn', ...
                'setappdata(get(gco,''Parent''),''SelectedLine'',gco);dualcursor2(''selectline'',[],[],[],get(gco,''Parent''))');

            %set erasemode to xor.  This speeds things up a ton with large data sets
            %             erasemode = get(lineh,'EraseMode');
            %             set(lineh,'EraseMode','xor');

            lineh = min(lineh);           %Lets just use one line.  This was probably added first
            setappdata(axh,'SelectedLine',lineh); %The currently selected line.

            %Why the last line?  Because it is the first one added
            xl = get(lineh,'XData');
            yl = get(lineh,'YData');
            ynew=yl;
        end;

        %Find nearest value on the line
        [xv1,yv1] = local_nearest(x_init(1),xl,yl);
        [xv2,yv2] = local_nearest(x_init(2),xl,yl);

        %Build the string for the data label
        textstring1 = feval(datalabelformatfcnh,xv1,yv1);
        textstring2 = feval(datalabelformatfcnh,xv2,yv2);

        %Add the data label
        th1 = text(xv1,yv1,textstring1,'FontSize',8,'Tag','CursorText','Parent',axh);
        th2 = text(xv2,yv2,textstring2,'FontSize',8,'Tag','CursorText','Parent',axh);

        %For R13 or higher (MATLAB 6.5), use a background color on the text string
        v=ver('MATLAB');
        v=str2num(v.Version(1:3));
        if v>=6.5
            set(th1,'BackgroundColor','y');
            set(th2,'BackgroundColor','y');
        end;

        yl = ylim(axh);
        lim = localObjbounds(axh);
        lim = lim(3:4);     %y values only
        yl(isinf(yl)) = lim(isinf(yl));

        %Add the cursors
        ph1 = line([xv1 xv1 xv1],[yl(1) yv1 yl(2)], ...
            'Color',color, ...
            'Marker',marker, ...
            'Tag','Cursor', ...
            'UserData',[lineh th1], ...
            'LineStyle','-', ...
            'Parent',axh);
        ph2 = line([xv2 xv2 xv2],[yl(1) yv2 yl(2)], ...
            'Color',color, ...
            'Marker',marker, ...
            'Tag','Cursor', ...
            'UserData',[lineh th2], ...
            'LineStyle','-', ...
            'Parent',axh);

        %Add context menu to the cursors
%         cmenu = uicontextmenu('Parent',get(axh,'Parent'));
        cmenu = uicontextmenu('Parent',hFig);
        set([ph1 ph2],'UIContextMenu',cmenu);

        % Define the context menu items
        item1 = uimenu(cmenu, 'Label', 'Export region to workspace', ...
            'Callback', 'dualcursor2(''exportws'',[],[],[],get(gco,''Parent''))');
        item2 = uimenu(cmenu, 'Label', 'Export region to new figure', ...
            'Callback', 'dualcursor2(''exportfig'',[],[],[],get(gco,''Parent''))');
        item3 = uimenu(cmenu, 'Label', 'Export cursor data to workspace', ...
            'Callback', 'dualcursor2(''exportcursor'',[],[],[],get(gco,''Parent''));');
        item4 = uimenu(cmenu, 'Label', 'Turn cursors off', ...
            'Callback', 'dualcursor2(''off'',[],[],[],get(gco,''Parent''))', ...
            'Separator','on');

        setappdata(th1,'Coordinates',[xv1 yv1]);
        setappdata(ph1,'Coordinates',[xv1 yv1]);
        setappdata(th2,'Coordinates',[xv2 yv2]);
        setappdata(ph2,'Coordinates',[xv2 yv2]);
        setappdata(th1,'CursorNumber',1);
        setappdata(ph1,'CursorNumber',1);
        setappdata(th2,'CursorNumber',2);
        setappdata(ph2,'CursorNumber',2);
        setappdata(th1,'FormatFcnH',datalabelformatfcnh);
        setappdata(th2,'FormatFcnH',datalabelformatfcnh);
        setappdata(th1,'Offset',[0 0]);     %Offset for the text label from the data value
        setappdata(th2,'Offset',[0 0]);     %Offset for the text label from the data value

        %     mh = uicontextmenu('Tag','DeleteObject', ...
        %         'Callback','ud = get(gco,''UserData'');delete([gco ud(2)]);');
        %     set([th1 th2],'UIContextMenu',mh);

        set(th1,'UserData',[lineh ph1]);     %Store handle to line
        set(th2,'UserData',[lineh ph2]);     %Store handle to line

        %Calculate Difference
        dx = xv2 - xv1;
        dy = yv2 - yv1;
        xpos1=find(xl==xv1);
        xpos2=find(xl==xv2);
        avg=mean(ynew(xpos1:xpos2));
        slope=dy/dx;
        stddev=std2(ynew(xpos1:xpos1));

        if round2(avg,.01)==0
            acc={'%+.3e'};
        elseif avg/1000>1
            acc={'%+.3e'};
        else
            acc={'%+.3f'};
        end

        %Add display for cursor deltas
        deltah(1) = text(deltalabelpos(1,1),deltalabelpos(1,2),['({x}_{2}-{x}_{1}) = ' num2str(dx,acc{1})], ...
            'Units','Normalized', ...
            'HorizontalAlignment','left', ...
            'Tag','CursorDeltaText', ...
            'FontSize',8, ...
            'Interpreter','tex', ...
            'Parent',axh);      %dx
        deltah(2) = text(deltalabelpos(2,1),deltalabelpos(2,2),['{y}_{2}-{y}_{1} = ' num2str(dy,acc{1})], ...
            'Units','Normalized', ...
            'HorizontalAlignment','left', ...
            'Tag','CursorDeltaText', ...
            'FontSize',8, ...
            'Interpreter','tex', ...
            'Parent',axh);      %dy

        avgh=text(deltalabelpos(3,1),deltalabelpos(3,2),['avg= ' num2str(avg,acc{1})], ...
            'Units','Normalized', ...
            'HorizontalAlignment','left', ...
            'Tag','CursorDeltaText', ...
            'FontSize',8, ...
            'Interpreter','tex', ...
            'Parent',axh);      %dy
        slopeh=text(deltalabelpos(4,1),deltalabelpos(4,2),['slope= ' num2str(slope,acc{1})], ...
            'Units','Normalized', ...
            'HorizontalAlignment','left', ...
            'Tag','CursorDeltaText', ...
            'FontSize',8, ...
            'Interpreter','tex', ...
            'Parent',axh);      %dy
        stddevh=text(deltalabelpos(5,1),deltalabelpos(5,2),['stddev= ' num2str(stddev,'%+.3e')], ...
            'Units','Normalized', ...
            'HorizontalAlignment','left', ...
            'Tag','CursorDeltaText', ...
            'FontSize',8, ...
            'Interpreter','tex', ...
            'Parent',axh);      %dy

        %Set Application Data.
        setappdata(axh,'Delta_Handle',deltah);
        setappdata(axh,'Avgh_Handle',avgh);
        setappdata(axh,'Slope_Handle',slopeh);
        setappdata(axh,'Stddev_Handle',stddevh);
        setappdata(axh,'Marker',marker);
        setappdata(axh,'Color',color);

        set(hFig,'WindowButtonDownFcn','dualcursor2(''down'',[],[],[],get(gco,''Parent''))')
        set(hFig,'DoubleBuffer','on');       %eliminate flicker

    case 'down'      % Execute the WindowButtonDownFcn
        htype = get(gco,'Type');
        tag = get(gco,'Tag');
        marker = getappdata(axh,'Marker');
        color  = getappdata(axh,'Color');

        %If it's a movable object (Cursor, CursorText label, or Cursor Delta Text
        %label), make it movable.
        if strcmp(tag,'CursorText') | strcmp(tag,'Cursor') | strcmp(tag,'CursorDeltaText')
            set(gco,'EraseMode','xor')
            set(gcf,'WindowButtonMotionFcn','dualcursor2(''move'',[],[],[],get(gco,''Parent''))', ...
                'WindowButtonUpFcn','dualcursor2(''up'',[],[],[],get(gco,''Parent''))');
        end;

        if strcmp(tag,'Cursor') %If clicked on a cursor
            %Label the cursor we are moving.  Add a text label just above the plot
            CursorNumber = getappdata(gco,'CursorNumber');
            cnstr = num2str(CursorNumber);

            xl = xlim(axh);
            lim = localObjbounds(axh);
            lim = lim(1:2);     %x values only
            xl(isinf(xl)) = lim(isinf(xl));
            %             yl = lim(3:4);

            xv = getappdata(gco,'Coordinates');
            xv = xv(1);

            %Convert to normalized
            xn = (xv - xl(1))/(xl(2) - xl(1));
            yn = 1.05;

            CNh = text(xn,yn,cnstr, ...
                'Units','Normalized', ...
                'HorizontalAlignment','center', ...
                'Parent',axh);      %dx

            %For R13 or higher (MATLAB 6.5), use a background color on the text string
            v=ver('MATLAB');
            v=str2num(v.Version);
            if v>=6.5
                set(CNh,'BackgroundColor','y');
            end;
            setappdata(gco,'CNh',CNh);

        end

    case 'move'          % Execute the WindowButtonMotionFcn
        htype = get(gco,'Type');
        tag = get(gco,'Tag');
        if ~isempty(gco)

            %Special (simple) case - just repositioning the delta labels
            if strcmp(tag,'CursorDeltaText')    %The cursor delta labels.
                cp = get(axh,'CurrentPoint');
                pt = cp(1,[1 2]);

                %Put into normalized units
                ax = axis;
                lim = localObjbounds(axh);
                ax(isinf(ax)) = lim(isinf(ax));

                pt(1) = (pt(1) - ax(1))/(ax(2)-ax(1));
                pt(2) = (pt(2) - ax(3))/(ax(4)-ax(3));

                set(gco,'Position', [pt 0])
                drawnow
                return
            end;

            %Is this the cursor or the text
            if strcmp(tag,'CursorText')             %The text
                th = gco;
                handles = get(gco,'UserData');
                ph = handles(2);
                slide = 0;      %Don't slide along line; just reposition text
            else                                    %The marker
                ph = gco;
                handles = get(gco,'UserData');
                th = handles(2);
                slide = 1;      %Slide along line to next data point
            end;

            offset = getappdata(th,'Offset');       %Offset from data value

            cp = get(axh,'CurrentPoint');
            pt = cp(1,[1 2]);

            %Constrain to Line
            lh = getappdata(get(th,'Parent'),'SelectedLine');
            %             lh = handles(1);        %Line

            x = cp(1,1);       %first xy values
            y = cp(1,2);       %first xy values

            if slide            %Move to new data value
                xl = get(lh,'XData');
                yl = get(lh,'YData');
                ynew=yl;

                %Get nearest value
                [xv,yv]=local_nearest(x,xl,yl);

                %If we are moving a cursor, must move the cursor number label, too
                if strcmp(tag,'Cursor')
                    %Move the Cursor Number label, too
                    CNh = getappdata(gco,'CNh');
                    pos = get(CNh,'Position');
                    xlm = xlim(axh);
                    lim = localObjbounds(axh);
                    lim = lim(1:2);     %x values only
                    xlm(isinf(xlm)) = lim(isinf(xlm));
                    xn = (xv - xlm(1))/(xlm(2)-xlm(1));
                    set(CNh,'Position',[xn pos(2:3)])
                end;

                yl = ylim(axh);
                lim = localObjbounds(get(lh,'Parent'));
                lim = lim(3:4);     %y values only
                yl(isinf(yl)) = lim(isinf(yl));
                datalabelformatfcnh = getappdata(th,'FormatFcnH');

                textstring = feval(datalabelformatfcnh,xv,yv);
                set(th,'Position', [xv yv 0] + [offset 0],'String',textstring)
                set(ph,'XData',[xv xv xv],'YData',[yl(1) yv yl(2)]);

                setappdata(ph,'Coordinates',[xv yv]);
                setappdata(th,'Coordinates',[xv yv]);

                %Update delta calculation
                cursors = findobj(axh,'Tag','Cursor');
                cn1 = getappdata(cursors(1),'CursorNumber');
                if cn1==2   %Switch order
                    temp = cursors(1);
                    cursors(1) = cursors(2);
                    cursors(2) = temp;
                end;

                deltah = getappdata(axh,'Delta_Handle');    %Handle to cursors
                avgh = getappdata(axh,'Avgh_Handle');
                slopeh = getappdata(axh,'Slope_Handle');
                stddevh= getappdata(axh,'Stddev_Handle');

                %Positions of two dualcursor2s
                xy1 = getappdata(cursors(1),'Coordinates');
                xy2 = getappdata(cursors(2),'Coordinates');

                %Calculate Difference
                dx = xy2(1) - xy1(1);
                dy = xy2(2) - xy1(2);
                xpos1=find(xl==xy1(1));
                xpos2=find(xl==xy2(1));
                avg=mean(ynew(xpos1:xpos2));
                slope=dy/dx;
                stddev=std2(ynew(xpos1:xpos2));

                if round2(avg,.01)==0
                    acc={'%+.3e'};
                elseif avg/1000>1
                    acc={'%+.3e'};
                else
                    acc={'%+.3f'};
                end

                set(deltah(1),'String',['({x}_{2}-{x}_{1}) = ' num2str(dx,acc{1})]);
                set(deltah(2),'String',['({y}_{2}-{y}_{1}) = ' num2str(dy,acc{1})]);
                set(avgh,'String',['avg= ' num2str(avg,acc{1})]);
                set(slopeh,'String',['slope= ' num2str(slope,acc{1})]);
                set(stddevh,'String',['stddev= ' num2str(stddev,'%+.3e')]);

            else                %Just move text around.
                set(th,'Position', [x y 0])
            end;
            drawnow
        end;

    case 'up'           % Execute the WindowButtonUpFcn
        htype = get(gco,'Type');
        tag = get(gco,'Tag');
        if strcmp(tag,'CursorText') | strcmp(tag,'Cursor') | strcmp(tag,'CursorDeltaText');
            set(gco,'EraseMode','Normal')
            set(gcf,'WindowButtonMotionFcn','')

            if strcmp(tag,'CursorText') %If the text label, record it's relative position
                cp = get(axh,'CurrentPoint');
                pt = cp(1,[1 2]);

                coords = getappdata(gco,'Coordinates');
                offset(1) = pt(1) - coords(1);
                offset = pt - coords;
                setappdata(gco,'Offset',offset);
            end;

            if strcmp(tag,'Cursor')        %Delete the temporary cursor number label
                CNh = getappdata(gco,'CNh');
                delete(CNh);
            end;

        end;

    case 'selectline'          % User selected a new line to be active
        %Make the selected line bold
        lh = getappdata(axh,'SelectedLine');
        lw = get(lh,'LineWidth');
        set(lh,'LineWidth',5*lw);
        drawnow

        %Update the cursors
        dualcursor2('update',deltalabelpos,marker_color,datalabelformatfcnh,axh);

        %Put the line back the way you found it!
        set(lh,'LineWidth',lw);

    case 'update'               %Update the cursor value
        %Find the position of the existing cursors

        cursors = findobj(axh,'Tag','Cursor');
        cd1 = getappdata(cursors(1),'Coordinates');
        cd2 = getappdata(cursors(2),'Coordinates');
        cn1 = getappdata(cursors(1),'CursorNumber');
        cn2 = getappdata(cursors(2),'CursorNumber');
        clear x
        x(cn1) = cd1(1);        %x value of cursor number cn1
        x(cn2) = cd2(1);        %x value of cursor number cn2

        lh = getappdata(axh,'SelectedLine');

        handles = get(cursors(cn1),'UserData');
        th1 = handles(2);
        handles = get(cursors(cn2),'UserData');
        th2 = handles(2);

        offset1 = getappdata(th1,'Offset');       %Offset from data value
        offset2 = getappdata(th2,'Offset');       %Offset from data value

        ylm = ylim(axh);
        lim = localObjbounds(get(lh,'Parent'));
        lim = lim(3:4);     %y values only
        ylm(isinf(ylm)) = lim(isinf(ylm));

        xl = get(lh,'XData');
        yl = get(lh,'YData');
        ynew=yl;

        %Get nearest value
        [xv1,yv1]=local_nearest(x(1),xl,yl);
        [xv2,yv2]=local_nearest(x(2),xl,yl);

        datalabelformatfcnh = getappdata(th1,'FormatFcnH');

        %         textstring1 = {['x=' num2str(xv1)];['y=' num2str(yv1)]};
        %         textstring2 = {['x=' num2str(xv2)];['y=' num2str(yv2)]};
        textstring1 = feval(datalabelformatfcnh,xv1,yv1);
        textstring2 = feval(datalabelformatfcnh,xv2,yv2);
        set(th1,'Position', [xv1 yv1 0] + [offset1 0],'String',textstring1)
        set(th2,'Position', [xv2 yv2 0] + [offset2 0],'String',textstring2)
        set(cursors(cn1),'XData',[xv1 xv1 xv1],'YData',[ylm(1) yv1 ylm(2)]);
        set(cursors(cn2),'XData',[xv2 xv2 xv2],'YData',[ylm(1) yv2 ylm(2)]);
        %Update delta calculation
        deltah = getappdata(axh,'Delta_Handle');    %Handle to delta calculation
        avgh = getappdata(axh,'Avgh_Handle');
        slopeh = getappdata(axh,'Slope_Handle');
        stddevh = getappdata(axh,'Stddev_Handle');

        %                 %Positions of two dualcursor2s
        xy1 = getappdata(cursors(1),'Coordinates');
        xy2 = getappdata(cursors(2),'Coordinates');

        %Calculate Difference
        dx = xv2 - xv1;
        dy = yv2 - yv1;
        xpos1=find(xl==xv1);
        xpos2=find(xl==xv2);
        avg=mean(ynew(xpos1:xpos2));
        slope=dy/dx;
        stddev=std2(ynew(xpos1:xpos2));

        if round2(avg,.01)==0
            acc={'%+.3e'};
        elseif avg/1000>1
            acc={'%+.3e'};
        else
            acc={'%+.3f'};
        end

        set(deltah(1),'String',['{x}_{2}-{x}_{1} ' num2str(dx,acc{1})]);
        set(deltah(2),'String',['{y}_{2}-{y}_{1} ' num2str(dy,acc{1})]);
        set(avgh,'String',['avg= ' num2str(avg,acc{1})]);
        set(slopeh,'String',['Slope= ' num2str(slope,acc{1})]);
        set(stddevh,'String',['stddev= ' num2str(stddev,'%+.3e')]);

    case 'exportws'     %Export selected region to workspace
        [xd,yd] = local_extractregion(axh);

        %If there's only one line, don't make the user hassle with cell arrays
        if length(xd)==1
            xd = xd{1};
            yd = yd{1};
        end;

        cursors.xd = xd;
        cursors.yd = yd;

        assignin('base','cursors',cursors);
        disp('Variable: cursors created in workspace');

    case 'exportfig'     %Export selected region to a new figure
        %         [xd,yd,proplist,props] = local_extractregion;
        [xd,yd,hgS] = local_extractregion(axh);

        %Create new plot
        fh = figure;
        ax = axes;
        struct2handle(hgS,ax);

        %Get labels from axes, too
        %Clunky, but it seems to work
        %First, create empty title, xlabel, ylabel
        newhandles(1) = title('');
        newhandles(2) = xlabel('');
        newhandles(3) = ylabel('');

        %Get strings from original plot
        props = {'Title','xlabel','ylabel'};
        vals = get(axh,props);
        handles = [vals{:}];
        str = get(handles,'String');

        %Set the new strings to match
        set(newhandles,{'String'},str)

    case 'exportcursor'
        % Get cursor coordinates
        cursordata = dualcursor2(axh);

        % Push the data into the base workspace.
        assignin('base','cursordata',cursordata);
        disp('Variable: cursordata created in workspace');

    case 'off'   % Unset the WindowButton...Fcns
        % set(get(axh,'Parent'),'WindowButtonDownFcn','','WindowButtonUpFcn','')
        set(hFig,'WindowButtonDownFcn','','WindowButtonUpFcn','')

        h1 = findobj(axh,'Tag','CursorText');       %All text
        h2 = findobj(axh,'Tag','Cursor');           %The cursors
        h3 = findobj(axh,'Tag','CursorDeltaText');           %The cursors

        lineh = local_findlines(axh);
        set(lineh,'ButtonDownFcn','');

        %         erasemode = getappdata(axh,'OriginalEraseMode');
        %         if isempty(erasemode), erasemode = 'normal'; end;   %handles first time
        %         set(local_findlines(axh),'EraseMode',erasemode);

        delete([h1;h2;h3]);

end %switch/case on action
end

% function hFig = figHandle(axh)
% % Get the handle to the parent figure
% hFig = ancestor(axh,'figure');
% end

function [xv,yv]=local_nearest(x,xl,yl)
%Inputs:
% x   Selected x value
% xl  Line Data (x)

%Find nearest value of [xl] to (x)
%Special Case: Line has a single non-singleton value
if sum(isfinite(xl))==1
    fin = find(isfinite(xl));
    xv = xl(fin);
    yv = yl(fin);
else
    %Normalize axes
    xlmin = min(xl);
    xlmax = max(xl);
    xln = (xl - xlmin)./(xlmax - xlmin);
    xn = (x - xlmin)./(xlmax - xlmin);

    %Find nearest x value only.
    c = abs(xln - xn);

    [junk,ind] = min(c);

    %Nearest value on the line
    xv = xl(ind);
    yv = yl(ind);
end;
end

function textstring = local_maketextstring(xv,yv)
textstring = {['x = ' num2str(xv,'%2g')];
    ['y = ' num2str(yv,'%2.2g')]};
end

function [xd,yd,hgS] = local_extractregion(axh)
val = dualcursor2(axh);

%Find all lines
lineh = local_findlines(axh);

Nl = length(lineh);

%Get all line properties, so we can reproduce the appearance
hgS = handle2struct(lineh);       %Get properties

xd = get(lineh,'XData');
yd = get(lineh,'YData');

%Figure out the index into these lines that val corresponds to
x1ind = zeros(Nl,1);
x2ind = zeros(Nl,1);

%Special handling for single line case
if Nl==1
    xd= {xd};
    yd = {yd};
end;

for ii = 1:Nl
    x1ind(ii) = max(find(xd{ii}<=val(1)));
    x2ind(ii) = max(find(xd{ii}<=val(3)));

    %Keep data from this region only
    xd{ii} = xd{ii}(x1ind(ii):x2ind(ii));
    yd{ii} = yd{ii}(x1ind(ii):x2ind(ii));

    %Update handle structure, to make plotting really easy
    hgS(ii).properties.XData = xd{ii};
    hgS(ii).properties.YData = yd{ii};

end;
end

function lineh = local_findlines(axh);
lineh = findobj(axh,'Type','line');        %Find a line to add cursor to
dots = findobj(axh,'Type','line','Tag','Cursor');  %Ignore existing cursors
lineh = setdiff(lineh,dots);

%Ignore lines with only one or two values - these are probably annotations of some
%sort
xdtemp = get(lineh,'XData');
linehtemp = lineh;
lineh=[];
if ~iscell(xdtemp)      %If there's only one line, force data into a cell array
    xdtemp = {xdtemp};
end;

for ii=1:length(xdtemp);
    if length(xdtemp{ii})>2
        lineh = [lineh; linehtemp(ii)];
    end;
end;
end

function lim = localObjbounds(axh);
% Get x limits of all data in axes axh
kids = get(axh,'Children');
xmin = Inf; xmax = -Inf;
ymin = Inf; ymax = -Inf;
for ii=1:length(kids)
    try % Pass through if can't get data.  hopefully we hit at least one
        xd = get(kids(ii),'XData');
        xmin = min([xmin min(xd(:))]);
        xmax = max([xmax max(xd(:))]);

        yd = get(kids(ii),'YData');
        ymin = min([ymin min(yd(:))]);
        ymax = max([ymax max(yd(:))]);
    end
end
% Nuclear option, in case things went really bad
xmin(xmin==Inf) = 0;  xmax(xmax==-Inf) = 1;
ymin(ymin==Inf) = 0;  ymax(ymax==-Inf) = 1;

lim = [xmin xmax ymin ymax];
end

author: Sable