text:
%TestSuiteTest Unit tests for TestSuite class

classdef TestSuiteTest < TestCaseInDir

   methods
      function self = TestSuiteTest(name)
         self = self@TestCaseInDir(name, ...
             fullfile(fileparts(which(mfilename)), 'helper_classes'));
      end

      function testClassNameIn(self)
         % Syntax check: TestSuite('classname')
         suite = TestSuite('TwoPassingTests');
         assertTrue(numel(suite.TestComponents) == 2, ...
            'TestSuite finds two test methods given class name');
      end

      function testCurrentDirectory(self)
         % See that the no-input syntax executes without error.
         % Not sure how to test this more effectively.
         suite = TestSuite();
      end

      function testNoTestMethods(self)
         % TestCase class containing no test methods
         suite = TestSuite('NoTestMethods');
         assertTrue(numel(suite.TestComponents) == 0, ...
            'No test cases when class contains no test methods');
      end

      function test_fromTestCaseClassName(self)
          suite = TestSuite.fromTestCaseClassName('TwoPassingTests');
          assertTrue(numel(suite.TestComponents) == 2);
          assertTrue(ismember(suite.TestComponents{1}.Name, ...
              {'testMethod1', 'testMethod2'}));
          assertTrue(ismember(suite.TestComponents{2}.Name, ...
              {'testMethod1', 'testMethod2'}));
      end

      function test_fromTestCaseClassName_badclass(self)
          assertExceptionThrown(@() TestSuite.fromTestCaseClassName('atan2'), ...
              'xunit:fromTestCaseClassName');
      end

      function test_fromName_TestCaseSubclass(self)
          suite = TestSuite.fromName('TwoPassingTests');
          assertTrue(numel(suite.TestComponents) == 2);
          assertEqual(suite.Name, 'TwoPassingTests');
      end

      function test_fromName_notTestCaseSubclass(self)
          suite = TestSuite.fromName('TestRunMonitor');
          assertTrue(isempty(suite.TestComponents));
          assertEqual(suite.Name, 'TestRunMonitor');
      end

      function test_fromName_simpleTest(self)
          suite = TestSuite.fromName('testSimple');
          assertEqual(numel(suite.TestComponents), 1);
          assertEqual(suite.Name, 'testSimple');
          assertEqual(suite.Location, which('testSimple'));
      end

      function test_fromName_subfunctions(self)
          suite = TestSuite.fromName('testFunctionHandlesA');
          assertEqual(numel(suite.TestComponents), 2);
          assertEqual(suite.Name, 'testFunctionHandlesA');
          assertEqual(suite.Location, which('testFunctionHandlesA'));
      end

      function test_fromName_bogus_name(self)
          suite = TestSuite.fromName('atan2');
          assertTrue(isempty(suite.TestComponents));
          assertEqual(suite.Name, 'atan2');
      end

      function test_fromName_with_filter_string(self)
          suite = TestSuite.fromName('testFunctionHandlesA:testA');
          assertEqual(numel(suite.TestComponents), 1);
          assertEqual(suite.TestComponents{1}.Name, 'testA');
          assertEqual(suite.Name, 'testFunctionHandlesA');
      end

      function test_fromName_with_nonmatching_filter_string(self)
          suite = TestSuite.fromName('testFunctionHandlesA:foobar');
          assertTrue(isempty(suite.TestComponents));
      end

      function test_fromName_with_dirname(self)
         xunit_test_dir = which('TestSuiteTest');
         xunit_test_dir = fileparts(xunit_test_dir);
         cwd_test_dir = fullfile(xunit_test_dir, 'cwd_test');
         suite = TestSuite.fromName(cwd_test_dir);

         assertEqual(suite.Name, 'cwd_test');
         assertEqual(suite.Location, cwd_test_dir);
         assertEqual(numel(suite.TestComponents), 3);
      end

      function test_fromPwd(self)
          % Verify that the fromPwd method returns a nonempty TestSuite object
          % from the helper_classes directory, with the correct number of
          % test components.
          suite = TestSuite.fromPwd();
          assertTrue(isa(suite, 'TestSuite'));
          assertTrue(numel(suite.TestComponents) == 16);
      end

   end

end

author: ricedsp


text:
function [ dat ] = prep_filter( dat, varargin )
% prep_filter (Pre-processing procedure):
%
% Description:
%     This function filters the data within specified frequency band
%
% Example:
%    EEG.data=prep_filter(EEG.data, {'frequency', [7 13];'fs',100});
%
% Input:
%     dat       - EEG data structure
% Option:
%     frequency - Frequency range that you want to filter
%     fs        - Sampling frequency
% Return:
%     dat       - Spectrally filtered data
%
% Min-ho Lee
% mhlee@image.korea.ac.kr
%

if iscell(varargin{:})
    opt=opt_cellToStruct(varargin{:});
elseif isstruct(varargin{:}) % already structure(x-validation)
    opt=varargin{:}
end

