text:
function talDat = loadTalairachCoordinates;

%loadTalairachCoordinates() Load Tal. coordinates from file.
%
%This function is supposed to prevent multiple hacking of the same coordinates
%instead write to a pure textfile and read them in for each subject.
%The file is assumed to reside at the base directory for the antomies. Hoever deviations
%in file location are posible
%the structure of the text-file is fixed and assumed to be the following:
%x-coord ycoord zcoord descripition-string color
%The file may consist of several lines each specifiying one set of coordinates
%Fields within a line are separated by whitespace(s). The field color is a valid matlab
%colorstring such as 'b' 'r' 'y' etc.
%
% Authors BD and JR 09/18/01

% ugly but conforms with the conventions
rootPath=which(mfilename);
rootPath = fullfile(fileparts(rootPath),'Coords');
pathstr = getPathStrDialog(rootPath,'Choose Talairach-coordinates file','*.txt;*.TXT');

if isempty(pathstr)
    talDat = [];
end

talDat = dlmreadMixedDataType(pathstr,'n1 n2 n3 s4 s5');


text:
function [weight] = bcrp_kernel(data, weight_o)
% BCRP for ps, kernel file
% Input
%     data    -- market historical prices
%     weight_o -- portfolio weight for last trading day
%
% Output
%     weight  -- weight column vector
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

return;
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [pnt, tri] = read_off(fn)

% READ_OFF reads vertices and triangles from a OFF format triangulation file
%
% [pnt, tri] = read_off(filename)
%
% See also READ_TRI, READ_BND

% Copyright (C) 1998, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(fn, 'rt');

% scan the file type
[s, count] = fscanf(fid, '%s\n', 1);
if ~strcmp(s,'OFF')
  msg = sprintf('wrong file type %s', s);
  ft_error(msg);
end

% read the number of vertex points and triangles
[val, count] = fscanf(fid, '%d', 3);
Npnt = val(1)
Ntri = val(2)

% read the vertex points
pnt  = fscanf(fid, '%f', [3, Npnt]);
pnt  = pnt(1:3,:)';

% read the triangles
tri = fscanf(fid, '%d', [4, Ntri]);
tri = (tri(2:4,:)+1)';
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Alpha, Trf] = bSSFP_GetSeq( VaryAlpha, FixTrf, VaryTrf, FixAlpha )

Alpha =  zeros(length(VaryTrf)+length(VaryAlpha), 1);
Trf   =  Alpha;

kk=1;
for ii = 1:length(VaryAlpha)
    Alpha(kk) = VaryAlpha(ii);
    Trf(kk) = FixTrf;
    kk = kk+1;
end

for jj = 1:length(VaryTrf)
    Alpha(kk) = FixAlpha;
    Trf(kk) = VaryTrf(jj);
    kk = kk+1;
end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [idi, rt, rn, bb] = idis(pt, pn, y, z, t, varargin)
%IDIS Integrated Discrimination Improvement between two models
%
%  Description
%    [IDI,RT,RN,BB] = IDIS(PT,PN,Y,Z,T,OPTIONS) Returns Integrated
%    Discrimination Improvement (IDI) given two vectors of event probabilities
%    at time T, PT for traditional model and PN for new model, and the
%    observed event or censoring times Y and the corresponding censoring
%    indicators Z (0=event, 1=censored). With these inputs, the estimator
%    attributed to Pencina et al. in the reference is used, with
%    restriction to time T.
%
%    [IDI,RT,RN,BB] = IDIS(PT,PN,OPTIONS) Returns Integrated
%    Discrimination Improvement (IDI) given two vectors of event probabilities
%    at time T, PT for traditional model and PN for new model. Here, the
%    model-based estimator is used (the "new estimator" in the reference).
%
%    Ouputs RT and RN are the R^2 statistics for the two models. BB are
%    Bayesian bootstrap samples of the IDI distribution.
%
%    OPTIONS is optional parameter-value pair
%      rsubstream - number of a random stream to be used for
%                   simulating dirrand variables. This way same
%                   simulation can be obtained for different models.
%                   See doc RandStream for more information.
%
%  Reference
%    L. E. Chambless, C. P. Cummiskey, and G. Cui (2011). Several
%    methods to assess improvement in risk prediction models:
%    Extension to survival analysis. Statistics in Medicine
%    30(1):22-38.

