text:
function [handles] = spm_uitab(hparent,labels,callbacks,tag,active,height,tab_height)
% Create tabs in a figure
% FORMAT [handles] = spm_uitab(hparent,labels,callbacks,tag,active,height,tab_height)
% This function creates tabs in the SPM graphics window.
% These tabs may be associated with different sets of axes and uicontrol,
% through the use of callback functions linked to the tabs.
% Inputs:
%   hparent    - the handle of the parent of the tabs (can be the SPM
%                graphics windows, or the handle of the uipanel of a former
%                spm_uitab...)
%   labels     - a cell array of string containing the labels of the tabs
%   callbacks  - a cell array of strings which will be evaluated using the
%                'eval' function when clicking on a tab [default: {[]}]
%   tag        - a string which is the tag associated with the tabs
%                (useful for finding them in a window...) [default: '']
%   active     - the index of the active tab when creating the uitabs
%                [default: 1, ie the first tab is active]
%   height     - the relative height of the tab panels within its parent
%                spatial extent [default: 1]
%   tab_height - the relative height of the tabs within its parent spatial
%                extent [default: 0.025]
% Output:
%   handles    - a structure of handles for the differents tab objects.
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau
% $Id: spm_uitab.m 6067 2014-06-26 15:33:30Z guillaume $

Ntabs = length(labels);

if ~exist('callbacks','var') || isempty(callbacks)
    for i=1:Ntabs
        callbacks{i} = [];
    end
end
if  ~exist('tag','var') || isempty(tag)
    tag = '';
end
if  ~exist('active','var') || isempty(active)
    active = 1;
end
if  ~exist('height','var') || isempty(height)
    height = 1;
end
if  ~exist('tab_height','var') || isempty(tab_height)
    tab_height = 0.025;
end
if ~isequal(get(hparent,'type'),'figure')
    set(hparent,'units','normalized')
    POS = get(hparent,'position');
    pos1 = [POS(1)+0.02,POS(2)+0.01,POS(3)-0.04,POS(4)-(tab_height+0.035)];
    dx = 0.095*(POS(3)-0.04)./0.98;
    dx2 = [0.04,0.93]*(POS(3)-0.04)./0.98;
else
    pos1 = [0.01 0.005 0.98 1-(tab_height+0.01)];
    dx = 0.095;
    dx2 = [0.04,0.93];
end
pos1(4) = pos1(4).*height;

COLOR = 0.95*[1 1 1];

handles.hp = uipanel(...
    'parent',hparent,...
    'position',pos1,...
    'BorderType','beveledout',...
    'BackgroundColor',COLOR,...
    'tag',tag);
set(handles.hp,'units','normalized');

xl = pos1(1);
yu = pos1(2) +pos1(4);
ddx = 0.0025;
ddy = 0.005;
dy = tab_height;

if Ntabs > 9
    handles.hs(1) = uicontrol(...
        'parent',hparent,...'enable','off',...
        'style','pushbutton',...
        'units','normalized','position',[xl yu dx2(1) dy],...
        'SelectionHighlight','off',...
        'BackgroundColor',COLOR,...
        'callback',@doScroll,...
        'value',0,'min',0,'max',Ntabs-9,...
        'string','<',...
        'tag',tag,...
        'BusyAction','cancel',...
        'Interruptible','off');

    handles.hs(2) = uicontrol(...
        'parent',hparent,...
        'style','pushbutton',...
        'units','normalized','position',[xl+dx2(2) yu 0.05 dy],...
        'SelectionHighlight','off',...
        'BackgroundColor',COLOR,...
        'callback',@doScroll,...
        'value',1,'min',1,'max',Ntabs-9,...
        'string','>',...
        'tag',tag,...
        'BusyAction','cancel',...
        'Interruptible','off');
    set(handles.hs,'units','normalized')
    xl = xl + dx2(1);
end

for i =1:min([Ntabs,9])
    pos = [xl+dx*(i-1) yu dx dy];
    handles.htab(i) = uicontrol(...
        'parent',hparent,...
        'style','pushbutton',...
        'units','normalized','position',pos,...
        'SelectionHighlight','off',...
        'string',labels{i},...
        'BackgroundColor',COLOR,...
        'tag',tag);
    set(handles.htab(i),'units','normalized')
    pos = [xl+dx*(i-1)+ddx yu-ddy dx-2*ddx 2*ddy];
    handles.hh(i) = uicontrol(...
        'parent',hparent,...
        'style','text',...
        'units','normalized','position',pos,...
        'BackgroundColor',COLOR,...
        'tag',tag);
    set(handles.hh(i),'units','normalized')
end
try
    set(handles.hh(active),'visible','on')
catch
    active = 1;
    set(handles.hh(active),'visible','on')
end
others = setdiff(1:min([Ntabs,9]),active);
set(handles.htab(active),...
    'FontWeight','bold');
set(handles.hh(others),'visible','off');
set(handles.htab(others),...
    'ForegroundColor',0.25*[1 1 1]);
ud.handles = handles;
ud.Ntabs = Ntabs;
for i =1:min([Ntabs,9])
    ud.ind = i;
    ud.callback = callbacks{i};
    set(handles.htab(i),'callback',@doChoose,'userdata',ud,...
        'BusyAction','cancel',...
        'Interruptible','off');
    if i > 9
        set(handles.htab(i),'visible','off');
    end
end

if Ntabs > 9
    UD.in = [1:9];
    UD.Ntabs = Ntabs;
    UD.h = handles;
    UD.active = active;
    UD.who = -1;
    UD.callbacks = callbacks;
    UD.labels = labels;
    set(handles.hs(1),'userdata',UD,'enable','off');
    UD.who = 1;
    set(handles.hs(2),'userdata',UD);
end

%==========================================================================
% doChoose
%==========================================================================
function doChoose(o1,o2)
ud = get(o1,'userdata');
% Do nothing if called tab is current (active) tab
if ~strcmp(get(ud.handles.htab(ud.ind),'FontWeight'),'bold')
    spm('pointer','watch');
    set(ud.handles.hh(ud.ind),'visible','on');
    set(ud.handles.htab(ud.ind),...
        'ForegroundColor',0*[1 1 1],...
        'FontWeight','bold');
    others = setdiff(1:length(ud.handles.hh),ud.ind);
    set(ud.handles.hh(others),'visible','off');
    set(ud.handles.htab(others),...
        'ForegroundColor',0.25*[1 1 1],...
        'FontWeight','normal');
    if ud.Ntabs >9
        UD = get(ud.handles.hs(1),'userdata');
        UD.active = UD.in(ud.ind);
        UD.who = -1;
        set(ud.handles.hs(1),'userdata',UD);
        UD.who = 1;
        set(ud.handles.hs(2),'userdata',UD);
    end
    drawnow
    if ~isempty(ud.callback)
        if isa(ud.callback, 'function_handle')
            feval(ud.callback);
        else
            eval(ud.callback);
        end
    end
    drawnow
    spm('pointer','arrow');
end

%==========================================================================
% doScroll
%==========================================================================
function doScroll(o1,o2)
ud = get(o1,'userdata');
% active = ud.in(ud.active);
ud.in = ud.in + ud.who;
if min(ud.in) ==1
    set(ud.h.hs(1),'enable','off');
    set(ud.h.hs(2),'enable','on');
elseif max(ud.in) ==ud.Ntabs
    set(ud.h.hs(1),'enable','on');
    set(ud.h.hs(2),'enable','off');
else
    set(ud.h.hs,'enable','on');
end
UD.handles = ud.h;
UD.Ntabs = ud.Ntabs;
for i = 1:length(ud.in)
    UD.ind = i;
    UD.callback = ud.callbacks{ud.in(i)};
    set(ud.h.htab(i),'userdata',UD,...
        'string',ud.labels{ud.in(i)});
    if ismember(ud.active,ud.in)
        ind = find(ud.in==ud.active);
        set(ud.h.hh(ind),'visible','on');
        set(ud.h.htab(ind),...
            'ForegroundColor',0*[1 1 1],...
            'FontWeight','bold');
        others = setdiff(1:9,ind);
        set(ud.h.hh(others),'visible','off');
        set(ud.h.htab(others),...
            'ForegroundColor',0.25*[1 1 1],...
            'FontWeight','normal');
    else
        others = 1:9;
        set(ud.h.hh(others),'visible','off');
        set(ud.h.htab(others),...
            'ForegroundColor',0.25*[1 1 1],...
            'FontWeight','normal');
    end
end
ud.who = -1;
set(ud.h.hs(1),'userdata',ud)
ud.who = 1;
set(ud.h.hs(2),'userdata',ud)


text:
function Y = spherharm(theta,phi,l,m);
% function Y = spherharm(theta,phi,l,m);

mp = abs(m);			% positive m

Y = sqrt((2*l + 1)/(4*pi) * (prod(1:(l-mp))/prod(1:(l+mp)))) * ...
    legendrex(cos(theta),mp,l) .* exp(sqrt(-1)*mp*phi);

if(m < 0),
  Y = (-1)^mp * conj(Y);
end
return

author: brainstorm-tools


text:
function [ind_Mp, ind_Md] = create_masks(ge, gridnd)
% Generating the indices to mask is better than generating a diagonal mask
% matrix, because multiplying the mask matrix cannot make NaN or Inf elements
% zeros.

chkarg(istypesizeof(ge, 'GT'), '"ge" should be instance of GT.');  % ge: grid type for the E-field
chkarg(istypesizeof(gridnd, 'Grid2d') || istypesizeof(gridnd, 'Grid3d'), '"gridnd" should be instance of Grid2d or Grid3d.');

v = Axis.x;
if istypesizeof(gridnd, 'Grid2d')
	v = Dir.h;
end

%% Get the shape.
N = gridnd.N;
bc = gridnd.bc;
ind = cell(1, v.count);

% Mask matrices
mask_p = cell(1, v.count);
for w = v.elems  % mask w-component of field on primary grid
	mask = false(N);
	us = setdiff(v.elems, w);

	for u = us
		if (ge == GT.prim && bc(u) == BC.e) || (ge == GT.dual && bc(u) == BC.m)
			[ind{:}] = deal(':');  % ind = {':', ':', ':'} for grid3d
			ind{u} = 1;
			mask(ind{:}) = true;
		end
	end
	mask_p{w} = mask(:);
end

ind_Mp = cell2mat(mask_p);
ind_Mp = ind_Mp(:);

mask_d = cell(1, v.count);
for w = v.elems  % mask w-component of field on dual grid
	mask = false(N);

	if (ge == GT.prim && bc(w) == BC.e) || (ge == GT.dual && bc(w) == BC.m)
		[ind{:}] = deal(':');  % ind = {':', ':', ':'} for grid3d
		ind{w} = 1;
		mask(ind{:}) = true;
	end
	mask_d{w} = mask(:);
end

ind_Md = cell2mat(mask_d);
ind_Md = ind_Md(:);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
% SOSDEMO4 --- Matrix Copositivity
% Section 3.4 of SOSTOOLS User's Manual
%

clear; echo on;
syms x1 x2 x3 x4 x5;
vartable = [x1; x2; x3; x4; x5];

% The matrix under consideration
J = [1 -1  1  1 -1;
    -1  1 -1  1  1;
     1 -1  1 -1  1;
     1  1 -1  1 -1;
    -1  1  1 -1  1];

% =============================================
% First, initialize the sum of squares program

prog = sosprogram(vartable);     % No decision variables.

% =============================================
% Next, define SOSP constraints

% Constraint : r(x)*J(x) - p(x) = 0
J = [x1^2 x2^2 x3^2 x4^2 x5^2]*J*[x1^2; x2^2; x3^2; x4^2; x5^2];
r = x1^2 + x2^2 + x3^2 + x4^2 + x5^2;

