text:
% RANGE = showIm (MATRIX, RANGE, ZOOM, LABEL, NSHADES )
%
% Display a MatLab MATRIX as a grayscale image in the current figure,
% inside the current axes.  If MATRIX is complex, the real and imaginary
% parts are shown side-by-side, with the same grayscale mapping.
%
% If MATRIX is a string, it should be the name of a variable bound to a
% MATRIX in the base (global) environment.  This matrix is displayed as an
% image, with the title set to the string.
%
% RANGE (optional) is a 2-vector specifying the values that map to
% black and white, respectively.  Passing a value of 'auto' (default)
% sets RANGE=[min,max] (as in MatLab's imagesc).  'auto2' sets
% RANGE=[mean-2*stdev, mean+2*stdev].  'auto3' sets
% RANGE=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th percentile
% value of the sorted MATRIX samples, and p2 is the 90th percentile
% value.
%
% ZOOM specifies the number of matrix samples per screen pixel.  It
% will be rounded to an integer, or 1 divided by an integer.  A value
% of 'same' or 'auto' (default) causes the zoom value to be chosen
% automatically to fit the image into the current axes.  A value of
% 'full' fills the axis region (leaving no room for labels).  See
% pixelAxes.m.
%
% If LABEL (optional, default = 1, unless zoom='full') is non-zero, the range
% of values that are mapped into the gray colormap and the dimensions
% (size) of the matrix and zoom factor are printed below the image.  If label
% is a string, it is used as a title.
%
% NSHADES (optional) specifies the number of gray shades, and defaults
% to the size of the current colormap.

% Eero Simoncelli, 6/96.

%%TODO: should use "newplot"

function range = showIm( im, range, zoom, label, nshades );

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (nargin < 1)
  error('Requires at least one input argument.');
end

MLv = version;

if isstr(im)
  if (strcmp(MLv(1),'4'))
    error('Cannot pass string arg for MATRIX in MatLab version 4.x');
  end
  label = im;
  im = evalin('base',im);
end

if (exist('range') ~= 1)
  range = 'auto1';
end

if (exist('nshades') ~= 1)
  nshades = size(colormap,1);
end
nshades = max( nshades, 2 );

if (exist('zoom') ~= 1)
  zoom = 'auto';
end

if (exist('label') ~= 1)
  if strcmp(zoom,'full')
    label = 0;				% no labeling
  else
    label = 1;				% just print grayrange & dims
  end
end

%------------------------------------------------------------

%% Automatic range calculation:
if (strcmp(range,'auto1') | strcmp(range,'auto'))
  if isreal(im)
    [mn,mx] = range2(im);
  else
    [mn1,mx1] = range2(real(im));
    [mn2,mx2] =  range2(imag(im));
    mn = min(mn1,mn2);
    mx = max(mx1,mx2);
  end
  if any(size(im)==1)
    pad = (mx-mn)/12;			% MAGIC NUMBER: graph padding
    range = [mn-pad, mx+pad];
  else
    range = [mn,mx];
  end

elseif strcmp(range,'auto2')
  if isreal(im)
    stdev = sqrt(var2(im));
    av = mean2(im);
  else
    stdev = sqrt((var2(real(im)) + var2(imag(im)))/2);
    av = (mean2(real(im)) + mean2(imag(im)))/2;
  end
  range = [av-2*stdev,av+2*stdev]; 	% MAGIC NUMBER: 2 stdevs

elseif strcmp(range, 'auto3')
  percentile = 0.1;			% MAGIC NUMBER: 0<p<0.5
  [N,X] = histo(im);
  binsz = X(2)-X(1);
  N = N+1e-10;  % Ensure cumsum will be monotonic for call to interp1
  cumN = [0, cumsum(N)]/sum(N);
  cumX = [X(1)-binsz, X] + (binsz/2);
  ctrRange = interp1(cumN,cumX, [percentile, 1-percentile]);
  range = mean(ctrRange) + (ctrRange-mean(ctrRange))/(1-2*percentile);

elseif isstr(range)
  error(sprintf('Bad RANGE argument: %s',range))

end

if ((range(2) - range(1)) <= eps)
  range(1) = range(1) - 0.5;
  range(2) = range(2) + 0.5;
end

if isreal(im)
  factor=1;
else
  factor = 1+sqrt(-1);
end

xlbl_offset = 0; % default value

if (~any(size(im)==1))
  %% MatLab's "image" rounds when mapping to the colormap, so we compute
  %%      (im-r1)*(nshades-1)/(r2-r1) + 1.5
  mult = ((nshades-1) / (range(2)-range(1)));
  d_im = (mult * im) + factor*(1.5 - range(1)*mult);
end

if isreal(im)
  if (any(size(im)==1))
    hh = plot( im);
    axis([1, prod(size(im)), range]);
  else
    hh = image( d_im );
    axis('off');
    zoom = pixelAxes(size(d_im),zoom);
  end
else
  if (any(size(im)==1))
    subplot(2,1,1);
    hh = plot(real(im));
    axis([1, prod(size(im)), range]);
    subplot(2,1,2);
    hh = plot(imag(im));
    axis([1, prod(size(im)), range]);
  else
    subplot(1,2,1);
    hh = image(real(d_im));
    axis('off'); zoom = pixelAxes(size(d_im),zoom);
    ax = gca; orig_units = get(ax,'Units');
    set(ax,'Units','points');
    pos1 = get(ax,'Position');
    set(ax,'Units',orig_units);
    subplot(1,2,2);
    hh = image(imag(d_im));
    axis('off'); zoom = pixelAxes(size(d_im),zoom);
    ax = gca; orig_units = get(ax,'Units');
    set(ax,'Units','points');
    pos2 = get(ax,'Position');
    set(ax,'Units',orig_units);
    xlbl_offset = (pos1(1)-pos2(1))/2;
  end
end

if ~any(size(im)==1)
  colormap(gray(nshades));
end

if ((label ~= 0))
  if isstr(label)
    title(label);
    h = get(gca,'Title');
    orig_units = get(h,'Units');
    set(h,'Units','points');
    pos = get(h,'Position');
    pos(1:2) = pos(1:2) + [xlbl_offset, -3]; % MAGIC NUMBER: y pixel offset
    set(h,'Position',pos);
    set(h,'Units',orig_units);
  end

  if (~any(size(im)==1))
    if (zoom > 1)
      zformat = sprintf('* %d',round(zoom));
    else
      zformat = sprintf('/ %d',round(1/zoom));
    end
    if isreal(im)
      format=[' Range: [%.3g, %.3g] \n Dims: [%d, %d] ', zformat];
        else
      format=['Range: [%.3g, %.3g]  ----  Dims: [%d, %d]', zformat];
    end
    xlabel(sprintf(format, range(1), range(2), size(im,1), size(im,2)));
    h = get(gca,'Xlabel');
    set(h,'FontSize', 9); 		% MAGIC NUMBER: font size!!!

    orig_units = get(h,'Units');
    set(h,'Units','points');
    pos = get(h,'Position');
    pos(1:2) = pos(1:2) + [xlbl_offset, 10]; % MAGIC NUMBER: y offset in points
    set(h,'Position',pos);
    set(h,'Units',orig_units);

    set(h,'Visible','on');		% axis('image') turned the  xlabel  off...
  end
end

return;

repo: SelfExSR


text:
function mv = mv_sortByVoxR(mv,plotFlag);
%
% mv = mv_sortByVoxR([mv],[plotFlag]);
%
% For MultiVoxel UI, sorts voxels by the "voxel reliability"
% metric, which is computed in mv_reliability, and stored
% in the 'voxR' field of the anal struct that code produces.
%
% Images the voxels sorted by this metric (rows at the top
% will have a higher voxel reliability), and produces
% a performance curve by stepping through progressively
% larger subsets of the data, running the WTA classifier
% on those subsets, to estimate if this metric finds voxels
% that reliability discriminate objects.
%
% if plotFlag is set to 0, the code doesn't visualize the
% results, just returns the results in a 'voxRSorting'
% substruct.
%
% ras, 05/05
if ieNotDefined('mv')
    mv = get(gcf,'UserData');
end

if ieNotDefined('plotFlag')
    plotFlag = 1;
end

% params
sel = find(tc_selectedConds(mv)) - 1;
% sel = setdiff(sel-1, 13); % GUM FOR E-R ANALYSES
nConds = length(sel);
nRuns = length(unique(mv.trials.run));
names = mv.trials.condNames(sel); % ignore null

%%%%% recompute voxR for current event-related parameters
A1 = mv_amps(mv, 1:2:nRuns);
A2 = mv_amps(mv, 2:2:nRuns);
nVoxels = size(A1,1);
for v = 1:nVoxels
    [R p] = corrcoef(A1(v,sel), A2(v,sel));
    voxR(v) = R(1,2);
end

% grab unsorted voxel data from all runs
amps = mv_amps(mv);

% get new row index for voxel data
[ignore rank] = sort(voxR);
rank = fliplr(rank); % high rank at the top

% re-sort vox data by this ranking
sortedAmps = amps(rank,:);
A1 = A1(rank,:);
A2 = A2(rank,:);

% plot results if selected
if plotFlag==1
    oldaxes = findobj('Type','axes','Parent',gcf);
    delete(oldaxes)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % display a regression, color-coding each point
    % by its approximate rank
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(221);
    regressPlot(mv.wta.mu1,mv.wta.mu2,'x');
    hold on

    binSz = ceil(nVoxels/256); % # voxels per color bin
    ranked1 = repmat(NaN,[binSz 256]);
    ranked2 = repmat(NaN,[binSz 256]);
    ranked1(1:nVoxels) = mv.wta.mu1(rank);
    ranked2(1:nVoxels) = mv.wta.mu2(rank);
    plot(fliplr(ranked1),fliplr(ranked2),'x');
    xlabel('Mean Response, subset 1',...
        'FontName',mv.params.font,'FontSize',10);
    ylabel('Mean Response, subset 2',...
        'FontName',mv.params.font,'FontSize',10);
    setLineColors(jet(256));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % show sorted voxels
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(223);
    imagesc(sortedAmps'); % sorted
    set(gca,'YTick',1:nConds,'YTickLabel',names);
    title('Sorted')
    xlabel Voxels

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Test subplots of the sorted data and plot results
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(122);
    mv.voxRSorting = mv_wtaCurve(mv, A1, A2, 1);
    set(gca, 'ButtonDownFcn', 'mv_selectSubset([],[],''voxR'');');
    axis square
else
    mv.voxRSorting = mv_wtaCurve(mv, A1, A2, 0);
end

% add the sorted amps
mv.voxRSorting.sortedAmps = sortedAmps;
mv.voxRSorting.ampType = mv.params.ampType;
mv.voxRSorting.A1 = A1;
mv.voxRSorting.A2 = A2;
mv.voxRSorting.rank = rank;
mv.voxRSorting.metric = mv.wta.voxR;

% if a UI exists, set as user data
if isfield(mv.ui,'fig') & ishandle(mv.ui.fig)
    set(mv.ui.fig,'UserData',mv);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [patch_loc] = patchLocation(imagesize, size_patch, size_skip)

y = 1:size_skip(1):imagesize(1)-size_patch(1)+1;
x = 1:size_skip(2):imagesize(2)-size_patch(2)+1;

y = [y imagesize(1)-size_patch(1)+1];
x = [x imagesize(2)-size_patch(2)+1];

[Y,X] = meshgrid(y,x);
[dY,dX] = meshgrid(0:size_patch(1)-1,0:size_patch(2)-1);
dY = repmat(reshape(Y,[1 1 size(Y(:),1)]), [size_patch(1) size_patch(2) 1]) + repmat(dY, [1 1 size(Y(:),1)]);
dX = repmat(reshape(X,[1 1 size(X(:),1)]), [size_patch(1) size_patch(2) 1]) + repmat(dX, [1 1 size(X(:),1)]);
patch_loc = dY+(dX-1)*imagesize(1);

author: csjcai

repo: RealSR

save_path: github-repos/MATLAB/csjcai-RealSR


text:
%%  compares two groups of DTI images using tensor statistics and FDR.
%
% Probably originated with Armin Schwartzman and Bob.
% Might be deprecated or fixed up.
%
% Armin (c) Stanford VISTASOFT Team 2006???

% Load data
clear
dataPath = pwd;
dataSet = 'SIRL_DTI_sampleData';
[dataSet,dataPath] = uigetfile('*.mat','Select the sample data file.',fullfile(dataPath, dataSet));
if(isnumeric(dataSet)), error('Data file required.'); end
load(fullfile(dataPath, dataSet));
g1 = dti.groups{1};
g2 = dti.groups{2};

%% Test Statistics

% Log transformation
[vec,val] = dtiEig(dti.dt6);
logVal = log(val);
logDt6 = dtiEigComp(vec,logVal);

% Each of the following tests summarizes first the data, and then
% calls the test using the summaries only.
% This avoids having to load the data multiple times and saves memory.

% FA test
fa = dtiComputeFA(val);
[M1, S1, N1] = deal(mean(fa(:,g1),2), std(fa(:,g1),0,2), length(g1));
[M2, S2, N2] = deal(mean(fa(:,g2),2), std(fa(:,g2),0,2), length(g2));
[T, DISTR, df, M, S] = dtiTTest(M1, S1, N1, M2, S2, N2);

% Two-sample test of first eigenvector
vec1 = squeeze(vec(:,:,1,:));
[M1, S1, N1, Sbar1] = dtiDirMean(vec1(:,:,g1));
[M2, S2, N2, Sbar2] = dtiDirMean(vec1(:,:,g2));
[T, DISTR, df, M, S] = dtiDirTest(Sbar1, N1, Sbar2, N2);

% Two-sample test of frame of eigenvectors
[M1, S1, N1] = dtiLogTensorMean(logDt6(:,:,g1));
[M2, S2, N2] = dtiLogTensorMean(logDt6(:,:,g2));
[T, DISTR, df, M, S] = dtiLogTensorTest('vec', M1, S1, N1, M2, S2, N2);

% Display for any of the above tests
Timg = dtiIndToImg(T, dti.maskWM, NaN);
figure, imagesc(Timg(:,:,38)), axis image xy off, colormap('hot'), colorbar
pVal = -log10(1 - cdf(DISTR, T, df(1), df(2)));
pValImg = dtiIndToImg(pVal, dti.maskWM, NaN);
figure, imagesc(pValImg(:,:,38)), axis image xy off, colormap('hot'), colorbar

%--------------------------------------------------------------------
% FDR Analysis

% Simple FDR computation for given p-value threshold
pVal = 1 - cdf(DISTR, T, df(1), df(2));
thresh = 1e-4;
FDR = thresh ./ (sum(pVal < thresh) / length(T));

%--------------------------------------------------------------------
% FDR Analysis with empirical null

% Quantile transformation
switch DISTR,
case 't',
    T = norminv(cdf(DISTR, T, df(1), df(2)));
    DISTR = 'norm';
case 'f',
    T = chi2inv(cdf(DISTR, T, df(1), df(2)), df(1));
    DISTR = 'chi2';
end

% FDR analysis
dt = 0.2;
level = 0.05;
switch DISTR,
case 'norm',
    % Empirical null
    theoNull = fdrEmpNull(T, 'norm', dt, {});
    empNull = fdrEmpNull(Z, 'norm', dt, {'mu','s'});

    % FDR analysis (right tail)
    [fdrTNull, t] = fdrCurve(theoNull.fit, 'tail', 1); % use -1 for left tail
    [fdrENull, t] = fdrCurve(empNull.fit, 'tail', 1);
    thrTheoNull = fdrThresh(fdrTheoNull_R(:,1), t, level, 1);
    thrEmpNull = fdrThresh(fdrEmpNull_R(:,1), t, level, 1);

case 'chi2',
    % Empirical null
    theoNull = fdrEmpNull(T, 'chi2', dt, {}, df(1));
    empNull = fdrEmpNull(T, 'chi2', dt, {'a','nu'}, df(1));

    % FDR analysis (right tail)
    [fdrTNull, t] = fdrCurve(theoNull.fit, 'tail', 1);
    [fdrENull, t] = fdrCurve(empNull.fit, 'tail', 1);
    thrTheoNull = fdrThresh(fdrTNull(:,1), t, level, 1);
    thrEmpNull = fdrThresh(fdrENull(:,1), t, level, 1);
end

%-----------------------------------------------------------------------
% FDR Plots

% Histogram of test stats
figure, set(gcf, 'name', 'Histograms'), hold on
h = bar(theoNull.H.x, theoNull.H.hist, 1, 'w');
h0 = plot(theoNull.H0.x, theoNull.H0.hist, 'b');
h1 = plot(empNull.fit.x, empNull.fit.yhat, 'r');
hold off, legend([h0 h1], 'theo null','emp null',1)
xlabel('T'); ylabel('voxel count');
switch DISTR,
case 'norm',
    a=axis; axis([-4 4 a(3:4)]);
case 'chi2',
    a=axis; axis([0 prctile(T,99) a(3:4)]);
end

% FDR curves
figure,	set(gcf, 'name', 'FDR'), hold on
plot(t, fdrTNull(:,1), 'b')
plot(t, fdrENull(:,1), 'r')
legend('theo null','emp null')
plot(t, fdrTNull(:,2), 'b:', t, fdrTNull(:,3), 'b:')
plot(t, fdrENull(:,2), 'r:', t, fdrENull(:,3), 'r:')
hold off, axis([0 prctile(T,99.99) 0 1]);
xlabel('threshold'); ylabel('FDR');

%-----------------------------------------------------------------------
% Significant voxels

Timg = dtiIndToImg(T, dti.maskWM, NaN);
sgn = 1; % use -1 for left tail
fdrVol = (sgn * Timg > sgn * thrEmpNull);

figure, set(gcf, 'name', 'Significant Voxels')
bg = makeMontage(isfinite(Timg) & ~(fdrVol));
fg = makeMontage(fdrVol);
mont = cat(3, fg, 0.5*fg, bg);
image(mont); axis xy image off


text:
function prob_test148 ( )

%*****************************************************************************80
%
%% TEST148 tests UNIFORM_NSPHERE_SAMPLE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST148\n' );
  fprintf ( 1, '  For the Uniform PDF on the N-Sphere:\n' );
  fprintf ( 1, '  UNIFORM_NSPHERE_SAMPLE samples.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dimension N of sphere =       %6d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Points on the sphere:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ x, seed ] = uniform_nsphere_sample ( n, seed );
    fprintf ( 1, '  %6d  %14f  %14f  %14f\n', i, x(1:n) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = vl_nnnoffset(x, param, dzdy)
%VL_NNNOFFSET CNN norm-dependent offset.
%   Y = VL_NNNOFFSET(X, PARAM) subtracts from each element of X the
%   weighted norm of the feature channels:
%
%     X(i,j,k) = X(i,j,k) - PARAM(1) * L(i,j) ^ PARAM(2)
%
%   where
%
%     L(i,j) = sum_K X(i,j,k)^2
%
%   DZDX = VL_NNNOFFSET(X, PARAM, DZDY) computes the derivative of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y respectively.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

L = sum(x.^2,3) ;
L = max(L, single(1e-8)) ;
param = single(param) ;

if nargin <= 2
  y = bsxfun(@minus, x, param(1)*L.^param(2)) ;
else
  y = dzdy - bsxfun(@times, (2*param(1)*param(2))* x, sum(dzdy,3) .* (L.^(param(2)-1))) ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
% t_meshEdit
% Script (tutorial) to create, edit and visualize a mesh.
%
% See also T_MESHCREATE, MRGSAVECLASSWITHGRAY.
%
% Stanford VISTA
%

%% See t_meshCreate for a little more on creating meshes
dataD = mrvDataRootPath;
fName = fullfile(dataD,'anatomy','anatomyNIFTI','t1_class.nii.gz');

% Run the build code
msh = meshBuildFromClass(fName,[],'left');
msh = meshSmooth(msh);
msh = meshColor(msh);

% Visualize the coarse, unshaded mesh
meshVisualize(msh);

%% Set up parameters, smooth and visualize the mesh
msh = meshSet(msh,'smooth_relaxation',1);
msh = meshSet(msh,'smooth_sinc_method',0);
msh = meshSet(msh,'smooth_iterations',200);
msh2 = meshSmooth(msh);
meshVisualize(msh2);

%% Shade the mesh using local curvature
msh3 = meshColor(msh2,[],.25);

% Visualize the mesh smoothed and colored with the curvature
meshID = 2;
meshVisualize(msh3,meshID);

%% Now, visualize a mesh with some gray matter added to it
% We will add two gray layers
nLayers = 2;

% We use mrgGrowGray to read the WM and create new nodes/edges and white
% matter class data.
[nodes,edges,classData] = mrgGrowGray(fName,nLayers);

% mrgDisplayGrayMatter(nodes,edges,80,[120 140 120 140]);

% Add the gray matter to the white matter prior to creating the mesh.
wm = uint8( (classData.data == classData.type.white) ...
    | (classData.data == classData.type.gray) );

% We should add to meshBuildFromClass function an argument that specifies
% how many gray layers to add.  Then we could skip this step and just start
% below.

%% Now do the same sequence as above
mmPerVox = 0.7*[1,1,1];
msh4 = meshBuildFromClass(wm,mmPerVox);
msh4 = meshSmooth(msh4);
msh4 = meshColor(msh4);
meshVisualize(msh4,2);

%% END


text:
% Section 4.5.4: Frobenius norm diagonal scaling (GP)
% Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 01/29/06
% Updated to use GP mode by Almir Mutapcic 02/08/06
%
% Given a square matrix M, the goal is to find a vector (with dii > 0)
% such that ||DMD^{-1}||_F is minimized, where D = diag(d).
% The problem can be cast as an unconstrained geometric program:
%           minimize sqrt( sum_{i,j=1}^{n} Mij^2*di^2/dj^2 )
%

rs = randn( 'state' );
randn( 'state', 0 );

% matrix size (M is an n-by-n matrix)
n = 4;
M = randn(n,n);

% formulating the problem as a GP
cvx_begin gp
  variable d(n)
  minimize( sqrt( sum( sum( diag(d.^2)*(M.^2)*diag(d.^-2) ) ) ) )
  % Alternate formulation: norm( diag(d)*abs(M)*diag(1./d), 'fro' )
cvx_end

% displaying results
D = diag(d);
disp('The matrix D that minimizes ||DMD^{-1}||_F is: ');
disp(D);
disp('The minimium Frobenius norm achieved is: ');
disp(norm(D*M*inv(D),'fro'));
disp('while the Frobunius norm of the original matrix M is: ');
disp(norm(M,'fro'));

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [boxes, labels] = imBoundingBox(img, varargin)
% Bounding box of regions within a 2D or 3D binary or label image.
%
%   BOX = imBoundingBox(IMG)
%   Compute the bounding boxes of the particles in labeled image IMG. If
%   the image is binary, one box, corresponding to the foreground (i.e.
%   the pixels with value 1) will be computed.
%
%   The result is a N-by-4 array BOX = [XMIN XMAX YMIN YMAX], containing
%   coordinates of the box extent.
%
%   The same result could be obtained with the regionprops function. The
%   advantage of using imBoundingBox is that equivalent boxes can be
%   obtained in one call.
%
%   BOX = imBoundingBox(IMG3D)
%   If input image is a 3D array, the result is a N-by-6 array, containing
%   the maximal coordinates in the X, Y and Z directions:
%   BOX = [XMIN XMAX YMIN YMAX ZMIN ZMAX].
%
%   [BOX, LABELS] = imBoundingBox(...)
%   Also returns the labels of the regions for which a bounding box was
%   computed. LABELS is a N-by-1 array with as many rows as BOX.
%
%   [...] = imBoundingBox(IMG, LABELS)
%   Specifies the labels of the regions whose bounding box need to be
%   computed.
%
%
%   Example
%   % Draw a complex particle together with its bounding box
%     img = imread('circles.png');
%     imshow(img); hold on;
%     boxes = imBoundingBox(img);
%     drawBox(boxes)
%
%   % Compute and display the bounding box of several particles
%     img = imread('rice.png');
%     img2 = img - imopen(img, ones(30, 30));
%     lbl = bwlabel(img2 > 50, 4);
%     boxes = imBoundingBox(lbl);
%     imshow(img); hold on;
%     drawBox(boxes, 'linewidth', 2, 'color', 'g');
%
%   See also
%     regionprops, drawBox, imOrientedBox, imEquivalentEllipse

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2011-03-30,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% History
% 2013-03-29 add support for 3D images

%% Process input arguments

% default spatial calibration
nd = ndims(img);
spacing = ones(1, nd);
origin  = ones(1, nd);
calib   = false;

% extract spacing
if ~isempty(varargin) && all(size(varargin{1}) == [1 nd])
    spacing = varargin{1};
    varargin(1) = [];
    calib = true;
    origin = zeros(1, nd);

    % extract origin
    if ~isempty(varargin) && all(size(varargin{1}) == [1 nd])
        origin = varargin{1};
        varargin(1) = [];
    end
end

% check if labels are specified
labels = [];
if ~isempty(varargin) && size(varargin{1}, 2) == 1
    labels = varargin{1};
end

%% Initialisations

% extract the set of labels, without the background
if isempty(labels)
    labels = imFindLabels(img);
end

if nd == 2
    %% Process planar case
    props = regionprops(img, 'BoundingBox');
    bb = reshape([props.BoundingBox], [4 length(props)])';
    bb = bb(labels, :);

    % convert to MatImage convention
    boxes = [bb(:, 1) bb(:, 1)+bb(:, 3) bb(:, 2) bb(:, 2)+bb(:, 4)];

    % spatial calibration
    if calib
        boxes = bsxfun(@plus, bsxfun(@times, (boxes - 1), spacing([1 1 2 2])), origin([1 1 2 2]));
    end

elseif nd == 3
    %% Process 3D case
    stats = regionprops3(img, 'BoundingBox');
    bb = reshape([stats.BoundingBox], [6 size(stats, 1)])';
    bb = bb(labels, :);

    % convert to MatImage convention
    boxes = [bb(:, 1) bb(:, 1)+bb(:, 4) bb(:, 2) bb(:, 2)+bb(:, 5) bb(:, 3) bb(:, 3)+bb(:, 6)];

    % spatial calibration
    if calib
        boxes = bsxfun(@plus, bsxfun(@times, (boxes - 1), spacing([1 1 2 2 3 3])), origin([1 1 2 2 3 3]));
    end
else
    error('Image dimension must be 2 or 3');
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% Compute empirical p-vals under the null hypothesis that observed samples
% come from a given surrogate distribution. P-values for Type I error in
% rejecting the null hypothesis are obtained by finding the proportion of
% samples in the distribution that
% (a) are larger than the observed sample (one-sided test)
% (b) are larger or smaller than the observed sample (two-sided test).
%
% This function is based on Arnaud Delorme's statcond:COMPUTE_PVALS
% function from EEGLAB
%
% Inputs:
%
%   distribution:   [d1 x d2 x ... x dM x N] matrix of surrogate samples.
%                   distribution(i,j,k,...,:) is a collection of N samples
%                   from a surrogate distribution.
%   observed:       [d1 x d2 x ... x dM] matrix of observations.
%   tail:           can be 'one' or 'both' indicating a one-tailed or
%                   two-tailed test
% Outputs:
%
%   pvals:          [d1 x d2 x ... x dM] matrix of p-values specifying
%                   probability of Type I error in rejecting the null
%                   hypothesis
%
% Author: Tim Mullen and Arnaud Delorme, SCCN/INC/UCSD

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

function pvals = stat_surrogate_pvals(distribution,observed,tail)
numDims = ndims(distribution);
if iscolumn(distribution)
	numDims = 1;
end

n = size(distribution, numDims);
% once support for matlab <= R2016b is dropped:
% pvals = sum(distribution >= observed, numDims) / n;
pvals = sum(bsxfun(@ge, distribution, observed), numDims) / n;

if any(strcmpi(tail, {'right', 'one'}))
	% nothing to be done
	return;
end

p_left = 1 - pvals + sum(bsxfun(@eq, distribution, observed), numDims) / n;

if strcmpi(tail, 'both')
	pvals = 2 * min(pvals, p_left);
elseif strcmpi(tail, 'left')
	pvals = p_left;
else
	error('invalid value for tail: "%s", should be left, right, one or both', tail);
end

end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% --------------------------------------------------------------------------- %
% MarrRevisited - Surface Normal Estimation
% Copyright (c) 2016 Adobe Systems Incorporated and Carnegie Mellon University.
% All rights reserved.[see LICENSE for details]
% -------------------------------------------------------------------------- %

% Written by Aayush Bansal. Please contact ab.nsit@gmail.com
% demo code to use the surface normal mode --
clc; clear all;

%
conv_cache = ['./cachedir/demo_results/'];
if(~isdir(conv_cache))
        mkdir(conv_cache);
end

% initialize caffe
NET_FILE_PATH = ['./cachedir/surface_normal_models/'];
net_file     = [NET_FILE_PATH, 'best_model.caffemodel'];
DEPLOY_FILE_PATH = ['./net/conv/'];
deploy_file  = [DEPLOY_FILE_PATH, 'deploy.prototxt'];

% set the gpu --
% if not using GPU, set it to CPU mode.
gpu_id = 0;
caffe.reset_all;
caffe.set_device(gpu_id);
caffe.set_mode_gpu;
net = caffe.Net(deploy_file, net_file, 'test');

cnn_input_size = 224;
crop_height = 224; crop_width = 224;
image_mean = cat(3,  103.9390*ones(cnn_input_size),...
		     116.7700*ones(cnn_input_size),...
		     123.6800*ones(cnn_input_size));

% read the image set for NYU
img_data = {'demo/img_000001.jpg', 'demo/img_000002.jpg'};

% for each image in the img_set
for i = 1:length(img_data)

	display(['Image : ', img_data{i}]);
	ith_Img = im2uint8(imread(img_data{i}));

	%
        save_file_name = [conv_cache, strrep(img_data{i}, '.jpg', '')];
        if(exist([save_file_name, '.mat'], 'file'))
                continue;
        end

        j_ims = single(ith_Img(:,:,[3 2 1]));
        j_tmp = imresize(j_ims, [cnn_input_size, cnn_input_size], ...
                           'bilinear', 'antialiasing', false);
        j_tmp = j_tmp - image_mean;
        ims(:,:,:,1) = permute(j_tmp, [2 1 3]);

        snd(:,:,:,1) = (1/sqrt(3))*ones(cnn_input_size+200, cnn_input_size+200,3);
        depd(:,:,:,1) = zeros(cnn_input_size+200, cnn_input_size+200);

        %
        net.blobs('data0').reshape([crop_height+200, crop_width+200, 3, 1]);
        net.blobs('data2').reshape([crop_height+200, crop_width+200, 1, 1]);
        net.blobs('data1').reshape([crop_height+200, crop_width+200, 3, 1]);

        input_data = zeros(crop_height+200,crop_width+200,3,1);
        input_data(101:crop_width+100, 101:crop_width+100, :, 1) = ims;
        depd(101:crop_width+100, 101:crop_width+100, :, 1) = 1;
        net.blobs('data0').set_data(input_data);
        net.blobs('data2').set_data(depd);
        net.blobs('data1').set_data(snd);

        net.forward_prefilled();
        out = net.blobs('fc8_hcol').get_data();

        %%
        f2 = out';
        f2 = reshape(f2, [224, 224,3]);
        f2 = permute(f2, [2,1,3]);

        % normalize
        nx = f2(:,:,1); ny = f2(:,:,2); nz = f2(:,:,3);
        N = (nx.^2 + ny.^2 + nz.^2).^0.5 + eps;
        nx = nx./N; ny = ny./N; nz = nz./N;

        predns = cat(3, nx, ny, nz);
        predns = imresize(predns,...
                 [size(ith_Img,1), size(ith_Img,2)]);
        pred_N = (predns(:,:,1).^2 + predns(:,:,2).^2 +...
                                 predns(:,:,3).^2).^0.5 + eps;
        predns(:,:,1) = predns(:,:,1)./pred_N;
        predns(:,:,2) = predns(:,:,2)./pred_N;
        predns(:,:,3) = predns(:,:,3)./pred_N;

        predns_vis = uint8(255*(max(min(predns,1),-1)+1)/2);

        % dump the nx/ny/nz
        predns = single(predns);
        imwrite(predns_vis, [save_file_name, '.png']);
        save([save_file_name, '.mat'], 'predns')

end

% reset caffe
caffe.reset_all;

author: aayushbansal

repo: MarrRevisited


text:
function [d, dict, idict] = seg2dmat(im, outformat, res)
% SEG2DMAT  Local neighbourhood distance matrix between segmentation voxels
%
% D = SEG2DMAT(IM)
%
%   D is a sparse matrix where D(i,j) gives the distance between the i-th
%   and j-th voxels in the binary segmentation IM.
%
%   The index values i, j are computed with SUB2IND(), in the usual Matlab
%   way if you reshape IM into a vector, IM(:).
%
%   A 26-neighbourhood is assumed. That is, a voxel is only connected to
%   the 26 voxels that form a cube around it. That's why the sparse matrix
%   representation is convenient.
%
%   This function is fully vectorized and is thus very fast.
%
%   This function can be used with 2D images instead of 3D volumes too,
%   although some of the intermediate steps are not as efficient
%   memory-wise as they could be.
%
% [D, DICT, IDICT] = SEG2DMAT(IM, OUTFORMAT)
%
%   OUTFORMAT is a string. By default, D(i,j) gives the distance between
%   voxels i and j in IM. However, for large image volumes, it can be more
%   convenient to compute a smaller D where D(m, n) is the distance between
%   the m-th and n-th voxels of the segmentation:
%
%     'im' (default): Indices correspond to the whole image volume.
%
%     'seg': Indices correspond only to the segmentation (smaller matrix).
%
%   DICT and IDICT are column vectors used to convert to whole image and
%   segmentation indices. If OUTFORMAT='im', then i==j, so there's no need
%   for conversion and DICT and IDICT are returned empty.
%
%     DICT(i) is the matrix index for voxel i in the image
%     IDICT(i) is the image index for matrix index i.
%
% ... = SEG2DMAT(IM, OUTFORMAT, RES)
%
%    RES is a 3-vector with the voxel size given as [row, col, slice]. By
%    default, RES=[1 1 1].
%
% See also: im2imat.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2011 University of Oxford
% Version: 0.2.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(1, 3, nargin, 'struct'));
error(nargoutchk(0, 3, nargout, 'struct'));

% defaults
if (nargin < 2 || isempty(outformat))
    outformat = 'im';
end
if (nargin < 3 || isempty(res))
    res = [1 1 1];
end

% total number of voxels in the image
N = numel(im);

% image size
sz = size(im);
if (length(sz) == 2)
    sz(3) = 1;
end

% get linear indices of segmented voxels
idx0 = find(im);

% create dictionary vector that will allow us to convert between indices
% referred to the whole image, and indices of segmented voxels
%
% note that if you do dict(i) == 0 means that voxel i is not segmented
% dict = zeros(numel(im),1);
% dict(idx0) = (1:length(idx0))';
dict = sparse(idx0, ones(length(idx0), 1), (1:length(idx0)));

% convert to r, c, s indices
idx = idx0;
[r, c, s] = ind2sub(sz, idx);

% compute neighbourhood with connectivity 26 around origin
[gr, gc, gs] = ndgrid(-1:1, -1:1, -1:1);

% compute distances from each point to the origin
dlocal = sqrt((res(1)*gr).^2 + (res(2)*gc).^2 + (res(3)*gs).^2);

% convert volume of distances and coordinates into vectors
dlocal = dlocal(:);
gr = gr(:);
gc = gc(:);
gs = gs(:);

% each row has the coordinates of a neighbourhood around a segmented voxel
r = repmat(gr', length(idx), 1) + repmat(r, 1, 27);
c = repmat(gc', length(idx), 1) + repmat(c, 1, 27);
s = repmat(gs', length(idx), 1) + repmat(s, 1, 27);

% convert out of range subscripts into NaNs, so that they are not
% considered later when computing distances
r(r < 1 | r > sz(1)) = nan;
c(c < 1 | c > sz(2)) = nan;
s(s < 1 | s > sz(3)) = nan;

% reduce the r, c, s matrices to a single matrix with linear indices for
% the neighbours. Now, we have vector idx and matrix nn. For example:
% idx(5) = 34864
% nn(5,:) = [NaN NaN NaN NaN NaN NaN NaN NaN NaN 34395 34396 34397 ...
%            34863 34864 34865 35331 35332 35333 108807 108808 ...
%            108809 109275 109276 109277 109743 109744 109745]
% means that voxel 34864 has only 18 neighbours (voxel 34864 is in the
% first slice, so the neighbourhood is constrained that way)
nn = sub2ind(sz, r, c, s);

% we are going to remove the central column, because we don't need to
% connect a voxel to itself with distance 0
nn = nn(:, [1:13 15:end]);
dlocal = dlocal([1:13 15:end]);

% replicate local distance and index vectors so that it corresponds with
% the nn matrix
dlocal = repmat(dlocal', length(idx), 1);
idx = repmat(idx, 1, 26);

% find the out of range connections
ok = ~isnan(nn);

% remove them
idx = idx(ok);
nn = nn(ok);
dlocal = dlocal(ok);

% find indices larger than the largest one in the segmentation
ok = (nn <= max(idx0));

% remove them
idx = idx(ok);
nn = nn(ok);
dlocal = dlocal(ok);

% find connections to voxels that are not part of the segmentation
ok = (dict(nn) ~= 0);

% remove them
idx = idx(ok);
nn = nn(ok);
dlocal = dlocal(ok);

switch outformat
    case 'im'
        % create sparse matrix for distances between all voxels in the
        % image
        d = sparse(idx, nn, dlocal, N, N);

        % in this case, the i-th position in the matrix corresponds to the
        % i-th voxel in the image, so no dictionary is necessary
        dict = [];
        idict = [];
    case 'seg'
        % create sparse matrix for distances between all voxels in the
        % segmentation
        d = sparse(dict(idx), dict(nn), dlocal);

        % create dictionary to translate from image linear indices to
        % distance matrix indices. For example, distances to voxel im(3)
        % are dict(dict(3), :)
        dict = sparse(idx0, ones(length(idx0), 1), (1:length(idx0)));

        % compute inverse dictionary
        idict = find(dict);
    otherwise
        error('Unrecognized output format string')
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [v1,v2,v3]=vehicle()

% [v1,v2,v3]=vehicle;
% creates 3 different vehicle configurations, which vary one
% from another only by the distance between prow and wings.

% -----------------------------------------------------------------------
% configuration # 1

% phisical world variables
v1.g_e=[0; 0; 9.81];	% Gravitational acceleration
v1.rho=1033;			% Marine water density

% mass
v1.m=1444.0;			% mass (kg)

% Rigid body mass matrix [Kg Kg*m; Kg*m Kg*m^2]

v1.Mrb=[	1444	0	0	0	252.7	0	;
			0	1444	0	-252.7	0	-2527	;
			0	0	1444	0	2527	0	;
			0	-252.7	0	142.8	0	0	;
			252.7	0	2527	0	2796	0	;
			0	-2527	0	0	0	2778	];

% Added mass matrix [Kg Kg*m; Kg*m Kg*m^2]

v1.Ma=[	55.7	0	0	0	0	0	;
			0	1460	0	0	0	-102	;
			0	0	1460	0	102	0	;
			0	0	0	83.4	0	0	;
			0	0	102	0	3401	0	;
			0	-102	0	0	0	3401	];

% inverse total mass matrix
v1.iM=inv(v1.Mrb+v1.Ma);

% geometric variables

v1.l=3.5;			% Fuselage length (m)
v1.d=0.7;			% Fuselage diameter (m)
v1.vol=1.39787;		% Fuselage volume (m^3)

% vectors

v1.P_b=[0; 0; 0];		% Pole wrt B (m)
v1.G_b=[-1.75; 0; 0.175];	% Center of mass wrt B (m)
v1.B_b=[-1.75; 0; 0];	% Center of buoyancy wrt B (m)

% wings

v1.sw=0.2;			% Wing surface (m^2)
v1.cw=0.4;			% Wing width (m)
v1.bw=0.5;			% Wing length (m)

% Position of wings wrt B along x (m) in config 1
lw=-1.65;

% distance of wing middle point from prow and from axis in config 1
dpfw =    -lw +0.25*v1.cw;
dafw = v1.d/2 +0.43*v1.bw;

% wings middle point positions wrt B in config 1
v1.P1_b = [-dpfw; dafw;    0];
v1.P2_b = [-dpfw;    0; dafw];
v1.P3_b = [-dpfw;-dafw;    0];
v1.P4_b = [-dpfw;    0;-dafw];

% tails

v1.st=0.2;			% Tail surface (m^2)
v1.ct=0.4;			% Tail width (m)
v1.bt=0.5;			% Tail length (m)

lt=-3.1;				% Position of tails along x wrt B (m)

% distance of tail middle point from prow and from axis
dpft =    -lt +0.25*v1.ct;
daft = v1.d/2 +0.43*v1.bt;

% tails middle point positions wrt B
v1.P5_b = [-dpft; daft;    0];
v1.P6_b = [-dpft;    0; daft];
v1.P7_b = [-dpft;-daft;    0];
v1.P8_b = [-dpft;    0;-daft];

% -----------------------------------------------------------------------
% configuration # 2

v2=v1;

% Position of wings wrt B along x (m) in config 2
lw=-1.15;

% distance of wing middle point from prow and from axis in config 2
dpfw =    -lw +0.25*v2.cw;
dafw = v2.d/2 +0.43*v2.bw;

% wings middle point positions wrt B in config 2
v2.P1_b = [-dpfw; dafw;    0];
v2.P2_b = [-dpfw;    0; dafw];
v2.P3_b = [-dpfw;-dafw;    0];
v2.P4_b = [-dpfw;    0;-dafw];

% -----------------------------------------------------------------------
% configuration # 3

v3=v1;

% Position of wings wrt B along x (m) in config 3
lw=-0.35;

% distance of wing middle point from prow and from axis in config 3
dpfw =    -lw +0.25*v3.cw;
dafw = v3.d/2 +0.43*v3.bw;

% wings middle point positions wrt B in config 2
v3.P1_b = [-dpfw; dafw;    0];
v3.P2_b = [-dpfw;    0; dafw];
v3.P3_b = [-dpfw;-dafw;    0];
v3.P4_b = [-dpfw;    0;-dafw];

author: Sable

repo: mcbench-benchmarks


text:
function out = cfg_run_gunzip_files(job)
% Run gunzip on a list of files.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_gunzip_files.m 7499 2018-11-28 11:57:34Z guillaume $

rev = '$Rev: 7499 $'; %#ok

if isempty(job.outdir) || isempty(job.outdir{1})
    out = gunzip(job.files);
else
    out = gunzip(job.files, job.outdir{1});
end
if ~job.keep
    delete(job.files{:});
end

author: spm


text:
function [t, x, u] = snd_mpc(snd, output_data, varargin)
    warning off all;
%%
    snd.net_load = snd.load *2 - snd.PV/3 - snd.wind; %Power Balance

    %   Solve the optimal control problem
    t_Start = tic;
    [snd.u0, V_current, exitflag, output] = snd_solveOptimalControlProblem( snd );
    t_Elapsed = toc( t_Start );
    %V_current

    if ( snd.iprint >= 1 )
        printSolution(snd.printClosedloopData, snd.mpciter, snd.xmeasure, snd.u0, snd.iprint, exitflag, t_Elapsed, output_data);
    end

    %Save current information:
    x_estimated = snd_computeOpenloopSolution( snd , snd.u0 );

    t = output;
    x = x_estimated;
    u = snd.u0;

end

author: juchengquan


text:
function net = rbfsetfw(net, scale)
%RBFSETFW Set basis function widths of RBF.
%
%	Description
%	NET = RBFSETFW(NET, SCALE) sets the widths of the basis functions of
%	the RBF network NET. If Gaussian basis functions are used, then the
%	variances are set to the largest squared distance between centres if
%	SCALE is non-positive and SCALE times the mean distance of each
%	centre to its nearest neighbour if SCALE is positive.  Non-Gaussian
%	basis functions do not have a width.
%
%	See also
%	RBFTRAIN, RBFSETBF, GMMEM
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Set the variances to be the largest squared distance between centres
if strcmp(net.actfn, 'gaussian')
   cdist = dist2(net.c, net.c);
   if scale > 0.0
      % Set variance of basis to be scale times average
      % distance to nearest neighbour
      cdist = cdist + realmax*eye(net.nhidden);
      widths = scale*mean(min(cdist));
   else
      widths = max(max(cdist));
   end
   net.wi = widths * ones(size(net.wi));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function drawboard(board)

for i = 1:numel(board)
    patch(board(i).xy(:,1), board(i).xy(:,2), [0.5 0.5 0.5]);
end

axis image off


text:
report_this_filefun(mfilename('fullpath'));

l=storedcat(:,3) > mati+ 2.5/365 &  storedcat(:,3) <  mati+ 3/365 + 7/365;
b=storedcat.subset(l);

cd /home2/stefan/srisk/

load lat
load lon
lon = -lon;

le = length(lat);
Y0 = zeros(le,1);

for i = 1:length(b)
    di2 = deg2km((distance(lat,lon,repmat(b(i,2),le,1),repmat(b(i,1),le,1))));
    R = di2;
    r = sqrt(R.^2 + 5.57^2);
    M = b(i,6);
    Y = -0.136 + 0.229*(M-6) - 0.778 * log10(r) ;
    Y = 10.^Y;
    c = [Y , Y0];
    mapga = max(c');
    Y0 = mapga';
end

mapga = mapga';

cd /home2/stefan/srisk/

save hpga2 mapga -ascii

do = [' ! /home2/stefan/srisk/myriskobs2 ' num2str(max(mapga)/2.4,2)  '  ' num2str(max(mapga)/10,2) ]; eval(do)

% do = [' ! /home2/stefan/srisk/myriskobs 0.04 0.008' ]; eval(do)

author: CelsoReyes


text:
function [channels, skel] = acclaimLoadChannels(fileName, skel)

% ACCLAIMLOADCHANNELS Load the channels from an AMC file.
% FORMAT
% DESC loads channels from an acclaim motion capture file.
% ARG fileName : the file name to load in.
% ARG skel : the skeleton structure for the file.
% RETURN channels : the channels read in for the file.
% RETURN skel : the skeleton for the file.
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% SEEALSO : acclaimReadSkel

% MOCAP

fid=fopen(fileName, 'rt');

lin = getline(fid);
lin = strtrim(lin);
while ~strcmp(lin,':DEGREES')
  lin = getline(fid);
  lin = strtrim(lin);
end

counter = 0;
lin = getline(fid);

while lin~=-1
  lin = strtrim(lin);
  parts = tokenise(lin, ' ');
  if length(parts)==1
    frameNo = str2num(parts{1});
    if ~isempty(frameNo)
      counter = counter + 1;
      if counter ~= frameNo
        error('Unexpected frame number.');
      end
    else
      error('Single bone name  ...');
    end
  else
    ind = skelReverseLookup(skel, parts{1});
    for i = 1:length(skel.tree(ind).channels)
      bone{ind}{frameNo}(i) = str2num(parts{i+1});
    end
  end
  lin = getline(fid);
end
fclose(fid);
numFrames = counter;
numChannels = 0;
for i = 1:length(skel.tree)
  numChannels = numChannels + length(skel.tree(i).channels);
end
channels = zeros(numFrames, numChannels);

endVal = 0;
for i =1:length(skel.tree)
  if length(skel.tree(i).channels)>0
    startVal = endVal + 1;
    endVal = endVal + length(skel.tree(i).channels);
    channels(:, startVal:endVal)= reshape([bone{i}{:}], ...
                                          length(skel.tree(i).channels), numFrames)';

  end
  [skel.tree(i).rotInd, skel.tree(i).posInd] = ...
      resolveIndices(skel.tree(i).channels, startVal);
end

channels = smoothAngleChannels(channels, skel);

function [rotInd, posInd] = resolveIndices(channels, startVal)

% RESOLVEINDICES Get indices from the channels.

baseChannel = startVal - 1;
rotInd = zeros(1, 3);
posInd = zeros(1, 3);
for i = 1:length(channels)
  switch channels{i}
   case 'Xrotation'
    rotInd(1, 1) = baseChannel + i;
   case 'Yrotation'
    rotInd(1, 2) = baseChannel + i;
   case 'Zrotation'
    rotInd(1, 3) = baseChannel + i;
   case 'Xposition'
    posInd(1, 1) = baseChannel + i;
   case 'Yposition'
    posInd(1, 2) = baseChannel + i;
   case 'Zposition'
    posInd(1, 3) = baseChannel + i;
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [choices, feedbacks, simulation]=demo_QlearningSimulation()
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [choices, feedbacks] = demo_QlearningSimulation()
% Demo of Q-learning simulation. This demonstrate in particular how to use
% the toolbox to simulate data with a dynamic feedack (behaviour at times t
% can serve directly or indirectly as an input at time t+1)
%
% This is based on simple example of reinforcement learning algorithm.
% This demo first simulates 150 choices of a Q-learning agent when faced to
% a sequence binary alternatives:
% - at each trial, the agent choose one of two actions
% - the action is rewarded according to the probabilistic contingency rule
% - every 25 trials, the contingencies are reversed.
%
% /////////////////////////////////////////////////////////////////////////

f_fname = @f_Qlearning; % evolution function (Q-learning)
g_fname = @g_QLearning; % observation function (softmax mapping)

% Create the feedback rule for the simulation
% =========================================================================

% define which action should be rewarded at each trial (contingencies)
% -------------------------------------------------------------------------
% probability of a positive reward following a 'correct' action
    probRewardGood = 75/100;
% draw 25 random feedbacks
contBloc = +(rand(1,25) < probRewardGood);
% create 6 blocs with reversals
contingencies = [contBloc, 1-contBloc, ...
                 contBloc, 1-contBloc, ...
                 contBloc, 1-contBloc] ;

% create feedback structure for the simulation with VBA
% -------------------------------------------------------------------------
% feedback function. Return 1 if action follow contingencies.
h_feedback = @(yt,t,in) +(yt == contingencies(t));
% feedback structure for the VBA
fb = struct( ...
    'h_fname', h_feedback, ... % feedback function
    'indy', 1, ... % where to store simulated choice
    'indfb', 2, ... % where to store simulated feedback
    'inH', struct() ...
   );

% Simulate choices for the given feedback rule
% =========================================================================

% define parameteters of the simulated agent
% -------------------------------------------------------------------------
% learning rate
theta = VBA_sigmoid(0.65,'inverse',true); % 0.65, once sigm transformed
% inverse temperature
phi = log(2.5); % will be exp transformed
% initial state
x0 = [.5; .5];

% options for the simulation
% -------------------------------------------------------------------------
% number of trials
n_t = numel(contingencies);
% fitting binary data
options.sources.type = 1;
% Normally, the expected first observation (choice) is g(x1), ie. after
% a first iteratition x1 = f(x0). The skipf flag will prevent this evolution
% and thus set x1 = x0
options.skipf = [1 zeros(1,n_t)];

% simulate choices
% -------------------------------------------------------------------------
[y,x,x0,eta,e,u] = VBA_simulate ( ...
    n_t+1, ... number of trials
    f_fname, ... evolution function
    g_fname, ... observation function
    theta, ... evolution parameters (learning rate)
    phi, ... observation parameters,
    nan(2,n_t), ... dummy inputs
    Inf, Inf, ... deterministic evolution and observation
    options, ... options
    x0, ... initial state
    fb ... feedback rule
   );

% plot simulated choices
% -------------------------------------------------------------------------
hf = figure( ...
    'name', 'Simulated Q-learning behaviour', ...
    'color','w' ...
   );
ha = axes('parent',hf,'nextplot','add');
plot(ha,y,'kx')
plot(ha,y-e,'r')
legend(ha,{'y: agent''s choices','p(y=1|theta,phi,m): behavioural tendency'})

% Return simulated choices, feedbacks, and parameters used for the
% simulation
% =========================================================================
choices = u(1,2:end);
feedbacks = u(2,2:end);
simulation = struct( ...
    'state', x(:,1:end-1), ...
    'initial', x0, ...
    'evolution', theta, ...
    'observation', phi ...
    );

author: MBB-team

repo: VBA-toolbox


text:
function wfdb2mat(varargin)
%
% wfdm2mat(recordName,signaList,N,N0)
%
%    Wrapper to WFDB WFDB2MAT:
%         http://physionet.org/physiotools/wag/wfdb2m-1.htm
%
% Converts a WFDB-compatible signal file to MATLAB/Octave *.mat file.
% The output files are recordNamem.mat and recordNamem.hea. The standard output of
% WFDB2MAT will be saved in a file recordNamem.info.
%
% Required Parameters:
%
% recorName
%       String specifying the name of the record in the WFDB path or
%       in the current directory.
%
% Optional Parameters are:
%
% signalList
%       A Mx1 array of integers. Read only the signals (columns)
%       named in the signalList (default: read all signals).
% N
%       A 1x1 integer specifying the sample number at which to stop reading the
%       record file (default read all the samples = N).
% N0
%       A 1x1 integer specifying the sample number at which to start reading the
%       record file (default 1 = first sample).
%
%
%
% NOTE:
%       You can use the WFDB2MAT command in order to convert the record data into a *.mat file,
%       which can then be loaded into MATLAB/Octave's workspace using the LOAD command.
%       This will load the signal data in raw units (use RDMAT to load the signal in physical units).
%
%
% Written by Ikaro Silva, 2014
% Last Modified: September 15, 2014
% Version 0.1
%
% Since 0.9.7
%
% %Example:
% wfdb2mat('mitdb/200')
% [tm,signal,Fs,labels]=rdmat('200m');
%
%
% See also RDSAMP, RDMAT

%endOfHelp

[javaWfdbExec,config]=getWfdbClass('wfdb2mat');

%Set default pararameter values
inputs={'recordName','signalList','N','N0'};
signalList=[];
N=[];
N0=1;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

wfdb_argument={'-r',recordName,'-f',['s' num2str(N0-1)]};

if(~isempty(N))
    wfdb_argument{end+1}='-t';
    wfdb_argument{end+1}=['s' num2str(N)];
end

if(~isempty(signalList))
    wfdb_argument{end+1}='-s ';
    %-1 is necessary because WFDB is 0 based indexed.
    for sInd=1:length(signalList)
        wfdb_argument{end+1}=[num2str(signalList(sInd)-1)];
    end
end

data=javaWfdbExec.execToStringList(wfdb_argument);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function blas1_z_test16 ( )

%*****************************************************************************80
%
%% TEST16 tests ZSIGN2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST16\n' );
  fprintf ( 1, '  ZSIGN2 ( C1, C2 ) transfers the sign of complex C2\n' );
  fprintf ( 1, '  to the ZABS2 magnitude of C1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '           C1                    C2                    C3\n' );
  fprintf ( 1, ...
    '  --------------------  --------------------  --------------------\n' );
  fprintf ( 1, '\n' );

  for j = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );
    c1 = 5.0 * c1;
    [ c2, seed ] = c8_uniform_01 ( seed );
    c2 = 5.0 * c2;
    c3 = zsign2 ( c1, c2 );

    fprintf ( 1, '  %10f  %10f  %10f  %10f  %10f  %10f\n', ...
      real ( c1 ), imag ( c1 ), real ( c2 ), imag ( c2 ), ...
      real ( c3 ), imag ( c3 ) );

  end

  return
end


text:
function netOutputs = fromDagNN(dag, customFn)
%FROMDAGNN Converts a DagNN/SimpleNN object to the AutoNN framework
%   OUTPUTS = Layer.fromDagNN(DAG) converts a MatConvNet DagNN object, DAG,
%   into the AutoNN framework (i.e., a set of recursively nested Layer
%   objects). SimpleNN is also supported.
%
%   Returns a cell array of Layer objects, each corresponding to an output
%   of the network. These can be composed with other layers, or compiled
%   into a Net object for training/evaluation.
%
%   To convert a SimpleNN, first use DagNN.fromSimpleNN and then use this
%   function.
%
%   Example:
%     layers = Layer.fromDagNN(myDag) ;
%     net = Net(layers{:}) ;
%     net.eval({'images', randn(5,5,1,3), 'labels', 1:3}) ;
%
%   OUTPUTS = Layer.fromDagNN(DAG, @CUSTOMFN) adds support for custom DagNN
%   blocks, by calling function CUSTOMFN to convert unknown DagNN blocks.
%
%   The function's signature is customFn(block, inputs, params), where
%   block is the DagNN layer object, inputs is a cell array with the AutoNN
%   input layers, and likewise for params.
%
%   The following example would convert a DagNN block of class myClass,
%   with 1 input and 1 parameter, to a AutoNN layer that calls function
%   myfunc:
%     function obj = customFn(block, inputs, params)
%       if isa(block, 'myClass')
%         f = Layer.fromFunction(@myfunc) ;
%         obj = f(inputs{1}, params{1}, 'alpha', block.alpha) ;
%       end
%     end

% Copyright (C) 2016-2017 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  if nargin < 2
    customFn = [] ;
  end

  % handle struct inputs (SimpleNN or serialized DagNN)
  if isstruct(dag)
    assert(isfield(dag, 'layers'), 'Invalid struct (must be DagNN or SimpleNN).') ;
    if iscell(dag.layers)  % SimpleNN
      dag = dagnn.DagNN.fromSimpleNN(dag, 'CanonicalNames', true) ;
    else
      % serialized DagNN
      % compatibility: first remove 'exBackprop' property, which does not
      % exist in recent versions
      for l = 1:numel(dag.layers)
        if isfield(dag.layers(l).block, 'exBackprop')
          dag.layers(l).block = rmfield(dag.layers(l).block, 'exBackprop') ;
        end
      end
      % from struct to DagNN object
      dag = dagnn.DagNN.loadobj(dag) ;
    end
  end

  assert(isa(dag, 'dagnn.DagNN'), 'Input must be a DagNN or SimpleNN.') ;

  % update DagNN variable indexes
  dag.rebuild() ;

  % like initParams, but does not overwrite them if already defined
  initUninitializedParams(dag) ;

  assert(all(cellfun(@isscalar, {dag.layers.outputIndexes})), ...
    'Layers with 0 or multiple outputs are not supported yet.') ;

  % one-to-out relationship between each layer and its output var
  idx = [dag.layers.outputIndexes] ;
  var2layer = zeros(size(dag.vars)) ;
  var2layer(idx) = 1:numel(idx) ;

  % Layer handle corresponding to each DAG layer, empty if not processed
  layers = cell(1, numel(dag.layers)) ;

  % create all Param objects
  allParams = cell(size(dag.params)) ;
  for k = 1:numel(dag.params)
    p = dag.params(k) ;
    allParams{k} = Param('name', p.name, 'value', p.value, ...
      'learningRate', p.learningRate, 'weightDecay', p.weightDecay, ...
      'trainMethod', p.trainMethod) ;
  end

  % create Input objects for network inputs, corresponding to sources in
  % the graph
  netInputsIdx = find([dag.vars.fanin] == 0) ;
  netInputs = cell(1, numel(netInputsIdx));
  for k = 1:numel(netInputsIdx)
    netInputs{k} = Input('name', dag.vars(netInputsIdx(k)).name) ;
    var2layer(netInputsIdx(k)) = numel(layers) + k ;  % Inputs are at the end of the layers list
  end
  layers = [layers, netInputs] ;

  % find sinks in the graph, these are the network outputs
  netOutputsIdx = find([dag.vars.fanout] == 0) ;

  % recursively process the layers, starting with the network outputs
  netOutputs = cell(1, numel(netOutputsIdx)) ;
  for o = 1:numel(netOutputsIdx)
    netOutputs{o} = convertLayer(dag, var2layer(netOutputsIdx(o))) ;
  end

  % copy meta properties to one of the Layers
  netOutputs{1}.meta = dag.meta ;

  % process a single layer, and recurse on its inputs
  function obj = convertLayer(dag, layerIdx)
    % if this layer has already been processed, return its handle
    if ~isempty(layers{layerIdx})
      obj = layers{layerIdx} ;
      return
    end

    layer = dag.layers(layerIdx) ;
    params = allParams(layer.paramIndexes) ;
    block = layer.block ;

    % recurse on inputs; they must be defined before this Layer
    inputVars = layer.inputIndexes ;
    inputs = cell(size(inputVars)) ;
    for i = 1:numel(inputVars)
      inputs{i} = convertLayer(dag, var2layer(inputVars(i))) ;
    end

    % now create a Layer with those inputs and parameters

    switch class(block)
    case 'dagnn.Conv'
      if isscalar(params), params{2} = [] ; end  % no bias

      % ommit arguments with default values, for conciseness
      args = {} ;
      if ~all(block.stride == 1)
        args(end+1:end+2) = {'stride', block.stride} ;
      end
      if ~all(block.pad == 0)
        args(end+1:end+2) = {'pad', block.pad} ;
      end
      if ~all(block.dilate == 1)
        args(end+1:end+2) = {'dilate', block.dilate} ;
      end

      obj = vl_nnconv(inputs{1}, params{1}, params{2}, args{:}, block.opts{:}) ;

    case 'dagnn.ConvTranspose'
      if isscalar(params), params{2} = [] ; end  % no bias

      obj = vl_nnconvt(inputs{1}, params{1}, params{2}, ...
        'upsample', block.upsample, 'crop', block.crop, ...
        'numGroups', block.numGroups, block.opts{:}) ;

    case 'dagnn.BatchNorm'
      % make sure the Params are not empty, but scalar
      defaults = single([0, 0, 1]) ;
      for i = 1:3
        if isempty(params{i}.value)
          params{i}.value = defaults(i) ;
        end
      end

      obj = vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
        'moments', params{3}, 'epsilon', block.epsilon) ;

    case 'dagnn.Pooling'
      % ommit arguments with default values, for conciseness
      args = {} ;
      if ~all(block.stride == 1)
        args(end+1:end+2) = {'stride', block.stride} ;
      end
      if ~all(block.pad == 0)
        args(end+1:end+2) = {'pad', block.pad} ;
      end

      obj = vl_nnpool(inputs{1}, block.poolSize, ...
       'method', block.method, args{:}, block.opts{:}) ;

    case 'dagnn.ReLU'
      % ommit arguments with default values, for conciseness
      args = {} ;
      if block.leak ~= 0
        args(end+1:end+2) = {'leak', block.leak} ;
      end

      obj = vl_nnrelu(inputs{1}, args{:}, block.opts{:}) ;

    case 'dagnn.DropOut'
      obj = vl_nndropout(inputs{1}) ;

    case 'dagnn.Loss'
      obj = vl_nnloss(inputs{1}, inputs{2}, 'loss', block.loss, block.opts{:}) ;

    case 'dagnn.LRN'
      obj = vl_nnnormalize(inputs{1}, block.param) ;

    case 'dagnn.NormOffset'
      obj = vl_nnnoffset(inputs{1}, block.param) ;

    case 'dagnn.SpatialNorm'
      obj = vl_nnspnorm(inputs{1}, block.param) ;

    case 'dagnn.Sigmoid'
      obj = vl_nnsigmoid(inputs{1}) ;

    case 'dagnn.SoftMax'
      obj = vl_nnsoftmax(inputs{1}) ;
    case 'dagnn.GlobalPooling'
      obj = vl_nnglobalpool(inputs{1});
    case 'dagnn.Axpy'
      obj = vl_nnaxpy(inputs{1},inputs{2});
    case 'dagnn.Concat'
      obj = cat(block.dim, inputs{:}) ;

    case 'dagnn.Sum'
      obj = inputs{1} ;
      for i = 2:numel(inputs)
        obj = obj + inputs{i} ;
      end

    case 'dagnn.Crop'
      % replicate Crop block's functionality using simple operations
      A = inputs{1} ;
      B = inputs{2} ;

      sizeA1 = size(A,1);  % reuse this layer
      sizeA2 = size(A,2);
      cropv = sizeA1 - size(B,1) ;
      cropu = sizeA2 - size(B,2) ;
      cropv1 = max(0, cropv - block.crop(1)) ;
      cropu1 = max(0, cropu - block.crop(2)) ;

      crop = [cropv - cropv1, cropv1, cropu - cropu1, cropu1] ;
      obj = vl_nncrop(A, crop) ;

    otherwise
      if ~isempty(customFn)
        % handle custom blocks
        obj = customFn(block, inputs, params) ;
      else
        error(['Unknown block type ''' class(block) '''.']) ;
      end
    end

    obj.name = layer.name ;
    layers{layerIdx} = obj ;
  end
end

function initUninitializedParams(obj)
  % Same as INITPARAMS, but does not replace parameters that were already
  % initialized

  for l = 1:numel(obj.layers)
    p = obj.getParamIndex(obj.layers(l).params) ;
    params = obj.layers(l).block.initParams() ;
    switch obj.device
      case 'cpu'
        params = cellfun(@gather, params, 'UniformOutput', false) ;
      case 'gpu'
        params = cellfun(@gpuArray, params, 'UniformOutput', false) ;
    end

    for i = 1:numel(params)
      if isempty(obj.params(p(i)).value)
        obj.params(p(i)).value = params{i} ;
      end
    end
  end
end

repo: MFT


text:
function Z = and(X,Y)
%AND Logical AND (&) for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Z = tenfun(@and,X,Y);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function vw = addROIsinglePoint(vw,sgn)
%
% vw = addROIsinglePoint(vw,[sgn])
%
% AUTHOR:  Wandell
% DATE:   07.26.00
% PURPOSE:
%   The user clicks on a point and all of the gray matter within
% a distance, d, is selected for an ROI.  The calculation of the
% disk is managed by mrManifoldDistance.
%
% We find gray matter that is within some distance of the selected point.
% This is always done using gray matter graph in the VOLUME and coords data set
% via the routine mrManifoldDistance.  It appears that the proper nodes and
% edges are in Gray/coords.mat
%
% Hence, we have to
%   (a) make sure the gray graph is loaded,
%   (b) find the selected point in the gray matter,
%   (c) compute the disk around the point
%       return the corresponding values in whatever view we are in.
% If sgn~=0, the routine adds user-specified disk to the selected ROI in
% current slice. If sgn==0, the disk is removed from the current ROI.
%
% Follows logic from:  djh, 7/98
%
% If you change this function make parallel changes in:
%   all addROI*.m functions

% error if no current ROI.  This is usually created by the
% callback, so it should be here.
if vw.selectedROI == 0
  myErrorDlg('No current ROI');
  return
end

% Default is to add the data to the current ROI
if ~exist('sgn','var')
  sgn = 1;
end

% Get current ROI coords
curCoords = getCurROIcoords(vw);

% Save prevSelpts for undo
vw.prevCoords = curCoords;

% Get curSlice
curSlice = viewGet(vw, 'Current Slice');

% Get a single point from the user.  We leave the designation
% rgn, rather than pnt, for now.
% (rgn is short for region)
figure(vw.ui.figNum);
rgn = round(ginput(1));

% Note: ginput hands them back in x, y order (1st col is x and
% 2nd col is y).  But we use them in the opposite order (row,col), so that
% we want (y,x).  So we flip 'em. Hmmm....BW, but, but ...
%
rgn = fliplr(rgn);

% Check if outside image
dims=size(vw.ui.image);
if (min(rgn(:,1))< 1 | max(rgn(:,1))>dims(1) | ...
      min(rgn(:,2))< 1 | max(rgn(:,2))>dims(2))
  myWarnDlg('Selected point is outside image boundaries');
  return;
end

% Compute new coordinates
newCoords(1,:) = rgn(1);
newCoords(2,:) = rgn(2);
newCoords(3,:) = curSlice;

% Convert coords to canonical frame of reference
% Do an (inverse) rotation if necessary
if (strcmp(vw.viewType,'Flat'))
    newCoords=(rotateCoords(vw,newCoords,1));
end

% for VOLUME vw.  It does nothing for other views.
newCoords = curOri2CanOri(vw,newCoords);

% Merge/remove coordinates
if sgn
  coords = mergeCoords(curCoords,newCoords);
else
  coords = removeCoords(newCoords,curCoords);
end
vw.ROIs(vw.selectedROI).coords = coords;

return;

author: vistalab


text:
% Demo4: Fault detection and fault diagnosis for TE process using KPCA
% X: training samples
% Y: test samples

%%%%%%%%%%%%%%%%%%%%%%%%%%     KPCA       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Improve the performance by adjusting the following parameters:
% 1. sigma  % kernel width
% 2. pcr    % principal contribution rate
% ---------------------------------------------------------------------%

clc
clear all
close all
addpath(genpath(pwd))

% load TE process data
load train
load test

% Normalization
[X, Y] = normalize(train,test);

% Train KPCA model
model = kpca_train(X,'type',1,'sigma',800,'fd',1);

% Test a new sample Y (vector of matrix)
model = kpca_test(model,Y);

% Plot the result
plotResult(model.SPE_limit,model.SPE_test);
plotResult(model.T2_limit,model.T2_test);

% Fault diagnosis
[CPs_T2_test_s, CPs_SPE_test_s] = CPsKPCA(X,Y,model, ...
    'start_time',300,'end_time',500,'theta',0.7);

% Plot Contribution Plots
plotCPs(CPs_SPE_test_s)
plotCPs(CPs_T2_test_s)


text:
function dist=KLDiv(P,Q)
%  dist = KLDiv(P,Q) Kullback-Leibler divergence of two discrete probability
%  distributions
%  P and Q  are automatically normalised to have the sum of one on rows
% have the length of one at each
% P =  n x nbins
% Q =  1 x nbins or n x nbins(one to one)
% dist = n x 1

if size(P,2)~=size(Q,2)
    error('the number of columns in P and Q should be the same');
end

if sum(~isfinite(P(:))) + sum(~isfinite(Q(:)))
   error('the inputs contain non-finite values!')
end

% normalizing the P and Q
if size(Q,1)==1
    Q = Q ./sum(Q);
    P = P ./repmat(sum(P,2),[1 size(P,2)]);
    dist =  sum(P.*log(P./repmat(Q,[size(P,1) 1])),2);

elseif size(Q,1)==size(P,1)

    Q = Q ./repmat(sum(Q,2),[1 size(Q,2)]);
    P = P ./repmat(sum(P,2),[1 size(P,2)]);
    dist =  sum(P.*log(P./Q),2);
end

% resolving the case when P(i)==0
dist(isnan(dist))=0;


text:
function pbma_io_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests PBMA_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2008
%
%  Author:
%
%    John Burkardt
%
  file_name = 'pbma_io_test02.ascii.pbm';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  PBMA_READ reads an ASCII PBMA file.\n' );

  pbma_write_test ( file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PBMA_WRITE_TEST created data and wrote it to a file.\n' );
%
%  Now try to read the data.
%
  b = pbma_read ( file_name );

  [ nrow, ncol ] = size ( b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data file was read by PBMA_READ.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of rows of data =    %d\n', nrow );
  fprintf ( 1, '  Number of columns of data = %d\n', ncol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample data:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : 30
    i = floor ( ( ( 30 - k ) * 1 + ( k - 1 ) * nrow ) / ( 30 - 1 ) );
    j = floor ( ( ( 30 - k ) * 1 + ( k - 1 ) * ncol ) / ( 30 - 1 ) );
    fprintf ( 1, '%4d  %4d  %6d\n', i, j, b(i,j) );
  end

  return
end

author: johannesgerer


text:
function model = mogUpdatePrior(model)

% MOGUPDATEPRIORML Update the priors of an MOG model.
% FORMAT
% DESC updates the prior probabilities of a mixtures of
% Gaussians model.
% ARG model : the model which is to be updated.
% RETURN model : the model with updated priors.
%
% SEEALSO : mogCreate, mogUpdateMean, mogUpdateCovariance, mogEstep
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2008

% MLTOOLS

if model.isInfinite
  % First compute expectations of v.
  sumS = sum(model.posterior);
  a0bar = model.a0 + sumS; % Posterior value for a_0.
  a1bar = model.a1 + cumsum(sumS); % Posterior value for a_1.
  model.v = a0bar./(a0bar+a1bar);
  tmp = cumprod(1-model.v);
  model.prior = model.v;
  model.prior(2:end) = model.prior(2:end).*tmp(1:end-1);
else
  model.prior = mean(model.posterior);
  model.prior(find(model.prior==0))=1e-100;
end


text:
function sparse_display_test ( )

%*****************************************************************************80
%
%% SPARSE_DISPLAY_TEST tests the SPARSE_DISPLAY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_DISPLAY_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SPARSE_DISPLAY library.\n' );

  sparse_display_test01 ( );
  sparse_display_test02 ( );
  sparse_display_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_DISPLAY_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function M = getmassmat3(node,elem2dof,volume,type,K)
%% GETMASSMAT Get mass matrix of the finite element space
%
% M = GETMASSMAT(node,elem2dof,volume,type,K) get mass matrix of the finite element
% space specified by elemType.
%
% The type can be:
% - 'P1': full mass matrix for P1 element
% - 'lump': lumped mass matrix for P1 element

N = size(node,1);
NT = length(volume);
Ndof = double(max(elem2dof(:)));

%% Coefficients and default type
if ~exist('type','var'), type = 'P1'; end
if ~exist('K','var'), K = []; end

%% Assembling
n = size(elem2dof,2);
switch n
    case 4  % P1 element
    if strcmp(type,'lump')
        %% Assemble the mass matrix by the mass lumping
        M = accumarray([elem2dof(:,1);elem2dof(:,2);elem2dof(:,3);elem2dof(:,4)],...
                       [volume;volume;volume;volume]/4,[N,1]);
        if exist('K','var') && ~isempty(K) && ~isnumeric(K) % K is a function
            M = K(node).*M;
        elseif exist('K','var') && ~isempty(K) && isnumeric(K) && size(K,1) == N
            M = K.*M;
        end
        M = spdiags(M,0,N,N);
    else
        %% Assemble the full mass matrix
        if exist('K','var') && ~isempty(K) && ~isnumeric(K) % K is a function
            center = (node(elem2dof(:,1),:) + node(elem2dof(:,2),:) + ...
                      node(elem2dof(:,3),:) + node(elem2dof(:,4),:))/4;
            volume = K(center).*volume;
        elseif exist('K','var') && ~isempty(K) && isnumeric(K) && size(K,1) == NT
            volume = K.*volume;
        end
        M = sparse(N,N);
        for i = 1:4
            for j = i:4
                ii = double(elem2dof(:,i));
                jj = double(elem2dof(:,j));
                if (j==i)
                    M = M + sparse(ii,jj,volume/10,N,N);
                else
                    M = M + sparse([ii;jj],[jj;ii],[volume/20; volume/20],N,N);
                end
            end
        end
    end
    case 10 % P2 element
        %% Assemble the full mass matrix using nodal basis
        % indexing follows Poisson3P2
        [lambda, w] = quadpts3(2);
        nQuad = size(lambda,1);
        ii = zeros(55*NT,1); % 55: # of upper triangular entries
        jj = zeros(55*NT,1);
        sM = zeros(55*NT,1);
        phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
        phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
        phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
        phi(:,4) = lambda(:,4).*(2*lambda(:,4)-1);
        phi(:,5) = 4*lambda(:,1).*lambda(:,2);
        phi(:,6) = 4*lambda(:,1).*lambda(:,3);
        phi(:,7) = 4*lambda(:,1).*lambda(:,4);
        phi(:,8) = 4*lambda(:,2).*lambda(:,3);
        phi(:,9) = 4*lambda(:,2).*lambda(:,4);
        phi(:,10)= 4*lambda(:,3).*lambda(:,4);

        index = 0;
        for i = 1:10
            for j = i:10
                Mij = 0;
                for p = 1:nQuad; Mij = Mij + w(p)*phi(p,i).*phi(p,j); end
                Mij = Mij.*volume;
                ii(index+1:index+NT) = double(elem2dof(:,i));
                jj(index+1:index+NT) = double(elem2dof(:,j));
                sM(index+1:index+NT) = Mij;
                index = index + NT;
            end
        end
        diagIdx = (ii == jj);   upperIdx = ~diagIdx;
        M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),Ndof,Ndof);
        MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),Ndof,Ndof);
        M = M + MU + MU';

end

repo: ifem


text:
function [chansel, trlsel, cfg] = rejectvisual_trial(cfg, data)

% SUBFUNCTION for ft_rejectvisual

% determine the initial selection of trials
ntrl = length(data.trial);
if isequal(cfg.trials, 'all') % support specification like 'all'
  cfg.trials = 1:ntrl;
elseif isempty(cfg.trials)
  cfg.trials = 1:ntrl;
elseif isnumeric(cfg.trials)
  % use the selection as it is
end
trlsel = false(1, ntrl);
trlsel(cfg.trials) = true;

% determine the initial selection of channels
nchan = length(data.label);
cfg.channel = ft_channelselection(cfg.channel, data.label); % support specification like 'all'
chansel = false(1, nchan);
chansel(match_str(data.label, cfg.channel)) = true;

% compute the sampling frequency from the first two timepoints
fsample = 1/mean(diff(data.time{1}));

% compute the offset from the time axes
offset = zeros(ntrl,1);
for i=1:ntrl
  offset(i) = time2offset(data.time{i}, fsample);
end

if (isfield(cfg, 'preproc') && ~isempty(cfg.preproc))
  ft_progress('init', cfg.feedback, 'filtering data');
  for i=1:ntrl
    ft_progress(i/ntrl, 'filtering data in trial %d of %d\n', i, ntrl);
    [data.trial{i}, label, time, cfg.preproc] = preproc(data.trial{i}, data.label, data.time{i}, cfg.preproc);
  end
  ft_progress('close');
end

% select the specified latency window from the data
% this is done AFTER the filtering to prevent edge artifacts
if ischar(cfg.latency)
  cfg.latency(1) = min(cellfun(@min, data.time));
  cfg.latency(2) = max(cellfun(@max, data.time));
end

% the info structure will be attached to the figure
% and passed around between the callback functions
info.trlop   = 1;
info.chanlop = 1;
info.quit    = 0;
info.ntrl    = ntrl;
info.nchan   = nchan;
info.data    = data;
info.cfg     = cfg;
info.offset  = offset;
info.chansel = chansel;
info.trlsel  = trlsel;

tmpcfg = [];
tmpcfg.layout = 'ordered';
tmpcfg.skipscale = 'yes';
tmpcfg.skipcomnt = 'yes';
tmpcfg.channel = data.label;
info.layout = ft_prepare_layout(tmpcfg);

minx = min(info.layout.pos(:,1) - info.layout.width/2);
maxx = max(info.layout.pos(:,1) + info.layout.width/2);
miny = min(info.layout.pos(:,2) - info.layout.height/2);
maxy = max(info.layout.pos(:,2) + info.layout.height/2);

h = create_figure([minx maxx miny maxy], info);

% Maintain focus while GUI is open
while ishandle(h)
  drawnow
  info = guidata(h);
  if info.quit == 0
    uiwait;
  else
    chansel = info.chansel;
    trlsel = info.trlsel;
    delete(h);
    break
  end
end % while ishandle(h)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = create_figure(limits, info)
% Creates the GUI and plots the data for the first time.
% All additional changes are done by manipulating the properties of the
% specific graphical objects
%
% Input:
% limits - limits of the area where the trials are plotted
% info   - the struct that is attched to figure using guidata
h = open_figure(keepfields(info.cfg, {'figure', 'position', 'visible', 'renderer', 'figurename', 'title'}));
axis('off')
axis(limits);

% Set buttons and callback functions
info.ui.quit        = uicontrol(h,'units','pixels','position',[  5 5 40 18],'String','quit','Callback',@stop);
info.ui.prev        = uicontrol(h,'units','pixels','position',[ 50 5 25 18],'String','<','Callback',{@step_trial, -1});
info.ui.next        = uicontrol(h,'units','pixels','position',[ 75 5 25 18],'String','>','Callback',{@step_trial, 1});
info.ui.prev10      = uicontrol(h,'units','pixels','position',[105 5 25 18],'String','<<','Callback',{@step_trial, -10});
info.ui.next10      = uicontrol(h,'units','pixels','position',[130 5 25 18],'String','>>','Callback',{@step_trial, 10});
info.ui.exclude     = uicontrol(h,'units','pixels','position',[160 5 70 18],'String','exclude','Callback',@markexclude);
info.ui.include     = uicontrol(h,'units','pixels','position',[230 5 70 18],'String','include','Callback',@markinclude);
info.ui.excludenext = uicontrol(h,'units','pixels','position',[310 5 70 18],'String','exclude >','Callback',@markexclude_next);
info.ui.includenext = uicontrol(h,'units','pixels','position',[380 5 70 18],'String','include >','Callback',@markinclude_next);
set(h, 'WindowButtonUpFcn', @button);
set(h, 'KeyPressFcn', @key);

% Determine the y-axis limits for this channel for all included trials
[ymin, ymax] = get_ylim(info);

% Plot individual channels
hold('on')
color = 'k';
tim = info.data.time{info.trlop};
for chanindx=1:info.nchan
  dat = info.data.trial{info.trlop}(chanindx,:);
  info.h_chan(chanindx) = ft_plot_vector(tim, dat, 'hpos', info.layout.pos(chanindx,1), 'vpos', info.layout.pos(chanindx,2), 'width', info.layout.width(chanindx), 'height', info.layout.height(chanindx), 'vlim', [ymin ymax], 'box', istrue(info.cfg.box), 'color', color, 'label', info.layout.label{chanindx}, 'labelpos', 'lowerleft');
end

% Only blank excluded channels if the trial is included
if info.trlsel(info.trlop)
  set(info.h_chan(~info.chansel), 'Color', 'none');
else % If channel is excluded blank out everything
  set(info.h_chan, 'Color', 'none');
end
hold('off')

% Enable or disable buttons as appropriate
check_buttons_status(info);

% Add title
title(description_trial(info),'interpreter','none');

% Update the GUI
guidata(h, info);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function step_trial(h, event, step)
info = guidata(h);
if info.trlop == 1 && step < 0
  fprintf('At first trial\n');
elseif info.trlop == info.ntrl && step > 0
  fprintf('At last trial\n');
else
  trls = (1:info.ntrl);
  info.trlop = trls(nearest(trls, info.trlop + step));
  % Update plots
  new_plots(info);
end
guidata(h, info);
uiresume;
end

function markexclude_next(h, event)
markexclude(h, event);
step_trial(h, event, 1);
end

function markinclude_next(h, event)
markinclude(h, event);
step_trial(h, event, 1);
end

function markinclude(h, event)
info = guidata(h);
% Include trial if it's not already included
if ~info.trlsel(info.trlop)
  info.trlsel(info.trlop) = 1;
  fprintf(description_trial(info));
  title(description_trial(info));
  % Color in black all channels that are included (excluded remain blank)
  set(info.h_chan(info.chansel), 'Color', 'k');
  guidata(h,info);
end
uiresume;
end

function markexclude(h, event)
info = guidata(h);
% Exclude channel if it's not already excluded
if info.trlsel(info.trlop)
  info.trlsel(info.trlop) = 0;
  fprintf(description_trial(info));
  title(description_trial(info));
  % Blank all channels
  set(info.h_chan(info.chansel), 'Color', 'none');
  guidata(h,info);
end
uiresume;
end

function key(h, event)

switch lower(event.Key)
  case 'rightarrow'
    step_trial(h, event, 1);
  case 'leftarrow'
    step_trial(h, event, -1);
  case 'g'
    markinclude(h, event);
  case 'b'
    markexclude(h, event);
  case 'q'
    stop(h);
  otherwise
    fprintf('unknown key pressed\n');
end
end

function button(h, event)
% Find selected trial
pos = get(gca, 'CurrentPoint');
x = pos(1,1);
y = pos(1,2);
info = guidata(h);
dx = info.layout.pos(:,1) - x;
dy = info.layout.pos(:,2) - y;
dd = sqrt(dx.^2 + dy.^2);
[d, i] = min(dd);
threshold = median(sqrt(info.layout.width.^2 + info.layout.height.^2)/2);
if d<threshold && i<=info.nchan
  % Toggle
  info.chansel(i) = ~info.chansel(i); % toggle
  info.chanlop = i;
  guidata(h,info);

  % Print to command window
  fprintf(description_channel(info));

  % Update plots
  update_plots(info);

  % Update trial color
  if strcmp(info.h_chan(i).Color, 'none') % channel is toggled to include
    info.h_chan(i).Color = 'k';
  else % channel is toggled to exclude
    info.h_chan(i).Color = 'none';
  end
else
  fprintf('button clicked\n');
end
uiresume;
end

function stop(h, event)
info = guidata(h);
info.quit = 1;
guidata(h,info);
uiresume;
end

function str = description_channel(info)
if info.chansel(info.chanlop)
  str = sprintf('channel %s marked to INCLUDE\n', info.data.label{info.chanlop});
else
  str = sprintf('channel %s marked to EXCLUDE\n', info.data.label{info.chanlop});
end
end

function str = description_trial(info)
if info.trlsel(info.trlop)
  str = sprintf('trial %d marked to INCLUDE\n', info.trlop);
else
  str = sprintf('trial %d marked to EXCLUDE\n', info.trlop);
end
end

function [ymin, ymax] = get_ylim(info)
% Determines the y-axis limits for one trial for all included channels

if ~isnumeric(info.cfg.ylim)
  ymin = nan(1,info.nchan);
  ymax = nan(1,info.nchan);
  for chanindx = find(info.chansel)
    dat = info.data.trial{info.trlop}(chanindx,:);
    ymin(chanindx) = min(dat);
    ymax(chanindx) = max(dat);
  end
  ymin = min(ymin);
  ymax = max(ymax);
  if strcmp(info.cfg.ylim, 'maxabs') % handle maxabs, make y-axis center on 0
    ymax = max(abs(ymax), abs(ymin));
    ymin = -ymax;
  elseif strcmp(info.cfg.ylim, 'zeromax')
    ymin = 0;
  elseif strcmp(info.cfg.ylim, 'minzero')
    ymax = 0;
  end
else
  ymin = info.cfg.ylim(1);
  ymax = info.cfg.ylim(2);
end
end

function update_plots(info)
% Updates the trial plots after inclusion/exclusion of trials
%
% Input:
% info     - info struct with all the data attched to the GUI

[ymin, ymax] = get_ylim(info);

% Get positions and heights
v_centers = info.layout.pos(:, 2);
heights   = info.layout.height;

for chanindx = find(info.chansel)
  % Shift the vertical axis to zero
  vdat = info.data.trial{info.trlop}(chanindx, :);
  % Scale to length 1 of the new data
  vdat = vdat ./ (ymax - ymin);
  % Scale by availabe plot area
  vdat = vdat .* heights(chanindx);
  % Shift to the vertical position
  vdat = vdat + v_centers(chanindx);
  % Update
  info.h_chan(chanindx).YData = vdat;
end
end

function new_plots(info)
% Updates the trial plots after channel change
%
% Input:
% info - info struct with all the data attched to the GUI

[ymin, ymax] = get_ylim(info);

% Get positions and heights
v_centers = info.layout.pos(:, 2);
h_centers = info.layout.pos(:, 1);
heights   = info.layout.height;
widths    = info.layout.width;

% Update data of ALL channels, also excluded channels
xmin = min(info.data.time{info.trlop});
xmax = max(info.data.time{info.trlop});
for chanindx = 1:info.nchan
  % Shift the horizontal axis to zero
  hdat = info.data.time{info.trlop} - (xmin + xmax)/2;
  % Scale to length 1 of the new data
  hdat = hdat ./ (xmax - xmin);
  % Scale by availabe plot area
  hdat = hdat .* widths(chanindx);
  % Shift to the horizontal position
  hdat = hdat + h_centers(chanindx);

  % Shift the vertical axis to zero
  vdat = info.data.trial{info.trlop}(chanindx, :);
  % Scale to length 1 of the new data
  vdat = vdat ./ (ymax - ymin);
  % Scale by availabe plot area
  vdat = vdat .* heights(chanindx);
  % Shift to the vertical position
  vdat = vdat + v_centers(chanindx);
  % Update
  set(info.h_chan(chanindx), {'XData', 'YData'}, {hdat, vdat});
end

% Blank excluded channels
if info.trlsel(info.trlop) % channel is included
  set(info.h_chan(info.chansel), 'Color', 'k');
  set(info.h_chan(~info.chansel), 'Color', 'none');
else % if channel is excluded blank out everything
  set(info.h_chan, 'Color', 'none');
end

% Enable or disable buttons as appropriate
check_buttons_status(info);

% Update title
title(description_trial(info),'interpreter','none');
end

function check_buttons_status(info)
% Updates buttons status if required
%
% Input:
% info - info struct with all the data attched to the GUI

if info.trlop == 1
  info.ui.prev.Enable   = 'off';
  info.ui.prev10.Enable = 'off';
else
  info.ui.prev.Enable   = 'on';
  info.ui.prev10.Enable = 'on';
end
if info.trlop == info.ntrl
  info.ui.next.Enable   = 'off';
  info.ui.next10.Enable = 'off';
else
  info.ui.next.Enable   = 'on';
  info.ui.next10.Enable = 'on';
end
if info.trlop == info.ntrl
  info.ui.excludenext.Enable = 'off';
  info.ui.includenext.Enable = 'off';
else
  info.ui.excludenext.Enable = 'on';
  info.ui.includenext.Enable = 'on';
end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function  [W, U, mu, UtU] = update_params(mu, W, U, dWUtot, nspikes)

[Nchan, Nfilt, Nrank] = size(U);

dWUtotCPU = gather_try(dWUtot);
ntot = sum(nspikes,2);

for k = 1:Nfilt
    if ntot(k)>5

        [Uall, Sv, Vall] = svd(gather_try(dWUtotCPU(:,:,k)), 0);
        Sv = diag(Sv);
        sumSv2 = sum(Sv(1:Nrank).^2).^.5;
        for irank = 1:Nrank
            [~, imax] = max(abs(Uall(:,irank)), [], 1);
            W(:,k,irank) = - Uall(:,irank) * sign(Uall(imax,irank)) * Sv(irank)/sumSv2;
            U(:,k,irank) = - Vall(:,irank) * sign(Uall(imax,irank));
        end
        mmax = max(abs(U(:,k,1)));
        Usize = squeeze(abs(U(:,k,:)));
        Usize = Usize .* repmat(Sv(1:Nrank)'/Sv(1), Nchan, 1);
        ibad = max(Usize, [], 2) < .1 * mmax;

        U(ibad,k,:) = 0;
    end
end

% mu = zeros(Nfilt,1, 'single');
for k = 1:Nfilt
    if ntot(k)>5
        wu = squeeze(W(:,k,:)) * squeeze(U(:,k,:))';
        mu(k) = sum(sum(wu.*squeeze(dWUtotCPU(:,:,k))));
    end
end

for k = 1:Nfilt
    if ntot(k)>5
        wu = squeeze(W(:,k,:)) * squeeze(U(:,k,:))';
        newnorm = sum(wu(:).^2).^.5;
        W(:,k,:) = W(:,k,:)/newnorm;
    end
end

% compute adjacency matrix UtU
U(isnan(U)) = 0;
U0 = gpuArray(U);
utu = gpuArray.zeros(Nfilt, 'single');
for irank = 1:Nrank
    utu = utu + (U0(:,:,irank)' * U0(:,:,irank));
end

UtU = logical(utu);

author: cortex-lab

repo: KiloSort

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


text:
%%DEMOGRAVCODENOCOMPILE This file demonstrates the functions for spherical
%                       harmonic synthesis of gravitational and terrain
%                       data. Unlike the function DemoGravCode, this
%                       function does not make any plots, but does run fast
%                       even when the mex files have not ben compiled. The
%                       gravitational code library (that this function is
%                       part of) must have been added to Matlab's search
%                       path for this function to work.
%
%Rather than using advanced spherical harmonic models at high orders, this
%function loads spherical harmonic coefficients only up to degree and order
%50 for a simple ellipsoidal Earth model. The gravitational potential and
%acceleration due to gravity from the ellipsoidal model as generated by the
%coefficients are shown to match (within reasonable precision bounds) the
%results obtained from closed-form functions.
%
%When run, the error in the potential should be 5.9474e-14% and the error
%in the magnitude of the gravitational acceleration should be 3.3790e-12%.
%
%February 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=50;%The maximum degree and order of the model.
%Generate coefficients for an ellipsoidal Earth model.
[C,S,a,c]=ellipsGravCoeffs(M);

lambda=-20*pi/180;%Longitude
theta=80*pi/180;%Ellipsoidal latitude
%The zero altitude places the point on the reference ellipsoid.
pointEllips=[theta;lambda;0];
pointSpher=ellips2Sphere(pointEllips);%Convert to spherical coordinates.

%Compute the potential and acceleration due to gravity.
[V,gradV]=spherHarmonicEval(C,S,pointSpher,a,c);

%Add in the centripetal effects of the Earth's rotation. The formulae for
%the extra terms in the potential and acceleration due to gravity can be
%found in Chapters 2.1 and 2.7 of
%B. Hofmann-Wellenhof and H. Moritz, Physical Geodesy, 2nd ed.
%SpringerWienNewYork, 2006.
CartPoint=ellips2Cart(pointEllips);%Convert to Cartesian coordinates.
omega=Constants.WGS84EarthRotationRate;
U=V+0.5*omega^2*(CartPoint(1)^2+CartPoint(2)^2);%The gravity potential
%Acceleration due to gravity with the Earth's rotation.
g=gradV+omega^2*[CartPoint(1);CartPoint(2);0];

%Get the exact values of U and g from explicit formulas for the ellipsoidal
%Earth approximation.
[U0,g0]=ellipsParam2Grav(CartPoint);

disp('The geodetic latitude and longitude of the point on the surface of')
disp('the WGS-84 reference ellipsoid where the values are computed are (in degrees)')
latitude=theta*(180/pi)
longitude=lambda*(180/pi)

disp('The gravity potential from the spherical harmonic solution is (in m^2/s^2)')
U

disp('The gravitational acceleration vector (with centripetal acceleration)')
disp('from the spherical harmonic solution is (in m/s^2)')
g

disp('The percent error in the potential from the spherical harmonic')
disp('approximation and from the explicit solution is')
percentErrorPotential=100*(U-U0)/U0

disp('The percent error in the gravitational acceleration magnitude from')
disp('the spherical harmonic approximation and from the explicit solution is')
percentErrorGravAccel=100*(norm(g)-norm(g0))/norm(g0)

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Test file for chebtech2/chebpts.m

function pass = test_chebpts(varargin)

% Set a tolerance (pref.chebfuneps doesn't matter)
tol = 10*eps;

% Test that n = 0 returns empty results:
[x, w, v] = chebtech2.chebpts(0);
pass(1) = ( isempty(x) && isempty(w) && isempty(v) );

% Test n = 1:
[x, w, v] = chebtech2.chebpts(1);
pass(2) = ( x == 0 && w == 2 && v == 1);

% Test that n = 2 returns [-1 ; 1]:
x = chebtech2.chebpts(2);
pass(3) = ( all(size(x) == [2, 1]) && all(x == [-1 ; 1]) );
[x, w, v] = chebtech2.chebpts(2);
pass(4) = ( all(size(w) == [1, 2]) && all(w == [1, 1]) );
pass(5) = ( all(size(v) == [2, 1]) && all(v == .5*[-1 ; 1]) );

% Test that n = 3 returns [-1 ; 0 ; 1]:
x = chebtech2.chebpts(3);
pass(6) = ( all(size(x) == [3, 1]) && all(x == [-1 ; 0 ; 1]) );
[x, w, v] = chebtech2.chebpts(3);
pass(7) = ( all(size(w) == [1, 3]) && norm( w - ([0, 1, 0] + 1/3), inf) < tol );
pass(8) = ( all(size(v) == [3, 1]) && norm( v - ([.5 ; -1 ; .5]) , inf) < tol );

% Test that n = 129 returns vectors of the correct size:
n = 129;
[x, w, v] = chebtech2.chebpts(n);
pass(9) = ( all([size(x) == [n, 1], size(w) == [1, n], size(v) == [n, 1]]) );
% and that the nodes are symmetric:
pass(10) = ( norm(x(1:(n-1)/2) + x(n:-1:(n+3)/2), inf) == 0 );
pass(11) =  ( x((n+1)/2) == 0 );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% FindInl    find inliers in joint image matrix
%			 by pairwise epipolar geometry
%
% function IdMatIn = findinl(Ws,IdMat,tol)
% Ws ... 3MxN joint image matrix
% IdMat ... MxN ... 0 -> no point detected
%                   1 -> point detected
% tol ... [pixels] tolerance for the epipolar geometry
%         the point are accpted as outliers only if they
%         are closer to the epipolar line than tol

% $Author: svoboda $
% $Revision: 2.1 $
% $Id: findinl.m,v 2.1 2003/07/30 10:28:29 svoboda Exp $
% $State: Exp $

function IdMatIn = findinl(Ws,IdMat,tol)

NoCams = size(IdMat,1);

% fill the array of structures not_used denoted as 0
% allocate the array of structures for used
for i=1:NoCams,
  not_used(i).pts = sum(IdMat(i,:));
  used(i).pts	  = -1;
end

% allocate IdMat for outliers
IdMatIn = zeros(size(IdMat));

while (sum([not_used.pts])>1-NoCams),
  [buff, id.cam_max]  = max([not_used.pts]);
  used	   = add(used, id.cam_max, not_used(id.cam_max).pts);
  not_used = remove(not_used, id.cam_max);
  Mask	   = repmat(IdMat(id.cam_max,:),NoCams,1);
  Corresp  = Mask & IdMat;
  Corresp(id.cam_max,:) = 0;
  [buff, id.cam_to_pair] = max(sum(Corresp')); % find the camera with most correspondences
  idx.corr_to_pair = find(sum(IdMat([id.cam_max,id.cam_to_pair],:))==2);
  % used	   = add(used, id.cam_to_pair, not_used(id.cam_to_pair).pts);
  % not_used = remove(not_used, id.cam_to_pair);
  if size(idx.corr_to_pair,2)<8,
	error('Not enough points to compute epipolar geometry in RANSAC validation')
  end
  Wspair   = [];
  Wspair   = Ws(id.cam_max*3-2:id.cam_max*3, idx.corr_to_pair);
  Wspair   = [Wspair; Ws(id.cam_to_pair*3-2:id.cam_to_pair*3, idx.corr_to_pair)];
  % id
  [F, inls] = rEG(Wspair,tol,tol,0.99);
  IdMatIn(id.cam_max, idx.corr_to_pair(inls)) = 1;
  IdMat(id.cam_max, :)						  = 0;
  IdMat(id.cam_max, idx.corr_to_pair(inls))	  = 1;
end

function list = add(list, id, value)
list(id).pts = value;
return

function list = remove(list, id)
list(id).pts = -1;
return

author: strawlab

repo: MultiCamSelfCal


text:
function value = r8vec_diff_norm ( n, a, b )

%*****************************************************************************80
%
%% R8VEC_DIFF_NORM returns the L2 norm of the difference of R8VEC's.
%
%  Discussion:
%
%    The vector L2 norm is defined as:
%
%      value = sqrt ( sum ( 1 <= I <= N ) A(I)**2 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in A.
%
%    Input, real A(N), B(N), the vectors.
%
%    Output, real VALUE, the L2 norm of A - B.
%
  value = norm ( a(1:n) - b(1:n) );

  return
end

repo: jburkardt-m


text:
% PROPAGATE A GAUSSIAN ERROR DISTRIBUTION OF EACH JOINT TO AN ERROR IN
% END EFFECTORS' POSITION
% AN ELLIPSE IS DRAWN TO REPRESENT THE ERROR IN POSITION

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

%load arm parameters
robot= load_robot('example','scara')

%find errors around pose
q=[pi/4 3*pi/4 0 0];

%Matriz de errores en las articulaciones
% sigmaq1=sigmaq2=0.0017 rad, sigmaq3=0.01 m.
sigmaq1=0.017;%rad, 0.01 grados
sigmaq2=0.017;%rad
sigmaq3=0.01;% m
Rq=[sigmaq1^2 0 0;
    0 sigmaq2^2 0;
    0 0 sigmaq3^2]

teta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alfa = eval(robot.DH.alpha);

Jq = eval(robot.J)

Rp=Jq*Rq*Jq'

T=directkinematic(robot, q);

drawrobot3d(robot,q), hold on
draw_ellipse([T(1,4),T(2,4)], Rp, 'r')

repo: ARTE


text:
%  Figure 3.18      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.18
%% fig3_18.m
clf;
num=[1];          % form transfer function
den=[1 2*.2 1];
t=0:.01:30;       % define time vector
sys=tf(num,den);  % define system
y=impulse(sys,t); % compute impulse response
plot(t,y)
hold on
% compute e^(-.2t)
ne=1;
de=[1 .2];
syse=tf(ne,de);
ye=impulse(syse,t);
plot(t,ye,':r')
plot(t,-ye,':r'), grid
axis([0 30 -1 1])
hold off
title('Figure 3.18')
xlabel('Time (sec)')
ylabel('h(t)')
text(7,.3,'e^{-\sigma t}');
text(7,-.3,'-e^{-\sigma t}');