text:
classdef TestImageCoordinateUtilities < CoreTest
    % TestImageCoordinateUtilities. Tests for the MimImageCoordinateUtilities class.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = TestImageCoordinateUtilities
            mock_reporting = CoreMockReporting;

            obj.TestGetDimensionPermutationVectorFromDicomOrientation(mock_reporting);
            obj.TestGetDimensionIndicesAsArray(mock_reporting);
            obj.TestGetDimensionIndicesFromOrientations(mock_reporting);
        end

        function TestGetDimensionPermutationVectorFromDicomOrientation(obj, mock_reporting)

            % Test GetDimensionPermutationVectorFromDicomOrientation
            obj.Assert(isequal(MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([1 0 0 0 1 0]', mock_reporting), [1, 2, 3]), 'Expected result');
            obj.Assert(isequal(MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([1 0 0 0 0 1]', mock_reporting), [3, 2, 1]), 'Expected result');
            obj.Assert(isequal(MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([0 1 0 0 0 1]', mock_reporting), [2, 3, 1]), 'Expected result');
            obj.Assert(isequal(MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation([0.9964, -0.0851, 0, 0, 0, -1.0000]', mock_reporting), [3, 2, 1]), 'Expected result');

        end

        function TestGetDimensionIndicesAsArray(obj, mock_reporting)

            % Test with rows
            obj.Assert(isequal([1,2,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([1,0,0], [0,1,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,1,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,1,0], [1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,1,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,1], [1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([1,3,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([1,0,0], [0,0,1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,3,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,1,0], [0,0,1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,2,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,1], [0,1,0], mock_reporting)), 'Expected result');

            % Test with columns
            obj.Assert(isequal([1,2,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([1,0,0]', [0,1,0]', mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,1,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,1,0]', [1,0,0]', mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,1,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,1]', [1,0,0]', mock_reporting)), 'Expected result');
            obj.Assert(isequal([1,3,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([1,0,0]', [0,0,1]', mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,3,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,1,0]', [0,0,1]', mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,2,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,1]', [0,1,0]', mock_reporting)), 'Expected result');

            % Test with negatives
            obj.Assert(isequal([1,2,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([-1,0,0], [0,-1,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,1,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,-1,0], [-1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,1,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,-1], [-1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([1,3,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([-1,0,0], [0,0,-1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,3,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,-1,0], [0,0,-1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,2,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,-1], [0,-1,0], mock_reporting)), 'Expected result');

            obj.Assert(isequal([1,2,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([-1,0,0], [0,1,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,1,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,-1,0], [1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,1,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,-1], [1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([1,3,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([-1,0,0], [0,0,1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,3,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,-1,0], [0,0,1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,2,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,-1], [0,1,0], mock_reporting)), 'Expected result');

            obj.Assert(isequal([1,2,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([1,0,0], [0,-1,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,1,3],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,1,0], [-1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,1,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,1], [-1,0,0], mock_reporting)), 'Expected result');
            obj.Assert(isequal([1,3,2],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([1,0,0], [0,0,-1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([2,3,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,1,0], [0,0,-1], mock_reporting)), 'Expected result');
            obj.Assert(isequal([3,2,1],TestImageCoordinateUtilities.GetDimensionIndicesAsArray([0,0,1], [0,-1,0], mock_reporting)), 'Expected result');
        end

        function TestGetDimensionIndicesFromOrientations(obj, mock_reporting)
            obj.Assert(isequal([1,2,3,0,0,1],TestImageCoordinateUtilities.GetCombinedDimensionPermutationVectorFromDicomOrientation([1,0,0,0,1,0], mock_reporting)), 'Expected result');
        end
    end

    methods (Static, Access = private)
        function dims_array = GetDimensionIndicesAsArray(v1, v2, reporting)
            [d1, d2, d3] = MimImageCoordinateUtilities.GetDimensionIndicesFromOrientations(v1, v2, reporting);
            dims_array = [d1, d2, d3];
        end

        function combined_array = GetCombinedDimensionPermutationVectorFromDicomOrientation(orientation, reporting)
            [permutation_vector, flip] = MimImageCoordinateUtilities.GetDimensionPermutationVectorFromDicomOrientation(orientation, reporting);
            combined_array = [permutation_vector, flip];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function fx = p49_fun ( n, x )

%*****************************************************************************80
%
%% P49_FUN evaluates the integrand for problem 49.
%
%  Discussion:
%
%    The function is singular at two internal points, 1 and sqrt(2).
%
%  Interval:
%
%    0 <= x <= 3
%
%  Integrand:
%
%    x^3 * log ( abs ( ( x^2 - 1 ) * ( x^2 - 2 ) ) )
%
%  Exact Integral:
%
%    61 log ( 2 ) + (77/4) log ( 7 ) - 27 = 52.7408...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 104.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  x = x ( : );

  fx(1:n,1) = abs ( ( x(1:n,1).^2 - 1.0 ) .* ( x(1:n,1).^2 - 2.0 ) );
  i = ( fx == 0.0 );
  fx = x.^3 .* log ( fx );
  fx(i) = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = comp_idtwfb(c,nodes,dualnodes,Lc,rangeLoc,rangeOut,ext,do_complex)

if do_complex
   % Split the coefficients
   c1 = cellfun(@(cEl1,cEl2) (cEl1 + cEl2)/2, c(1:end/2),c(end:-1:end/2+1),...
             'UniformOutput',0);
   c2 = cellfun(@(cEl1,cEl2) (-1i*cEl1 + 1i*cEl2)/2, c(1:end/2),c(end:-1:end/2+1),...
             'UniformOutput',0);
else
   c1 = cellfun(@real,c,'UniformOutput',0);
   c2 = cellfun(@imag,c,'UniformOutput',0);
end

f1 = comp_iwfbt(c1,nodes,Lc,rangeLoc,rangeOut,ext);
f = f1 + comp_iwfbt(c2,dualnodes,Lc,rangeLoc,rangeOut,ext);

if ~do_complex
   f = real(f);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% CELLTOMAT - convert cell array to matrix
%
% Usage: >> M = celltomat( C );
%
% Author: Arnaud Delorme, CNL / Salk Institute, Jan 25 2002
%
% Note: This function overloads the neuralnet toolbox function CELLTOMAT,
% but does not have all its capacities. Delete this version if you have
% the toolbox.

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

function M = celltomat( C, varargin );

if nargin < 1
	help celltomat;
	return;
end

M = zeros(size(C));
for i=1:size(C,1)
    for j=1:size(C,2)
        M(i,j) = C{i,j};
    end
end;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function phi=sinCosConformLat2EllipsLat(sinChi,cosChi,f)
%%SINCOSCONFORMLAT2ELLIPSLAT Given the sine and cosine of a conformal
%           latitude, determine an ellipsoidal (geodetic) latitude. A
%           conformal latitude is a latitude in terms of a sphere that is
%           conformal with regard to the ellipsoid. That is, the
%           transformation from the ellipsoid to the sphere preserves
%           angles. This is such that the angle of intersection  between
%           any two lines on the ellipsoid is the same as the corresponding
%           angle on the sphere.
%
%INPUTS: sinChi The real sine of the conformal latitude. This can be a
%               scalar or a matrix of values to convert.
%        cosChi The real cosine of the conformal latitude. This can be a
%               scalar or a matrix of values to convert. This has to be the
%               same size as sinChi.
%             f The flattening factor of the reference ellipsoid. If this
%               argument is omitted, the value in Constants.WGS84Flattening
%               is used. This function is only accurate for values from 0
%               to about 0.999. Typically, this will be close to zero, as
%               in Constants.WGS84Flattening.
%
%OUTPUTS: phi The ellipsoidal latitudes in radians. This is the same size
%             as sinChi and cosChi.
%
%Conformal latitudes are discussed in Chapter 3 of [1]. This function
%implements the fixed-point iteration that is described in Section 2.9 of
%[2]. Smaller values of f lead to faster convergence. The maximum number of
%iterations in the function, which is typically never obtained, is set to
%be able to convert unusually large values of f up to about 0.999. A
%different fixed point iteration when given the conformal latitude itself
%is given in [3].
%
%REFERENCES:
%[1] J. P. Snyder, "Map projections- a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%[2] Office of Geomatics, "National geospatial-intelligence agency
%    standardization document: Implementation practice: The universal
%    grids and the transverse mercator and polar stereographic
%    map projections," 25 Mar. 2014. [Online]. Available:
%    http://earth-info.nga.mil/GandG/publications/NGA_SIG_0012_2_0_0_UTMUPS/NGA.SIG.0012_2.0.0_UTMUPS.pdf
%[3] Weisstein, Eric W. "Conformal Latitude." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/ConformalLatitude.html
%
%July 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(f))
    f=Constants.WGS84Flattening;
end

%The first numerical eccentricity of the ellipsoid.
e=sqrt(2*f-f^2);

N=numel(sinChi);
phi=zeros(size(sinChi));

for curPoint=1:N
    s=sinChi(curPoint);
    %It should converge in 7 or fewer iterations for
    %f=Constants.WGS84Flattening. The upper bound of 5e7 should be enough for
    %f=0.999.
    for curIter=1:3e7
        P=exp(e*atanh(e*s));
        term1=(1+sinChi(curPoint))*P*P;
        term2=(1-sinChi(curPoint));

        sNew=(term1-term2)/(term1+term2);

        if(abs(s-sNew)<=eps(s))
            break
        end
        s=sNew;
    end
    P=exp(e*atanh(e*s));
    c=(1/2)*((1+s)/P+(1-s)*P)*cosChi(curPoint);
    phi(curPoint)=atan2(s,c);
end

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:
function model = yalmip2scs(interfacedata)

%% From internal
data.A = -interfacedata.F_struc(:,2:end);
data.b = full(interfacedata.F_struc(:,1));
data.c =  interfacedata.c;
cones = [];
cones.z = interfacedata.K.f; %scs renamed the cone from f to z
cones.l = interfacedata.K.l;
cones.q = interfacedata.K.q;
cones.s = interfacedata.K.s;
cones.ep =interfacedata.K.e;
param = interfacedata.options.scs;
param.verbose = interfacedata.options.verbose;

%% Extract lower diagonal form for new SCS format
if ~isempty(cones.s) && any(cones.s)
    sdpA = data.A(1+cones.l + cones.z+sum(cones.q):end,:);
    sdpb = data.b(1+cones.l + cones.z+sum(cones.q):end,:);
    expA = data.A(end-3*cones.ep+1:end,:);
    expb = data.b(end-3*cones.ep+1:end,:);
    data.A = data.A(1:cones.l + cones.z+sum(cones.q),:);
    data.b = data.b(1:cones.l + cones.z+sum(cones.q),:);
    top = 1;
    for i = 1:length(cones.s)
        A = sdpA(top:top + cones.s(i)^2-1,:);
        b = sdpb(top:top + cones.s(i)^2-1,:);
        n = cones.s(i);
        ind = find(speye(n));
        b(ind) = b(ind)/sqrt(2);
        A(ind,:) = A(ind,:)/sqrt(2);
        ind = find(tril(ones(n)));
        A = A(ind,:);
        b = b(ind);
        data.A = [data.A;A];
        data.b = [data.b;b];
        top = top  + cones.s(i)^2;
    end
    data.A = [data.A;expA];
    data.b = [data.b;expb];
end

%% Collect in one structure
model.data = data;
model.cones = cones;
model.param = param;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Fast linear interpolation of ordered values (one dimendionnal)
%
% This matlab version is made to illustrate the behavior of the C version.
% Do not expect speed from this Matlab version  !
% Also, if this function is updated, its C version has to be updated
% accordingly.
%
% Note tha the C version makes absolutely NO check at all.
% Thus, it can crash very easily.
%
% Build the C version using
%  $ mex interp1ordered.c
%
% Input
%  x   ; [Nx1] Ordered reference abscissa
%  y   ; [Nx1] Ordered reference ordinates
%  xi  : [Mx1] Ordered abscissa where to interpolate
%  yid : [Mx1] Default out of bounds value
%
% Output
%  yi  : Interpolated output values
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>

function yi = interp1ordered(x, y, xi, yid)

    yi = yid*ones(size(xi));

    i=1;
    j=1;

    while xi(i)<x(1)
        yi(i) = yid;
        i=i+1;
    end

    while i<=length(xi)
        while j<length(x) && x(j+1)<xi(i);
            j=j+1;
        end
        if j==length(x); break; end

        g = (xi(i)-x(j))./(x(j+1)-x(j));

        yi(i) = y(j,:)*(1-g) + y(j+1,:)*g;

        i=i+1;
    end

    if i<=length(xi)
        yi(i:end) = yid;
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [fP] = callback_GetPercentile(fPercentile, vDistribution, fValue)
    % get percentile
    fP = abs(prctile(vDistribution, fPercentile) - fValue);
end

author: CelsoReyes


text:
function title = p02_title ( )

%*****************************************************************************80
%
%% P02_TITLE returns the title of problem p02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Product Peak';

  return
end


text:
%%%%% Author Muhammad Adil Farooqi %%%%
%%%%% E-mail: engr.adilfarooqi@gmail.com %%%%%
%%%%% OCR_ANN GUI %%%%%

%%
function varargout = OCR_ANN(varargin)
% OCR_ANN M-file for OCR_ANN.fig
%
%
%
% % Last Modified by GUIDE v2.5 28-May-2011 04:54:03

% Begin initialization code -
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @OCR_ANN_OpeningFcn, ...
                   'gui_OutputFcn',  @OCR_ANN_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code -

%%
% --- Executes just before OCR_ANN is made visible.
%%
function OCR_ANN_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to OCR_ANN (see VARARGIN)

% Choose default command line output for OCR_ANN
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes OCR_ANN wait for user response (see UIRESUME)
% uiwait(handles.figure1);

%%
% --- Outputs from this function are returned to the command line.
function varargout = OCR_ANN_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

%%
% --- Executes on button press in plot1_pushbutton.
function plot1_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to plot1_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%%%% step# 01 Image Acquisition %%%%%%

clc;
[filename, pathname] = uigetfile({'*.jpg';'*.bmp';'*.gif';'*.*'}, 'Pick an Image File');
TestIm = imread([pathname,filename]);

axes(handles.axes1);
imshow(TestIm);

save('TestIm','TestIm');
clear all;

%%
% --- Executes on button press in Axes2_pushbutton.
function Axes2_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Axes2_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%%% step# 02 Preprocessing %%%%
load TestIm

%Convert to gray scale
if size(TestIm,3)==3 %RGB image
    TestIm = rgb2gray(TestIm);
end
%Convert to binary image
threshold = graythresh(TestIm);
TestIm =~im2bw(TestIm,threshold);

%Remove Salt and paper noise
TestIm = medfilt2(TestIm);

axes(handles.axes2);
imshow(TestIm)

save('TestIm','TestIm');
clc; clear all;

%%
% --- Executes on button press in axes3_pushbutton.
function axes3_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to axes3_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%%%% step# 03 Segmentation %%%%%
load TestIm

%Label and Count Connected components
[L Ne] = bwlabel(TestIm);
glyphs = []  %%%Initilize glyphs matrix

   for n=1:Ne
    [r c] =find(L==n);

    glyph = TestIm(min(r):max(r),min(c):max(c));

    glyph = imresize(glyph,[35 35]);

    % Again convert to binaray image
    glyph = double(glyph);
    thresh = graythresh(glyph);
    glyph = im2bw(glyph,thresh);

    glyphs = [glyphs glyph]

    axes(handles.axes3);
    imshow(glyph)

   end
   [m n] = size(glyphs)
   set(handles.text6,'string',n/35)
   save('glyphs','glyphs');
   clear all;
%%
% --- Executes on button press in axes4_pushbutton.
function axes4_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to axes4_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%%%% step# 04 Classification %%%%

load glyphs;
load net;
[r c] = size(glyphs);

    glyphsclass = [];performance = [];%%initilize matrix
    for n = 0:35:(c-35)
        glyph = glyphs(1:35,1+n:(n+35));
        glyph = mat2vect(glyph);
        [glyphclass,PF,AF,E,Perf] = sim(net,glyph);

        performance = [performance Perf]
        axes(handles.axes4)
       plot(performance,'-.b*','LineWidth',2,'MarkerSize',10)
        glyphsclass = [glyphsclass glyphclass]

    end

save('glyphsclass','glyphsclass');
clear all;

%%
% --- Executes on button press in pushbutton11.
function pushbutton11_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton11 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%%% step# 04 Output %%%%%

load glyphsclass

[r c] = size(glyphsclass);
alphabets = [];space = ' ';haroof = [] %%%initilize matrix

for n=1:c

    glyphclass = glyphsclass(1:r,n) %%% individiual alphabet

    [i j] = max(glyphclass); %%% alphabet???

if j==5 | j==6 | j== 7 | j== 8 | j==9 | j==10 | j==11 | j==12 | j==13 | j==14 | j==15 | j==16 |...
        j==17 | j==18 | j==19 |j==20 | j==21 | j==22 | j==23 | j==24 | j==25

if j==5
    alphabet = '0627'
    harf='alif      '
end
if j==6
    alphabet = '0628'
    harf = 'bay       ';
end

if j==7
    alphabet = '062D'
    harf = 'bari he   ';
end

if j==8
    alphabet = '062F'
    harf = 'dal       ';
end

if j==9
    alphabet = '0631'
    harf = 're        ';
end

if j==10
    alphabet = '0633'
    harf = 'sin       ';
end

if j==11
    alphabet = '0635'
    harf = 'suad      ';
end

if j==12
    alphabet = '0637';
    harf = 'toe       ';
end

if j==13
    alphabet = '0639';
    harf = 'ain       ';
end

if j==14
    alphabet = '0641';
    harf = 'fe        ';
end
if j==15
    alphabet = '0642';
    harf = 'qaf       ';
end
if j==16
    alphabet = '06A9';
    harf = 'kaf       ';
end

if j==17
    alphabet = '0644';
    harf = 'lam       ';
end
if j==18
    alphabet ='0645';
    harf = 'mim       ';
end
if j==19
    alphabet = '0646';
    harf = 'nun       ';
end

if j==20
    alphabet = '0648';
    harf = 'vao       ';
end
if j==21
    alphabet = '06C1';
    harf = 'choti he  ';
end
if j==22
    alphabet = '06BE';
    harf = 'DoCashmiHe';
end
if j==23
    alphabet = '0621';
    harf = 'hamzah    ';
end
if j==24
    alphabet = '06CC';
    harf = 'choti ye  ';
end
if j==25
    alphabet = '06D2';
    harf = 'bari ye   ';
end

alphabets = [alphabets space alphabet]
haroof = [haroof space harf]
end

%%%If glyph is Diacritics then...
[x y] = size(alphabets);
[m n] = size(haroof);

    if j==1 & alphabet=='0628' %%%%bay
        alphabets(1,(y-3):y) ='0628'%%%each alphabet is of length 4
        haroof(1,(n-9):n) = 'bay       '%%%each harf is of length 10
    end
    if j==2 & alphabet=='0628'
        alphabets(1,(y-3):y) = '062A'%%te
        haroof(1,(n-9):n) = 'te        '
    end
    if j==1 & alphabet== '062A'%%te
        alphabets(1,(y-3):y) = '062B'%%se
        haroof(1,(n-9):n) = 'se        '
    end
    if j==3 & alphabet =='0628'%%bay
        alphabets(1,(y-3):y) = '0679'%%te
        haroof(1,(n-9):n) = 'te        '

    end
    if j==1 & alphabet == '062D' %%%%bar? he
        alphabets(1,(y-3):y) = '062C'%%j?m
        haroof(1,(n-9):n) = 'jim       '
    end
    if j==2 & alphabet =='062D'
        alphabets(1,(y-3):y) = '0686'%%chay
        haroof(1,(n-9):n) = 'chay      '
    end
    if j==1 & alphabet =='062F'  %%%d?l
        alphabets(1,(y-3):y) = '0630'%%z?l
        haroof(1,(n-9):n) = 'zal       '
     end
    if j==3 & alphabet=='062F'
        alphabets(1,(y-3):y) = '0688'%%dd?l
        haroof(1,(n-9):n) = 'ddal      '
     end
    if j == 1 & alphabet =='0631' %%%raa
        alphabets(1,(y-3):y) = '0632'%%ze
        haroof(1,(n-9):n) = 'ze        '
     end
    if j==2 & alphabet =='0631'
        alphabets(1,(y-3):y) = '0698'%%zhe
        haroof(1,(n-9):n) = 'zhe       '
    end
    if j== 3 & alphabet =='0631'
        alphabets(1,(y-3):y) ='0691'%%%rre
        haroof(1,(n-9):n) = 'rre       '
     end
    if j==1 & alphabet == '0635'%%%soad
        alphabets(1,(y-3):y) = '0636'%%%zu'?d
        haroof(1,(n-9):n) = 'zuad      '
     end
    if j==2 & alphabet== '0633'%%%seen
        alphabets(1,(y-3):y) = '0634'%%%sh?n
        haroof(1,(n-9):n) = 'shin      '
     end
    if j==1 & alphabet== '0637'%%%toa
        alphabets(1,(y-3):y)='0638'%%%zo'e
        haroof(1,(n-9):n) = 'zoe       '

    end
     if j ==1 & alphabet== '0639'%%%aaen
        alphabets(1,(y-3):y) = '063A'%%ghain
        haroof(1,(n-9):n) = 'ghain     '
     end
     if j ==4 & alphabet== '06A9'%%%kaaf
        alphabets(1,(y-3):y) = '06AF'%%g?f
        haroof(1,(n-9):n) = 'gaf       '
     end

end
[m n] = size(alphabets)
set(handles.text5,'string',n/5)
set(handles.text4,'string',haroof)

save('alphabets','alphabets')

%%
%--- Executes on button press in pushbutton12.
function pushbutton12_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton12 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%%%% step# 05 Export Alphabets %%%
load alphabets

fid = fopen('text.RTF','w');

fprintf(fid,alphabets);
fclose(fid)

    winopen('text.RTF')

function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

%%
% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double

%%
% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit5_Callback(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit5 as text
%        str2double(get(hObject,'String')) returns contents of edit5 as a double

%%
% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%%
% --- Executes on selection change in listbox1.
function listbox1_Callback(hObject, eventdata, handles)
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns listbox1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from listbox1

%%
% --- Executes during object creation, after setting all properties.
function listbox1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Sable

repo: mcbench-benchmarks


text:
clear all;
% path to db files
addpath('/media/carsen/DATA1/2P/dbfiles/')
compile_dbs;

% path to saved F files (and root directory)
ops0.RootDir                = '//zserver.cortexlab.net/Data/Subjects/';
ops0.ResultsSavePath        = '/media/carsen/DATA1/2P/F/';

% choose experiment
iexp = 31;
db = db0(iexp);

% build ops
ops = build_ops3(db, ops0);
root = ops.ResultsSavePath;

% load mean image and stats for each plane in recording
stat = [];
xL = [];
yL = [];
clear mimg;
for iplane = 1:ops.nplanes
    fname  = sprintf('F_%s_%s_plane%d.mat', db.mouse_name, db.date, iplane);
    dat = load(fullfile(root, fname));

    mimg{iplane} = dat.ops.mimg1(dat.ops.yrange, dat.ops.xrange);

    stat = cat(2, stat, dat.stat);
    yL(iplane) = numel(dat.ops.yrange);
    xL(iplane) = numel(dat.ops.xrange);
end

% compute outlines of cells
clear img;
ij = 1;
clf;
for iplane = 1:ops.nplanes
    mimg{iplane} = mimg{iplane} - min(mimg{iplane}(:));
    mimg{iplane} = mimg{iplane} / max(mimg{iplane}(:));

    % background of mean image
    img{iplane} = repmat(mimg{iplane}*2,1,1,3);
    img{iplane} = rgb2hsv(img{iplane});
    img{iplane} = reshape(img{iplane},[],3);

    colormap('gray');
    hold all;
    while stat(ij).iplane == iplane
        if stat(ij).iscell
            % find pixels that are exterior to the cell
            idist  = sqrt(bsxfun(@minus, stat(ij).xpix', stat(ij).xpix).^2 + ...
                bsxfun(@minus, stat(ij).ypix', stat(ij).ypix).^2);
            idist  = idist - diag(NaN*diag(idist));
            extpix = sum(idist <= sqrt(2)) <= 6;
            xext = stat(ij).xpix(extpix);
            yext = stat(ij).ypix(extpix);

            % color the exterior pixels
            ipix = sub2ind([yL(iplane) xL(iplane)], yext, xext);
            img{iplane}(ipix, 1) = rand;
            img{iplane}(ipix, 2) = 1;
            img{iplane}(ipix, 3) = 1;

        end

        ij = ij+1;

        if ij > numel(stat)
            break;
        end
    end

end

%%
clf;
set(gcf,'color','w');

% plot of one plane with circled ROIs
iplane = 5;

imgj = img{iplane};
imagesc(hsv2rgb(reshape(imgj,yL(iplane),xL(iplane),3)))
axis off;
axis square;

author: cortex-lab


text:
% lapplot() -  Compute the discrete laplacian of EEG scalp distribution(s)
%
% Usage:
%   >> laplace = lapplot(map,eloc_file,draw)
%
% Inputs:
%    map        - Activity levels, size (nelectrodes,nmaps)
%    eloc_file	- Electrode location filename (.loc file)
%                 For format, see  >> topoplot example
%    draw       - If defined, draw the map(s) {default: no}
%
% Output:
%    laplace    - Laplacian map, size (nelectrodes,nmaps)
%
% Note: uses del2()
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1998
%
% See also: topoplot(), gradplot()

% Copyright (C) Scott Makeig, SCCN/INC/UCSD, La Jolla, 1998
%
% 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

% 01-25-02 reformated help & license, added links -ad

function [laplac] = lapplot(map,filename,draw)

if nargin < 2
	help lapplot;
	return;
end;

MAXCHANS = size(map,1);
GRID_SCALE = 2*MAXCHANS+5;
MAX_RADIUS = 0.5;

% ---------------------
% Read the channel file
% ---------------------
if isstr( filename )
	fid = fopen(filename);
	locations = fscanf(fid,'%d %f %f %s',[7 MAXCHANS]);
	fclose(fid);
	locations = locations';
	Th = pi/180*locations(:,2);   % convert degrees to rads
	Rd = locations(:,3);
	ii = find(Rd <= MAX_RADIUS); % interpolate on-scalp channels only
	Th = Th(ii);
	Rd = Rd(ii);
	[x,y] = pol2cart(Th,Rd);
else
	x = real(filename);
	y = imag(filename);
end;

% ---------------------------------------------------
% Locate nearest position of an electrode in the grid
% ---------------------------------------------------
xi = linspace(-0.5,0.5,GRID_SCALE);   % x-axis description (row vector)
yi = linspace(-0.5,0.5,GRID_SCALE);   % y-axis description (row vector)
for i=1:MAXCHANS
   [useless_var horizidx(i)] = min(abs(y(i) - xi));    % find pointers to electrode
   [useless_var vertidx(i)] = min(abs(x(i) - yi));     % positions in Zi
end;

% -----------------
% Compute laplacian
% -----------------
for i=1:size(map,2)
   	[Xi,Yi,Zi] = griddata(y,x,map(:,i),yi',xi, 'v4');   % interpolate data

   	laplac2D = del2(Zi);
	positions = horizidx + (vertidx-1)*GRID_SCALE;
	laplac(:,i) = laplac2D(positions(:));

	% ------------------
	% Draw laplacian map
	% ------------------
	if exist('draw');
        mask = (sqrt(Xi.^2+Yi.^2) <= MAX_RADIUS);
        laplac2D(find(mask==0)) = NaN;

		subplot(ceil(sqrt(size(map,2))), ceil(sqrt(size(map,2))), i);
		contour(laplac2D);
		title( int2str(i) );

% %%% Draw Head %%%%
ax = axis;
width = ax(2)-ax(1);
axis([ax(1)-width/3 ax(2)+width/3 ax(3)-width/3 ax(4)+width/3])
steps = 0:2*pi/100:2*pi;
basex = .18*MAX_RADIUS;
tip = MAX_RADIUS*1.15;
base = MAX_RADIUS-.004;
EarX = [.497 .510 .518 .5299 .5419 .54 .547 .532 .510 .489];
EarY = [.0555 .0775 .0783 .0746 .0555 -.0055 -.0932 -.1313 -.1384 -.1199];

HCOLOR = 'k';
HLINEWIDTH = 1.8;

% Plot Head, Ears, Nose
hold on
plot(1+width/2+cos(steps).*MAX_RADIUS*width,...
     1+width/2+sin(steps).*MAX_RADIUS*width,...
    'color',HCOLOR,'Linestyle','-','LineWidth',HLINEWIDTH); % head

plot(1+width/2+[.18*MAX_RADIUS*width;0;-.18*MAX_RADIUS*width],...
     1+width/2+[base;tip;base]*width,...
    'Color',HCOLOR,'LineWidth',HLINEWIDTH);                 % nose

plot(1+width/2+EarX*width,...
     1+width/2+EarY*width,...
           'color',HCOLOR,'LineWidth',HLINEWIDTH)           % l ear
plot(1+width/2-EarX*width,...
     1+width/2+EarY*width,...
           'color',HCOLOR,'LineWidth',HLINEWIDTH)           % r ear

hold off
axis off

	end;
end;

return;

author: PatternRecognition

repo: OpenBMI


text:
function legendre_set_test ( )

%*****************************************************************************80
%
%% LEGENDRE_SET_TEST tests LEGENDRE_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_SET_TEST\n' );
  fprintf ( 1, '  LEGENDRE_SET returns points and weights of \n' );
  fprintf ( 1, '  Gauss-Legendre quadrature rules.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N               1             X^4           Runge\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 10
    [ x, w ] = legendre_set ( n );
    e1 = sum ( w(1:n,1) );
    e2 = w(1:n,1)' * x(1:n,1).^4;
    e3 = w(1:n,1)' * ( 1.0 ./ ( 1.0 + 25.0 * x(1:n,1).^2 ) );
    fprintf ( 1, '  %2d  %14.6g  %14.6g  %14.6g\n', n, e1, e2, e3 );
  end

  return
end


text:
function pf = loadPoleFigure_scintag(fname,varargin)
% import data fom scintag ascii file
%
% Syntax
%   pf = loadPoleFigure_scintag(fname)
%
% Input
%  fname    - filename
%
% Output
%  pf - vector of @PoleFigure
%
% See also
% ImportPoleFigureData loadPoleFigure

fid = fopen(fname,'r');
d = fread(fid);
fclose(fid);

try
  d(d==13) = [];
  d = char(d)';

  file_id = regexpi(d(1:200),'(?<=File ID:(.*?)\n)(.*?)\n','match');
  comment = file_id{1};
  h = string2Miller(comment);

  ranges = regexpsplit(d,'Range');
  axis = regexpsplit(ranges{1},'Axis');

  % extract the stepsize
  f = regexpi(axis(2:end),'(?<=Step Size  :)(.*?)\n','match');
  f = [f{:}]; f = [f{:}];
  f = str2num(f);
  step = f(find(f>0,1,'first'));

  % process the ranges
  for k=2:numel(ranges)
    range = ranges{k};

    po = sscanf(range,'%d')*step;

    low = sscanf(range(41:49),'%f');
    high = sscanf(range(87:95),'%f');

    bg = (low+high)/2;

    t = find(double(range(97:end))==10,1,'first');
    f = numel(sscanf(range(97:97+t),'%f'));

    data = sscanf(range(97:end),'%f');

    if f==3 % counts per second
      az = data(1:3:end);
      I(:,k-1) = data(2:3:end);
    else    % counts in time
      az = data(1:4:end);
      I(:,k-1) = data(2:4:end)*1000./data(3:4:end);
    end

    % background correction
    I(:,k-1) = I(:,k-1) - bg;

    % specimen direction
    r(:,k-1) = vector3d('theta',po*degree,'rho',az*degree);

  end

  pf = PoleFigure( h,r,I,varargin{:});
catch
  interfaceError(fname);
end

author: mtex-toolbox

repo: mtex

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


text:
% Setting An Initial Pressure Gradient Example
%
% This example demonstrates how to set an initial pressure gradient using
% kspaceSecondOrder. It builds on the Comparison Of Modelling Functions
% Example.
%
% author: Bradley Treeby
% date: 28th October 2010
% last update: 17th October 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 64;            % number of grid points in the x (row) direction
Ny = 64;            % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
k = kgrid.k;

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.alpha_power = 1.5;   % [dB/(MHz^y cm)]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]

% set plot_frames to true to produce the plots given in the documentation
% set plot_frames to false to visualise the propagation of the pressure
% field
plot_frames = false;
if plot_frames
    % calculate and plot the pressure at a particular value of t
    kgrid.t_array = [0 1000]*1e-9;
else
    % calculate the plot the progression of the pressure field with t
    dt = 5e-9;
    t_end = 5e-6;
    kgrid.t_array = 0:dt:t_end;
end

% create source distribution using makeDisc
disc_magnitude = 4; % [Pa]
disc_x_pos = 25;    % [grid points]
disc_y_pos = 40;    % [grid points]
disc_radius = 4;    % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; % [Pa]
disc_x_pos = 40;    % [grid points]
disc_y_pos = 25;    % [grid points]
disc_radius = 3;    % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source_distribution = disc_1 + disc_2;

% define a single sensor point
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/2, Ny/2) = 1;

% define the input arguments
input_args = {'PlotFrames', plot_frames, 'MeshPlot', true, 'PlotScale', [0 3], 'ExpandGrid', true};

% assign the source distribution to the initial pressure
source.p0 = source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});

% assign the source distribution to the initial pressure gradient
source = rmfield(source, 'p0');
source.dp0dt = 5e6*source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});

% assign the source distribution to both the initial pressure and the
% initial pressure gradient
source.p0 = source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});

repo: Aero-matlab


text:
function mat = eulerAnglesToRotation3d(phi, theta, psi, varargin)
%EULERANGLESTOROTATION3D Convert 3D Euler angles to 3D rotation matrix.
%
%   MAT = eulerAnglesToRotation3d(PHI, THETA, PSI)
%   Creates a rotation matrix from the 3 euler angles PHI THETA and PSI,
%   given in degrees, using the 'XYZ' convention (local basis), or the
%   'ZYX' convention (global basis). The result MAT is a 4-by-4 rotation
%   matrix in homogeneous coordinates.
%
%   PHI:    rotation angle around Z-axis, in degrees, corresponding to the
%       'Yaw'. PHI is between -180 and +180.
%   THETA:  rotation angle around Y-axis, in degrees, corresponding to the
%       'Pitch'. THETA is between -90 and +90.
%   PSI:    rotation angle around X-axis, in degrees, corresponding to the
%       'Roll'. PSI is between -180 and +180.
%   These angles correspond to the "Yaw-Pitch-Roll" convention, also known
%   as "Tait–Bryan angles".
%
%   The resulting rotation is equivalent to a rotation around X-axis by an
%   angle PSI, followed by a rotation around the Y-axis by an angle THETA,
%   followed by a rotation around the Z-axis by an angle PHI.
%   That is:
%       ROT = Rz * Ry * Rx;
%
%   MAT = eulerAnglesToRotation3d(ANGLES)
%   Concatenates all angles in a single 1-by-3 array.
%
%   ... = eulerAnglesToRotation3d(ANGLES, CONVENTION)
%   CONVENTION specifies the axis rotation sequence.
%   Supported conventions are: 'ZYX', 'ZYZ'. Default is 'ZYX'.
%
%   Example
%   [n e f] = createCube;
%   phi     = 20;
%   theta   = 30;
%   psi     = 10;
%   rot = eulerAnglesToRotation3d(phi, theta, psi);
%   n2 = transformPoint3d(n, rot);
%   drawPolyhedron(n2, f);
%
%   See also
%   transforms3d, createRotationOx, createRotationOy, createRotationOz
%   rotation3dAxisAndAngle
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-07-22,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY
%   2011-06-20 rename and use degrees

p = inputParser;
validStrings = {'ZYX','ZYZ'};
addOptional(p,'convention','ZYX',@(x) any(validatestring(x,validStrings)));
parse(p,varargin{:});
convention=p.Results.convention;

% Process input arguments
if size(phi, 2) == 3
    % manages arguments given as one array
    psi     = phi(:, 3);
    theta   = phi(:, 2);
    phi     = phi(:, 1);
end

% create individual rotation matrices
k = pi / 180;

switch convention
    case 'ZYX'
        rot1 = createRotationOx(psi * k);
        rot2 = createRotationOy(theta * k);
        rot3 = createRotationOz(phi * k);
    case 'ZYZ'
        rot1 = createRotationOz(psi * k);
        rot2 = createRotationOy(theta * k);
        rot3 = createRotationOz(phi * k);
end

% concatenate matrices
mat = rot3 * rot2 * rot1;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [brs,delay] = ssr_brs_nfchoa(X,phi,xs,src,irs,conf)
%SSR_BRS_NFCHOA binaural room scanning (BRS) set for use with the SoundScape
%Renderer
%
%   Usage: [brs,delay] = ssr_brs_nfchoa(X,phi,xs,src,irs,conf)
%
%   Input parameters:
%       X       - listener position / m
%       phi     - azimuthal head orientation / rad
%       xs      - virtual source position [ys > Y0 => focused source] / m
%       src     - source type: 'pw' - plane wave
%                              'ps' - point source
%                              'fs' - focused source
%       irs     - impulse response data set for the second sources
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       brs     - conf.N x 2*nangles matrix containing all impulse responses (2
%                 channels) for every angles of the BRS set
%       delay    - delay added by driving function / s
%
%   SSR_BRS_NFCHOA(X,phi,xs,src,irs,conf) prepares a BRS set for a virtual
%   source at position xs for a virtual loudspeaker array driven by
%   nearfield compensated higher order Ambisonics (NFC-HOA) and the given
%   listener position. One way to use this BRS set is using the
%   SoundScapeRenderer (SSR), see http://spatialaudio.net/ssr/
%
%   See also: ir_generic, ir_nfchoa, driving_function_imp_nfchoa

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 6;
nargmax = 6;
narginchk(nargmin,nargmax);
isargposition(X);
isargxs(xs);
isargscalar(phi);
isargstruct(conf);

%% ===== Computation =====================================================
% Secondary sources
x0 = secondary_source_positions(conf);
% Calculate driving function
[d,~,delay] = driving_function_imp_nfchoa(x0,xs,src,conf);
% Calculate brs set
brs = ssr_brs(X,phi,x0,d,irs,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function [Z,Q1,R] = sparsenull(A)

[Q,R] = qr(A');
n = max(find(sum(abs(R),2)));
Q1 = Q(:,1:n);
R = R(1:n,:);
Z = Q(:,n+1:end); % New basis


text:
function unique_num = point_tol_unique_count ( m, n, a, tol )

%*****************************************************************************80
%
%% POINT_TOL_UNIQUE_COUNT counts the tolerably unique points.
%
%  Discussion:
%
%    The input data is an M x N array A, representing the M-dimensional
%    coordinates of N points.
%
%    This function uses a simple but expensive approach.  The first point
%    is accepted as unique.  Each subsequent point is accepted as unique
%    only if it is at least a tolerance away from all accepted unique points.
%    This means the expected amount of work is O(N^2).
%
%    The output is the number of unique points in the list.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 July 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows.
%
%    Input, integer N, the number of columns.
%
%    Input, real A(M,N), the array of N columns of data.
%
%    Input, real TOL, a tolerance.
%
%    Output, integer UNIQUE_NUM, the number of unique points.
%
  unique(1:n) = 1;
  unique_num = n;

  for i = 2 : n

    for j = 1 : i - 1
      if ( unique(j) )
        dist = sqrt ( sum ( ( a(1:m,i) - a(1:m,j) ).^2 ) );
        if ( dist <= tol )
          unique(i) = 0;
          unique_num = unique_num - 1;
          break
        end
      end
    end

  end

  return
end


text:
function [y] = cellshift(x, shift, dim, maxshift, flag)

% CELLSHIFT(X, SHIFT, DIM)

if nargin<5 || isempty(flag)
  flag = 'overlap'; % can also be 'zero' or 'nan'
  % overlap: only return the data samples that fully overlap for all shifts
  % zero: return 0's for the samples at which no data is present
  % nan: return nans for the samples at which no data is present
end

if nargin<4 || isempty(maxshift)
  maxshift = max(abs(shift));
else
  if any(abs(shift(:)))>abs(maxshift)
    error('the value for maxshift should be >= shift');
  end
end

if nargin<3 || isempty(dim)
  dim = find(size(x{1})>1, 1, 'first');
end

siz = size(shift);
if siz(1)==1 && siz(2)==1
  % this is OK, continue below
elseif siz(1)==1 && siz(2)>1
  % call the function recursively, and concatenate, same set of shifts for
  % all channels
  x = x(:)';
  y = cell(numel(shift),numel(x));
  for k = 1:numel(shift)
    %y(k,:) = cellshift(x, shift(k), dim, max(abs(shift)));
    y(k,:) = cellshift(x, shift(k), dim, [abs(min(shift)) abs(max(shift))], flag);
  end
  for k = 1:size(y,2)
    y{1,k} = cat(1,y{:,k});
  end
  y = y(1,:);
  return;
elseif siz(1)>1 && siz(2)>=1
  % different shifts, call the function recursively, and concatenate.
  % note that the order of the rows in the output is different from when
  % the same shift applies to all channels
  edges =[abs(nanmin(shift(:))) abs(nanmax(shift(:)))];
  y     = cell(siz(1),numel(x));
  for k = 1:siz(1)
    tmpshift = shift(k,:);
    tmpshift(~isfinite(tmpshift)) = [];
    if isempty(tmpshift)
      tmpshift = 0;
    end
    y(k,:) = cellshift(cellrowselect(x, k), tmpshift, dim, edges, flag);
  end
  for k = 1:size(y,2)
    y{1,k} = cat(1,y{:,k});
  end
  y = y(1,:);
  return;
end

%-----------------------------------------------------
maxshift = abs(maxshift);
if numel(maxshift)==1, maxshift = maxshift([1 1]); end

nx = size(x);
if ~iscell(x) || length(nx)>2 || all(nx>1)
  error('incorrect input for cellshift');
end

n    = numel(x);
nsmp = cellfun('size', x, dim);

switch flag
  case 'overlap'
    beg1 = ones(1,n) + shift + maxshift(1);
    end1 = nsmp      + shift - maxshift(2);

    beg2 = ones(size(beg1));
    end2 = end1-beg1+1;

    nsmp = nsmp - sum(maxshift);

  case {'zero' 'nan'}
    beg1 = ones(1,n);
    end1 = nsmp;

    beg2 = ones(1,n) - shift + maxshift(2);
    end2 = nsmp      - shift + maxshift(2);

    nsmp = nsmp + sum(maxshift);
end

usenans = strcmp(flag, 'nan');

switch dim
case 1
  y   = cell(size(x));
  for k = 1:n
    tmp = zeros(nsmp(k), size(x{k},2));
    if usenans, tmp(:) = nan; end
    tmp(beg2(k):end2(k),:) = x{k}(beg1(k):end1(k),:);
    y{k} = tmp;
  end
case 2
  y   = cell(size(x));
  for k = 1:n
    tmp = zeros(size(x{k},1), nsmp(k));
    if usenans, tmp(:) = nan; end
    tmp(:, beg2(k):end2(k)) = x{k}(:,beg1(k):end1(k));
    y{k} = tmp;
  end
otherwise
  error('dimensionality of >2 is not supported');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%GUI script for reading ctf data into eeglab.
%to run type ctf2eeglab on the command line.

clear all;
folder = uigetdir('*.ds','Data Set Finder');
[ctf2] = ctf_read_res4(folder,1);
sensors = {'meg','eeg','ref','vc'};
sensloc = 'ctf2.sensor.index';
sensorlen = '1:';
num2str(length(ctf2.sensor.index.meg_sens));

k = menu('Choose the channels you would like to use:','MEG','EEG','Reference','Virtual');
sensloc = strcat(sensloc,sensors(k));
sensorlen = strcat(sensorlen, num2str(length(eval(sensloc{1}))));
prompt   = {'Enter Sensor Numbers to Read (ie 2:7 or 1 3 7 23 65)'};
title    = 'Input for Sensor Range (Default is all sensors)';
lines = 1;
def     = {sensorlen};
sensnum   = inputdlg(prompt,title,lines,def);
for i = str2num(sensnum{1});
    if ~ismember(str2num(def{1}),i)
        errordlg('There is no such sensor','Sensor Range Error')
        exit(1);
    end
end

trialnums = '1:';
trialnums = strcat(trialnums, num2str(ctf2.setup.number_trials));
prompt2 = {'Enter the trials you would like to use(ie 1:3 or 2 4 5):'};
def2 = {trialnums};
trials = inputdlg(prompt2, 'Input for trial range', 1,def2);
for i = str2num(trials{1})
    if ~ismember(str2num(def2{1}),i)
        errordlg('Trial not in dataset','Trial Error')
        exit(1);
    end
end
button = questdlg('Would you like to use markers?','Use Markers?','Yes','No','Yes');
if strcmp(button,'Yes')
    [marker_info] = readmarkerfile(folder);
    markers = marker_info.marker_names;
    m = menu('Please choose which marker you would like to use', markers);
end
prompt   = {'Enter the start time (sec):', 'Enter the end time:'};
title    = 'Time Window';
lines = 1;
def     = {'0', '.15'};
def(1) = {num2str(ctf2.setup.start_sec)};
def(2) = {num2str(ctf2.setup.end_sec)};
markeranswer = inputdlg(prompt,title,lines,def);
wind = [0 0];
wind(1) = str2num(markeranswer{1});
wind(2) = str2num(markeranswer{2});

if strcmp(button,'Yes')
    [data]=readepochs(folder,'marker_info',marker_info,'ctf',ctf2,sensors{k}, str2num(sensnum{1}),'trials',str2num(trials{1}), 'markers',markers(m),'window', wind);
else
    [data]=readepochs(folder,'ctf',ctf2,sensors{k}, str2num(sensnum{1}),'trials',str2num(trials{1}), 'window', wind);
end

%wind = [-.1 .15];
%[data]=readepochs(folder,'markers',{'click'},'window',wind);
%[data]=readepochs(folder,'window',wind,'trials',[1:10],'megsens',[2:7]);

dat=permute(data.epochs{1},[2,1,3]);
size(dat)
[ALLEEG EEG CURRENTSET ALLCOM] = eeglab;

EEG = pop_importdata( 'nbchan', size(dat,1), 'dataformat', 'array', 'data', 'dat', 'pnts', size(dat,2), 'srate', data.setup.sample_rate, 'xmin', wind(1));
[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET,  'setname',data.setup.subject);
eeglab redraw;

for i=1:size(dat,1)
    if ~isempty(data.sensor.location)
        EEG.chanlocs(i).labels=char(data.sensor.label(i));
        EEG.chanlocs(i).X=data.sensor.location(1,i);
        EEG.chanlocs(i).Y=data.sensor.location(2,i);
        EEG.chanlocs(i).Z=data.sensor.location(3,i);
    end
end
EEG.chanlocs=pop_chanedit(EEG.chanlocs,  'convert', 'cart2topo', 'convert', 'cart2sph');
[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
eeglab redraw;
clear dat data ctf2 marker_info markers sensors;


text:
function setup(subnet, snum, enum, setupfile, chanmatch)
debug.printfunctionstack('>');
    [paths,PARAMS,subnets]=pf2PARAMS(setupfile);
    if ~exist('enum','var')
        enum = snum;
    end
    subnets(1).sites = get_closest_sites(subnets(1).source.longitude, subnets(1).source.latitude, subnets(1).radius, paths.DBMASTER, Inf, snum, enum+1, chanmatch);
    sites = subnets(1).sites;
	show_sites(sites);
    save sites.mat sites

%     %gismo_datasource = datasource('antelope', dbwfdata);
%     % Add response data (calib only?) for each scnl
%     for k=1:numel(subnets(1).sites)
%         %try
%                 subnets(1).sites(k).response = response_get_from_db(subnets(1).sites(k).channeltag.station, subnets(1).sites(k).channeltag.channel, snum, PARAMS.f, paths.DBMASTER);
%                 % for Sakurajima I get:
%                     % Error using response_get_from_db (line 59)
%                     % Database does not contain a sensor table.
%         %catch
%         %        subnets(1).sites(k).response.calib = NaN;
%         %end
%     end
    %save2mat(sprintf('pf/%s_%s.mat',subnet,datestr(snum,'yyyymmdd')), subnets, paths, PARAMS)

    % Replace this whole complicated response thing above with a simple
    % grab for calib value. Not much good for short periods, but who cares?
    % - implemented in closest stations

    save2mat(sprintf('pf/%s.mat',subnet), subnets, paths, PARAMS)
debug.printfunctionstack('<');
end

function [paths,PARAMS,subnets]=pf2PARAMS(setupfile)
    debug.printfunctionstack('>')

    % create pointer to main parameter file
    [dirname, filename, ext] = fileparts(setupfile);
    if exist(setupfile, 'file')
        setuppf = dbpf(sprintf('%s/%s',dirname,filename));

        % subnets
        subnet_tbl = pfget_tbl(setuppf, 'subnets');
        for c=1:numel(subnet_tbl)
            fields = regexp(subnet_tbl{c}, '\s+', 'split');
            subnets(c).name = fields{1};
            subnets(c).source.latitude = str2double(fields{2});
            subnets(c).source.longitude = str2double(fields{3});
            subnets(c).radius = str2double(fields{4});
            subnets(c).use = str2double(fields{5});
        end

        % Maximum number of scnls to display in a spectrogram
        PARAMS.max_number_scnls = pfget_num(setuppf, 'max_number_scnls');

        % Select channels to use according to this channel mask
        PARAMS.channel_mask = pfget(setuppf, 'channel_mask');

        % paths (removed from setup.pf file 2013/04/22)
        paths.DBMASTER = pfget(setuppf, 'dbmaster');
        paths.PFS = 'pf';
        paths.spectrogram_plots = 'spectrograms';

        % datasource
        try
            datasources = pfget_tbl(setuppf, 'datasources');
            for c=1:numel(datasources)
                fields = regexp(datasources{c}, '\s+', 'split');
                PARAMS.datasource(c).type = fields{1};
                PARAMS.datasource(c).path = fields{2};
                if numel(fields)>2
                    PARAMS.datasource(c).port = fields{3};
                end
            end
        end

        % archive_datasource
        PARAMS.switch_to_archive_after_days = pfget(setuppf, 'switch_to_archive_after_days');
        try
            archive_datasources = pfget_tbl(setuppf, 'archive_datasources');
            for c=1:numel(archive_datasources)
                fields = regexp(archive_datasources{c}, '\s+', 'split');
                PARAMS.archive_datasource(c).type = fields{1};
                PARAMS.archive_datasource(c).path = fields{2};
                if numel(fields)>2
                    PARAMS.archive_datasource(c).port = fields{3};
                end
            end
        end

        % waveform processing
        lowcut	 = pfget_num(setuppf,'lowcut');
        highcut	 = pfget_num(setuppf,'highcut');
        npoles	 = pfget_num(setuppf,'npoles');
        PARAMS.filterObj = filterobject('b',[lowcut highcut],npoles);

        % Spectrograms
        PARAMS.spectralobject = spectralobject( ...
            pfget_num(setuppf,'nfft'), ...
            pfget_num(setuppf,'overlap'), ...
            pfget_num(setuppf,'max_freq'), ...
            [ pfget_num(setuppf,'blue') pfget_num(setuppf,'red')] ...
        );

        % Derived data
        PARAMS.surfaceWaveSpeed = pfget_num(setuppf,'surfaceWaveSpeed');
        PARAMS.df = pfget_num(setuppf, 'df');
        PARAMS.f = 0:PARAMS.df:50;

        % Alarm system
        PARAMS.triggersForAlarmFraction = pfget_num(setuppf,'triggersForAlarmFraction');

        debug.print_debug(1, 'PARAMS setup OK')
    else
        error(sprintf('%s: parameter file %s.pf does not exist',mfilename, setupfile));
    end
    debug.printfunctionstack('<')
end

function save2mat(RUNTIMEMATFILE, subnets, paths, PARAMS)
    % write pf/tremor_runtime.mat, preserve current version if it already exists
    if exist(RUNTIMEMATFILE, 'file')
        system(sprintf('mv %s %s.%s',RUNTIMEMATFILE,RUNTIMEMATFILE,datestr(now,30)));
    end

    save(RUNTIMEMATFILE, 'subnets', 'paths', 'PARAMS');
end

author: geoscience-community-codes

repo: GISMO


text:
%%*****************************************************************
%% HKMcorr: corrector step for the HKM direction.
%%
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [dX,dy,dZ,resnrm,EinvRc] = HKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
            dX,dZ,coeff,L,X,Z);

    global matfct_options solve_ok

    printlevel = par.printlevel;
%%
    [rhs,EinvRc] = HKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ);
    m = length(rp); ncolU = size(coeff.mat12,2);
    rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
%%
    solve_ok = 1; resnrm = norm(rhs);
    if strcmp(matfct_options,'chol') | strcmp(matfct_options,'spchol') ...
       | strcmp(matfct_options,'ldl') | strcmp(matfct_options,'spldl')
       [xx,resnrm,solve_ok] = symqmr(coeff,rhs,L,[],[],printlevel);
       if (solve_ok<=0) & (printlevel)
          fprintf('\n  warning: symqmr fails: %3.1f.',solve_ok);
       end
    else
       [xx,resnrm,solve_ok] = mybicgstab(coeff,rhs,L,[],[],printlevel);
       if (solve_ok<=0) & (printlevel)
          fprintf('\n  warning: bicgstab fails: %3.1f.',solve_ok);
       end
    end
    if (printlevel>=3); fprintf('%2.0d ',length(resnrm)-1); end
%%
    [dX,dy,dZ] = HKMdirfun(blk,At,par,Rd,EinvRc,X,xx,m);
%%*****************************************************************

author: yu-jiang

repo: radpbook

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


text:
function ps = normpdf2( xs, m, C )
% Normal prob. density function (pdf) with arbitrary covariance matrix.
%
% Evaluate the multi-variate density with mean vector m and covariance
% matrix C for the input vector xs.  Assumes that the N datapoints are d
% dimensional.  Then m is dx1 or 1xd, C is dxd, and xs is dxN or NxD where
% N is the number of samples to be evaluated.
%
% USAGE
%  ps = normpdf2( xs, m, C )
%
% INPUTS
%  xs  - points to evaluated (Nxd or dxN)
%  m   - mean vector (dx1 or 1xd)
%  C   - Covariance matrix (dxd)
%
% OUTPUTS
%  ps  - probability density at each x (Nx1)
%
% EXAMPLE
%  ps = normpdf2( randn(10,2), [0 0], eye(2) )
%
% See also NORMPDF
%
% Piotr's Image&Video Toolbox      Version 2.30
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get dimensions of data
d=length(m);
if( size(xs,1)~=d ); xs=xs'; end
N=size(xs,2);

if( d==1 ) % fast special case
  ps = 1/sqrt(2*pi*C) * exp(-(xs-m).*(xs-m)/(2*C))';

elseif( rcond(C)<eps ) % if matrix is badly conditioned
  warning('normpdf2: Covariance matrix close to singular.'); %#ok<WNTAG>
  ps = zeros(N,1);

else % get probabilities
  xs = (xs-m(:)*ones(1,N))';
  denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
  mahal = sum( (xs/C).*xs, 2 );
  numer = exp(-0.5*mahal);
  ps = numer/denom;
end

author: zouchuhang

repo: LayoutNet


text:
function sigma = sigupdate(nrad,nphi,nits,vol,mu0,maxrad,polafm2,lastsigma)
%   sigupdate: update RG kernel parameters
% G.A. Reina 16 Jan 2007
% Modified from the C code provided by D. L. Jones and R. G. Baraniuk
% "An Adaptive Optimal-Kernel Time-Frequency Representation"
%   by D. L. Jones and R. G. Baraniuk, IEEE Transactions on Signal
%   Processing, Vol. 43, No. 10, pp. 2361--2371, October 1995.

sigma = lastsigma;

for ii=0:(nits-1),
    gradsum = 0.0;
    gradsum1 = 0.0;

    for i=0:(nphi-1),
        grad(i+1) = 0.0;

        ee1 = exp( - 1.0/(sigma(i+1).^2) );	% use Kaiser's efficient method
        ee2 = 1.0;
        eec = ee1*ee1;

        for j=1:(maxrad(i+1)-1)
            ee2 = ee1*ee2;
            ee1 = eec*ee1;
            grad(i+1) = grad(i+1) + (j.^3)*ee2*polafm2(j+1, i+1);
        end
        grad(i+1) = grad(i+1)/(sigma(i+1).^3);

        gradsum = gradsum + grad(i+1).^2;
        gradsum1 = gradsum1 + sigma(i+1)*grad(i+1);

    end

    gradsum1 = 2*gradsum1;

    if ( gradsum < 0.0000001 )
        gradsum = 0.0000001;
    end

    if ( gradsum1 < 0.0000001 )
        gradsum1 = 0.0000001;
    end

    mu = ( sqrt(gradsum1.^2 + 4.0*gradsum*vol*mu0) - gradsum1 ) / ( 2.0*gradsum );

    sigma = sigma + mu*grad;
    sigma(sigma < 0.5) = 0.5;
    tvol = sum(sigma.^2);

    volfac = sqrt(vol/tvol);

    sigma = volfac*sigma;

end


text:
classdef MimSetRotateTool < MimGuiPlugin
    % MimSetRotateTool. Gui Plugin for using a preset bone window/level
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimSetRotateTool is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Rotate'
        SelectedText = 'Rotate'
        ToolTip = 'Enables Matlab''s rotate tool'
        Category = 'Tools'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'rotate.png'
        Location = 14
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.SetControl('Rotate');
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && strcmp(gui_app.ImagePanel.Mode, MimModes.View3DMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = strcmp(gui_app.ImagePanel.SelectedControl, 'Rotate');
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%% Machine Learning Online Class
%  Exercise 7 | Principle Component Analysis and K-Means Clustering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     pca.m
%     projectData.m
%     recoverData.m
%     computeCentroids.m
%     findClosestCentroids.m
%     kMeansInitCentroids.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ================= Part 1: Find Closest Centroids ====================
%  To help you implement K-Means, we have divided the learning algorithm
%  into two functions -- findClosestCentroids and computeCentroids. In this
%  part, you shoudl complete the code in the findClosestCentroids function.
%
fprintf('Finding closest centroids.\n\n');

% Load an example dataset that we will be using
load('ex7data2.mat');

% Select an initial set of centroids
K = 3; % 3 Centroids
initial_centroids = [3 3; 6 2; 8 5];

% Find the closest centroids for the examples using the
% initial_centroids
idx = findClosestCentroids(X, initial_centroids);

fprintf('Closest centroids for the first 3 examples: \n')
fprintf(' %d', idx(1:3));
fprintf('\n(the closest centroids should be 1, 3, 2 respectively)\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ===================== Part 2: Compute Means =========================
%  After implementing the closest centroids function, you should now
%  complete the computeCentroids function.
%
fprintf('\nComputing centroids means.\n\n');

%  Compute means based on the closest centroids found in the previous part.
centroids = computeCentroids(X, idx, K);

fprintf('Centroids computed after initial finding of closest centroids: \n')
fprintf(' %f %f \n' , centroids');
fprintf('\n(the centroids should be\n');
fprintf('   [ 2.428301 3.157924 ]\n');
fprintf('   [ 5.813503 2.633656 ]\n');
fprintf('   [ 7.119387 3.616684 ]\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =================== Part 3: K-Means Clustering ======================
%  After you have completed the two functions computeCentroids and
%  findClosestCentroids, you have all the necessary pieces to run the
%  kMeans algorithm. In this part, you will run the K-Means algorithm on
%  the example dataset we have provided.
%
fprintf('\nRunning K-Means clustering on example dataset.\n\n');

% Load an example dataset
load('ex7data2.mat');

% Settings for running K-Means
K = 3;
max_iters = 10;

% For consistency, here we set centroids to specific values
% but in practice you want to generate them automatically, such as by
% settings them to be random examples (as can be seen in
% kMeansInitCentroids).
initial_centroids = [3 3; 6 2; 8 5];

% Run K-Means algorithm. The 'true' at the end tells our function to plot
% the progress of K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters, true);
fprintf('\nK-Means Done.\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============= Part 4: K-Means Clustering on Pixels ===============
%  In this exercise, you will use K-Means to compress an image. To do this,
%  you will first run K-Means on the colors of the pixels in the image and
%  then you will map each pixel on to it's closest centroid.
%
%  You should now complete the code in kMeansInitCentroids.m
%

fprintf('\nRunning K-Means clustering on pixels from an image.\n\n');

%  Load an image of a bird
A = double(imread('rainbow2.png'));

% If imread does not work for you, you can try instead
%   load ('bird_small.mat');

A = A / 255; % Divide by 255 so that all values are in the range 0 - 1

% Size of the image
img_size = size(A);

% Reshape the image into an Nx3 matrix where N = number of pixels.
% Each row will contain the Red, Green and Blue pixel values
% This gives us our dataset matrix X that we will use K-Means on.
X = reshape(A, img_size(1) * img_size(2), 3);

% Run your K-Means algorithm on this data
% You should try different values of K and max_iters here
K = 16;
max_iters = 10;

% When using K-Means, it is important the initialize the centroids
% randomly.
% You should complete the code in kMeansInitCentroids.m before proceeding
initial_centroids = kMeansInitCentroids(X, K);

% Run K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================= Part 5: Image Compression ======================
%  In this part of the exercise, you will use the clusters of K-Means to
%  compress an image. To do this, we first find the closest clusters for
%  each example. After that, we

fprintf('\nApplying K-Means to compress an image.\n\n');

% Find closest cluster members
idx = findClosestCentroids(X, centroids);

% Essentially, now we have represented the image X as in terms of the
% indices in idx.

% We can now recover the image from the indices (idx) by mapping each pixel
% (specified by it's index in idx) to the centroid value
X_recovered = centroids(idx,:);
round(centroids * 255)

% Reshape the recovered image into proper dimensions
X_recovered = reshape(X_recovered, img_size(1), img_size(2), 3);

% Display the original image
subplot(1, 2, 1);
imagesc(A);
title('Original');

% Display compressed image side by side
subplot(1, 2, 2);
imagesc(X_recovered)
title(sprintf('Compressed, with %d colors.', K));

fprintf('Program paused. Pres s enter to continue.\n');
pause;

author: benoitvallon

repo: coursera-machine-learning


text:
function colors = prtPlotUtilLightenColors(colors)
% Internal function,
% xxx Need Help xxx

colors = colors + 0.2;
colors(colors > 1) = 1;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function d =  testing(a,d)

 Kt=get_x(a.Xsv)*get_x(d)';
 Yt=get_y(d);

 Yest=((a.alpha'* Kt)+a.b0)';
 if a.algorithm.use_signed_output
   Yest=sign(Yest);
 end

 d=set_x(d,Yest);
 d=set_name(d,[get_name(d) ' -> ' get_name(a)]);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function b = isPerpendicular(v1, v2, varargin)
%ISPERPENDICULAR Check orthogonality of two vectors
%
%   B = isPerpendicular(V1, V2)
%   where V1 and V2 are two 1-by-2 row arrays, returns 1 if the vectors are
%   perpendicular, and 0 otherwise.
%
%   Also works when V1 and V2 are two N-by-2 arrays with same number of
%   rows. In this case, return a N-by-1 array containing 1 at the positions
%   of perpendicular vectors.
%
%   Also works when one of V1 or V2 is 1-by-2 and the other one is a N-by-2
%   array. In this case the result has size N-by-1.
%
%   B = isPerpendicular(V1, V2, ACCURACY)
%   specifies accuracy of numerical tests, default is 1e-14.
%
%
%   Example
%   isPerpendicular([1 2 1], [2 4 2])
%   ans =
%       1
%
%   isPerpendicular([1 2 1], [1 3 2])
%   ans =
%       0
%
%   See also
%   vectors2d, isParallel, lines2d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-04-25
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   2007-09-18 copy from isPerpendicular, adapt to any dimension, and add
%       psb to specify precision
%   2009-09-21 fix bug for array of 3 vectors
%   2011-01-20 replace repmat by ones-indexing (faster)
%   2017-08-31 use normalized vectors

% default accuracy
acc = 1e-14;
if ~isempty(varargin)
    acc = abs(varargin{1});
end

% normalize vectors
v1 = normalizeVector(v1);
v2 = normalizeVector(v2);

% adapt size of inputs
n1 = size(v1, 1);
n2 = size(v2, 1);
if n1 ~= n2
    if n1 == 1
        v1 = v1(ones(n2, 1), :);
    elseif n2==1
        v2 = v2(ones(n1, 1), :);
    else
        error('Inputs must either have same size, or one must be scalar');
    end
end

% performs test
b = abs(dot(v1, v2, 2)) < acc;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function determ = conference_determinant ( n )

%*****************************************************************************80
%
%% CONFERENCE_DETERMINANT returns the determinant of the CONFERENCE matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.  N-1 must be an odd prime,
%    or a power of an odd prime.
%
%    Output, real DETERM, the determinant.
%
  if ( mod ( n - 1, 4 ) == 1 )
    determ = - sqrt ( ( n - 1 )^n );
  else
    determ = + sqrt ( ( n - 1 )^n );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% (Internal) Display as formated strings an XML node array
%
%
% See also read_hl7a_format
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 10/05/2016
% Birthdate  : 10/05/2016
% Copyright 2008-2016
%
function strAux = disp_xml_nodes(thisNode)

    strAux = [];

    for jj = 0:(thisNode.getLength-1)

        element = thisNode.item(jj);

        strAux = sprintf( '%s <%s' , strAux, char(element.getTagName) );

        thisCode_att = element.getAttributes;

        for ii = 0:(thisCode_att.getLength-1)

            this_att = thisCode_att.item(ii);

            strAux = sprintf( '%s %s ' , strAux, char(this_att.toString) );

        end

        strAux = sprintf( '%s>\n' , strAux );

    end

author: marianux

repo: ecg-kit

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


text:
function [units, data_description] = segtype2units(unitCode)
   %segtype2units   retrieves unit and description based on the segtype code
   %'segtype' in antelope datasets indicate the natural units of the detector
   persistent codeKey
   if isempty(codeKey)
      codeKey = createCodeKey();
   end
   if codeKey.isKey(unitCode)
      details = codeKey(unitCode);
      units = details{1};
      data_description = details{2};
   else
      [units, data_description] = deal('null');
   end

   function SU = createCodeKey()
      %createSegUnits   creates a map  where SU(char) = {units, data_description}
      SU = containers.Map;
      SU('A') = {'nm / sec / sec','acceleration'};
      SU('B') = {'25 mw / m / m','UV (sunburn) index(NOAA)'};
      SU('D') = {'nm', 'displacement'};
      SU('H') = {'Pa','hydroacoustic'};
      SU('I') = {'Pa','infrasound'};
      SU('J') = {'watts','power (Joulses/sec) (UCSD)'};
      SU('K') = {'kPa','generic pressure (UCSB)'};
      SU('M') = {'mm','Wood-Anderson drum recorder'};
      SU('P') = {'mb','barometric pressure'};
      SU('R') = {'mm','rain fall (UCSD)'};
      SU('S') = {'nm / m','strain'};
      SU('T') = {'sec','time'};
      SU('V') = {'nm / sec','velocity'};
      SU('W') = {'watts / m / m', 'insolation'};
      SU('a') = {'deg', 'azimuth'};
      SU('b') = {'bits/ sec', 'bit rate'};
      SU('c') = {'counts', 'dimensionless integer'};
      SU('d') = {'m', 'depth or height (e.g., water)'};
      SU('f') = {'micromoles / sec / m /m', 'photoactive radiation flux'};
      SU('h') = {'pH','hydrogen ion concentration'};
      SU('i') = {'amp','electric curent'};
      SU('m') = {'bitmap','dimensionless bitmap'};
      SU('n') = {'nanoradians','angle (tilt)'};
      SU('o') = {'mg/l','diliution of oxygen (Mark VanScoy)'};
      SU('p') = {'percent','percentage'};
      SU('r') = {'in','rainfall (UCSD)'};
      SU('s') = {'m / sec', 'speed (e.g., wind)'};
      SU('t') = {'C','temperature'};
      SU('u') = {'microsiemens/cm','conductivity'};
      SU('v') = {'volts','electric potential'};
      SU('w') = {'rad / sec', 'rotation rate'};
      SU('-') = {'null','null'};
   end
end


text:
function [rss,bic,aic ] = ml_kmeans_eval(X,labels,C,dist,lambda)
%ML_KMEANS_EVAL
%
%   input ----------------------------------------------------------------
%
%       o   X    :   (N x D), dataset of N samples of dimension D.
%
%       o labels :   (N x 1), class labels
%
%       o   C    :   (K x D), K centroids of dimension D.
%
%       o dist   :   string,  distance metric used.
%
%       o lambda :   (1 x 1), regulariser for AIC
%
%   output ----------------------------------------------------------------
%
%       o rss   :   (1 x 1), residual sum of squares
%
%       o aic   :   (1 x 1), Akaike information criteria
%
%       o bic   :   (1 x 1), Bayesian information criterion
%
%   info ------------------------------------------------------------------
%
%   http://nlp.stanford.edu/IR-book/html/htmledition/cluster-cardinality-in-k-means-1.html
%

N       = size(X,1);
[K,D]   = size(C);

if ~exist('lambda','var'),lambda = 2; end

% compute RSS
rss       = 0;
for k=1:K
    idx     = k==labels;
    rss     = rss + sum(sqrt(ml_distfunc(X(idx,:), C(k,:), dist)).^2);
end

aic = ml_aic(-0.5*rss,D * K,lambda);
bic = ml_bic(-0.5*rss,D * K,N);

end

author: epfl-lasa


text:
function day_list_common ( y1, m1, d1, y2, m2, d2 )

%*****************************************************************************80
%
%% DAY_LIST_COMMON prints a list of days between two dates.
%
%  Discussion:
%
%    Given the dates of September 25, 2005 and October 2, 2005,
%    the routine should print out:
%
%    Sun, Sep 25 2005 -
%    Mon, Sep 26 2005 -
%    Tue, Sep 27 2005 -
%    Wed, Sep 28 2005 -
%    Thu, Sep 29 2005 -
%    Fri, Sep 30 2005 -
%    Sat, Oct 01 2005 -
%    Sun, Oct 02 2005 -
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, the first date.
%
%    Input, integer Y2, M2, D2, the second date.
%
  y = y1;
  m = m1;
  d = d1;
  f = 0.0;

  cmp = '<';

  while ( cmp ~= '>' )

    w = ymdf_to_weekday_common ( y, m, d, f );

    w_name = weekday_to_name_common3 ( w );

    m_name = month_to_month_name_common3 ( m );

    fprintf ( 1, '%3s, %3s %02d %4d -\n', w_name, m_name, d, y );

    [ y, m, d, f ] = ymdf_next_common ( y, m, d, f );

    cmp = ymdf_compare ( y, m, d, f, y2, m2, d2, f );

  end

  return
end


text:
function x = grid1 ( ndim, nstep, x1, x2 )

%*****************************************************************************80
%
%% GRID1 finds grid points between X1 and X2 in N dimensions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NDIM, the dimension of the points X1 and X2.
%
%    Input, integer NSTEP, the number of points to be generated.
%    NSTEP must be at least 2.
%
%    Input, real X1(NDIM), X2(NDIM), the first and last
%    points, between which the equally spaced points are
%    to be computed.
%
%    Output, real X(NDIM,NSTEP), the set of equally spaced
%    points.  Each column of X represents one point, with X(*,1) = X1
%    and X(*,NSTEP) = X2.
%
  if ( nstep <= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GRID1 - Fatal error!\n' );
    fprintf ( 1, '  NSTEP <= 1.\n' );
    fprintf ( 1, '  NSTEP = %d\n', nstep );
    error ( 'GRID1 - Fatal error!' );
  end

  for i = 1 : nstep
    x(1:ndim,i) = ...
      ( ( nstep - i     ) * x1(1:ndim)'   ...
      + (         i - 1 ) * x2(1:ndim)' ) ...
      / ( nstep     - 1 );
  end

  return
end


text:
function V = est_calcInvCovMatFourier(Rinv,E,foi,fs,M,p, verb)
%
% Obtain the frequency domain transform of the (inverse) covariance matrix
% of an M-variate VAR[p] process.
%
% Let
% This is needed for calculating the renormalized PDC [2]
%
% Inputs:
%
%     Rinv: inverse process covariance matrix obtained from est_calcInvCovMat
%     E:    noise covariance matrix
%     foi:  frequencies of interest (Hz)
%     fs:   sampling rate
%     M:    # chans
%     p:    model order
%     verb: verbosity level. 0 = no output, 1=text.
%
% Outputs:
%
%     V:    Frequency-domain transform of the inverse covariance matrix
%
% References:
%
% [1] Lutkepohl, H. (2007) New Introduction to Time Series Analysis. Springer.
% [2] Schelter et al, (2009). Assessing the strength of directed influences
% among neural signals using renormalized partial directed coherence.
% Journal of Neuroscience Methods. 179:121-130.
%
% See Also: est_calcInvCovMat()
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

DEBUG = 0;

%% extract the diagonal elements of H=Rinv
% structure of Hd is (e.g., for p=2):
% [diag(H(1,1)), diag(H(2,1)), diag(H(1,2)), diag(H(2,2))]
% where diag(H(u,v)) is the column vector formed by the diagonal of
% submatrix H(u,v) of H
Hd = zeros(M,p^2);
cnt=1;
for v=1:p
    for u=1:p
        Hd(:,cnt)=diag(Rinv((u-1)*M+1:u*M,(v-1)*M+1:v*M));
        cnt=cnt+1;
    end
end

if DEBUG
    for j=1:N
        try chol(reshape(Hd(j,:),[p p])); catch; fprintf('Hj<0: j=%d - ',j); keyboard; end;
    end
end

%% create u,v index vectors
% us = [1 2 ... p 1 2 ... p ... ]'   (p^2 length)
% vs = [1 1 ... 1 2 2 ... 2 ... ]'   (p^2 length)
us = repmat(1:p,1,p)';
vs = zeros(p^2,1);
for ii=1:p
    vs((ii-1)*p+1:(ii-1)*p+p)=ones(p,1)*ii;
end

%% construct the V matrix for all freqs
fi=0;
COS=zeros(p^2,4);
freqs=(2*pi*foi)/fs;
V = zeros(length(freqs),M,M,2,2); % NOTE: V will end up (M,M,freqs,2,2)
if verb, h=waitbar(0,'calculating V^-1...'); end
for f=freqs
    fi=fi+1;

    %% construct cosine matrix
    COS(:,1) = cos(us*f).*cos(vs*f);
    COS(:,2) = sin(us*f).*cos(vs*f);
    COS(:,3) = cos(us*f).*sin(vs*f);
    COS(:,4) = sin(us*f).*sin(vs*f);

    %% multiply Hd and COS matrices to get matrix where row j is
    % [sum_{u,v=1 : p} Hjj(u,v)COS_11(u,v), ...
    %  sum_{u,v=1 : p} Hjj(u,v)COS_21(u,v), ...
    %  sum_{u,v=1 : p} Hjj(u,v)COS_12(u,v), ...
    %  sum_{u,v=1 : p} Hjj(u,v)COS_22(u,v)]
    %
    % where COSab(u,v) is the a,bth element of the sine-cos transform
    % matrix evaluated at u,v.
    %
    % NOTE: reshaping row j to 2x2 yeilds a matrix proportional to V_ij(f)
    % for some specific i
    Vm = Hd*COS;

    % now we multiply in the variances of the i's (E(i,i)) to generate the
    % full V matrix
    Vm = kron(diag(E)',Vm);   % Vm = kron(diag(p^2*E)',Vm);

    % next, reshape to desired structure (NxNx2x2)
    V(fi,:,:,:,:) = permute(reshape(Vm',2,2,M,M),[3 4 1 2]);

    if verb, waitbar(fi/length(freqs),h); end
end % for freqs

% permute V: (chs,chs,freqs,2,2)
V = permute(V,[2 3 1 4 5]);
if verb, close(h); end

repo: REST


text:
classdef TestInterpolateCornersCharuco
    %TestInterpolateCornersCharuco

    methods (Static)
        function test_1
            [img, board] = get_image_markers();
            [corners, ids] = cv.detectMarkers(img, board{end});
            [charucoCorners, charucoIds, num] = cv.interpolateCornersCharuco(...
                corners, ids, img, board);
            validateattributes(charucoCorners, {'cell'}, {'vector'});
            cellfun(@(c) validateattributes(c, {'numeric'}, ...
                {'vector', 'numel',2}), charucoCorners);
            validateattributes(charucoIds, {'numeric'}, ...
                {'vector', 'integer', 'nonnegative'});
            assert(isequal(numel(charucoCorners),numel(charucoIds)));
            validateattributes(num, {'numeric'}, ...
                {'scalar', 'integer', 'nonnegative'});
            assert(isequal(num,numel(charucoIds)));
        end

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

end

function [img, board] = get_image_markers()
    % markers in a 5x7 charuco board
    board = {5, 7, 60, 40, '6x6_50'};
    img = cv.drawCharucoBoard(board, [340 460], 'MarginSize',20);
    img = repmat(img, [1 1 3]);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function gBest = gBest_get(AA)
% The gBest (global best) updating strategy of S-ECSO

%------------------------------- 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".
%--------------------------------------------------------------------------

    [M,~]     = size(AA.decs);
    r         = rand(1,size(AA.objs,2));
    r_matr    = repmat(r,M,1);
    f_gBest   = sum(r_matr.*AA.objs,2)/sum(r);
    [~,index] = min(f_gBest);
    B         = AA.decs;
    gBest     = B(index,:);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function SO3F = conj(SO3F)
% Construct the complex conjugate function $\overline{f}$ of an SO3Fun $f$.
%
% Syntax
%   SO3F = conj(F)
%
% Input
%  F - @SO3FunHarmonic
%
% Output
%  SO3F - @SO3FunHarmonic
%

bw = SO3F.bandwidth;
for n=0:bw
  ind = deg2dim(n)+1:deg2dim(n+1);
  ind2 = flip(ind);
  SO3F.fhat(ind,:) = SO3F.fhat(ind2,:);
end
SO3F.fhat = conj(SO3F.fhat);

end

author: mtex-toolbox

repo: mtex