text:
% The static class geometricTools encapsulates the most used functions for surface processing
% used to solve EEG forward or inverse problems.
%
% Author: Alejandro Ojeda, SCCN/INC/UCSD, 2012
%
% Contributors:
% nonrigid_version23    -> D.Kroon, University of Twente, August
%                             2010 (http://www.mathworks.com/matlabcentral/fileexchange/20057)
% getSurfaceLaplacian   -> Nelson Trujillo Barreto and Pedro Antonio Valdes
%                             Hernandez, Cuban Neuroscience Center
% iso2mesh dependencies -> Qianqian Fang, http://iso2mesh.sourceforge.net/cgi-bin/index.cgi

classdef geometricTools
    methods
        function obj = geometricTools()
        end
        %%
        function disp(obj)
            disp(['Static class: ' class(obj)])
            disp('Static classes work like a namespace, they are used to create data and functions that can be accessed without creating an instance of the class.')
            methods(obj);
        end
    end
    methods(Static)
        %%
        function Xcentered = correctOrigin(X)
            [m,n] = size(X);
            K = ones(m,n);
            B = pinv(K'*K)*K'*X;
            X0 = K*B;
            Aff = eye(4);
            Aff([1 2],4) = X0(1,[1 2])';
            Xcentered = Aff\[X ones(m,1)]';
            Xcentered = Xcentered(1:3,:)';
        end
        %%
        function [Aff,Sn, scale] = affineMapping(S,T)
            % S: source space
            % T: target space
            % S = [sx1 sy1 sz1; sx2 sy2 sz2; ... sxk syn szk]
            % T = [tx1 ty1 tz1; tx2 ty2 tz2; ... txk tyn tzk]
            %
            % d(Aff) = min frobenius(T -S* Aff')
            % Sn = S*Aff';

            [~,~,transform] = procrustes(T,S);
            scale = transform.b;
            Aff = [[transform.b*transform.T;transform.c(1,:)] [0;0;0;1]]';
            Sn = geometricTools.applyAffineMapping(S,Aff);
        end
        %%
        function T = applyAffineMapping(S,M)
            % [T 1] = [S 1]*M';
            T = [S ones(size(S,1),1)]*M';
            T(:,4) = [];
        end
        %%
        function [def,spacing,offset,SgridWarped] = bSplineMapping(S,T,Sgrid,options)
            if nargin < 4,
                options.Verbose = false;
                options.MaxRef = 5;
            end
            mn = min(Sgrid);
            Smn = bsxfun(@minus,S,mn);
            dim = max(Sgrid) - mn;
            Tmn = bsxfun(@minus,T,mn);
            [def,spacing,SgridWarped] = point_registration(dim,Smn,Tmn,options);
            offset = mn;
            SgridWarped = bsxfun(@plus,SgridWarped,offset);
        end
        %%
        function SgridWarped = applyBSplineMapping(def,spacing,offset,Sgrid)
            Smn = bsxfun(@minus,Sgrid,offset);
            SmnWarped = bspline_trans_points_double(def,spacing,Smn);
            SgridWarped = bsxfun(@plus,SmnWarped,offset);
        end
        %%
        function [neighbors,D,loc] = nearestNeighbor(vertices,T)
            if exist('DelaunayTri','file')
                 dt = DelaunayTri(vertices(:,1),vertices(:,2),vertices(:,3)); %#ok
            else dt = delaunayTriangulation(vertices(:,1),vertices(:,2),vertices(:,3));
            end
            try [loc,D] = nearestNeighbor(dt, T);
            catch
                loc = nearestNeighbor(dt, T);
                D = sqrt(sum((vertices(loc,:)-T).^2,2));
            end
            neighbors = vertices(loc,:);
        end
        %%
        function Yi = ridgeInterpolation(vertices,faces,elec,Y)
            L = geometricTools.getSurfaceLaplacian(vertices,faces);
            K = geometricTools.localGaussianInterpolator(vertices,elec,1);
            K = full(K);
            Yi = ridgeGCV(Y,K,L,100,0);
        end
        %%
        function J = simulateGaussianSource(X,X0,h)
            if nargin < 3, h = 0.1;end
            J = geometricTools.localGaussianInterpolator(X,X0,h)';
        end
        function W = localGaussianInterpolator(X,Xi,h,normalize)
            if nargin < 3, h = 0.1;end
            if nargin < 4, normalize = false;end
            N = size(Xi,1);
            M = size(X,1);
            W = zeros(N,M);
            for it=1:N
                d = sum(bsxfun(@minus,X,Xi(it,:)).^2,2);
                W(it,:) = exp(-d/(2*h^2));
            end
            if normalize, W = bsxfun(@rdivide,W,sum(W,2)+eps);end
        end
        %%
        function D = isInConvexHull(X,Xi)

            X = geometricTools.correctOrigin(X);
            X = bsxfun(@rdivide,X,sqrt(sum(X.^2,2)));

            [x,y,z] = sphere(64);

            figure;
            surf(x,y,z,'FaceColor','g','FaceAlpha',0.7,'EdgeColor','none');
            set(gca,'Projection','perspective','DataAspectRatio',[1 1 1]); hold on;axis tight;camlight

            plot3(X(:,1),X(:,2),X(:,3),'.')

            X = geometricTools.correctOrigin(X);
            X = bsxfun(@rdivide,X,sqrt(sum(X.^2,2)));

            N = size(Xi,1);
            M = size(X,1);
            D = zeros(M,N);
            for it=1:N
                d = bsxfun(@minus,X,Xi(it,:));
                D(:,it) = sqrt(sum(( d ).^2,2));
            end
        end
        %%
        function Yi = localGaussianScatterInterpolator(X,Y,Xi)
            n = size(Y,2);
            if n==1
                F = TriScatteredInterp(X(:,1),X(:,2),X(:,3),Y,'nearest');
                Yi = F(Xi(:,1),Xi(:,2),Xi(:,3));
            else
                Yi = zeros(size(Xi,1),n);
                for it=1:n
                    F = TriScatteredInterp(X(:,1),X(:,2),X(:,3),Y(:,it),'nearest');
                    Yi(:,it) = F(Xi(:,1),Xi(:,2),Xi(:,3));
                end
            end
            W = geometricTools.localGaussianInterpolator(Xi,Xi,16);
            Yi = W*Yi;
        end
        %%
        function [rVertices,rFaces] = resampleSurface(vertices,faces,decimationPercent)
            if nargin < 2, error('Not enough input arguments.');end
            if nargin < 3, decimationPercent = 0.1;end
            if isempty(which('meshresample')), error('This function uses Iso2Mesh toolbox, you can download it for free fom: http://iso2mesh.sourceforge.net');end
            [rVertices,rFaces]=meshresample(vertices,faces,decimationPercent);
        end
        %%
        function sVertices = smoothSurface(vertices,faces,lambda,method)
            if nargin < 2, error('Not enough input arguments.');end
            if nargin < 3, lambda = 0.2;end
            if nargin < 4, method = 'lowpass';end

            maxIter = 20;
            N = size(vertices,1);

            if isempty(which('meshresample'))
                warning('MoBILAB:noIso2Mesh','This function uses Iso2Mesh toolbox if is installed, you can download it for free fom: http://iso2mesh.sourceforge.net');
                sVertices = vertices;
                for it=1:N
                    ind = any(faces==it,2);
                    indices = faces(ind,:);
                    indices = indices(:);
                    indices(indices==it) = [];
                    W = geometricTools.localGaussianInterpolator(vertices(indices,:),vertices(it,:),10);
                    sVertices(it,:) = sum((1./W)*vertices(indices,:),1)./sum(1./W);
                end
                return;
            end
            conn = neighborelem(faces,size(vertices,1));
            for it=1:N
                tmp = faces(conn{it},:);
                conn{it} = unique(tmp(:)');
            end
            sVertices = smoothsurf(vertices,[],conn,maxIter,lambda,method);
        end
        %%
        function [rVertices,rFaces] = refineSurface(vertices,faces,decimationRate,maxIter)
            if nargin < 3, decimationRate = 0.5;end
            if nargin < 4, maxIter = 3;end
            if isempty(which('meshresample')), error('This function uses Iso2Mesh toolbox, you can download it for free fom: http://iso2mesh.sourceforge.net');end

            tmpVertices = vertices;
            tmpFaces = faces;
            for it=1:maxIter
                [tmpVertices,tmpFaces] = geometricTools.resampleSurface(tmpVertices,tmpFaces,decimationRate);
                tmpVertices = geometricTools.smoothSurface(tmpVertices,tmpFaces);
                disp(it)
            end
            rVertices = tmpVertices;
            rFaces = tmpFaces;
        end
        %%
        function verticesExt = repareIntersectedSurface(surfInt,surfOut,dmax)
            if nargin < 3, dmax = 8;end
            verticesInt = surfInt.vertices;
            verticesExt = surfOut.vertices;
            [nVerticesInt,d] = geometricTools.nearestNeighbor(verticesExt,verticesInt);
            I = d < dmax;
            while any(I)
                I2 = ismember(verticesExt,nVerticesInt(I,:),'rows');
                verticesExt(I2,:) = 1.005*verticesExt(I2,:);
                [nVerticesInt,d] = geometricTools.nearestNeighbor(verticesExt,verticesInt);
                I = d < dmax;
            end
            if any(verticesExt(:) ~= surfOut.vertices(:))
                verticesExt = geometricTools.smoothSurface(verticesExt,surfOut.faces);
            end
        end
        %%
        function [normals,faces] = getSurfaceNormals(vertices,faces,normalsIn)
            if nargin < 3, normalsIn = true;end
            h = figure('visible','off');
            h2 = patch('vertices',vertices,'faces',fliplr(faces));
            normals = get(h2,'vertexnormals');close(h);
            if isempty(normals)
                normals = vertices;
            end
            normals = normals./(sqrt(sum(normals.^2,2))*[1 1 1]);
            area1 = geometricTools.getSurfaceArea(vertices,faces);
            area2 = geometricTools.getSurfaceArea(vertices+normals,faces);
            if area2 < area1% && normalsIn
                faces = fliplr(faces);
                h = figure('visible','off');h2 = patch('vertices',vertices,'faces',fliplr(faces));
                normals = get(h2,'vertexnormals');close(h);
                normals = normals./(sqrt(sum(normals.^2,2))*[1 1 1]);
            end
            if normalsIn
                faces = fliplr(faces);
                h = figure('visible','off');h2 = patch('vertices',vertices,'faces',fliplr(faces));
                normals = get(h2,'vertexnormals');close(h);
                normals = normals./(sqrt(sum(normals.^2,2))*[1 1 1]);
            end
        end
        %%
        function [area,areas] = getSurfaceArea(vertices,faces)
            x1= vertices(faces(:,1),1);
            y1= vertices(faces(:,1),2);
            z1= vertices(faces(:,1),3);
            x2= vertices(faces(:,2),1);
            y2= vertices(faces(:,2),2);
            z2= vertices(faces(:,2),3);
            x3= vertices(faces(:,3),1);
            y3= vertices(faces(:,3),2);
            z3= vertices(faces(:,3),3);
            area = sqrt(((y2-y1).*(z3-z1)-(y3-y1).*(z2-z1)).^2+((z2-z1).*(x3-x1)-(z3-z1).*(x2-x1)).^2+...
                ((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)).^2)/2;
            areas = area;
            area = sum(area);
        end
        function L = getSurfaceLaplacian1(vertices,faces)
            % LAPLACES Calculates a Discrete Surface Laplacian Matrix
            %          for a triangulated surface

            Nv = size(vertices,1);
            Nf = size(faces,1);
            L = spalloc(Nv,Nv,3*Nf);
            %L = speye(Nvtx);
            for fi=1:Nf
                for k=1:3
                    kk = mod(k,3)+1;
                    L(faces(fi,k),faces(fi,kk)) = sqrt( sum((vertices(faces(fi,k),:)-vertices(faces(fi,kk),:)).^2,2));
                end
            end
            L(L>0) = 1./L(L>0);
            L = (L+L')/2;
            L = L - spdiags(sum(L,2),0,Nv,Nv);
        end
        function L = getSurfaceLaplacian(vertices,faces)
            % LAPLACES Calculates a Discrete Surface Laplacian Matrix
            %          for a triangulated surface
            %
            % Reference:
            % [1] Huiskamp, G., 1991, Difference formulas for the surface laplacian
            %     on a triangulated surface, Journal of Computational Physics 95,
            %     477-496.
            %
            % Nelson Trujillo Barreto
            % Pedro antonio Valdes Hernandez
            % Cuban Neuroscience Center

            Cortex.vertices = vertices;
            Cortex.faces = faces;
            vtx = Cortex.vertices;
            tri = Cortex.faces;
            [nei,nei_tri] = geometricTools.get_neis(Cortex);
            Nvtx = size(vtx,1);
            L = speye(Nvtx);
            for j=1:Nvtx,
                nei_tri_j = tri(nei_tri{j},:);
                nei_j = nei{j};
                PHI_jk = [];
                rj = vtx(j,:);
                Nj = length(nei_j);
                for k=1:Nj,

                    [indi,indj]=find(nei_tri_j==nei_j(k)); %#ok
                    nei_tri_jk = nei_tri_j(indi,:);
                    if size(nei_tri_jk,1) < 2,
                        break;
                    end
                    nei_k_lr = setxor(nei_tri_jk(1,:),nei_tri_jk(2,:));

                    rk2 = sum((rj-vtx(nei_j(k),:)).^2);
                    rl2 = sum((rj-vtx(nei_k_lr(1),:)).^2);
                    rr2 = sum((rj-vtx(nei_k_lr(2),:)).^2);
                    rkl2 = sum((vtx(nei_k_lr(1),:)-vtx(nei_j(k),:)).^2);
                    rkr2 = sum((vtx(nei_k_lr(2),:)-vtx(nei_j(k),:)).^2);

                    cos_phi_kl = (rk2+rl2-rkl2)./sqrt(rk2.*rl2)./2;
                    cos_phi_kr = (rk2+rr2-rkr2)./sqrt(rk2.*rr2)./2;
                    sin_phi_kl = sqrt(1-cos_phi_kl.^2);
                    sin_phi_kr = sqrt(1-cos_phi_kr.^2);

                    PHI_jk = [PHI_jk (1-cos_phi_kl)./(sin_phi_kl+eps)+(1-cos_phi_kr)./(sin_phi_kr+eps)]; %#ok
                end
                if ~isempty(PHI_jk)
                    rjk = sqrt(sum((vtx(nei_j,:)-repmat(rj,Nj,1)).^2,2));
                    rj_bar = mean(rjk);

                    theta_jk = 4*PHI_jk'./(rj_bar*sum(PHI_jk)*rjk);
                    L(j,nei_j) = theta_jk'; %#ok
                else
                    L(j,nei_j) = 0;     %#ok
                    L(j,j) = 1;         %#ok
                end
            end;

            L = L - speye(Nvtx,Nvtx);
            L = L - spdiags(sum(L,2),0,Nvtx,Nvtx);
            d = diag(L);
            ind = find(d==0);
            if ~isempty(ind), for it=1:length(ind), L(ind(it),ind(it)) = 1;end;end
            %th = prctile(nonzeros(L),[0.1 99.9]);
            %L(L>th(2)) = 0;
            %L(L<th(1)) = 0;
        end
        %%
        function [nei,nei_tri] = get_neis(P)
            % helper function for getSurfaceLaplacian
            % Nelson Trujillo Barreto
            % Pedro antonio Valdes Hernandez
            % Cuban Neuroscience Center
            n = size(P.vertices,1);
            nei_tri = cell(n,1);
            nei = cell(n,1);
            for i = 1:n
                [r,c] = find(P.faces == i); %#ok
                nei_tri{i} = r;
                nei{i} = setdiff(unique(P.faces(r,:)),i);
            end
        end
        function [Nei_faces,Nei_vertices] = get_neis1(P)
            n = size(P.vertices,1);
            Nei_faces = cell(n,1);
            Nei_vertices = cell(n,1);
            % hbar = waitbar(0,'calculating neigs...');
            for i = 1:n
                [r,c] = find(P.faces == i); %#ok
                tmp = P.faces(r,:)';
                tmp(tmp == i) = [];
                m = length(tmp)/2;
                Nei_faces{i} = reshape(tmp,2,m)';
                for j = 1:m
                    Nei_vertices{i}{j} = P.vertices(Nei_faces{i}(j,:),:);
                end
                waitbar(i/n,hbar);
            end
            close(hbar);
        end
       %%
        function [nVertices,nFaces] = openSurface(vertices,faces,rmIndices)
            nVertices = vertices;
            vertices(rmIndices,:) = [];
            [~,rm_1] = ismember(faces(:,1),rmIndices);
            [~,rm_2] = ismember(faces(:,2),rmIndices);
            [~,rm_3] = ismember(faces(:,3),rmIndices);
            rm_faces = rm_1 | rm_2 | rm_3;
            faces(rm_faces,:) = [];
            [~,J] = ismember(nVertices,vertices,'rows');
            nFaces = J(faces);
            nVertices = vertices;
        end
        %%
        function [nVertices,nFaces] = getSurfaceROI(vertices,faces,roiIndices)
            rmIndices = setdiff(1:size(vertices,1),roiIndices);
            [nVertices,nFaces] = geometricTools.openSurface(vertices,faces,rmIndices);
        end
        %%
        function yi = interpOnSurface(vertices,faces,elec,y,method)
            if nargin < 5, method = 'spline';end
            switch method
                case 'ridge'
                    yi = geometricTools.ridgeInterpolation(vertices,faces,elec,y);
                case 'linear'
                    W = geometricTools.localGaussianInterpolator(elec,vertices,32);
                    yi = W*y;
                case 'spline'
                    yi = geometricTools.spSplineInterpolator(elec,y,vertices);
                otherwise
                    yi = geometricTools.spSplineInterpolator(elec,y,vertices);
            end
        end
        %%
        function atlas = labelSurface(Surf,imgAtlasfile, txtAtlasLabel,maxColorValue)
            if nargin < 4, maxColorValue = 90;end
            % Atlas
            v =spm_vol(imgAtlasfile); % atlas
            A = spm_read_vols(v);
            A(A>maxColorValue) = 0;
            indNonZero = A(:)~=0;
            colorTable = A(indNonZero);
            [x,y,z] = ndgrid(1:v.dim(1),1:v.dim(2),1:v.dim(3));
            M = v.mat;
            X = [x(:) y(:) z(:) ones(numel(x),1)]*M';
            X = X(indNonZero,1:3);
            clear x y z
            F = TriScatteredInterp(X,colorTable,'nearest');
            n = size(Surf.vertices,1);
            labelsValue = F(Surf.vertices);
            colorTable = labelsValue;
            hwait = waitbar(0,'Atlas correction...');
            for it=1:n
                neigInd = any(Surf.faces == it,2);
                vertexInedex = Surf.faces(neigInd,:);
                vertexInedex = vertexInedex(:);
                [y,x] = hist(labelsValue(vertexInedex));
                [~,loc] = max(y);
                [~,loc] = min(abs(labelsValue(vertexInedex) - x(loc)));
                labelsValue(it) = colorTable(vertexInedex(loc));
                waitbar(it/n,hwait);
            end
            waitbar(1,hwait);
            close(hwait);
            atlas.colorTable = labelsValue;
            atlas.label = textfile2cell(txtAtlasLabel);
            atlas.label = atlas.label(1:max(atlas.colorTable));
            for it=1:length(atlas.label)
                ind = find(atlas.label{it} == ' ');
                atlas.label{it} = atlas.label{it}(ind(1)+1:ind(end)-1);
            end
        end
        %%
        function X = projectOntoUnitarySphere(X)
            [~,X(:,1),X(:,2),X(:,3)] = geometricTools.projectOnSphere(X(:,1),X(:,2),X(:,3));
            [azimuth,elevation,r] = cart2sph(X(:,1),X(:,2),X(:,3));
            [X(:,1),X(:,2),X(:,3)] = sph2cart(azimuth,elevation,elevation*0+1);
        end
        %%
        function [Yi,W] = spSplineInterpolator(X,Y,Xi,plotFlag)
            % Computes the spherical spline interpolator based on Perrin, F.,
            % Pernier, J., Bertrand, O., Echallier, J.F. (1990). Corrigenda
            % EEG 02274. Electroencephalography and Clinical Neurophysiology, 76, 565.

            if nargin < 4, plotFlag = false;end
            %X0 = mean(X);
            %X  = bsxfun(@minus,X,X0);
            %Xi  = bsxfun(@minus,Xi,X0);
            X  = geometricTools.projectOntoUnitarySphere(X);
            Xi = geometricTools.projectOntoUnitarySphere(Xi);
            %X  = bsxfun(@rdivide,X,sqrt(sum(X.^2,2)));
            %Xi = bsxfun(@rdivide,Xi,sqrt(sum(Xi.^2,2)));

            %--
            M = size(X,1);
            One = ones(size(Xi,1),1);
            %--

            % Solving eq. 4 of Perrin et al. (1989)
            COS_X  = geometricTools.cosines(X,X);
            COS_Xi = geometricTools.cosines(Xi,X);

            % Solving eq. 3 of Perrin et al. (1989)
            Gx  = geometricTools.sphericalSpline(COS_X);
            Gxi = geometricTools.sphericalSpline(COS_Xi);

            % Solving eq. 2 Perrin et al. (1989)
            [C,~,~,T] = ridgeGCV([Y;0],[Gx ones(M,1);ones(1,M) 0],eye(M+1));

            % Interpolating with the spherical harmonics
            Yi = [Gxi One]* C;

            W = [Gxi One]* T(:,1:end-1);

            % Plot the input & projected electrode positions on a sphere
            if plotFlag
                geometricTools.plot_on_sphere(X,Y,Xi,Yi);
            end
        end
        %%
        function Gx = sphericalSpline(x)
            % sphericalSpline solves eq. 3 of Perrin et al. (1989)
            % g(COS) = 1/4pi * sum[n=1:inf] (( (2*n+1)/( n^m * (n+1)^m ) ) * Pn(COS));

            m = 4;
            N = 16;    % gives accuracy of 10^-6

            P = cat(3, ones(size(x)), x);
            Gx = 3 / 2 ^ m * P(:, :, 2);
            for n = 2:N
                P(:, :, 3) = ((2 * n - 1) * x .* P(:, :, 2) - (n - 1) * P(:, :, 1)) / n;
                P = P(:,:,[2 3 1]);
                Gx = Gx + (2 * n + 1) / (n ^ m * (n + 1) ^ m) * P(:, :, 2);
            end
            Gx = Gx / (4 * pi);
        end
        %%
        function [r,x,y,z] = projectOnSphere(X,Y,Z,xo,yo,zo)
            % projectOnSphere - calculates projections of xyz positions
            % onto the unitary sphere
            %
            % Usage: [r,x,y,z] = projectOnSphere(X,Y,Z,xo,yo,zo,plotFlag)
            %
            % Notes:    The general formula for a sphere, with radius r is given by:
            %
            %           (x - xo)^2  +  (y - yo)^2  +  (z - zo)^2  =  r^2
            %
            %           This function takes arguments for cartesian co-ordinates
            %           of X,Y,Z (assume Z > 0) and the center of the sphere (xo,yo,zo).
            %           If (xo,yo,zo) is not provided a cnter at (0,0,0) is assumed.
            %
            %           Returned values are the fitted radius 'r' (constant)
            %           and the (x,y,z) Cartesian coordinates of the projected points
            %
            %
            % $Revision: 1.3 $ $Date: 2005/07/12 22:16:48 $
            % Licence:  GNU GPL, no express or implied warranties
            % History:  02/2002, Darren.Weber_at_radiology.ucsf.edu
            %                    adapted from elec_fit_sph
            %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % initialise centroid, unless input parameters defined
            if nargin < 4, xo = 0;end
            if nargin < 5, yo = 0;end
            if nargin < 6, zo = 0;end

            % Initialise r0 as a rough guess at the sphere radius
            rX = (max(X) - min(X)) / 2;
            rY = (max(Y) - min(Y)) / 2;
            rZ =  max(Z) - zo;
            r0 = mean([ rX rY rZ ]);

            % perform least squares estimate of spherical radius (r)
            options = optimset('fminsearch');
            r = fminsearch(@geometricTools.fit2sphere,r0, options, X, Y, Z, xo, yo, zo);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Find the projection point of X,Y,Z to the fitted sphere radius r

            % Convert Cartesian X,Y,Z to spherical (radians)
            theta = atan2( (Y-yo), (X-xo) );
            phi = atan2( sqrt( (X-xo).^2 + (Y-yo).^2 ), (Z-zo) );
            % do not recalc: r = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2);

            %   Recalculate X,Y,Z for constant r, given theta & phi.
            R = ones(size(phi)) * r;
            x = R .* sin(phi) .* cos(theta);
            y = R .* sin(phi) .* sin(theta);
            z = R .* cos(phi);
        end
    end
    methods(Static,Hidden=true)
        %%
        function f = fit2sphere(r, X, Y, Z, xo, yo, zo)
            S = (X-xo).^2  +  (Y-yo).^2  +  (Z-zo).^2  -  r^2;
            f = sum( S.^2 );
        end
        %%
        function Cos = cosines(A,B)
            Na = size(A,1);
            Nb = size(B,1);
            One1 = ones(1,Nb);
            One2 = ones(Na,1);
            Xe = A(:,1)*One1;
            Ye = A(:,2)*One1;
            Ze = A(:,3)*One1;
            Xf = One2*B(:,1)';
            Yf = One2*B(:,2)';
            Zf = One2*B(:,3)';
            Cos = (Xe-Xf).^2 + (Ye-Yf).^2 + (Ze-Zf).^2;
            Cos = 1-Cos/2;
            Cos(Cos > 1) = 1-eps;
            Cos(Cos < -1) = -1+eps;
        end
        function plot_on_sphere(X,Y,Xi,Yi)
            [~,X(:,1), X(:,2), X(:,3)]  = geometricTools.projectOnSphere(X(:,1), X(:,2), X(:,3));
            [~,Xi(:,1),Xi(:,2),Xi(:,3)] = geometricTools.projectOnSphere(Xi(:,1),Xi(:,2),Xi(:,3));
            X  = bsxfun(@rdivide,X,sqrt(sum(X.^2,2)));
            Xi = bsxfun(@rdivide,Xi,sqrt(sum(Xi.^2,2)));

            Xt = [X;Xi];
            Yt = [Y;Yi];
            Xt = bsxfun(@rdivide,Xt,sqrt(sum(Xt.^2,2)));
            Xi = bsxfun(@rdivide,Xi,sqrt(sum(Xi.^2,2)));
            Ne = size(X,1);
            Nf = 72;
            [Xs,Ys,Zs]=sphere(Nf);
            Xsp = [Xs(:) Ys(:) Zs(:)];
            Fsp = geometricTools.localGaussianInterpolator(Xt,Xsp,0.2);

            %[J,lambdaOpt,~,iFsp] = ridgeGCV(Yt,Fsp',eye(size(Xsp,1)),100,1);
            %J = iFsp*Yt;
            Ysp = Fsp*Yt;
            figure('NumberTitle','off','Name','Electrode Placements');
            set(gca,'Projection','perspective','DataAspectRatio',[1 1 1]); hold on
            %plot3(x,y,z,'b.');
            plot3(X(:,1),X(:,2),X(:,3),'ro');
            plot3(Xi(:,1),Xi(:,2),Xi(:,3),'k.')
            legend('input xyz','projected head','Location','BestOutside');
            surf(Xs,Ys,Zs,reshape(Ysp,[Nf Nf]+1),'specularstrength',0.1,'facealpha',0.9,'linestyle','none');
            camlight
            camlight headlight
            view(2); rotate3d;
            axis vis3d
        end
    end
end

repo: REST


text:
function delete_wp(pos, host, port)
% DELETE_WP(POS, HOST, PORT)
%
% Delete waypoint at position POS in stack of Flightgear over socket
% connection (HOST, PORT).
% Run Flightgear on HOST (e.g. '127.0.0.1' = localhost) with
% parameter "--props = PORT".

import java.net.Socket
import java.io.*

fg_socket = Socket(host, port);
out = PrintWriter(fg_socket.getOutputStream(), true);

out.println(sprintf(...
    'set /autopilot/route-manager/input @delete%d\r\n', pos-1));

out.close();
fg_socket.close();

end

author: Sable

repo: mcbench-benchmarks


text:
function N = numAtomsOfElementInFormula(formula, element, printLevel)
% Returns the number of atoms of a single element in a formula
%
% USAGE:
%
%    N = numAtomsOfElementInFormula(formula, element, printLevel)
%
% INPUTS:
%    formula:       formula in format Element then NumberOfAtoms with no spaces
%    element:       Abbreviation of element e.g. C or Mg
%    printLevel:    verbose level
%
% OUTPUT:
%    N:             number of atoms of this element in the formula provided
%
% .. Author:
%       - Ronan Fleming 9 March 09
%       - Ronan Fleming 21 July 09 handles formulae with same first letter for element e.g.C & Co in 'C55H80CoN15O11'
%       - Ronan Fleming 18 Sept 09 handles composite formulae like C62H90N13O14P.C10H11N5O3.Co
%       - Hulda SH 7 July 2011 Simplified and generalized code. Now handles most or all exceptional formulas.

if nargin < 3
    printLevel = 1;
end

if any(~isletter(element))
    disp(element)
    error('Element must not contain numbers')
end

if length(element)==1
    if ~strcmp(upper(element),element)
        disp(element)
        error('Single letter element must not be lower case')
    end
end

if ~isletter(formula(1))
    disp(formula)
    if printLevel>0
        warning('Formula format expected is element then number of elements, not the other way around')
    end
end

zero='0';
indZero=strfind(formula,zero);
if ~isempty(indZero)
    if isletter(formula(indZero-1))
        fprintf('%s%s\n','Formula before:, ', formula)
        warning('Formula contains a zero with a letter preceeding it, replacing with letter.')
        formula = strrep(formula, formula(indZero-1:indZero), formula(indZero-1));
        fprintf('%s%s\n','Formula after:, ', formula)
    end
end

%Count FULLR and FULLR2 groups as R groups (Human reconstruction)
formula = strrep(formula, 'FULLR3', 'R3');
formula = strrep(formula, 'FULLR2', 'R2');
formula = strrep(formula, 'FULLR', 'R');

element = strrep(element, 'FULLR3', 'R3');
element = strrep(element, 'FULLR2', 'R2');
element = strrep(element, 'FULLR', 'R');

if ischar(element) && ischar(formula)
    elementStart = regexp(formula, '[A-Z]', 'start'); % Get indices of all capital letters in the string formula. Treated as starting indices of elements.
    elementArray = cell(length(elementStart),2); % Initialize cell array for element symbols (Column 1) and numerical subscripts (Column 2)

    for n = 1:length(elementStart) % Loop for each element in formula
        if n < length(elementStart)
            splitFormula = formula(elementStart(n):(elementStart(n+1)-1)); % Extract section of formula from starting index of element n to starting index of element n+1

        else
            splitFormula = formula(elementStart(n):end);

        end

        if ~isempty(regexp(splitFormula, '[^a-z_A-Z]', 'once')) % If current section of formula contains non-alphabetic characters
            elementArray{n,1} = splitFormula(1:(regexp(splitFormula, '[^a-z_A-Z]', 'once')-1)); % Element symbol assumed to extend from beginning of current section to the first non-alphabetic character
            rest = splitFormula(regexp(splitFormula, '[^a-z_A-Z]', 'once'):end); % Rest of section after element symbol.

            if ~isempty(regexp(rest, '\W', 'once')) % Rest of section may contain word characters such as numbers and dots.
                if ~isempty(regexp(rest(1), '\d', 'once')) % If first character following element symbol is numeric it is assumed to represent the number of atoms of that element.
                    elementArray{n,2} = rest(1:(regexp(rest, '\W', 'once')-1)); % Extract element's numeric subscript.

                else
                    elementArray{n,2} = '1'; % If no number follows element symbol the numeric subscript is assumed to be 1.

                end

            else
                elementArray{n,2} = rest;

            end

        else
            elementArray{n,1} = splitFormula;
            elementArray{n,2} = '1';

        end
    end

    elementRows = strmatch(element, elementArray(:,1), 'exact'); % Element may appear in two different locations within formula.

    if ~isempty(elementRows)
        elementCount = zeros(length(elementRows),1);

        for m = 1:length(elementRows)
            elementCount(m,1) = str2double(elementArray{elementRows(m),2}); % Get numeric subscript by each instance of element in formula and convert from char to double

        end

        N = sum(elementCount);

    else
        N = 0;

    end

else
    if ~ischar(element)
        disp(element)
        error('Element must be given by a variable of class char')
    else
        disp(element)
        error('Formula must be given by a variable of class char')
    end
end

author: opencobra

repo: cobratoolbox


text:
function [net, info] = trainModel(varargin)

opts.batchNormalization = false ;
opts.network = [];
opts.networkType = 'dagnn' ;
opts.database = 'polyu1000';
opts.model = 'vgg16';
opts.downsample = 5;
opts.cdim = 4;
opts.ndim = 1024;
opts.warpsize = 8;
opts.trainNum = 1000;
opts.RoIRoD = 'RoIRoD';
opts.seed = 1;

[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.networkType ;
if opts.batchNormalization, sfx = [sfx '-bnorm'] ; end
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.expDir = fullfile('data',opts.database) ;
opts.train = struct() ;
opts.train.derOutputs = {'losscls', 1} ;
opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = [1]; end;

switch opts.RoIRoD
    case 'RoIOnly'
        opts.expDir = fullfile(opts.expDir,opts.model,'RoIOnly');
    case 'RoDOnly'
        opts.expDir = fullfile(opts.expDir,opts.model,'RoDOnly');
    case 'RoIRoD'
        opts.expDir = fullfile(opts.expDir,opts.model,'RoIRoD');
    otherwise
        error('Unknown RoIRoD choice.');
end

%% prepare the model and data
switch opts.model
    case 'vgg16'
        net = vgg16_intialization(opts.downsample,opts.warpsize,opts.cdim,opts.ndim, opts.RoIRoD);
    case 'resnet50'
        net = resnet50_intialization(opts.downsample,opts.warpsize,opts.cdim,opts.ndim, opts.RoIRoD);
end

opts.minScale = 256;
opts.imdbPath = fullfile(['imdb_' opts.database num2str(opts.minScale) '.mat']);
opts.expDir = fullfile(opts.expDir,[num2str(opts.trainNum) '_down' num2str(opts.downsample)...
              '_warp' num2str(opts.warpsize) 'cdim' num2str(opts.cdim) 'ndim' num2str(opts.ndim) 'seed' num2str(opts.seed)]);

if strcmp(opts.database,'GAIC')
    imdb = setup_database_GAIC(opts.trainNum,opts.seed);
end

gt_scores = cat(2,imdb.bbox.gt_scores{imdb.images.set==1});
gt_scores_means = mean(gt_scores);
gt_scores_stds = std(gt_scores);
for i = 1:length(imdb.bbox.gt_scores)
    imdb.bbox.gt_scores{i} = (imdb.bbox.gt_scores{i} - gt_scores_means) / gt_scores_stds;
end
imdb.gt_scores.means = gt_scores_means;
imdb.gt_scores.stds = gt_scores_stds;

imdb.images.dataMean = net.meta.normalization.averageImage;
net.meta.minScale = opts.minScale;

% --------------------------------------------------------------------
%  Train
% --------------------------------------------------------------------

[net, info] = cnn_train_dag(net, imdb, @(i,b) getBatch_rois(opts,i,b),...
  'expDir', opts.expDir, net.meta.trainOpts, opts.train, ...
  'val', find(imdb.images.set == 3),'solver',@solver.adam);

% --------------------------------------------------------------------
function inputs = getBatch_rois(opts, imdb, batch)
% --------------------------------------------------------------------

imo = vl_imreadjpeg(imdb.meta.img_path(batch),'numThreads',1,'Contrast',0.5,'Saturation',0.5);
num_samples = 64;
for i = 1:length(batch)

    [x1,x2,x3] = size(imo{i});
    imre = imresize(imo{i},opts.minScale/min([x1,x2]),'bilinear');

    [r1,r2,r3] = size(imre);
    r1 = 32*round(r1/32);
    r2 = 32*round(r2/32);
    imre = imresize(imre,[r1,r2],'bilinear');

    imre = bsxfun(@minus,imre,imdb.images.dataMean);

    sel = randperm(numel(imdb.bbox.gt_scores{batch(i)}),min(num_samples,numel(imdb.bbox.gt_scores{batch(i)})));
    rois{i} = imdb.bbox.boxes{batch(i)}(:,sel);

    scale1 = r1/x1;
    scale2 = r2/x2;
    rois{i}(1,:) = max(floor(rois{i}(1,:) * scale1),1);
    rois{i}(2,:) = max(floor(rois{i}(2,:) * scale2),1);
    rois{i}(3,:) = min(ceil(rois{i}(3,:) * scale1),r1);
    rois{i}(4,:) = min(ceil(rois{i}(4,:) * scale2),r2);

    rois{i} = [i*ones(1,length(sel));rois{i}];
    gt_scores{i} = imdb.bbox.gt_scores{batch(i)}(sel);
end

rois = cat(2,rois{:});
gt_scores = cat(2,gt_scores{:});

if numel(opts.train.gpus) > 0
  imre = gpuArray(imre) ;
  rois = gpuArray(single(rois)) ;
end

inputs = {'input', imre, 'label', gt_scores, 'rois', rois} ;

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
function lengths = grEdgeLength(nodes, edges, varargin)
%GREDGELENGTH  Compute length of edges in a geometric graph.
%
%   Deprecated, use 'grEdgeLengths' instead.
%
%   LENGTHS = grEdgeLength(NODES, EDGES)
%
%   Example
%   grEdgeLength
%
%   See also
%     grEdgeLenghts

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2014-01-13, using Matlab 7.9.0.529 (R2009b)
% Copyright 2014 INRA - Cepia Software Platform

warning('MatGeom:graphs:deprecated', ...
    'function grEdgeLength is obsolete, use grEdgeLengths instead');

nEdges = size(edges, 1);
lengths = zeros(nEdges, 1);

for iEdge = 1:nEdges
    ed = edges(iEdge, :);
    node1 = nodes(ed(1),:);
    node2 = nodes(ed(2),:);
    lengths(iEdge) = sqrt(sum((node1 - node2).^2));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function ok = validate(logger)
% Checks that an object of type Logger obeys certain rules that
% must always be true.
% Inputs:
%   logger: the object to be checked.
% Outputs:
%   ok: a boolean value indicating whether the object is valid.

assert(nargin==1)
assert(isa(logger,'Logger'))

ok = length(logger.fidlist)==length(logger.fidlevels);
ok = ok && length(logger.filenamelist)==length(logger.filenamefids);
ok = ok && length(logger.filenamelist)==length(logger.filenamelevels);
ok = ok && all(logger.fidlist > 0);
ok = ok && all(logger.filenamefids > 0);
ok = ok && all(logger.fidlevels >= 0);
ok = ok && all(logger.fidlevels < 256);
ok = ok && all(logger.filenamelevels >= 0);
ok = ok && all(logger.filenamelevels < 256);

author: nesl


text:
% RAY_MESH_INTERSECT  Find first hit (if it exists) for each ray.
%
% [flag, t, lambda] = ray_mesh_intersect(src, dir, V, F);
%
% Input:
%    src #rays by 3 list of 3D vector ray origins
%    dir #rays by 3 list of 3D vector ray directions
%    V  #V by 3 list of vertex positions
%    F  #F by 3 list of triangle indices
% Output:
%    id  #rays list of indices into F (0 if no hit)
%    t  #rays list of distances from the ray origin (inf if no hit)
%    lambda  #rays by 3 list of barycentric coordinate of hit
%

author: alecjacobson

repo: gptoolbox


text:
function [L, D] = LegendreGaussianIntegral(x, n)
% function [L, D] = legendreGaussianIntegral(x, n)
% Computes legendre gaussian integrals up to the order specified and the
% derivatives if requested
%
% The integral takes the following form, in Mathematica syntax,
%
% L[x, n] = Integrate[Exp[-x \mu^2] Legendre[2*n, \mu], {\mu, -1, 1}]
% D[x, n] = Integrate[Exp[-x \mu^2] (-\mu^2) Legendre[2*n, \mu], {\mu, -1, 1}]
%
% INPUTS:
%
% x should be a column vector of positive numbers, specifying the
% parameters of the gaussian
%
% n should be a non-negative integer, such that 2n specifies the maximum order
% of legendre polynomial
%
% The maximum value for n is 6.
%
% OUTPUTS:
%
% L will be a two-dimensional array with each row containing the
% legendre gaussian integrals of the orders 0, 2, 4, ..., to 2n for the
% parameter value at the corresponding row in x
%
% Note that the legendre gaussian integrals of the odd orders are zero.
%
% D will be the 1st order derivative of L
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

% Make sure n is not larger than 6
if n > 6
	error('The maximum value for n is 6, which corresponds to the 12th order Legendre polynomial');
end

%
% Computing the related exponent gaussian integrals
% I[x, n] = Integrate[Exp[-x \mu^2] \mu^(2*n), {\mu, -1, 1}]
% with the following recursion:
% 1) I[x, 0] = Sqrt[\pi]Erf[x]/Sqrt[x]
% 2) I[x, n+1] = -Exp[-x]/x + (2n+1)/(2x) I[x, n]
%
% This does not work well when x is small

exact = find(x>0.05);
approx = find(x<=0.05);
% Necessary to make matlab happy when x is a single value
exact = exact(:);
approx = approx(:);

if nargout > 1
	mn = n + 2;
else
	mn = n + 1;
end
I = zeros(length(x),mn);
sqrtx = sqrt(x(exact));
I(exact,1) = sqrt(pi)*erf(sqrtx)./sqrtx;
dx = 1.0./x(exact);
emx = -exp(-x(exact));
for i = 2:mn
    I(exact,i) = emx + (i-1.5)*I(exact,i-1);
    I(exact,i) = I(exact,i).*dx;
end

% Computing the legendre gaussian integrals for large enough x
L = zeros(length(x),n+1);
for i = 1:n+1
	if i == 1
		L(exact,1) = I(exact,1);
	elseif i == 2
		L(exact,2) = -0.5*I(exact,1) + 1.5*I(exact,2);
	elseif i == 3
		L(exact,3) = 0.375*I(exact,1) - 3.75*I(exact,2) + 4.375*I(exact,3);
	elseif i == 4
		L(exact,4) = -0.3125*I(exact,1) + 6.5625*I(exact,2) - 19.6875*I(exact,3) + 14.4375*I(exact,4);
	elseif i == 5
		L(exact,5) = 0.2734375*I(exact,1) - 9.84375*I(exact,2) + 54.140625*I(exact,3) - 93.84375*I(exact,4) + 50.2734375*I(exact,5);
    elseif i == 6
        L(exact,6) = -(63/256)*I(exact,1) + (3465/256)*I(exact,2) - (30030/256)*I(exact,3) + (90090/256)*I(exact,4) - (109395/256)*I(exact,5) + (46189/256)*I(exact,6);
    elseif i == 7
        L(exact,7) = (231/1024)*I(exact,1) - (18018/1024)*I(exact,2) + (225225/1024)*I(exact,3) - (1021020/1024)*I(exact,4) + (2078505/1024)*I(exact,5) - (1939938/1024)*I(exact,6) + (676039/1024)*I(exact,7);
	end
end

% Computing the legendre gaussian integrals for small x
x2=x(approx,1).^2;
x3=x2.*x(approx,1);
x4=x3.*x(approx,1);
x5=x4.*x(approx,1);
x6=x5.*x(approx,1);
for i = 1:n+1
	if i == 1
		L(approx,1) = 2 - 2*x(approx,1)/3 + x2/5 - x3/21 + x4/108;
	elseif i == 2
		L(approx,2) = -4*x(approx,1)/15 + 4*x2/35 - 2*x3/63 + 2*x4/297;
	elseif i == 3
		L(approx,3) = 8*x2/315 - 8*x3/693 + 4*x4/1287;
	elseif i == 4
		L(approx,4) = -16*x3/9009 + 16*x4/19305;
	elseif i == 5
		L(approx,5) = 32*x4/328185;
    elseif i == 6
		L(approx,6) = -64*x5/14549535;
    elseif i == 7
		L(approx,7) = 128*x6/760543875;
	end
end

if nargout == 1
	return;
end

% Computing the derivatives for large enough x
D = zeros(length(x),n+1);
for i = 1:n+1
	if i == 1
		D(exact,1) = -I(exact,2);
	elseif i == 2
		D(exact,2) = 0.5*I(exact,2) - 1.5*I(exact,3);
	elseif i == 3
		D(exact,3) = -0.375*I(exact,2) + 3.75*I(exact,3) - 4.375*I(exact,4);
	elseif i == 4
		D(exact,4) = 0.3125*I(exact,2) - 6.5625*I(exact,3) + 19.6875*I(exact,4) - 14.4375*I(exact,5);
	elseif i == 5
		D(exact,5) = -0.2734375*I(exact,2) + 9.84375*I(exact,3) - 54.140625*I(exact,4) + 93.84375*I(exact,5) - 50.2734375*I(exact,6);
    elseif i == 6
        D(exact,6) = (63/256)*I(exact,2) - (3465/256)*I(exact,3) + (30030/256)*I(exact,4) - (90090/256)*I(exact,5) + (109395/256)*I(exact,6) - (46189/256)*I(exact,7);
    elseif i == 7
        D(exact,7) = -(231/1024)*I(exact,2) + (18018/1024)*I(exact,3) - (225225/1024)*I(exact,4) + (1021020/1024)*I(exact,5) - (2078505/1024)*I(exact,6) + (1939938/1024)*I(exact,7) - (676039/1024)*I(exact,8);
	end
end

% Computing the derivatives for small x
for i = 1:n+1
	if i == 1
		D(approx,1) = -2/3 + 2*x(approx,1)/5 - x2/7 + x3/27 - x4/132;
	elseif i == 2
		D(approx,2) = -4/15 + 8*x(approx,1)/35 - 2*x2/21 + 8*x3/297 - 5*x4/858;
	elseif i == 3
		D(approx,3) = 16*x(approx,1)/315 - 8*x2/231 + 16*x3/1287 - 4*x4/1287;
	elseif i == 4
		D(approx,4) = -16*x2/3003 + 64*x3/19305 - 8*x4/7293;
	elseif i == 5
		D(approx,5) = 128*x3/328185 - 32*x4/138567;
    elseif i == 6
        D(approx,6) = -64*x4/2909907 + 128*x5/10140585;
    elseif i == 7
        D(approx,7) = 256*x5/253514625;
	end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
clear all; close all; clear classes; clc;

%% Set flags.
inspect_only = false;

%% Solve the system.
s = 220;
w = 20;
dL = 5;
dl = 1;
[E, H, obj_array, src_array, J] = maxwell_run(...
	'OSC', 1e-9, 200, ...
	'DOM', {'vacuum', 'none', 1.0}, [-310 310; -310 310; 0 dl], [dL dL dl], BC.p, [10*dL 10*dL 0], ...
	'OBJ', ...
		{'vacuum', 'none', 1.0}, Box([-w/2 w/2; -w/2 w/2; 0 dl], dl), ...
		{'Johnson/Au', 'y'}, ...
			PolygonalCylinder(Axis.z, dl, dl/2, [w/2 0; w/2+s*sqrt(3)/2 -s/2; w/2+s*sqrt(3)/2 s/2], dl), ...
			PolygonalCylinder(Axis.z, dl, dl/2, [-w/2 0; -w/2-s*sqrt(3)/2 s/2; -w/2-s*sqrt(3)/2 -s/2], dl), ...
	'SRCJ', PointSrc(Axis.x, [0 0 0]), ...
	inspect_only);

%% Visualize the solution.
figure
clear opts
opts.withinterp = false;
opts.withobjsrc = true;
opts.withabs = true;
% opts.withgrid = true;
opts.cscale = 1e-3;
z_location = 0;
vis2d(E{Axis.x}, Axis.z, z_location, obj_array, src_array, opts)
% vis2d(H{Axis.z}, Axis.z, z_location, obj_array, src_array, opts)

% %% Calculate the power emanating from the source.
% power = powerflux_box(E,H,[-10 10; -10 10; 0 1]);
% fprintf('power = %e\n', power);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [MStitch, result] = GrayListMain_Process(file)
for i = 1 : length(file)-1
    if i == 1
        file1 = file{i};
        im1 = imread(file1);
        im1 = rgb2gray(im1);
        MStitch.im1 = double(im1);

        [Pheight, Pwidth] = size(im1);

        MStitch.Pwidth = Pwidth;
        MStitch.Pheight = Pheight;

        MStitch.W_min = round(0.60*Pwidth);
        MStitch.W_max = round(0.83*Pwidth);
        MStitch.H_min = round(0.98*Pheight);
        MStitch.minval = 255;

        im1 = imread(file1);
        MStitch.imrgb1 = double(im1);
        im1 = rgb2gray(im1);
        MStitch.im1 = double(im1);
    else
        MStitch.im1 = double(result1);
    end
    file2 = file{i+1};
    im2 = imread(file2);
    MStitch.imrgb2 = double(im2);
    im2 = rgb2gray(im2);
    im2 = double(im2);
    [W_box, H_box, bdown, MStitch] = Fun_Match(im2, MStitch);
    [MStitch, result1] = Fun_Stitch(im2, W_box, H_box, bdown, MStitch);
end
result = result1;

author: decouples


text:
function out = sum(f)
%SUM   Definite integral of a SINGFUN on the interval [-1,1].
%   SUM(F) is the integral of F from -1 to 1.
%
% See also CUMSUM, DIFF.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Note for developers:
%
% The main algorithm:
%
% When the smoothPart of F is a CHEBTECH, that is, it can be written as
% Chebyshev sum, the integral
%
% I = \int_{-1}^{1} F dx = \sum_{0}^{n-1} c_r M_r,
%
% where M_r = \int_{-1}^{1} (1+x)^a(1-x)^b T_r(x) dx is the rth Jacobi moment.
%
% The computation of M_r is treated differently for different a and b:
%
% (I) when a == b, M_r are the Gegenbauer moments, which bear a closed-form
% solution as indicated in [2] and [4].
%
% (II) when a ~= b, M_r are the general Jacobi moments, which can be obtained
% using a three-term recursive relation discussed in [3].
%
% This way, all quadratures in SINGFUN, along with those in CHEBTECH are now
% entirely of Clenshaw-Curtis style.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Useful References:
%
% [1]. K. Xu and M. Javed, Singfun Working Note, August 2013
%
% [2]. Hunter, D., and Nikolov, G., Gaussian Quadrature of Chebyshev
% Polynomials, J. Comput. Appl. Math. 94 (1998), 123-131.
%
% [3]. Piessens, R., and Branders, M., The Evaluation and Application of Some
% Modified Moments, BIT 13 (1973), 443-450.
%
% [4]. Sommariva, A., Fast construction of Fejer and Clenshaw–Curtis rules for
% general weight functions, Computers & Mathematics with Applications 65
% (2012), 682-693.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
% Trivial cases:
if ( all(f.exponents == 0) )
    % If both the exponents are trivial, then compute the integral by calling
    % the sum in smoothfun.
    out = sum(f.smoothPart);
    return

elseif ( all(f.exponents <= -1) )
    % The integral is divergent or not a number when both the exponents are
    % non-zero.

    if ( isreal(f) )

        % The sign of the smoothPart of F at the end points:
        sl = sign(get(f.smoothPart, 'lval'));
        sr = sign(get(f.smoothPart, 'rval'));

        if  ( sl == sr )
            out = sl.*inf;
        else
            out = NaN;
        end

    else
        % The sign of the real part of the smoothPart of F at the end
        % points:
        rsl = sign(real(get(f.smoothPart, 'lval')));
        rsr = sign(real(get(f.smoothPart, 'rval')));

        % The sign of the real part of the smoothPart of F at the end
        % points:
        isl = sign(imag(get(f.smoothPart, 'lval')));
        isr = sign(imag(get(f.smoothPart, 'rval')));

        if ( (rsl == rsr) && (isl == isr) )
            out = Inf + 1i*Inf;
        else
            out = NaN;
        end
    end

    return

elseif ( any(f.exponents <= -1) )

    % The integral is divergent or not a number when one of the exponents
    % are non-zero.

    if ( isreal(f) )

        % The sign of the smoothPart of F at the end points:
        sl = sign(get(f.smoothPart, 'lval'));
        sr = sign(get(f.smoothPart, 'rval'));

        s = [sl sr];
        ind = ( f.exponents <= -1 );
        out = Inf*s(ind);

    else
        % The sign of the real part of the smoothPart of F at the end
        % points:
        rsl = sign(real(get(f.smoothPart, 'lval')));
        rsr = sign(real(get(f.smoothPart, 'rval')));

        % The sign of the real part of the smoothPart of F at the end
        % points:
        isl = sign(imag(get(f.smoothPart, 'lval')));
        isr = sign(imag(get(f.smoothPart, 'rval')));

        rs = [rsl rsr];
        is = [isl isr];
        ind = ( f.exponents <= -1 );

        % The real part:
        realPart = 0;
        if ( any(rs) )
            realPart = Inf*rs(ind);
        end

        % The imaginary part:
        imagPart = 0;
        if ( any(is) )
            imagPart = Inf*is(ind);
        end

        % The sum:
        out = realPart + imagPart;
    end

    return
end

%%
% The non-trivial case:

if ( isa(f.smoothPart, 'chebtech') )

    % If the smooth part of F is a CHEBTECH, then evaluate the integral by using
    % Clenshaw-Curtis-Jacobi moments and the Chebyshev coefficients:

    % Grab the number of points for the smooth part of f:
    n = length(f);

    % Grab the exponents:
    a = f.exponents(1);

    if ( diff(f.exponents) == 0 )
        % If the exponents at the endpoints are same, then compute the
        % appropriate modified moments for Gegenbauer weights.

        r = a + .5;
        m0 = gamma(r + .5)*sqrt(pi)/gamma(r + 1);
        k = 1:floor((n-1)/2);
        % Even modified moments for M_2k = \int_{-1}^1 (1-x)^a(1+x)^a T_2k(x) dx
        % and notice that the odd moments vanish due to parity.
        m = m0*[1, cumprod((k - r - 1)./(k + r))];
        % Form the modified moments vector:
        M(1:2:n) = m;
        M(2:2:n) = 0;

    else

        % The general case:
        b = f.exponents(2);

        % Common coefficient for the modified moments:
        c1 = a + 1;
        c2 = b + 1;
        c3 = a + b + 1;
        c4 = c1 + c2;
        c5 = a - b;
        c0 = (2^c3)*beta(c1, c2);

        % Compute the hypergeometric function related to the modified moments:
        M = zeros(1,n);
        M(1) = 1;
        if ( n > 1 )
            M(2) = c5/c4;
            % Sister Celine's three-term recurrence:
            for j = 3:n
                M(j) = (2*c5*M(j-1) + (j - 2 - c4)*M(j-2)) / (c3 + j - 1);
            end
        end
        % Compute the modified moments:
        M = c0*M;

    end

    % Chebyshev coefficients of the smooth part of F:
    coeffs = get(f, 'coeffs');

    % multiplication of weights and values
    out = M*coeffs;

else
    %%
    % If f.smoothPart is not a CHEBTECH, we evaluate the integral by using
    % Gauss-Jacobi points and weights.

    % Give a sufficiently large number:
    % [TODO]: This number needs to be determined in future when other 'tech's join.
    % [TODO]: Or perhaps compute iteratively until result doesn't change?
    n = 1000;

    [x, w] = jacpts(ceil(n/2) + 1, f.exponents(2), f.exponents(1));
    out = w*feval(f.smoothPart, x);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function rte = pred2path(P,s,t)
%PRED2PATH Convert predecessor indices to shortest paths from 's' to 't'.
%   rte = pred2path(P,s,t)
%     P = |s| x n matrix of predecessor indices (from DIJK)
%     s = FROM node indices
%       = [] (default), paths from all nodes
%     t = TO node indices
%       = [] (default), paths to all nodes
%   rte = |s| x |t| cell array of paths (or routes) from 's' to 't', where
%         rte{i,j} = path from s(i) to t(j)
%                  = [], if no path exists from s(i) to t(j)
%
% (Used with output of DIJK)

% Copyright (c) 1994-2006 by Michael G. Kay
% Matlog Version 9 13-Jan-2006 (http://www.ie.ncsu.edu/kay/matlog)

% Input Error Checking ****************************************************
narginchk(1,3);

[~,n] = size(P);

if nargin < 2 || isempty(s), s = (1:n)'; else, s = s(:); end
if nargin < 3 || isempty(t), t = (1:n)'; else, t = t(:); end

if any(P < 0 | P > n)
   error(['Elements of P must be integers between 1 and ',num2str(n)]);
elseif any(s < 1 | s > n)
   error(['"s" must be an integer between 1 and ',num2str(n)]);
elseif any(t < 1 | t > n)
   error(['"t" must be an integer between 1 and ',num2str(n)]);
end
% End (Input Error Checking) **********************************************

rte = cell(length(s),length(t));

[~,idxs] = find(P==0);

for i = 1:length(s)
%    if rP == 1
%       si = 1;
%    else
%       si = s(i);
%       if si < 1 | si > rP
%          error('Invalid P matrix.')
%       end
%    end
   si = find(idxs == s(i));
   for j = 1:length(t)
      tj = t(j);
      if tj == s(i)
         r = tj;
      elseif P(si,tj) == 0
         r = [];
      else
         r = tj;
         while tj ~= 0
            if tj < 1 || tj > n
               error('Invalid element of P matrix found.')
            end
            r = [P(si,tj) r];
            tj = P(si,tj);
         end
         r(1) = [];
      end
      rte{i,j} = r;
   end
end

if length(s) == 1 && length(t) == 1
   rte = rte{:};
end

%rte = t;
while 0%t ~= s
   if t < 1 || t > n || round(t) ~= t
      error('Invalid "pred" element found prior to reaching "s"');
   end
   rte = [P(t) rte];
   t = P(t);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function polygonal_surface_display ( node_file_name, face_file_name )

%*****************************************************************************80
%
%% POLYGONAL_SURFACE_DISPLAY displays a polygonal surface file.
%
%  Discussion:
%
%    This program reads two files defining the nodes and the faces
%    of a polygonal surface in 3D.
%
%    It displays the surface using MATLAB's graphics system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.
%
%    Input, string FACE_FILE_NAME, the name of the face file.
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGONAL_SURFACE_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the nodes and faces defining\n' );
  fprintf ( 1, '  a polygonal surface and display the object\n' );
  fprintf ( 1, '  as a MATLAB shape.\n' );
%
%  If at least one command line argument, it's the node file name.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGONAL_SURFACE_DISPLAY:\n' );
    node_file_name = input ( 'Enter the name of the node file:' );

  end
%
%  If at least two command line arguments, it's the face file name.
%
  if ( nargin < 2 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGONAL_SURFACE_DISPLAY:\n' );
    face_file_name = input ( 'Enter the name of the face file:' );

  end
%
%  Read the data.
%
  [ node_num, node_xyz ] = polygonal_surface_node_file_read ( node_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  r8mat_transpose_print_some ( 3, node_num, node_xyz, 1, 1, 3, 5, ...
    '  First five nodes:' );

  [ face_order, face_num, face_node ] = ...
    polygonal_surface_face_file_read ( face_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Face order =      %d\n', face_order );
  fprintf ( 1, '  Number of faces = %d\n', face_num );

  i4mat_transpose_print_some ( face_order, face_num, face_node, ...
    1, 1, face_order, 5, '  First five faces:' );
%
%  Display the image as a collection of polygons.
%
  handle = patch ( 'Vertices', node_xyz', 'Faces', face_node' );

  set ( handle, 'FaceColor', [0.5, 0.6, 0.8], 'EdgeColor', 'Black' );

  axis equal;
  grid on;

  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
  zlabel ( '--Z axis--' )
%
%  The TITLE function will interpret underscores in the title.
%  We need to unescape such escape sequences!
%
  title_string = s_escape_tex ( face_file_name );
  title ( title_string )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGONAL_SURFACE_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 10;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function [ face_order, face_num, face_node ] = ...
  polygonal_surface_face_file_read ( face_file_name )

%*****************************************************************************80
%
%% POLYGONAL_SURFACE_FACE_FILE_READ reads faces defining a surface.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Angel,
%    Interactive Computer Graphics,
%    a Top-Down Approach with OpenGL,
%    Addison-Wesley, 2000.
%
%  Parameters:
%
%    Input, string FACE_FILE_NAME, the name of the face file.
%
%    Output, integer FACE_ORDER, the order of the face.
%
%    Output, integer FACE_NUM, the number of face defined.
%
%    Output, real FACE_NODE(FACE_ORDER,FACE_NUM),
%    the nodes making up each face.
%
  [ face_order, face_num ] = i4mat_header_read ( face_file_name );

  face_node = i4mat_data_read ( face_file_name, face_order, face_num );

  return
end
function [ node_num, node_xyz ] = polygonal_surface_node_file_read ( ...
  node_file_name )

%*****************************************************************************80
%
%% POLYGONAL_SURFACE_NODE_FILE_READ reads node information from a polygonal surface node file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Angel,
%    Interactive Computer Graphics,
%    a Top-Down Approach with OpenGL,
%    Addison-Wesley, 2000.
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.
%
%    Output, integer NODE_NUM, the number of nodes defined.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the coordinates of nodes.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_file_name );

  node_xyz = r8mat_data_read ( node_file_name, dim_num, node_num );

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function asset_path_test ( )

%*****************************************************************************80
%
%% ASSET_PATH_TEST tests ASSET_PATH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASSET_PATH_TEST:\n' );
  fprintf ( 1, '  Demonstrate the simulation of an asset price path.\n' );

  s0 = 2.0;
  mu = 0.1;
  sigma = 0.3;
  n = 100;
  t0 = 0.0;
  t1 = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The asset price at time 0,     S0    = %f\n', s0 );
  fprintf ( 1, '  The asset expected growth rate MU    = %f\n', mu );
  fprintf ( 1, '  The asset volatility           SIGMA = %f\n', sigma );
  fprintf ( 1, '  The expiry date                T1    = %f\n', t1 );
  fprintf ( 1, '  The number of time steps       N     = %d\n', n );

  s = asset_path ( s0, mu, sigma, t1, n );
%
%  Plot.
%
  figure ( 1 )
  t = ( linspace ( t0, t1, n + 1 ) )';
  plot ( t, s )
  grid on
  xlabel ( '<-- Time -->' )
  ylabel ( '<-- Value --> ' )
  title ( 'Simulated asset path' )
%
%  Print a little.
%
  r8vec_print_part ( n + 1, s, 10, '  Partial results:' );
%
%  Write to a file.
%
  output_filename = 'asset_path.txt';

  r8vec_write ( output_filename, n + 1, s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Full results written to "%s".\n', output_filename );

  return
end


text:
function parPath = er_deconvolvedParfile(view, stim, params);
%
% parPath = er_deconvolvedParfile(view, <stim>, <params>);
%
% Create a parfile representing the deconvolved time courses
% from a GLM. The format the data are stored in is alternating
% beta values and residual variance from the GLM, each taking
% [nh] time points, where nh is the number of temporal frames in the
% estimated hemodynamic response. (This should correspond to the number of
% MR frames contained within the event-related 'timeWindow' parameter.
%
% stim: stim/'trials' struct loaded from er_concatParfiles for the
%       source scans from the deconvolution.
% params: event-related parameters from the source scans for the
%       deconvolution.
% <default: get stim / params from the view's current scan group.>
%
% ras, 02/2006.
if notDefined('view'),      view = getSelectedInplane;          end
if notDefined('stim'),      stim = er_concatParfiles(view);     end
if notDefined('params'),    params = er_getParams(view);        end

tr = params.framePeriod;
frameWindow = unique(round(params.timeWindow/tr));
nh = length(frameWindow);
nConds = sum(stim.condNums>=1);

% construct the par struct
par.cond = [1:nConds 0];
par.onset = 0:2*nh:2*nh*nConds;
par.label = stim.condNames([2:nConds+1 1]);
par.color = stim.condColors([2:nConds+1 1]);

% write out the parfile
mrGlobals
dt = existDataType('Deconvolved');
if dt==0,
    scan = 1;
else
    scan = length(dataTYPES(dt).scanParams); % assume last deconvolved scan
end
parFileName = sprintf('deconvolved_scan%i.par', scan);
parPath = fullfile(parfilesDir(view), parFileName);
writeParfile(par, parPath);

return


text:
function err = gzSeek(gzFp, skipBytes, origin)
%
% err = gzSeek(gzFp, skipBytes, origin)
%
% HISTORY:
% 2007.02.12 RFD: wrote it.

if isempty(gzFp) || ~isfield(gzFp,'inStream') || isempty(gzFp.inStream)
    disp('Can''t read file! Is it open?');
    err = -1;
    return;
end

if(~exist('origin','var')||isempty(origin))
    origin = 0;
end

if((isnumeric(origin)&&origin==-1) || strcmpi(origin,'bof'))
    if(gzFp.inStream.markSupported)
        gzFp.inStream.reset;
    elseif(skipBytes>=gzFp.curPos)
        % If we just want to go forward, then we can do that
        skipBytes = skipBytes - gzFp.curPos;
    else
        % Reset isn't supported
        disp('reset not supported- can''t skip from beginning.');
        err = -1;
        return;
    end
elseif((isnumeric(origin)&&origin==1) || strcmpi(origin,'eof'))
    disp('''eof'' option not supported.');
    err = -1;
    return;
end

if(skipBytes>0)
    gzFp.inStream.skip(int32(skipBytes));
    gzFp.curPos = gzFp.curPos+skipBytes;
elseif(skipBytes<0)
    disp('negative skip not supported.');
    err = -1;
    return;
end

return;


text:
% TERNAXES create ternary axis
%   HOLD_STATE = TERNAXES(MAJORS) creates a ternary axis system using the system
%   defaults and with MAJORS major tickmarks.

% Author: Carl Sandrock 20050211

% To Do

% Modifications

% Modifiers
% (CS) Carl Sandrock

function [hold_state, cax, next] = ternaxes(majors)

%TODO: Get a better way of offsetting the labels
xoffset = 0.04;
yoffset = 0.02;

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
ls = get(cax,'gridlinestyle');

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');

set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
    'DefaultTextFontName',   get(cax, 'FontName'), ...
    'DefaultTextFontSize',   get(cax, 'FontSize'), ...
    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
    'DefaultTextUnits','data')

% only do grids if hold is off
if ~hold_state
	%plot axis lines
	hold on;
	plot ([0 1 0.5 0],[0 0 sin(1/3*pi) 0], 'color', tc, 'linewidth',1,...
                   'handlevisibility','off');
	set(gca, 'visible', 'off');

    % plot background if necessary
    if ~isstr(get(cax,'color')),
       patch('xdata', [0 1 0.5 0], 'ydata', [0 0 sin(1/3*pi) 0], ...
             'edgecolor',tc,'facecolor',get(gca,'color'),...
             'handlevisibility','off');
    end

	% Generate labels
	majorticks = linspace(0, 1, majors + 1);
	majorticks = majorticks(1:end-1);
	labels = num2str(majorticks'*100);

    zerocomp = zeros(size(majorticks)); % represents zero composition

	% Plot right labels (no c - only b a)
	[lxc, lyc] = terncoords(1-majorticks, majorticks, zerocomp);
	text(lxc, lyc, [repmat('  ', length(labels), 1) labels]);

	% Plot bottom labels (no b - only a c)
	[lxb, lyb] = terncoords(majorticks, zerocomp, 1-majorticks); % fB = 1-fA
	text(lxb, lyb, labels, 'VerticalAlignment', 'Top');

	% Plot left labels (no a, only c b)
	[lxa, lya] = terncoords(zerocomp, 1-majorticks, majorticks);
	text(lxa-xoffset, lya, labels);

	nlabels = length(labels)-1;
	for i = 1:nlabels
        plot([lxa(i+1) lxb(nlabels - i + 2)], [lya(i+1) lyb(nlabels - i + 2)], ls, 'color', tc, 'linewidth',1,...
           'handlevisibility','off');
        plot([lxb(i+1) lxc(nlabels - i + 2)], [lyb(i+1) lyc(nlabels - i + 2)], ls, 'color', tc, 'linewidth',1,...
           'handlevisibility','off');
        plot([lxc(i+1) lxa(nlabels - i + 2)], [lyc(i+1) lya(nlabels - i + 2)], ls, 'color', tc, 'linewidth',1,...
           'handlevisibility','off');
	end;
end;

% Reset defaults
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits', fUnits );


text:
% SPHERE_READ Read a NIST SPHERE audio file
%
% Usages
%    [y, fs] = SPHERE_READ(filename, N)
%
%    sz = SPHERE_READ(filename, 'size')
%
% Input
%    filename (char): The name of the file to read.
%    N (int, optional): The number of samples to read.
%
% Output
%    y (numeric): The audio data in the file.
%    fs (numeric): The sample rate of the data.
%    sz (int): The number of samples in the file.
%
% Description
%    Analogous in use to the WAVREAD function, SPHERE_READ processes NIST
%    SPHERE files, such as those found in TIMIT. For more information, see
%    WAVREAD.
%
% See also
%    WAVREAD

function [y, fs] = sphere_read(filename, N)
	if nargin < 2
		N = [];
	end

	fid = fopen(filename);

	if fid == -1
		error('File not found.');
	end

	line = fgets(fid);

	if (isnumeric(line) && line == -1) || isempty(strfind(line, 'NIST_1A'))
		error('Invalid NIST SPHERE file (header not present).');
	end

	line = fgets(fid);

	header_length = str2num(line);

	if isempty(header_length)
		error('Invalid NIST SPHERE file (header length not present).');
	end

	header_length = round(header_length);

	fclose(fid);

	fid = fopen(filename,'rt','ieee-le');

	[data, count] = fread(fid, header_length, 'uchar');

	fclose(fid);

	if count ~= header_length
		error('Invalid NIST SPHERE file (unable to read header).');
	end

	header = char(data');

	lines = regexp(header,char(10),'split');

	fs = [];
	byte_format = [];
	count = [];

	for k = 1:length(lines)
		if strfind(lines{k}, 'sample_coding') == 1
			if ~isempty(strfind(lines{k}, 'shorten'))
				error('Unsupported coding scheme (''shorten'').');
			end
		elseif strfind(lines{k}, 'sample_rate') == 1
			fs = strread(lines{k}, 'sample_rate -i %d');
		elseif strfind(lines{k}, 'sample_count') == 1
			count = strread(lines{k}, 'sample_count -i %d');
		elseif strfind(lines{k}, 'channel_count') == 1
			if strread(lines{k}, 'channel_count -i %d') > 1
				error('Unsupported number of channels (>1).');
			end
		elseif strfind(lines{k}, 'sample_n_bytes') == 1
			if strread(lines{k}, 'sample_n_bytes -i %d') ~= 2
				error('Unsupported bit depth (not 16 bits).');
			end
		elseif strfind(lines{k}, 'sample_byte_format') == 1
			byte_format = strread(lines{k}, 'sample_byte_format -s2 %s');
		end
	end

	if ischar(N) && strcmp(N,'size')
		y = count;
		return;
	end

	if ~isempty(N) && isnumeric(N)
		count = N;
	end

	if isempty(byte_format)
		warning('''sample_byte_format'' not specified, assuming little-endian.');
		byte_format = '01';
	end

	if strcmp(byte_format, '01')
		endianness = 'le';
	elseif strcmp(byte_format, '10')
		endianness = 'be';
	else
		error('Unknown ''sample_byte_format''.');
	end

	fid = fopen(filename, 'r', ['ieee-' endianness]);

	temp = fread(fid, header_length, 'char');

	if isempty(count)
		count = Inf;
	end

	y = fread(fid, count, 'short');

	fclose(fid);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function setSrcComboUserData(hObject,handles,plotFrame)
    contents = cellstr(get(hObject,'String'));
    sel = strtrim(contents{get(hObject,'Value')});
    set(hObject,'UserData',sel);

    if(plotFrame && get(handles.playButton,'Value')==0)
        plotOneFrame(handles);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Whitener = bst_whitener(NoiseCov, ChannelFile, DataTypes, ChannelFlag)
% BST_WHITENER: Compute a whitener from a NoiseCov matrix
%
% USAGE:  Whitener = bst_whitener(NoiseCov);
%         Whitener = bst_whitener(NoiseCov, ChannelFile, DataTypes, ChannelFlag);
%         Whitener = bst_whitener(NoiseCov,     Channel, DataTypes, ChannelFlag);

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

% If channels definition is not provided
if (nargin < 4)
    % Use all channels
    iChan = 1:length(NoiseCov);
else
    % If channel file provided: load ChannelFile
    if ischar(ChannelFile)
        ChannelMat = in_bst_channel(ChannelFile, 'Channel');
        Channel = ChannelMat.Channel;
    % If Channel structure provided
    elseif isstruct(ChannelFile)
        Channel = ChannelFile;
    end
    % Get good channels
    iChan = good_channel(Channel, ChannelFlag, DataTypes);
end

% Detect the rows with only zero values
iZeroRow = find(sum(NoiseCov .^ 2) == 0);
% Remove those channels from list of valid channels
iChan = setdiff(iChan, iZeroRow);

% Initialize output matrix
Whitener = zeros(size(NoiseCov));
% Decomposition
[U,S] = svd(NoiseCov(iChan,iChan));
% Check matrix rank
m = length(iChan);
r = sum(diag(S) > m * S(1) * eps('double'));
if (r < m)
    error(['You have deficient data. Please:' 10 10 ...
           '1) Look for bad channels in your recordings' 10 ...
           '2) Tag them bad channels in all your recordings' 10 ...
           '3) Try automatic detection of flat channels: Right click > Good/bad channels > ...' 10 ...
           '4) Recompute the noise covariance matrix' 10, ...
           '5) Restart this proces']);
end
% Create whitener
Whitener(iChan,iChan) = pinv(U * sqrt(S));

author: brainstorm-tools


text:
function bnet = mk_qmr_bnet(G, inhibit, leak, prior, tabular_findings, onodes)
% MK_QMR_BNET Make a QMR model
% bnet = mk_qmr_bnet(G, inhibit, leak, prior)
%
% G(i,j) = 1 iff there is an arc from disease i to finding j
% inhibit(i,j) = inhibition probability on i->j arc
% leak(j) = inhibition prob. on leak->j arc
% prior(i) = prob. disease i is on
% tabular_findings = 1 means multinomial leaves (ignores leak/inhibit params)
%   = 0 means noisy-OR leaves (default = 0)

if nargin < 5, tabular_findings = 0; end

[Ndiseases Nfindings] = size(inhibit);
N = Ndiseases + Nfindings;
finding_node = Ndiseases+1:N;
ns = 2*ones(1,N);
dag = zeros(N,N);
dag(1:Ndiseases, finding_node) = G;
if nargin < 6, onodes = finding_node; end
bnet = mk_bnet(dag, ns, 'observed', onodes);

for d=1:Ndiseases
    CPT = [1-prior(d) prior(d)];
    bnet.CPD{d} = tabular_CPD(bnet, d, CPT');
end

for i=1:Nfindings
    fnode = finding_node(i);
    ps = parents(G, i);
    if tabular_findings
        bnet.CPD{fnode} = tabular_CPD(bnet, fnode);
    else
        bnet.CPD{fnode} = noisyor_CPD(bnet, fnode, leak(i), inhibit(ps, i));
    end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% Histogram Comparison
%
% In this demo, we show how to:
%
% * Use the function |cv.compareHist| to get a numerical parameter that
%   express how well two histograms match with each other
% * Use different metrics to compare histograms
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d8/dc8/tutorial_histogram_comparison.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/Histograms_Matching/compareHist_Demo.cpp>
%

%% Theory
%
% To compare two histograms ($H_{1}$ and $H_{2}$), first we have to choose a
% _metric_ ($d(H_{1}, H_{2})$) to express how well both histograms match.
%
% OpenCV implements the function |cv.compareHist| to perform a comparison. It
% also offers 4 different metrics to compute the matching:
%
% * *Correlation*:
%
% $$d(H_1,H_2) =  \frac{\sum_I (H_1(I) - \bar{H_1}) (H_2(I) - \bar{H_2})}
%                      {\sqrt{\sum_I(H_1(I) - \bar{H_1})^2
%                             \sum_I(H_2(I) - \bar{H_2})^2}}$$
%
% where
%
% $$\bar{H_k} =  \frac{1}{N} \sum _J H_k(J)$$
%
% and $N$ is the total number of histogram bins.
%
% * *Chi-Square*:
%
% $$d(H_1,H_2) =  \sum _I  \frac{\left(H_1(I)-H_2(I)\right)^2}{H_1(I)}$$
%
% * *Intersection*:
%
% $$d(H_1,H_2) =  \sum _I  \min (H_1(I), H_2(I))$$
%
% * *Bhattacharyya*:
%
% $$d(H_1,H_2) =  \sqrt{1 - \frac{1}{\sqrt{\bar{H_1} \bar{H_2} N^2}}
%                           \sum_I \sqrt{H_1(I) \cdot H_2(I)}}$$
%

%% Code
%
% This program:
%
% * Loads a _base image_ and 2 _test images_ to be compared with it.
% * Generate 1 image that is the lower half of the _base image_
% * Convert the images to HSV format
% * Calculate the H-S histogram for all the images and normalize them in order
%   to compare them.
% * Compare the histogram of the _base image_ with respect to the 2 test
%   histograms, the histogram of the lower half base image and with the same
%   base image histogram.
% * Display the numerical matching parameters obtained.
%

%%
% Load source images (base image and the two other images to compare)
im = {
    'https://docs.opencv.org/3.3.1/Histogram_Comparison_Source_0.jpg'
    'https://docs.opencv.org/3.3.1/Histogram_Comparison_Source_1.jpg'
    'https://docs.opencv.org/3.3.1/Histogram_Comparison_Source_2.jpg'
};
src = cell(3,1);
for i=1:3
    [~,name,ext] = fileparts(im{i});
    fname = fullfile(mexopencv.root(), 'test', [name ext]);
    if exist(fname, 'file') ~= 2
        disp('Downloading image...')
        urlwrite(im{i}, fname);
    end
    src{i} = cv.imread(fname, 'Color',true);
end

%%
% also create an image of half the base image
src{4} = src{1}(end/2:end,:,:);
src = src([1 4 2 3]);

%%
% show images
names = {'Base', 'Half', 'Test1', 'Test2'};
for i=1:numel(src)
    subplot(2,2,i), imshow(src{i}), title(names{i})
end

%%
% convert images to HSV color space
hsv = cell(size(src));
for i=1:numel(src)
    hsv{i} = cv.cvtColor(src{i}, 'RGB2HSV');
end

%%
% calculate H-S 2D histograms
%ranges = {linspace(0,180,50+1), linspace(0,256,60+1)};
ranges = {[0 180], [0 256]};
hsizes = [50, 60];
histo = cell(size(hsv));
for i=1:numel(hsv)
    histo{i} = cv.calcHist(hsv{i}(:,:,1:2), ranges, 'HistSize',hsizes, 'Uniform',true);
    histo{i} = cv.normalize(histo{i}, 'NormType','MinMax', 'Alpha',0, 'Beta',1);
end

%%
% compare histogram of the base image against the other histograms
% using four different metrics
algs = {'Correlation', 'ChiSquare', 'Intersection', 'Bhattacharyya'};
D = zeros(numel(algs),numel(histo));
for j=1:numel(histo)
    for i=1:numel(algs)
        D(i,j) = cv.compareHist(histo{1}, histo{j}, 'Method',algs{i});
    end
end

%%
% The first image is the base (to be compared to the others), the other 2 are
% the test images. We also compare the first image with respect to itself and
% with respect of half the base image.
%
% We should expect a perfect match when we compare the base image histogram
% with itself. Also, compared with the histogram of half the base image, it
% should present a high match since both are from the same source. For the
% other two test images, we can observe that they have very different lighting
% conditions, so the matching should not be very good.
%

%%
% Here are the numeric results
if ~mexopencv.isOctave()
    t = array2table(D, 'VariableNames',names, 'RowNames',algs);
    disp(t);
else
    %HACK: use cell array instead of table
    t = [{''}, names; algs(:), arrayfun(@num2str,D,'UniformOutput',false)];
    t = t';
    fprintf('%13s %9s %9s %9s %9s\n',t{:});
end

%%
% For the _Correlation_ and _Intersection_ methods, the higher the metric, the
% more accurate the match. As we can see, the match _base-base_ is the highest
% of all as expected. Also we can observe that the match _base-half_ is the
% second best match (as we predicted). For the other two metrics, the less the
% result, the better the match. We can observe that the matches between the
% test 1 and test 2 with respect to the base are worse, which again, was
% expected.
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function quadrule_test10 ( )

%*****************************************************************************80
%
%% TEST10 tests FEJER2_RULE_COMPUTE and FEJER2_RULE_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10\n' );
  fprintf ( 1, '  FEJER2_RULE_COMPUTE computes\n' );
  fprintf ( 1, '  a Fejer type 2 quadrature rule.\n' );
  fprintf ( 1, '  FEJER2_RULE_SET sets\n' );
  fprintf ( 1, '  a Fejer type 2 quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compare:\n' );
  fprintf ( 1, '    (X1,W1) from FEJER2_RULE_SET\n' );
  fprintf ( 1, '    (X2,W2) from FEJER2_RULE_COMPUTE\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '     Order        W1            W2            X1            X2\n' );
  fprintf ( 1, '\n' );

  for order = 1 : 3 : 10

    [ x1, w1 ] = fejer2_rule_set ( order );
    [ x2, w2 ] = fejer2_rule_compute ( order );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %8d\n', order );

    for i = 1 : order
      fprintf ( 1, '            %12f  %12f  %12f  %12f\n', ...
        w1(i), w2(i), x1(i), x2(i) );
    end

  end

  return
end


text:
function [grad_X, grad_W, grad_b] = B_weighting(prev_layer, curr_layer, future_layers)

input = prev_layer{1}.a;
weight = curr_layer.W;
nWeight = length(weight);
D = size(input,1)/nWeight;

grad_X = 0; grad_W=0; grad_b=0;
future_grad = GetFutureGrad(future_layers, curr_layer);
if 0
    for i=1:length(future_layers)
        if curr_layer.update
            for j=1:size(input,2)
                curr_input = reshape(input(:,j), D, nWeight);
                grad_W = grad_W + curr_input' * future_layers{i}.grad(:,j);
            end
            if isfield(curr_layer, 'updateBias')== 0 || curr_layer.updateBias==1
                grad_b = grad_b + sum(future_layers{i}.grad,2);
            end
        end
        for j=1:size(input,2)
            grad_X = grad_X + future_layers{i}.grad(:,j)*weight';
        end
    end
else
    if curr_layer.update
        for j=1:size(input,2)
            curr_input = reshape(input(:,j), D, nWeight);
            grad_W = grad_W + curr_input' * future_grad(:,j);
        end
        if isfield(curr_layer, 'updateBias')== 0 || curr_layer.updateBias==1
            grad_b = grad_b + sum(future_grad,2);
        end
    end
    for j=1:size(input,2)
        grad_X = grad_X + future_grad(:,j)*weight';
    end
end
grad_X = reshape(grad_X, prod(size(grad_X)), 1);

end

author: singaxiong

repo: SignalGraph


text:
function [recall, rankloss]= testNet(db, net, opts, ID, qFeat, dbFeat)

    relja_display('testNet: %s %s', opts.sessionID, ID);

    rankloss= testCoreRank(db, qFeat, dbFeat, opts.margin, opts.nNegChoice, 'nTestSample', opts.nTestRankSample);
    recall= testCore(db, qFeat, dbFeat, 'nTestSample', opts.nTestSample, 'recallNs', opts.recallNs);
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function [normal, d] = plot_line(p1, p2, p3)

% This function plots a line from three points.
% I/P arguments:
%   p1, p2, p3 eg, p1 = [x y z]
%
%
% O/P is:
% normal: it contains a,b,c coeff , normal = [a b c]
% d : coeff

normal = cross(p1 - p2, p1 - p3);
d = p1(1)*normal(1) + p1(2)*normal(2) + p1(3)*normal(3);
d = -d;
x = -100:100; y = -100:100;
[X,Y] = meshgrid(x,y);
Z = (-d - (normal(1)*X) - (normal(2)*Y))/normal(3);
mesh(X,Y,Z)

author: Sable


text:
file = 'PnP_mex.cpp';

INCLUDE = {'/usr/local/include' , '/usr/include/glog', '/usr/include/eigen3',...
	   '/home/htaira/Desktop/zips/ceres-solver/build/config/'};

LIBFOLDER = {'/home/htaira/Desktop/zips/ceres-solver/build/lib/', '/usr/lib/x86_64-linux-gnu/'};

LIBNAME = {'ceres', 'ceres', 'glog', 'gflags', 'suitesparseconfig',...
	   'cholmod', 'lapack', 'blas', 'spqr', 'gomp', 'cxsparse'};

include = strjoin(strcat('-I', INCLUDE));
libfolder = strjoin(strcat('-L', LIBFOLDER));
libname = strjoin(strcat('-l', LIBNAME));

eval(['mex' ' ' file ' ' include ' '  libfolder ' ' libname]);

author: HajimeTaira


text:
%FILLCONVEXPOLY  Fills a convex polygon
%
%     img = cv.fillConvexPoly(img, pts)
%     [...] = cv.fillConvexPoly(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __img__ Image.
% * __pts__ Polygon vertices, stored in numeric array (Nx2/Nx1x2/1xNx2) or
%   cell-array of 2-element vectors `{[x,y], [x,y], ...}`. Supports integer
%   `int32` 2D points.
%
% ## Output
% * __img__ Output image, same size and type as input `img`.
%
% ## Options
% * __Color__ 3-element floating-point vector specifying polygon color.
%   default zeros
% * __LineType__ Type of the polygon boundaries. One of:
%   * __4__ 4-connected line
%   * __8__ 8-connected line (default)
%   * __AA__ anti-aliased line
% * __Shift__ Number of fractional bits in the vertex coordinates. default 0
%
% The function cv.fillConvexPoly draws a filled convex polygon. This function
% is much faster than the function cv.fillPoly. It can fill not only convex
% polygons but any monotonic polygon without self-intersections, that is, a
% polygon whose contour intersects every horizontal line (scan line) twice at
% the most (though, its top-most and/or the bottom edge could be horizontal).
%
% See also: cv.line, cv.fillPoly
%


text:
classdef prtDataSetClassReshape < prtDataSetClass & prtDataInterfaceReshape
    % An helper class that interfaces from prtDataSetClass and prtDataInterfaceReshape
    %
    % There are no additional methods. This just combines prtDataSetClass and the
    % mixin prtDataInterfaceReshape
    %
    %  Example:
    %     ds = prtDataSetClassReshape(randn(20,100), [], 'observationSize',[10 10]);
    %
    %     size(ds.getObservations(1))
    %     size(ds.getObservations)
    %     size(ds.getObservationsAsMat(1))
    %     size(ds.getObservationsAsMat)
    %     size(ds.getObservationsAsMat([1 3]))
    %     size(ds.getObservationsAsCell([1 3]))
    %     size(ds.getObservationsAsCell())
    %     cell2mat(cellfun(@size,ds.getObservationsAsCell([1 3]),'uniformOutput',false))
    methods
        function self = prtDataSetClassReshape(varargin)
            self = self@prtDataSetClass(varargin{:});
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT