text:
function [apexPoints,signVals,exitCode]=trace2EarthMagApex(zCart,modSelParam,a,f,RelTol,AbsTol,maxSteps)
%TRACE2EARTHMAGAPEX Given a point in the in the International Terrestial
%           Reference System (ITRS), a type of Earth-Centered Earth-Fixed
%           (ECEF) system trace along the magnetic field line of the Earth
%           until the apex is found. The apex is the highest point of the
%           magnetic field line above the reference ellipsoid. Such tracing
%           is used for determining apex coordinates and quasi-dipole
%           coordinates, both of which play a role in ionospheric analyses.
%
%INPUTS: zCart One or more points given in Cartesian coordinates in the
%              ITRS with units of meters. zCart is a 3XN matrix with each
%              column having the format [x;y;z]. One would generally
%              assume that points are not underground.
%  modSelParam An optional parameter selecting or providing the magnetic
%              field model to use. If omitted or an empty matrix is
%              passed, the International Geomagnetic Reference Field
%              (IGRF) at the latest epoch of the model is used via the
%              function getIGRFCoeffs. Possible other values are
%              1) year The parameters of the IGRF for the specified epoch
%                 year are used. The year is in the Gregorian calendar and
%                 is specified as noted in the comments to the
%                 getIGRFCoeffs  function.
%              2) modSelParam is a structure where modSelParam.algorithm
%                 selects the algorithm to use Possible values for
%                 modSelParam.algorithm are 'IGRF', 'WMM' and 'preloaded'.
%                 If modSelParam.algorithm is 'IGRF or 'WMM', then the IGRF
%                 or the World Magnetic Model (WMM) is used at the
%                 fractional year specified by modSelParam.year (using the
%                 function getIGRFCoeffs or getWMMCoeffs). If
%                 modSelParam.algorithm is preloaded, then the parameters
%                 for the model are given by the modSelParam.C,
%                 modSelParam.S, modSelParam.a, and modSelParam.c, where
%                 the parameters have the same format as the fully
%                 normalized outputs of getIGRFCoeffs or getWMMCoeffs. See
%                 comments below for the use of custom coeffients.
%            a The semi-major axis of the reference ellipsoid. The apex is
%              defined in terms of a maximum magnetic field line height
%              above the reference ellipsoid. If this argument is omitted
%              or an empty matrix is passed, the value in
%              Constants.WGS84SemiMajorAxis is used.
%            f The flattening factor of the reference ellipsoid. If this
%              argument is omitted or an empty matrix is passed, the value
%              in Constants.WGS84Flattening is used.
%       RelTol The maximum relative error tolerance allowed for the
%              Runge-Kutta algorithm to trace the path towards the apex.
%              If omitted or an empty matrix is passed, the default value
%              of 1e-6 is used.
%       AbsTol The absolute error tolerance allowed, a positive scalar.
%              If omitted or an empty matrix is passed, the default value
%              of 1e-9 is used. This is used in the Runge-Kutta algorithm
%              to trace along the path towards the apex, and it is used
%              for the fminbnd function to find the apex location once it
%              has been bounded.
%     maxSteps The maximum allowable number of steps to perform the
%              adaptive Runge-Kutta integration along a magnetic field
%              line to find the apex. If omitted or an empty matrix is
%              passed, the default of 1024 is used.
%
%OUTPUS: apexPoints A 3XN matrix of the apex points corresponding to the
%               points in zCart. The apex points are obtained by tracing the
%               magnetic field line starting from zCart until reaching the
%               point farthest from the reference ellipsoid. Near the
%               geomagnetic poles, these points will be very far from the
%               Earth. Thus, if the point ends up being 1e37 m or more away
%               from the Earth, then integration is stopped and the
%               components of apexPoints are just set to Inf (rather than
%               returning an error). If  exitCode~=0, for any point, an
%               empty matrix is returned.
%      signVals An NX1 vector of values (+1, or -1) indicating whether the
%               direction to the apex point along the field line was
%               obtained by tracing in the direction of the magnetic field
%               or opposite to it. This can be used to determine which
%               magnetic hemisphere the points in zCart are in, because
%               +1 means that one starts closer to the magnetic North pole
%               (geographic South pole), and -1 means that one starts
%               closer to the magnetic South pole (geographic North pole).
%               If an error occurs, an empty matrix is returned.
%      exitCode A NX1 vector of codes indicating whether an error occurred.
%               If an error occurs, then the algorithm is halted at the
%               point causing the error.
%               0: Integration was successful.
%               1: Unable to get a small enough step size.
%               2: Apex not found within the maximum number of iterations.
%               3: Non-finite number encountered.
%               4: Problems performing a line search to the peak point.
%
%Basically, starting at the given point, one traces along a magnetic field
%line until reaching the highest point above the reference ellipsoid as
%defined in [2]. The tracing is done using an order 5(4) adaptive
%Runge-Kutta method until the direction of the field line is descending, in
%which case the peak is bounded between the last two points on the field
%line. An interpolation routine is then used with fminbnd to find the peak.
%This assumes that there is only one peak, which is safe to assume with the
%IMM and WMM. On the other hand, if using a custom magnetic field model,
%there might be more peaks and this function will not return the correct
%solution.
%
%Note that in [1], the definition of the apex is given in terms of a
%maximum height or a magnetic field line above the geoid, not above the
%reference ellipsoid (which is an approximation of the geoid). The
%definition with respect to the geoid seems to not be used in the
%literature. Moreover, the fast undulations of the geoid would make a
%magnetic field line have more than one peak, requiring that one find every
%peak to determine where the apex is.
%
%REFERENCES:
%[1] T. E. VanZandt, W. L. Clark, and J. M. Warnock, "Magnetic apex
%    coordinates: A magnetic coordinate system for the ionospheric f2
%    layer," Journal of Geophysical Research, vol. 77, no. 13, pp. 2406-
%    2411, 1 May 1972.
%[2] A. D. Richmond, "Ionospheric electrodynamics using magnetic apex
%    coordinates," Journal of Geomagnetism and Geoelectricity, vol. 47,
%    no. 2, pp. 191-212, 1995.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(modSelParam))
    %If no magnetic model is specified, use the IGRF at the reference
    %epoch.
    [C,S,aMagMod,cMagMod]=getIGRFCoeffs([],true);
elseif(isa(modSelParam,'double'))
    %If just a scalar is passed, assume that it is the epoch year for the
    %IGRF.
    [C,S,aMagMod,cMagMod]=getIGRFCoeffs(modSelParam,true);
else
    %Otherwise, a structrue should be passed indicating the model and the
    %year.
    switch(modSelParam.algorithm)
        case 'WMM'%World Magnetic Model
            [C,S,aMagMod,cMagMod]=getWMMCoeffs(modSelParam.year,true);
        case 'IGRF'%Intergnational Geomagnetic Reference Field
            [C,S,aMagMod,cMagMod]=getIGRFCoeffs(modSelParam.year,true);
        case 'preloaded'%Custom coefficients.
            C=modSelParam.C;
            S=modSelParam.S;
            aMagMod=modSelParam.a;
            cMagMod=modSelParam.c;
        otherwise
            error('Unknown model specified')
    end
end

if(nargin<3||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

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

if(nargin<5||isempty(RelTol))
    RelTol=1e-6;
end

if(nargin<6||isempty(AbsTol))
    AbsTol=1e-9;
end

if(nargin<7||isempty(maxSteps))
    maxSteps=1024;
end

%The parameters selecting the Runge-Kutta algorithm to use for the
%numerical integration.
order=5;
solutionChoice=0;

numPoints=size(zCart,2);
apexPoints=zeros(3,numPoints);
signVals=zeros(numPoints,1);
exitCode=zeros(numPoints,1);
for curPoint=1:numPoints
    zITRS=zCart(:,curPoint);

    %The initial stepsize is arbitrarily set to 10km.
    deltaS=10e3;
    %The maximum stepsize is arbitrarily set very high, so that it can be
    %quickly detected when the apex is extremely far away, because the
    %chosen point is very near the magnetic pole.
    deltaSMaxMag=1e36;

    xCur=zITRS;%The initial Cartesian location.
    sCur=0;%Distance traveled from the initial point is zero.

    %First, we have to decide which direction to go. We want to go in the
    %direction of increasing ellipsoidal height.
    dxdsCur=dxds(xCur,1,C,S,aMagMod,cMagMod);%The positive direction.

    uVert=getEllipsVert(zITRS,a,f);

    %If the angle between the direction and the local vertical is >90
    %degrees, then it is a decreasing direction. If it is greater, then it
    %is an increasing direction. This is the same as checking whether the
    %dot product is negative (the dot product is proportional to the cosine
    %of the angle).
    if(dot(uVert,dxdsCur)<0)
        signVal=-1;%Go the other way.
    else
        signVal=1;
    end
    signVals(curPoint)=signVal;

    %The path direction function.
    derivFun=@(x,s)dxds(x,signVal,C,S,aMagMod,cMagMod);

    %Now, we integrate along the path until the next point takes us to a
    %location where the direction along the magnetic field line is no longer
    %going in a direction of increasing ellipsoidal height. That means that we
    %have passed the apex. Thus, the apex height can be interpolated beteen
    %those final two points.
    apexPoint=[];
    for curStep=1:maxSteps
        xPrev=xCur;
        sPrev=sCur;

        %Arbitrary minimum step size.
        deltaSMinMag=2^4*eps(sCur);

        [deltaS,xCur,sCur,k,dxdsCur,exitCode(curPoint)]=performOneAdaptiveRKStep(xCur,sCur,derivFun,deltaS,deltaSMinMag,deltaSMaxMag,dxdsCur,order,solutionChoice,AbsTol,RelTol);

        if(exitCode(curPoint)~=0)
            signVals=[];
            apexPoints=[];
            return;
        end

        %Determine whether the new point is still a point of increasing
        %ellipsoidal height.
        uVert=getEllipsVert(xCur,a,f);

        %If the peak has been reached.
        if(dot(uVert,dxdsCur)<0)
            %Interpolate to find the peak.

            %Get a Hermite interpolating polynomial over the step. Note that the
            %interpolating polynomial takes a parameter from 0->1 indicating the
            %fraction of the distance traveled between sPrev and sCur.
            [interpPolyA,interpPolyC]=RKInterpPolys(xPrev,sPrev,xCur,sCur,derivFun,order,solutionChoice,k);

            %We will now do a simple line search to find where the direction of
            %the magnetic field vector is 90 degrees offset.
            costFun=@(sFracEst)levelCostFun(sFracEst,interpPolyA,interpPolyC,a,f);
            try
                options=optimset('TolX',AbsTol);
                minFrac=fminbnd(costFun,0,1,options);
            catch%If some error in fminbnd occurred.
                apexPoints=[];
                exitCode(curPoint)=4;
                return
            end

            apexPoint=polyValNewton(minFrac,interpPolyA,interpPolyC);

            break;
        elseif(Cart2Ellipse(xCur,[],a,f)>=1e37)
            %If the ellipsoidal height exceeds 1e37m, then just declare it
            %to be at infinity. This helps near the magnetic poles.
            apexPoint=[Inf;Inf;Inf];
        end
    end

    if(isempty(apexPoint))
       %If the maximum number of iterations was achieved without hitting the
       %apex point, then return with an error.
        apexPoints=[];
        signVals=[];
        exitCode(curPoint)=2;
        return;
    end

    %We have the apexPoint. We can now compute the apex coordinates.
    apexPoints(:,curPoint)=apexPoint;
    exitCode(curPoint)=0;
end

end

function uVert=getEllipsVert(xCart,a,f)
%Get a vector pointing in the direction of the ellipsoidal vertical. This
%requires the position to be converted into ellipsoidal coordinates.

    ellipsCur=Cart2Ellipse(xCart,[],a,f);
    justVertical=true;
    uVert=getENUAxes(ellipsCur,justVertical,a,f);
end

function costVal=levelCostFun(sFracEst,interpPolyA,interpPolyC,a,f)
%The cost function to determine when the angle between the local vertical
%and the magnetic field is 90 degrees. We are trying to drive the dot
%product to zero.

    xValCart=polyValNewton(sFracEst,interpPolyA,interpPolyC);
    dxdsFracCart=polyDerValNewton(sFracEst,interpPolyA,interpPolyC);

    uVert=getEllipsVert(xValCart,a,f);

    costVal=abs(dot(uVert,dxdsFracCart));
end

function derivVal=dxds(xCart,signVal,C,S,aMagMod,cMagMod)
%The function providing the derivative of position with respect to
%arclength for tracing out the field line.

    pointSpher=Cart2Sphere(xCart);
    [~,gradV]=spherHarmonicEval(C,S,pointSpher,aMagMod,cMagMod);
    B=-gradV;%The magnetic flux

    %Integration is performed in space, so the normalized B vector is the
    %position derivative for the differential equation. The problem of running
    %into a point where B=0 should not arise.
    derivVal=signVal*B/norm(B);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Shear Waves And Critical Angle Reflection Example
%
% This example illustrates Snell's law for elastic media using a weakly
% focused ultrasound transducer incident on a soft-tissue / bone interface.
% It builds on the Explosive Source In A Layered Medium and Snell's Law And
% Critical Angle Reflection examples.
%
% author: Bradley Treeby
% date: 14th Nov 2013
% last update: 15th February 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

clear all;

% =========================================================================
% SIMULATION PARAMETERS
% =========================================================================

% change scale to 2 to reproduce the higher resolution figures used in the
% help file
scale = 1;

% create the computational grid
PML_size = 10;               % size of the PML in grid points
Nx = 128*scale - 2*PML_size; % number of grid points in the x direction
Ny = 192*scale - 2*PML_size; % number of grid points in the y direction
dx = 0.1e-3/scale;           % grid point spacing in the x direction [m]
dy = 0.1e-3/scale;           % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the medium properties
cp1 = 1540;         % compressional wave speed [m/s]
cs1 = 0;            % shear wave speed [m/s]
rho1 = 1000;        % density [kg/m^3]
alpha0_p1 = 0.1;    % compressional wave absorption [dB/(MHz^2 cm)]
alpha0_s1 = 0.1;    % shear wave absorption [dB/(MHz^2 cm)]

cp2 = 3000;         % compressional wave speed [m/s]
cs2 = 1400;         % shear wave speed [m/s]
rho2 = 1850;        % density [kg/m^3]
alpha0_p2 = 1;      % compressional wave absorption [dB/(MHz^2 cm)]
alpha0_s2 = 1;      % shear wave absorption [dB/(MHz^2 cm)]

% create the time array
cfl   = 0.1;
t_end = 12e-6;
kgrid.t_array= makeTime(kgrid, cp1, cfl, t_end);

% define position of heterogeneous slab
slab = zeros(Nx, Ny);
slab(Nx/2:Nx, :) = 1;

% define the source properties
source_freq = 2e6;              % [Hz]
source_strength = 1e6;          % [Pa]
source_cycles = 3;              % number of tone burst cycles
source_focus_dist = 5*scale;    % position of focus inside slab [grid points]
source_slab_dist = 5*scale;     % distance between source and slab [grid points]
source_mask = makeCircle(Nx, Ny, Nx/2 + source_focus_dist, Ny/2, 50*scale, pi/3);
source_mask(Nx/2 - source_slab_dist:end, :) = 0;

% define the sensor to record the maximum particle velocity everywhere
sensor.record = {'u_max_all'};

% set the input arguments
input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ...
    'PMLInside', false, 'PlotScale', [-1, 1]*source_strength, ...
    'DisplayMask', 'off', 'DataCast', 'single'};

% =========================================================================
% FLUID SIMULATION
% =========================================================================

% assign the medium properties
medium.sound_speed            = cp1*ones(Nx, Ny);
medium.sound_speed(slab == 1) = cp2;
medium.density                = rho1*ones(Nx, Ny);
medium.density(slab == 1)     = rho2;
medium.alpha_coeff            = alpha0_p1*ones(Nx, Ny);
medium.alpha_coeff(slab == 1) = alpha0_p2;
medium.alpha_power            = 2;

% assign the source
source.p_mask = source_mask;
source.p = source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);

% run the fluid simulation
sensor_data_fluid = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% ELASTIC SIMULATION
% =========================================================================

% define the medium properties
clear medium
medium.sound_speed_compression            = cp1*ones(Nx, Ny);
medium.sound_speed_compression(slab == 1) = cp2;
medium.sound_speed_shear                  = cs1*ones(Nx, Ny);
medium.sound_speed_shear(slab == 1)       = cs2;
medium.density                            = rho1*ones(Nx, Ny);
medium.density(slab == 1)                 = rho2;
medium.alpha_coeff_compression            = alpha0_p1*ones(Nx, Ny);
medium.alpha_coeff_compression(slab == 1) = alpha0_p2;
medium.alpha_coeff_shear                  = alpha0_s1*ones(Nx, Ny);
medium.alpha_coeff_shear(slab == 1)       = alpha0_s2;

% assign the source
clear source
source.s_mask = source_mask;
source.sxx = -source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);
source.syy = source.sxx;

% run the elastic simulation
sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% VISUALISATION
% =========================================================================

% define plot vector
x_vec = kgrid.x_vec(1 + PML_size:end - PML_size)*1e3;
y_vec = kgrid.y_vec(1 + PML_size:end - PML_size)*1e3;

% calculate square of velocity magnitude
u_e = sensor_data_elastic.ux_max_all.^2 + sensor_data_elastic.uy_max_all.^2;
u_f = sensor_data_fluid.ux_max_all.^2 + sensor_data_fluid.uy_max_all.^2;

% plot layout
figure;
imagesc(y_vec, x_vec, double(source_mask | slab));
xlabel('y [mm]');
ylabel('x [mm]');
axis image;
colormap(flipud(gray));

% plot beam patterns
figure;
subplot(2, 1, 1);
imagesc(y_vec, x_vec, 20*log10(u_f./max(u_f(:))));
xlabel('y [mm]');
ylabel('x [mm]');
axis image;
colorbar;
caxis([-50, 0]);
title('Fluid Model');

subplot(2, 1, 2);
imagesc(y_vec, x_vec, 20*log10(u_e./max(u_e(:))));
xlabel('y [mm]');
ylabel('x [mm]');
axis image;
colorbar;
caxis([-50, 0]);
title('Elastic Model');
colormap(jet(256));

repo: Aero-matlab


text:
function [count, value] = mrAnatHistogramSmooth(img, nbins, fwhm)
%
% [count, value] = mrAnatHistogramSmooth(img, nbins, fwhm)
%
% Returns a smoothed image histogram.
%
% * nbins should be even (defaults to 256)
% * fwhm is the full-width,half-max of the gaussian smoothing kernel,
% expressed as a proportion of nbins. (defaults to 0.1)
%
% HISTORY:
% 2005.02.05 RFD: wrote it.

if(~exist('nbins','var') | isempty(nbins))
    nbins = 256;
end
if(~exist('fwhm','var') | isempty(fwhm))
    fwhm = 0.1;
end
[count,value] = hist(img(:), nbins);
a = fwhm/2.3548;
x = [-(nbins-1)/2:(nbins-1)/2]/(nbins/2);
g = a*sqrt(2*pi) * exp(-0.5 * (x./a).^2) ./ (sqrt(2*pi) .* a);
count = conv(count, g);
count = count(nbins/2:end-nbins/2);
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function cPoly=composePolys(a,b)
%%COMPOSEPOLYS Compose two polynomials given as power series. This means,
%              for polynomials A(x) and B(x) in a scalar x, evaluate the
%              composition C(x)=A(B(x)). C(x) itself is a polynomial. Thus,
%              the output of the function is the set of coefficients of C
%              (subject to finite precision limitations in the recursion).
%
%INPUTS: a Coefficients for a power series of the form
%          y(z)=a(end)+a(end-1)*z+a(end-2)*z^2+...
%          into which the power series in b is to be substituted. The
%          ordering of the coefficients is the same as in Matlab's polyval
%          function.
%        b Coefficients for the power series that is substitued into a.
%
%OUTPUTS: c Coefficients for the power series that represents the
%           substitution of the polynomial b into a.
%
%The substitution is performed by explicitly evaluating power of b (via
%repeated convolutions), multiplying the powers by the appropriate value in
%a and summing the results.
%
%February 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

orderA=length(a)-1;
orderB=length(b)-1;

%The order of the output.
prodOrder=orderA*orderB;

%Zero-pad b so that it is the same length as the length of the output. This
%simplifies assignments in the loop.
b=[zeros(prodOrder-orderB,1);b(:)];

%bProd will hold powers of b. Right now, it is just initialized to 1 (b^0).
bProd=[zeros(prodOrder,1);1];

%Initialize the output with the zeroth-order term.
cPoly=zeros(prodOrder+1,1);
cPoly=cPoly+a(end)*bProd;
for curPow=1:orderA
    %Get the polynomial for the next power of b. The convolution operation
    %is equivalent to polynomial multiplication. The extra leading terms
    %(which are all zeros) are just thrown out.
    bProd=conv(b,bProd);
    bProd=bProd((prodOrder+1):end);

    %Add in the next power term.
    cPoly=cPoly+a(end-curPow)*bProd;
end
end

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


text:
%% DEMO_MixedTetHexMeshing
% Below is a demonstration of how to create a mixed mesh consisting of
% (linear) hexahedral and tetrahedral elements. The hexahedral mesh is
% regular while the tetrahedral mesh is derived using TetGen.

%%

clear; close all; clc;

%%
% Plot settings for the examples below
fontSize=20;
faceAlpha1=1;
faceAlpha2=0.3;
plotColors=gjet(4);

%%

searchRadius=6;

%% CONVERTING A TRIANGULATED SURFACE TO AN IMAGE WITH DESIRED SIZE, VOXEL SIZE AND ORIGIN
% Defining an example triangulated surface model

[Fs,Vs]=stanford_bunny;

%%
% Setting control parameters

% Defining the full set of possible control parameters
voxelSize=6; % The output image voxel size.
imOrigin=min(Vs,[],1)-voxelSize;
imMax=max(Vs,[],1)+voxelSize;
imSiz=round((imMax-imOrigin)/voxelSize);
imSiz=imSiz([2 1 3]); %Image size (x, y corresponds to j,i in image coordinates, hence the permutation)

% Using |triSurf2Im| function to convert patch data to image data
[M,~]=triSurf2Im(Fs,Vs,voxelSize,imOrigin,imSiz);

%%
% Plotting the results

hf1=cFigure;
subplot(1,2,1);
title('Closed triangulated surface','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;
gpatch(Fs,Vs,'g','none');
axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
camlight('headlight'); lighting phong;
set(gca,'fontSize',fontSize);

subplot(1,2,2);
title('Boundary, intertior and exterior image','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fs,Vs,'g','none',faceAlpha2);

L_plot=false(size(M));
L_plot(:,:,round(size(M,3)/2))=1;
[Fm,Vm,Cm]=ind2patch(L_plot,double(M),'sk');
[Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3));
Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
gpatch(Fm,Vm,Cm,'k');

L_plot=false(size(M));L_plot(round(size(M,1)/2),:,:)=1;
[Fm,Vm,Cm]=ind2patch(L_plot,M,'si');
[Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3));
Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
gpatch(Fm,Vm,Cm,'k');

L_plot=false(size(M));L_plot(:,round(size(M,2)/2),:)=1;
[Fm,Vm,Cm]=ind2patch(L_plot,M,'sj');
[Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3));
Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
gpatch(Fm,Vm,Cm,'k');

colormap(gray(3)); caxis([0 2]);
hc=colorbar;
set(hc,'YTick',[1/3 1 5/3]);
set(hc,'YTickLabel',{'Exterior','Boundary','Intertior'});
set(hc,'fontSize',fontSize);
axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
set(gca,'fontSize',fontSize);
drawnow;

%% GET HEXAHEDRAL ELEMENT SET

L_model=(M==2); %Interior&Boundary choosen here

%Defining erosion/dilation kernel
k=3;
p=k-round(k./2);
hb=zeros(3,3);
hb(2,2,2)=1;
hb(2,2,1)=1;
hb(2,2,3)=1;
hb(1,2,2)=1;
hb(3,2,2)=1;
hb(2,3,2)=1;
hb(2,1,2)=1;

L_model_rep=zeros(size(L_model)+(2.*p));
L_model_rep(p+(1:size(L_model,1)),p+(1:size(L_model,2)),p+(1:size(L_model,3)))=L_model;
L_model_blur = convn(double(L_model_rep),hb,'valid');
L_model=L_model_blur>=(sum(hb(:)));

[E_hex,V_hex,C_hex]=ind2patch(L_model,M,'hu');

% Convert Coordinates
[V_hex(:,1),V_hex(:,2),V_hex(:,3)]=im2cart(V_hex(:,2),V_hex(:,1),V_hex(:,3),voxelSize*ones(1,3));
V_hex=V_hex+imOrigin(ones(size(V_hex,1),1),:);

% Use element2patch to get patch data to plot the model
[F_hex_cut1,C_hex_F]=element2patch(E_hex,C_hex);

%Pass through unique_patch to reduce "weight" of plot
[Fp,Vp,~,~,~,F_count]=unique_patch(F_hex_cut1,V_hex,[],5);
logicUni=F_count==1; %Logic for boundary faces

Fq=Fp(logicUni,:);
Vq=Vp;
[Fq,Vq,~]=patchCleanUnused(Fq,Vq);

[Ft,Vt]=quad2tri(Fq,Vq,'b');

%%
% Plotting the results

hf1=cFigure;
title('Visualizing internal voxels=hexahedral elements','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fs,Vs,0.5*ones(1,3),'none',faceAlpha2);
gpatch(Fq,Vq,plotColors(2,:),'k');

camlight('headlight'); lighting flat;
axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
drawnow;

%%

%Joining surface sets
F=[Fs;Ft+size(Vs,1)];
V=[Vs;Vt];
C_tet=[ones(size(Fs,1),1);2*ones(size(Ft,1),1)]; %Surface marker colors

%% Get hole point
[V_hole]=getInnerPoint(Ft,Vt,searchRadius,voxelSize/2,0);
plotV(V_hole,'r.','MarkerSize',25);

%% Get region point

L_in=(M==1);

[indInternal]=getInnerVoxel(double(L_in),searchRadius,0);

[I_in,J_in,K_in]=ind2sub(size(L_in),indInternal); %Convert to subscript coordinates
[X_in,Y_in,Z_in]=im2cart(I_in,J_in,K_in,voxelSize*ones(1,3));
V_in=[X_in Y_in Z_in];
V_in=V_in+imOrigin(ones(size(V_in,1),1),:);

plotV(V_in,'k.','MarkerSize',25);

%%
% DEFINE FACE BOUNDARY MARKERS
faceBoundaryMarker=C_tet;

%%
% Define region points
V_regions=[V_in];

%%
% Define hole points
V_holes=[V_hole];

%%
% Regional mesh parameters
[edgeLengths]=patchEdgeLengths(F,V);
edgeLengthsMean=mean(edgeLengths);
meanProposedVolume=edgeLengthsMean^3./(6*sqrt(2)); %For regular tetrahedron
regionA=meanProposedVolume;

%%
% CREATING THE SMESH STRUCTURE, meshing without the surface constraints
% imposed by the -Y this time.

stringOpt='-pq1.2AaY';

modelName='tetGenModel';

meshStruct.stringOpt=stringOpt;
meshStruct.Faces=F;
meshStruct.Nodes=V;
meshStruct.holePoints=V_holes;
meshStruct.faceBoundaryMarker=faceBoundaryMarker; %Face boundary markers
meshStruct.regionPoints=V_regions; %region points
meshStruct.regionA=regionA;
meshStruct.minRegionMarker=2; %Minimum region marker
meshStruct.modelName=modelName;

%%
% Mesh model using tetrahedral elements using tetGen (see:
% <http://wias-berlin.de/software/tetgen/>)

[meshOutput]=runTetGen(meshStruct); %Run tetGen

%%
% Access model element and patch data
F_tet_cut1=meshOutput.faces;
V_tet=meshOutput.nodes;
C_tet=meshOutput.faceMaterialID;
E_tet=meshOutput.elements;

indBoundary=meshOutput.facesBoundary(meshOutput.boundaryMarker==1,:);
indBoundary=unique(indBoundary(:));

%% MERGING NODE SETS

V=[V_tet;V_hex];
E_hex=E_hex+size(V_tet,1);

[~,V,ind1,ind2]=mergeVertices(F_tet_cut1,V);

E_tet=ind2(E_tet);
E_hex=ind2(E_hex);
indBoundary=ind2(indBoundary);

E={E_tet, E_hex};

%%
% Visualizing mesh

cFigure;
title('Mixed TET/HEX mesh','FontSize',fontSize);

%Selecting half of the model to see interior
X=V(:,2); XE=mean(X(E{1}),2);
L=XE>mean(X);
[F_tet_cut1,~]=element2patch(E{1}(L,:),C_tet(L));

%Selecting half of the model to see interior
X=V(:,2); XE=mean(X(E{2}),2);
L=XE>mean(X);
[F_hex_cut1,~]=element2patch(E{2}(L,:),C_hex(L));

gpatch(F_tet_cut1,V,plotColors(1,:),'k');
gpatch(F_hex_cut1,V,plotColors(2,:),'k');
gpatch(Fs,Vs,0.5*ones(1,3),'none',faceAlpha2);

axisGeom(gca,fontSize);
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

%%
% Smoothing meshes

cPar.Method='LAP';
cPar.n=25;
cPar.RigidConstraints=indBoundary;

[F1,~]=element2patch(E{1},[]);
[F1,~,~]=uniqueIntegerRow(F1);
[F2,~]=element2patch(E{2},[]);
[F2,~,~]=uniqueIntegerRow(F2);
[~,IND_V1,~]=tesIND(F1,V,0);
[~,IND_V2,~]=tesIND(F2,V,0);
IND_V=[IND_V1 IND_V2];

[VS]=tesSmooth([],V,IND_V,cPar);

%%
% Visualizing mesh

hf=cFigure;
title('Mixed TET/HEX mesh','FontSize',fontSize);

%Selecting half of the model to see interior
X=V(:,2); XE=mean(X(E{1}),2);
L=XE>mean(X);
[F_tet_cut1,~]=element2patch(E{1}(L,:),C_tet(L));

%Selecting half of the model to see interior
X=V(:,2); XE=mean(X(E{2}),2);
L=XE>mean(X);
[F_hex_cut1,~]=element2patch(E{2}(L,:),C_hex(L));

hp1=gpatch(F_tet_cut1,VS,plotColors(1,:),'k');
hp2=gpatch(F_hex_cut1,VS,plotColors(2,:),'k');
gpatch(Fs,Vs,0.5*ones(1,3),'none',faceAlpha2);

axisGeom(gca,fontSize);
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

%%
% Set up animation
nSteps=25; %Number of animation steps
X=V(:,2);
XE1=mean(X(E{1}),2);
XE2=mean(X(E{2}),2);

animStruct.Time=linspace(0,1,nSteps); %Time vector
cutLevel=linspace(min(X(:)),max(X(:)),nSteps); %Property to set

for q=1:1:nSteps %Step through time
    cutLevelNow=cutLevel(q); %The current cut level

    L1=XE1>cutLevelNow;
    [F_tet_cut1,~]=element2patch(E{1}(L1,:));

    L2=XE2>cutLevelNow;
    [F_hex_cut1,~]=element2patch(E{2}(L2,:));

    %Set entries in animation structure
    animStruct.Handles{q}=[hp1 hp2]; %Handles of objects to animate
    animStruct.Props{q}={'Faces','Faces'}; %Properties of objects to animate
    animStruct.Set{q}={F_tet_cut1,F_hex_cut1}; %Property values for to set in order to animate
end

%Add animation layer
anim8(hf,animStruct);

%%
%
% <<gibbVerySmall.gif>>
%
% GIBBON
%
% Kevin M. Moerman (kevinmoerman@hotmail.com)

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function score = nmi_onmf(true_labels, cluster_labels)
%NMI Compute normalized mutual information (NMI) using the true and cluster
%   labels and return the value in 'score'.
%
%   Input    : true_labels    : N-by-1 vector containing true labels
%              cluster_labels : N-by-1 vector containing cluster labels
%
%   Output   : score          : NMI value
%
%   Reference: Shi Zhong, 2003.
%              http://www.cse.fau.edu/~zhong/software/textclust.zip

% Compute the confusion matrix 'cmat', where
%   col index is for true label (CAT),
%   row index is for cluster label (CLS).
n = length(true_labels);
cat = spconvert([(1:n)' true_labels ones(n,1)]);
cls = spconvert([(1:n)' cluster_labels ones(n,1)]);
cls = cls';
cmat = full(cls * cat);

n_i = sum(cmat, 1); % Total number of data for each true label (CAT), n_i
n_j = sum(cmat, 2); % Total number of data for each cluster label (CLS), n_j

% Calculate n*n_ij / n_i*n_j
[row, col] = size(cmat);
product = repmat(n_i, [row, 1]) .* repmat(n_j, [1, col]);
index = find(product > 0);
n = sum(cmat(:));
product(index) = (n*cmat(index)) ./ product(index);
% Sum up n_ij*log()
index = find(product > 0);
product(index) = log(product(index));
product = cmat .* product;
score = sum(product(:));
% Divide by sqrt( sum(n_i*log(n_i/n)) * sum(n_j*log(n_j/n)) )
index = find(n_i > 0);
n_i(index) = n_i(index) .* log(n_i(index)/n);
index = find(n_j > 0);
n_j(index) = n_j(index) .* log(n_j(index)/n);
denominator = sqrt(sum(n_i) * sum(n_j));

% Check if the denominator is zero
if denominator == 0
  score = 0;
else
  score = score / denominator;
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function area = polygonArea3d(poly, varargin)
%POLYGONAREA3D Area of a 3D polygon
%
%   AREA = polygonArea3d(POLY)
%   POLY is given as a N-by-3 array of vertex coordinates. The resulting
%   area is positive.
%   Works also for polygons given as a cell array of polygons.
%
%   Example
%     % area of a simple 3D square
%     poly = [10 30 20;20 30 20;20 40 20;10 40 20];
%     polygonArea3d(poly)
%     ans =
%        100
%
%     % Area of a 3D mesh
%     [v f] = createCubeOctahedron;
%     polygons = meshFacePolygons(v, f);
%     areas = polygonArea3d(polygons);
%     sum(areas)
%     ans =
%         18.9282
%
%   See also
%     polygons3d, triangleArea3d, polygonArea, polygonCentroid3d

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

%   HISTORY
%   2013-08-20 add support for multiple polygons

% Check multiple polygons
if iscell(poly) || sum(sum(isnan(poly))) > 0
    % split the polygons into a cell array
    polygons = splitPolygons3d(poly);
    nPolys = length(polygons);

    % compute area of each polygon
    area = zeros(nPolys, 1);
    for  i = 1:nPolys
        area(i) = polygonArea3d(polygons{i});
    end

    return;
end

% put the first vertex at origin (reducing computation errors for polygons
% far from origin)
v0 = poly(1, :);
poly = bsxfun(@minus, poly, v0);

% indices of next vertices
N = size(poly, 1);
iNext = [2:N 1];

% compute cross-product of each elementary triangle formed by origin and
% two consecutive vertices
cp = cross(poly, poly(iNext,:), 2);

% choose one of the triangles as reference for the normal direction
vn = vectorNorm(cp);
[tmp, ind] = max(vn); %#ok<ASGLU>
cpRef = cp(ind,:);

% compute the sign of the area of each triangle
% (need to compute the sign explicitely, as the norm of the cross product
% does not keep orientation within supporting plane)
sign_i = sign(dot(cp, repmat(cpRef, N, 1), 2));

% compute area of each triangle, using sign correction
area_i = vectorNorm3d(cp) .* sign_i;

% sum up individual triangles area
area = sum(area_i) / 2;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function cs_install (do_pause)
%CS_INSTALL: compile and install CSparse for use in MATLAB.
%   Your current working directory must be CSparse/MATLAB in order to use this
%   function.
%
%   The directories
%
%       CSparse/MATLAB/CSparse
%       CSparse/MATLAB/Demo
%       CSparse/MATLAB/UFget
%
%   are added to your MATLAB path (see the "pathtool" command to add these to
%   your path permanently, for future MATLAB sessions).
%
%   Next, the MATLAB CSparse demo program, CSparse/MATLAB/cs_demo is executed.
%   To run the demo with pauses so you can see the results, use cs_install(1).
%   To run the full MATLAB test programs for CSparse, run testall in the
%   Test directory.
%
%   Example:
%       cs_install          % install and run demo with no pauses
%       cs_install(1)       % install and run demo with pauses
%
%   See also: cs_demo
%
%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

fprintf ('Compiling and installing CSparse\n') ;
if (nargin < 1)
    do_pause = 0 ;
end

if (~isempty (strfind (computer, '64')))
    error ('64-bit version not supported; use CXSparse instead') ;
end

if (do_pause)
    input ('Hit enter to continue: ') ;
end
addpath ([pwd filesep 'CSparse']) ;
addpath ([pwd filesep 'Demo']) ;

v = getversion ;
if (v >= 7.0)
    addpath ([pwd filesep 'UFget']) ;
else
    fprintf ('UFget not installed (MATLAB 7.0 or later required)\n') ;
end

cd ('CSparse') ;
cs_make (1) ;
cd ('../Demo') ;
cs_demo (do_pause)

%-------------------------------------------------------------------------------
function v = getversion
% determine the MATLAB version, and return it as a double.
v = sscanf (version, '%d.%d.%d') ;
v = 10.^(0:-1:-(length(v)-1)) * v ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Y, spectrum] = sllemap(G, d, sch)
%SLLEMAP Solves Laplacian Eigenmap Embedding
%
% $ Syntax $
%   - Y = sllemap(G, d)
%   - Y = sllemap(G, d, sch)
%   - [Y, spectrum] = sllemap(...)
%
% $ Arguments $
%   - G:        The affinity graph (in any acceptable form): n x n
%   - d:        The embedding dimension
%   - sch:      The scheme to use
%   - Y:        The solved embedded coordinates (d x n)
%
% $ Description $
%   - Y = sllemap(G, d) uses the default scheme to solve the Laplacian
%     Eigenmap embedding in a d-dimensional space.
%
%   - Y = sllemap(G, d, sch) uses the specified scheme to solve the
%     Laplacian Eigenmap embedding in a d-dimensional space.
%
%     Three schemes are implemented to resolve the problem, they are
%     under three different formulations:
%       (1) 'minLI':
%           objective: minimize sum_ij w_ij ||y_i - y_j||^2
%                      s.t. forall i, ||y_i||^2 = 1
%           in matrix form, it is expressed as:
%               minimize y^T * L * y,  s.t. y^T * y = 1
%       (2) 'minLD':
%           objective: minimize sum_ij w_ij ||y_i - y_j||^2
%                      s.t. forall i, d_ii ||y_i||^2 = 1
%           in matrix form, it is expressed as:
%               minimize y^T * L * y, s.t. y^T * D * y = 1
%           This is the original formulation in many papers in the fields
%           of spectral learning, clustering and manifold learning.
%       (3) 'maxWD': (default)
%           objective: maximize sum_ij w_ij <y_i, y_j>
%                      s.t. forall i, d_ii ||y_i||^2 = 1
%           in matrix form, it is expressed as:
%               maximize y^T * W * y, s.t. y^T * D * y = 1
%           In theory, this objective is equivalent to 'minLD'. However,
%           due to that its implementation is based on finding the
%           eigenvectors corresponding to the largest eigenvalues instead
%           of the smallest ones, thus it is much more efficient and
%           numerically stable. Hence, it is selected as the default
%           scheme.
%
%   - [Y, spectrum] = sllemap(...) additionally outputs the spectrum
%     of the embedding. The definition of the spectrum varies for different
%     schemes:
%     'minLI': the spectrum is the eigenvalues of L, in ascending order
%     'minLD': the spectrum is the eigenvalues of D^(-1/2) * L * D^(-1/2),
%              in ascending order
%     'maxWD': the spectrum is the eigenvalues of D^(-1/2) * W * D^(-1/2),
%              in descending order.
%
% $ Remarks $
%   - If the input graph does not have edge values or it is logical,
%     it just assume 1 between neighboring samples and 0 for other pairs.
%
%   - The embedding dimension d should be strictly less than the number
%     of samples n.
%
% $ History $
%   - Created by Dahua Lin, on Sep 12nd, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('sllemap', 2);
end

gi = slgraphinfo(G, {'square'});
n = gi.n;
if strcmp(gi, 'adjmat')
    if isnumeric(G)
        W = G;
    else
        W = double(G);
    end
else
    W = sladjmat(G);
end

if d >= n
    error('sltoolbox:invalidarg', ...
        'The embeded dimension d should be strictly less than n');
end

if nargin < 3 || isempty(sch)
    sch = 'maxWD';
end

%% main delegation

% L = Di + Dj - Wij - Wji
% we let
%   W = Wij + Wji
%   D = Di + Dj = diag(diag(W))
%   L = D - W
W = W + W';

switch sch
    case 'maxWD'
        [Y, spectrum] = solve_maxWD(W, d);
    case 'minLD'
        [Y, spectrum] = solve_minLD(W, d);
    case 'minLI'
        [Y, spectrum] = solve_minLI(W, d);
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid scheme for solving eigenmap: %s', sch);
end

%% Core functions

function [Y, spectrum] = solve_maxWD(W, d)

vD = calcDv(W);
W = calcNormalizeMat(W, vD);

[spectrum, Y] = slsymeig(W, d+1, 'descend');
spectrum = spectrum(2:d+1);
Y = Y(:, 2:d+1)';
Y = denormY(Y, vD);

function [Y, spectrum] = solve_minLD(W, d)

vD = calcDv(W);
L = calcL(vD, W);
L = calcNormalizeMat(L, vD);

[spectrum, Y] = slsymeig(L, d+1, 'ascend');
spectrum = spectrum(2:d+1);
Y = Y(:, 2:d+1)';
Y = denormY(Y, vD);

function [Y, spectrum] = solve_minLI(W, d)

vD = calcDv(W);
L = calcL(vD, W);

[spectrum, Y] = slsymeig(L, d+1, 'ascend');
spectrum = spectrum(2:d+1);
Y = Y(:, 2:d+1)';

%% Computation routines

function vD = calcDv(W)

vD = sum(W, 1)';

function L = calcL(vD, W)

n = length(vD);
if issparse(W)
    D = sparse((1:n)', (1:n)', vD, n, n, n);
    L = D - W;
else
    L = -W;
    dinds = (1:n)'*(n+1) - n;
    L(dinds) = L(dinds) + vD;
end

function Mn = calcNormalizeMat(M, vD)

vD(vD < eps) = eps;
cv = 1 ./ sqrt(vD);

if issparse(M)
    n = size(M,1);
    Mn = M;
    for i = 1 : n
        Mn(:,i) = Mn(:,i) * cv(i);
    end
    for i = 1 : n
        Mn(i,:) = Mn(i,:) * cv(i);
    end
else
    rv = cv';
    Mn = slmulrowcols(M, rv, cv);
end

function Y = denormY(Y, vD)

vD = vD';
vD(vD < eps) = eps;
Y = slmulvec(Y, 1 ./ sqrt(vD), 2);

author: lmthang


text:
% readMNIST by Siddharth Hegde
%
% Description:
% Read digits and labels from raw MNIST data files
% File format as specified on http://yann.lecun.com/exdb/mnist/
% Note: The 4 pixel padding around the digits will be remove
%       Pixel values will be normalised to the [0...1] range
%
% Usage:
% [imgs labels] = readMNIST(imgFile, labelFile, readDigits, offset)
%
% Parameters:
% imgFile = name of the image file
% labelFile = name of the label file
% readDigits = number of digits to be read
% offset = skips the first offset number of digits before reading starts
%
% Returns:
% imgs = 20 x 20 x readDigits sized matrix of digits
% labels = readDigits x 1 matrix containing labels for each digit
%
function [imgs labels] = readMNIST(imgFile, labelFile, readDigits, offset)

    % Read digits
    fid = fopen(imgFile, 'r', 'b');
    header = fread(fid, 1, 'int32');
    if header ~= 2051
        error('Invalid image file header');
    end
    count = fread(fid, 1, 'int32');
    if count < readDigits+offset
        error('Trying to read too many digits');
    end

    h = fread(fid, 1, 'int32');
    w = fread(fid, 1, 'int32');

    if offset > 0
        fseek(fid, w*h*offset, 'cof');
    end

    imgs = zeros([h w readDigits]);

    for i=1:readDigits
        for y=1:h
            imgs(y,:,i) = fread(fid, w, 'uint8');
        end
    end

    fclose(fid);

    % Read digit labels
    fid = fopen(labelFile, 'r', 'b');
    header = fread(fid, 1, 'int32');
    if header ~= 2049
        error('Invalid label file header');
    end
    count = fread(fid, 1, 'int32');
    if count < readDigits+offset
        error('Trying to read too many digits');
    end

    if offset > 0
        fseek(fid, offset, 'cof');
    end

    labels = fread(fid, readDigits, 'uint8');
    fclose(fid);

    % Calc avg digit and count
    imgs = trimDigits(imgs, 4);
    imgs = normalizePixValue(imgs);
    %[avg num stddev] = getDigitStats(imgs, labels);

end

function digits = trimDigits(digitsIn, border)
    dSize = size(digitsIn);
    digits = zeros([dSize(1)-(border*2) dSize(2)-(border*2) dSize(3)]);
    for i=1:dSize(3)
        digits(:,:,i) = digitsIn(border+1:dSize(1)-border, border+1:dSize(2)-border, i);
    end
end

function digits = normalizePixValue(digits)
    digits = double(digits);
    for i=1:size(digits, 3)
        digits(:,:,i) = digits(:,:,i)./255.0;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function y = log_normcdf( x, approx ) %#ok

%LOG_NORMCDF   Logarithm of the cumulative normal distribution.
%   Y = LOG_NORMCDF(X) is the logarithm of the CDF of the normal
%   distribution at the point X.
%
%                                1    / x
%       LOG_NORMCDF(X) = LOG( ------- |   exp(-t^2/2) dt )
%                             sqrt(2) / -Inf
%
%   For numeric X, LOG_NORMCDF(X) is computed using the equivalent
%   expression LOG(0.5*ERFC(-X*SQRT(0.5))). When X is a CVX variable, a
%   a piecewise quadratic *approximation* is employed instead. This
%   approximation gives good results when -4 <= x <= 4, and will be
%   improved in future releases of CVX.
%
%   For array values of X, the LOG_NORMCDF returns an array of identical
%   size with the calculation applied independently to each element.
%
%   X must be real.
%
%   Disciplined convex programming information:
%       LOG_NORMCDF is concave and nondecreasing in X. Therefore, when used
%       in CVX specifications, X must be concave.

narginchk(1,2);
if ~isreal( x ),
    error( 'Argument must be real.' );
end
if nargin > 1,
    % For debugging purposes only
    y = cvx_constant(log_normcdf(cvx(x)));
else
    y = log(0.5*erfc(-x*sqrt(0.5)));
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function prepare_voc12_imgs(img_set, output_img_dir, opts)
% PREPARE_VOC12_IMGS
%   prepare voc12 images (cropped from ground truth bboxes with jittering)
% input:
%   img_set: 'train' or 'val'
%   output_img_dir: output folder names, final structure is <output_img_dir>/<category>/<imgs>
%   opts: matlab struct with flip, aug_n, jitter_IoU, difficult, truncated, occluded fields
%       if flip is 1, images will be flipped.
%       if ang_n>1, images will be augmented by jittering bbox. jitter_IoU==1 means normal crop
%       difficult, truncated, occluded \in {0,1}, where 0 indicated that we do not
%           want images with that property (e.g. 0,0,0 means we want easy images only)
% output:
%   cropped images according to ground-truth bounding boxes (with jittering) and image filelists
%

addpath(fullfile(mfilename('fullpath'), '../../'));
global_variables;
cls_names = g_cls_names;

% paths
annotation_path = fullfile(g_pascal3d_root_folder, 'Annotations');
image_path = fullfile(g_pascal3d_root_folder, 'Images');
addpath(fullfile(g_pascal3d_root_folder, 'VDPM'));
addpath(fullfile(g_pascal3d_root_folder, 'PASCAL/VOCdevkit/VOCcode'));

% read ids of train/val set images
VOCinit;
ids = textread(sprintf(VOCopts.imgsetpath, img_set), '%s');
M = numel(ids);

% avoid duplication
assert(exist(output_img_dir,'dir')==0);
mkdir(output_img_dir);
for cls_idx = 1:numel(cls_names)
    mkdir([output_img_dir '/' cls_names{cls_idx}]);
end

for cls_idx = 1:numel(cls_names)
    cls = cls_names{cls_idx};
    labelfile = fopen(fullfile(output_img_dir, sprintf('%s.txt',cls)),'w');

    for i = 1:M
        if mod(i,100) == 0
            fprintf('%s: %d/%d\n', cls, i, M);
        end

        anno_filename = fullfile(annotation_path, sprintf('%s_pascal/%s.mat', cls, ids{i}));
        if ~exist(anno_filename,'file')
            continue;
        end
        anno = load(anno_filename);
        objects = anno.record.objects;

        for k = 1:length(objects)
            obj = objects(k);

            if ~isempty(obj.viewpoint) && strcmp(obj.class, cls)
                try
                    % write view annotation
                    azimuth = mod(round(obj.viewpoint.azimuth), 360);
                    elevation = mod(round(obj.viewpoint.elevation), 360);
                    tilt = mod(round(obj.viewpoint.theta), 360);
                    truncated = obj.truncated;
                    occluded = obj.occluded;
                    difficult = obj.difficult;
                    % skip un-annotated image
                    if azimuth == 0 && elevation == 0 && theta == 0
                        fprintf('skip %s...', ids{i});
                        continue;
                    end
                    % skip unwanted image
                    if (difficult==1 && opts.difficult==0) || (truncated==1 && opts.truncated==0) || (occluded==1 && opts.occluded==0)
                        fprintf('fliter skip %s...', ids{i});
                        continue;
                    end

                    img_filename = fullfile(image_path, sprintf('%s_pascal/%s.jpg', cls, ids{i}));
                    im = imread(img_filename);
                    box = obj.bbox;
                    w = box(3)-box(1)+1;
                    h = box(4)-box(2)+1;
                    rect = [box(1),box(2), w, h];

                    for aug_i = 1:opts.aug_n
                        % write cropped img
                        if aug_i == 1
                            cropped_im = imcrop(im, rect);
                        else
                            cropped_im = jitter_imcrop(im, rect, opts.jitter_IoU);
                        end
                        cropped_im_filename = sprintf('%s/%s/%s_%s_%s_%s.jpg', output_img_dir, cls, cls, ids{i}, num2str(k), num2str(aug_i));
                        imwrite(cropped_im, cropped_im_filename);
                        fprintf(labelfile, '%s %d %d %d %d\n', cropped_im_filename, cls_idx-1, azimuth, elevation, tilt);
                        if opts.flip
                            cropped_im_flip = fliplr(cropped_im); % flip the image horizontally
                            cropped_im_flip_filename = sprintf('%s/%s/%s_%s_%s_%s_%s.jpg', output_img_dir, cls, cls, ids{i}, num2str(k), num2str(aug_i), 'flip');
                            imwrite(cropped_im_flip, cropped_im_flip_filename);
                            fprintf(labelfile, '%s %d %d %d %d\n', cropped_im_flip_filename, cls_idx-1, mod(360-azimuth,360), elevation, mod(-1*tilt,360));
                        end
                    end
                catch
                end
            end
        end
    end
    fclose(labelfile);
end

author: ShapeNet

repo: RenderForCNN


text:
function r=rotqr2ro(q)
%ROTQR2RO converts a real quaternion to a 3x3 rotation matrix
% Inputs:
%
%     Q(4,1)   real-valued quaternion (with magnitude = 1)
%
% Outputs:
%
%     R(3,3)   Input rotation matrix
%              Plots a diagram if no output specified
%
% In the quaternion representation of a rotation, and q(1) = cos(t/2)
% where t is the angle of rotation in the range 0 to 2pi
% and q(2:4)/sin(t/2) is a unit vector lying along the axis of rotation
% a positive rotation about [0 0 1] takes the X axis towards the Y axis.
%
%      Copyright (C) Mike Brookes 2007
%      Version: $Id: rotqr2ro.m,v 1.5 2008/12/03 09:53:15 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent a b c d e f g
if isempty(a)
    a=[1 5 9];
    b=[11 16 6];
    c=[16 6 11];
    d=[4 8 3];
    e=[10 15 14];
    f=[4 2 3];
    g=[2 6 7];
end
p=2*(q*q.')/(q.'*q);            % force normalized
r=zeros(3,3);
r(a)=1-p(b)-p(c);
r(d)=p(e)-p(f);
r(g)=p(e)+p(f);
if ~nargout
    % display rotated pyramid
    cla % clear current axis
%     vv=[0,0,0;1,0,0;0,1,0;0,0,1]*r';  % pyramid
%     ff=[1 2 4; 2 1 3; 3 1 4; 4 2 3];
%     cc=[0 1 0; 0 0 1; 1 0 0; 1 1 0];
    vv=[0,0,0;1,0,0;0,1,0;0,0,1;0 1 1; 1 0 1; 1 1 0; 1 1 1]*r';    % cube
    ff=[1 2 6 4; 2 7 8 6; 7 3 5 8; 4 5 3 1; 3 7 2 1; 6 8 5 4];
    cc=[0 1 0; 1 0 0; 0 1 0; 1 0 0; 0 0 1; 0 0 1];
    pa=patch('Vertices',vv,'Faces',ff,'FaceVertexCData',cc,'FaceColor','Flat');
    xlabel('x axis');
    ylabel('y axis');
    zlabel('z axis');
    title(sprintf('qr = [%.2f, %.2f, %.2f, %.2f]''    initial xyz=0 are rgb',q))
    axis([-1 1 -1 1 -1 1 0 1]*sqrt(3));
    grid on
    view(3);
    axis equal;
end

author: decouples


text:
filename='Tests_Triangle_Linear';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeterInterior'};
weights = [1 1];
constraint = {'volumeConstraint'};
%optimizer = 'DualNestedInPrimal';
optimizer = 'AlternatingPrimalDual';

optimizerUnconstrained = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'PDE';

nsteps = 1;
Vfrac_final = 0.3;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

plotting = false;
printing = false;
monitoring = false;

maxiter = 2;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% [psi] = dai_readfg (filename)
%
%    INPUT:  filename   = filename of a .fg file
%
%    OUTPUT: psi        = linear cell array containing the factors
%                         (psi{i} is a structure with a Member field
%                         and a P field, like a CPTAB).

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function out = EN_wentropy(y,whaten,p)
% EN_wentropy   Entropy of time series using wavelets.
%
% Uses the wentropy function from Matlab's Wavelet toolbox.
%
%--INPUTS:
% y, the input time series
% whaten, the entropy type:
%               'shannon',
%               'logenergy',
%               'threshold' (with a given threshold),
%               'sure' (with a given parameter).
%               (see the wentropy documentaiton for information)
% p, the additional parameter needed for threshold and sure entropies
%
%---NOTE:
% It seems likely that this implementation of wentropy is nonsense.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check that a Wavelet Toolbox license is available:
% ------------------------------------------------------------------------------
BF_CheckToolbox('wavelet_toolbox')

% ------------------------------------------------------------------------------
% Check inputs
% ------------------------------------------------------------------------------
if nargin < 2 || isempty(whaten)
    whaten = 'shannon'; % default
end

N = length(y); % time-series length

switch whaten
	case 'shannon' % Shannon entropy
		out = wentropy(y,'shannon')/N; % scales with N for large N

	case 'logenergy' % Log Energy entropy
		out = wentropy(y,'log energy')/N; % scales with N for large N

    case 'threshold' % Magnitude of the signal greater than some value
        out = wentropy(y,'threshold',p)/N;

    case 'sure'
        % Equivalent to threshold entropy?
        out = wentropy(y,'sure',p)/N;

    otherwise
        error('Unknown entropy type ''%s''', whaten);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function mse = EvalMSE(NN,tNN,sqi,HRVparams,windows_all)

% mse = EvalMSE(NN,tNN,sqi,HRVparams,windows_all)
%
%   OVERVIEW:   This function returns MultiScale Entropy calculated on
%               input NN intervals for each window.
%
%   INPUT:      MANDATORY:
%               NN             : a single row of NN (normal normal) interval
%                                data in seconds
%               tNN            : the time indices of the NN interval data
%                                (seconds)
%               sqi            : (Optional )Signal Quality Index; Requires
%                                a matrix with at least two columns. Column
%                                1 should be timestamps of each sqi measure,
%                                and Column 2 should be SQI on a scale from 0 to 1.
%               HRVparams      : struct of settings for hrv_toolbox analysis
%               windows_all    : vector containing the starting time of each
%                                windows (in seconds)
%
%   OUTPUT:
%               mse            : vector of [max_tau, 1] doubles for each
%                                window
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Giulia Da Poian
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%

% Verify input arguments
if nargin < 4
    error('Wrong number of input parameters');
end
if nargin < 5 || isempty(windows_all)
    windows_all = 0;
end
if isempty(sqi)
     sqi(:,1) = tNN;
     sqi(:,2) = ones(length(tNN),1);
end
% Set Defaults

if isempty(HRVparams.MSE.windowlength)
    windowlength = length(NN);
else
    windowlength = HRVparams.MSE.windowlength*3600;
end

threshold1 = HRVparams.sqi.LowQualityThreshold;
threshold2 = HRVparams.RejectionThreshold;

m = HRVparams.MSE.patternLength;
r = HRVparams.MSE.RadiusOfSimilarity;
maxTau = HRVparams.MSE.maxCoarseGrainings;

cg_moment = HRVparams.MSE.moment;
cg_method = HRVparams.MSE.method;
constant_r = HRVparams.MSE.constant_r;

% Preallocate arrays (all NaN) before entering the loop
mse = nan(maxTau,length(windows_all));

%Analyze by Window

% Loop through each window of RR data
for i_win = 1:length(windows_all)
    % Check window for sufficient data
    if ~isnan(windows_all(i_win))
        % Isolate data in this window
        idx_NN_in_win = find(tNN >= windows_all(i_win) & tNN < windows_all(i_win) + windowlength);
        idx_sqi_win = find(sqi(:,1) >= windows_all(i_win) & sqi(:,1) < windows_all(i_win) + windowlength);

        sqi_win = sqi(idx_sqi_win,:);
        nn_win = NN(idx_NN_in_win);

        % Analysis of SQI for the window
        lowqual_idx = find(sqi_win(:,2) < threshold1);

        % If enough data has an adequate SQI, perform the calculations
        if numel(lowqual_idx)/length(sqi_win(:,2)) < threshold2
            mse(:,i_win) = ComputeMultiscaleEntropy(nn_win, m, r, maxTau, [], cg_moment,cg_method,constant_r);
        end

    end % end check for sufficient data

end % end of loop through window

end % end of function

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
% This code just simply run the SVM on the example data set "heart_scale",
% which is scaled properly. The code divides the data into 2 parts
% train: 1 to 200
% test: 201:270
% Then plot the results vs their true class. In order to visualize the high
% dimensional data, we apply MDS to the 13D data and reduce the dimension
% to 2D

clear
clc
close all

% read the data set
[heart_scale_label, heart_scale_inst] = libsvmread('heart_scale');
[N,D] = size(heart_scale_inst);

% Determine the train and test index
trainIndex = zeros(N,1); trainIndex(1:200) = 1;
testIndex = zeros(N,1); testIndex(201:N) = 1;
trainData = heart_scale_inst(trainIndex==1,:);
trainLabel = heart_scale_label(trainIndex==1,:);
testData = heart_scale_inst(testIndex==1,:);
testLabel = heart_scale_label(testIndex==1,:);

% Train the SVM
model = svmtrain(trainLabel, trainData, '-c 1 -g 0.07 -b 1');
% Use the SVM model to classify the data
[predict_label, accuracy, prob_values] = svmpredict(testLabel, testData, model, '-b 1'); % run the SVM model on the test data

% ================================
% ===== Showing the results ======
% ================================

% Assign color for each class
% colorList = generateColorList(2); % This is my own way to assign the color...don't worry about it
colorList = prism(100);

% true (ground truth) class
trueClassIndex = zeros(N,1);
trueClassIndex(heart_scale_label==1) = 1;
trueClassIndex(heart_scale_label==-1) = 2;
colorTrueClass = colorList(trueClassIndex,:);
% result Class
resultClassIndex = zeros(length(predict_label),1);
resultClassIndex(predict_label==1) = 1;
resultClassIndex(predict_label==-1) = 2;
colorResultClass = colorList(resultClassIndex,:);

% Reduce the dimension from 13D to 2D
distanceMatrix = pdist(heart_scale_inst,'euclidean');
newCoor = mdscale(distanceMatrix,2);

% Plot the whole data set
x = newCoor(:,1);
y = newCoor(:,2);
patchSize = 30; %max(prob_values,[],2);
colorTrueClassPlot = colorTrueClass;
figure; scatter(x,y,patchSize,colorTrueClassPlot,'filled');
title('whole data set');

% Plot the test data
x = newCoor(testIndex==1,1);
y = newCoor(testIndex==1,2);
patchSize = 80*max(prob_values,[],2);
colorTrueClassPlot = colorTrueClass(testIndex==1,:);
figure; hold on;
scatter(x,y,2*patchSize,colorTrueClassPlot,'o','filled');
scatter(x,y,patchSize,colorResultClass,'o','filled');
% Plot the training set
x = newCoor(trainIndex==1,1);
y = newCoor(trainIndex==1,2);
patchSize = 30;
colorTrueClassPlot = colorTrueClass(trainIndex==1,:);
scatter(x,y,patchSize,colorTrueClassPlot,'o');
title('classification results');

author: ravenxrz

repo: Mathematical-Modeling


text:
function EUI_data=riceVST_EUI(Efz_data,sigma_data,VST_ABC)
% Applies exact unbiased inverse of the variance-stabilizing transformation f
% --------------------------------------------------------------------------------------------
%
% SYNTAX
% ------
% EUI_data = riceVST_EUI ( Efz_data , sigma_data , VST_ABC )
%
%
% OUTPUT
% ------
% EUI_data    :  exact unbiased inverse of the input Efz_data
%
%
% INPUTS
% ------
% Efz_data    :  Filtered variance-stabilized data (e.g, output of the denoising filter)
% sigma_data  :  standard-deviation parameter of the Rice distribution;  the method assumes
%                that the data z_data before stabilization is distributed according to
%                z_data ~ Rice(nu,sigma_data), nu being some unknown noise-free signal.
%                EUI_data should coincide with nu if the denoising of fz_data is perfect.
% VST_ABC     :  file-selector for the variance-stabilizing transformation  (default = 'A')
%
%
% --------------------------------------------------------------------------------------------
%
% The software implements the method published in the paper:
%
%  A. Foi, "Noise Estimation and Removal in MR Imaging: the Variance-Stabilization Approach",
%  in Proc. 2011 IEEE Int. Sym. Biomedical Imaging, ISBI 2011, Chicago (IL), USA, April 2011.
%  doi:10.1109/ISBI.2011.5872758
%
% --------------------------------------------------------------------------------------------
%
%
% author:                Alessandro Foi
%
% web page:              http://www.cs.tut.fi/~foi/RiceOptVST
%
% contact:               firstname.lastname@tut.fi
%
% --------------------------------------------------------------------------------------------
% Copyright (c) 2010-2012 Tampere University of Technology.
% All rights reserved.
% This work should be used for nonprofit purposes only.
% --------------------------------------------------------------------------------------------
%
% Disclaimer
% ----------
%
% Any unauthorized use of these routines for industrial or profit-oriented activities is
% expressively prohibited. By downloading and/or using any of these files, you implicitly
% agree to all the terms of the TUT limited license (included in the file Legal_Notice.txt).
% --------------------------------------------------------------------------------------------
%

%% Defaults
if nargin<2
    sigma_data=1;
end
if isempty(sigma_data)
    sigma_data=1;
end

if nargin<3
    VST_ABC='A';
end
if isempty(VST_ABC)
    VST_ABC='A';
end

if ischar(VST_ABC)
    if numel(VST_ABC)==1
        Rice_VST_matFile=['Rice_VST_',VST_ABC];
    else
        Rice_VST_matFile=VST_ABC;
    end
end

%% load variance-stabilizing transformation and its exact unbiased inverse from file
load(Rice_VST_matFile,'Efz','nu','z','f')

%% scale data before applying exact unbiased inverse of the variance-stabilizing transformation (see riceVST.m)
sigma_data_scaling=sigma_data;
a=f(end)-sqrt(sigma_data_scaling^2*max(z)^2/sigma_data^2-1/2);
nu=nu*sigma_data_scaling;
%% apply exact unbiased inverse of the variance-stabilizing transformation
EUI_data=interp1(Efz,nu,Efz_data,'linear','extrap');
%% small values (see Equation 11 in the ISBI2011 paper)
EUI_data(Efz_data<=min(Efz))=min(nu);
%% apply asymptotical exact unbiased inverse of the variance-stabilizing transformation (used only for large values)
maxEfz=max(Efz);
if 0
    EUI_data(Efz_data>maxEfz)=sigma_data*sqrt((Efz_data(Efz_data>maxEfz)-a).^2+0.5);  %% this is the asymptotic inversion to Ez (identical to the above theta_or_Ez==1 case)
    EUI_data(Efz_data>maxEfz)=EUI_data(Efz_data>maxEfz).*(1-0.5*(sigma_data./EUI_data(Efz_data>maxEfz)).^2);  %% this is the asymptotic correction between nu and Ez
else
    EUI_data(Efz_data>maxEfz)=sigma_data*((Efz_data(Efz_data>maxEfz)-a).^2)./sqrt((Efz_data(Efz_data>maxEfz)-a).^2+0.5);   %% this is the result of the two lines above
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% GAINSANDREFERENCES compute gains matrices, references and regularization
%                    parameters for the controller for a specific robot

%% --- Initialization ---

% references for the CoM (sine trajectory)
Config.referencesCoM.amplitude = 0.02; %[m]
Config.referencesCoM.frequency = 0.25; %[Hz]
Config.referencesCoM.direction = [0;1;0];

% pseudoinverse tolerances
Config.Reg.pinv_tol        = 1e-4;
Config.Reg.pinvDampBaseVel = 1e-4;

% Gains matrices

% CoM gains
Config.Gain.KP_CoM    = diag([50;50;50]);
Config.Gain.KD_CoM    = 2*sqrt(Config.Gain.KP_CoM);

% Feet pose Gains
Config.Gain.KP_feet   = diag([50;50;50]);
Config.Gain.KD_feet   = 2*sqrt(Config.Gain.KP_feet);

% Postural task gains         % torso    % lArm       % rArm       % lLeg             % rLeg
Config.Gain.KP_joints = diag([20;20;20; 10;10;10;10; 10;10;10;10; 30;30;20;30;40;40; 30;30;20;30;40;40]);
Config.Gain.KD_joints = 2*sqrt(Config.Gain.KP_joints);

author: robotology

repo: whole-body-controllers


text:
function [p2_closest, ip2_closest, dmin] = nearest_point(p2, p1, ith_closest)
%
% Usage:
%
%    [p2_closest ip2_closest] = nearest_point(p2,p1)
%
% Description:
%
%    Find the points in p2 that are nearest to the points in p1.
%
% AUTHOR: Jay Dubb (jdubb@nmr.mgh.harvard.edu)
% DATE:   11/13/2008
%
% Updates:
%
% 01/21/2010 - Jay Dubb changed variable names and modified comment
%              to better describe the way this function works.
%

% Output arguments
p2_closest  = [];
ip2_closest = [];
dmin        = 0;

if isempty(p1) || isempty(p2)
    return;
end
if ~exist('ith_closest','var') | isempty(ith_closest)
    ith_closest = 1;
end

% Figure out the dimentions of p1 and p2 to see if they're compatible
sz1 = size(p1);
sz2 = size(p2);
ndim = max(sz1(sz1==sz2));
if isempty(ndim)
    if length(sz1)>2 || length(sz1)~=length(sz2)
        MenuBox('nearest_point error: p1 and p2 number of dimensions incompatible','OK');
        return;
    end
    if length(sz1)>2 || length(sz1)~=length(sz2)
        MenuBox('nearest_point error: p1 and p2 number of dimensions incompatible','OK');
        return;
    end
    if ismember(1,sz1) && ismember(1,sz2)
        ndim=1;
    end
end
if ndim~=3 && ndim~=1
    MenuBox('nearest_point warning: currently only handles arguments in 1 and 3 dimensions.','OK');
    return;
end

m=size(p1,1);
n=size(p2,1);

if ith_closest>n
    fprintf('Error: %dth closest element is greater than size of first argument (%d)\n', ith_closest, n);
    return;
end

if ndim==3
    p2_closest  = zeros(m,3);
    ip2_closest = zeros(m,1);
    dmin = zeros(m,1);
    for k=1:m
        d = sqrt((p2(:,1)-p1(k,1)).^2+(p2(:,2)-p1(k,2)).^2+(p2(:,3)-p1(k,3)).^2);
        [d2, j] = sort(d);
        dmin(k) = d2(ith_closest);

        ip2_closest(k)  = j(ith_closest);
        p2_closest(k,:) = p2(ip2_closest(k),:);
    end
elseif ndim==1
    p2_closest  = zeros(m,1);
    ip2_closest = zeros(m,1);
    dmin = zeros(m,1);
    for k=1:m
        d = sqrt((p2(:)-p1(k)).^2);
        [d2, j] = sort(d);
        dmin(k) = d2(ith_closest);

        ip2_closest(k) = j(ith_closest);
        p2_closest(k)  = p2(ip2_closest(k));
    end
end

author: BUNPC


text:
function deform_GUI(V,F,C,W)
%%%%%%%%%%
% you dont  need

fprintf( ...
    ['Linear Blend Skinning: \n' ...
    '- CLICK handles to visualize weights \n', ...
    '- DRAG a handle to move\n', ...
    '- SHIFT+DRAG a handle to rotate\n'] ...
    );

P = 1:size(C,1); % set default point handles
if(size(C,2) == 3)
C = C(:,1:2); % Be sure that control vertices are in 2D
end
CW = W; % Weights used for contours
WVW = W;  % Weights used for weight visualization
np = numel(P);  % number of point handles

DeformOBJ.new_C = [];
% keep track of rotations stored at each control point, for 2D this is a m
% by 1 list of angles
DeformOBJ.R = zeros(np,1);
DeformOBJ.update_positions = @update_positions;

subplot(1,2,1)
DeformOBJ.tsh = tsurf(F,V);
axis equal
hold on
C_plot = scatter3( ...
    C(:,1),C(:,2),0.1+0*C(:,1), ...
    'o','MarkerFaceColor',[0.9 0.8 0.1], 'MarkerEdgeColor','k',...
    'LineWidth',2,'SizeData',100, ...
    'ButtonDownFcn',@oncontrolsdown);
hold off;
% axis manual

subplot(1,2,2)
% plot the original mesh
DeformOBJ.wvsh = tsurf(F,V);
view(2);
axis equal
% axis manual

% keep track of window xmin, xmax, ymin, ymax
win_min = min([C(:,1:2); V(:,1:2)]);
win_max = max([C(:,1:2); V(:,1:2)]);
% keep track of down position
down_pos = [];
% keep track of last two drag positions
drag_pos = [];
last_drag_pos = [];
% keep track of mesh vertices at mouse down
down_V = [];
% keep track of index of selected control point
ci = [];
% type of click ('left','right')
down_type  = '';

% Callback for mouse down on control points
  function oncontrolsdown(src,ev)
    % get current mouse position, and remember old one
    down_pos=get(gca,'currentpoint');
    down_pos=[down_pos(1,1,1) down_pos(1,2,1)];
    last_drag_pos=down_pos;
    drag_pos=down_pos;
    % keep track of control point positions at mouse down
    DeformOBJ.new_C = [get(C_plot,'XData')' get(C_plot,'YData')'];
    % get index of closest control point
    [minD,ci] =  ...
      min(sum((DeformOBJ.new_C(:,1:2) - ...
      repmat(down_pos,size(DeformOBJ.new_C,1),1)).^2,2));
    % keep track of mesh vertices at mouse down
    down_V = get(DeformOBJ.tsh,'Vertices');
    down_V = down_V(:,1:2);

    % tell window that drag and up events should be handled by controls
    set(gcf,'windowbuttonmotionfcn',@oncontrolsdrag)
    set(gcf,'windowbuttonupfcn',@oncontrolsup)
    set(gcf,'KeyPressFcn',@onkeypress)
    if(strcmp('normal',get(gcf,'SelectionType')))
      % left-click
      down_type = 'left';
    else
      % other (right) click
      down_type = 'right';
    end

    % try to find ci in list of point handles
    [found, iP] = ismember(ci,P);
    if(found)
      % set color of mesh plot to weights of selected
      %set(tsh,'CData',W(:,iP));
      % change weights in weight visualization
      set(DeformOBJ.wvsh,'CData',WVW(:,iP));
    end

  end

  % Callback for mouse drag on control points
  function oncontrolsdrag(src,ev)
    % keep last drag position
    last_drag_pos=drag_pos;
    % get current mouse position
    drag_pos=get(gca,'currentpoint');
    drag_pos=[drag_pos(1,1,1) drag_pos(1,2,1)];
    if(strcmp('left',down_type))
      % move selected control point by drag offset
      DeformOBJ.new_C(ci,:) = ...
        DeformOBJ.new_C(ci,:) + drag_pos-last_drag_pos;
    else
      [found, iP] = ismember(ci,P);
      if(found)
        DeformOBJ.R(iP) = ...
          DeformOBJ.R(iP) + 2*pi*(drag_pos(1)-last_drag_pos(1))/100;
      end
    end
    update_positions();
  end

  function update_positions()
    % update display positions
    set(C_plot,'XData',DeformOBJ.new_C(:,1));
    set(C_plot,'YData',DeformOBJ.new_C(:,2));

    % USING LINEAR BLEND SKINNING
    % get transformations stored at each point and bone handle
    TR = ...
    skinning_transformations(C,P,[],DeformOBJ.new_C,DeformOBJ.R);

    % linear blend skinning
    [new_V] = linear_blend_skinning(V(:,1:2),TR,W);

    % update mesh positions
    set(DeformOBJ.tsh,'Vertices',new_V(:,1:2));
  end

  % Callback for mouse release of control points
  function oncontrolsup(src,ev)
    % Tell window to handle drag and up events itself
    set(gcf,'windowbuttonmotionfcn','');
    set(gcf,'windowbuttonupfcn','');
    cur_V = get(DeformOBJ.tsh,'Vertices');
    cur_V = cur_V(:,1:2);

    % scale window to fit
    win_min = min([win_min; cur_V]);
    win_max = max([win_max; cur_V]);
    axis(reshape([win_min;win_max],1,2*size(cur_V,2)))
  end

  function onkeypress(src,ev)
    if(strcmp(ev.Character,'r'))
      DeformOBJ.new_C = C;
      DeformOBJ.R = zeros(np,1);
      update_positions();
    elseif(strcmp(ev.Character,'u'))
      update_positions();
    end
  end
end

author: odedstein

repo: sgi-introduction-course


text:
function g = grad(SO3F,varargin)
% right-sided gradient of an SO3Fun
%
% Syntax
%   G = SO3F.grad % compute the gradient
%   g = SO3F.grad(rot) % evaluate the gradient in rot
%
%   % go 5 degree in direction of the gradient
%   ori_new = exp(rot,5*degree*normalize(g))
%
% Input
%  SO3F - @SO3FunCBF
%  rot  - @rotation / @orientation
%
% Output
%  G - @SO3VectorField
%  g - @vector3d
%
% Description
% general formula:
%
% $$s(g1_i) = sum_j c_j DRK(<g h_j,r_j>) g h_j x r_j $$
%

% fallback to generic method
if check_option(varargin,'check') || nargin == 1 || ~isa(varargin{1},'rotation')
  g = grad@SO3Fun(SO3F,varargin{:});
  return
end

ori = varargin{1};

% symmetrise - only crystal symmetry
[h,l] = symmetrise(SO3F.h.normalize,'unique');
r = repelem(SO3F.r.normalize,l);
w = repelem(SO3F.weights./l,l);

g = vector3d.zeros(size(ori));
for i = 1:length(h)
  g = g + w(i) * SO3F.psi.grad(dot(ori*h(i),r(i),'noSymmetry'),'polynomial') .* ...
      cross(h(i),inv(ori) * r(i));
end
end

function test

cs = crystalSymmetry('1');
odf = fibreODF(Miller(0,0,1,cs),vector3d.Z);
omega = linspace(-20,20)*degree;
omega = 15 *degree;
ref = orientation.byAxisAngle(vector3d(1,0,10),omega,cs)

g1 = odf.grad(ref)
g2 = odf.grad(ref,'check','delta',0.05*degree)

plot(omega./degree,[g1.x,g2.x])

omega2 = linspace(-5,5)*degree;
ori1 = ref * rotation.byAxisAngle(g1,omega2);
ori2 = ref * rotation.byAxisAngle(g2,omega2);
ori3 = ref * rotation.byAxisAngle(normalize(g1+g2),omega2);
%ori4 = ref * rotation.byAxisAngle(normalize(g1-g2),'angle',omega2);
ori4 = ref * rotation.byAxisAngle(vector3d(-1,2,0),omega2);

plot(omega2./degree,[odf.eval(ori1(:)),odf.eval(ori2(:)),odf.eval(ori3(:)),odf.eval(ori4(:))])

end

function test2

  cs = crystalSymmetry('321');
  odf = fibreODF(Miller(1,2,3,cs),vector3d(-1,3,2));

  ref = orientation.rand(1000,cs)

  g1 = odf.grad(ref)
  g2 = odf.grad(ref,'check','delta',0.05*degree)

  hist(norm(g1-g2)./degree)

  ref = orientation.id(cs) * cs(5);

  f = S2FunHarmonic.quadrature(@(r) ...
    odf.eval(rotation.byAxisAngle(r,5*degree)*ref));

  plot(f,'lower')

  annotate(odf.grad(ref))
end

author: mtex-toolbox

repo: mtex

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


text:
function w = faddeeva(z,N)
% FADDEEVA   Faddeeva function
%   W = FADDEEVA(Z) is the Faddeeva function, aka the plasma dispersion
%   function, for each element of Z. The Faddeeva function is defined as:
%
%     w(z) = exp(-z^2) * erfc(-j*z)
%
%   where erfc(x) is the complex complementary error function.
%
%   W = FADDEEVA(Z,N) can be used to explicitly specify the number of terms
%   to truncate the expansion (see (13) in [1]). N = 16 is used as default.
%
%   Example:
%       x = linspace(-10,10,1001); [X,Y] = meshgrid(x,x);
%       W = faddeeva(complex(X,Y));
%       figure;
%       subplot(121); imagesc(x,x,real(W)); axis xy square; caxis([-1 1]);
%       title('re(faddeeva(z))'); xlabel('re(z)'); ylabel('im(z)');
%       subplot(122); imagesc(x,x,imag(W)); axis xy square; caxis([-1 1]);
%       title('im(faddeeva(z))'); xlabel('re(z)'); ylabel('im(z)');
%
%   Reference:
%   [1] J.A.C. Weideman, "Computation of the Complex Error Function," SIAM
%       J. Numerical Analysis, pp. 1497-1518, No. 5, Vol. 31, Oct., 1994
%       Available Online: http://www.jstor.org/stable/2158232

if nargin<2, N = []; end
if isempty(N), N = 16; end

w = zeros(size(z)); % initialize output

%%%%%
% for purely imaginary-valued inputs, use erf as is if z is real
idx = real(z)==0; %
w(idx) = exp(-z(idx).^2).*erfc(imag(z(idx)));

if all(idx), return; end
idx = ~idx;

%%%%%
% for complex-valued inputs

% make sure all points are in the upper half-plane (positive imag. values)
idx1 = idx & imag(z)<0;
z(idx1) = conj(z(idx1));

M = 2*N;
M2 = 2*M;
k = (-M+1:1:M-1)'; % M2 = no. of sampling points.
L = sqrt(N/sqrt(2)); % Optimal choice of L.

theta = k*pi/M;
t = L*tan(theta/2); % Variables theta and t.
f = exp(-t.^2).*(L^2+t.^2);
f = [0; f]; % Function to be transformed.
a = real(fft(fftshift(f)))/M2; % Coefficients of transform.
a = flipud(a(2:N+1)); % Reorder coefficients.

Z = (L+1i*z(idx))./(L-1i*z(idx));
p = polyval(a,Z); % Polynomial evaluation.
w(idx) = 2*p./(L-1i*z(idx)).^2 + (1/sqrt(pi))./(L-1i*z(idx)); % Evaluate w(z).

% convert the upper half-plane results to the lower half-plane if necesary
w(idx1) = conj(2*exp(-z(idx1).^2) - w(idx1));


text:
function [ v, lambda ] = l1nd_eigen ( n, h )

%*****************************************************************************80
%
%% L1ND_EIGEN returns eigeninformation for the 1D ND Laplacian.
%
%  Discussion:
%
%    The grid points are assumed to be evenly spaced by H.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real H, the spacing between points.
%
%    Output, real V(N,N), the eigenvectors.
%
%    Output, real LAMBDA(N), the eigenvalues.
%
  v = zeros ( n, n );
  lambda = zeros ( n, 1 );

  for j = 1 : n
    theta = pi * ( j - 0.5 ) / ( 2.0 * n + 1.0 );
    lambda(j) = 4.0 * ( sin ( theta ) / h ) ^ 2;
    for i = 1 : n
      theta = pi * ( i - 0.5 ) * ( 2.0 * j - 1.0 ) / ( 2.0 * n + 1.0 );
      v(i,j) = sqrt ( 2.0 / ( n + 0.5 ) ) * cos ( theta );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [M, L, Y, C] = lmnn(X, labels)
%LMNN Learns a metric using large-margin nearest neighbor metric learning
%
%   [M, L, Y, C] = lmnn(X, labels)
%
% The function uses large-margin nearest neighbor (LMNN) metric learning to
% learn a metric on the data set specified by the NxD matrix X and the
% corresponding Nx1 vector labels. The metric is returned in M.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    % Initialize some variables
    [N, D] = size(X);
    assert(length(labels) == N);
    [lablist, ~, labels] = unique(labels);
    K = length(lablist);
    label_matrix = false(N, K);
    label_matrix(sub2ind(size(label_matrix), (1:length(labels))', labels)) = true;
    same_label = logical(double(label_matrix) * double(label_matrix'));
    M = eye(D);
    C = Inf; prev_C = Inf;

    % Set learning parameters
    min_iter = 50;          % minimum number of iterations
    max_iter = 1000;        % maximum number of iterations
    eta = .1;               % learning rate
    mu = .5;                % weighting of pull and push terms
    tol = 1e-3;             % tolerance for convergence
    best_C = Inf;           % best error obtained so far
    best_M = M;             % best metric found so far
    no_targets = 3;         % number of target neighbors

    % Select target neighbors
    sum_X = sum(X .^ 2, 2);
    DD = bsxfun(@plus, sum_X, bsxfun(@plus, sum_X', -2 * (X * X')));
    DD(~same_label) = Inf; DD(1:N + 1:end) = Inf;
    [~, targets_ind] = sort(DD, 2, 'ascend');
    targets_ind = targets_ind(:,1:no_targets);
    targets = false(N, N);
    targets(sub2ind([N N], vec(repmat((1:N)', [1 no_targets])), vec(targets_ind))) = true;

    % Compute pulling term between target neigbhors to initialize gradient
    slack = zeros(N, N, no_targets);
    G = zeros(D, D);
    for i=1:no_targets
        G = G + (1 - mu) .* (X - X(targets_ind(:,i),:))' * (X - X(targets_ind(:,i),:));
    end

    % Perform main learning iterations
    iter = 0;
    while (prev_C - C > tol || iter < min_iter) && iter < max_iter

        % Compute pairwise distances under current metric
        XM = X * M;
        sum_X = sum(XM .* X, 2);
        DD = bsxfun(@plus, sum_X, bsxfun(@plus, sum_X', -2 * (XM * X')));

        % Compute value of slack variables
        old_slack = slack;
        for i=1:no_targets
            slack(:,:,i) = ~same_label .* max(0, bsxfun(@minus, 1 + DD(sub2ind([N N], (1:N)', targets_ind(:,i))), DD));
        end

        % Compute value of cost function
        prev_C = C;
        C = (1 - mu) .* sum(DD(targets)) + ...  % push terms between target neighbors
                 mu  .* sum(slack(:));          % pull terms between impostors

        % Maintain best solution found so far (subgradient method)
        if C < best_C
            best_C = C;
            best_M = M;
        end

        % Perform gradient update
        for i=1:no_targets

            % Add terms for new violations
            [r, c] = find(slack(:,:,i) > 0 & old_slack(:,:,i) == 0);
            G = G + mu .* ((X(r,:) - X(targets_ind(r, i),:))' * ...
                           (X(r,:) - X(targets_ind(r, i),:)) - ...
                           (X(r,:) - X(c,:))' * (X(r,:) - X(c,:)));

            % Remove terms for resolved violations
            [r, c] = find(slack(:,:,i) == 0 & old_slack(:,:,i) > 0);
            G = G - mu .* ((X(r,:) - X(targets_ind(r, i),:))' * ...
                           (X(r,:) - X(targets_ind(r, i),:)) - ...
                           (X(r,:) - X(c,:))' * (X(r,:) - X(c,:)));
        end
        M = M - (eta ./ N) .* G;

        % Project metric back onto the PSD cone
        [V, L] = eig(M);
        V = real(V); L = real(L);
        ind = find(diag(L) > 0);
        if isempty(ind)
            warning('Projection onto PSD cone failed. All eigenvalues were negative.'); break
        end
        M = V(:,ind) * L(ind, ind) * V(:,ind)';
        if any(isinf(M(:)))
            warning('Projection onto PSD cone failed. Metric contains Inf values.'); break
        end
        if any(isnan(M(:)))
            warning('Projection onto PSD cone failed. Metric contains NaN values.'); break
        end

        % Update learning rate
        if prev_C > C
            eta = eta * 1.01;
        else
            eta = eta * .5;
        end

        % Print out progress
        iter = iter + 1;
        no_slack = sum(slack(:) > 0);
        if rem(iter, 10) == 0
            [~, sort_ind] = sort(DD, 2, 'ascend');
            disp(['Iteration ' num2str(iter) ': error is ' num2str(C ./ N) ...
                  ', nearest neighbor error is ' num2str(sum(labels(sort_ind(:,2)) ~= labels) ./ N) ...
                  ', number of constraints: ' num2str(no_slack)]);
        end
    end

    % Return best metric and error
    M = best_M;
    C = best_C;

    % Compute mapped data
    [L, S, ~] = svd(M);
    L = bsxfun(@times, sqrt(diag(S)), L);
    Y = X * L;
end

function x = vec(x)
    x = x(:);
end

repo: cluster


text:
function model= dbnFit(X, numhid, y, varargin)
%fit a DBN to bianry data in X

%INPUTS:
%X              ... data. should be binary, or in [0,1] interpreted as
%               ... probabilities
%numhid         ... list of numbers of hidden units
%y              ... List of discrete labels

%OUTPUTS:
%model          ... A cell array containing models from all RBM's

%varargin may contain options for the RBM's of this DBN, in row one by one
%for example:
%dbnFit(X, [500,400], opt1, opt2) uses opt1 for 500 and opt2 for 400
%dbnFit(X, [500,400], opt1) uses opt1 only for 500, and defaults for 400

numopts=length(varargin);
H=length(numhid);
model=cell(H,1);
cant_clases = length(unique(y));

if H>=2

    %train the first RBM on data
    if(numopts>=1)
        model{1}= rbmBB(X, numhid(1),varargin{1});
    else
        model{1}= rbmBB(X, numhid(1));
    end

%     %extra-debug
%     visualize(model{1}.W);
%     drawnow;

    %train all other RBM's on top of each other
    for ii=2:H
        if(numopts>=ii)
            model{ii}=rbmBB(model{ii-1}.top, numhid(ii), varargin{ii});
        else
            model{ii}=rbmBB(model{ii-1}.top, numhid(ii));
        end

%         %extra-debug
%         visualize(model{ii}.W);
%         drawnow;

    end

%     %the last RBM has access to labels too
%     if(numopts>=H)
%         model{H}= rbmFit(model{H-1}.top, numhid(end), y, varargin{H});
%     else
%         model{H}= rbmFit(model{H-1}.top, numhid(end), y);
%     end
    model{H}.Wc = 0.1*randn(cant_clases, size(model{H}.W,2) );
    model{H}.cc = 0.1*randn(1, cant_clases);
    model{H}.labels = 1:cant_clases;

    %fine-tuning of the DBN

    targets = zeros(length(y), cant_clases);
    for ii = 1:cant_clases
        targets(y == ii, ii) = 1;
    end

    for epoch = 1:200

        %%%%%%%%%%%%%%% PERFORM CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        max_iter=3;

        if epoch<6  % First update top-level weights holding other weights fixed.
%             N = size(data,1);
%             XX = [data ones(N,1)];
%             w1probs = 1./(1 + exp(-XX*w1)); w1probs = [w1probs  ones(N,1)];
%             w2probs = 1./(1 + exp(-w1probs*w2)); w2probs = [w2probs ones(N,1)];
%             w3probs = 1./(1 + exp(-w2probs*w3)); %w3probs = [w3probs  ones(N,1)];

            testdata = X;
            for ii=1:length(model)
                testdata= rbmVtoH(model{ii}, testdata);
            end

            VV = colvec([ model{H}.Wc'; model{H}.cc ]);
            Dim = [ numhid(end); cant_clases];
            XX = minimize(VV,'CG_CLASSIFY_INIT',max_iter,Dim,testdata,targets);
            aux = reshape(XX, numhid(end)+1, cant_clases);
            model{H}.Wc = aux(1:end-1,:)';
            model{H}.cc = aux(end,:);

        else

%             VV = [w1(:)' w2(:)' w3(:)' w_class(:)']';
%             Dim = [l1; l2; l3; l4; l5];
            VV = [];
            for ii=1:H
                VV = [VV; colvec([model{ii}.W; model{ii}.b ]) ];
            end
            VV = [VV; colvec([model{ii}.Wc'; model{ii}.cc ])];

            Dim = [ size(X,2); colvec(numhid); cant_clases];

            XX = minimize(VV,'CG_CLASSIFY',max_iter,Dim,X,targets);

            xxx = 0;
            for ii=1:length(model)
                aux = reshape(XX(xxx+1:xxx+(Dim(ii)+1)*Dim(ii+1)),Dim(ii)+1,Dim(ii+1));
                model{ii}.W = aux(1:end-1,:);
                model{ii}.b = aux(end,:);
                xxx = xxx+(Dim(ii)+1)*Dim(ii+1);
            end
            aux = reshape(XX(xxx+1:xxx+(Dim(ii+1)+1)*Dim(ii+2)),Dim(ii+1)+1,Dim(ii+2));
            model{ii}.Wc = aux(1:end-1,:)';
            model{ii}.cc = aux(end,:);

%             w1 = reshape(X(1:(l1+1)*l2),l1+1,l2);
%             xxx = (l1+1)*l2;
%             w2 = reshape(X(xxx+1:xxx+(l2+1)*l3),l2+1,l3);
%             xxx = xxx+(l2+1)*l3;
%             w3 = reshape(X(xxx+1:xxx+(l3+1)*l4),l3+1,l4);
%             xxx = xxx+(l3+1)*l4;
%             w_class = reshape(X(xxx+1:xxx+(l4+1)*l5),l4+1,l5);

        end
        %%%%%%%%%%%%%%% END OF CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        prediction = dbnPredict(model, X);
        errors = sum(prediction~=y);
        fprintf(1, 'Errores %d\n', errors);

    end

%     %extra-debug
%     visualize(model{H}.W);
%     drawnow;

else

    %numhid is only a single layer... but we should work anyway
    if (numopts>=1)
        model{1}= rbmFit(X, numhid(1), y, varargin{1});
    else
        model{1}= rbmFit(X, numhid(1), y);
    end
end

author: marianux

repo: ecg-kit

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


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrpmht
%TFRPMHT Unit test for the time frequency representation TFRPMH.

%       O. Lemoine - March 1996.

% We test each property of the corresponding TFR :

N=128;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrpmh(sig1);
tfr2=tfrpmh(sig2);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrpmh test 1 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrpmh(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrpmh test 2 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrpmh(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrpmh test 3 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrpmh(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrpmh test 4 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(N/4,1);noisecg(N/2);zeros(N/4,1)];
tfr=tfrpmh(sig);
if sum(any(abs(tfr(:,1:N/4-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(3*N/4+1):N))>sqrt(eps))),
 error('tfrpmh test 5 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrpmh(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrpmh test 6 failed');
end;

% A PMHD with a constant window is a MHD
sig=fmlin(N).*amgauss(N);
tfr1=tfrpmh(sig,1:N,N,ones(2*N+1,1));
tfr2=tfrmh(sig);
if max(max(abs(tfr1-tfr2)))>1e-5,
 error('tfrpmh test 7 failed');
end

N=131;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrpmh(sig1);
tfr2=tfrpmh(sig2);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrpmh test 8 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrpmh(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrpmh test 9 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrpmh(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrpmh test 10 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrpmh(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrpmh test 11 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(round(N/4),1);noisecg(round(N/2));zeros(round(N/4),1)];
tfr=tfrpmh(sig);
if sum(any(abs(tfr(:,1:round(N/4)-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(round(3*N/4)+2):N))>sqrt(eps))),
 error('tfrpmh test 12 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrpmh(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrpmh test 13 failed');
end;

% A PMHD with a constant window is a MHD
sig=fmlin(N).*amgauss(N);
tfr1=tfrpmh(sig,1:N,N,ones(2*N+1,1));
tfr2=tfrmh(sig);
if max(max(abs(tfr1-tfr2)))>1e-5,
 error('tfrpmh test 14 failed');
end

author: HeLiangHIT