prog = sosineq(prog,r*J);

% =============================================
% And call solver
prog = sossolve(prog);

% =============================================
% Program is feasible. The matrix J is copositive.

echo off

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [C,S]=getEarth2014TerrainCoeffs(M,coeffType)
%%GETEARTH2014TERRAINCOEFFS Get fully normalized spherical harmonic
%                   coefficients for the radial distance offset of the
%                   terrain out from the surface of the GRS80 reference
%                   ellipsoid in meters under various parts of the degree
%                   2160 Earth2014 model.
%
%INPUTS: M The integer maximum order of the spherical harmonic
%          coefficients obtained. This is a value between 0 and 2160. If
%          this parameter is omitted, the default value is 2160.
% coeffType An optional parameter selecting the type of coefficients to
%          load. There are five variants of this model. Possible values
%          are
%          0 (The default if omitted or an empty matrix is passed) The
%            Earth's surface including water in lakes and major ice
%            sheets. oceans are placed on the GRS80 ellipsoid surface.
%          1 The Earth's bedrock.
%          2 The Earth's toporgraphy, bedrock plus major ice sheets.
%          3 The rock-equivalent topography of the Earth (ice and water
%            masses are condensed to layers of rock)
%          4 Major ice sheets. Everything else is set to the surface of the
%            GRS80 ellipsoid.
%
%OUTPUTS: C An array holding the coefficient terms that are multiplied by
%           cosines in the spherical harmonic expansion. If given to a
%           CountingClusterSet class, C(n+1,m+1) is the coefficient of
%           degree n and order m. When a maximum degree of M is used, all C
%           have values for all n from 0 to M and for all m from 0 to n for
%           each n. The coefficients are unitless.
%         S An array holding the coefficient terms that are multiplied by
%           sines in the spherical harmonic expansion. The format of S is
%           the same as that of C.
%
%The model is described in [1]. This function is only for values up to
%degree and order 2160 as the spherical harmonic synthesis routine in the
%Tracker Component Library cannot currently handle the full degree 10,800
%model, because it does not use extended precision arithmetic.
%
%The Earth 2014 data can be downloaded from
%http://ddfe.curtin.edu.au/models/Earth2014/
%To use the data for all of the models, place the files
%Earth2014.SUR2014.degree2160.bshc, Earth2014.BED2014.degree2160.bshc
%Earth2014.TBI2014.degree2160.bshc, Earth2014.RET2014.degree2160.bshc
%and Earth2014.ICE2014.degree2160.bshc in the data folder that is in the
%same folder as this function.
%
%EXAMPLE:
%The fact that this only stores radial distance offsets can be a little
%confusing when one wants to find the location of a point on the Earth's
%surface in the WGS-84 coordinate system/ International Terrestrial
%Reference System (ITRS). Here, we give an example:
% %WGS-84 latitude and longitude in radians.
% latLon=[19.4721;-155.5922]*(pi/180);
% [C,S]=getEarth2014TerrainCoeffs();
% %To use the terrain coefficients, the latitude must be changed from
% %geodetic to geocentric. We will assume that the latitude is given in
% %WGS-84 coordinates, as is standard with GPS.
% latSpher=ellipsLat2SpherLat(latLon(1));
% azElSphere=[latLon(2);latSpher];
% terHeight=spherHarmonicEval(C,S,azElSphere);
% %Add the ellipsoidal radius of the GRS80 reference ellipsoid at the
% %point. Note that this actually changes the ellipsoidal latitude
% %slightly.
% a=Constants.GRS80SemiMajorAxis;
% f=Constants.GRS80Flattening;
% terRad=terHeight+ellipsoidalRadius(latSpher,0,a,f);
% %Convert from spherical to Cartesian coordinates to find the point.
% CartLoc=spher2Cart([terRad;azElSphere])
%
%REFERENCES:
%[1] C. Hirt and M. Rexer, "Earth2014: 1arc-min shape, topography,
%    bedrock and ice-sheet models - available as gridded data and degree-
%    10,800 spherical harmonics," International Journal of Applied Earth
%    Observation and Geoinformation, vol. 39, pp. 103-112, Jul. 2015.
%
%June 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(coeffType))
     coeffType=0;
end

%The Earth2014 terrain coefficient data file, should be located in a data
%folder that is in the same folder as this file. This finds the path to
%this file.
ScriptPath=mfilename('fullpath');
ScriptFolder = fileparts(ScriptPath);

switch(coeffType)
    case 0
        fileID=fopen([ScriptFolder,'/data/Earth2014.SUR2014.degree2160.bshc'],'rb');
    case 1
        fileID=fopen([ScriptFolder,'/data/Earth2014.BED2014.degree2160.bshc'],'rb');
    case 2
        fileID=fopen([ScriptFolder,'/data/Earth2014.TBI2014.degree2160.bshc'],'rb');
    case 3
        fileID=fopen([ScriptFolder,'/data/Earth2014.RET2014.degree2160.bshc'],'rb');
    case 4
        fileID=fopen([ScriptFolder,'/data/Earth2014.ICE2014.degree2160.bshc'],'rb');
    otherwise
        error('Unknown coefficient type specified.')
end

data = fread(fileID,Inf,'double');
fclose(fileID);

%Note that data(1) should be zero. data(1) is the degree of the lowest
%degree coefficient. data(2) is the degree of the maximum degree
%coefficient.
maxDeg = data(2);
maxNumCoeffs=(maxDeg+1)*(maxDeg+2)/2;

if(nargin<1||isempty(M))
    M=maxDeg;
end

totalNumCoeffs=(M+1)*(M+2)/2;

%Extract the data. The additive 2 skips the first two entries, which
%indicate the lowest and highest coefficient degrees.
C=data((2+1):(2+totalNumCoeffs));
S=data((2+maxNumCoeffs+1):(2+maxNumCoeffs+totalNumCoeffs));
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Test Jacobian primitives

p560

jacob0(p560, qz)
jacob0(p560, qr)
jacob0(p560, qn)

jacobn(p560, qz)
jacobn(p560, qr)
jacobn(p560, qn)

t = fkine(p560, qn)
tr2jac(t)

author: petercorke

repo: robotics-toolbox-matlab


text:
 function diff_im = anisodiff2D(im, num_iter, delta_t, kappa, option)
%ANISODIFF2D Conventional anisotropic diffusion
%   DIFF_IM = ANISODIFF2D(IM, NUM_ITER, DELTA_T, KAPPA, OPTION) perfoms
%   conventional anisotropic diffusion (Perona & Malik) upon a gray scale
%   image. A 2D network structure of 8 neighboring nodes is considered for
%   diffusion conduction.
%
%       ARGUMENT DESCRIPTION:
%               IM       - gray scale image (MxN).
%               NUM_ITER - number of iterations.
%               DELTA_T  - integration constant (0 <= delta_t <= 1/7).
%                          Usually, due to numerical stability this
%                          parameter is set to its maximum value.
%               KAPPA    - gradient modulus threshold that controls the conduction.
%               OPTION   - conduction coefficient functions proposed by Perona & Malik:
%                          1 - c(x,y,t) = exp(-(nablaI/kappa).^2),
%                              privileges high-contrast edges over low-contrast ones.
%                          2 - c(x,y,t) = 1./(1 + (nablaI/kappa).^2),
%                              privileges wide regions over smaller ones.
%
%       OUTPUT DESCRIPTION:
%                DIFF_IM - (diffused) image with the largest scale-space parameter.
%
%   Example
%   -------------
%   s = phantom(512) + randn(512);
%   num_iter = 15;
%   delta_t = 1/7;
%   kappa = 30;
%   option = 2;
%   ad = anisodiff2D(s,num_iter,delta_t,kappa,option);
%   figure, subplot 121, imshow(s,[]), subplot 122, imshow(ad,[])
%
% See also anisodiff1D, anisodiff3D.

% References:
%   P. Perona and J. Malik.
%   Scale-Space and Edge Detection Using Anisotropic Diffusion.
%   IEEE Transactions on Pattern Analysis and Machine Intelligence,
%   12(7):629-639, July 1990.
%
%   G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz.
%   Nonlinear Anisotropic Filtering of MRI Data.
%   IEEE Transactions on Medical Imaging,
%   11(2):221-232, June 1992.
%
%   MATLAB implementation based on Peter Kovesi's anisodiff(.):
%   P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing.
%   School of Computer Science & Software Engineering,
%   The University of Western Australia. Available from:
%   <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
%
% Credits:
% Daniel Simoes Lopes
% ICIST
% Instituto Superior Tecnico - Universidade Tecnica de Lisboa
% danlopes (at) civil ist utl pt
% http://www.civil.ist.utl.pt/~danlopes
%
% May 2007 original version.

% Convert input image to double.
im = double(im);

% PDE (partial differential equation) initial condition.
diff_im = im;

% Center pixel distances.
dx = 1;
dy = 1;
dd = sqrt(2);

% 2D convolution masks - finite differences.
hN = [0 1 0; 0 -1 0; 0 0 0];
hS = [0 0 0; 0 -1 0; 0 1 0];
hE = [0 0 0; 0 -1 1; 0 0 0];
hW = [0 0 0; 1 -1 0; 0 0 0];
hNE = [0 0 1; 0 -1 0; 0 0 0];
hSE = [0 0 0; 0 -1 0; 0 0 1];
hSW = [0 0 0; 0 -1 0; 1 0 0];
hNW = [1 0 0; 0 -1 0; 0 0 0];