switch isstruct(dat)
    case true %struct
        if ~isfield(dat, 'x') && ~isfield(dat, 'fs')
            warning('Parameter is missing: dat.x or dat.fs');
        end

        if ~isfield(opt,'fs')
            if isfield(dat, 'fs')
                opt.fs=dat.fs;
            else
                error('Parameter is missing: fs');
            end
        end
        tDat=dat.x;
        if ndims(tDat)==3   %smt
            [nD nT nC]=size(tDat);
            tDat=reshape(tDat, [nD*nT,nC]);
            band=opt.frequency;
            [b,a]= butter(5, band/(opt.fs/2),'bandpass');
            tDat(:,:)=filter(b, a, tDat(:,:));
            tDat=reshape(tDat, [nD, nT,nC]);
            dat.x=tDat;
        elseif ndims(tDat)==2  %cnt
            band=opt.frequency;
            [b,a]= butter(5, band/(opt.fs/2),'bandpass');
            tDat(:,:)=filter(b, a, tDat(:,:));
            dat.x=tDat;
        end

    case false
        % add if dat is not struct
        tDat=dat;
        if ndims(tDat)==3   %smt
            [nD nT nC]=size(tDat);
            tDat=reshape(tDat, [nD*nT,nC]);
            band=opt.frequency;
            [b,a]= butter(5, band/(opt.fs/2),'bandpass');
            tDat(:,:)=filter(b, a, tDat(:,:));
            tDat=reshape(tDat, [nD, nT,nC]);
            dat.x=tDat;
        elseif ndims(tDat)==2  %cnt
            band=opt.frequency;
            [b,a]= butter(5, band/(opt.fs/2),'bandpass');
            %% TODO: It' doesn't work when the sampling rate is 1000 Hz;
            tDat(:,:)=filter(b, a, tDat(:,:));
            fld='x';
            dat=tDat;
        end

        % History
        if isfield(dat,'stack')
            c = mfilename('fullpath');
            c = strsplit(c,'\');
            dat.stack{end+1}=c{end};
        end
end

author: PatternRecognition

repo: OpenBMI


text:
function [x,param,values] = coarserefine(f,intrv,theta,N,alpha)
% 1-D adaptive residual subsampling method for radial basis function
% interpolation
%
% f: function defined on interval intrv = [a b].
% intrv: interval [a b] where f is defined.
% theta: threshold interval [thetar thetac] where thetac < thetar.
%        thetar: refinement threshold.
%        thetac: coarsening threshold.
% N: Initial equally-spaced N centers.
% alpha: global multiplier of multiquadric parameters
%
% Example 1 :  f = @(x) abs(x+.04)
%              coarserefine(f,[-1 1],[5e-5 5e-7],11,0.75)
%
% Example 2 :  coarserefine(@(x)myfun(x,5),[-1 1],[5e-5 5e-7],11,0.75)
%              %----------------------%
%              function y = myfun2(x,c)
%              y = 1./(1 + (c*x).^2);
%              %----------------------%
% For reference, see:
% Adaptive residual subsampling methods for radial basis function
% interpolation and collocation problems. submitted to Computers Math. Appl
%
% Tobin A. Driscoll and Alfa R.H. Heryudono    05/14/2006
% MATLAB 7 is recommended.

% Initial points
x = linspace(intrv(1),intrv(2),N)';
N = length(x); dx = diff(x);
epsilon = alpha*min([Inf;1./dx],[1./dx;Inf]);
y = x(1:N-1) + 0.5*dx;

refine = true;
while any(refine)
  A = zeros(N); B = zeros(N-1,N);
  for j=1:N
    A(:,j) = mq(x,x(j),epsilon(j));
    B(:,j) = mq(y,x(j),epsilon(j));
  end
  lambda = A\feval(f,x);
  resid = abs(B*lambda - feval(f,y));

  refine = resid > theta(1);
  fprintf('Adding %i centers.',sum(refine))
  vals = sortrows([[x;y(refine)] [feval(f,x);B(refine,:)*lambda]]);
  x = vals(:,1); vals(:,1)=[];
  dx = diff(x); epsilon = alpha*min([Inf;1./dx],[1./dx;Inf]);

  coarsen = resid(1:N-2) < theta(2) & resid(2:N-1) < theta(2);
  coarsen = 1+find(coarsen); x(coarsen) = []; vals(coarsen) = []; epsilon(coarsen)=[];
  fprintf(' Removing %i centers.\n',length(coarsen))

  N = length(x);
  y = x(1:N-1) + 0.5*diff(x);
end

if nargout > 1
    param = epsilon;
    if nargout > 2
        values = vals;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function DECODED = RS_E_E_DEC(received, erasures,n,k,t,h,g,field);

%Check for decoding failures
%Previous decoder RS_E_E_DEC4

%****************
%*** Decoding ***
%****************

%syndrome calculation
S = [];
%Subtitute alpha^i in received polynomial - Lin + Costello p.152 eq. 6.13
for ii = 1:2*t
    S(ii)= -Inf;
    for cc = 1:n
        S(ii) = gfadd(S(ii),gfmul(received(cc),gfpow(ii,cc-1,n),field),field); %Sum all the terms
    end
end
%S

%Test if syndrome  = 0, if syndrome equals 0, assume that no errors occured
for i = 1:2*t
    test_pol(i) = -Inf;
end

if all (S == test_pol)

    message = received;

    for i = 1:n
        if message(i) < 0
            message(i) = -Inf;
        end
    end

else

    %Compute the erasure locator polynomial:
    erasures_pos = erasures - 1;
    num_erasures = length(erasures);

    %Compute the erasure-locator polynomial
    erasure_loc_pol = 0;
    for i = 1:length(erasures_pos)
        erasure_loc_pol = gfconv(erasure_loc_pol, [0 erasures_pos(i)],field);
    end

    %Compute modified syndrome polynomial:
    S_pol = [-inf S];
    dividend = gfconv(erasure_loc_pol,gfadd(0,S_pol,field),field);
    dividend = gfadd(dividend,0,field);

    divisor = [];
    for i = 1:2*t+2
        divisor(i) = -Inf;
    end
    divisor(2*t+2) = 0;

    [q,mod_syn] = gfdeconv(dividend,divisor,field);

    while length(mod_syn) < h+1
        mod_syn = [mod_syn -Inf];
    end

    S_M = [];
    for i = 1:h - num_erasures
        S_M(i) = mod_syn(i + num_erasures + 1);
    end

    flag = 0;
    if isempty(S_M) == 1
        flag = 0;
    else
        for i = 1:length(S_M)
            if (S_M(i) ~= -Inf)
                flag = 1;     %Other errors occured in conjunction with erasures
            end
        end
    end

    %Find error-location polynomial sigma (Berlekamp's iterative algorithm -
    %sigma = [0 7 4 6]
    if (flag == 1)
        %sigma = M_B2(n,k,length(S_M) - 1,S_M,field);

        num_iter = t - num_erasures/2;

        sigma = massey_berlekamp_M3(n,k,num_iter,S_M,field);

        %Chien search
        %step 3 from Lin + Costello p.175
        %the error locating polynomial have a maximum of t entries
        error_loc = [];
        kk = 0;

        for ii = 0:n-1
            error_r = -Inf;
            for cc = 1:length(sigma)
                error_r = gfadd(error_r,gfmul(sigma(cc),gfpow(ii,cc-1,n),field),field); %Sum all the terms
            end
            if error_r == -Inf
                kk = kk + 1;
                error_loc(kk) = ii;
            end
        end

        % Test if the roots are distinct
        % Form a test polynomial by multiplying the roots of error_loc with each other
        % Divide the error_loc pol by test pol
        % if the degree of the quotient exceeds a constant, then the roots are
        % not distinct

        test_pol = 0;
        for ii = 1:length(error_loc)
            test_pol = gfconv(test_pol,[error_loc(ii) 0],field);
        end

        %test_pol
        %error_loc
        %sigma

        [QQ,RR] = gfdeconv(sigma,test_pol,field);
        if length(QQ) > 1
            DECODED = received;
            return
        end

        comp_error_locs = [];
        %Take reciprocals of elements in error_loc - error location numbers
        for ii = 1:length(error_loc)
            comp_error_locs(ii) = gfdiv(0,error_loc(ii),field);
        end
        %error_loc_p %places where errors occur
    else
        sigma = 0;
        comp_error_locs = [];
    end

    %Calculate error magnitudes - Forney algorithm?
    %Step 4. Lin and Costello - This program uses another algorithm from:
    %            drake.ee.washington.edu/~adina/rsc/slide/node9.html
    %            http://www.ee.ucla.edu/~matache/rsc/slide.html
    %Compute the error magnitude polynomial:
    %1.  Form the function [1 + S(x)]

    SS(1) = 0;
    for ii = 1: 2*t
        SS(ii+1) = S(ii);
    end

    %SS

    %2. form the product of SS and the KEY Equation
    %OMEGA = gfconv(SS,sigma,field);
    OMEGA = gfconv(sigma,gfadd(0,mod_syn,field),field);

    %3. OMEGA = (SS * sigma)mod(x^(2t+1))
    %3.1. Form a function := x^(2t+1)
    for ii = 1: (2*t)
        DIV(ii)= -Inf;
    end
    DIV(2*t+1) = 0;

    %3.2.  OMEGA = (SS * sigma)mod(x^(2t+1))
    [DUMMY, OMEGA] = gfdeconv(OMEGA,DIV,field);
    %OMEGA

    %4. Differentiate the key equation with respect to x
    %sigma_diff = gfdiff(sigma);
    tsi = gfconv(sigma,erasure_loc_pol,field);
    tsi_diff = gfdiff(tsi);

    e_e_places = [erasures_pos comp_error_locs];

    %Calculate the error magnitudes
    %Substitute the inverse into sigma_diff
    for ii = 1:length(e_e_places)
        %error_loc_p(ii)
        ERR_DEN = gfsubstitute(tsi_diff,gfdiv(0,e_e_places(ii),field),length(tsi_diff),n,field);
        ERR_NUM = gfsubstitute(OMEGA,gfdiv(0,e_e_places(ii),field),length(OMEGA),n,field);
        ERR_NUM = gfmul(ERR_NUM,e_e_places(ii),field);

        if ERR_DEN == -Inf
            DECODED = received;
            %display('Decoding Failure XXX')
            return
        end

        ERR(ii) = gfmul(ERR_NUM,gfdiv(0,ERR_DEN,field),field);
    end

    %error_loc_p
    %ERR

    %Determine introduced error
    for ii = 1:n
        ERR_p(ii) = -Inf;
    end

    %Error -  t must be substituted by amount of errors
    for ii = 1:length(e_e_places)
        pp = e_e_places(ii);
        ERR_p(pp+1) = ERR(ii);
    end

    %ERR_p

    message = gfadd(received,ERR_p,field);

end

DECODED = message;


text:
close all;
clear all;
clc;

data_file_path = 'bin/mmv_phase_transition_snr_30db_s_4.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'mmv_snr_30_db_s_4';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'MMV, SNR=30dB s=4';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function Study = readLeksellStudyFile(filename)
%"readLeksellStudyFile"
%   Reads a Leksell Study file using decodeLeksellData and places the
%   fields into a datastructure with meaningful names.  Many of the values
%   in Leksell files are a mystery at the moment, and are stored in
%   variables called "mystery_value_x" where x is a number.  If anyone
%   discovers the meaning of these values please report them on the CERR
%   webpage at radium.wustl.edu/cerr or to jalaly@radium.wustl.edu.
%
%JRA 6/13/05
%
%Usage:
%   Study = readLeksellStudyFile(filename)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

fid = fopen(filename, 'r', 'b');

data = decodeLeksellData(fid);

fclose(fid);

%If data does not have at least two elements it is effectively empty.
if length(data) < 2
    Study = [];
    return;
end

for i=1:length(data) - 1
    rawStudy = data{i};

	Study(i).modality           = rawStudy{1};
    Study(i).mystery_value_1    = rawStudy{2}; %always a small integer
    %Read in the original transM from the file
    originalTransM = reshape(rawStudy{3}, [4 4]);
    %Take the transpose and the inverse to put into the format used by CERR
    modifiedTransM = inv(originalTransM');
    %Divide the translations by 10 since the matrix was originally for mm
    %(converts the translations to cm)
    modifiedTransM(1,4) = modifiedTransM(1,4)/10;
    modifiedTransM(2,4) = modifiedTransM(2,4)/10;
    modifiedTransM(3,4) = modifiedTransM(3,4)/10;
    Study(i).rcsToxyzTransM = modifiedTransM;
    Study(i).pixelIntensityRange = rawStudy{4}; % appears to be the image intensity range, but does not need to be used for CERR import
    Study(i).registration_value = rawStudy{5};

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout = gammaincinv_a(varargin)

switch class(varargin{1})

    case 'double'
        varargout{1} = gammaincinv(varargin{2},varargin{1});

    case 'char'

        operator = CreateBasicOperator('increasing','positive','callback');
        operator.range = [0 inf];
        operator.domain = [0 inf];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout = samplinthist(varargin)
% VL_SAMPLINTHIST  Sample integral histogram
%   HISTS = VL_SAMPLINTHIST(INTHIST, BOXES) samples the integral
%   histogram INTHIST to obtain the histograms of the specified
%   BOXES.
%
%   INTHIST is a MxNxK array, where M x N are ``spatial'' dimensions,
%   and K is the number of histogram bins. INTHIST may be of class
%   UINT32 or DOUBLE.
%
%   Each box is a four dimensional vector [IMIN JMIN IMAX JMAX]' of
%   class UINT32 and correspond to the index set [IMIN, IMAX] x [JMIN,
%   JMAX]. To specify an empty box, let IMIN > IMAX.
%
%   HISTS stores one histogram per column (one for each box) and has K
%   rows, one for each histogram bin. HIST is of the same class of
%   INTHIST.
%
%   See also: VL_INTHIST(), VL_IMINTEGRAL(), VL_HELP().
[varargout{1:nargout}] = vl_samplinthist(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function bb = bb_correct(img,bb)
% Change

K  = 2;
K2 = 1;
while 1

    clf; imshow(img); bb_rectangle(bb,'EdgeColor','y');
    %maximize();
    [x,y,ch] = ginput(1);
    if isempty(ch), break; end
    switch ch
        case 29
            bb = [bb(1)+K; bb(2); bb(3)+K; bb(4)];
        case 28
            bb = [bb(1)-K; bb(2); bb(3)-K; bb(4)];
        case 30
            bb = [bb(1); bb(2)-K; bb(3); bb(4)-K];
        case 31
            bb = [bb(1); bb(2)+K; bb(3); bb(4)+K];
        case 46 % bigger
            bb = [bb(1)-K2; bb(2)-K2; bb(3)+K2; bb(4)+K2];
        case 44 % smaller
            bb = [bb(1)+K2; bb(2)+K2; bb(3)-K2; bb(4)-K2];

    end

end

repo: OpenTLD


text:
% AVEREF - convert common-reference EEG data to average reference
%            Note that this old function is not being used in EEGLAB. The
%            function used by EEGLAB is REREF.
%
% Usage:
%   >> data = averef(data);
%   >> [data_out W_out S_out meandata] = averef(data,W);
%
% Inputs:
%   data - 2D data matrix (chans,frames*epochs)
%   W    - ICA weight matrix
%
% Outputs:
%   data_out - Input data converted to average reference.
%   W_out    - ICA weight matrix converted to average reference
%   S_out    - ICA sphere matrix converted to EYE
%   meandata - (1,dataframes) mean removed from each data frame (point)
%
% Note: If 2 args, also converts the weight matrix W to average reference:
%         If ica_act = W*data, then data = inv(W)*ica_act;
%         If R*data is the average-referenced data,
%         R*data=(R*inv(W))*ica_act and W_out = inv(R*inv(W));
%         The average-reference ICA maps are the columns of inv(W_out).
%
% Authors: Scott Makeig and Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 1999
%
% See also: REREF

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

% 12/16/99 Corrected denomiator on the suggestion of Ian Nimmo-Smith, Cambridge UK
% 01-25-02 reformated help & license -ad

function [data, W, S, meandata] = averef(data, W, S)

if nargin<1
  help averef
  return
end
chans = size(data,1);
if chans < 2
  help averef
  return
end

% avematrix = eye(chans)-ones(chans)*1/chans;
% data = avematrix*data; % implement as a matrix multiply
% else (faster?)

meandata = sum(data)/chans;
data = data - ones(chans,1)*meandata;

% treat optional ica parameters
if nargin == 2
	winv  = pinv(W);
    size1 = size(winv,1);
	avematrix = eye(size1)-ones(size1)*1/size1;
	W = pinv(avematrix*winv);
end
if nargin >= 3
	winv = pinv(W*S);
    size1 = size(winv,1);
	avematrix = eye(size1)-ones(size1)*1/size1;
	W = pinv(avematrix*winv);
	S = eye(chans);
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%ISEMPTY Dataset overload
%
%	I = ISEMPTY(A,FIELD)
%
% INPUT
%  A     Dataset
%  FIELD Dataset field, default 'data'
%
% OUTPUT
%  I     Flag, 1 if field is empty, 0 otherwise.
%
% DESCRIPTION
% Dataset overload for ISEMPTY. This is particulary useful for
% ISEMPTY(A) to test on an empty dataset, and
% ISEMPTY(A,'prior') to test on an undefined PRIOR field
%
% SEE ALSO (PRTools Guide)
% PRDATASET, SETPRIOR, GETPRIOR

author: marianux

repo: ecg-kit

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


text:
function [intl] = embedbp(I,txt,b)
% EMBEDBP Embeds the data in the bth bit plane of the image
%
% [INTL] = EMBEDBP(I,TXT,B) embeds the string TXT in the Bth bit-plane of
% the image I and returns the watermarked image INTL. If B is not speified,
% it is taken as 1.
%
% See also RECOVERBP

if nargin == 2
    b=1;
end
N = 8*numel(txt);
S = numel(I);
if N > S
    warning('Text truncated to be within size of image')
    txt = txt(1:floor(S/8));
    N = 8*numel(txt);
end
p = 2^b;
h = 2^(b-1);
I1 = reshape(I,1,S);
addl = S-N;
dim = size(I);
I2 = round(abs(I1(1:N)));
si = sign(I1(1:N));
for k = 1:N
    if si(k) == 0
        si(k) = 1;
    end
    I2(k) = round(I2(k));
    if mod((I2(k)),p) >= h
        I2(k) = I2(k) - h;
    end
end
bt = dec2bin(txt,8);
bint = reshape(bt,1,N);
d = h*48;
bi = (h*bint) - d;
I3 = double(I2) + bi;
binadd = [bi zeros(1,addl)];
I4 = double(si).*double(I3);
I5 = [I4 I1(N+1:S)];
intl = reshape(I5,dim);
return


text:
function varargout = setupSolarized(varargin) %#ok<STOUT>
%===============================================================================
% FUNCTION:
%   setupSolarized.m
%
% PURPOSE:
%   Loads the Solarized color scheme, either light or dark as specified by the
%   user. More information about the Solarzied project can be found online:
%
%       http://ethanschoonover.com/solarized
%
% RESULTS OF SCRIPT:
%   o The script loads the color pallet as defined by the solarized project
%   o Optional string inputs are 'light', ''dark', or 'default' (see examples
%   below for proper usage. The default option will restore the Matlab default
%   color scheme.
%
% EXAMPLE:
%   o setupSolarized('light');
%   o setupSolarized('dark');
%   o setupSolarized('default');
%===============================================================================

%===============================================================================
% REVISION HISTORY:
%   o Add colors setup to return to default
%       - Wang Yudong / 2015-01-05
%   o Utilities Toolbox:
%       - B. Hager / 2013-05-31 <- Original and debugged.
%
% COPYRIGHT:
%   This work and derivatives of this work are copyrighted by Honda
%   Aircraft Company. Any use without permission is strictly
%   prohibited.
%===============================================================================
%% Initial Setup
if (nargin==1)
    slrzd = varargin{1};
elseif (nargin==0)
    slrzd = 'default';
end

% Setup the color pallet           % SOLARIZED HEX     16/8 TERMCOL  XTERM/HEX   L*A*B      RGB         HSB
% divide by 256 for Matlab         % --------- ------- ---- -------  ----------- ---------- ----------- -----------
sol.base03  = [  0  43  54] / 256; % base03    #002b36  8/4 brblack  234 #1c1c1c 15 -12 -12   0  43  54 193 100  21
sol.base02  = [  7  54  66] / 256; % base02    #073642  0/4 black    235 #262626 20 -12 -12   7  54  66 192  90  26
sol.base01  = [ 88 110 117] / 256; % base01    #586e75 10/7 brgreen  240 #585858 45 -07 -07  88 110 117 194  25  46
sol.base00  = [101 123 131] / 256; % base00    #657b83 11/7 bryellow 241 #626262 50 -07 -07 101 123 131 195  23  51
sol.base0   = [131 148 150] / 256; % base0     #839496 12/6 brblue   244 #808080 60 -06 -03 131 148 150 186  13  59
sol.base1   = [147 161 161] / 256; % base1     #93a1a1 14/4 brcyan   245 #8a8a8a 65 -05 -02 147 161 161 180   9  63
sol.base2   = [238 232 213] / 256; % base2     #eee8d5  7/7 white    254 #e4e4e4 92 -00  10 238 232 213  44  11  93
sol.base3   = [253 246 227] / 256; % base3     #fdf6e3 15/7 brwhite  230 #ffffd7 97  00  10 253 246 227  44  10  99
sol.yellow  = [181 137   0] / 256; % yellow    #b58900  3/3 yellow   136 #af8700 60  10  65 181 137   0  45 100  71
sol.orange  = [203  75  22] / 256; % orange    #cb4b16  9/3 brred    166 #d75f00 50  50  55 203  75  22  18  89  80
sol.red     = [220  50  47] / 256; % red       #dc322f  1/1 red      160 #d70000 50  65  45 220  50  47   1  79  86
sol.magenta = [211  54 130] / 256; % magenta   #d33682  5/5 magenta  125 #af005f 50  65 -05 211  54 130 331  74  83
sol.violet  = [108 113 196] / 256; % violet    #6c71c4 13/5 brmagenta 61 #5f5faf 50  15 -45 108 113 196 237  45  77
sol.blue    = [ 38 139 210] / 256; % blue      #268bd2  4/4 blue      33 #0087ff 55 -10 -45  38 139 210 205  82  82
sol.cyan    = [ 42 161 152] / 256; % cyan      #2aa198  6/6 cyan      37 #00afaf 60 -35 -05  42 161 152 175  74  63
sol.green   = [133 153   0] / 256; % green     #859900  2/2 green     64 #5f8700 60 -20  65 133 153   0  68 100  60
mat.black   = [  0   0   0] / 256;
mat.white   = [255 255 255] / 256;
mat.blue    = [  0   0 255] / 256;
mat.green   = [ 34 139  34] / 256;
mat.violet  = [160  32 240] / 256;
mat.red178  = [178   0   0] / 256;
mat.red230  = [230   0   0] / 256;
mat.golden  = [178 140   0] / 256;
mat.orange  = [255 148   0] / 256;
mat.chamois = [238 225 180] / 256;
mat.iceblue = [202 232 232] / 256;
mat.persian = [  0 163 163] / 256;
mat.yellow  = [252 252 220] / 256;

%% Setup Color Scheme for Light/Dark Options
if strcmp(slrzd,'light')
    sycb = 0;           % Set 'Use system colors' checkbox to False
    txc  = sol.base00;  % Set 'Text' color
    bgc  = sol.base3;   % Set 'Background' color
    kwd  = sol.blue;    % Set 'Keywords' color
    cmt  = sol.cyan;    % Set 'Comments' color - spec calls for base1, not sure about that...
    str  = sol.magenta; % Set 'Strings' color
    ustr = sol.orange;  % Set 'Unterminated strings' color
    scmd = sol.yellow;  % Set 'System commands' color
    errs = sol.red;     % Set 'Errors' color
    hyp  = sol.blue;    % Set 'Hyperlinks' color
    warn = sol.orange;  % Set 'Warnings' color
    afhb = 0;           % Set 'Autofix highlight' checkbox to False
    afh  = sol.base2;   % Set 'Autofix highlight' color
%     ahib = 0;           % Set 'Automatically highlight' checkbox to False
%     ahi  = sol.violet;  % Set 'Automatically highlight' color
    vwss = sol.green;   % Set 'Variables with shared scope' color
    hsb  = 0;           % Set 'Highlight sections' checkbox to False
    hsbc = mat.yellow;  % Set 'Highlight sections' color to default yellow
    hclb = 1;           % Set 'Highlight current line' checkbox to True
    hcl  = sol.base2;   % Set 'Highlight current line' color
    slnb = 1;           % Set 'Show line numbers' checkbox to True
    rtlb = 1;           % Set 'Show line' checkbox in Right-hand text limit to True
elseif strcmp(slrzd,'dark')
    sycb = 0;           % Set 'Use system colors' checkbox to False
    txc  = sol.base0;   % Set 'Text' color
    bgc  = sol.base03;  % Set 'Background' color
    kwd  = sol.blue;    % Set 'Keywords' color
    cmt  = sol.cyan;    % Set 'Comments' color - spec calls for base1, not sure about that...
    str  = sol.magenta; % Set 'Strings' color
    ustr = sol.orange;  % Set 'Unterminated strings' color
    scmd = sol.yellow;  % Set 'System commands' color
    errs = sol.red;     % Set 'Errors' color
    hyp  = sol.blue;    % Set 'Hyperlinks' color
    warn = sol.orange;  % Set 'Warnings' color
    afhb = 0;           % Set 'Autofix highlight' checkbox to False
    afh  = sol.base00;  % Set 'Autofix highlight' color
%     ahib = 0;           % Set 'Automatically highlight' checkbox to False
%     ahi  = sol.violet;  % Set 'Automatically highlight' color
    vwss = sol.green;   % Set 'Variables with shared scope' color
    hsb  = 0;           % Set 'Highlight sections' checkbox to False
    hsbc = sol.base02;  % Set 'Highlight sections' color to default yellow
    hclb = 1;           % Set 'Highlight current line' checkbox to True
    hcl  = sol.base02;  % Set 'Highlight current line' color
    slnb = 1;           % Set 'Show line numbers' checkbox to True
    rtlb = 1;           % Set 'Show line' checkbox in Right-hand text limit to True
elseif strcmp(slrzd,'default')
    sycb = 1;           % Set 'Use system colors' checkbox to True
    txc  = mat.black;   % Set 'Text' color
    bgc  = mat.white;   % Set 'Background' color
    kwd  = mat.blue;    % Set 'Keywords' color
    cmt  = mat.green;   % Set 'Comments' color - spec calls for base1, not sure about that...
    str  = mat.violet;  % Set 'Strings' color
    ustr = mat.red178;  % Set 'Unterminated strings' color
    scmd = mat.golden;  % Set 'System commands' color
    errs = mat.red230;  % Set 'Errors' color
    hyp  = mat.blue;    % Set 'Hyperlinks' color
    warn = mat.orange;  % Set 'Warnings' color
    afhb = 1;           % Set 'Autofix highlight' checkbox to True
    afh  = mat.chamois; % Set 'Autofix highlight' color
%     ahib = true;        % Set 'Automatically highlight' checkbox to True
%     ahi  = mat.iceblue; % Set 'Automatically highlight' color
    vwss = mat.persian; % Set 'Variables with shared scope' color
    hsb  = 1;           % Set 'Highlight sections' checkbox to True
    hsbc = mat.yellow;  % Set 'Highlight sections' color to default yellow
    hclb = 0;           % Set 'Highlight current line' checkbox to False
    hcl  = sol.base2;   % Set 'Highlight current line' color
    slnb = 1;           % Set 'Show line numbers' checkbox to True
    rtlb = 1;           % Set 'Show line' checkbox in Right-hand text limit to True
end

%% Desktop tool colors
com.mathworks.services.Prefs.setBooleanPref('ColorsUseSystem',sycb); clear('sycb')

com.mathworks.services.Prefs.setColorPref('ColorsText',java.awt.Color(txc(1), txc(2), txc(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsText'); clear('txc')

com.mathworks.services.Prefs.setColorPref('ColorsBackground',java.awt.Color(bgc(1), bgc(2), bgc(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsBackground'); clear('bgc')

%% Setup MATLAB syntax highlighting colors
com.mathworks.services.Prefs.setColorPref('Colors_M_Keywords',java.awt.Color(kwd(1), kwd(2), kwd(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Keywords'); clear('kwd')

com.mathworks.services.Prefs.setColorPref('Colors_M_Comments',java.awt.Color(cmt(1), cmt(2), cmt(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Comments'); clear('cmt')

com.mathworks.services.Prefs.setColorPref('Colors_M_Strings',java.awt.Color(str(1), str(2), str(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Strings'); clear('str')

com.mathworks.services.Prefs.setColorPref('Colors_M_UnterminatedStrings',java.awt.Color(ustr(1), ustr(2), ustr(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_UnterminatedStrings'); clear('ustr')

com.mathworks.services.Prefs.setColorPref('Colors_M_SystemCommands',java.awt.Color(scmd(1), scmd(2), scmd(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_SystemCommands'); clear('scmd')

com.mathworks.services.Prefs.setColorPref('Colors_M_Errors',java.awt.Color(errs(1), errs(2), errs(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Errors');clear('errs')

%% Other colors
com.mathworks.services.Prefs.setColorPref('Colors_HTML_HTMLLinks',java.awt.Color(hyp(1), hyp(2), hyp(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_HTML_HTMLLinks'); clear('hyp')

%% Code analyzer colors
com.mathworks.services.Prefs.setColorPref('Colors_M_Warnings',java.awt.Color(warn(1), warn(2), warn(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Warnings'); clear('warn')

com.mathworks.services.Prefs.setBooleanPref('ColorsUseMLintAutoFixBackground',afhb);
com.mathworks.services.Prefs.setColorPref('ColorsMLintAutoFixBackground',java.awt.Color(afh(1), afh(2), afh(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsMLintAutoFixBackground'); clear('afhb','afh')

%% Variable and function colors
% com.mathworks.services.Prefs.setBooleanPref('Editor.VariableHighlighting.Automatic',ahib);
% com.mathworks.services.Prefs.setColorPref('Editor.VariableHighlighting.Color',java.awt.Color(ahi(1), ahi(2), ahi(3)));
% com.mathworks.services.ColorPrefs.notifyColorListeners('Editor.VariableHighlighting.Color'); clear('ahib','ahi')

% *** need to add checkbox option for variables with shared scope here
com.mathworks.services.Prefs.setColorPref('Editor.NonlocalVariableHighlighting.TextColor',java.awt.Color(vwss(1), vwss(2), vwss(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Editor.NonlocalVariableHighlighting.TextColor');clear('vwss')

%% Section display options
com.mathworks.services.Prefs.setBooleanPref('EditorCodepadHighVisible',hsb);
com.mathworks.services.Prefs.setColorPref('Editorhighlight-lines', java.awt.Color(hsbc(1), hsbc(2), hsbc(3)));
clear('hsb','hsbc');

%% Editor/Debugger General display options
com.mathworks.services.Prefs.setBooleanPref('Editorhighlight-caret-row-boolean',hclb);
com.mathworks.services.Prefs.setColorPref('Editorhighlight-caret-row-boolean-color',java.awt.Color(hcl(1), hcl(2), hcl(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Editorhighlight-caret-row-boolean-color'); clear('hclb','hcl')

com.mathworks.services.Prefs.setBooleanPref('EditorShowLineNumbers',slnb); clear('slnb')

com.mathworks.services.Prefs.setBooleanPref('EditorRightTextLineVisible',rtlb);
%EditorRightTextLineLimit=I80
%EditorRightTextLimitLineWidth=I1
clear('rtlb')

%% Cleanup
clear('sol')
clear('mat')

author: benhager

repo: solarized-matlab


text:
function pass = test_plus( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e4*pref.techPrefs.chebfuneps;

% Example 1
f = ballfun(ones(21,20,22));
V1 = ballfun.coeffs2vals(f.coeffs);
g = f+f;
V2 = ballfun.coeffs2vals(g.coeffs);
pass(1) = ( norm(V2(:)-2*V1(:),inf) < tol );

% Example 2
f = ballfun(@(x,y,z)x.^2.*cos(y)-1);
g = f+1;
exact = ballfun(@(x,y,z)x.^2.*cos(y));
pass(2) = norm( g - exact ) < tol;

% Example 3
f = ballfun(@(x,y,z)y.*sin(z));
g = 3+f;
exact = ballfun(@(x,y,z)y.*sin(z)+3);
pass(3) = norm( g - exact ) < tol;

% Example 4
f = ballfun(@(x,y,z)x.*sin(z).^2.*cos(y));
g = 3+f+2;
exact = ballfun(@(x,y,z)x.*sin(z).^2.*cos(y)+5);
pass(4) = norm( g - exact ) < tol;

% Example 5
f = ballfun(@(r,lam,th)1);
g = f+f;
exact = ballfun(@(r,lam,th)2);
pass(5) = norm( g - exact ) < tol;

% Example 6
f = ballfun(0) + 1;
exact = ballfun(1);
pass(6) = norm( f - exact ) < tol;

if (nargout > 0)
    pass = all(pass(:));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ A, rhs ] = boundary ( nx, ny, x, y, A, rhs )

%*****************************************************************************80
%
%% BOUNDARY sets up the matrix and right hand side at boundary nodes.
%
%  Discussion:
%
%    For this simple problem, the boundary conditions specify that the solution
%    is 10 on the left size, 100 on the right side, and 0 on the top and bottom.
%
%    Nodes are assigned a single index K, which increases as:
%
%    (NY-1)*NX+1  (NY-1)*NX+2  ...  NY * NX
%           ....         ....  ...    .....
%           NX+1         NX+2  ...   2 * NX
%              1            2  ...       NX
%
%    The index K of a node on the lower boundary satisfies:
%      1 <= K <= NX
%    The index K of a node on the upper boundary satisfies:
%      (NY-1)*NX+1 <= K <= NY * NX
%    The index K of a node on the left boundary satisfies:
%      mod ( K, NX ) = 1
%    The index K of a node on the right boundary satisfies:
%      mod ( K, NX ) = 0
%
%    If we number rows from bottom I = 1 to top I = NY
%    and columns from left J = 1 to right J = NX, then the relationship
%    between the single index K and the row and column indices I and J is:
%      K = ( I - 1 ) * NX + J
%    and
%      J = 1 + mod ( K - 1, NX )
%      I = 1 + ( K - J ) / NX
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of grid points in X and Y.
%
%    Input, real X(NX), Y(NY), the coordinates of grid lines.
%
%    Input, real sparse A(N,N), the system matrix, with the entries for the
%    interior nodes filled in.
%
%    Input, real RHS(N), the system right hand side, with the entries for the
%    interior nodes filled in.
%
%    Output, real sparse A(N,N), the system matrix, with the entries for all
%    nodes filled in.
%
%    Output, real RHS(N), the system right hand side, with the entries for
%    all nodes filled in.
%

%
%  Left boundary.
%
  j = 1;
  for i = 1 : ny
    kc = ( i - 1 ) * nx + j;
    xc = x(j);
    yc = y(i);
    A(kc,kc) = 1.0;
    rhs(kc,1) = 10.0;
  end
%
%  Right boundary.
%
  j = nx;
  for i = 1 : ny
    kc = ( i - 1 ) * nx + j;
    xc = x(j);
    yc = y(i);
    A(kc,kc) = 1.0;
    rhs(kc,1) = 100.0;
  end
%
%  Lower boundary.
%
  i = 1;
  for j = 1 : nx
    kc = ( i - 1 ) * nx + j;
    xc = x(j);
    yc = y(i);
    A(kc,kc) = 1.0;
    rhs(kc,1) = 0.0;
  end
%
%  Upper boundary.
%
  i = ny;
  for j = 1 : nx
    kc = ( i - 1 ) * nx + j;
    xc = x(j);
    yc = y(i);
    A(kc,kc) = 1.0;
    rhs(kc,1) = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mrAnatComputeVanatSpatialNorm(vAnatFileName, outFileName)
% Computes an SPM2 spatial normalization and its inverse for a vAnatomy.
%
% mrAnatComputeVanatSpatialNorm(vAnatFileName, [outFileName])
%
%
% HISTORY:
% 2005.08.05 RFD: wrote it.
%

if(~exist('vAnatFileName','var') | isempty(vAnatFileName))
    vAnatFileName = mrvSelectFile('r','dat',[],'Select a vAnatomy file...');
    if isempty(vAnatFileName), return;
    elseif ~exist(vAnatFileName,'file'), error(sprintf('%s not found.\n',vAnatFileName));
    end
end
if(~exist('outFileName','var'))
    [p,f,e] = fileparts(vAnatFileName);
    outFileName = fullfile(p,[f '_sn.mat']);
end

[img, mm] = readVolAnat(vAnatFileName);
img = uint8(img);
xform = mrAnatXform(mm,size(img),'vanatomy2acpc');

%
% Compute the normalization params.
%
disp('Computing spatial norm params...');
[sn, Vtemplate] = mrAnatComputeSpmSpatialNorm(img, xform);

%
% Convert the SN params to a deformation field.
%
disp('Computing inverse spatial norm...');
% Build the bounding box needed to sample the deformation field. We want it
% to capture the entire template, so we use the template vox-to-physical xform
% (sn.VG) to generate it.
%
% We first extract scales from the template xform. We want the
% defomation field to be sampled at this same scale. That makes pulling
% values from the inverse deformation easier and saves space, since it
% isn't necessary to sample at a higher resolution than the template.
mmTemplate = sqrt(sum(sn.VG.mat(1:3,1:3).^2));
origin  = sn.VG.mat\[0 0 0 1]';
origin  = origin(1:3)';
bb = [-mmTemplate.*(origin-1) ; mmTemplate.*(sn.VG.dim(1:3)-origin)];
% Now compute the deformation:
d = single(mrAnatSnToDeformation(sn, mmTemplate, bb));

% Invert the deformation. The inverse deformation is essentially a look-up
% table that maps voxel coords to the physical space of the template. This
% is exactly what we want.
%
% 4th arg is 4x4 xform from mm to voxels in the coordinate frame of the inverse deformation field
% 5th arg is xform from voxels to mm in the coordinate frame of the forward deformation field
[defX,defY,defZ] = spm_invdef(d(:,:,:,1), d(:,:,:,2), d(:,:,:,3), sn.VF.dim(1:3), ...
                             inv(sn.VF.mat), sn.VG.mat);

% We can save it as int8, since real brain templates never extend beyond
% 128mm from the origin. But, we'll check, just to be sure.
if(all(abs(defX(:))<128) & all(abs(defY(:))<128) & all(abs(defZ(:))<128))
    voxToTemplateLUT = zeros([3 size(defX)], 'int8');
    voxToTemplateLUT(1,:,:,:) = int8(round(defX));
    voxToTemplateLUT(2,:,:,:) = int8(round(defY));
    voxToTemplateLUT(3,:,:,:) = int8(round(defZ));
else
    voxToTemplateLUT = zeros([3 size(defX)], 'int16');
    voxToTemplateLUT(1,:,:,:) = int16(round(defX));
    voxToTemplateLUT(2,:,:,:) = int16(round(defY));
    voxToTemplateLUT(3,:,:,:) = int16(round(defZ));
end

fprintf('Saving to %s...\n', outFileName);
save(outFileName, 'sn', 'voxToTemplateLUT');
disp('done.');
return;

% DEBUGGING CODE:

mmPerVox = [1 1 1];

% Get a bounding box in ac-pc space that captures the whole template image.
templateMmPerVox = sqrt(sum(sn.VG(1).mat(1:3,1:3).^2));
if det(sn.VG.mat(1:3,1:3))<0, mmPerVox(1) = -mmPerVox(1); end;
templateOrigin  = sn.VG.mat\[0 0 0 1]';
templateOrigin  = templateOrigin(1:3)';
bb = [-templateMmPerVox .* (templateOrigin-1) ; templateMmPerVox.*(sn.VG(1).dim(1:3)-templateOrigin)];
% The code in mrAnatResliceSpm needs to know how to transform the image
% coords from the template's ac-pc space to actual deformation-image space.
% The following was gleaned from spm_reslice_sn. I only vaguely understand
% it, but it seems to work.
og  = -templateMmPerVox .* templateOrigin;
M1  = [templateMmPerVox(1) 0 0 og(1) ; 0 templateMmPerVox(2) 0 og(2) ; 0 0 templateMmPerVox(3) og(3) ; 0 0 0 1];
outMmPerVox = [1 1 1];
of  = -outMmPerVox.*(round(-bb(1,:)./outMmPerVox)+1);
M2  = [outMmPerVox(1) 0 0 of(1) ; 0 outMmPerVox(2) 0 of(2) ; 0 0 outMmPerVox(3) of(3) ; 0 0 0 1];
d.inMat = inv(sn.VG(1).mat*inv(M1)*M2);

dField = mrAnatSnToDeformation(sn, [1 1 1], bb);
d.deformX = dField(:,:,:,1);
d.deformY = dField(:,:,:,2);
d.deformZ = dField(:,:,:,3);
d.outMat = inv(xform);
[imgSn, newXform] = mrAnatResliceSpm(img, d, bb, mmPerVox, [7 7 7 0 0 0]);
imgSn(isnan(imgSn)) = 0;

figure; imagesc(makeMontage(imgSn)); axis image; colormap(gray)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [imdepthMin, imdepthMax, imdepthCol, contact, x, y, z] = makeDepthImages(bndinfo)

global DO_DISPLAY;
DO_DISPLAY = 0;

load '/home/dhoiem/src/cmu/iccv07Final/data/contactdt.mat';

[tmp, glabels] = max(bndinfo.result.geomProb, [], 2);
glabels((glabels>=2) & (glabels<=4)) = 2;
glabels(glabels==5) = 3;

lab = bndinfo.edges.boundaryType;

lab = lab(1:end/2) + 2*lab(end/2+1:end);

[imdepthMin, imdepthMax, imdepthCol, contact, x, y, z] = ...
    getDepthRangeForDisplay(bndinfo, glabels, lab, contactdt, 0.5);

author: Cloud-CV

repo: object-proposals


text:
function [y,ysup] = erfc_rnd(x,rnd)
% input x real non-negative column vector
% rnd  -1  y = lower bound for erf(x)
%       1  y = upper bound for erf(x)
%      []  [y,ysup] inclusion of erf(x)
% rounding may be altered after leaving erfc_rnd

% written  05/30/13     S.M. Rump
%

  xmax = hex2num('403b369a6244e684');   % ~27.2: erfc(x)<subrealmin fuer x>=xmax

  y = x;
  if isempty(rnd)
    ysup = x;
  end

  index = ( x<-0.5 );                   % Use 2-erfc(-x)
  if any(index)                         % x in [-inf,-0.5)
    if isempty(rnd)
      [yindex,ysupindex] = erfc_rnd(-x(index),-rnd);
      setround(-1)
      y(index) = 2 - ysupindex;
      setround(1)
      ysup(index) = 2 - yindex;
    else
      yindex = erfc_rnd(-x(index),-rnd);
      setround(rnd)
      y(index) = 2 - yindex;
    end
  end
  Index = index;                        % store finished indices

  index = ( ~Index ) & ( x<0 );         % first method
  if any(index)                         % x in [-0.5,0)
    if isempty(rnd)
      [yindex,ysupindex] = erf_rnd1(-x(index),rnd,6);
      setround(-1)
      y(index) = 1 + yindex;
      setround(1)
      ysup(index) = 1 + ysupindex;
    else
      yindex = erf_rnd1(-x(index),rnd,6);
      setround(rnd)
      y(index) = 1 + yindex;
    end
  end
  Index = Index | index;                % store finished indices

  index = ( ~Index ) & ( x<=0.5 );      % first method
  if any(index)                         % x in [0,0.5]
    if isempty(rnd)
      [yindex,ysupindex] = erf_rnd1(x(index),rnd,6);
      setround(-1)
      y(index) = 1 - ysupindex;
      setround(1)
      ysup(index) = 1 - yindex;
    else
      yindex = erf_rnd1(x(index),-rnd,6);
      setround(rnd)
      y(index) = 1 - yindex;
    end
  end
  Index = Index | index;                % store finished indices

  index = ( ~Index ) & ( x<=7 );        % second method
  if any(index)                         % x in (0.5,7]
    yindex = erfc_rnd2(x(index));
    if rnd==-1
      y(index) = yindex.inf;
    elseif rnd==1
      y(index) = yindex.sup;
    else
      y(index) = yindex.inf;
      ysup(index) = yindex.sup;
    end
  end
  Index = Index | index;                % store finished indices

  index = ( ~Index ) & ( x<=10 );       % third method
  if any(index)                         % x in (7,10]
    if isempty(rnd)
      [y(index),ysup(index)] = erfc_rnd3(x(index),rnd,16);
    else
      y(index) = erfc_rnd3(x(index),rnd,16);
    end
  end
  Index = Index | index;                % store finished indices

  index = ( ~Index ) & ( x<=xmax );     % third method
  if any(index)                         % x in (10,xmax]
    if isempty(rnd)
      [y(index),ysup(index)] = erfc_rnd3(x(index),rnd,10);
    else
      y(index) = erfc_rnd3(x(index),rnd,10);
    end
  end

  index = ( x>xmax );
  if any(index)                         % x in (xmax,inf]
    if isempty(rnd)                     % inclusion [y,ysup] of erfc(x)
      y(index) = 0;
      ysup(index) = hex2num('0000000000000001');    % subrealmin
    elseif rnd==1                       % y upper bound for erfc(x)
      y(index) = hex2num('0000000000000001');    % subrealmin
    else                                % y lower bound for erfc(x)
      y(index) = 0;
    end
  end

repo: OpenMAS


text:
% Demo for spm_rar; robust autogregressive modelling

close all

secs=3;
ns=128;
t=[1/ns:1/ns:secs];
N=length(t);

% Generate Gaussian-Mixture Noise
m=2;
mix.m=m;
mix.state(1).prior=0.9;
mix.state(2).prior=0.1;
mix.state(1).m=0;
mix.state(2).m=0;
mix.state(1).C=1;
mix.state(2).C=100;
[noise,gamma_true]=spm_samp_mix(mix,N);
new_index=randperm(N);
noise=noise(new_index);
gamma_true=gamma_true(new_index);

% Convolve noise with AR filter
a_true=[-1.8517,1.3741,0.1421,-0.6852,0.3506];
p_true=length(a_true);
y=filter(1,[1,a_true],noise);
y=y(1:N);

% Fit AR model
ar=spm_ar(y,p_true);

% Fit Robust AR model
[rar,yclean] = spm_rar(y,p_true,m);

rar3=spm_rar(y,p_true,3);
rar4=spm_rar(y,p_true,4);

% Get model evidences
fm=[ar.fm,rar.fm,rar3.fm,rar4.fm];
fm=fm-mean(fm);

figure
subplot(2,2,1);
plot(t,y);
title('Data');
subplot(2,2,3);
hist(noise,20);
title('Noise histogram');
subplot(2,2,2);
[tmp,outlier]=min(rar.pi);
standard=m+1-outlier;
plot(gamma_true);
title('Outliers');
axis([0 N -0.1 1.1]);
subplot(2,2,4);
bar(fm);
title('Model evidence');
xlabel('m');

disp(' ');
disp('OUTLIER DETECTION:');
pos_prob=rar.gamma(outlier,find(gamma_true==1));
sens=length(find(pos_prob>0.5))/length(pos_prob);

neg_prob=rar.gamma(standard,find(gamma_true==0));
spec=length(find(neg_prob>0.5))/length(neg_prob);

disp(sprintf('Proportion of outliers correctly detected = %1.2f',sens));
disp(sprintf('Proportion of standards correctly detected = %1.2f',spec));

disp(' ');
disp('ACCURACY OF AR COEFFICIENT ESTIMATION:');
d_ar=norm(ar.a_mean-a_true');
d_rar=norm(rar.posts.a_mean-a_true');
disp(sprintf('Error for AR=%1.3f',d_ar));
disp(sprintf('Error for RAR=%1.3f',d_rar));
disp(sprintf('Ratio E_RAR/E_AR=%1.3f',d_rar/d_ar));

% figure;plot(y);hold on;plot(yclean,'r');
% legend('Original','Clean');

author: spm


text:
function cm = render_on_surface(obj, surface_handles, varargin)
% Map voxels in obj to isosurface (patch) handles in han and change surface colors according to values in obj
%
% :Usage:
% ::
%
%     cm = render_on_surface(obj, han, [optional inputs])
%
% - Object can be thresholded or unthresholded statistic_image, or other image_vector object
% - Uses only the first image in obj
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2020 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **obj:**
%       - image_vector object or subclasses
%       - e.g., fmri_data, thresholded or unthresholded statistic_image obj
%       - Uses only the first image in obj
%
%   **han:**
%       - A vector of one or more handles to patch objects (from surfaces)
%       - e.g., from patch(isosurface()), patch(isocaps()),
%       fmri_data.surface(), fmri_data.isosurface()
%       - Try: surface_handles = addbrain('left_cutaway');
%       - Try: surface_handles = addbrain('coronal_slabs_4');
%
% :Optional Inputs:
%   **'clim':**
%        Followed by [lowerlimit upperlimit] values to define color limits
%
%   **'colormap':**
%        Followed by the name of a Matlab colormap string
%        e.g., 'hot' (default), 'summer', 'winter'
%        - Creates a split colormap with gray values where object values
%        are zero, and colors where object values are nonzero
%
%   **'indexmap':**
%       Interpret data as indices into colormap. You must specify
%       'colormap' argument as well, otherwise this does nothing and a
%       warning is thrown.
%
% :Outputs:
%
%   **renders colors on surfaces:**
%   Note that it is not always possible to reset the surfaces to their
%   original colors completely; you may have to redraw them if you change
%   the threshold/voxels that should be colored.
%
% :Examples:
% ::
%
% % Load a standard dataset (emotion regulation test data)
% % Perform a t-test, and render the results on a series of coronal slabs
%
% imgs = load_image_set('emotionreg', 'noverbose');
% t = ttest(imgs, .01, 'unc');
% figure; han = addbrain('coronal_slabs_4');
%
% % When the t map has positive and negative values, creates a special
% % bicolor split colormap that has warm colors for positive vals, and cool colors
% % for negative vals:
% render_on_surface(t, han);
%
% % You can set the color limits (here, in t-values)
% render_on_surface(t, han, 'clim', [-4 4]);
%
% % You can set the colormap to any Matlab colormap:
% render_on_surface(t, han, 'colormap', 'summer');
%
% % If your image is positive-valued only (or negative-valued only), a
% % bicolor split map will not be created:
%
% t = threshold(t, [2 Inf], 'raw-between');
% render_on_surface(t, han, 'colormap', 'winter', 'clim', [2 6]);
%
% Note:
% To erase rendered blobs, use:
% sh = addbrain('eraseblobs', sh);

% Programmers' Notes:
% Tor Wager, Jan 2020
% desirable features:
% - Speed
% - control of colormap, can update
% - consistent value->color mapping across surfaces
% - need to know meaning of values, as they often represent statistic values
% - can re-run without changing colors
% - need to plot on isosurface or isocaps objects that use colormap to
% render gray-scale images.
% This last part makes it a bit more complicated.
% Rendering only on solid-color patch objects from isosurfaces is easier:
%     c = isocolors(mesh_struct.X, mesh_struct.Y, mesh_struct.Z, mesh_struct.voldata, han(i));
%     set(han(i), 'FaceColor', 'interp')
%     c_orig = get(han(i), 'FaceVertexCData');
%     han(i).FaceVertexCData = c;
%
% But to render on either solid-color or gray scale colormapped patches, we want to use
% a split colormap.
%
% Different surface handles (to patches) will have FaceColor is an rgb
% triplet (for solid-color patches) or 'interp' for colormapped patches
% (like isocaps). This will be set to interp for all patches, but we have
% to preserve the existing colormapped grayscale values for interp patches
% by using the split colormap.
%
% Sept 2020: add solid-color 'color' option
% revise surface() method to use this

if any(~ishandle(surface_handles))
    error('Some surface_handles are not valid handles');
end

% -------------------------------------------------------------------------
% OPTIONAL INPUTS
% -------------------------------------------------------------------------

% cm = [];
nvals = 256;                % number of values for each segment of colormap
colormapname = 'hot';       % default colormap name (all Matlab defined colormaps OK)
custom_colormap = false;
pos_colormap = [];
neg_colormap = [];
clim = [];
axis_handle = get(surface_handles, 'Parent');          % axis handle to apply colormap to; can be altered with varargin
dolegend = true;
doindexmap = false;

allowable_keyword_value_pairs = {'clim' 'color' 'colormap' 'colormapname' 'axis_handle' 'pos_colormap' 'neg_colormap'};

% optional inputs with default values - each keyword entered will create a variable of the same name

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'color'
                % Special instructions for solid colors
                mycolor = varargin{i+1};
                validateattributes(mycolor, {'numeric'}, {'>=', 0, '<=', 1, 'size', [1 3]})

                nvals = 256;
                pos_colormap = repmat(mycolor, nvals, 1);

            case 'colormap'

                colormapname = varargin{i+1}; varargin{i+1} = [];
                custom_colormap = true;

            case 'cmaprange'
                % option to match region.montage method

                clim =  varargin{i+1}; varargin{i+1} = [];

            case 'nolegend'

                dolegend = false;

            case 'indexmap'
                if ~contains('colormap',varargin(cellfun(@ischar,varargin))),...
                        warning('render_on_surface() given ''indexmap'' argument, but no ''colormap'' argument. ''indexmap'' doesn''t do anything without a ''colormap'' argument.');
                end

                doindexmap = true;

            case allowable_keyword_value_pairs

                eval([varargin{i} ' = varargin{i+1}; varargin{i+1} = [];']);

                % eliminate this here because we may have passed many irrelevant values in in surface() call
                %otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% There are 2 ways to enter custom colormaps.
% 1. By name
% 'colormap', [name of Matlab colormap]
%
% 2. Enter colormaps generated externally
% 'pos_colormap', followed by [n x 3] color vectors and/or
% 'neg_colormap', followed by [n x 3] color vectors
%
% If you enter 'pos_colormap' / 'neg_colormap', these will be stacked and
% added to a gray colormap

% Get range of values in object
% -----------------------------------------------------------------------
% These will be used to map color limits and create legend
% we will map values into the range of [1 nvals] where nvals = colormap length
% Deal with possibility of multiple images in object
% Exclude zeros

[datvec, clim] = get_data_range(obj, clim);

% -------------------------------------------------------------------------
% Define colormap
% -------------------------------------------------------------------------

if ~isempty(pos_colormap) ||  ~isempty(neg_colormap)

    if any(strcmp(varargin, 'colormap'))
        error('Entering ''colormap'' argument and name cannot be done when entering custom colormaps in ''pos_colormap''/''neg_colormap''');
    end

    if ~isempty(pos_colormap) && ~isempty(neg_colormap) && size(pos_colormap, 1) ~= size(neg_colormap, 1)
        error('Lengths of positive-color and negative-color colormaps must match.')
    end

    if ~isempty(pos_colormap)
        validateattributes(pos_colormap, {'numeric'}, {'>=', 0, '<=', 1, 'size', [NaN 3]})
        nvals = size(pos_colormap, 1);
    end

    if ~isempty(neg_colormap)
        validateattributes(neg_colormap, {'numeric'}, {'>=', 0, '<=', 1, 'size', [NaN 3]})
        nvals = size(neg_colormap, 1);
    end

    if isempty(pos_colormap)
        pos_colormap = hot(nvals);
    end

    if isempty(neg_colormap)
        neg_colormap = cool(nvals);
    end

    custom_colormap = true;

    % Build the colormap
    % -----------------------------------------------------------------------
    % Skip colormap generator, already
    % found the range of colors for pos/neg values to split around 0 here.
    % this sets the colormap for axis_handle
    [cm, kpos, kneg] = hotcool_split_colormap(nvals, clim, axis_handle, pos_colormap(1, :), pos_colormap(end, :), neg_colormap(1, :), neg_colormap(end, :));

    % colormapname = [neg_colormap; pos_colormap];   % colormapname is either name or [nvals x 3] matrix
    % nvals = size(colormapname, 1);                 % needs to match for color and gray maps to work right

else

    % Build the colormap
    % -----------------------------------------------------------------------
    % See notes below

    if doindexmap

        nvals = min([256, length(unique(obj.dat))]);

    end

    if custom_colormap

        [cm, kpos, kneg] = split_colormap(nvals, colormapname, axis_handle); % colormapname is either name or [nvals x 3] matrix

    elseif diff(sign(clim))

        % Default colormap for objects with - and + values
        [cm, kpos, kneg] = hotcool_split_colormap(nvals, clim, axis_handle);

    else
        [cm, kpos, kneg] = split_colormap(nvals, colormapname, axis_handle);
    end

end % custom posneg or other colormap

% -------------------------------------------------------------------------
% Change colors
% -------------------------------------------------------------------------

% Map object data values to indices in split colormap
% -----------------------------------------------------------------------
% Define mapping function from image values -> colormap indices (1-512).
% We will apply this later to vertex color data (only colored)
% Defining this here preserves the same mapping across different surfaces
% clim(1) mapped to lowest color, clim(2) to highest color
% map_function = @(c) 1 + nvals + (c - clim(1)) ./ range(clim) .* nvals;

map_function = @(c, x1, x2, y1, y2)  y1 + (c - x1) * (y2 - y1) ./ (x2 - x1);

% problem with above is that border gets mapped to low color when interpolating, which is
% blue for split colormap, creating blue outline around orange blobs. SO
% we need to map separately for pos and neg colored blobs
% ***

% problem with below is that vertices near empty (zero) voxels get interpolated down to zero
%map_function = @(c) 1 + nvals + (c - 0) ./ range(clim) .* nvals;

% Reconstruct volume data
% -----------------------------------------------------------------------
% Needed to map to vertices later
[~, ~, mesh_struct] = reconstruct_image(obj);                % get volume data for slices

% Deal with edge interpolation effects
% problem is that vertices near empty (zero) voxels get interpolated down to zero
% solution: map all non-zero voxels to lowest limit value
% Need to preserve signs as well.
% Then replace empty voxels with gray
% lower lim to interp to 10% of the way towards zero,
% preserving sign
% lowerclimit = linspace(clim(1), 0, 10);
% lowerclimit = lowerclimit(2);
% mesh_struct.voldata(mesh_struct.voldata == 0) = lowerclimit;
% lowerclimit = 0; ... wh = c == 0 | isnan(c);

% Kludgy fix for interpolation issues around border of blobs
% Changing c or c_colored won't do it because of how gray and colored
% colormaps are stacked. With one color, we interpolate into the gray zone,
% into about the top 256-150 values.
% for i = 256:-1:150
%     cm(i, :) = [.5 .5 .5];
% end
% colormap(cm)

for i = 1:length(surface_handles)
    % For each surface handle entered
    % -----------------------------------------------------------------------

    % Get vertex colors for image (obj) to map
    % -----------------------------------------------------------------------
    % Use isocolors to get colormapped values
    % Convert to index color in split colormap by adding nvals to put in
    % colored range
    % isocolors returns nans sometimes even when no NaNs in data
    c = isocolors(mesh_struct.X, mesh_struct.Y, mesh_struct.Z, mesh_struct.voldata, surface_handles(i));
    if doindexmap, c = round(c); end

    % doesn't work to fix interpolation error.
    %     border_percent = prctile(abs(c(c ~= 0)), 80);
    %     whlow = abs(c) < border_percent | isnan(c);

    c_colored = c;

    whpos = c > 0;
    %     kpos = 61;   % which block of 256 colors; depends on colormap
    cpos = map_function(c(whpos), 0, clim(2), (kpos-1)*nvals+1, kpos*nvals); % map into indices in hot cm range of colormap
    c_colored(whpos) = cpos;

    whneg = c < 0;
    %     kneg = 55;   % which block of 256 colors
    cneg = map_function(c(whneg), clim(1), 0, (kneg-1)*nvals+1, kneg*nvals); % map into indices in cool cm range of colormap
    c_colored(whneg) = cneg;

    wh = c == 0 | isnan(c);                    % save these to replace with gray-scale later

    %c_colored = map_function(c);    % Map to colormap indices (nvals = starting range, nvals elements)

    % FIX: rescale range so we don't map into gray range at edges due to
    % interpolation
    %     xx = c_colored(~wh); % in-blob vertices
    %     xx = 1 + nvals + max(xx) .* (xx - min(xx)) ./ (range(xx) + nvals);
    %     c_colored(~wh) = xx;
    %cm(256, :) = [1 0 1];

    % Get the original grayscale values
    % -----------------------------------------------------------------------
    % Exclude colored values

    c_gray = get(surface_handles(i), 'FaceVertexCData');

    if isempty(c_gray)

        % solid
        % c_gray = repmat(round(nvals ./ 2), size(get(surface_handles(i), 'Vertices'), 1), 1);
        c_gray = repmat(0 .* nvals + .5 .* round(nvals), size(get(surface_handles(i), 'Vertices'), 1), 1);

    else
        %orig: c_gray(~wh) = (c_gray(~wh) - min(c_gray(~wh))) ./ range(c_gray(~wh)) .* nvals;

        if range(c_gray(wh)) == 0
            % Solid surface color
            % do nothing - skip rescaling
            % -----------------------------------------------------------------------
            %c_colored(wh) = c_gray(wh);

        else
            % rescale to gray part of colormap
            c_gray(wh) = (c_gray(wh) - min(c_gray(wh))) ./ range(c_gray(wh)) .* nvals;
            % ***

            % Now merge graycsale and object-mapped colors
            % -----------------------------------------------------------------------
            % Everything that was originally zero in isocolors c
            % gets replaced with its original grayscale values in c_gray

            %c_colored(wh) = c_gray(wh);

        end
    end

    c_colored(wh) = c_gray(wh);

    % Set object handle properties
    % -----------------------------------------------------------------------
    % Change from 'scaled' to 'direct' mode

    % Save FaceVertexCData f or later
    prevdata = get(surface_handles(i), 'FaceVertexCData');
    set(surface_handles(i), 'UserData', prevdata);

    set(surface_handles(i), 'FaceVertexCData', c_colored);  % dot indexing sometimes works, sometimes doesn't...depends on handle type
    if ~doindexmap
        set(surface_handles(i), 'FaceColor', 'interp');
    else
        set(surface_handles(i), 'FaceColor', 'flat');
    end
    set(surface_handles(i), 'CDataMapping', 'direct')
    set(surface_handles(i), 'EdgeColor', 'none');

end

% Colorbars - legend
% -----------------------------------------------------------------------

if ~dolegend, return, end

if any(datvec > 0)

    bar1axis = axes('Position', [.55 .55 .38 .4]);
    colormap(bar1axis, cm(1+(kpos-1)*nvals:kpos*nvals, :))
    colorbar_han = colorbar(bar1axis);
    set(bar1axis, 'Visible', 'off');

    minpos = min(datvec(datvec > 0));
    set(colorbar_han, 'YTick', [0 1], 'YTickLabel', round([minpos clim(2)] * 100)/100, 'FontSize', 18);

end

if any(datvec < 0)

    bar2axis = axes('Position', [.55 .1 .38 .4]);
    colormap(bar2axis, cm(1+(kneg-1)*nvals:kneg*nvals, :))
    colorbar_han = colorbar(bar2axis);
    set(bar2axis, 'Visible', 'off');

    maxneg = max(datvec(datvec < 0));
    set(colorbar_han, 'YTick', [0 1], 'YTickLabel', round([clim(1) maxneg] * 100)/100, 'FontSize', 18);

end

end % main function

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
% Subfunctions
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function [datvec, clim] = get_data_range(obj, clim)

if isa(obj, 'statistic_image')

    sig = logical(obj.sig);
    dat = obj.dat(sig);

else

    wh = obj.dat ~= 0 & ~isnan(obj.dat);
    dat = obj.dat(wh);

end

datvec = dat(:);

if any(isinf(datvec))
    warning('Some image values are Inf. Expect erratic behavior/errors.');
    whinf = isinf(datvec);
    datvec(whinf) = sign(datvec(whinf)) .* max(abs(datvec(~whinf)));
end

if isempty(clim)

    %clim = [min(datvec) max(datvec)];  % clim: data values for min and max, should become min/max colors

    if any(datvec < 0) && any(datvec > 0) % split colormap

        clim = double([prctile(datvec(datvec < 0), 10) prctile(datvec(datvec > 0), 90)]); % Match defaults for montage in render_blobs

    else
        % may need adjustment
        clim = [min(datvec) max(datvec)];

    end

    % if no variance, we have constant data values - special case.
    % reducing lower clim(1) will effectively map all data to max color value
    if abs(diff(clim)) < 100 * eps
        clim(1) = clim(2) - 1;
    end

end

end % function

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
% Colormap subfunctions
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% - This will include nvals (generally 256) gray-scale and nvals colored
% values. If nvals = 256, the first 256 colors will be grayscale, and the
% next 256 color values.
% - If the object has positive and negative values, assume we want a split
% colormap showing negative values in cool colors, positive values in hot
% colors. The split colormap will also have different ranges for negative and
% positive values
% - If object is positive-valued or negative-valued only, assume we want a unipolar colormap

function [cm, kpos, kneg] = hotcool_split_colormap(nvals, clim, axis_handle, varargin)
%
% cm = hotcool_split_colormap(nvals, clim, [lowhot hihot lowcool hicool]), each is [r g b] triplet
%
% Get hotcool split colormap (around 0)
% note: colormaps with these colors match spm_orthviews_hotcool_colormap used in orthviews()
% [.4 .4 .3], [1 1 0] (hot/pos) and [0 0 1], [.4 .3 .4] (cool/neg)

% Split blue/yellow
% lowhot = [.4 .4 .3];
% hihot = [1 1 0];
% lowcool = [0 0 1];
% hicool = [.4 .3 .4];

% Default addblobs - orange/pink
hihot = [1 1 0]; % max pos, most extreme values
lowhot = [1 .4 .5]; % [.8 .3 0]; % min pos
hicool = [0 .8 .8]; % [.3 .6 .9]; % max neg
lowcool = [0 0 1]; % min neg, most extreme values

graybuffer = 20;  % for border - fade to gray in lowest k values

if ~isempty(varargin)
    if length(varargin) < 4
        error('Enter no optional arguments or 4 [r g b] color triplets');
    end
    lowhot = varargin{1};
    hihot = varargin{2};
    lowcool = varargin{3};
    hicool = varargin{4};
end

cmgray = gray(nvals);

thr = 0;    % Threshold to split colormap into two colors

% vec = linspace(clim(1), clim(2), nvals);
% pos = find(vec > thr);
% np = length(pos);

[hotcm, coolcm] = deal([]);

if graybuffer

    hotcm = [colormap_tor([.5 .5 .5], lowhot, 'n', graybuffer); ...
        colormap_tor(lowhot, hihot, 'n', nvals-graybuffer)];

else

    hotcm = colormap_tor(lowhot, hihot, 'n', nvals); %np);

end

%hotcm(1, :) = [.5 .5 .5]; % first element is gray
% end

% neg = find(vec < -thr);
% nn = length(neg);

% if nn > 0
%coolcm(end, :) = [.5 .5 .5]; % last element is gray
% end

if graybuffer

    coolcm = [colormap_tor(lowcool, hicool, 'n', nvals-graybuffer); ...
        colormap_tor(hicool, [.5 .5 .5], 'n', graybuffer)];

else

    coolcm = colormap_tor(lowcool, hicool, 'n', nvals); % nn);

end

% for interpolation issues
% 1:256 is gray-scale anat; 257:512 is buffer; 513:768 is negcm; 769:1024
% is buffer; 1025:1280 is poscm
buffercm = repmat([.5 .5 .5], nvals, 1);

% which block of nvals (usually 256) has colored values
kpos = 61;
kneg = 55;

%cm = [coolcm; buffercm; buffercm; cmgray; hotcm];
cm = [cmgray; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    buffercm; buffercm; buffercm; buffercm; coolcm; ...
    buffercm; buffercm; buffercm; buffercm; buffercm; ...
    hotcm];

% Kludgy fix for interpolation issues
% cm(128+3*256:4*256, :) = .5;
% cm(1+3*256:50+3*256, :) = .5;

% cm(128+0*256:1*256, :) = .5;
% cm(1+0*256:50+0*256, :) = .5;

% cm = [cmgray; cm];

% Change colormap(s)
if isempty(axis_handle)
    colormap(cm)

elseif iscell(axis_handle)

    for i = 1:length(axis_handle)
        colormap(axis_handle{i}, cm);
    end

else
    colormap(axis_handle, cm);
end

end

function [cm, kpos, kneg] = split_colormap(nvals, cmcolor, axis_handle)
% Get split colormap with gray then colored values
% e.g., cm = split_colormap(256, [2 8], 'hot')
%
% colormapname is either name or [nvals x 3] matrix

%validateattributes(cmcolor, {'numeric' 'char'})

% which block of nvals (usually 256) has colored values
kpos = 12;
kneg = 12;

cmgray = gray(nvals);

if ischar(cmcolor)
    % Named colormap function
    cmcolor = eval(sprintf('%s(%d)', cmcolor, nvals));
end

validateattributes(cmcolor, {'numeric'}, {'>=', 0, '<=', 1, 'size', [NaN 3]})

buffercm = repmat([.5 .5 .5], nvals, 1);

cm = [cmgray; repmat(buffercm, 10, 1); cmcolor];

% Change colormap(s)
if isempty(axis_handle)
    colormap(cm)

elseif iscell(axis_handle)

    for i = 1:length(axis_handle)
        colormap(axis_handle{i}, cm);
    end

else
    colormap(axis_handle, cm);
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [soln,eqn,info] = StokesP2P0(node,elem,bdFlag,pde,option)
%% STOKESP2P0 Stokes equation: P2-P0 elements.
%
%   [soln,eqn,info] = STOKESP2P0(node,elem,bdFlag,pde) use quadratic and piceswise
%   constant elements to approximate velocity u and pressure p, repectively.
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%   with
%       Dirichlet boundary condition        u = g_D  on \Gamma_D,
%       Neumann boundary condition du/dn - np = g_N  on \Gamma_N.
%
%   It is a choice of option.fem in Stokes. Please read Stokes for more
%   information on the input and output.
%
% See also Stokes, Poisson, StokesP2P1
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end

%% Construct Data Structure
[elem2dof,edge,bdDof] = dofP2(elem);
N = size(node,1);  NT = size(elem,1);  Nu = N+size(edge,1);   Np = NT;

t = cputime;
%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix for Laplace operator
% generate sparse pattern
ii = zeros(21*NT,1); jj = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        index = index + NT;
    end
end
% quadrature points
if ~isfield(pde,'nu'), pde.nu = []; end
if ~isfield(option,'quadorder')
    % constant viscosity
    option.quadorder = 2;        % default order
    if ~isempty(pde.nu) && isnumeric(pde.nu) % numerical viscosity
        option.quadorder = 3;    % exact for linear diffusion coefficient
    end
end
[lambda, w] = quadpts(option.quadorder);
nQuad = size(lambda,1);
% compute non-zeros
sA = zeros(21*NT,nQuad);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    index = 0;
    for i = 1:6
        for j = i:6
            Aij = 0;
            if isempty(pde.nu) || isnumeric(pde.nu)
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
            else
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
            end
            if ~isempty(pde.nu) && (pde.nu~=1)
                Aij = pde.nu*Aij;
            end
            Aij = Aij.*area;
            sA(index+1:index+NT,p) = Aij;
            index = index + NT;
        end
    end
end
sA = sum(sA,2);
% assemble the matrix
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Nu,Nu);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Nu,Nu);
A = A + AU + AU';
A = blkdiag(A,A);
clear Aij ii jj sA Dphip

%% Assemble matrix for divergence operator
% Since Dphi is linear, 1-pt quadrature is exact for the integral
% int(div(phi)). We evaluate each basis at the barycenter.
Dphic(:,:,6) = 4/3*(Dlambda(:,:,1) + Dlambda(:,:,2));
Dphic(:,:,1) = 1/3*Dlambda(:,:,1); % (4*lambda(p,1)-1).*Dlambda(:,:,1);
Dphic(:,:,2) = 1/3*Dlambda(:,:,2);
Dphic(:,:,3) = 1/3*Dlambda(:,:,3);
Dphic(:,:,4) = 4/3*(Dlambda(:,:,2) + Dlambda(:,:,3));
Dphic(:,:,5) = 4/3*(Dlambda(:,:,3) + Dlambda(:,:,1));
clear Dlambda
% divergence matrix
Dx = sparse(Np,Nu);
Dy = sparse(Np,Nu);
for i = 1:6
    Dx = Dx + sparse(1:NT,double(elem2dof(:,i)),Dphic(:,1,i).*area,Np,Nu);
    Dy = Dy + sparse(1:NT,double(elem2dof(:,i)),Dphic(:,2,i).*area,Np,Nu);
end
B = -[Dx Dy]; %#ok<*NASGU>
clear Dphi Dx Dy

%% Assemble right hand side
f1 = zeros(Nu,1);
f2 = zeros(Nu,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,weight] = quadpts(option.fquadorder);
    % basis values at quadrature points
    phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
    phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
    phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
    phi(:,4) = 4*lambda(:,2).*lambda(:,3);
    phi(:,5) = 4*lambda(:,3).*lambda(:,1);
    nQuad = size(lambda,1);
    ft1 = zeros(NT,6);
    ft2 = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        % function values at quadrature points
        fp = pde.f(pxy);
        for j = 1:6
            ft1(:,j) = ft1(:,j) + fp(:,1).*phi(p,j)*weight(p);
            ft2(:,j) = ft2(:,j) + fp(:,2).*phi(p,j)*weight(p);
        end
    end
    ft1 = ft1.*repmat(area,1,6);
    ft2 = ft2.*repmat(area,1,6);
    f1 = accumarray(elem2dof(:),ft1(:),[Nu 1]);
    f2 = accumarray(elem2dof(:),ft2(:),[Nu 1]);
end

%% Boundary condition
[AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesP2P0;

%% Record assembeling time
assembleTime = cputime - t;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(ufreeDof), return; end
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if length(f)+length(g) <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else          % Multigrid-type  solver for large size systems
        option.solver = 'asmg';
    end
end
solver = option.solver;

%% Solver
switch solver
    case 'direct'
        t = cputime;
        bigA = [AD, BD'; ...
                BD, sparse(Np,Np)];
        bigF = [f; g];
        bigu = [u; p];
        bigFreeDof = [ufreeDof; 2*Nu+pDof];
        bigu(bigFreeDof) = bigA(bigFreeDof,bigFreeDof)\bigF(bigFreeDof);
        u = bigu(1:2*Nu);
        p = bigu(2*Nu+1:end);
        residual = norm(bigF - bigA*bigu);
        info = struct('solverTime',cputime - t,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
        option.solver  = 'WCYCLE';
        [u(ufreeDof),p,info] = mgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                        u(ufreeDof),p,elem,ufreeDof,option);
    case 'asmg'
        [u(ufreeDof),p,info] = asmgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                          u,p,node,elem,bdFlag,ufreeDof,option);
end

%% Post-process
if length(pDof) ~= Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(p.*area)/sum(area);
    p = p - c;
end

%% Output
soln = struct('u',u,'p',p);
eqn = struct('A',AD,'B',BD,'Lap',A,'f',f,'g',g,...
             'edge',edge,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdStokesP2P0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesP2P0
    %% Boundary condition of Stokes equation: P2-P0 elements

    %% Initial set up
%     f = [f1; f2];    % set in Neumann boundary condition
    g = zeros(Np,1);
    u = zeros(2*Nu,1);
    p = zeros(Np,1);
    ufreeDof = (1:Nu)';
    pDof = (1:Np)';

    if ~exist('bdFlag','var'), bdFlag = []; end
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Find Dirichlet dof and modify the matrix
    % Find Dirichlet boundary dof: fixedDof and pDof
    isFixedDof = false(Nu,1);
    if ~isempty(bdFlag)       % case: bdFlag is not empty
        elem2edge = elem2dof(:,4:6)-N;
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;   % nodes of all D-edges
        isFixedDof(N + find(isDirichlet')) = true;% dof on D-edges
        fixedDof = find(isFixedDof);
        ufreeDof = find(~isFixedDof);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        fixedDof = bdDof;
        isFixedDof(fixedDof) = true;
        ufreeDof = find(~isFixedDof);
    end
    if isempty(fixedDof) % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        fixedDof = 1;
        ufreeDof = 2:Nu;    % eliminate the kernel by enforcing u(1) = 0;
    end

    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedDof,fixedDof)=I, AD(fixedDof,ufreeDof)=0, AD(ufreeDof,fixedDof)=0.
    % BD(:,fixedDof) = 0 and thus BD'(fixedDof,:) = 0.
    bdidx = zeros(2*Nu,1);
    bdidx(fixedDof) = 1;
    bdidx(Nu+fixedDof) = 1;
    Tbd = spdiags(bdidx,0,2*Nu,2*Nu);
    T = spdiags(1-bdidx,0,2*Nu,2*Nu);
    AD = T*A*T + Tbd;
    BD = B*T;

    %% Part 2: Find boundary edges and modify the right hand side f and g
    % Find boundary edges: Neumann and Robin
    Neumann = []; Robin = []; %#ok<*NASGU>
    if ~isempty(bdFlag)
        isNeumann(elem2edge((bdFlag(:)==2)|(bdFlag(:) == 3))) = true;
        isRobin(elem2edge(bdFlag(:)==3)) = true;
        Neumannidx = find(isNeumann);
        Neumann   = edge(isNeumann,:);
        Robin     = edge(isRobin,:);
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        Neumann = edge(bdDof>N,:);
        if ~isempty(pde.g_R)
            Robin = Neumann;
        end
    end

    % Neumann boundary condition
    if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        [lambda,w] = quadpts1(3);
        nQuad = size(lambda,1);
        % quadratic bases (1---3---2)
        bdphi(:,1) = (2*lambda(:,1)-1).*lambda(:,1);
        bdphi(:,2) = (2*lambda(:,2)-1).*lambda(:,2);
        bdphi(:,3) = 4*lambda(:,1).*lambda(:,2);
        % length of edge
        ve = node(Neumann(:,1),:) - node(Neumann(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        % update RHS
        gex = zeros(size(Neumann,1),2);
        gey = zeros(size(Neumann,1),2);
        for pp = 1:nQuad
            pxy = lambda(pp,1)*node(Neumann(:,1),:)+lambda(pp,2)*node(Neumann(:,2),:);
            gp = pde.g_N(pxy);
            gex(:,1) = gex(:,1) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,1);
            gex(:,2) = gex(:,2) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,2);
            gey(:,1) = gey(:,1) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,1);
            gey(:,2) = gey(:,2) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,2);
            f1(N+Neumannidx) = f1(N+Neumannidx) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,3); % interior bubble
            f2(N+Neumannidx) = f2(N+Neumannidx) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,3); % interior bubble
        end
        f1(1:N) = f1(1:N) + accumarray(Neumann(:), gex(:),[N,1]);
        f2(1:N) = f2(1:N) + accumarray(Neumann(:), gey(:),[N,1]);
    end
    f = [f1; f2];
    % The case non-empty Neumann but g_N=[] corresponds to the zero flux
    % boundary condition on Neumann edges and no modification is needed.

    % Dirichlet boundary conditions
    if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
        u1 = zeros(Nu,1);
        u2 = zeros(Nu,1);
        idx = (fixedDof > N);              % index of edge dof
        uD = pde.g_D(node(fixedDof(~idx),:));  % bd value at vertex dofs
        u1(fixedDof(~idx)) = uD(:,1);
        u2(fixedDof(~idx)) = uD(:,2);
        bdEdgeIdx = fixedDof(idx)-N;
        bdEdgeMid = (node(edge(bdEdgeIdx,1),:)+node(edge(bdEdgeIdx,2),:))/2;
        uD = pde.g_D(bdEdgeMid);         % bd values at middle points of edges
        u1(fixedDof(idx)) = uD(:,1);
        u2(fixedDof(idx)) = uD(:,2);
        u = [u1; u2]; % Dirichlet bd condition is built into u
        f = f - A*u;  % bring affect of nonhomgenous Dirichlet bd condition to
        g = g - B*u;  % the right hand side
        g = g - mean(g);
        f(fixedDof) = u1(fixedDof);
        f(fixedDof+Nu) = u2(fixedDof);
    end
    % The case non-empty Dirichlet but g_D=[] corresponds to the zero Dirichlet
    % boundary condition and no modification is needed.

    % modfiy pressure dof for pure Dirichlet
    if isempty(Neumann)
        pDof = (1:Np-1)';
    end

    ufreeDof = [ufreeDof; Nu+ufreeDof];
    end % end of function getbdStokesP2P0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% TODO: impose compatible condition int(pde.g_D*n)=0 numerically.
% NOTE: the data should be compatible, which is not easy. One
% special case is the enclose flow i.e. pde.g_D*n=0 everywhere.

repo: ifem


text:
function [A,b,C] = update_spatial_components(Y,C,f,A_,P,options)

% update spatial footprints and background through Basis Pursuit Denoising
% for each pixel i solve the problem
%   [A(i,:),b(i)] = argmin sum(A(i,:))
%       subject to || Y(i,:) - A(i,:)*C + b(i)*f || <= sn(i)*sqrt(T);
% for each pixel the search is limited to a few spatial components

% INPUTS:
% Y:    raw data
% C:    temporal components
% f:    temporal background
% A_:   current estimate of spatial footprints (used for determining search locations only)
% P:    dataset parameters (used for noise values and interpolated entries)

% options    parameter struct (for noise values and other parameters)

% OUTPUTS:
% A:    new estimate of spatial footprints
% b:    new estimate of spatial background
% C:    temporal components (updated only when spatial components are completely removed)

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

warning('off', 'MATLAB:maxNumCompThreads:Deprecated');
memmaped = isobject(Y);
if memmaped
    sizY = size(Y,'Y');
    d = prod(sizY(1:end-1));
    T = sizY(end);
else
    [d,T] = size(Y);
end
if nargin < 6 || isempty(options); options = []; end
if ~isfield(options,'d1') || isempty(options.d1); d1 = input('What is the total number of rows? \n'); options.d1 = d1; else d1 = options.d1; end          % # of rows
if ~isfield(options,'d2') || isempty(options.d2); d2 = input('What is the total number of columns? \n'); options.d2 = d2; else d2 = options.d2; end          % # of columns
if ~isfield(options,'show_sum'); show_sum = 0; else show_sum = options.show_sum; end            % do some plotting while calculating footprints
if ~isfield(options,'interp'); Y_interp = sparse(d,T); else Y_interp = options.interp; end      % identify missing data
if ~isfield(options,'use_parallel'); use_parallel = ~isempty(which('parpool')); else use_parallel = options.use_parallel; end % use parallel toolbox if present
if ~isfield(options,'search_method'); method = []; else method = options.search_method; end     % search method for determining footprint of spatial components

if nargin < 2 || (isempty(A_) && isempty(C))  % at least either spatial or temporal components should be provided
    error('Not enough input arguments')
else
    if ~isempty(C); K = size(C,1); elseif islogical(A_); K = size(A_,2); else K = size(A_2,2) - options.nb; end
end

if nargin < 5 || isempty(P); P = preprocess_data(Y,1); end  % etsimate noise values if not present
if nargin < 4 || isempty(A_);
    IND = ones(d,size(C,1));
else
    if islogical(A_)     % check if search locations have been provided, otherwise estimate them
        IND = A_;
        if isempty(C)
            INDav = double(IND)/diag(sum(double(IND)));
            px = (sum(IND,2)>0);
            f = mean(Y(~px,:));
            b = max(Y*f',0)/norm(f)^2;
            C = max(INDav'*Y - (INDav'*b)*f,0);
        end
    else
        IND = determine_search_location(A_(:,1:K),method,options);
    end
end

options.sn = P.sn;
if ~memmaped
    Y(P.mis_entries) = NaN; % remove interpolated values
end

Cf = [C;f];

if use_parallel         % solve BPDN problem for each pixel
    Nthr = max(2*maxNumCompThreads,round(d*T/2^24));
    siz_row = [floor(d/Nthr)*ones(Nthr-1,1);d-floor(d/Nthr)*(Nthr-1)];
    indeces = [0;cumsum(siz_row)];
    if ~memmaped
        Ycell = mat2cell(Y,siz_row,T);
    else
        Ycell = cell(Nthr,1);
    end
    INDc =  mat2cell(IND,siz_row,K);
    Acell = cell(Nthr,1);
    Psnc = mat2cell(options.sn(:),siz_row,1);
    Yf = cell(Nthr,1);
    parfor nthr = 1:Nthr
        Acell{nthr} = zeros(siz_row(nthr),size(Cf,1));
        if memmaped
            Ytemp = double(Y.Yr(indeces(nthr)+1:indeces(nthr+1),:));
        else
            Ytemp = Ycell{nthr};
        end
        Yf{nthr} = Ytemp*f';
        for px = 1:siz_row(nthr)
            fn = ~isnan(Ytemp(px,:));       % identify missing data
            ind = find(INDc{nthr}(px,:));
            if ~isempty(ind);
                ind2 = [ind,K+(1:size(f,1))];
                %[~, ~, a, ~] = lars_regression_noise(Ycell{nthr}(px,fn)', Cf(ind2,fn)', 1, Psnc{nthr}(px)^2*T);
                [~, ~, a, ~] = lars_regression_noise(Ytemp(px,fn)', Cf(ind2,fn)', 1, Psnc{nthr}(px)^2*T);
                Acell{nthr}(px,ind2) = a';
            end
        end
    end
    A = cell2mat(Acell);
    Yf = cell2mat(Yf);
else
    A = [zeros(d,K),zeros(d,size(f,1))];
    sA = zeros(d1,d2);
    Yf = Y*f';
    for px = 1:d   % estimate spatial components
        fn = ~isnan(Y(px,:));       % identify missing data
        ind = find(IND(px,:));
        if ~isempty(ind);
            ind2 = [ind,K+(1:size(f,1))];
            [~, ~, a, ~] = lars_regression_noise(Y(px,fn)', Cf(ind2,fn)', 1, options.sn(px)^2*T);
            A(px,ind2) = a';
            sA(px) = sum(a);
        end
        if show_sum
            if mod(px,d1) == 0;
               figure(20); imagesc(sA); axis square;
               title(sprintf('Sum of spatial components (%i out of %i columns done)',round(px/d1),d2)); drawnow;
            end
        end
    end
end

A(isnan(A))=0;
A = sparse(A);
A = threshold_components(A,options);  % post-processing of components

fprintf('Updated spatial components \n');

ff = find(sum(A(:,1:K))==0);           % remove empty components
if ~isempty(ff)
    K = K - length(ff);
    A(:,ff) = [];
    C(ff,:) = [];
end

b = max((Yf - A(:,1:K)*(C(1:K,:)*f'))/(f*f'),0);
A = A(:,1:K);

author: zhoupc


text:
function value = p18_param ( action, name, value )

%*****************************************************************************80
%
%% P18_PARAM defaults, sets, or gets the parameters for problem 18.
%
%  Discussion:
%
%    In typical use, the user calls
%      p18_param ( 'DEFAULT' )
%    Then, once the variables have been given default values, any
%    value can be retrieved, as for example by:
%      value = p18_param ( 'GET', 'DELTA' )
%    or set to a new value, as in:
%      p18_param ( 'SET', 'DELTA', 7 )
%
%    DELTA is the only parameter.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION:
%    'DEFAULT' sets the variables to default values.
%    'SET' sets the variable chosen by NAME to have the value VALUE.
%    'GET' gets the variable chosen by NAME.
%
%    Input, string NAME, the name of the variable, which should
%    be 'DELTA'.  NAME is ignored if the ACTION is 'DEFAULT'.
%
%    Input, real VALUE, the value to be assigned the variable NAME,
%    if the value of ACTION was 'SET'.  Otherwise VALUE is ignored.
%
%    Output, real VALUE, depends on the input value of ACTION:
%    'DEFAULT', returns the number of variables;
%    'SET', returns the input value of VALUE;
%    'GET', returns the value of the variable NAME.
%
  persistent delta;

  if ( s_eqi ( action, 'DEFAULT' ) )

    delta = 0.5;

    value = 1;

  elseif ( s_eqi ( action, 'SET' ) )

    if ( s_eqi ( name, 'DELTA' ) )
      delta = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P18_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized variable name: "%s".\n', name );
      error ( 'P18_PARAM - Fatal error!' );
    end

  elseif ( s_eqi ( action, 'GET' ) )

    if ( s_eqi ( name, 'DELTA' ) )
      value = delta;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P18_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized variable name: "%s".\n', name );
      error ( 'P18_PARAM - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P18_PARAM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action: "%s".\n', action );
    error ( 'P18_PARAM - Fatal error!' );

  end

  return
end


text:
% vgg_vec  De-/vectorization of a matrix.
%
% For a matrix X, vgg_vec(X) = X(:).
% For a N^2-vector x, vgg_vec(x) = reshape(x,N,N).
%
% Classical matrix re-arrangement operator, see book Magnus-Neudecker.
% Trivial function, included mainly for consistency with notation in literature.
%
% Useful for rearranging matrix equations. Matrix from the middle of a product can be put to the right as
%
%    vgg_vec(A*B*C) = kron(C',A)*vgg_vec(B)
%
% See also vgg_vec_swap, vgg_commut_matrix, vgg_duplic_matrix.

function v = vgg_vec(A)

if any(size(A)==1)
  v = reshape(A,[1 1]*sqrt(prod(size(A))));
else
  v = A(:);
end

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function compute_movie_file(A, filename, options)

% compute_movie_file - create an avi or gif file
%
%   compute_movie_file(A, filename, options);
%
%   A can be a 3D array or a cell array of 2D matrices.
%   Each A{k} or A(:,:,k) is one frame in the animation.
%
%   Output avi or gif file depending on the extension of the filename.
%   (you can also use options.format='avi' or 'gif')
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
fps = getoptions(options, 'fps', 10);
format = getoptions(options, 'format', []);

if isempty(format)
    % find extension
    i = strfind(filename,'.');
    if isempty(i)
        error('Does not recognize extension, you should specify options.format.');
    end
    format = filename(i+1:end);
end

if not(iscell(A)) && not(isstr(A))
    B = A; A = {};
    for i=1:size(B,3)
        A{end+1} = B(:,:,i);
    end
    clear B;
end

% close all;
% figure;

if iscell(A)
    nb = length(A);
else
    nb = length( dir([A '*.png']) );
end

if strcmp(format, 'avi')
    % 'Indeo3', 'Indeo5', 'Cinepak', 'MSVC', 'RLE' or 'None'
    compressor = 'Cinepak';
    compressor = 'none';
    qual = 50;
    aviobj = avifile(filename, 'compression', compressor, 'quality', qual, 'fps', fps );
    % aviobj = avifile(filename, 'fps', fps );
end

imap = getoptions(options, 'colormap', jet(256));
loopCount = getoptions(options, 'loopcount', 1);

for j=1:nb
    progressbar(j,nb);
    if iscell(A)
        B = A{j};
    else
        % A is a filename
        B = load_image( [A '-' num2string_fixeddigit(j-1,3)] );
    end
    B = round(rescale( B, 1,256));
    switch format
        case 'avi'
            if size(B,3)==3
                frame = im2frame( rescale(B) );
            else
                frame = im2frame(B,gray(256));
            end
            % imageplot(A{j});
            % frame = getframe ( gca );
            aviobj = addframe( aviobj, frame );
        case 'gif'
            warning off;
            if j==1
                imwrite(B,imap,filename,'gif',...
                    'DelayTime',1/fps, 'LoopCount',loopCount,...
                    'WriteMode','overwrite')
            else
                imwrite(B,imap,filename,'gif','WriteMode','append',...
                    'DelayTime',1/fps);
            end
            warning on;
        otherwise
            error('Unknown format');
    end
end

if strcmp(format, 'avi')
    warning off;
    aviobj = close( aviobj );
    warning on;
end

author: gpeyre

repo: matlab-toolboxes


text:
function [Z]=rplus(X,Y)
%
% Replicating addition
%
% Does element by element operations on X and Y where non-same sized
% dimensions are implicity wrapped round to match the size of the larger
% to give a result matrix Z with size max(size(X),size(Y));
%
% In this case returns double array with X+Y
%
% See also repops, plus, minus
%
% Copyright 2006- by Jason D.R. Farquhar (jdrf@zepler.org)
% Inspired by code by Douglas M. Schwarz & Aki Vehtari.
Z=repop(X,'+',Y);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function K_mat=Kernel_func(x,y,mu,lamda)
xdatasize=size(x,1);
ydatasize=size(y,1);
for i=1:xdatasize
        for j=1:ydatasize
       K_mat(i,j)=mu*exp(-norm(x(i)-y(j))^2/(2*lamda^2));
        end
    end
end

author: epfl-lasa