% Copyright (C) 2012 Ernesto Ulloa, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
if nargin < 3 || ischar(y)
    % model-based estimator
    model_based_estimator = true;

    ip.addRequired('pt', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
    ip.addRequired('pn', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
    ip.addParamValue('rsubstream', 0, @(x) isreal(x) && isscalar(x) && isfinite(x) && x>0)
    if nargin > 2
        % if more than 2 input arguments, there must be four as only a
        % single optional parameter is implemented
        if nargin == 4
            ip.parse(pt, pn, y, z);
        else
            error('Invalid number of arguments.');
        end
    else
        ip.parse(pt, pn);
    end
else
    % Pencina et al. estimator
    model_based_estimator = false;

    ip.addRequired('pt', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
    ip.addRequired('pn', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
    ip.addRequired('y', @(x) isreal(x) && all(isfinite(x(:))))
    ip.addRequired('z', @(x) isreal(x) && all(isfinite(x(:))))
    ip.addRequired('t', @(x) isreal(x) && isscalar(x) && ~isnan(x))
    ip.addParamValue('rsubstream', 0, @(x) isreal(x) && isscalar(x) && isfinite(x) && x>0)
    ip.parse(pt, pn, y, z, t, varargin{:})
end
rsubstream=ip.Results.rsubstream;

if nargout < 4
    % without boostrap
    if model_based_estimator
        rt = rsqr(pt);
        rn = rsqr(pn);
    else
        rt = rsqr(pt, y, z, t);
        rn = rsqr(pn, y, z, t);
    end
else
    % with boostrap
    if model_based_estimator
        if rsubstream == 0
            [rt, bbt] = rsqr(pt);
            [rn, bbn] = rsqr(pn);
        else
            [rt, bbt] = rsqr(pt, 'rsubstream', rsubstream);
            [rn, bbn] = rsqr(pn, 'rsubstream', rsubstream);
        end
    else
        if rsubstream == 0
            [rt, bbt] = rsqr(pt, y, z, t);
            [rn, bbn] = rsqr(pn, y, z, t);
        else
            [rt, bbt] = rsqr(pt, y, z, t, 'rsubstream', rsubstream);
            [rn, bbn] = rsqr(pn, y, z, t, 'rsubstream', rsubstream);
        end
    end
    bb = bbn - bbt;
end

idi = rn - rt;

end


text:
% LFDispVidCirc - visualize a 4D light field animating a circular path through two dimensions
%
% Usage:
%     FigureHandle = LFDispVidCirc( LF )
%     FigureHandle = LFDispVidCirc( LF, [RenderOptions], [ScaleFactor] )
%
%
% Inputs:
%
%     LF : a colour or single-channel light field, and can a floating point or integer format. For
%          display, it is converted to 8 bits per channel. If LF contains more than three colour
%          channels, as is the case when a weight channel is present, only the first three are used.
%
% Optional Inputs:
%
%     RenderOptions : struct controlling rendering
%        .PathRadius_percent : radius of the circular path taken by the viewpoint. Values that are
%                              too high can result in collision with the edges of the lenslet image,
%                              while values that are too small result in a less impressive
%                              visualization. The default value is 60%.
%        .FrameDelay         : delay between frames, in seconds. Default 1/30.
%        .NumCycles          : how many times the circular path should repeat. Default inf.
%        .FramesPerCycle     : how many frames along a single cycle. Default 20.
%
%     ScaleFactor : Adjusts the size of the display -- 1 means no change, 2 means twice as
%                   big, etc. Integer values are recommended to avoid scaling artifacts.
%
% Outputs:
%
%     FigureHandle
%
%
% User guide: LFToolbox.pdf
% See also:  LFDisp, LFDispMousePan, LFDispLawnmower, LFDispProj, LFDispTiles

% Copyright (c) 2013-2020 Donald G. Dansereau

function FigureHandle = LFDispVidCirc( LF, RenderOptions, ScaleFactor )

%---Defaults---
ScaleFactor = LFDefaultVal('ScaleFactor', 1);
RenderOptions = LFDefaultField( 'RenderOptions', 'PathRadius_percent', 60 );
RenderOptions = LFDefaultField( 'RenderOptions', 'FrameDelay', 1/30 );
RenderOptions = LFDefaultField( 'RenderOptions', 'NumCycles', inf );
RenderOptions = LFDefaultField( 'RenderOptions', 'FramesPerCycle', 20 );

%---Check for mono and clip off the weight channel if present---
Mono = (ndims(LF) == 4);
if( ~Mono )
    LF = LF(:,:,:,:,1:3);
end

%---Rescale for 8-bit display---
if( isfloat(LF) )
    LF = uint8(LF ./ max(LF(:)) .* 255);
else
    LF = uint8(LF.*(255 / double(intmax(class(LF)))));
end

%---Setup the motion path---
[TSize,SSize, ~,~] = size(LF(:,:,:,:,1));
TCent = (TSize-1)/2 + 1;
SCent = (SSize-1)/2 + 1;

t = 0;
RotRad = TCent*RenderOptions.PathRadius_percent/100;
NumFrames = RenderOptions.NumCycles * RenderOptions.FramesPerCycle;

%---Setup the display---
TIdx = round(TCent + RotRad);
SIdx = round(SCent);
[ImageHandle,FigureHandle] = LFDispSetup( squeeze(LF(TIdx,SIdx,:,:,:)), ScaleFactor );

while(1)
    TVal = TCent + RotRad * cos( 2*pi*t/RenderOptions.FramesPerCycle );
    SVal = SCent + RotRad * sin( 2*pi*t/RenderOptions.FramesPerCycle );

    SIdx = round(SVal);
    TIdx = round(TVal);

    CurFrame =  squeeze(LF( TIdx, SIdx, :,:,: ));
    set(ImageHandle,'cdata', CurFrame );

    pause(RenderOptions.FrameDelay)
    t = t + 1;

	if( t > NumFrames )
		break;
	end
end

repo: LFToolbox


text:
function tran03_values_test ( )

%*****************************************************************************80
%
%% TRAN03_VALUES_TEST demonstrates the use of TRAN03_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRAN03_VALUES_TEST:\n' );
  fprintf ( 1, '  TRAN03_VALUES stores values of\n' );
  fprintf ( 1, '  the transportation function of order 3.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X               FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = tran03_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%                               manageOptionsForBenchmark.m
%
% This file consists of functions used to check the validity of
% user-provided options, provideds default values
% for those unspecified options and raises warnings for those unnecessary
% but user-provided fields.
%
% manageOptionsForBenchmark invokes helper functions getExtraOpts,
% applyOpts, and warnExtraOpts in the folder util/.
% It is used in the general benchmark interface benchmarkPR.m.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

% This function integrates the user specified values and default values for
% fields in opts and raises warning for those unnecessary but user-provided
% fields.
% params is as defined in benchmarkPR.m. See its header or User Guide for
% details.
function params = manageOptionsForBenchmark(dataSet, params)

    % Obtain and apply default options recognized by a dataSet
    defaults = getDefaultOptsForBenchmack(dataSet);
    extras = getExtraOpts(params, defaults);
    params = applyOpts(params, defaults, false);

    % Warn the user regarding unrecognized options
    if ~strcmp(dataSet, 'custom')
        warnExtraOpts(extras);
    end
end

% Return a struct that consists of all the special field(initialized to default values)
% for the dataSet
function params = getDefaultOptsForBenchmack(dataSet)
    % General options
    params.verbose = true;           % If true, each trial's result will be reported.
    params.plotType = 'auto';    % The type of the plot generated. It
                                     % currently supports ['semilogy', 'linear','auto'].
    params.numTrials = 1;             % The number of trials each algorithm/dataset
                                     % combination will run.
    params.policy = 'median';        % How to compute the final yvalue used
                                     % for ploting from the yvalues one get
                                     % by running numTrials trials. It
                                     % currently supports
                                     % ['median','average','best',
                                     % 'successRate']
    params.successConstant = 1e-5;   % If the yvalue of the current trial
                                     % is less than it,
                                     % the trial will be counted as a
                                     % success. This parameter will only be
                                     % used when policy='successRate'.
    params.maxTime = 300;            % max time allowed for a single
                                     % algorithm/dataset trial.

    params.recordSignals = false;    % whether record the recovered signal
                                     % at each trial.

    % Specific options
    switch lower(dataSet)
        case '1dgaussian'
            params.n = 10;                        % Length of the unknown
                                                  % signal
            params.m = 80;                        % Number of measurements
            params.isComplex = true;              % If the signal is complex
            params.isNonNegativeOnly = false;     % If the signal is real
                                                  % and non-negative
            params.SNR = inf;                     % Signal-to-noise ratio
        case '2dimage'
            params.imagePath = 'data/shapes.png'; % Image path
            params.L = 12;                        % Number of of masks
                                                  % created
        case 'transmissionmatrix'
            params.n=256;
            params.m=256*20;
            params.isComplex = true;              % If the signal is complex
            params.isNonNegativeOnly = false;     % If the signal is real
                                                  % and non-negative
        otherwise
            error('unknown dataset: %s\n',dataSet);
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
function fMeasCreateLUT_main(sFilename)
%FMEASCREATELUT Creates look-up table from measdata file(s).
%   FMEASCREATELUT(SFILENAME) parses the given Siemens measdata
%   file and creates a look-up table of the file's MDH entries and saves it
%   to a file '%sFilename%.mat'. This information can be used to load
%   measdata from that file later on using fMeasRead way faster than
%   usually. Additionally, only MDHs with specific loop counter values can
%   be read.
%
%   FMEASCREATELUT looks for '.dat'-files in the current working directory
%   and creates look-up tables/files for all Siemens measdata files.
%
%   Class support for input sFilename:
%      string
%
%   See also FMEASINFO, FMEASREAD
%
%   Copyright 2014 Christian Wuerslin, University of Tuebingen, Germany
%   christian.wuerslin@med.uni-tuebingen.de
%   Modifications 2014 Thomas Kuestner, University of Tuebingen, Germany
%   thomas.kuestner@med.uni-tuebingen.de

% csLCS = {'Smp', 'Cha', 'Lin', 'Acq', 'Sli', 'Par', 'Eco', 'Pha', 'Rep', 'Set', 'Seg', 'FreePara1', 'FreePara2', 'FreePara3', 'FreePara4'};
% Loop Counter names
%
% (c) Christian Wuerslin, Thomas Kuestner
% ---------------------------------------------------------------------

if nargin

    % ---------------------------------------------------------------------
    % Filename give => parse file

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Create the involved filenames
    [sPath, sName] = fileparts(sFilename);
    if isempty(sPath)
        sFullName = which(sFilename);
        if isempty(sFullName), error('File ''%s'' not on path!', sFilename); end
        [sPath, sName, sExt] = fileparts(sFullName);
        sFilename = [sPath, filesep, sName, sExt];
    else
        if ~exist(sFilename, 'file'), error('File ''%s'' doesn''t exist!', sFilename); end
    end
    sFilenameXML = [sPath, filesep, sName, '.xml'];
    sFilenameMat = [sPath, filesep, sName, '.mat'];
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Get filesize and extract the header
    SDir = dir(sFilename);
    iFileSize = SDir(1).bytes;
    fidmeas = fopen(sFilename, 'r');
    iHeaderLength = fread(fidmeas, 1, 'uint32');
    fprintf(1, '\n*** Parsing %s ***\n', sFilename);
    fprintf(1, 'File Size            : %u  (%u MB)\n', iFileSize, round(iFileSize/(1024^2)));
    fprintf(1, 'Length of the Header : %u\n', iHeaderLength);

    sXML = fread(fidmeas, iHeaderLength, '*char');
    fidXML = fopen(sFilenameXML, 'w');
    fwrite(fidXML, sXML);
    fclose(fidXML);
    fclose(fidmeas);
    fprintf(1, 'Header exported to file ''%s''\n', sFilenameXML);
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Get the MDHs
    [dMDH, iCount] = fMeasCreateLUT_mex(sFilename);
    iCount = uint32(iCount);
    dMDH = dMDH(:,1:iCount)';
    fprintf(1, 'Number of ADCs: %u\n', iCount);
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Get the DrecksMDH as struct
    lDrecksMDH = dMDH(:, 4) == 20;
    if(nnz(lDrecksMDH) == 0)
        lVer = -1;
    else
        lVer = bitget(uint16(dMDH(find(lDrecksMDH, 1, 'first'), 6)), 16);
    end
    if  lVer == -1 % no DrecksMDH
        iLCDrecks = [];
        iLCPositions = dMDH(~lDrecksMDH, 32:38);
    elseif lVer == 0 % or dMDH(1,6) ?, bitget(uint16(dMDH(1, 1)), 16) == 0
        % DrecksMDH 1.0
        iLCDrecks = [dMDH(lDrecksMDH,4:19),zeros(nnz(lDrecksMDH),4),dMDH(lDrecksMDH,30),zeros(nnz(lDrecksMDH),1),dMDH(lDrecksMDH,20:29),dMDH(lDrecksMDH,31:37),zeros(nnz(lDrecksMDH),1),dMDH(lDrecksMDH,38)];
        iLCPositions = dMDH(~lDrecksMDH, 31:end);
    else
        % DrecksMDH 2.0
        dDrecksMDH = dMDH(lDrecksMDH, :);
        iNChannels = dDrecksMDH(1, 5);
        dDrecksMDH = dDrecksMDH(1:iNChannels:end, :);
        dInt = dDrecksMDH(:, 6:28)';
        dFloat = dDrecksMDH(:, 29:30)';
        dDrecksMDH = [dInt(:); dFloat(:)];
        SDrecksMDH = fReadNewDrecksMDH(dDrecksMDH);
    end
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    iSP = uint64(dMDH(~lDrecksMDH,1));
    iEvalInfoMask = uint32(dMDH(~lDrecksMDH, 2:3));
    iLC = uint16(dMDH(~lDrecksMDH, [4:14,26:29]));
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    if (lVer == -1)
        save(sFilenameMat, 'iLC', 'iSP', 'iEvalInfoMask', 'iLCPositions');
        fMeasInfo(sFilename);
    elseif (lVer == 0)   % or dMDH(1,6) ?, bitget(uint16(dMDH(1, 1)), 16) == 0   %
        % DrecksMDH 1.0
        save(sFilenameMat, 'iLC', 'iSP', 'iEvalInfoMask', 'iLCDrecks', 'iLCPositions');
        fMeasInfo(sFilename);
    else
        % DrecksMDH 2.0
        dLCPositions = dMDH(~lDrecksMDH, 31:end);
        SPos.dSliceDataPosition = dLCPositions(:,1:3);
        if SDrecksMDH.Seq.Is3D
            SPos.dSliceDataPosition = unique(SPos.dSliceDataPosition, 'rows');
        end
%         SPos.dTablePosition     = unique(-((dLCPositions(:,4) +5)/10)); % wrong: that is shift
        SPos.dTablePosition     = unique(-dLCPositions(:,8)/10);
        SPos.dRotMatrix         = quat2dcm(dLCPositions(1,4:7));
        SPos.dRotMatrix(:,2:3)  = -SPos.dRotMatrix(:,2:3); % due to negative orientation of y and z axis (compared to standard cartesian system)
        [yaw, pitch, roll]      = dcm2angle(SPos.dRotMatrix);
        SPos.dRotAngles         = [yaw, pitch, roll];

        save(sFilenameMat, 'iLC', 'iSP', 'iEvalInfoMask', 'SDrecksMDH', 'SPos');
        fMeasInfo2(sFilename);       % Print loop counter summary
    end
else
    % No input argument: Execute on all *.dat files in folder (batch processing)
    fprintf(1, '\n***********************\n');
    fprintf(1, '* Entering batch-mode *\n');
    fprintf(1, '*   Get some coffee!  *\n');
    fprintf(1, '***********************\n\n');
    SDir = dir('*.dat');
    for iI = 1:length(SDir)
        fMeasCreateLUT_main(SDir(iI).name);
    end
end

author: thomaskuestner


text:
function [label, centroid, dis] = fkmeans(X, k, options)
% FKMEANS Fast K-means with optional weighting and careful initialization.
% [L, C, D] = FKMEANS(X, k) partitions the vectors in the n-by-p matrix X
% into k (or, rarely, fewer) clusters by applying the well known batch
% K-means algorithm. Rows of X correspond to points, columns correspond to
% variables. The output k-by-p matrix C contains the cluster centroids. The
% n-element output column vector L contains the cluster label of each
% point. The k-element output column vector D contains the residual cluster
% distortions as measured by total squared distance of cluster members from
% the centroid.
%
% FKMEANS(X, C0) where C0 is a k-by-p matrix uses the rows of C0 as the
% initial centroids instead of choosing them randomly from X.
%
% FKMEANS(X, k, options) allows optional parameter name/value pairs to
% be specified. Parameters are:
%
%   'weight' - n-by-1 weight vector used to adjust centroid and distortion
%              calculations. Weights should be positive.
%   'careful' - binary option that determines whether "careful seeding"
%               as recommended by Arthur and Vassilvitskii is used when
%               choosing initial centroids. This option should be used
%               with care because numerical experiments suggest it may
%               be counter-productive when the data is noisy.
%
% Notes
% (1) The careful seeding procedure chooses the first centroid at random
% from X, and each successive centroid from the remaining points according
% to the categorical distribution with selection probabilities proportional
% to the point's minimum squared Euclidean distance from the already chosen
% centroids. This tends to spread the points out more evenly, and, if the
% data is made of k well separated clusters, is likely to choose an initial
% centroid from each cluster. This can speed convergence and reduce the
% likelihood of getting a bad solution [1]. However, in experiments where
% 5% uniformly distributed noise data was added to such naturally clustered
% data the results were frequently worse then when centroids were chosen at
% random.
% (2) If, as is possible, a cluster is empty at the end of an iteration,
% then there may be fewer than k clusters returned. In practice this seems
% to happen very rarely.
% (3) Unlike the Mathworks KMEANS this implementation does not perform a
% final, slow, phase of incremental K-means ('onlinephase') that guarantees
% convergence to a local minimum.
%
% References
% [1] "k-means++: The Advantages of Careful Seeding", by David Arthur and
% Sergei Vassilvitskii, SODA 2007.

n = size(X,1);

% option defaults
weight = 0; % uniform unit weighting
careful = 0;% random initialization

if nargin == 3
    if isfield(options, 'weight')
        weight = options.weight;
    end
    if isfield(options,'careful')
        careful = options.careful;
    end
end

% If initial centroids not supplied, choose them
if isscalar(k)
    % centroids not specified
    if careful
        k = spreadseeds(X, k);
    else
        k = X(randsample(size(X,1),k),:);
    end
end

% generate initial labeling of points
[~,label] = max(bsxfun(@minus,k*X',0.5*sum(k.^2,2)));
k = size(k,1);

last = 0;

if ~weight
    % code defactoring for speed
    while any(label ~= last)
        % remove empty clusters
        [~,~,label] = unique(label);
        % transform label into indicator matrix
        ind = sparse(label,1:n,1,k,n,n);
        % compute centroid of each cluster
        centroid = (spdiags(1./sum(ind,2),0,k,k)*ind)*X;
        % compute distance of every point to each centroid
        distances = bsxfun(@minus,centroid*X',0.5*sum(centroid.^2,2));
        % assign points to their nearest centroid
        last = label;
        [~,label] = max(distances);
    end
    dis = ind*(sum(X.^2,2) - 2*max(distances)');
else
    while any(label ~= last)
        % remove empty clusters
        [~,~,label] = unique(label);
        % transform label into indicator matrix
        ind = sparse(label,1:n,weight,k,n,n);
        % compute centroid of each cluster
        centroid = (spdiags(1./sum(ind,2),0,k,k)*ind)*X;
        % compute distance of every point to each centroid
        distances = bsxfun(@minus,centroid*X',0.5*sum(centroid.^2,2));
        % assign points to their nearest centroid
        last = label;
        [~,label] = max(distances);
    end
    dis = ind*(sum(X.^2,2) - 2*max(distances)');
end
label = label';

% Code below this line reused from the file exchange submission K-means++
% (http://www.mathworks.com/matlabcentral/fileexchange/28901-k-means) in
% accordance with the license:
% vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
% Copyright (c) 2010, Michael Chen
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * 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 OWNER 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 D = sqrdistance(A, B)
% Square Euclidean distances between all sample pairs
% A:  n1 x d data matrix
% B:  n2 x d data matrix
% WB: n2 x 1 weights for matrix B
% D: n2 x n1 pairwise square distance matrix
%    D(i,j) is the squared distance between A(i,:) and B(j,:)
% Written by Michael Chen (sth4nth@gmail.com). July 2009.
n1 = size(A,1); n2 = size(B,2);
m = (sum(A,1)+sum(B,1))/(n1+n2);
A = bsxfun(@minus,A,m);
B = bsxfun(@minus,B,m);
D = full((-2)*(A*B'));
D = bsxfun(@plus,D,full(sum(B.^2,2))');
D = bsxfun(@plus,D,full(sum(A.^2,2)))';
end

function [S, idx] = spreadseeds(X, k)
% X: n x d data matrix
% k: number of seeds
% reference: k-means++: the advantages of careful seeding.
% by David Arthur and Sergei Vassilvitskii
% Adapted from softseeds written by Mo Chen (mochen@ie.cuhk.edu.hk),
% March 2009.
[n,d] = size(X);
idx = zeros(k,1);
S = zeros(k,d);
D = inf(n,1);
idx(1) = ceil(n.*rand);
S(1,:) = X(idx(1),:);
for i = 2:k
    D = min(D,sqrdistance(S(i-1,:),X));
    idx(i) = find(cumsum(D)/sum(D)>rand,1);
    S(i,:) = X(idx(i),:);
end
end

end

author: Sable

repo: mcbench-benchmarks


text:
% SB2_PREPROCESSBASIS  Do basis matrix pre-processing for SPARSEBAYES
%
% [BASIS, SCALES] = SB2_PREPROCESSBASIS(BASIS)
%
% OUTPUT ARGUMENTS:
%
%	BASIS	NxM matrix of basis vectors appropriately pre-processed
%			(scaled to unit length per column)
%
%	SCALES	Vector of scaling factors applied
%
%
% INPUT ARGUMENTS:
%
%	BASIS	NxM matrix of basis vectors (one column per basis function)
%
% NOTES:
%
% Simply normalises the basis vectors to unit length, returning the
% original lengths so that the weights can be rescaled appropriately
% before returning to the user.
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the Free
% Software Foundation; either version 2 of the License, or (at your option)
% any later version.
%
% SPARSEBAYES is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
% FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
% more details.
%
% You should have received a copy of the GNU General Public License along
% with SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function [BASIS, Scales] = SB2_PreProcessBasis(BASIS)

%
[N,M]	= size(BASIS);
%
% Compute "lengths" of basis vectors (columns of BASIS)
%
Scales	= sqrt(sum(BASIS.^2));
%
% Work-around divide-by-zero inconvenience
%
Scales(Scales==0)	= 1;
%
% Normalise each basis vector to "unit length"
%
for m=1:M
  BASIS(:,m)	= BASIS(:,m) / Scales(m);
end

repo: REST


text:
function [y,feature_names,cache] = ComputeHarmonicWindowFeatures(x,varargin)

x = x(:)';
N = numel(x);
y = nan(0,N);
feature_names = {};

%% default parameters

[windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii] = ...
  SetDefaultWindowParameters();

% use all transformation types by default
%trans_types = 'all';
trans_types = uint8(15);

% for debugging purposes
SANITY_CHECK = true;

% whether to cache results
DOCACHE = true;

% initialize empty cache
cache = InitializeCache();

% initialize feature_types already computed to empty
feature_types = {};

% harmonic features
num_harmonic = 1;

% relative params
relativeParams = [];
%% parse parameters

[...
  windows,...
  window_radii,window_offsets,...
  min_window_radius,max_window_radius,nwindow_radii,...
  feature_types,...
  trans_types,...
  SANITY_CHECK,...
  DOCACHE,...
  cache,...
  num_harmonic,...
  relativeParams...
  ] = myparse(varargin,...
  'windows',windows,...
  'window_radii',window_radii,'window_offsets',window_offsets,...
  'min_window_radius',min_window_radius,'max_window_radius',max_window_radius,'nwindow_radii',nwindow_radii,...
  'feature_types',feature_types,...
  'trans_types',trans_types,...
  'sanitycheck',SANITY_CHECK,...
  'docache',DOCACHE,...
  'cache',cache,...
  'num_harmonic',num_harmonic,...
  'relativeParams',relativeParams); %#ok<ASGLU>

%% whether we've specified to use all trans types by default
%if ischar(trans_types) && strcmpi(trans_types,'all'),
%  trans_types = {'none','abs','flip','relative'};
%end

%% select default windows from various ways of specifying windows

[windows,window_radii,windowi2radiusi,nradii] = ...
  SetWindowParameters(...
  windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii);

%% main computation

%if ismember('relative',trans_types)
if bitand(8,trans_types)
  if DOCACHE && ~isempty(cache.relX)
    modX = cache.relX;
  else
    modX = convertToRelative(x,relativeParams);
    cache.relX = modX;
  end
end

for radiusi = 1:nradii,
  r = window_radii(radiusi);
  w = 2*r+1;

  for num_harmonic_curr = 1:num_harmonic,

    % can't split into more segments than frames
    if num_harmonic_curr >= w,
      continue;
    end

    res = HarmonicWindowCore(x,w,num_harmonic_curr);

    %if ismember('relative',trans_types),
    if bitand(8,trans_types),
       resRel = HarmonicWindowCore(modX,w,num_harmonic_curr);
    end

    % all offsets for this radius
    windowis = find(windowi2radiusi == radiusi);
    for windowi = windowis',
      off = windows(windowi,2);
      res1 = padgrab2(res,nan,1,1,1+off,N+off);

      %if ismember('none',trans_types),
      if bitand(1,trans_types),
        y(end+1,:) = res1; %#ok<*AGROW>
        feature_names{end+1} = {'stat','harmonic','trans','none','radius',r,'offset',off,'num_harmonic',num_harmonic_curr};
      end

      %if ismember('abs',trans_types),
      if bitand(2,trans_types),
        y(end+1,:) = abs(res1);
        feature_names{end+1} = {'stat','harmonic','trans','abs','radius',r,'offset',off,'num_harmonic',num_harmonic_curr};
      end

      %if ismember('relative',trans_types),
      if bitand(8,trans_types),
        resRel1 = padgrab2(resRel,nan,1,1,1+off,N+off);
        y(end+1,:) = resRel1;
        feature_names{end+1} = {'stat','harmonic','trans','relative','radius',r,'offset',off,'num_harmonic',num_harmonic_curr};
      end

      if SANITY_CHECK,
        extraStr = sprintf('num_harmonic = %d ',num_harmonic_curr);
        funcType = 'harmonic';

        %if ismember('none',trans_types),
        if bitand(1,trans_types),
          fastY = res1; %#ok<*AGROW>
          res_dumb = nan(1,N);
          for n_dumb = 1:N,
            r_dumb = min([r,n_dumb+off-1,N-n_dumb-off]);
            if r_dumb < 1,
              continue;
            end
            w_dumb = 2*r_dumb+1;
            fil_dumb = cos(linspace(0,pi*num_harmonic_curr,w_dumb))/w_dumb*(num_harmonic_curr+1);
            res_dumb(n_dumb) = sum(fil_dumb.*x(n_dumb+off-r_dumb:n_dumb+off+r_dumb));
          end
          checkSanity(fastY,res1,r,off,funcType,'none',extraStr);
        end

        %if ismember('abs',trans_types),
        if bitand(2,trans_types),
          fastY = abs(res1);
          res_dumb = nan(1,N);
          for n_dumb = 1:N,
            r_dumb = min([r,n_dumb+off-1,N-n_dumb-off]);
            if r_dumb < 1,
              continue;
            end
            w_dumb = 2*r_dumb+1;
            fil_dumb = cos(linspace(0,pi*num_harmonic_curr,w_dumb))/w_dumb*(num_harmonic_curr+1);
            res_dumb(n_dumb) = abs(sum(fil_dumb.*x(n_dumb+off-r_dumb:n_dumb+off+r_dumb)));
          end
          checkSanity(fastY,res1,r,off,funcType,'abs',extraStr);
        end
      end

    end
  end

end

author: kristinbranson

repo: JAABA


text:
function timer_tictoc_test ( )

%*****************************************************************************80
%
%% TIMER_TICTOC_TEST uses TIC and TOC as the timer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2006
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TIMER_TICTOC\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate the use of the TIC routine.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  TIC and TOC are built-in MATLAB routines.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To measure and print elapsed wallclock time:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    tic;\n' );
  fprintf ( 1, '    THINGS TO BE TIMED;\n' );
  fprintf ( 1, '    toc;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To store elapsed wallclock time:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    tic;\n' );
  fprintf ( 1, '    THINGS TO BE TIMED;\n' );
  fprintf ( 1, '    wallclock = toc;\n' );

  timer_tictoc_test01 ( );
  timer_tictoc_test02 ( );
  timer_tictoc_test03 ( );
  timer_tictoc_test04 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TIMER_TICTOC\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function ecgpuwave(varargin)
%
% ecgpuwave(recordName,annFileName,startTime,stopTime,qrsAnn,pflag,signalList)
%
%    Wrapper to the ECGPUWAVE binary written by Pablo Laguna (laguna@posta.unizar.es), Raimon Jan[�e], Eudald Bogatell, and David Vigo Anglada:
%         http://www.physionet.org/physiotools/wag/ecgpuw-1.htm
%
% ECGPUWAVE analyses an ECG signal from the specified record 'recordName', detecting the QRS complexes
% and locating the beginning, peak, and end of the P, QRS, and ST-T waveforms. The output of ecgpuwave
% is written as a standard WFDB-format annotation file associated
% with the specified annotator 'outputAnnFileName'.
%
% The QRS detector is based on the algorithm of Pan and Tompkins (reference 1) with some improvements that
% make use of slope information (reference 2). Optionally, QRS annotations can be provided as input (see option 'qrsAnn'),
% permitting the use of external QRS detectors such as sqrs(1) or manually-edited annotations
% (which can be created using wave(1) ). The waveform limit locator is based on the algorithm
% described in reference 3 and evaluated in references 3 and 4.
%
% The output annotation file contains PWAVE ("p") and TWAVE ("t") annotations
% that indicate the P- and T-wave peaks, as well as QRS annotations (NORMAL ("N") if generated
% by the built-in QRS detector, or copies of the input QRS annotations if these were supplied).
% ECGPUWAVE classifies each T wave as type 0 (normal), 1 (inverted), 2 (positive monophasic),
% 3 (negative monophasic), 4 (biphasic negative-positive), or 5 (biphasic positive-negative);
% this numeric classification is written into the num field of each TWAVE annotation.
% The P, QRS, and T waveform onsets and ends are marked in the output annotation file using
% WFON ("(") and WFOFF (")") annotations. The num field of each WFON and WFOFF annotation designates
% the type of waveform with which it is associated: 0 for a P wave, 1 for a QRS complex, or 2 for a T wave.
%
% WARNING:
%   If ECGPUWAVE is used without providing an annotator file for the R peaks it relies on a Pan-Tompkins detector.
%   This detector may not be optimized for your particular signal, and thus may need to be tuned, yielding anempty annotation
%   file. In this case, you may want to first run a QRS detector to generate an appropriate QRS annotation file,
%   which you can then check manually prior to calling ECGPUWAVE with the
%   'qrsAnn' option described below.
%
% IMPORTANT NOTE:
%     A patch to the original Fortran source files was provided by Roberto Sassi in
%     order to compile ECGPUWAVE in 64-bit arch. This patch requires
%     compiling the source code with gfortran (g77 is not supported in gcc 4.x), which can
%     yield slighly different results. For more information please see:
%           http://www.dti.unimi.it/~sassi/software/wfdb64HowTo.htm
%
%
%
% References:
%1. Pan J and Tompkins WJ. A Real-Time QRS Detection Algorithm. IEEE Transactions on Biomedical Engineering 32(3):230-236, 1985.
%2. Laguna P. New Electrocardiographic Signal Processing Techniques: Application to Long-term Records. Ph. D. dissertation, Science Faculty, University of Zaragoza, 1990.
%3. Laguna P, Jan[�e] R, Caminal P. Automatic Detection of Wave Boundaries in Multilead ECG Signals: Validation with the CSE Database. Computers and Biomedical Research 27(1):45-60, 1994.
%4. Jan[�e] R, Blasi A, Garc[�i]a J, and Laguna P. Evaluation of an automatic threshold based detector of waveform limits in Holter ECG with the QT database. Computers in Cardiology 24:295-298 (1997; available at http://www.physionet.org/physiobank/database/qtdb/eval/ )
%
% Required Parameters:
%
% recordName
%       String specifying the name of the record to annotate with ECGPUWAVE.
%
% annFileName
%       String specifying the output file name that ECPUWAVE will store its
%       output.
%
% ecgpuwave(recordName,outputAnnFileName,startTime,stopTime,qrsAnn,pflag,signalList)
%
% Optional Parameters are:
%
% startTime
%      Either WFDB Time String or 1x1 integer specifying the sample number
%      specifying the starting location.
%
% stopTime
%      Either WFDB Time String or 1x1 integer specifying the sample number
%      specifying the ending location.
%
% qrsAnn
%       A String represengint the QRS annotation file that ECGPUWAVE should
%       use (must be in the WFDP PATH or within the current directory).
%
% pflag
%      A boolean (default = 0), if true all beats are processed, if false,
%      only normal beats are used.
%
% signalList
%     A Nx1 vector of integers specifying which signals in 'recordName' are
%     to be analyzed (default = [0].
%
% Wrapper written by Ikaro Silva, 2013
% Last Modified: May 29 , 2014
% Version 0.0.1
%
%
% %Example - Will go into the WFDB directory of examples to load data
%curdir=pwd;
%[~,config]=wfdbloadlib;eval(['cd ' config.WFDB_JAVA_HOME filesep 'example'])
%ecgpuwave('100s','test');
%[signal,Fs,tm]=rdsamp('100s');
%pwaves=rdann('100s','test',[],[],[],'p');
%plot(tm,signal(:,1));hold on;grid on
%plot(tm(pwaves),signal(pwaves),'or')
%cd(curdir) %CD back to current directory
%
%
% Since 0.9.5
%
%
% See also WFDBDESC, PHYSIONETDB, RDANN, WRANN, SQRS, WQRS, WFDBEXEC

%endOfHelp

javaWfdbExec=getWfdbClass('ecgpuwave');

%Set default pararamter values
inputs={'recordName','annFileName','startTime','stopTime','qrsAnn','pflag','signalList'};
startTime=[];
stopTime=[];
qrsAnn=[];
pflag=[];
signalList=[];
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

wfdb_argument={'-r',recordName,'-a',annFileName};

if(~isempty(startTime))
    wfdb_argument{end+1}='-f';
    wfdb_argument{end+1}=num2str(startTime-1);
end
if(~isempty(stopTime))
    wfdb_argument{end+1}='-t';
    wfdb_argument{end+1}=num2str(stopTime-1);
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

if(~isempty(qrsAnn))
    wfdb_argument{end+1}='-i';
    wfdb_argument{end+1}=qrsAnn;
    if(~isempty(pflag))
        wfdb_argument{end+1}='-n';
        wfdb_argument{end+1}=[num2str(pflag)];
    end
end

javaWfdbExec.execToStringList(wfdb_argument);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function [ Temp,e,Geopot,P,longrid,latgrid] = aps_weather_model_nan_check( Temp,e,Geopot,P,longrid,latgrid);
% Function which fills in the gaps in the weather model data with the
% nearest neighbor value, or from the level above.
%
%     Copyright (C) 2016  Bekaert David
%     davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By David Bekaert
% April 2016

% getting the pressure
Pressure_level = squeeze(P(1,1,:));
lon = double(longrid(:,:,1));
lat = double(latgrid(:,:,1));

% step_levels
% define it in such a way we loop from the upper atmopshere to the lower atmosphere
step_level = [length(Pressure_level):-1:1]';
if Pressure_level(1)<Pressure_level(end) % first level is top APS
    step_level = [1:length(Pressure_level)]';
end

% looping though and use TEMP as inital check for nan
n_pixels = size(P,1)*size(P,2);
for k=1:length(Pressure_level)
    ix_nan = isnan(Temp(:,:,step_level(k)));
    n_nan = sum(sum(ix_nan));

    % check if there are any NaN values
    if n_nan<n_pixels-3
        temp_data = Temp(:,:,step_level(k));
        Temp(:,:,step_level(k)) = griddata(reshape(lon(~ix_nan),[],1),reshape(lat(~ix_nan),[],1),double(reshape(temp_data(~ix_nan),[],1)),lon,lat,'nearest');

        temp_data = e(:,:,step_level(k));
        e(:,:,step_level(k)) = griddata(reshape(lon(~ix_nan),[],1),reshape(lat(~ix_nan),[],1),double(reshape(temp_data(~ix_nan),[],1)),lon,lat,'nearest');

    % check if all are NaN values
    elseif n_nan==n_pixels
        if k==1
            error('Weird seems like the top of the atmopshere has no data')
        end
        Temp(:,:,step_level(k))=Temp(:,:,step_level(k-1));
        e(:,:,step_level(k))=e(:,:,step_level(k-1));
    else
        Temp(:,:,step_level(k))=nanmean(nanmean(Temp(:,:,step_level(k))));
        e(:,:,step_level(k))=nanmean(nanmean(e(:,:,step_level(k))));

    end
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function [Coh] = SLC_cov(slcstack,SHP)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file is part of TomoSAR.
%
%   TomoSAR is distributed in the hope that it will be useful,
%   but without warranty of any kind; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%   See the Apache License for more details.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author : Dinh Ho Tong Minh (INRAE) and Yen Nhi Ngo, Jan. 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nlines,nwidths,n_slc]=size(slcstack);

% Normalize
slcstack(slcstack~=0) = slcstack(slcstack~=0)./abs(slcstack(slcstack~=0));

CalWin =SHP.CalWin;
RadiusRow=(CalWin(1)-1)/2;
RadiusCol=(CalWin(2)-1)/2;

mlistack = abs(slcstack);

% Coherence matrix estimate
Coh=zeros(n_slc,n_slc,nlines,nwidths,'single');
for ii=1:n_slc
    m1_intial = mlistack(:,:,ii);
    for ss = ii+1:n_slc
        m2_intial = mlistack(:,:,ss);
        Dphi = exp(1i*(angle(slcstack(:,:,ii).*conj(slcstack(:,:,ss)))));
        Intf= sqrt(m1_intial.*m2_intial).*Dphi;

        % Padding at edge
        m1 = padarray(m1_intial,[RadiusRow RadiusCol],'symmetric');
        m2 = padarray(m2_intial,[RadiusRow RadiusCol],'symmetric');
        Intf= padarray(Intf,[RadiusRow RadiusCol],'symmetric');
        nu = zeros(nlines,nwidths,'single');
        de1=nu;
        de2=nu;
        num=1;
        for jj = 1:nwidths
            for kk= 1:nlines
                x_global  = jj+RadiusCol;
                y_global  = kk+RadiusRow;
                MasterValue= m1(y_global-RadiusRow:y_global+RadiusRow,x_global-RadiusCol:x_global+RadiusCol);
                SlaveValue = m2(y_global-RadiusRow:y_global+RadiusRow,x_global-RadiusCol:x_global+RadiusCol);
                InterfValue= Intf(y_global-RadiusRow:y_global+RadiusRow,x_global-RadiusCol:x_global+RadiusCol);
                MasterValue= MasterValue(SHP.PixelInd(:,num));
                SlaveValue = SlaveValue(SHP.PixelInd(:,num));
                InterfValue= InterfValue(SHP.PixelInd(:,num));
                nu(kk,jj)  = sum(InterfValue);
                de1(kk,jj) = sum(MasterValue);
                de2(kk,jj) = sum(SlaveValue);
                num=num+1;
            end
        end
        Coh(ii,ss,:,:) = nu./sqrt(de1.*de2);
    end
end

% Make mirror operator
temp = ones(1,n_slc) ;
for jj = 1:nwidths
    for kk= 1:nlines
        W = Coh(:, :,kk,jj) ;
        Coh(:, :,kk,jj) = W + (W - diag(temp))';
    end
end

return

author: DinhHoTongMinh

repo: TomoSAR


text:
function x=pcmu2lin(p,s)
%PCMU2LIN Convert Mu-law PCM to linear X=(P,S)
%	lin = pcmu2lin(pcmu) where pcmu contains a vector
%	of mu-law values in the range 0 to 255.
%	No checking is performed to see that numbers are in this range.
%
%	Output values are divided by the scale factor s:
%
%		   s		Output Range
%
%		   1		+-8031	(integer values)
%		4004.2	+-2.005649 (default)
%		8031		+-1
%		8159		+-0.9843118 (+-1 nominal full scale)
%
%	The default scaling factor 4004.189931 is equal to
%	sqrt((2207^2 + 5215^2)/2) this follows ITU standard G.711.
%	The sine wave with PCM-Mu values [158 139 139 158 30 11 11 30]
%	has a mean square value of unity corresponding to 0 dBm0.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: pcmu2lin.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2
  t=9.98953613E-4;
else
  t=4/s;
end

m=15-rem(p,16);
q=floor(p/128);
e=(127-p-m+128*q)/16;
x=(q-0.5).*(pow2(m+16.5,e)-16.5)*t;

author: decouples


% pivot_setup.m // Jon Lee
% set up a pivoting example

clear all;

global A b c m n beta eta

try
    reply = input('NAME of NAME.m file holding input? [pivot_input]: ', 's');
    if isempty(reply)
       reply = 'pivot_input';
       end
catch
    reply = 'pivot_input'; % need catch for execution on MathWorks Cloud
end
eval(reply);

if (size(b) ~= m)
    display('size(b) does not match number of rows of A')
    return
end
if (size(c) ~= n)
    display('size(c) does not match number of columns of A')
    return
end
if(size(setdiff(beta,1:n)) > 0)
    display('beta has elements not in 1,2,...,n')
    return
end
if(size(setdiff(eta,1:n)) > 0)
    display('eta has elements not in 1,2,...,n')
    return
end
if (size(beta) ~= m)
    display('size(beta) does not match number of rows of A')
    return
end
if (size(eta) ~= n-m)
    display('size(eta) does not match number of cols minus number of rows of A')
    return
end

display('Available quantites: A, b, c, m, n, beta, eta');
display('Seems like a good time to run pivot_algebra');


text:
%TROTY Rotation about Y axis
%
% T = TROTY(THETA) is a homogeneous transformation (4x4) representing a rotation
% of THETA radians about the y-axis.
%
% T = TROTY(THETA, 'deg') as above but THETA is in degrees.
%
% Notes::
% - Translational component is zero.
%
% See also ROTY, TROTX, TROTZ, TROT2.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function T = troty(t, varargin)
	T =    [roty(t, varargin{:}) [0 0 0]'; 0 0 0 1];

author: Allopart

repo: rbpf-gmapping


text:
function [MriFileReg, errMsg, fileTag, sMriReg] = mri_coregister(MriFileSrc, MriFileRef, Method, isReslice, isAtlas)
% MRI_COREGISTER: Compute the linear transformations on both input volumes, then register the first on the second.
%
% USAGE:  [MriFileReg, errMsg, fileTag, sMriReg] = mri_coregister(MriFileSrc, MriFileRef, Method, isReslice)
%            [sMriReg, errMsg, fileTag, sMriReg] = mri_coregister(sMriSrc,    sMriRef, ...)
%
% INPUTS:
%    - MriFileSrc : Relative path to the Brainstorm MRI file to register
%    - MriFileRef : Relative path to the Brainstorm MRI file used as a reference
%    - sMriSrc    : Brainstorm MRI structure to register (fields Cube, Voxsize, SCS, NCS...)
%    - sMriRef    : Brainstorm MRI structure used as a reference
%    - Method     : Method used for the coregistration of the volume: 'spm', 'mni', 'vox2ras'
%    - isReslice  : If 1, reslice the output volume to match dimensions of the reference volume
%    - isAtlas    : If 1, perform only integer/nearest neighbors interpolations (MNI and VOX2RAS registration only)
%
% OUTPUTS:
%    - MriFileReg : Relative path to the new Brainstorm MRI file (containing the structure sMriReg)
%    - errMsg     : Error messages if any
%    - fileTag    : Tag added to the comment/filename
%    - sMriReg    : Brainstorm MRI structure with the registered volume

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2016-2023

% ===== LOAD INPUTS =====
% Parse inputs
if (nargin < 5) || isempty(isAtlas)
    isAtlas = 0;
end
% Initialize returned variables
MriFileReg = [];
errMsg = [];
fileTag = '';
sMriReg = [];
% Progress bar
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'MRI register', 'Loading input volumes...');
end
% USAGE: mri_coregister(sMriSrc, sMriRef, ...)
if isstruct(MriFileSrc)
    sMriSrc = MriFileSrc;
    sMriRef = MriFileRef;
    MriFileSrc = [];
    MriFileRef = [];
% USAGE: mri_coregister(MriFileSrc, MriFileRef, ...)
elseif ischar(MriFileSrc)
    % Get the default MRI for this subject
    if isempty(MriFileRef)
        sSubject = bst_get('MriFile', MriFileSrc);
        MriFileRef = sSubject.Anatomy(sSubject.iAnatomy).FileName;
    end
    sMriSrc = in_mri_bst(MriFileSrc);
    sMriRef = in_mri_bst(MriFileRef);
else
    error('Invalid call.');
end
% % Not available for multiple volumes
% if (size(sMriRef.Cube, 4) > 1) || (size(sMriSrc.Cube, 4) > 1)
%     errMsg = 'The input files cannot contain multiple volumes.';
%     return;
% end
% Inialize various variables
isUpdateScs = 0;
isUpdateNcs = 0;

% ===== REGISTER VOLUMES =====
switch lower(Method)

    % ===== METHOD: SPM =====
    case 'spm'
        % Initialize SPM
        [isInstalled, errMsg] = bst_plugin('Install', 'spm12');
        if ~isInstalled
            if ~isProgress
                bst_progress('stop');
            end
            return;
        end
        bst_plugin('SetProgressLogo', 'spm12');

        % === SAVE FILES IN TMP FOLDER ===
        bst_progress('text', 'Saving temporary files...');
        % Get temporary folder
        TmpDir = bst_get('BrainstormTmpDir', 0, 'spmcoreg');
        % Save source MRI in .nii format
        NiiSrcFile = bst_fullfile(TmpDir, 'spm_src.nii');
        out_mri_nii(sMriSrc, NiiSrcFile);
        % Save reference MRI in .nii format
        NiiRefFile = bst_fullfile(TmpDir, 'spm_ref.nii');
        out_mri_nii(sMriRef, NiiRefFile);

        % === CALL SPM COREGISTRATION ===
        bst_progress('text', 'Calling SPM batch...');
        % Code initially coming from Olivier David's ImaGIN_anat_spm.m function
        % Initial translation according to centroids
        % Reference volume
        Vref = spm_vol([NiiRefFile, ',1']);
        [Iref,XYZref] = spm_read_vols(Vref);
        Iindex = find(Iref>max(Iref(:))/6);
        Zindex = find(max(XYZref(3,:))-XYZref(3,:)<200);
        index = intersect(Iindex,Zindex);
        CentroidRef = mean(XYZref(:,index),2);
        % Volume to register
        V2 = spm_vol([NiiSrcFile, ',1']);
        [I2,XYZ2] = spm_read_vols(V2);
        Iindex = find(I2>max(I2(:))/6);
        Zindex = find(max(XYZ2(3,:))-XYZ2(3,:)<200);
        index = intersect(Iindex,Zindex);
        Centroid2 = mean(XYZ2(:,index),2);
        % Apply translation
        B = [CentroidRef'-Centroid2' 0 0 0 1 1 1 0 0 0];
        M = spm_matrix(B);
        Mat = spm_get_space(V2.fname);
        spm_get_space(V2.fname, M*Mat);

        % Create coregistration batch
        if isReslice
            % Coreg: Estimate and reslice
            matlabbatch{1}.spm.spatial.coreg.estwrite.ref      = {[NiiRefFile, ',1']};
            matlabbatch{1}.spm.spatial.coreg.estwrite.source   = {[NiiSrcFile, ',1']};
            matlabbatch{1}.spm.spatial.coreg.estwrite.other    = {''};
            matlabbatch{1}.spm.spatial.coreg.estwrite.eoptions = spm_get_defaults('coreg.estimate');
            matlabbatch{1}.spm.spatial.coreg.estwrite.woptions = spm_get_defaults('coreg.write');
            matlabbatch{1}.spm.spatial.coreg.estwrite.woptions.outdir = TmpDir;
            % Output file
            NiiRegFile = bst_fullfile(TmpDir, 'rspm_src.nii');
        else
            % Coreg: Estimate
            matlabbatch{1}.spm.spatial.coreg.estimate.ref      = {[NiiRefFile, ',1']};
            matlabbatch{1}.spm.spatial.coreg.estimate.source   = {[NiiSrcFile, ',1']};
            matlabbatch{1}.spm.spatial.coreg.estimate.other    = {''};
            matlabbatch{1}.spm.spatial.coreg.estimate.eoptions = spm_get_defaults('coreg.estimate');
            % Output file
            NiiRegFile = NiiSrcFile;
        end
        % Run SPM batch
        spm_jobman('initcfg');
        % spm_jobman('interactive', matlabbatch)
        spm_jobman('run',matlabbatch)
        % Read output volume
        sMriReg = in_mri(NiiRegFile, 'ALL', 0, 0);
        % If an error occurred in SPM
        if isempty(sMriReg)
            errMsg = 'An unknown error occurred while executing SPM. See the logs in the command window.';
            if ~isProgress
                bst_progress('stop');
            end
            return;
        end
        % Delete the temporary files
        file_delete(TmpDir, 1, 1);
        % Output file tag
        fileTag = '_spm';
        % Remove logo
        bst_plugin('SetProgressLogo', []);

        % === UPDATE FIDUCIALS ===
        if isReslice
            % Use the reference SCS coordinates
            if isfield(sMriRef, 'SCS')
                sMriReg.SCS = sMriRef.SCS;
            end
            % Use the reference NCS coordinates
            if isfield(sMriRef, 'NCS')
                sMriReg.NCS = sMriRef.NCS;
            end
        else
            isUpdateScs = 1;
            isUpdateNcs = 1;
        end

        % ===== COMPUTE TRANSFORMATION MATRIX =====
        % Get transformations MRI=>WORLD (vox2ras) for original volume
        vox2ras_src = cs_convert(sMriSrc, 'mri', 'world');
        vox2ras_reg = cs_convert(sMriReg, 'mri', 'world');
        % If there are vox2ras transformation matrices available
        if ~isempty(vox2ras_src) && ~isempty(vox2ras_reg)
            % Compute the transformation from the original to the registered volume
            TransfRegSpm = vox2ras_reg * inv(vox2ras_src);
            % Save in the registered MRI structure (this transformation must be applied on world coordinates)
            sMriReg.InitTransf(end+1,[1 2]) = {'reg', TransfRegSpm};
        end

    % ===== METHOD: MNI =====
    case 'mni'
        % === COMPUTE MNI TRANSFORMATIONS ===
        % Source MRI
        if ~isfield(sMriSrc, 'NCS') || ~isfield(sMriSrc.NCS, 'R') || ~isfield(sMriSrc.NCS, 'T') || isempty(sMriSrc.NCS.R) || isempty(sMriSrc.NCS.T)
            [sMriSrc,errMsg] = bst_normalize_mni(sMriSrc, 'maff8');
        end
        % Reference MRI
        if ~isfield(sMriRef, 'NCS') || ~isfield(sMriRef.NCS, 'R') || ~isfield(sMriRef.NCS, 'T') || isempty(sMriRef.NCS.R) || isempty(sMriRef.NCS.T)
            [sMriRef,errMsg] = bst_normalize_mni(sMriRef, 'maff8');
        end
        % Handle errors
        if ~isempty(errMsg)
            if ~isProgress
                bst_progress('stop');
            end
            return;
        end
        % Get MNI transformations
        TransfSrc = [sMriSrc.NCS.R, sMriSrc.NCS.T; 0 0 0 1];
        TransfRef = [sMriRef.NCS.R, sMriRef.NCS.T; 0 0 0 1];

        % === RESLICE VOLUME ===
        if isReslice
            % Reslice the volume
            [sMriReg, errMsg] = mri_reslice(sMriSrc, sMriRef, TransfSrc, TransfRef, isAtlas);
        else
            % Save the original input volume
            sMriReg = sMriSrc;
            isUpdateScs = 1;
            isUpdateNcs = 0;
        end
        % Output file tag
        fileTag = '_mni';

    % ===== VOX2RAS =====
    case 'vox2ras'
        % Nothing to do, just reslice if needed
        if isReslice
            % Reslice the volume
            [sMriReg, errMsg] = mri_reslice(sMriSrc, sMriRef, 'vox2ras', 'vox2ras', isAtlas);
            % Output file tag
            fileTag = '_reg';
        else
            % Save the original input volume
            sMriReg = sMriSrc;
            isUpdateScs = 1;
            isUpdateNcs = 1;
            % Output file tag
            fileTag = '';
        end
end
% Handle errors
if ~isempty(errMsg)
    if ~isProgress
        bst_progress('stop');
    end
    return;
end

% ===== REMOVE NEW NEGATIVE VALUES =====
% If some negative values appeared just because of the registration/reslicing: remove them
if any(sMriReg.Cube(:) < 0) && ~any(sMriSrc.Cube(:) < 0)
    sMriReg.Cube(sMriReg.Cube < 0) = 0;
end

% ===== UPDATE FIDUCIALS =====
if isUpdateScs || isUpdateNcs
    % Get vox2ras transformations for all volumes
    mri2world_reg = cs_convert(sMriReg, 'mri', 'world');
    mri2world_ref = cs_convert(sMriRef, 'mri', 'world');
    % Transfer fiducials from reference MRI to registered MRI
    if ~isempty(mri2world_reg) && ~isempty(mri2world_ref)
        % Apply transformation: reference MRI => SPM RAS/world => registered MRI
        Transf = inv(mri2world_reg) * (mri2world_ref);
        % Convert to millimeters, just like the SCS and NCS transformation matrices
        Transf(1:3,4) = Transf(1:3,4) .* 1000;
        % SCS coordinates
        if isUpdateScs && isfield(sMriRef, 'SCS')
            SCS = sMriRef.SCS;
            % Update fiducials coordinates
            if all(isfield(SCS, {'NAS','LPA','RPA'})) && ~isempty(SCS.NAS) && ~isempty(SCS.LPA) && ~isempty(SCS.RPA)
                NAS = (Transf * [SCS.NAS, 1]')';
                LPA = (Transf * [SCS.LPA, 1]')';
                RPA = (Transf * [SCS.RPA, 1]')';
                sMriReg.SCS.NAS = NAS(1:3);
                sMriReg.SCS.LPA = LPA(1:3);
                sMriReg.SCS.RPA = RPA(1:3);
            end
            % Update SCS transformation
            if all(isfield(SCS, {'T','R'})) && ~isempty(SCS.R) && ~isempty(SCS.T)
                Tscs = [SCS.R, SCS.T; 0 0 0 1] * inv(Transf);
                sMriReg.SCS.R = Tscs(1:3,1:3);
                sMriReg.SCS.T = Tscs(1:3,4);
            end
        end
        % NCS coordinates
        if isUpdateNcs && isfield(sMriRef, 'NCS')
            NCS = sMriRef.NCS;
            % Update NCS fiducials
            if all(isfield(NCS, {'AC','PC','IH'})) && ~isempty(NCS.AC) && ~isempty(NCS.PC) && ~isempty(NCS.IH)
                AC = (Transf * [NCS.AC, 1]')';
                PC = (Transf * [NCS.PC, 1]')';
                IH = (Transf * [NCS.IH, 1]')';
                sMriReg.NCS.AC = AC(1:3);
                sMriReg.NCS.PC = PC(1:3);
                sMriReg.NCS.IH = IH(1:3);
            end
            % Update NCS transformation
            if all(isfield(NCS, {'T','R'})) && ~isempty(NCS.R) && ~isempty(NCS.T)
                Tncs = [NCS.R, NCS.T; 0 0 0 1] * inv(Transf);
                sMriReg.NCS.R = Tncs(1:3,1:3);
                sMriReg.NCS.T = Tncs(1:3,4);
            end
        end
    end
end

% ===== SAVE NEW FILE =====
% Add file tag
if isReslice
    fileTag = [fileTag, '_reslice'];
end
sMriReg.Comment = [sMriSrc.Comment, fileTag];
% Save output
if ~isempty(MriFileSrc)
    bst_progress('text', 'Saving new file...');
    % Get subject
    [sSubject, iSubject, iMri] = bst_get('MriFile', MriFileSrc);
    % Update comment
    sMriReg.Comment = file_unique(sMriReg.Comment, {sSubject.Anatomy.Comment});
    % Add history entry
    sMriReg.History = sMriSrc.History;
    sMriReg = bst_history('add', sMriReg, 'resample', ['MRI co-registered on default file (' Method '): ' MriFileRef]);
    % Save new file
    MriFileRegFull = file_unique(strrep(file_fullpath(MriFileSrc), '.mat', [fileTag '.mat']));
    MriFileReg = file_short(MriFileRegFull);
    % Save new MRI in Brainstorm format
    sMriReg = out_mri_bst(sMriReg, MriFileRegFull);

    % Register new MRI
    iAnatomy = length(sSubject.Anatomy) + 1;
    sSubject.Anatomy(iAnatomy) = db_template('Anatomy');
    sSubject.Anatomy(iAnatomy).FileName = MriFileReg;
    sSubject.Anatomy(iAnatomy).Comment  = sMriReg.Comment;
    % Update subject structure
    bst_set('Subject', iSubject, sSubject);
    % Refresh tree
    panel_protocols('UpdateNode', 'Subject', iSubject);
    panel_protocols('SelectNode', [], 'anatomy', iSubject, iAnatomy);
    % Save database
    db_save();
else
    % Return output structure
    MriFileReg = sMriReg;
end
% Close progress bar
if ~isProgress
    bst_progress('stop');
end

author: brainstorm-tools


text:
function obj = c_func_hand_eye_new(in1,in2,in3,in4)
coef_f1_q_sym1 = in2(:,1);
coef_f2_q_sym1 = in3(:,1);
coef_f3_q_sym1 = in4(:,1);
coef_f1_q_sym11 = in2(:,11);
coef_f2_q_sym11 = in3(:,11);
coef_f3_q_sym11 = in4(:,11);
obj = [-coef_f1_q_sym1-coef_f1_q_sym11;-coef_f2_q_sym1-coef_f2_q_sym11;-coef_f3_q_sym1-coef_f3_q_sym11];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function scan_data = p10rd(ufilename,plotmode,verbose)
%P10RD Read Tencor P-10 surface profilometer data file
%
% SCAN_DATA = P10RD(UFILENAME,PLOTMODE,VERBOSE)
%
%P10RD reads the data from file UFILENAME.txt, and returns a 5-column
% array SCAN_DATA:
%  (:,1) is the horizontal position data in microns
%  (:,2) is the raw vertical profile data in microns
%  (:,3) is the levelled vertical profile data in microns
%  (:,4) is the Wav/A data (unscaled)
%  (:,5) is the Rough/A data (unscaled)
%
% Options:
% If VERBOSE == 0, no messages will be displayed.
% If VERBOSE == 1, the data summary will be displayed. (default)
% If VERBOSE == 2, the file header and data summary will be displayed.
% If PLOTMODE == 0, the data will not be displayed.
% if PLOTMODE == 1, the levelled data will be displayed in a new figure. (default)
% If PLOTMODE == 2, a second data plot with the raw data will be shown.
%
% If no input parameters are specified, the program is interactive, and the
% user will be prompted to enter a filename.
%
% NOTE: The function can only read ASCII text data files, not binary data files.
%  The name of the data file must have the extension ".txt", but you do not
%  need include the extension when specifying the file name.
%
% Examples:
%
%  If the data is stored in a text file called "mydata.txt", then:
%
%      p10rd mydata;
%
%  will load the file, plot the data, and print a summary. You can also use:
%
%      p10rd;
%
%  and you will be prompted for the file name. The command:
%
%      profile_data=p10rd('mydata',0,1);
%
%  will open the file mydata.txt, plot the data,
%   and save the data to the variable profile_data in the workspace
%   without printing anything in the command window.
%
%
% M.A. Hopcroft
%  mhopeng at gmail dot com
%
% The license governing use of this code is in the accompanying file
%  "license.txt".
%
versionstr='p10rd 1.3';

% MH Jul2013
% v1.3  add uifileget
% MH Sep2012
% v1.2  updates for modern matlab
% MH MAY2007
% v1.1  re-order input (verbose last)
%       filename w/ or w/o .txt
% MH OCT2006
% v1.0 out of beta
%
%#ok<*NBRAK>

% handle input arguments - set defaults
% default to verbose mode 1 (summary only)
if nargin < 3
    verbose = 1;
end

% default to plot mode 1 (levelled data only)
if nargin < 2
    plotmode = 1;
end

% get interactive if not being called by another program
if (nargout == 0)
    if (verbose >= 1 && nargin == 0)
        clc;
        fprintf(1, '\n\n  %s\n  %s\n\n', versionstr, 'Read in Tencor P-10 data files');
    end
end

%%%%
% open user's file
if nargin >= 1 % if the user already specified a file, use it
    userfilename=ufilename;

    % check to see if file exists
    cfid=fopen(userfilename);
    if cfid== -1
        % try with ".txt"
        userfilename=strcat(userfilename,'.txt');
        cfid=fopen(userfilename);
        if cfid== -1
            disp(pwd);
            fprintf(1,'ERROR: data file ["%s"] not found in current directory',userfilename);
            scan_data=[];
            return
        end
        fclose(cfid);
    end

else % if user did not specify the file name, ask for one
    [ufilename, ufilepath]=uigetfile('.txt','Select the Tencor P-10 data file');
    if isequal(ufilename,0)
        scan_data=[];
        return
    end
    userfilename=fullfile(ufilepath,ufilename);
    %userfilename=input(' Enter the name of the Dektak XT .csv file: ','s');
end

% open the file
userfile=fopen(userfilename);
[fpath, fname, fext]=fileparts(userfilename);

%%%%
% process the file header
% the first 5 lines are header information
% read them in and get some data: scan length, number of points, etc.
if verbose == 2
    fprintf(1,'\n\n %s\n %s %s\n\n %s\n\n', versionstr, 'Filename:', fname, 'File header:');
end

% additional linein statements should handle unix/DOS ASCII/ANSI problem
%  [note: investigate this further]

% line 1: filename and date
linein=fgetl(userfile);
if isempty(linein), linein=fgetl(userfile); end % spurious carriage returns?
if verbose == 2
    fprintf(1,' %s\n', linein);
end
scan_date=sscanf(linein,'* %*s %s',1);

% line 2: number of data points
linein=fgetl(userfile);
if isempty(linein), linein=fgetl(userfile); end % spurious carriage returns?
if verbose == 2
    fprintf(1,' %s\n', linein);
end
% %  because the bastards use tabs instead of spaces, we have to
% %   parse the line to isolate the number at the start of it
% lspace=0; % initialize
% for i=1:8
%     if isspace(linein(i))
%         lspace=i;
%         break
%     end
% end
% num_points=str2num(linein(1:lspace-1));
num_points=sscanf(linein,'%d',1);

% line 3: length of scan (um)
linein=fgetl(userfile); % line 3: length of scan (um)
if isempty(linein), linein=fgetl(userfile); end % spurious carriage returns?
if verbose == 2
    fprintf(1,' %s\n', linein);
end
scan_length=sscanf(linein,'%d',1);

%line 4: "A per Y count"
% (I believe this represents the minimum vertical resolution for the scan)
linein=fgetl(userfile);
if isempty(linein), linein=fgetl(userfile); end %spurious carriage returns?
if verbose == 2
    fprintf(1,' %s\n', linein);
end

%line 5: column headers
linein=fgetl(userfile);
if isempty(linein), linein=fgetl(userfile); end
if verbose == 2
    fprintf(1,' %s\n', linein);
end

%the position of the data points is now known
% the P-10 appears to add spurious data points, i.e. 700um scan, 702 data points
%  ignore step
% stepsize=scan_length/num_points;
%rstep=round(stepsize);
%if rstep-stepsize <= 0.005
%    stepsize=rstep; %the scan length is effectively lengthed slightly
%end

%%%%
% get data from file
%  data is formatted in 5 columns
%  column 1 is horizontal position in um
%  column 2 is the "raw data" in Angstroms
%  column 3 is the levelled data in Angstroms
%  column 4 is "Wav/A" ??
%  column 5 is "Rough/A" ?? surface roughness?

d=textscan(userfile,'%f %f %f %f %f',inf);
scan_data(:,1)=d{1};
scan_data(:,2)=d{2};
scan_data(:,3)=d{3};
scan_data(:,4)=d{4};
scan_data(:,5)=d{5};
count=size(scan_data,1);
sd=diff(scan_data(:,1));
if all(sd==sd(1))
    stepsize=sd(1);
else
    stepsize=scan_length/num_points;
end

% close the file
fclose(userfile);

%adjust data point number to be position
%scan_data(:,1)=scan_data(:,1).*stepsize;
%scale data to um from Angstroms
scan_data(:,2:5)=scan_data(:,2:5)./10000;

% plot
% plot the raw data
if plotmode >= 1

    figure;
    plot(scan_data(:,1),scan_data(:,2),'.-b');
    title(['p10rd: ' fname],'Interpreter','none','FontSize',16,'FontWeight','bold');
    xlabel('Scan Length (\mum)','FontSize',14);
    ylabel('Profile (\mum)','FontSize',14);
    grid on;

end

% plot the levelled data
if plotmode == 2

    hold on
    plot(scan_data(:,1),scan_data(:,3),'.-g');
    legend('Raw Data','Levelled Data');

end

%print summary data
if verbose >= 1

    fprintf(1,'\n');
    fprintf(1,' %s\n', 'Data Summary:');
    fprintf(1,'\n');
    fprintf(1,' %s %s\n', 'File name:', [fname fext]);
    if ~isempty(scan_date)
        fprintf(1,' %s %s\n','Scan Date:',scan_date);
    end

    fprintf(1,'\n');

    fprintf(1,' %s %i\n', 'Number of data points read:', count);
    fprintf(1,' %s %g %s\n', 'Scan Length:', scan_length, 'um');
    fprintf(1,' %s %d %s\n %s %d %s\n %s %g %s\n',...
        'Scan Start:', scan_data(1,1), 'um', 'Scan End:',...
        scan_data(end,1), 'um', 'Scan Step:', stepsize, 'um');

    fprintf(1,'\n');

end


text:
classdef TestUndistort
    %TestUndistort

    methods (Static)
        function test_1
            im = uint8([...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
            ]);
            result = cv.undistort(im, eye(3), []);
            validateattributes(result, {class(im)}, {'size',size(im)});
            assert(isequal(result, im));
        end

        function test_error_argnum
            try
                cv.undistort();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function bemcp_example
% Simple function to test/demonstrate how the Boundary element functions are
% used in combination with Fildtrip/Forwinv routines.
%
% 1. A model is created as 3 concentric meshed spheres (using FT's
% icosahedron routines),
% 2. then random electrodes are placed on the upper part of the outer
% sphere.
% 3. the model is then "prepared" with 'ft_prepare_bemmodel', this bits
% takes most time as it requires LOTS of calculation.
% 4. sensors and volumes are plugged together by 'forwinv_prepare_vol_sens'
% 5. Finally the leadfiled for 3 orthogonal sources placed at one location
% is calculated with 'forwinv_compute_leadfield.m'
% 6. Display the 3 leadfields
%
% NOTE:
% this bit of code needs access to low level fieldtrip/forwinv routines
% which have been copy/pasted here under.
% Be aware that this way of programming is generally NOT advisable!
% I used it only to ensure a quick & dirty check of the BEM module...

% Christophe Phillips
% $Id: bemcp_example.m 4793 2012-07-19 11:35:24Z spm $

% create volume conductor starting from unit sphere
[pnt, tri] = icosahedron162;

vol = [];
vol.cond = [1 1/80 1];
vol.source = 1; % index of source compartment
vol.skin_surface   = 3; % index of skin surface
% inner_skull_surface
vol.bnd(1).pnt = pnt*88;
vol.bnd(1).tri = tri;
% outer_skull_surface
vol.bnd(2).pnt = pnt*92;
vol.bnd(2).tri = tri;
% skin_surface
vol.bnd(3).pnt = pnt*100;
vol.bnd(3).tri = tri;

% create the BEM system matrix
cfg = [];
cfg.method = 'bemcp';
vol1 = ft_prepare_bemmodel(cfg, vol);

% create some random electrodes
pnt = randn(200,3);
pnt = pnt(pnt(:,3)>0, :);  % only those on the upper half
sens = [];
for i=1:size(pnt,1)
  sens.pnt(i,:) = pnt(i,:) / norm(pnt(i,:)); % scale towards the skin surface
  sens.label{i} = sprintf('%02d', i);
end

% prepare the sensor array and volume conduction, i.e. set up the linear
% interpolation from vertices to electrodes
[vol2, sens] = forwinv_prepare_vol_sens(vol1, sens);

lf = forwinv_compute_leadfield([0 0 50], sens, vol2);

figure; triplot(sens.pnt, [], lf(:,1)); colorbar
figure; triplot(sens.pnt, [], lf(:,2)); colorbar
figure; triplot(sens.pnt, [], lf(:,3)); colorbar

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Subfunctions from FieldTrip
function [pnt, dhk] = icosahedron162

% ICOSAHEDRON162 creates a 2-fold refined icosahedron

% Copyright (C) 2003, Robert Oostenveld
%
% $Log: icosahedron162.m,v $
% Revision 1.3  2003/11/28 09:40:12  roberto
% added a single help line
%
% Revision 1.2  2003/03/04 21:46:18  roberto
% added CVS log entry and synchronized all copyright labels
%

[pnt, dhk] = icosahedron;
[pnt, dhk] = refine(pnt, dhk);
[pnt, dhk] = refine(pnt, dhk);

pnt = pnt ./ repmat(sqrt(sum(pnt.^2,2)), 1,3);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pnt, dhk] = icosahedron

% ICOSAHEDRON creates an icosahedron
%
% [pnt, dhk] = icosahedron
% creates an icosahedron with 12 vertices and 20 triangles
%
% See also OCTAHEDRON, ICOSAHEDRON42, ICOSAHEDRON162, ICOSAHEDRON642, ICOSAHEDRON2562

% Copyright (C) 2002, Robert Oostenveld
%
% $Log: icosahedron.m,v $
% Revision 1.4  2006/07/26 11:03:38  roboos
% added "see also octahedron"
%
% Revision 1.3  2003/03/11 15:35:20  roberto
% converted all files from DOS to UNIX
%
% Revision 1.2  2003/03/04 21:46:18  roberto
% added CVS log entry and synchronized all copyright labels
%

dhk = [
   1   2   3
   1   3   4
   1   4   5
   1   5   6
   1   6   2
   2   8   3
   3   9   4
   4  10   5
   5  11   6
   6   7   2
   7   8   2
   8   9   3
   9  10   4
  10  11   5
  11   7   6
  12   8   7
  12   9   8
  12  10   9
  12  11  10
  12   7  11
];

pnt = zeros(12, 3);

rho=0.4*sqrt(5);
phi=2*pi*(0:4)/5;

pnt( 1, :) = [0 0  1];			% top point

pnt(2:6, 1) = rho*cos(phi)';
pnt(2:6, 2) = rho*sin(phi)';
pnt(2:6, 3) = rho/2;

pnt(7:11, 1) = rho*cos(phi - pi/5)';
pnt(7:11, 2) = rho*sin(phi - pi/5)';
pnt(7:11, 3) = -rho/2;

pnt(12, :) = [0 0 -1];			% bottom point

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pntr, dhkr] = refine(pnt, dhk, method, varargin)

% REFINE a 3D surface that is described by a triangulation
%
% Use as
%   [pnt, tri] = refine(pnt, tri)
%   [pnt, tri] = refine(pnt, tri, 'updown', numtri)
%
% The default method is to refine the mesh globally by inserting a vertex at
% each edge according to the algorithm described in Banks, 1983.
%
% The alternative 'updown' method refines the mesh a couple of times
% using Banks' algorithm, followed by a downsampling using the REDUCEPATCH
% function.

% The Banks method is a memory efficient implementation which remembers
% the previously inserted vertices. The refinement algorithm executes in
% linear time with the number of triangles.

% Copyright (C) 2002-2005, Robert Oostenveld
%
% $Log: refine.m,v $
% Revision 1.4  2005/05/06 08:54:31  roboos
% added 'updown method, using matlab reducepatch function
%
% Revision 1.3  2003/03/11 15:35:20  roberto
% converted all files from DOS to UNIX
%
% Revision 1.2  2003/03/04 21:46:19  roberto
% added CVS log entry and synchronized all copyright labels
%

if nargin<3
  method = 'banks';
end

switch lower(method)
case 'banks'
  npnt   = size(pnt,1);
  ndhk   = size(dhk,1);
  insert = spalloc(3*npnt,3*npnt,3*ndhk);

  dhkr  = zeros(4*ndhk,3);		% allocate memory for the new triangles
  pntr  = zeros(npnt+3*ndhk,3);		% allocate memory for the maximum number of new vertices
  pntr(1:npnt,:) = pnt;			% insert the original vertices
  current = npnt;

  for i=1:ndhk

    if ~insert(dhk(i,1),dhk(i,2))
      current = current + 1;
      pntr(current,:) = (pnt(dhk(i,1),:) + pnt(dhk(i,2),:))/2;
      insert(dhk(i,1),dhk(i,2)) = current;
      insert(dhk(i,2),dhk(i,1)) = current;
      v12 = current;
    else
      v12 = insert(dhk(i,1),dhk(i,2));
    end

    if ~insert(dhk(i,2),dhk(i,3))
      current = current + 1;
      pntr(current,:) = (pnt(dhk(i,2),:) + pnt(dhk(i,3),:))/2;
      insert(dhk(i,2),dhk(i,3)) = current;
      insert(dhk(i,3),dhk(i,2)) = current;
      v23 = current;
    else
      v23 = insert(dhk(i,2),dhk(i,3));
    end

    if ~insert(dhk(i,3),dhk(i,1))
      current = current + 1;
      pntr(current,:) = (pnt(dhk(i,3),:) + pnt(dhk(i,1),:))/2;
      insert(dhk(i,3),dhk(i,1)) = current;
      insert(dhk(i,1),dhk(i,3)) = current;
      v31 = current;
    else
      v31 = insert(dhk(i,3),dhk(i,1));
    end

    % add the 4 new triangles with the correct indices
    dhkr(4*(i-1)+1, :) = [dhk(i,1) v12 v31];
    dhkr(4*(i-1)+2, :) = [dhk(i,2) v23 v12];
    dhkr(4*(i-1)+3, :) = [dhk(i,3) v31 v23];
    dhkr(4*(i-1)+4, :) = [v12 v23 v31];

  end

  % remove the space for the vertices that was not used
  pntr = pntr(1:current, :);

case 'updown'
  ndhk = size(dhk,1);
  while ndhk<varargin{1}
    % increase the number of triangles by a factor of 4
    [pnt, dhk] = refine(pnt, dhk, 'banks');
    ndhk = size(dhk,1);
  end
  % reduce number of triangles using Matlab function
  [dhkr, pntr] = reducepatch(dhk, pnt, varargin{1});

otherwise
  error(['unsupported method: ' method]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hs, hc, contour] = triplot(pnt, tri, val, mode, levels)

% TRIPLOT make 2D or 3D plot of triangulated surface and interpolated values
% The surface can be displayed with linear interpolated values or with
% linear interpolated contours of a potential distribution.
%
% Use as
%   triplot(pnt, tri, value)
%   triplot(pnt, tri, value, mode)
%   triplot(pnt, tri, value, mode, levels)
% This will make a plot of value on the surface described by triangles
% tri with vertices pnt. The matrix tri can be [], in which case a
% it will be computed using a delaunay triangulation.
%
% The visualization mode can be
%   'surface'     make interpolated plot of value on surface (default)
%   'faces'       plot white triangles only        (value can be [])
%   'faces_red'   plot red triangles only          (value can be [])
%   'faces_blue'  plot blue triangles only         (value can be [])
%   'faces_skin'  plot skin-colored triangles only (value can be [])
%   'face_index'  plot index of each triangle      (value can be [])
%   'nodes'       plot black vertices only         (value can be [])
%   'node_index'  plot index of each vertex        (value can be [])
%   'node_label'  plot label of each vertex        (value should be cell array)
%   'edges'       plot black edges only            (value can be [])
%   'contour'     make interpolated contour plot of value on surface
%   'contour_bw'  make interpolated black-white contour plot
%   'contour_rb'  make interpolated contour plot with red-blue
%
% With the optional levels, you can specify the levels at which contours will
% be plotted
%
% See also PATCH, COLORMAP, VIEW (general Matlab commands)

% updated on Mon Jul 23 12:41:44 MET DST 2001
% updated on Fri Jan 31 11:47:28 CET     2003

% Copyright (C) 2001=2006, Robert Oostenveld
%
% $Log: triplot.m,v $
% Revision 1.7  2008/06/24 13:37:51  roboos
% added option faces_blue
% always return handle to objects that were plotted
%
% Revision 1.6  2007/01/03 17:00:35  roboos
% updated documentation, changed layout of code and comments
%
% Revision 1.5  2006/09/19 16:11:35  roboos
% added support for line segments, removed "axis equal" at end
%
% Revision 1.4  2006/05/02 19:15:28  roboos
% added 'faces_red' style
%
% Revision 1.3  2004/06/28 07:51:39  roberto
% improved documentation, added faces_skin
%
% Revision 1.2  2003/03/17 10:37:29  roberto
% improved general help comments and added copyrights
%

% start with empty return values
hs      = [];
hc      = [];
contour = [];

% everything is added to the current figure
holdflag = ishold;
hold on

% check the input variables
if ~isempty(val)
  val = val(:);
end

if nargin<4
  mode = 'surface';
end

if isempty(tri) & ~strcmp(mode, 'nodes')
  % no triangulation was specified but a triangulation is needed
  if size(pnt,2)==2
    % make a 2d triangulation of the points using delaunay
    tri = delaunay(pnt(:,1), pnt(:,2));
  else
    % make a 2d triangulation of the projected points using delaunay
    prj = elproj(pnt);
    tri = delaunay(prj(:,1), prj(:,2));
  end
end

if size(tri,2)==2
  % lines are specified instead of triangles, convert to triangles
  tri(:,3) = tri(:,2);
end

if nargin<3
  warning('only displaying triangle edges')
  mode='edges';
  val = [];
elseif nargin<4
  % warning('default displaying surface')
  mode='surface';
elseif nargin<5
  % determine contour levels
  if ~isempty(val) & ~iscell(val)
    absmax = max(abs([min(val) max(val)]));
    levels = linspace(-absmax,absmax,21);
  else
    levels = [];
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute contours for 2D or 3D triangulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(mode, 'contour') || strcmp(mode, 'contour_bw') || strcmp(mode, 'contour_rb')
  triangle_val = val(tri);
  triangle_min = min(triangle_val, [], 2);
  triangle_max = max(triangle_val, [], 2);

  for cnt_indx=1:length(levels)
    cnt = levels(cnt_indx);
    use = cnt>=triangle_min & cnt<=triangle_max;
    counter = 0;
    intersect1 = [];
    intersect2 = [];

    for tri_indx=find(use)'
      pos  = pnt(tri(tri_indx,:), :);
      v(1) = triangle_val(tri_indx,1);
      v(2) = triangle_val(tri_indx,2);
      v(3) = triangle_val(tri_indx,3);
      la(1) = (cnt-v(1)) / (v(2)-v(1));	% abcissa between vertex 1 and 2
      la(2) = (cnt-v(2)) / (v(3)-v(2));	% abcissa between vertex 2 and 3
      la(3) = (cnt-v(3)) / (v(1)-v(3));	% abcissa between vertex 1 and 2
      abc(1,:) = pos(1,:) + la(1) * (pos(2,:) - pos(1,:));
      abc(2,:) = pos(2,:) + la(2) * (pos(3,:) - pos(2,:));
      abc(3,:) = pos(3,:) + la(3) * (pos(1,:) - pos(3,:));
      counter = counter + 1;
      sel     = find(la>=0 & la<=1);
      intersect1(counter, :) = abc(sel(1),:);
      intersect2(counter, :) = abc(sel(2),:);
    end

    % remember the details for external reference
    contour(cnt_indx).level = cnt;
    contour(cnt_indx).n     = counter;
    contour(cnt_indx).intersect1 = intersect1;
    contour(cnt_indx).intersect2 = intersect2;
  end

  % collect all different contourlevels for plotting
  intersect1 = [];
  intersect2 = [];
  cntlevel   = [];
  for cnt_indx=1:length(levels)
    intersect1 = [intersect1; contour(cnt_indx).intersect1];
    intersect2 = [intersect2; contour(cnt_indx).intersect2];
    cntlevel   = [cntlevel; ones(contour(cnt_indx).n,1) * levels(cnt_indx)];
  end

  X = [intersect1(:,1) intersect2(:,1)]';
  Y = [intersect1(:,2) intersect2(:,2)]';
  C = [cntlevel(:)     cntlevel(:)]';

  if size(pnt,2)>2
    Z = [intersect1(:,3) intersect2(:,3)]';
  else
    Z = zeros(2, length(cntlevel));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the desired detail
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch lower(mode)

  case 'faces'
    % plot the faces of the 2D or 3D triangulation
    hs = patch('Vertices', pnt, 'Faces', tri);
    set(hs, 'FaceColor', 'white');
    set(hs, 'EdgeColor', 'none');

  case 'faces_skin'
    % plot the faces of the 2D or 3D triangulation
    skin_surface   = [255 213 119]/255;
    inner_skull_surface  = [202 100 100]/255;
    cortex = [255 213 119]/255;
    hs = patch('Vertices', pnt, 'Faces', tri);
    set(hs, 'FaceColor', skin_surface);
    set(hs, 'EdgeColor', 'none');
    lighting gouraud
    material shiny
    camlight

  case 'faces_red'
    % plot the faces of the 2D or 3D triangulation
    hs = patch('Vertices', pnt, 'Faces', tri);
    set(hs, 'FaceColor', [1 0 0]);
    set(hs, 'EdgeColor', 'none');
    lighting gouraud
    material shiny
    camlight

  case 'faces_blue'
    % plot the faces of the 2D or 3D triangulation
    hs = patch('Vertices', pnt, 'Faces', tri);
    set(hs, 'FaceColor', [0 0 1]);
    set(hs, 'EdgeColor', 'none');
    lighting gouraud
    material shiny
    camlight

  case 'face_index'
    % plot the triangle indices (numbers) at each face
    for face_indx=1:size(tri,1)
      str = sprintf('%d', face_indx);
      tri_x = (pnt(tri(face_indx,1), 1) +  pnt(tri(face_indx,2), 1) +  pnt(tri(face_indx,3), 1))/3;
      tri_y = (pnt(tri(face_indx,1), 2) +  pnt(tri(face_indx,2), 2) +  pnt(tri(face_indx,3), 2))/3;
      tri_z = (pnt(tri(face_indx,1), 3) +  pnt(tri(face_indx,2), 3) +  pnt(tri(face_indx,3), 3))/3;
      h   = text(tri_x, tri_y, tri_z, str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
      hs  = [hs; h];
    end

  case 'edges'
    % plot the edges of the 2D or 3D triangulation
    hs = patch('Vertices', pnt, 'Faces', tri);
    set(hs, 'FaceColor', 'none');
    set(hs, 'EdgeColor', 'black');

  case 'nodes'
    % plot the nodes (vertices) only as points
    if size(pnt, 2)==2
      hs = plot(pnt(:,1), pnt(:,2), 'k.');
    else
      hs = plot3(pnt(:,1), pnt(:,2), pnt(:,3), 'k.');
    end

  case 'nodes_blue'
    % plot the nodes (vertices) only as points
    if size(pnt, 2)==2
      hs = plot(pnt(:,1), pnt(:,2), 'b.', 'MarkerSize', 20);
    else
      hs = plot3(pnt(:,1), pnt(:,2), pnt(:,3), 'b.', 'MarkerSize', 20);
    end

  case 'nodes_red'
    % plot the nodes (vertices) only as points
    if size(pnt, 2)==2
      hs = plot(pnt(:,1), pnt(:,2), 'r.', 'MarkerSize', 20);
    else
      hs = plot3(pnt(:,1), pnt(:,2), pnt(:,3), 'r.', 'MarkerSize', 20);
    end

  case 'node_index'
    % plot the vertex indices (numbers) at each node
    for node_indx=1:size(pnt,1)
      str = sprintf('%d', node_indx);
      if size(pnt, 2)==2
        h   = text(pnt(node_indx, 1), pnt(node_indx, 2), str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
      else
        h   = text(pnt(node_indx, 1), pnt(node_indx, 2), pnt(node_indx, 3), str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
      end
      hs  = [hs; h];
    end

  case 'node_label'
    % plot the vertex indices (numbers) at each node
    for node_indx=1:size(pnt,1)
      str = val{node_indx};
      if ~isempty(str)
        if size(pnt, 2)==2
          h   = text(pnt(node_indx, 1), pnt(node_indx, 2), str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
        else
          h   = text(pnt(node_indx, 1), pnt(node_indx, 2), pnt(node_indx, 3), str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
        end
      else
        h = -1;
      end
      hs  = [hs; h];
    end

  case 'surface'
    % plot a 2D or 3D triangulated surface with linear interpolation
    if length(val)==size(pnt,1)
      hs = patch('Vertices', pnt, 'Faces', tri, 'FaceVertexCData', val, 'FaceColor', 'interp');
    else
      hs = patch('Vertices', pnt, 'Faces', tri, 'CData', val, 'FaceColor', 'flat');
    end
    set(hs, 'EdgeColor', 'none');

  case 'contour_bw'
    % make black-white contours
    hc = [];
    for i=1:length(cntlevel)
      if cntlevel(i)>0
        linestyle = '-';
        linewidth = 1;
      elseif cntlevel(i)<0
        linestyle = '--';
        linewidth = 1;
      else
        linestyle = '-';
        linewidth = 2;
      end
      h1 = patch('XData', X(:,i), 'Ydata', Y(:,i), ...
        'ZData', Z(:,i), 'CData', C(:,i), ...
        'facecolor','none','edgecolor','black', ...
        'linestyle', linestyle, 'linewidth', linewidth, ...
        'userdata',cntlevel(i));
      hc = [hc; h1];
    end

  case 'contour_rb'
    % make red-blue contours
    hc = [];
    for i=1:length(cntlevel)
      if cntlevel(i)>0
        edgecolor = 'red';
      elseif cntlevel(i)<0
        edgecolor = 'blue';
      else
        edgecolor = 'black';
      end
      h1 = patch('XData', X(:,i), 'Ydata', Y(:,i), ...
        'ZData', Z(:,i), 'CData', C(:,i), ...
        'facecolor','none','edgecolor',edgecolor, ...
        'linestyle', '-', 'linewidth', 3, ...
        'userdata',cntlevel(i));
      hc = [hc; h1];
    end

  case 'contour'
    % make full-color contours
    hc = [];
    for i=1:length(cntlevel)
      h1 = patch('XData', X(:,i), 'Ydata', Y(:,i), ...
        'ZData', Z(:,i), 'CData', C(:,i), ...
        'facecolor','none','edgecolor','flat',...
        'userdata',cntlevel(i));
      hc = [hc; h1];
    end

end	% switch

axis off
axis vis3d
axis equal

if nargout==0
  clear contour hc hs
end

if ~holdflag
  hold off
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [proj] = elproj(pos, method);

% ELPROJ makes a azimuthal projection of a 3D electrode cloud
%  on a plane tangent to the sphere fitted through the electrodes
%  the projection is along the z-axis
%
%  [proj] = elproj([x, y, z], 'method');
%
% Method should be one of these:
%	  'gnomic'
%	  'stereographic'
%	  'ortographic'
%	  'inverse'
%	  'polar'
%
% Imagine a plane being placed against (tangent to) a globe. If
% a light source inside the globe projects the graticule onto
% the plane the result would be a planar, or azimuthal, map
% projection. If the imaginary light is inside the globe a Gnomonic
% projection results, if the light is antipodal a Sterographic,
% and if at infinity, an Orthographic.
%
% The default projection is a polar projection (BESA like).
% An inverse projection is the opposite of the default polar projection.

% Copyright (C) 2000-2008, Robert Oostenveld
%
% $Log: elproj.m,v $
% Revision 1.4  2008/05/15 10:54:24  roboos
% updated documentation
%
% Revision 1.3  2007/03/20 10:29:35  roboos
% renamed method 'default' into 'polar'
%
% Revision 1.2  2003/03/17 10:37:28  roberto
% improved general help comments and added copyrights
%

x = pos(:,1);
y = pos(:,2);
if size(pos, 2)==3
  z = pos(:,3);
end

if nargin<2
  method='polar';
end

if nargin<3
  secant=1;
end

if strcmp(method, 'orthographic')
  % this method compresses the lowest electrodes very much
  % electrodes on the bottom half of the sphere are folded inwards
  xp = x;
  yp = y;
  num = length(find(z<0));
  str = sprintf('%d electrodes may be folded inwards in orthographic projection\n', num);
  if num
    warning(str);
  end
  proj = [xp yp];

elseif strcmp(method, 'gnomic')
  % the lightsource is in the middle of the sphere
  % electrodes on the equator are projected at infinity
  % electrodes below the equator are not projected at all
  rad = mean(sqrt(x.^2 + y.^2 + z.^2));
  phi = cart2pol(x, y);
  th  = atan(sqrt(x.^2 + y.^2) ./ z);
  xp  = cos(phi) .* tan(th) .* rad;
  yp  = sin(phi) .* tan(th) .* rad;
  num = length(find(th==pi/2 | z<0));
  str = sprintf('removing %d electrodes from gnomic projection\n', num);
  if num
    warning(str);
  end
  xp(find(th==pi/2 | z<0)) = NaN;
  yp(find(th==pi/2 | z<0)) = NaN;
  proj = [xp yp];

elseif strcmp(method, 'stereographic')
  % the lightsource is antipodal (on the south-pole)
  rad = mean(sqrt(x.^2 + y.^2 + z.^2));
  z   = z + rad;
  phi = cart2pol(x, y);
  th  = atan(sqrt(x.^2 + y.^2) ./ z);
  xp  = cos(phi) .* tan(th) .* rad * 2;
  yp  = sin(phi) .* tan(th) .* rad * 2;
  num = length(find(th==pi/2 | z<0));
  str = sprintf('removing %d electrodes from stereographic projection\n', num);
  if num
    warning(str);
  end
  xp(find(th==pi/2 | z<0)) = NaN;
  yp(find(th==pi/2 | z<0)) = NaN;
  proj = [xp yp];

elseif strcmp(method, 'inverse')
  % compute the inverse projection of the default angular projection
  [th, r] = cart2pol(x, y);
  [xi, yi, zi] = sph2cart(th, pi/2 - r, 1);
  proj = [xi, yi, zi];

else
  % use default angular projection
  [az, el, r] = cart2sph(x, y, z);
  [x, y] = pol2cart(az, pi/2 - el);
  proj = [x, y];
end

author: spm


text:
function [om2,T2,dom2dom,dom2dT,dT2dom,dT2dT] = inverse_motion(om,T);

% This function computes the inverse motion corresponding to (om,T)

om2 = -om;
dom2dom = -eye(3);
dom2dT = zeros(3,3);

[R,dRdom] = rodrigues(om);
Rinv = R';
dRinvdR = zeros(9,9);
dRinvdR([1 4 7],[1 2 3]) = eye(3);
dRinvdR([2 5 8],[4 5 6]) = eye(3);
dRinvdR([3 6 9],[7 8 9]) = eye(3);
dRinvdom = dRinvdR * dRdom;

Tt = Rinv * T;
[dTtdRinv,dTtdT] = dAB(Rinv,T);

T2 = -Tt;

dT2dom = - dTtdRinv * dRinvdom;
dT2dT = - dTtdT;

return;

% Test of the Jacobians:

om = randn(3,1);
T = 10*randn(3,1);
[om2,T2,dom2dom,dom2dT,dT2dom,dT2dT] = inverse_motion(om,T);

dom = randn(3,1) / 100000;
dT  = randn(3,1) / 100000;

[om3r,T3r] = inverse_motion(om+dom,T+dT);

om3p = om2 + dom2dom*dom +  dom2dT*dT;
T3p  =  T2 + dT2dom*dom  +  dT2dT*dT;

%norm(om3r - om2) / norm(om3r - om3p)  %-> Leads to infinity, since the opreation is linear!
norm(T3r - T2) / norm(T3r - T3p)

repo: ekfmonoslam