% Anisotropic diffusion.
for t = 1:num_iter

        % Finite differences. [imfilter(.,.,'conv') can be replaced by conv2(.,.,'same')]
        nablaN = imfilter(diff_im,hN,'conv');
        nablaS = imfilter(diff_im,hS,'conv');
        nablaW = imfilter(diff_im,hW,'conv');
        nablaE = imfilter(diff_im,hE,'conv');
        nablaNE = imfilter(diff_im,hNE,'conv');
        nablaSE = imfilter(diff_im,hSE,'conv');
        nablaSW = imfilter(diff_im,hSW,'conv');
        nablaNW = imfilter(diff_im,hNW,'conv');

        % Diffusion function.
        if option == 1
            cN = exp(-(nablaN/kappa).^2);
            cS = exp(-(nablaS/kappa).^2);
            cW = exp(-(nablaW/kappa).^2);
            cE = exp(-(nablaE/kappa).^2);
            cNE = exp(-(nablaNE/kappa).^2);
            cSE = exp(-(nablaSE/kappa).^2);
            cSW = exp(-(nablaSW/kappa).^2);
            cNW = exp(-(nablaNW/kappa).^2);
        elseif option == 2
            cN = 1./(1 + (nablaN/kappa).^2);
            cS = 1./(1 + (nablaS/kappa).^2);
            cW = 1./(1 + (nablaW/kappa).^2);
            cE = 1./(1 + (nablaE/kappa).^2);
            cNE = 1./(1 + (nablaNE/kappa).^2);
            cSE = 1./(1 + (nablaSE/kappa).^2);
            cSW = 1./(1 + (nablaSW/kappa).^2);
            cNW = 1./(1 + (nablaNW/kappa).^2);
        end

        % Discrete PDE solution.
        diff_im = diff_im + ...
                  delta_t*(...
                  (1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...
                  (1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE + ...
                  (1/(dd^2))*cNE.*nablaNE + (1/(dd^2))*cSE.*nablaSE + ...
                  (1/(dd^2))*cSW.*nablaSW + (1/(dd^2))*cNW.*nablaNW );

        % Iteration warning.
        %fprintf('\rIteration %d\n',t);
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function lobatto_set_test ( )

%*****************************************************************************80
%
%% LOBATTO_SET_TEST tests LOBATTO_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LOBATTO_SET_TEST\n' );
  fprintf ( 1, '  LOBATTO_SET sets a Lobatto rule;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         I      X             W\n' );

  for n = 4 : 3 : 12

    [ x, w ] = lobatto_set ( n );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %8d  %12f  %12f\n', i, x(i), w(i) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef prtClassMap < prtClass
 %prtClassMap  Maximum a Posteriori classifier
 %
 %    CLASSIFIER = prtClassMap returns a Maximum a Posteriori classifier
 %
 %    CLASSIFIER = prtClassMap(PROPERTY1, VALUE1, ...) constructs a
 %    prtClassMAP object CLASSIFIER with properties as specified by
 %    PROPERTY/VALUE pairs.
 %
 %    A prtClassMap object inherits all properties from the abstract class
 %    prtClass. In addition is has the following property:
 %
 %    rvs    - A prtRv object. This property describes the random variable
 %             model used for Maximum a Posteriori classification.
 %
 %    A prtClassMap object inherits inherits the TRAIN, RUN, CROSSVALIDATE
 %    and KFOLDS methods from prtClass.
 %
 %    Example:
 %
 %    Test` = prtDataGenUnimodal;       % Create some test and
 %    TrainingDataSet = prtDataGenUnimodal;   % training data
 %    classifier = prtClassMap;               % Create a classifier
 %    classifier = classifier.train(TrainingDataSet);    % Train
 %    classified = run(classifier, TestDataSet);         % Test
 %
 %    subplot(2,1,1); classifier.plot;  % Plot results
 %    subplot(2,1,2); prtScoreRoc(classified,TestDataSet);
 %    set(get(gca,'Children'), 'LineWidth',3)
 %
 %    See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
 %    prtClassMap, prtClassFld, prtClassBinaryToMaryOneVsAll, prtClassDlrt,
 %    prtClassPlsda, prtClassFld, prtClassRvm, prtClassGlrt,  prtClassSvm,
 %    prtClassTreeBaggingCap, prtClassKmsd, prtClassKnn

    properties (SetAccess=private)
        % Required by prtAction
        name = 'Maximum a Posteriori'   % Maximum a Posteriori
        nameAbbreviation = 'MAP'        % MAP
        isNativeMary = true;            % True
    end

    properties
        rvs = prtRvMvn; % Random variable object containing mean and variance
    end
    properties (Hidden)
        runLogLikelihoods = false
    end

    methods
        % Constructor
        function self = prtClassMap(varargin)

            self.classTrain = 'prtDataInterfaceCategoricalTargets';
            self.classRun = 'prtDataSetBase';
            self.classRunRetained = false;

            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
        % Set function
        function self = set.rvs(self,val)
            if ~(isa(val, 'prtRv') || isa(val, 'prtBrv'))
                error('prtClassMAP:rvs','Rvs parameter must be of class prtRv');
            else
                self.rvs = val;
            end
        end
    end

    methods (Access = protected, Hidden = true)

        function self = trainAction(self,ds)

            % Repmat the rv objects to get one for each class
            % The pattern of supplied RVs repeats, if 2 rvs are suplied for
            % a three class problem, the third class gets the same as class
            % 1, if there were a fourth class if would be the same as class
            % 2.
            self.rvs = repmat(self.rvs(:), (ds.nClasses - length(self.rvs)+1),1);
            self.rvs = self.rvs(1:ds.nClasses);

            % Get the ML estimates of the RV parameters for each class
            for iY = 1:ds.nClasses
                self.rvs(iY) = train(self.rvs(iY), ds.retainClassesByInd(iY));
            end
        end

        function ds = runAction(self,ds)

            % We call run for each RV
            % Typically this is the loglikelihood, but it might not be.
            logLikelihoods = zeros(ds.nObservations, length(self.rvs));
            for iY = 1:length(self.rvs)
                logLikelihoods(:,iY) = getObservations(run(self.rvs(iY), ds));
            end

            if ~self.runLogLikelihoods
                % If we don't want loglikelihoods transform them to
                % probabilities
                logLikelihoods = exp(bsxfun(@minus, logLikelihoods, prtUtilSumExp(logLikelihoods.').'));
            end

            if ~isa(ds,'prtDataSetClass')
                ds = ds.toPrtDataSetClassNoData();
            end
            ds.X = logLikelihoods;
        end

        function xOut = runActionFast(self,xIn,ds) %#ok<INUSD>

            logLikelihoods = zeros(size(xIn,1),length(self.rvs));
            for iY = 1:length(self.rvs)
                logLikelihoods(:,iY) = runFast(self.rvs(iY), xIn);
            end

            if ~self.runLogLikelihoods
                % If we don't want loglikelihoods transform them to
                % probabilities
                xOut = exp(bsxfun(@minus, logLikelihoods, prtUtilSumExp(logLikelihoods.').'));
            else
                xOut = logLikelihoods;
            end
        end
    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% function output = test_add_overlay(input)
% Test image overlay creation and plotting of 3D-RGB-image via montage
%
%   output = test_add_overlay(input)
%
% IN
%
% OUT
%
% EXAMPLE
%   test_add_overlay
%
%   See also tapas_uniqc_add_overlay

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-25
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

clear;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load example data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pathExamples    = tapas_uniqc_get_path('examples');
fileTest        = fullfile(pathExamples, 'nifti', 'rest', 'meanfmri.nii');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create Example overlay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X               = MrImage(fileTest);
edgeX           = edge(X);

imageMatrix     = X.data;
overlayMatrix   = edgeX.data;

rgbMatrix       = tapas_uniqc_add_overlay(imageMatrix, overlayMatrix, 'jet');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot resulting overlay via montage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X.plot('overlayImages', edgeX)

stringTitle = 'imageMatrix';
figure('Name', stringTitle);
montage(permute(imageMatrix, [1 2 4 3]), 'DisplayRange', []);
title(tapas_uniqc_str2label(stringTitle));

stringTitle = 'overlayMatrix';
figure('Name', stringTitle);
montage(permute(overlayMatrix, [1 2 4 3]), 'DisplayRange', []);
title(tapas_uniqc_str2label(stringTitle));
colormap hot;

stringTitle = 'test_add_overlay Montage Image';
figure('Name', stringTitle);
montage(rgbMatrix, 'DisplayRange', []);
title(tapas_uniqc_str2label(stringTitle));

author: translationalneuromodeling

repo: tapas


text:
% std_dipplot() - Commandline function to plot cluster component dipoles. Dipoles for each
%                 named cluster is displayed in a separate figure. To view all the clustered
%                 components in the STUDY on the same figure (in a separate subplot), all
%                 STUDY clusters must be requested.
%                 To visualize dipoles, they first must be stored in the EEG dataset structures
%                 using dipfit(). Only components that have dipole locations will be displayed,
%                 along with the cluster mean dipole (in red).
% Usage:
%                 >> [STUDY] = std_dipplot(STUDY, ALLEEG, clusters);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the dataset(s) included in
%                the STUDY. ALLEEG for a STUDY set is typically created using load_ALLEEG().
%
% Optional inputs:
%   'clusters' - [numeric vector | 'all']  -> specific cluster numbers to plot.
%                'all'  -> plot all clusters in STUDY.
%                {default: 'all'}.
%   'comps'    - [numeric vector]  -> indices of the cluster components to plot.
%                'all'  -> plot all the components in the cluster
%                {default: 'all'}.
%   'mode'     - ['together'|'apart'] Display all requested cluster on one
%                figure ('together') or separate figures ('apart').
%                'together'-> plot all 'clusters' in one figure (without the gui).
%                'apart'   -> plot each cluster in a separate figure. Note that
%                this parameter has no effect if the 'comps' option (above) is used.
%                {default: 'together'}
%   'figure'   - ['on'|'off'] plots on a new figure ('on')  or plots on current
%                figure ('off'). If 'figure','off' does not display gui controls,
%                Useful for incomporating a cluster dipplot into a complex figure.
%                {default: 'on'}.
% Outputs:
%   STUDY      - the input STUDY set structure modified with plotted cluster
%                mean dipole, to allow quick replotting (unless cluster means
%                already exists in the STUDY).
%   Example:
%   >> [STUDY] = std_dipplot(STUDY,ALLEEG, 'clusters', 5, 'mode', 'apart', 'figure', 'off');
%                % Plot cluster-5 component dipoles (in blue), plus ther mean dipole (in red),
%                % on an exisiting (gui-less) figure.
%
%  See also  pop_clustedit(), dipplot()
%
% Authors:  Hilit Serby, Arnaud Delorme, Scott Makeig, SCCN, INC, UCSD, June, 2005

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, June 08, 2005, hilit@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function STUDY = std_dipplot(STUDY, ALLEEG, varargin)

% Set default values
cls = []; % plot all clusters in STUDY
figureon = 1; % plot on a new figure
mode = 'apart';
opt_dipplot = {'projlines','off', 'normlen', 'on', 'pointout', 'on', 'verbose', 'off', 'dipolelength', 0,'spheres','on'};
%, 'spheres', 'on'
for k = 3:2:nargin
    switch varargin{k-2}
        case 'clusters'
            if isnumeric(varargin{k-1})
                cls = varargin{k-1};
                if isempty(cls)
                    cls = 2:length(STUDY.cluster);
                end
            else
                if isstr(varargin{k-1}) & strcmpi(varargin{k-1}, 'all')
                    cls = 2:length(STUDY.cluster);
                else
                    error('std_dipplot: ''clusters'' input takes either specific clusters (numeric vector) or keyword ''all''.');
                end
            end
            if length(cls) == 1, mode = 'apart'; else mode = 'together'; end;
        case 'comps'
            STUDY = std_plotcompdip(STUDY, ALLEEG,  cls, varargin{k-1});
            return;
        case 'plotsubjects', % do nothing
        case 'mode', mode = varargin{k-1};
        case 'figure'
            if strcmpi(varargin{k-1},'off')
                opt_dipplot{end + 1} = 'gui';
                opt_dipplot{end + 1} = 'off';
                figureon = 0;
            else
                opt_dipplot{end + 1} = 'projimg';
                opt_dipplot{end + 1} = 'on';
            end
    end
end

% select clusters to plot
% -----------------------
if isempty(cls)
    tmp =[];
    cls = 2:length(STUDY.cluster); % plot all clusters in STUDY
    for k = 1: length(cls)
        % don't include 'Notclust' clusters
        if ~strncmpi('Notclust',STUDY.cluster(cls(k)).name,8) & ~strncmpi('ParentCluster',STUDY.cluster(cls(k)).name,13)
            tmp = [tmp cls(k)];
        end
    end
    cls = tmp;
end;

if strcmpi(mode, 'apart')  % case each cluster on a separate figure
    for clus = 1: length(cls) % For each cluster requested
        if length(STUDY.cluster(cls(clus)).comps) > 0  % check there are comps in cluster
            max_r = 0;
            clear cluster_dip_models;
            len = length(STUDY.cluster(cls(clus)).comps);
            ndip = 0;
            dip_ind = [];
            if ~isfield(STUDY.cluster(cls(clus)),'dipole')
                STUDY = std_centroid(STUDY,ALLEEG, cls(clus) , 'dipole');
            elseif isempty(STUDY.cluster(cls(clus)).dipole)
                STUDY = std_centroid(STUDY,ALLEEG, cls(clus) , 'dipole');
            end
            for k = 1:len
                abset   = STUDY.datasetinfo(STUDY.cluster(cls(clus)).sets(1,k)).index;
                subject = STUDY.datasetinfo(STUDY.cluster(cls(clus)).sets(1,k)).subject;
               if ~isfield(ALLEEG(abset), 'dipfit')
                   warndlg2(['No dipole information available in dataset ' ALLEEG(abset).filename ' , abort plotting'], 'Aborting plot dipoles');
                   return;
               end
               comp = STUDY.cluster(cls(clus)).comps(k);
               cluster_dip_models(k).posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
               cluster_dip_models(k).momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
               cluster_dip_models(k).rv = ALLEEG(abset).dipfit.model(comp).rv;
               if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                   if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                       load('-mat', ALLEEG(abset).dipfit.hdmfile);
                       max_r = max(max_r, max(vol.r));
                   else % old version of dipfit
                       max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
                   end
               end
               comp_to_disp{k} = [subject  ', ' 'IC' num2str(comp) ];
               if ~isempty(cluster_dip_models(k).posxyz)
                   ndip = ndip +1;
                   dip_ind = [dip_ind k];
               end
            end % finished going over cluster comps
            STUDY.cluster(cls(clus)).dipole = computecentroid(cluster_dip_models);
            cluster_dip_models(end + 1) = STUDY.cluster(cls(clus)).dipole;

           % additional options
           % ------------------
           dip_color = cell(1,ndip+1);
           dip_color(1:ndip) = {'b'};
           dip_color(end) = {'r'};
           options = opt_dipplot;
           options{end+1} =  'mri';
           options{end+1} =  ALLEEG(abset).dipfit.mrifile;
           options{end+1} =  'coordformat';
           options{end+1} =  ALLEEG(abset).dipfit.coordformat;
           options{end+1} =  'dipnames';
           options{end+1} = {comp_to_disp{dip_ind } [STUDY.cluster(cls(clus)).name ' mean']};
           options{end+1} = 'color';
           options{end+1} = dip_color;
           if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
               options{end+1} = 'sphere';
               options{end+1} = max_r;
           else
               options{end+1} = 'meshdata';
               options{end+1} = ALLEEG(abset).dipfit.hdmfile;
           end
           if ndip < 6 && strcmpi(options{1}, 'projlines') && length(cls) == 1 % less than 6 dipoles, project lines
               options{2} = 'on';
           end

           if figureon
               dipplot(cluster_dip_models, options{:});
               fig_h = gcf;
               set(fig_h,'Name', [STUDY.cluster(cls(clus)).name ' - ' num2str(length(unique(STUDY.cluster(cls(clus)).sets(1,:)))) ...
                       ' sets - ' num2str(length(STUDY.cluster(cls(clus)).comps)) ' components (' num2str(ndip) ' dipoles)' ],'NumberTitle','off');
           else
               dipplot(cluster_dip_models, options{:},'view', [0.5 -0.5 0.5],'spheres', 'on', 'projlines', 'on');
               for gind = 1:length(options) % remove the 'gui' 'off' option
                   if isstr(options{gind})
                       if strfind(options{gind}, 'gui')
                           break;
                       end
                   end
               end
               options(gind:gind+1) = [];
               dipinfo.dipmod =  cluster_dip_models;
               dipinfo.op = options;
               diptitle = [STUDY.cluster(cls(clus)).name ', ' num2str(length(unique(STUDY.cluster(cls(clus)).sets(1,:)))) ' sets -' ...
                   num2str(length(STUDY.cluster(cls(clus)).comps)) ' components (' num2str(ndip) ' dipoles)' ];
               dipinfo.title = diptitle;
               set(gcf, 'UserData', dipinfo);
               set(gca,'UserData', dipinfo);
               rotate3d off;
               axcopy(gca, ['dipinfo = get(gca, ''''UserData''''); dipplot(dipinfo.dipmod, dipinfo.op{:}); set(gcf, ''''Name'''', dipinfo.title,''''NumberTitle'''',''''off''''); ']);
           end
        end % finished the if condition that cluster isn't empty
    end % finished going over requested clusters
end

if strcmpi(mode, 'together')  % case all clusters are plotted in the same figure (must be a new figure)
    N = length(cls);
    rowcols(2) = ceil(sqrt(N)); % Number of rows in the subplot figure.
    rowcols(1) = ceil(N/rowcols(2));
    fig_h = figure;
    orient tall
    set(fig_h,'Color', 'black');
    set(fig_h,'Name', 'All clusters dipoles','NumberTitle','off');
    set(fig_h, 'resize','off');
    for l = 1:N
        len = length(STUDY.cluster(cls(l)).comps);
        max_r = 0;
        clear cluster_dip_models;
        if ~isfield(STUDY.cluster(cls(l)),'dipole')
            STUDY = std_centroid(STUDY,ALLEEG, cls(l), 'dipole');
        elseif isempty(STUDY.cluster(cls(l)).dipole)
            STUDY = std_centroid(STUDY,ALLEEG, cls(l), 'dipole');
        end
        for k = 1: len
            abset = STUDY.datasetinfo(STUDY.cluster(cls(l)).sets(1,k)).index;
           if ~isfield(ALLEEG(abset), 'dipfit')
               warndlg2(['No dipole information available in dataset ' num2str(abset) ' , abort plotting'], 'Aborting plot dipoles');
               return;
           end
           comp = STUDY.cluster(cls(l)).comps(k);
           cluster_dip_models(k).posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
           cluster_dip_models(k).momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
           cluster_dip_models(k).rv = ALLEEG(abset).dipfit.model(comp).rv;
           if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                   load('-mat', ALLEEG(abset).dipfit.hdmfile);
                   max_r = max(max_r, max(vol.r));
               else % old version of dipfit
                   max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
               end
           end
        end % finished going over cluster comps
        STUDY.cluster(cls(l)).dipole = computecentroid(cluster_dip_models);
        cluster_dip_models(end + 1) = STUDY.cluster(cls(l)).dipole;
        dip_color = cell(1,length(cluster_dip_models));
        dip_color(1:end-1) = {'b'};
        dip_color(end) = {'r'};
        options = opt_dipplot;
        options{end + 1} =  'gui';
        options{end + 1} =  'off';
        options{end+1} =  'mri';
        options{end+1} =  ALLEEG(abset).dipfit.mrifile;
        options{end+1} =  'coordformat';
        options{end+1} =  ALLEEG(abset).dipfit.coordformat;
        options{end+1} = 'color';
        options{end+1} = dip_color;
        if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
            options{end+1} = 'sphere';
            options{end+1} = max_r;
        else
            options{end+1} = 'meshdata';
            options{end+1} = ALLEEG(abset).dipfit.hdmfile;
        end
        subplot(rowcols(1),rowcols(2),l) ,
        dipplot(cluster_dip_models, options{:});
        title([ STUDY.cluster(cls(l)).name ' (' num2str(length(unique(STUDY.cluster(cls(l)).sets(1,:)))) ' Ss, '  num2str(length(STUDY.cluster(cls(l)).comps)),' ICs)'],'color','white');
        %diptitle = [STUDY.cluster(cls(l)).name ', ' num2str(length(unique(STUDY.cluster(cls(l)).sets(1,:)))) 'Ss'];
        %title(diptitle, 'Color', 'white');
        % Complex axcopy
        %if l == 1
        %    for gind = 1:length(options) % remove the 'gui' 'off' option
        %        if isstr(options{gind})
        %            if strfind(options{gind}, 'gui')
        %                break;
        %            end
        %        end
        %    end
        %    options(gind:gind+1) = [];
        %end
        %dipinfo.dipmod =  cluster_dip_models;
        %dipinfo.op = options;
        %dipinfo.title = diptitle;
        %set(gcf, 'UserData', dipinfo);
        %set(gca,'UserData', dipinfo);
        %axcopy(gcf, ['dipinfo = get(gca, ''''UserData''''); dipplot(dipinfo.dipmod, dipinfo.op{:}); set(gcf, ''''Name'''', dipinfo.title,''''NumberTitle'''',''''off'''');']);
   end %finished going over all clusters
   set(fig_h, 'resize','on');
end % finished case of 'all' clusters
% std_plotcompdip() - Commandline function, to visualizing cluster components dipoles.
%                   Displays the dipoles of specified cluster components with the cluster mean
%                   dipole on separate figures.
%                   To visualize dipoles they first must be stored in the EEG dataset structures
%                   using dipfit(). Only components that have a dipole locations will be displayed,
%                   along with the cluster mean dipole in red.
% Usage:
%                   >> [STUDY] = std_plotcompdip(STUDY, ALLEEG, cluster, comps);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the dataset(s) included in the STUDY.
%                     ALLEEG for a STUDY set is typically created using load_ALLEEG().
%   cluster     - single cluster number.
%
% Optional inputs:
%   comps      - [numeric vector]  -> indices of the cluster components to plot.
%                       'all'                       -> plot all the components in the cluster {default: 'all'}.
%
% Outputs:
%   STUDY    - the input STUDY set structure modified with plotted cluster
%                     dipole mean, to allow quick replotting (unless cluster mean
%                     already existed in the STUDY).
%
%   Example:
%                         >> cluster = 4; comps= 1;
%                         >> [STUDY] = std_plotcompdip(STUDY,ALLEEG, cluster, comps);
%                    Plots component 1 dipole in blue with the cluster 4 mean dipole in red.
%
%  See also  pop_clustedit, dipfit, std_dipplot
%
% Authors:  Hilit Serby, Arnaud Delorme, Scott Makeig, SCCN, INC, UCSD, June, 2005

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, June 08, 2005, hilit@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function STUDY = std_plotcompdip(STUDY, ALLEEG, cls, varargin)
if ~exist('cls')
    error('std_plotcompdip: you must provide a cluster number as an input.');
end
if isempty(cls)
   error('std_plotcompdip: you must provide a cluster number as an input.');
end
if nargin == 3 % no components indices were given
    % Default plot all components of the cluster
    [STUDY] = std_dipplot(STUDY, ALLEEG, 'clusters', cls);
    return
else
    comp_ind = varargin{1};
end

for ci = 1:length(comp_ind)
    abset = STUDY.datasetinfo(STUDY.cluster(cls).sets(1,comp_ind(ci))).index;
    comp = STUDY.cluster(cls).comps(comp_ind(ci));
    subject = STUDY.datasetinfo(STUDY.cluster(cls).sets(1,comp_ind(ci))).subject;
    if ~isfield(ALLEEG(abset), 'dipfit')
        warndlg2(['No dipole information available in dataset ' num2str(abset) ' , abort plotting'], 'Aborting plot dipoles');
        return;
    end
    if length(comp_ind) == 1 & isempty(ALLEEG(abset).dipfit.model(comp).posxyz)
        warndlg2(strvcat('There is no dipole information available in', ...
                       [ 'dataset ' num2str(abset) ' for this component, abort plotting']), 'Aborting plot dipoles');
        return;
    end;
    if ~isfield(STUDY.cluster(cls),'dipole')
        STUDY = std_centroid(STUDY,ALLEEG, cls , 'dipole');
    elseif isempty(STUDY.cluster(cls).dipole)
        STUDY = std_centroid(STUDY,ALLEEG, cls , 'dipole');
    end
    comp_to_disp = [subject  ' / ' 'IC' num2str(comp) ];
    cluster_dip_models.posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
    cluster_dip_models.momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
    cluster_dip_models.rv = ALLEEG(abset).dipfit.model(comp).rv;
    cluster_dip_models(2) = STUDY.cluster(cls).dipole;
    if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
        if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
            load('-mat', ALLEEG(abset).dipfit.hdmfile);
            max_r = max(vol.r);
        else
            max_r = max(ALLEEG(abset).dipfit.vol.r);
        end
        dipplot(cluster_dip_models, 'sphere', max_r, 'mri', ALLEEG(abset).dipfit.mrifile,'coordformat', ALLEEG(abset).dipfit.coordformat , ...
           'normlen' ,'on', 'pointout' ,'on','color', {'b', 'r'}, 'dipnames', {comp_to_disp [ STUDY.cluster(cls).name ' mean' ] },...
            'spheres', 'on', 'verbose', 'off', 'projlines', 'on');
    else
       dipplot(cluster_dip_models, 'meshdata', ALLEEG(abset).dipfit.hdmfile, 'mri', ALLEEG(abset).dipfit.mrifile,'coordformat', ALLEEG(abset).dipfit.coordformat , ...
          'normlen' ,'on', 'pointout' ,'on','color', {'b', 'r'}, 'dipnames', {comp_to_disp [STUDY.cluster(cls).name ' mean']}, ...
          'spheres', 'on', 'verbose', 'off', 'projlines', 'on');
    end
    fig_h = gcf;
    set(fig_h,'Name', [subject ' / ' 'IC' num2str(comp) ', ' STUDY.cluster(cls).name],'NumberTitle','off');
end

% -----------------------
% load all dipoles and
% compute dipole centroid
% DEVELOPMENT: this function
% should be the only one to
% access dipole information
% -----------------------
function STUDY = std_centroid(STUDY,ALLEEG, clsind, tmp);

    for clust = 1:length(clsind)
        max_r = 0;
        len = length(STUDY.cluster(clsind(clust)).comps);
        tmppos = [ 0 0 0 ];
        tmpmom = [ 0 0 0 ];
        tmprv = 0;
        ndip = 0;
        for k = 1:len
            fprintf('.');
            comp  = STUDY.cluster(clsind(clust)).comps(k);
            abset = STUDY.cluster(clsind(clust)).sets(1,k);
            if ~isfield(ALLEEG(abset), 'dipfit')
               warndlg2(['No dipole information available in dataset ' num2str(abset) ], 'Aborting compute centroid dipole');
               return;
            end
            if ~isempty(ALLEEG(abset).dipfit.model(comp).posxyz)
                ndip   = ndip +1;
                posxyz = ALLEEG(abset).dipfit.model(comp).posxyz;
                momxyz = ALLEEG(abset).dipfit.model(comp).momxyz;
                if size(posxyz,1) == 2
                    if all(posxyz(2,:) == [ 0 0 0 ])
                        posxyz(2,:) = [];
                        momxyz(2,:) = [];
                    end;
                end;
                tmppos = tmppos + mean(posxyz,1);
                tmpmom = tmpmom + mean(momxyz,1);
                tmprv = tmprv + ALLEEG(abset).dipfit.model(comp).rv;
                if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                   if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                       load('-mat', ALLEEG(abset).dipfit.hdmfile);
                       max_r = max(max_r, max(vol.r));
                   else % old version of dipfit
                       max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
                   end
               end
            end
        end
        centroid{clust}.dipole.posxyz =  tmppos/ndip;
        centroid{clust}.dipole.momxyz =  tmpmom/ndip;
        centroid{clust}.dipole.rv =  tmprv/ndip;
        if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical') & (~isfield(ALLEEG(abset).dipfit, 'hdmfile')) %old dipfit
            centroid{clust}.dipole.maxr = max_r;
        end
        STUDY.cluster(clsind(clust)).dipole = centroid{clust}.dipole;
    end
    fprintf('\n');

% --------------------------------
% new function to compute centroid
% was programmed to debug the function
% above but is now used in the code
% --------------------------------
function dipole = computecentroid(alldipoles)

        max_r = 0;
        len = length(alldipoles);
        dipole.posxyz = [ 0 0 0 ];
        dipole.momxyz = [ 0 0 0 ];
        dipole.rv = 0;
        ndip = 0;
        count = 0;
        warningon = 1;
        for k = 1:len
            if size(alldipoles(k).posxyz,1) == 2
                if all(alldipoles(k).posxyz(2,:) == [ 0 0 0 ])
                    alldipoles(k).posxyz(2,:) = [];
                    alldipoles(k).momxyz(2,:) = [];
                end;
            end;
            if ~isempty(alldipoles(k).posxyz)
                dipole.posxyz = dipole.posxyz + mean(alldipoles(k).posxyz,1);
                dipole.momxyz = dipole.momxyz + mean(alldipoles(k).momxyz,1);
                dipole.rv     = dipole.rv     + alldipoles(k).rv;
                count = count+1;
            elseif warningon
                disp('Some components do not have dipole information');
                warningon = 0;
            end;
        end
        dipole.posxyz = dipole.posxyz/count;
        dipole.momxyz = dipole.momxyz/count;
        dipole.rv     = dipole.rv/count;
        if isfield(alldipoles, 'maxr')
            dipole.maxr = alldipoles(1).max_r;
        end;

author: PatternRecognition

repo: OpenBMI


text:
function dudt = burgers(t,u,epsilon,D1,D2)

dudt = epsilon*D2*u - u.*(D1*u);

author: Sable

repo: mcbench-benchmarks


text:
%% The Class MLSSolver
%
%%
%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
pb=1e-6;
data=[];
while pb<0.5
    data=[data;[pb,geterror(pb)]];
    pb=pb*1.5;
end
plot(log(data(:,1))/log(10),log(data(:,2))/log(10));


text:
function dn = dnmult(dn0,dn1)

% DMULT dual number multiplication
%
%   DN = DNMULT(DN0,DN1) returns the dual number DN which is the dual number
%     multiplication of the dual numbers DN0 and DN1
%      - DN0 (resp. DN1) is a dual number (DN0 = a0 + eps*b0, eps^2 = 0).
%         It is a 2-vector or a 2*N array (column i represents dual number
%         i) where N is the number of dual numbers. DN0 and DN1 must have
%         the same size.
%      - DN is a dual number. It is a 2*N array (each column is the dual
%          multiplication of the corresponding columns in DN0 and DN1)

s0 = size(dn0);
s1 = size(dn1);
if s0 == [1 2]
    dn0 = dn0.';
    s0 = size(dn0);
end
if s1 == [1 2]
    dn1 = dn1.';
    s1 = size(dn1);
end

% wrong format
if s0(1) ~= 2 || s1(1) ~= 2
    error('DualQuaternion:Dmult:wrongsize',...
        '%d rows in the DN0 array and %d rows in the DN1 array. It should be 2 for both.',...
        s0(1),s1(1));
end

% sizes do not match
n1 = s0(2);
n2 = s1(2);
if n1 ~= n2
    error('DualQuaternion:Dmult:wrongsize',...
        '%d dual numbers in DN0 array and %d dual numbers in DN1 array.They should be equal.',...
        n1,n2);
end

dn = sym(zeros(2,n1));
dn(1,:) = dn0(1,:).*dn1(1,:);
dn(2,:) = dn0(1,:).*dn1(2,:)+dn0(2,:).*dn1(1,:);


text:
function data = memmap_file_sequence(foldername)

% read a sequence of stacked tiff arrays, reshapes it to 2d array and saves it a mat file that can be memory mapped.
% The file is saved both in its original format and as a reshaped 2d matrix

% INPUTS
% foldername:     path to tiff folder containing a series of tiff files

% OUTPUT
% data:         object with the data containing:
%   Yr:         reshaped data file
% sizY:         dimensions of original size
%   nY:         minimum value of Y (not subtracted from the dataset)

% Author: Eftychios A. Pnevmatikakis, Simons Foundation, 2016

nY = Inf;
files = dir([foldername,'/*.tif']);
T = 0;
data = matfile([foldername,'/',files(1).name(1:end-8),'.mat'],'Writable',true);
tt1 = tic;
for i = 1:length(files)
    filename = [foldername,'/',files(i).name];
    Y = tiff_reader(filename);
    sizY = size(Y);
    Yr = reshape(Y,prod(sizY(1:end-1)),[]);
    if i == 1; data.Yr(size(Yr,1),size(Yr,2)*length(files)) = single(0); end
    nY = min(min(Yr(:)),nY);
    if length(sizY) == 3
        if i == 1; data.Y(sizY(1),sizY(2),sizY(3)*length(files)) = single(0); end
        data.Y(1:sizY(1),1:sizY(2),T+(1:size(Yr,2))) = single(Y);
    elseif legnth(sizY) == 4
        if i == 1; data.Y(sizY(1),sizY(2),sizY(3),sizY(4)*length(files)) = single(0); end
        data.Y(1:sizY(1),1:sizY(2),1:sizY(3),T+(1:size(Yr,2))) = single(Y);
    end
    data.Yr(1:prod(sizY(1:end-1)),T+(1:size(Yr,2))) = single(Yr);
    T = sizY(end) + T;
    toc(tt1);
end
sizY(end) = T;
data.sizY = sizY;
data.nY = nY;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
% Remove unwanted periods in positions, recreate time and move spike times accordingly
%
% This function can be used to remove some of position samples keeping spike timestamps
% in sync. One application is to remove positions that correspond to slow speed.
%
%  USAGE
%   [newPos, spkInd, newCounter] = general.trimPositionAndSpikes(pos, spikeCounter, <options>)
%   pos             NxM matrix of position samples in form [t x y ...]. The very first
%                   column must contain time.
%   spikeCounter    Nx1 vector, i.e. it contains the same number of samples as pos.
%                   Each element in spikeCounter indicates how many spikes fall to that
%                   particular position sample. spikeCounter(3) = 2 indicates that 2 spikes
%                   correspond to position sample 3. And spikeCounter(5) = 0 indicates that
%                   no spikes correspond to position sample 5.
%   <options>       Optional list of property-value pairs (see table below)
%
%   ==============================================================================================
%    Properties     Values
%   ----------------------------------------------------------------------------------------------
%    'trim'         Nx1 vector. This is a logical vector that indicates which position samples should
%                   be removed. If not provided, then trim = isnan(pos(:, 2)), i.e. trim vector
%                   will be indices of NaN values in second column of pos matrix.
%   ==============================================================================================
%   newPos          JxM matrix of trimmed positions in form [t x y ...]. Time column is recreated
%                   to be linear with sampling rate extracted from pos.
%   spkInd          Linear indices of spikes in newPos matrix, size Kx1.
%   newCounter      Trimmed version of spikeCounter, size Jx1.
%
function [newPos, spkInd, spikeCounter] = trimPositionAndSpikes(pos, spikeCounter, varargin)
    inp = inputParser;
    deafaultTrim = [];

    addRequired(inp, 'pos');
    addRequired(inp, 'spikeCounter');
    addParameter(inp, 'trim', deafaultTrim);
    parse(inp, pos, spikeCounter, varargin{:});
    badPos = inp.Results.trim;
    if isempty(badPos)
        badPos = isnan(pos(:, 2));
    else
        if numel(badPos) ~= size(pos, 1)
            error('BNT:arg:length', 'Parameter ''trim'' should have the same number of elements as number of rows in ''pos'' matrix');
        end
    end

    sampleTime = mean(diff(pos(:, bntConstants.PosT)));
    if isnan(sampleTime)
        sampleTime = data.sampleTime('sec');
    end

    newPos = pos;
    newPos(badPos, :) = [];
    spikeCounter(badPos, :) = [];

    % make new timestamps
    newPos(:, 1) = 0:sampleTime:(size(newPos, 1) - 1) * sampleTime;
    %spikes = repelem(newPos(:, bntConstants.PosT), spikeCounter);
    spkInd = repelem(1:size(newPos, 1), spikeCounter);
end

author: kavli-ntnu


text:
% PURPOSE: shows pop function ending statement (located at msg2end.txt)
%
% Format:
%
% msg2end
%
% See also readmsg2end.m msg2end.txt
%
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2011

function msg2end
[msg2show x colormsg] = readmsg2end;
try
        cprintf(colormsg, msg2show');
catch
        fprintf(msg2show);
end
fprintf('\n');

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%% patchEdgeCrossProduct
% Below is a demonstration of the features of the |patchEdgeCrossProduct| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=patchEdgeCrossProduct(F,V);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%--- help for dsge_solver_first_order_autoregress_1 ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named dsge_solver_first_order_autoregress_1
%
%       dsge/dsge_solver_first_order_autoregress_1
%

author: jmaih


text:
function weig = get_w(supVec,method)
%
%                w = get_w()
%
% Returns value of w from svm
%
% AE: i have changed this procedure to include non linear feature
% selection with RFE. When non-linear, it gives  a score vector for each
% feature (used in RFE) which is not the margin but is sorted as the margin

kerTemp = sv.child.kerparam;
if strcmp(sv.child.ker,'linear')
  weig=sv.alpha'*get_x(sv.Xsv);
  return;
end

if strcmp(sv.child.ker,'weighted_linear')
  numEx=get_dim(sv.Xsv);
  weig=sv.alpha'* (get_x(sv.Xsv) .*  repmat(kerTemp,numEx,1));
  return;
end
if strcmp(sv.child.ker,'rbf')  | strcmp(sv.child.ker,'poly'),
    xTemp = get_x(sv.Xsv);
    alphaTemp = sv.alpha;

    if strcmp(sv.child.ker,'rbf'),

      % compute the kernel matrix for all components
      K = xTemp*xTemp';
      Kdn = sum(xTemp.^2,2);
      Kn = sum(xTemp.^2,2);
      K = ones(size(xTemp,1),1)*Kdn' + Kn*ones(1,size(xTemp,1)) - 2*K;
      K = exp(-K/(2*kerTemp));
      % compute the margin when one component is removed
       for i = 1:size(xTemp,2),
           Ki = xTemp(:,i)*ones(1,size(xTemp,1)) - ones(size(xTemp,1),1)*xTemp(:,i)';
           Ki = Ki.^2;
           Ki = Ki/(2*kerTemp);
           Ki = exp(Ki);
           weig(i) = (alphaTemp'*(K.*Ki)*alphaTemp);
       end;
    elseif strcmp(sv.child.ker,'poly'),
        % compute the margin when one component is removed
        Ktmp = xTemp*xTemp';
        for i = 1:size(xTemp,2),
           Ki = xTemp(:,i)*xTemp(:,i)';
           K_i = (Ktmp - Ki+1).^(kerTemp);
           weig(i) = (alphaTemp'*K_i*alphaTemp);
    	end;
    end;% if strcmp(...,'rbf')
    weig=max(weig)-weig;  %% make largest the smallest --- wrong way round!
    return;
end;% if strcmp(... | strcmp(...,'poly')
%%% if this point is reached then the kernel is a non classic one
ker=supVec.child.ker;
if strcmp(ker,'custom'),
    error('Get w not implemented for CUSTOM kernels.');
    return;
end;
%%% if tpoly kernel -> used for nfe
if strcmp(ker,'tpoly'),
    [numEx,vDim,oDim]=get_dim(supVec.Xsv);
    xTemp = get_x(supVec.Xsv);
    alphaTemp = supVec.alpha;
    weig = ones(length(supVec.child.param,vDim));

    for j=1:length(supVec.child.kerparam),
        ktmp = supVec.child;
        kerp = ktmp.kerparam;
        for i = 1:length(kerp{j}),
            kerptmp=kerp;
            kerptmp{j} = [kerp{j}(1:i-1),kerp{j}(i+1:length(kerp{j}))];
            ktmp.kerparam=kerptmp;
            K_i = get_kernel(ktmp,supVec.Xsv,supVec.Xsv);
            weig(j,kerp{j}(i)) = (alphaTemp'*K_i*alphaTemp);
        end;
    end;
    weig = max(max(weig))-weig;
    return;
end;
for i = 1:size(xTemp,2),
    dtmp = data('tmp',[xTemp(:,1:(i-1)),xTemp(:,(i+1):size(xTemp,2))],[]);
    Ki = get_kernel(supVec.child,dtmp,dtmp);
    weig(i) = (alphaTemp'*Ki*alphaTemp);
end;
weig=max(weig)-weig;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function tf = isTearDownString(str)
%isTearDownString True if string looks like the name of a teardown function
%   tf = isTearDownString(str) returns true if the string str looks like the
%   name of a teardown function.  If str is a cell array of strings, then
%   isTearDownString tests each string in the cell array, returning the results
%   in a logical array with the same size as str.

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

setup_exp = '^[tT]ear[dD]own';
tf = mtest.utils.containsRegexp(str, setup_exp);

author: ricedsp


text:
classdef CEC2017_F20 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{12}.o;
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Y = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = Y.^2;
            PopObj = sum(0.5+(sin(sqrt(Y(:,1:end-1)+Y(:,2:end))).^2-0.5)./(1+0.001*sqrt(Y(:,1:end-1)+Y(:,2:end))).^2,2) + ...
                     0.5+(sin(sqrt(Y(:,end)+Y(:,1))).^2-0.5)./(1+0.001*sqrt(Y(:,end)+Y(:,1))).^2;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Y = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = sum(Y,2);
            PopCon(:,1) = cos(Y).^2 - 0.25*cos(Y) - 0.125;
            PopCon(:,2) = exp(cos(Y)) - exp(0.25);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% INTERSECT - SPHERE AND ROTATED CUBOID (OBB)
function [intersectFlag,separation_OBB]  = CheckOBBSphereIntersection(centerA,radiusA,centerB,cuboidVerticesB)
% This function computes the intersection between a sphere and a
% rotated bounding box (OBB). This function assumes that the
% problem is being resolved in the axes of A (i.e) 'R' is the
% rotation matrix taking the cube from its own frame to the
% axes of A.

% RELATIVE POSITION RAY
[BAray] = ray(centerB,(centerA-centerB));       % Vector between
% If we define A to be sphere, its projection on the separation
% vector will always be its radius
projection_A = radiusA;
% GET THE PROJECTION OF ALL THE VERTICES ON THE SEPARATING VECTOR
projection_B = zeros(size(cuboidVerticesB,1),1);
for i = 1:size(cuboidVerticesB,1)
    % Its projection on the separation vector
    [projection_B(i)] = GetPointProjectionOnRay(BAray,cuboidVerticesB(i,:)');
end
% THE MAXIMAL PROJECTION TOWARDS THE SPHERE
projection_B = max(projection_B);
% CHECK IF THIS EXCEEDS THE SEPARATION OF THE TWO OBJECTS
separation_OBB = norm(centerB-centerA) - (projection_A + projection_B);
% EVALUATE INTERSECTION
intersectFlag = 0;
if separation_OBB < 0
    separation_OBB = 0;
    intersectFlag = 1;
end
end

repo: OpenMAS


text:
function sparse_test06 ( )

%*****************************************************************************80
%
%% SPARSE_TEST06 times the zeroing out operation for full and sparse matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_TEST06:\n' );
  fprintf ( 1, '  Zero out a portion of a matrix.\n' );
  fprintf ( 1, '  Compare the time required when using full or sparse storage.\n' );
  fprintf ( 1, '  The sparse matrix takes longer to modify, and it takes longer\n' );
  fprintf ( 1, '  when there are more elements to remove.\n' );

  n = 2000;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix is of size %d by %d\n', n, n );
  fprintf ( 1, '  Initial number of nonzeros is %d\n', n * n );

  for amax = [ 1, n * n / 2, n * n - 1 ]

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Zero out all entries greater than %d\n', amax );

    A = magic ( n );
    tic
    A(A>amax) = 0.0;
    t = toc;
    fprintf ( 1, '  Full storage matrix required %g seconds\n', t );

    A = magic ( n );
    A = sparse ( A );
    tic
    A(A>amax) = 0.0;
    t = toc;
    fprintf ( 1, '  Sparse storage matrix required %g seconds\n', t );
    fprintf ( 1, '  Number of nonzeros is %d\n', nnz ( A ) );

  end

  return
end


text:
alpha=1;
beta=0.5;
gamma=0.00001;
lambda=0.00001;
numK=10;
maxIter=100;
flag=1;
inputPath='D:/chengxh/matlab/TLLibrary64/CRA/data/img/inputData.mat';
% CRA_enterFunc(alpha,beta,gamma,lambda,numK,maxIter,flag,inputPath)

alpha=1;
beta=0.5;
gamma=0.00001;
numK=10;
numC=2;
maxIter=100;
flag=1;
TrainXPath='D:/chengxh/matlab/TLLibrary/TLDA/data/img/TrainData.mat';
TrainYPath='D:/chengxh/matlab/TLLibrary/TLDA/data/img/TrainLabel.mat';
TestXPath='D:/chengxh/matlab/TLLibrary/TLDA/data/img/TestData.mat';
TestYPath='D:/chengxh/matlab/TLLibrary/TLDA/data/img/TestLabel.mat';
InputPath='TLDA/data/img/inputData.mat'
% TLDA_enterFunc(alpha,beta,gamma,numK,maxIter,flag,InputPath)

alpha=1;
beta=0.5;
gamma=0.00001;
lambda=0.00001;
numK=10;
sourceSize='166 200 185';
maxIter=10;
flag=1;
TrainXPath='D:/chengxh/matlab/TLLibrary/CRA/data/img/TrainData.mat';
TrainYPath='D:/chengxh/matlab/TLLibrary/CRA/data/img/TrainLabel.mat';
TestXPath='D:/chengxh/matlab/TLLibrary/CRA/data/img/TestData.mat';
TestYPath='D:/chengxh/matlab/TLLibrary/CRA/data/img/TestLabel.mat';
%CRA_enterFunc(alpha,beta,gamma,lambda,numK,sourceSize,maxIter,flag,TrainXPath,TrainYPath,TestXPath,TestYPath)

alpha=2.4;
beta=2.4;
numCluster=15;
maxIter=200;
TrainXPath='MTrick/data/TrainData.mat';
TrainYPath='MTrick/data/TrainLabel.mat';
TestXPath='MTrick/data/TestData.mat';
TestYPath='MTrick/data/TestLabel.mat';
inputPath='data.mat';
% result = MTrick_enterFunc(alpha,beta,numCluster,maxIter,inputPath);

numIdentical=20;
numAlike=20;
numDistinct=10;
numIter=10;
inputPath='data.mat';
HIDC_enterFunc(numIdentical,numAlike,numDistinct,numIter,inputPath)

numIdentical=20;
numAlike=20;
numDistinct=10;
numIter=100;
sourceSize = '1976';
targetSize = '1977';
TrainXPath='TriTL/data/TrainData.mat';
TrainYPath='TriTL/data/TrainLabel.mat';
TestXPath='TriTL/data/TestData.mat';
TestYPath='TriTL/data/TestLabel.mat';
InputPath='data.mat';
% TriTL_enterFunc(numIdentical,numAlike,numDistinct,numIter,InputPath)

numCluster=64;
maxIter=10;
inputPath='data.mat';
% CDPLSA_enterFunc(numCluster,maxIter,inputPath)

numIdentical=20;
numAlike=20;
numDistinct=10;
numIter=500;
inputPath='data.mat';
% HIDC_enterFunc(numIdentical,numAlike,numDistinct,numIter,inputPath)

alpha=0.4;
beta=15;
gamma=0.12;
threshold=0.1;
inputPath='D:/chengxh/matlab/TLLibrary64/IHR/data/inputData.mat';
neighborPath='D:/chengxh/matlab/TLLibrary64/IHR/data/KneighborFile.txt';
%IHR_enterFunc(alpha,beta,gamma,threshold,inputPath,neighborPath)

inputPath='D:/chengxh/matlab/TLLibrary64/MAMUDA/data/inputData.mat';
ite1=4;
ite2=4;
itermediateD=50;
reducedD=50;
sharedD=50;
% meanF1=MAMUDA_enterFunc(inputPath,ite1,ite2,itermediateD,reducedD,sharedD);

inputPath='D:/chengxh/matlab/TLLibrary64/CSL_MTMV/data/inputData.mat';
ite=20;
h=15;
alpha=0.0039;
beta=0.0078;
gamma=0.0039;
% meanF1=CASO_MVMT_enterFunc(inputPath,ite,h,alpha,beta,gamma);

inputPath='CCR3/data/inputData.mat';
gamma=150;
%results=CCR3_enterFunc(gamma, inputPath);

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%-------------------------------------------------------------------------------------------------------------
% This is an implementation of the MCWNNM algorithm for real color image denoising.
% Author:  Jun Xu, csjunxu@comp.polyu.edu.hk
%              The Hong Kong Polytechnic University
%
% Please refer to the following paper if you use this code:
%
% @article{MCWNNM,
% 	author = {Jun Xu and Lei Zhang and David Zhang and Xiangchu Feng},
% 	title = {Multi-channel Weighted Nuclear Norm Minimization for Real Color Image Denoising},
% 	journal = {ICCV},
% 	year = {2017}
% }
%
% Please see the file License.txt for the license governing this code.
%-------------------------------------------------------------------------------------------------------------

clear;
Original_image_dir  =    'kodak_color/';
fpath = fullfile(Original_image_dir, '*.png');
im_dir  = dir(fpath);
im_num = length(im_dir);

nSig = [5 30 15];

Par.nSig      =   nSig;                                 % Variance of the noise image
Par.win =   20;                                   % Non-local patch searching window
Par.Constant         =  2 * sqrt(2);                              % Constant num for the weight vector
Par.Innerloop =   2;                                    % InnerLoop Num of between re-blockmatching
Par.ps       =   6;                            % Patch size
Par.step        =   5;
Par.Iter          =   6;                            % total iter numbers
Par.display = true;

% Par.method = 'WNNM_ADMM'
Par.method = 'MCWNNM_ADMM'
Par.maxIter = 10;

Par.model = '2';

Par.delta     =   0.1;                                  % Parameter between each iter
Par.lambda = 0.75;
Par.mu = 1.001;
Par.rho = 0.05;
% record all the results in each iteration
Par.PSNR = zeros(Par.Iter, im_num, 'single');
Par.SSIM = zeros(Par.Iter, im_num, 'single');
for i = 1:im_num
    Par.image = i;
    Par.nSig0 = nSig;
    Par.nlsp        =   70;   % Initial Non-local Patch number
    Par.I =  double( imread(fullfile(Original_image_dir, im_dir(i).name)) );
    S = regexp(im_dir(i).name, '\.', 'split');
    [h, w, ch] = size(Par.I);
    Par.nim = zeros(size(Par.I));
    for c = 1:ch
        randn('seed',0);
        Par.nim(:, :, c) = Par.I(:, :, c) + Par.nSig0(c) * randn(size(Par.I(:, :, c)));
    end
    fprintf('%s :\n',im_dir(i).name);
    PSNR =   csnr( Par.nim, Par.I, 0, 0 );
    SSIM      =  cal_ssim( Par.nim, Par.I, 0, 0 );
    fprintf('The initial value of PSNR = %2.4f, SSIM = %2.4f \n', PSNR,SSIM);
    %
    time0 = clock;
    if Par.model == '1'
        [im_out, Par] = MCWNNM_ADMM1_Denoising( Par.nim, Par.I, Par );
    elseif Par.model == '2'
        [im_out, Par] = MCWNNM_ADMM2_Denoising( Par.nim, Par.I, Par );
    else
        [im_out, Par] = MCWNNM_ADMM_Denoising( Par.nim, Par.I, Par );
    end
    fprintf('Total elapsed time = %f s\n', (etime(clock,time0)) );
    im_out(im_out>255)=255;
    im_out(im_out<0)=0;
    % calculate the PSNR
    Par.PSNR(Par.Iter, Par.image)  =   csnr( im_out, Par.I, 0, 0 );
    Par.SSIM(Par.Iter, Par.image)      =  cal_ssim( im_out, Par.I, 0, 0 );
    imname = sprintf([Par.method '_nSig' num2str(nSig(1)) num2str(nSig(2)) num2str(nSig(3)) '_' Par.model '_Oite' num2str(Par.Iter) '_Iite' num2str(Par.maxIter) '_rho' num2str(Par.rho) '_mu' num2str(Par.mu) '_lambda' num2str(Par.lambda) '_' im_dir(i).name]);
    imwrite(im_out/255, imname);
    fprintf('%s : PSNR = %2.4f, SSIM = %2.4f \n',im_dir(i).name, Par.PSNR(Par.Iter, Par.image),Par.SSIM(Par.Iter, Par.image)     );
end
mPSNR=mean(Par.PSNR,2);
[~, idx] = max(mPSNR);
PSNR =Par.PSNR(idx,:);
SSIM = Par.SSIM(idx,:);
mSSIM=mean(SSIM,2);
fprintf('The best PSNR result is at %d iteration. \n',idx);
fprintf('The average PSNR = %2.4f, SSIM = %2.4f. \n', mPSNR(idx),mSSIM);
name = sprintf([Par.method '_' Par.model '_nSig' num2str(nSig(1)) num2str(nSig(2)) num2str(nSig(3)) '_' Par.model '_Oite' num2str(Par.Iter) '_Iite' num2str(Par.maxIter) '_rho' num2str(Par.rho) '_mu' num2str(Par.mu) '_lambda' num2str(Par.lambda) '.mat']);
save(name,'nSig','PSNR','SSIM','mPSNR','mSSIM');

author: csjunxu


text:
%% gaborFilter
% Below is a demonstration of the features of the |gaborFilter| function

%%
clear; close all; clc;

%% Syntax
% |GB=gaborFilter(Vr,S,f,tagShift,d,siz,p);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.


text:
classdef dl_list_node < handle

% Copyright (C) 2012 Jean-Pierre de la Croix
% see the LICENSE file included with this software

    properties
        key_
        next_
        prev_
    end

    methods
        function obj = dl_list_node(key)
            obj.key_ = key;
            obj.next_ = [];
            obj.prev_ = obj.next_;
        end
    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
%% Example 8.15: Duffing van der Pol oscillator
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%% Duffing van der Pol

  % Time-span
  tspan = 0:2^-5:20;

  % Parameters
  alpha = 1;

  % Define arrow (for visalization)
  arrow1 = [-1 1 0 -1; -.5 -.5 2 -.5]';

  % The model
  f = @(x,t) [x(2,:); x(1,:).*(alpha - x(1,:).^2)-x(2,:)];
  L = @(x,t) [zeros(1,size(x,2)); x(1,:)];

%% ODE

  figure(1); clf; hold on

  for j=1:10
    %x = rk4(f,tspan,[-2-.2*j; 0]);
    x = rk4simple(f,tspan,[-2-.2*j; 0]);
    %[~,x] = ode45(@(t,x) f(x,t),tspan,[-2-.2*j; 0]); x = x';

    % Plot trajectory
    plot(x(1,:),x(2,:),'-k','LineWidth',.25)

    % Plot direction
    uv = f(x,[]); ind = 10;
    newquiver(x(1,ind),x(2,ind),uv(2,ind),-uv(1,ind), ...
      'X',arrow1,'scale',.04*[1 14.8/8.8])

  end

  % Set axis limits
  axis([-4.4 4.4 -4.8 10]), axis square
  %set(gca,'XTick',-4:2:4,'YTick',-4:2:8)
  xlabel('$x_1$'); ylabel('$x_2$')
  box on

%% SDE

  figure(2); clf; hold on
  figure(3); clf; hold on

  for j=1:10

    figure(2);

    % Lock random seed
    if exist('rng') % Octave doesn't have rng
      rng(3,'twister')
    else
      randn('state',2);
      rand('state',2);
    end

    % Use the strong order 1.0 method
    x = srkS10scalarnoise(f,L,tspan,[-2-.2*j; 0],.5^2);

    % Plot trajectory
    plot(x(1,:),x(2,:),'-k','LineWidth',.25)

    % Plot direction
    uv = f(x,[]); ind = 10;
    newquiver(x(1,ind),x(2,ind),uv(2,ind),-uv(1,ind), ...
      'X',arrow1,'scale',.04*[1 14.8/8.8])

    figure(3);
    plot(tspan,x(1,:),'-k')
    plot(tspan,x(2,:),'-','Color',[.7 .7 .7])

  end

  % Set axis limits
  figure(2)
  axis([-4.4 4.4 -4.8 10]), axis square
  %set(gca,'XTick',-4:2:4,'YTick',-4:2:8)
  xlabel('$x_1$'); ylabel('$x_2$')
  box on

  % Set axis limits
  figure(3)
  xlim([0 20])
  xlabel('Time, $t$'); ylabel('$x$')
  legend('$x_1(t)$','$x_2(t)$')
  box on

%% Weak approximation

  % Time-span
  tspan = 0:2^-4:20;

  % Reset random seed
  % Lock random seed
  if exist('rng') % Octave doesn't have rng
    rng(1,'twister')
  else
    randn('state',2);
    rand('state',2);
  end

  % Initial point
  x0 = [-3;0];

  % Number of smaples
  x = zeros(2,10000);

  % Iterate
  for j=1:size(x,2)

    % Weak SRK scheme
    foo = srkW20(f,L,tspan,x0,.5^2,true);

    % Store
    x(:,j) = foo(:,end);

    % Report
    if rem(j,100)==0,
      figure(4); clf
        hist(x(1,1:j),ceil(sqrt(j)))
        drawnow
      j
    end

  end

%% Histogram

  % Bins for histogram
  t = linspace(min(x(1,:)),max(x(1,:)),64);

  figure(5); clf

    % Show solution w2.0
    n = histc(x(1,:),t);
    stairs(t-(t(2)-t(1))/2,n/size(x,2),'-k')

    % Label
    xlabel('$x_1$')

    % Ticks
    box off
    xlim([-2.2 2.2])
    %set(gca,'XTick',0:.2:1.2)

  figure(6); clf

    % Show solution w2.0
    plot(x(1,:),x(2,:),'.k')

    % Label
    xlabel('$x_1$')
    ylabel('$x_2$')

    % Ticks
    box on

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
% test for texture synthesis

path(path, 'images/');
path(path, '../images/');

name = 'grass';
name = 'turbulence1';
name = 'fingerprint3';
name = 'dunes';

n = 128;
M = load_image(name, n);
M = crop(M,n);

s = size(M,3);
n1 = 128;
M1 = randn(n1,n1,s);

options.niter_synthesis = 5;
options.verb = 1;
disp('--> Steerable synthesis.');
options.synthesis_method = 'steerable';
Msteer = perform_wavelet_matching(M1,M,options);
disp('--> Wavelet orthogonal synthesis.');
options.synthesis_method = 'wavelets-ortho';
Morth = perform_wavelet_matching(M1,M,options);
disp('--> Quincunx synthesis.');
options.synthesis_method = 'quincunx-ti';
Mquin = perform_wavelet_matching(M1,M,options);
disp('--> Wavelets TI synthesis.');
options.synthesis_method = 'wavelets-ti';
Mwavti = perform_wavelet_matching(M1,M,options);

clf;
imageplot({M Msteer Morth Mquin Mwavti},{'Original', 'Synth. steerable', 'Synth. wav. ortho', 'Synth. quincunx', 'Synth. wav. ti'});

author: gpeyre

repo: matlab-toolboxes


text:
%% Calculating optical flow of a sequence of images or from Video file
%% Author Gurkirt Singh
%% Date 09 June 2011 at Vision Lab IIT DELHI
%% Reading a Video file
clear all
clc
Object=mmreader('1.avi');
lastframe=read(Object, inf);
NumFrames = Object.NumberOfFrames;
%% Initialization of paramEters
%  Alpha
alpha=10;
% Number of iterations.
iterations=3;
% Do you Want to reduce size of Image for less computation? If yes DoYou=1;
DoYou=1;
% Width and Height
U=cell(NumFrames/2-1,1);
V=cell(NumFrames/2-1,1);
n=0;
% Normal Intensity flow

Normal=cell(NumFrames-1,1);

Window=[1/12 1/6 1/12;1/6 0 1/6;1/12 1/6 1/12];

%% Main program Starts here
hsize=[5 5];
sigma=1;
I=rgb2gray(read(Object,1));
h = fspecial('gaussian',hsize,sigma);
% I = imfilter(I,h,'replicate');
PreviousFrame= imfilter(I,h,'replicate');
if DoYou==1
PreviousFrame=impyramid(PreviousFrame,'reduce');
PreviousFrame=impyramid(PreviousFrame,'reduce');
end
[height width]=size(PreviousFrame);
%%
for k = 2:1:NumFrames

    I=rgb2gray(read(Object,k));
%     I=imfilter(I,h,'replicate');
    CurrentFrame=imfilter(I,h,'replicate');
    if DoYou==1;
    CurrentFrame=impyramid(CurrentFrame,'reduce');
    CurrentFrame=impyramid(CurrentFrame,'reduce');
    end

    [TempU TempV TempNormal]=Opticalflow(CurrentFrame,PreviousFrame,alpha,iterations,height,width);
    n=n+1;
    U{n,1}=single(TempU);
    V{n,1}=single(TempV);
%     Normal{n,1}=single(TempNormal);
    PreviousFrame=CurrentFrame;

end
%% Last Part
% myObj = VideoWriter('Ite3.wmv');
% myObj.FrameRate = 10;
% open(myObj);

% for i=1:49-1
% %     u=U{i,1};
% %     v=V{i,1};
% %     temp=mat2gray(sqrt(v.^2+u.^2));
% temp=mat2gray(abs(Normal{i,1}));
%     imshow(temp);
% %     BW =im2bw(I,0.12);
%     currFrame = getframe;
%     writeVideo(myObj,currFrame);
% end
% close(myObj)


text:
function obj = as_channeltag(obj)
   switch(class(obj))
      case 'ChannelTag' %good. do nothing.
      case 'scnlobject'
         obj = get(obj, 'channeltag');
      otherwise
         obj = ChannelTag(obj); %attempt natural conversion
         % should be able to handle 'N.S.L.C'
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function [acc, rt] = FDDL_wrapper(Y_train, label_train, Y_test , label_test, ...
                            k, lambda1, lambda2)
% function [acc, rt] = FDDL_wrapper(Y_train, label_train, Y_test , ...
%       label_test, k, lambda1, lambda2)
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    if nargin == 0 % test mode
        dataset = 'myYaleB';
        N_train = 10;
        [~, Y_train, Y_test, label_train, label_test] = ...
            train_test_split(dataset, N_train);
        k = 8;
        lambda = 0.001;
        eta = 0.01;
    end
    C                = max(label_train);
    k0               = 0;
    opts.k           = k;
    opts.k0          = 0;
    opts.show_cost   = 0;
    opts.lambda1     = lambda1;
    opts.lambda2     = lambda2;
    opts.lambda3     = 0;
    opts.D_range     = k*(0:C);
    opts.D_range_ext = [opts.D_range k*C+k0];
    opts.initmode    = 'normal';
    opts.max_iter    = 100;
    opts             = initOpts(opts);
    opts.verbose      = true;
    opts.tol         = 1e-8;
    %% Train
    [D, ~, ~, ~, CoefM, ~, opts, rt] = ...
                    LRSDL(Y_train, label_train, opts);
    Y_range = label_to_range(label_train);
    C = max(label_train);
    opts.verbose = 0;
    opts.weight = 0.1;
    acc = [];
    for vgamma = [0.0001, 0.001, 0.01, 0.1]
        opts.gamma = vgamma;
        pred = FDDL_pred(Y_test, D, CoefM, opts);
        acc1 = double(numel(find(pred == label_test)))/...
            numel(label_test);
        fprintf('gamma = %.4f, acc = %.4f\n', vgamma, acc1);
        acc = [acc acc1];
    end
    acc = max(acc);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [Lrgb bgm] = watershed_seg(ip_wshed,ip_img)
% Watershed segmentation algorithm

se = strel('disk', 5); % disk = 20
%Opening
Io = imopen(ip_img, se);

%Opening-by-reconstruction (Iobr)
Ie = imerode(ip_img, se);
Iobr = imreconstruct(Ie, ip_img);

%Opening-closing (Ioc)
Ioc = imclose(Io, se);

%Opening-closing by reconstruction (Iobrcbr)
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);

%Regional maxima of opening-closing by reconstruction (fgm)
fgm = imregionalmax(Iobrcbr);

%Regional maxima superimposed on original image (I2)
I2 = ip_img;
I2(fgm) = 255;

se2 = strel(ones(3, 3)); % ones(5,5)
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

fgm4 = bwareaopen(fgm3, 20); %(fgm3, 20)

%Modified regional maxima superimposed on original image (fgm4)
I3 = ip_img;
I3(fgm4) = 255;

%Thresholded opening-closing by reconstruction (bw)
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));

%Watershed ridge lines (bgm)
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;

gradmag2 = imimposemin(ip_wshed, bgm | fgm4);

L = watershed(gradmag2);

%Markers and object boundaries superimposed on original image (I4)
I4 = ip_img;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;

Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');

author: Sable


text:
function varargout = efficency(X,P)
% function [con_eff,hrf_eff] = efficency(X,PARAMS)
% Tor Wager, 9/30/02
%
% a useful summary function for a design, incorporating tests of contrasts and of HRF shape estimation
%
% PARAMS structure contains inputs
% Optional fields:
%   contrasts   required for testing contrast efficiency; see calcEfficiency.m
%   Vi          intrinsic autocorrelation matrix; empty assumes identity matrix (independence)
%   S           smoothing matrix; empty for no smoothing
%   HRF         canonical hemodynamic response function; empty assumes SPM's canonical HRF
%               in this function, FIR estimation length fixed at 12 s.
%   nonlint     saturation threshold (ceiling) for predictors; simple nonlinearity model
%
% Required fields for condition function (vector) input:
%   ISI         sampling resolution of design vector, as time between elements in s
%   TR          sampling resolution of final design matrix (volumes acquired), in s
%
% Functions Called: (from OptimizeDesign GA toolbox)
% designvector2model.m
% tor_make_deconv_mtx2.m
% calcEfficiency.m

hrfsec = 30;

if ~isfield(P,'contrasts'), P.contrasts = []; end

if ~isfield(P,'S'), P.S = []; end

if any(size(X) == 1)    % then it's a condition function

    if ~isfield(P,'ISI'), error('ISI is required field in PARAMS input'), end
    if ~isfield(P,'TR'), error('TR is required field in PARAMS input'), end
    if ~isfield(P,'Vi'), P.Vi = eye(size(des,1)); end

    if ~isfield(P,'HRF'), P.HRF = spm_hrf(.1); end
    if ~isfield(P,'nonlint'), P.nonlint = []; end

    numsamps = ceil(length(X)*P.ISI/P.TR);
    X = designvector2model(X,P.ISI,P.HRF,P.TR,numsamps,P.nonlint,P.S);

    delta = [];
    for i = 1:max(X)
        delta(:,i) = (X == i);
    end
else                    % then it's a design matrix
end

svi = P.S * P.Vi;

if ~isempty(P.contrasts)
        % -------------------------------------------------------------------------------------------------
		% * efficiency
		% -------------------------------------------------------------------------------------------------

        	xtxitx = pinv(X);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)

            contrastweights = ones(1,size(xtxitx,1)-1);
			[dummy,varargout{1}] = calcEfficiency(contrastweights,P.contrasts,pinv(X),svi);

else

        [eff, eff_vector] = calcEfficiency([], [], pinv(X), svi);
        varargout{1} = eff_vector;
end

if isfield(P,'ISI') & isfield(P,'TR') & nargout > 1 & exist('delta', 'var')

        % -------------------------------------------------------------------------------------------------
		% * HRF shape estimation efficiency
		% -------------------------------------------------------------------------------------------------

			[X2] = tor_make_deconv_mtx2(delta,round(hrfsec / P.TR),P.TR / P.ISI);
            if ~isempty(P.S), X2 = P.S * X2;,end

        	xtxitx = pinv(X2);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
			[dummy,varargout{2}] = calcEfficiency([],[],xtxitx,svi);

elseif nargout > 1
            disp(['Missing ISI or TR field, or already-constructed model entered: HRF efficiency not calculated.'])
            varargout{2} = [];
end

return


text:
function [tcpipClient] = createTcpIpClient(varargin)
%createTcpIpClient Summary of this function goes here
%   Detailed explanation goes here
    global rHost;

    if(isempty(varargin))
        port = 8282;
        role = 'Client';
        if(isempty(rHost))
            rHost = 'localhost';
        end
    else
        port = varargin{1};
        role = varargin{2};
        try
            rHost = varargin{3};
        catch
            if(isempty(rHost))
                rHost = 'localhost';
            end
        end
    end

    tcpipClient = tcpip(rHost,port,'NetworkRole',role);

    tcpipClient.ByteOrder = 'bigEndian';
    set(tcpipClient, 'InputBufferSize',120480);
    set(tcpipClient, 'OutputBufferSize',120480);
    set(tcpipClient, 'Timeout',10);
    set(tcpipClient, 'ReadAsyncMode', 'continuous');
    set(tcpipClient, 'BytesAvailableFcnMode', 'byte');
    set(tcpipClient, 'BytesAvailableFcnCount', 120240);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef MatrixVectorizedInverter_3x3 < MatrixVectorizedInverter_Interface

    methods (Access = public)

        function B = computeInverse(obj,A)
            d = obj.computeAdjointsTranspose(A);

            det = obj.computeDeterminant(A);

            B = zeros(size(A));
            for i = 1:3
                for j = 1:3
                    B(i,j,:) = squeeze(d(i,j,:))./det;
                end
            end
        end

        function detA = computeDeterminant(obj,A)
            d    = obj.computeAdjointsTranspose(A);
            d = permute(d,[2 1 3]);
            detA = squeeze(sum(A.*d,[1 2]))/3;
        end

    end

    methods (Access = private, Static)

        function d = computeAdjointsTranspose(A)
            d = zeros(size(A));

            A11 = A(1,1,:);
            A12 = A(1,2,:);
            A13 = A(1,3,:);
            A21 = A(2,1,:);
            A22 = A(2,2,:);
            A23 = A(2,3,:);
            A31 = A(3,1,:);
            A32 = A(3,2,:);
            A33 = A(3,3,:);

            d(1,1,:) = A22.*A33-A23.*A32;
            d(1,2,:) = A32.*A13-A33.*A12;
            d(1,3,:) = A12.*A23-A13.*A22;

            d(2,1,:) = A23.*A31-A21.*A33;
            d(2,2,:) = A33.*A11-A31.*A13;
            d(2,3,:) = A13.*A21-A11.*A23;

            d(3,1,:) = A21.*A32-A22.*A31;
            d(3,2,:) = A31.*A12-A32.*A11;
            d(3,3,:) = A11.*A22-A12.*A21;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan