text:
function varargout = hopfieldNetwork(varargin)
% HOPFIELDNETWORK M-file for hopfieldNetwork.fig
%      HOPFIELDNETWORK, by itself, creates a new HOPFIELDNETWORK or raises the existing
%      singleton*.
%
%      H = HOPFIELDNETWORK returns the handle to a new HOPFIELDNETWORK or the handle to
%      the existing singleton*.
%
%      HOPFIELDNETWORK('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in HOPFIELDNETWORK.M with the given input arguments.
%
%      HOPFIELDNETWORK('Property','Value',...) creates a new HOPFIELDNETWORK or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before hopfieldNetwork_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to hopfieldNetwork_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help hopfieldNetwork

% Last Modified by GUIDE v2.5 21-Jan-2007 15:45:38

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

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

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

% Choose default command line output for hopfieldNetwork
handles.output = hObject;
N = str2num(get(handles.imageSize,'string'));
handles.W = [];
handles.hPatternsDisplay = [];

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes on button press in reset.
function reset_Callback(hObject, eventdata, handles)
% cleans all data and enables the change of the number of neurons used
    for n=1 : length(handles.hPatternsDisplay)
        delete(handles.hPatternsDisplay(n));
    end
    handles.hPatternsDisplay = [];
    set(handles.imageSize,'enable','on');
    handles.W = [];
    guidata(hObject, handles);

function imageSize_Callback(hObject, eventdata, handles)
% hObject    handle to imageSize (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    num = get(hObject,'string');
    n = str2num(num);
    if isempty(n)
        num = '32';
        set(hObject,'string',num);
    end
    if n > 32
        warndlg('It is strongly recomended NOT to work with networks with more then 32^2 neurons!','!! Warning !!')
    end

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

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

% --- Executes on button press in loadIm.
function loadIm_Callback(hObject, eventdata, handles)
    [fName dirName] = uigetfile('*.bmp;*.tif;*.jpg;*.tiff');
    if fName
        set(handles.imageSize,'enable','off');
        cd(dirName);
        im = imread(fName);
        N = str2num(get(handles.imageSize,'string'));
        im = fixImage(im,N);
        imagesc(im,'Parent',handles.neurons);
        colormap('gray');
    end

% --- Executes on button press in train.
function train_Callback(hObject, eventdata, handles)

    Npattern = length(handles.hPatternsDisplay);
    if Npattern > 9
        msgbox('more then 10 paterns isn''t supported!','error');
        return
    end

    im = getimage(handles.neurons);
    N = get(handles.imageSize,'string');
    N = str2num(N);
    W = handles.W;  %weights vector
    avg = mean(im(:));   %removing the cross talk part
    if ~isempty(W)
        %W = W +( kron(im,im))/(N^2);
        W = W + ( kron(im-avg,im-avg))/(N^2)/avg/(1-avg);
    else
        % W = kron(im,im)/(N^2);
        W = ( kron(im-avg,im-avg))/(N^2)/avg/(1-avg);
    end
    % Erasing self weight
    ind = 1:N^2;
    f = find(mod(ind,N+1)==1);
    W(ind(f),ind(f)) = 0;

    handles.W = W;

    % Placing the new pattern in the figure...
    xStart = 0.01;
    xEnd = 0.99;
    height = 0.65;
    width = 0.09;
    xLength = xEnd-xStart;
    xStep = xLength/10;
    offset = 4-ceil(Npattern/2);
    offset = max(offset,0);
    y = 0.1;

    if Npattern > 0
        for n=1 : Npattern
            x = xStart+(n+offset-1)*xStep;
            h = handles.hPatternsDisplay(n);
            set(h,'units','normalized');
            set(h,'position',[x y width height]);
        end
        x = xStart+(n+offset)*xStep;
        h = axes('units','normalized','position',[x y width height]);
        handles.hPatternsDisplay(n+1) = h;
        imagesc(im,'Parent',h);
    else
        x = xStart+(offset)*xStep;
        h = axes('units','normalized','position',[x y width height]);
        handles.hPatternsDisplay = h;
    end

    imagesc(im,'Parent',h);
    set(h, 'YTick',[],'XTick',[],'XTickMode','manual','Parent',handles.learnedPaterns);

    guidata(hObject, handles);

% --- Executes on button press in addNoise.
function addNoise_Callback(hObject, eventdata, handles)
    im = getimage(handles.neurons);
    % N = get(handles.imageSize,'string');
    % N = floor(str2num(N)/2)+1;
    noisePercent = get( handles.noiseAmount, 'value' );
    N = round( length(im(:))* noisePercent );
    N = max(N,1);   %minimum change one neuron
    ind = ceil(rand(N,1)*length(im(:)));
%    im(ind) = -1*im(ind); %!!!!
    im(ind) = ~im(ind);
    imagesc(im,'Parent',handles.neurons);
    colormap('gray');

% --- Executes on button press in run.
function run_Callback(hObject, eventdata, handles)
    im = getimage(handles.neurons);
    [rows cols] = size(im);
    if rows ~= cols
        msgbox('I don''t support non square images','error');
        return;
    end
    N = rows;
    W = handles.W;
    if isempty(W)
        msgbox('No train data - doing nothing!','error');
        return;
    end
    %figure; imagesc(W)
    mat = repmat(im,N,N);
    mat = mat.*W;
    mat = im2col(mat,[N,N],'distinct');
    networkResult = sum(mat);
    networkResult = reshape(networkResult,N,N);
    im = fixImage(networkResult,N);
    imagesc(im,'Parent',handles.neurons);

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

function im = fixImage(im,N)
%    if isrgb(im)
	if length( size(im) ) == 3
        im = rgb2gray(im);
    end
    im = double(im);
    m = min(im(:));
    M = max(im(:));
    im = (im-m)/(M-m);  %normelizing the image
    im = imresize(im,[N N],'bilinear');
    %im = (im > 0.5)*2-1;    %changing image values to -1 & 1
    im = (im > 0.5);    %changing image values to 0 & 1

% --- Executes on slider movement.
function noiseAmount_Callback(hObject, eventdata, handles)
% hObject    handle to noiseAmount (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    percent = get(hObject,'value');
    percent = round(percent*100);
    set(handles.noisePercent,'string',num2str(percent));

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

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

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

author: Sable

repo: mcbench-benchmarks


text:
function geoEastOfNorth=magHeading2Geog(points,PEastOfNorth,coeffType,year,a,f)
%MAGHEADING2GEOG Convert a heading in terms of radians clockwise
%                (East) of MAGNETIC North to a heading in terms of radians
%                clockwise from GEOGRAPHIC North as defined on a particular
%                reference ellipsoid. The model for the Earth's magnetic
%                field can be selected. In the rare spots (i.e. near the
%                magnetic poles), where the magnetic field vector points
%                directly towards the reference ellipsoid, the geographic
%                heading will be undefined and a NaN will be returned.
%
%INPUTS: points One or more points given in geodetic latitude and
%               longitude, in radians, and height, in meters where the
%               magnetic headings apply. To convert N headings, points is a
%               3XN matrix with each column having the format
%               [latitude;longitude; height]. All points should be
%               associated with the same time. At the geographic poles, the
%               longitude value determines the orientation of the local
%               coordinate system axes. Thus, geographic headings ARE
%               defined at the poles.
%  PEastOfNorth An NX1 array of N magnetic headings in radians clockwise
%               from North that should be turned into geographic headings.
%     coeffType This specifies the coefficient model for the coefficients.
%               If one wishes to explcitely pass a model, then this is a
%               structure with members C, S, aH, and cH, which are defined
%               in the same manner as the return values from the function
%               getWMMCoeffs. Otherwise, a string can be passed. Possible
%               values are:
%               'WMM' (The default if omitted or an empty matrix is
%                     passed). Use the World Magnetic Model via the
%                     function getWMMCoeffs.
%               'IGRF' Use the International Geomagnetic Reference Field
%                     model via the function getIGRFCoeffs.
%          year This is only used if coeffType is a string. A decimal
%               number indicating a year in the Gregorian calendar as
%               specified by UTC. For example, halfway through the year
%               2012 would be represented as 2012.5. The precision of the
%               model is not sufficiently fine that leap seconds matter. If
%               this parameter is omitted, then the last reference epoch of
%               the geomagnetic model is used.
%             a The semi-major axis of the reference ellipsoid. If this
%               argument is omitted, the value in
%               Constants.WGS84SemiMajorAxis is used.
%             f The flattening factor of the reference ellipsoid. If this
%               argument is omitted, the value in
%               Constants.WGS84Flattening is used.
%
%OUTPUT: geoEastOfNorth The headings converted to radians clockwise of
%                       geographic North on the reference ellipse.
%
%First, spherical harmonic coefficients for the geomagnetic model are
%obtained at the desired time and points. Then, the magnetic flex vector B
%is determined at each of the points. Strictly speaking, the direction of
%magnetic North would be the component of the magnetic flux vector
%projected onto the local gravitationally horizontal plane at each point.
%However, since B is already rather imprecise, nothing is really lost by
%using a projection onto the local tangent plane of the reference
%ellipsoid of the Earth. Thus, rotating B about the negative of the local
%vertical to the reference ellipsoid by PEastOfNorth provides a vector
%whose projection into the local tangent plane points in the desired
%geographic heading.
%
%This function makes use of the functions getWMM2010Coeffs and
%spherHarmonicEval to determine the magnetic flux vector at the specified
%points as well as the function getENUAxes to determine the local
%East-North-Up coordinate axes.
%
%January 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

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

if(nargin<3||isempty(coeffType))
    coeffType='WMM';
end

if(nargin<4)
    year=[];
end

%Get the fully normalized spherical harmonic coefficients for the
%geomagnetic model.
if(isstruct(coeffType))
    %If the user explicitly provided the coefficients.
    C=coeffType.C;
    S=coeffType.S;
    aH=coeffType.aH;
    cH=coeffType.cH;
else
    switch(coeffType)
        case 'WMM'
            [C,S,aH,cH]=getWMMCoeffs(year);
        case 'IGRF'
            [C,S,aH,cH]=getIGRFCoeffs(year);
        otherwise
            error('Unknown coefficients selected')
    end
end

[~,gradV]=spherHarmonicEval(C,S,ellips2Sphere(points,a,f),aH,cH,true);
B=-gradV;
%B is now a matrix of vectors of the magnetic flux of the Earth's field at
%the points. The direction of B in the local tangent plane to the reference
%ellipsoid approximately defined magnetic north.

numPoints=size(points,2);
geoEastOfNorth=zeros(numPoints,1);%Allocate space
for curPoint=1:numPoints
    %-B(:,curPoint) defines the direction of magnetic North at the current
    %point. We want the direction vector PEastOfNorth(curPoint). To go East
    %of North, one must rotate that number of degrees about the local
    %"Down" vector.
    u=getENUAxes(points(:,curPoint),false,a,f);

    %The negative of u(:,3) is a unit vector in the down direction.
    vRot=rotateVectorAxisAng(B(:,curPoint),-u(:,3),PEastOfNorth(curPoint));

    %Given the rotated field vector, only the components in the local
    %tangent plane matter.
    vEast=dot(vRot,u(:,1));
    vNorth=dot(vRot,u(:,2));
    %Find the angle East of North.
    geoEastOfNorth(curPoint)=atan2(vEast,vNorth);
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ undx, xdnu ] = r8col_tol_undex ( m, n, a, unique_num, tol )

%*****************************************************************************80
%
%% R8COL_TOL_UNDEX indexes tolerably unique entries of an R8COL.
%
%  Discussion:
%
%    An R8COL is an M x N array of R8 values, regarded as N columns
%    each of M R8 values.
%
%    The goal of this routine is to determine a vector UNDX,
%    which points to the unique elements of A, in sorted order,
%    and a vector XDNU, which identifies, for each entry of A, the index of
%    the unique sorted element of A.
%
%    This is all done with index vectors, so that the elements of
%    A are never moved.
%
%    The first step of the algorithm requires the indexed sorting
%    of A, which creates arrays INDX and XDNI.  (If all the entries
%    of A are unique, then these arrays are the same as UNDX and XDNU.)
%
%    We then use INDX to examine the entries of A in sorted order,
%    noting the unique entries, creating the entries of XDNU and
%    UNDX as we go.
%
%    Once this process has been completed, the object A could be
%    replaced by a compressed object XU, containing the unique entries
%    of A in sorted order, using the formula
%
%      XU(1:UNIQUE_NUM) = A(UNDX(1:UNIQUE_NUM)).
%
%    We could then, if we wished, reconstruct the entire vector A, or
%    any element of it, by index, as follows:
%
%      A(I) = XU(XDNU(I)).
%
%    We could then replace A by the combination of XU and XDNU.
%
%    Later, when we need the I-th entry of A, we can locate it as
%    the XDNU(I)-th entry of XU.
%
%    Here is an example of a vector A, the sort and inverse sort
%    index vectors, and the unique sort and inverse unique sort vectors
%    and the compressed unique sorted vector.
%
%      I    A   Indx  Xdni      XU   Undx  Xdnu
%    ----+-----+-----+-----+--------+-----+-----+
%      1 : 11.     1     1 :    11.     1     1
%      2 : 22.     3     5 :    22.     2     2
%      3 : 11.     6     2 :    33.     4     1
%      4 : 33.     9     8 :    55.     5     3
%      5 : 55.     2     9 :                  4
%      6 : 11.     7     3 :                  1
%      7 : 22.     8     6 :                  2
%      8 : 22.     4     7 :                  2
%      9 : 11.     5     4 :                  1
%
%    INDX(2) = 3 means that sorted item(2) is A(3).
%    XDNI(2) = 5 means that A(2) is sorted item(5).
%
%    UNDX(3) = 4 means that unique sorted item(3) is at A(4).
%    XDNU(8) = 2 means that A(8) is at unique sorted item(2).
%
%    XU(XDNU(I))) = A(I).
%    XU(I)        = A(UNDX(I)).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the dimension of the data values.
%
%    Input, integer N, the number of data values.
%
%    Input, real A(M,N), the data values.
%
%    Input, integer UNIQUE_NUM, the number of unique values in A.
%    This value is only required for languages in which the size of
%    UNDX must be known in advance.
%
%    Input, real TOL, a tolerance for equality.
%
%    Output, integer UNDX(UNIQUE_NUM), the UNDX vector.
%
%    Output, integer XDNU(N), the XDNU vector.
%
  undx = zeros ( unique_num, 1 );
  xdnu = zeros ( n, 1 );
%
%  Implicitly sort the array.
%
  indx = r8col_sort_heap_index_a ( m, n, a );
%
%  Consider entry I = 1.
%  It is unique, so set the number of unique items to K.
%  Set the K-th unique item to I.
%  Set the representative of item I to the K-th unique item.
%
  i = 1;
  k = 1;
  undx(k) = indx(i);
  xdnu(indx(i)) = k;
%
%  Consider entry I.
%
%  If it is unique, increase the unique count K, set the
%  K-th unique item to I, and set the representative of I to K.
%
%  If it is not unique, set the representative of item I to a
%  previously determined unique item that is close to it.
%
  for i = 2 : n

    unique = 1;

    for j = 1 : k
      diff = max ( abs ( a(1:m,indx(i)) - a(1:m,undx(j)) ) );
      if ( diff <= tol )
        unique = 0;
        xdnu(indx(i)) = j;
        break
      end
    end

    if ( unique )
      k = k + 1;
      undx(k) = indx(i);
      xdnu(indx(i)) = k;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% PTrc2sp()

% from emuflight

function [RCcommand_dps] = PTrc2sp(RCcommand, rateSystem)
%function [RCcommand_dps] = PTrc2sp(RCcommand)
%   converts rc command to deg/s

BetaflightRates = 1;
BFActualRates = 2;

EmuflightRates = 2;

    rcCommandf    = RCcommand / 500.0;
    rcCommandfAbs = abs(rcCommandf);

          %  if rc_expo
                expof = rc_expo / 100;
                rcCommandf = rcCommandf * Math.pow(rcCommandfAbs, RC_EXPO_POWER) * expof + rcCommandf * (1-expof);
          %  end

            rcRate = rc_rates / 100.0;
          %  if (rcRate > 2.0)
                rcRate += RC_RATE_INCREMENTAL * (rcRate - 2.0);
          %  end

            angleRate = 200.0 * rcRate * rcCommandf;
         %   if rates
                rcSuperfactor = 1.0 / ((1.0 - (rcCommandfAbs * rates / 100.0)));
                angleRate = angleRate * rcSuperfactor;
         %   end

 return calculateSetpointRate(axis, value);

%%%%% From betaflight configurator RC3

 const minRc = 1000;
const midRc = 1500;
const maxRc = 2000;

const RateCurve = function (useLegacyCurve) {
    this.useLegacyCurve = useLegacyCurve;
    this.maxAngularVel = null;

    this.constrain = function (value, min, max) {
        return Math.max(min, Math.min(value, max));
    };

    this.rcCommand = function (rcData, rcRate, deadband) {
        const tmp = Math.min(Math.max(Math.abs(rcData - midRc) - deadband, 0), 500);

        let result = tmp * rcRate;

        if (rcData < midRc) {
            result = -result;
        }

        return result;
    };

    this.getBetaflightRates = function (rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo, superExpoActive, limit) {
        let angularVel;

        if (rcRate > 2) {
            rcRate = rcRate + (rcRate - 2) * 14.54;
        }

        let expoPower;
        let rcRateConstant;

        if (semver.gte(FC.CONFIG.apiVersion, "1.20.0")) {
            expoPower = 3;
            rcRateConstant = 200;
        } else {
            expoPower = 2;
            rcRateConstant = 205.85;
        }

        if (rcExpo > 0) {
            rcCommandf =  rcCommandf * Math.pow(rcCommandfAbs, expoPower) * rcExpo + rcCommandf * (1-rcExpo);
        }

        if (superExpoActive) {
            const rcFactor = 1 / this.constrain(1 - rcCommandfAbs * rate, 0.01, 1);
            angularVel = rcRateConstant * rcRate * rcCommandf; // 200 should be variable checked on version (older versions it's 205,9)
            angularVel = angularVel * rcFactor;
        } else {
            angularVel = (((rate * 100) + 27) * rcCommandf / 16) / 4.1; // Only applies to old versions ?
        }

        angularVel = this.constrain(angularVel, -1 * limit, limit); // Rate limit from profile

        return angularVel;
    };

    this.getRaceflightRates = function (rcCommandf, rate, rcRate, rcExpo) {
        let angularVel = ((1 + 0.01 * rcExpo * (rcCommandf * rcCommandf - 1.0)) * rcCommandf);
        angularVel = (angularVel * (rcRate + (Math.abs(angularVel) * rcRate * rate * 0.01)));
        return angularVel;
    };

    this.getKISSRates = function (rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo) {
        const kissRpy = 1 - rcCommandfAbs * rate;
        const kissTempCurve = rcCommandf * rcCommandf;
        rcCommandf = ((rcCommandf * kissTempCurve) * rcExpo + rcCommandf * (1 - rcExpo)) * (rcRate / 10);
        return ((2000.0 * (1.0 / kissRpy)) * rcCommandf);
    };

    this.getActualRates = function (rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo) {
        let angularVel;
        const expof = rcCommandfAbs * ((Math.pow(rcCommandf, 5) * rcExpo) + (rcCommandf * (1 - rcExpo)));

        angularVel = Math.max(0, rate-rcRate);
        angularVel = (rcCommandf * rcRate) + (angularVel * expof);

        return angularVel;
    };

    this.getQuickRates = function (rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo) {
        rcRate = rcRate * 200;
        rate = Math.max(rate, rcRate);

        let angularVel;
        const superExpoConfig = (((rate / rcRate) - 1) / (rate / rcRate));
        const curve = Math.pow(rcCommandfAbs, 3) * rcExpo + rcCommandfAbs * (1 - rcExpo);

        angularVel = 1.0 / (1.0 - (curve * superExpoConfig));
        angularVel = rcCommandf * rcRate * angularVel;

        return angularVel;
    };

};

RateCurve.prototype.rcCommandRawToDegreesPerSecond = function (rcData, rate, rcRate, rcExpo, superExpoActive, deadband, limit) {
    let angleRate;

    if (rate !== undefined && rcRate !== undefined && rcExpo !== undefined) {
        let rcCommandf = this.rcCommand(rcData, 1, deadband);
        if (semver.gte(FC.CONFIG.apiVersion, API_VERSION_1_43)) {
            rcCommandf = rcCommandf / (500 - deadband);
        } else {
            rcCommandf = rcCommandf / 500;
        }

        const rcCommandfAbs = Math.abs(rcCommandf);

        switch(TABS.pid_tuning.currentRatesType) {
            case TABS.pid_tuning.RATES_TYPE.RACEFLIGHT:
                angleRate=this.getRaceflightRates(rcCommandf, rate, rcRate, rcExpo);

                break;

            case TABS.pid_tuning.RATES_TYPE.KISS:
                angleRate=this.getKISSRates(rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo);

                break;

            case TABS.pid_tuning.RATES_TYPE.ACTUAL:
                angleRate=this.getActualRates(rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo);

                break;

            case TABS.pid_tuning.RATES_TYPE.QUICKRATES:
                angleRate=this.getQuickRates(rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo);

                break;

            // add future rates types here
            default: // BetaFlight
                angleRate=this.getBetaflightRates(rcCommandf, rcCommandfAbs, rate, rcRate, rcExpo, superExpoActive, limit);

                break;
        }
    }

    return angleRate;
};

RateCurve.prototype.getMaxAngularVel = function (rate, rcRate, rcExpo, superExpoActive, deadband, limit) {
    let maxAngularVel;
    if (!this.useLegacyCurve) {
        maxAngularVel = this.rcCommandRawToDegreesPerSecond(maxRc, rate, rcRate, rcExpo, superExpoActive, deadband, limit);
    }

    return maxAngularVel;
};

RateCurve.prototype.setMaxAngularVel = function (value) {
    this.maxAngularVel = Math.ceil(value/200) * 200;
    return this.maxAngularVel;

};

RateCurve.prototype.draw = function (rate, rcRate, rcExpo, superExpoActive, deadband, limit, maxAngularVel, context) {
    if (rate !== undefined && rcRate !== undefined && rcExpo !== undefined) {
        const height = context.canvas.height;
        const width = context.canvas.width;

        if (this.useLegacyCurve) {
            this.drawLegacyRateCurve(rate, rcRate, rcExpo, context, width, height);
        } else {
            this.drawRateCurve(rate, rcRate, rcExpo, superExpoActive, deadband, limit, maxAngularVel, context, width, height);
        }
    }
};

repo: PIDtoolbox


text:
function CPT = CPD_to_CPT(CPD)
% Compute the big CPT for an HHMM Q node (including F parents)
% by combining internal transprob and startprob
% function CPT = CPD_to_CPT(CPD)

Qsz = CPD.Qsz;

if ~isempty(CPD.Fbelow_ndx)
  if ~isempty(CPD.Fself_ndx) % general case
    error('not implemented')
  else % no F from self, hence no startprob (top level)
    nps = length(CPD.dom_sz)-1; % num parents
    CPT = 0*myones(CPD.dom_sz);
    % when Fself=1, the CPT(i,j) = delta(i,j) for all k
    for k=1:prod(CPD.Qpsizes)
      Qps_vals = ind2subv(CPD.Qpsizes, k);
      ndx = mk_multi_index(nps+1, [CPD.Fbelow_ndx CPD.Qps_ndx], [1 Qps_vals]);
      CPT(ndx{:}) = eye(Qsz); % CPT(:,2,k,:) or CPT(:,k,2,:) etc
    end
    ndx = mk_multi_index(nps+1, CPD.Fbelow_ndx, 2);
    CPT(ndx{:}) = CPD.transprob; % we assume transprob is in topo order
  end
else % no F signal from below
  if ~isempty(CPD.Fself_ndx) % bottom level
    nps = length(CPD.dom_sz)-1; % num parents
    CPT = 0*myones(CPD.dom_sz);
    ndx = mk_multi_index(nps+1, CPD.Fself_ndx, 1);
    CPT(ndx{:}) = CPD.transprob;
    ndx = mk_multi_index(nps+1, CPD.Fself_ndx, 2);
    CPT(ndx{:}) = CPD.startprob;
  else % no F from self
    error('An hhmmQ node without any F parents is just a tabular_CPD')
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function c = t_project_coefficients ( n, f )

%*****************************************************************************80
%
%% T_PROJECT_COEFFICIENTS: function projected onto T(0:n,x).
%
%  Discussion:
%
%    It is assumed that the interval of definition is -1 <= x <= +1.
%
%    Over this interval, f(x) will be well approximated by
%
%      f(x) approx sum ( 0 <= i <= n ) c(i) * T(i,x)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the highest order polynomial to compute.
%
%    Input, function handle F, of the form
%      function v = f ( x )
%
%    Output, real C(N+1,1), the projection coefficients of f(x) onto
%    T(0,x) through T(n,x).
%
  for k = 1 : n + 1
    y = cos ( pi * ( k - 0.5 ) / ( n + 1 ) );
    d(k) = f ( y );
  end

  fac = 2.0 / ( n + 1 );
  for j = 1 : ( n + 1 )
    sum = 0.0;
    for k = 1 : ( n + 1 )
      sum = sum + d(k) * cos ( ( pi * ( j - 1 ) ) * ( ( k - 0.5 ) / ( n + 1 ) ) );
    end
    c(j) = fac * sum;
  end

  c(1) = c(1) / 2.0;

  return
end


text:
classdef ShFunc_Chomog_fraction < ShFunc_Chomog

    properties (Access = private)
        invChAA
        invChBB
        invChAB
        invChBA
    end

    properties (Access = private)
        alpha
        beta
    end

    methods (Access = public)

        function obj=ShFunc_Chomog_fraction(cParams)
            obj.initChomog(cParams);
            obj.alpha = cParams.alpha/norm(cParams.alpha);
            obj.beta  = cParams.beta/norm(cParams.beta);
        end

        function computeFunctionValue(obj)
            obj.computeInvChProyections();
            obj.value = obj.invChAB/obj.invChAA + obj.invChBA/obj.invChBB;
        end

        function computeGradientValue(obj)
            obj.computeChDerivative();
            a = obj.alpha;
            b = obj.beta;
            beta1 = obj.invChAA*b - obj.invChAB*a;
            beta2 = obj.invChBB*a - obj.invChBA*b;
            g1    = obj.computedChInv(obj.Chomog,a,beta1);
            g2    = obj.computedChInv(obj.Chomog,b,beta2);
            grad = g1/(obj.invChAA)^2 + g2/(obj.invChBB)^2;
            obj.gradient = grad;
        end

    end

    methods (Access = private)

        function computeInvChProyections(obj)
            invCh = inv(obj.Chomog);
            a = obj.alpha;
            b = obj.beta;
            obj.invChAB  = obj.projectTensor(invCh,a,b);
            obj.invChBA  = obj.projectTensor(invCh,b,a);
            obj.invChAA  = obj.projectTensor(invCh,a,a);
            obj.invChBB  = obj.projectTensor(invCh,b,b);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [y,m] = getPathsAsian(S0, r, sigma, dt, NSim, Nr)
% path retrieval for Asian options

if length(S0) > 1
    lenS = length(S0);             % for using this with different starting values
    R = exp((r - sigma^2/2) * dt ...
    + sigma * sqrt(dt) * randn(NSim,Nr,lenS));   % random noise + drift
    tmp = zeros(NSim,1,lenS);
    tmp(:,1,:) = repmat(log(S0)',NSim,1);
    y = exp(cumsum([tmp, R], 2));    % path set simple example
else
    R = exp((r - sigma^2/2) * dt ...
    + sigma * sqrt(dt) * randn(NSim,Nr));
    y = cumprod([S0*ones(NSim,1), R], 2); %S0 in Spalte 1 -> Indexverschiebung

end

tmp = repmat(1:Nr, NSim, 1); %Hilfsmatrix
m = [y(:,1) cumsum(y(:,2:end), 2) ./ tmp];

end


text:
function [ Z, H, dnorm, lambdas ] = ssdeep_seminmf ( X, layers, y, misc, varargin )

% Process optional arguments
pnames = { ...
    'z0' 'h0' 'maxiter' 'TolFun', 'verbose', 'lambdas', 'cache', ...
    'geometric'
};

num_of_layers = numel(layers);

Z = cell(1, num_of_layers);
H = cell(1, num_of_layers);
L = cell(1, num_of_layers);
W = cell(1, num_of_layers);
Di = cell(1, num_of_layers);

dflts  = {0, 0, 500, 1e-5, 1, 1, 1};

[z0, h0, maxiter, tolfun, verbose, lambdas, cache] = ...
        internal.stats.parseArgs(pnames,dflts,varargin{:});

nSmp = size(X, 2);

if  ~iscell(h0)
    for i_layer = 1:length(layers)
        if i_layer == 1
            V = X;
        else
            V = H{i_layer-1};
        end

        if verbose
            display(sprintf('Initialising Layer #%d with k=%d with size(V)=%s...', i_layer, layers(i_layer), mat2str(size(V))));
        end

        if numel(y{i_layer}) == 1
            options = [];
            options.WeightMode = 'Binary';
            options.k = y{i_layer};
            W{i_layer} = constructW(V',options)';
        else

            options = [];
            options.NeighborMode = 'Supervised';
            options.gnd = y{i_layer};
            options.WeightMode = 'HeatKernel';
            options.t = 1;
            W{i_layer} = constructW(V', options);
%             W{i_layer} = zeros(nSmp, nSmp);
%
%             for j=1:nSmp;
%                 W{i_layer}(j, :) =  y{i_layer} == y{i_layer}(j);
%                 W{i_layer}(j, j) = 0;
%             end
        end

        W{i_layer} = lambdas(i_layer) * W{i_layer};
        DCol = full(sum(W{i_layer},2));
        Di{i_layer} = spdiags(DCol,0,nSmp,nSmp);

        if ~iscell(z0)
%              [Z{i_layer}, H{i_layer}, ~] = ...
%              seminmf(V, ...
%                  layers(i_layer), ...
%                  'maxiter', 500, 'verbose', verbose,  'save', cache, 'fast', 1);
        [Z{i_layer}, H{i_layer}, ~, ~, ~] = ...
             ssnmf(V, y{i_layer}, ...
                 layers(i_layer), ...
                 'maxiter', 500, 'verbose', verbose, 'lambda', lambdas(i_layer), 'save', cache);
        else
            display('Using existing Z');
          [Z{i_layer}, H{i_layer}, ~, ~, ~] = ...
             ssnmf(V, y{i_layer}, ...
                 layers(i_layer), ...
                 'maxiter', 500, 'verbose', verbose, 'lambda', lambdas(i_layer), 'save', cache, 'z0', z0);
        end

         L{i_layer} = Di{i_layer} - W{i_layer};
    end

else
    Z=z0;
    H=h0;

    if verbose
        display('Skipping initialization, using provided init matrices...');
    end
end

eval(Z, H, misc);

dnorm0 = cost_function(X, Z, H, L);
dnorm = dnorm0;

if verbose
    display(sprintf('#%d error: %f', 0, dnorm0));
end

eval(Z, H, misc);

%% Error Propagation
if verbose
    display('Finetuning...');
end

for iter = 1:maxiter
%     H_err{numel(layers)} = H{numel(layers)};
%     for i_layer = numel(layers)-1:-1:1
%         H_err{i_layer} = Z{i_layer+1} * H_err{i_layer+1};
%     end

    for i = 1:numel(layers)
        try
            if i == 1
                Z{i} = X  * pinv(Z{2} * H{2});
            else
                Z{i} = pinv(D') * X * pinv(H{i});
            end
        catch
            display(sprintf('Convergance error %f. min Z{i}: %f. max %f', norm(Z{i}, 'fro'), min(min(Z{i})), max(max(Z{i}))));
            any(isnan(Z{i}(:)))
            any(isinf(Z{i}(:)))
        end

        if i == 1
            D = Z{1}';
        else
            D = Z{i}' * D;
        end

        A = D * X;
        Ap = max(A, 0);
        An = -min(A, 0);

        B = D * D';

        Bp = max(B, 0);
        Bn = -min(B, 0);

        H{i} = H{i} .* ((Ap + Bn * H{i} + H{i} * W{i} ) ./ max(An + Bp * H{i} +  H{i} * Di{i}, 1e-6)).^0.5;
    end

    assert(i == numel(layers));

    if mod(iter, 50) == 0
        eval(Z, H, misc);
    end

    dnorm = cost_function(X, Z, H, L);

    if verbose
        if mod(iter, 10) == 0
        display(sprintf('#%d error: %f', iter, dnorm));
        end
    end

%     assert(dnorm <= dnorm0 + 1, ...
%         sprintf('Rec. error increasing! From %f to %f. (%d)', ...
%         dnorm0, dnorm, iter) ...
%     );

%     if dnorm0-dnorm <= tolfun*max(1,dnorm0)
%         if verbose
%             display( ...
%                 sprintf('Stopped at %d: dnorm: %f, dnorm0: %f', ...
%                     iter, dnorm, dnorm0 ...
%                 ) ...
%             );
%         end
%         break;
%     end

    dnorm0 = dnorm;
end
end

function eval(Z, H, misc)
    names = {'Pose', 'Emotion', 'Identity'};

    for i = 1:3;
        fprintf('%s: ', names{i});
        ....
            for j = 1:3;
                try
                    %          mdl = train(misc.Y_train{i}, sparse(reshape(cell2mat(Htrain), 1200, [])'), '-q');
                    mdl = train(misc.Y_train{i}, sparse(H{j}'), '-q');
                    D = Z{1};
                    for k = 2:j
                        D = D * Z{k};
                    end

                    Hr = pinv(D) * misc.Xr;
                    [~, ac, ~] = predict(misc.Y_test{i}, sparse(Hr'), mdl, '-q');
                    fprintf(1, '%.2f | ', ac(1));
                catch

                end
            end
            fprintf(1, '\n');
    end
end

function eval(Z, H, misc)
    names = {'Pose', 'Emotion', 'Identity'};

    fprintf('\n');
    for i = 1:3;
        fprintf('%s: ', names{i});
        for j = 1:length(H);
            %mdl = train(misc.Y_train{i}, sparse(reshape(cell2mat(Htrain), 1200, [])'), '-q');
            mdl = train(misc.Y_train{i}, sparse(H{j}'), '-q');
            D = Z{1};
            for k = 2:j
                D = D * Z{k};
            end

            Hr = pinv(D) * misc.Xr;
            [~, ac, ~] = predict(misc.Y_test{i}, sparse(Hr'), mdl, '-q');

            fprintf(1, '%.2f | ', ac(1));
        end
        fprintf(1, '\n');
    end
end

function error = cost_function(X, Z, H, L)
    error = norm(X - reconstruction(Z, H), 'fro');

    for i = 1:length(H);
       error = error + trace(H{i} * L{i} * H{i}');
    end
end

function [ out ] = reconstruction( Z, H )

    out = H{numel(H)};

    for k = numel(H) : -1 : 1;
        out =  Z{k} * out;
    end

end

author: jwyang


text:
function plotLinCon(A,b,Aeq,beq,data)
%PLOTLINCON Plot Linear Constraints on the current figure
%   plotLinCon(A,b,Aeq,beq)

%   Copyright (C) 2011 Jonathan Currie (IPL)

if(nargin < 5 || isempty(data))
    xl = xlim; yl = ylim;
else
    xl = data.xl; yl = data.yl;
end

hold on;

%Plot Inequality Constraints
if(~isempty(A))
    A = full(A);
    %Check for 1D, fake second column as zeros
    if(size(A,2) == 1)
        A = [A zeros(size(A,1),1)];
    end

    for i = 1:length(b)
        if(any(A(i,:)))
            if(all(A(i,:)) && b(i)) %normal, full A + b
                %%disp('normal')
                c = b(i) / A(i,2);
                m = -(b(i) / A(i,2)) / (b(i) / A(i,1));
                xi = (yl - c) / m;
            elseif(all(A(i,:)) && ~b(i)) %b is zero, line through origin
                %%disp('b is zero, full A')
                c = 0;
                m = -A(i,1)/A(i,2);
                xi = (yl - c) / m;
            elseif(b(i)) %bounds, b scales
                %%disp('partially empty A, b scales')
                m = inf;
            else %bounds, b is zero
                %%disp('partially empty A, empty b');
                m = inf;
            end

            %Normal Line
            if(~isinf(m) && ~isnan(m))
                %Grab a point (+1,+1) or (-1,+1)
                if(m <= 0)
                    x = [xi(1)+1; yl(1)+1];
                else
                    x = [xi(1)-1; yl(1)+1];
                end
                %Evaluate point to check direction of patch
                in = A(i,:)*x;
                %Patch based on type + gradient
                if(in > b(i)) % <=
                    %%disp('leq')
                    if(m <= 0)
                        %%disp('neg grad')
                        patch([xi(1) xl(2) xl(2) xi(2)],[yl(1) yl(1) yl(2) yl(2)],'y','FaceAlpha',0.3)
                    else
                        %%disp('pos grad')
                        patch([xi(1) xl(1) xl(1) xi(2)],[yl(1) yl(1) yl(2) yl(2)],'y','FaceAlpha',0.3)
                    end
                else % >=
                    %%disp('geq')
                    if(m <= 0)
                        %%disp('neg grad')
                        patch([xi(1) xl(1) xl(1) xi(2)],[yl(1) yl(1) yl(2) yl(2)],'y','FaceAlpha',0.3)
                    else
                        %%disp('pos grad')
                        patch([xi(1) xl(2) xl(2) xi(2)],[yl(1) yl(1) yl(2) yl(2)],'y','FaceAlpha',0.3)
                    end
                end
            %Bounds
            else
                if(A(i,1) == 0)
                    %%disp('A1 empty')
                    a = A(i,2);
%                     if(xl(1) == -1), xxl = 0.1; else xxl = 1; end
%                     in = a*(xl(1)+xxl) > b(i);
                    i1 = a*(yl(1))-b(i); i2 = a*(yl(2))-b(i);
                    if(i1 < i2), in = false; else in = true; end
                    vbnd = b(i)/a;
                    bndtype = 1; %horizontal (y)
                else
                    %%disp('A2 empty')
                    a = A(i,1);
%                     if(yl(1) == -1), yyl = 0.1; else yyl = 1; end
%                     in = a*(yl(1)+yyl) > b(i);
                    i1 = a*(xl(1))-b(i); i2 = a*(xl(2))-b(i);
                    if(i1 < i2), in = false; else in = true; end
                    vbnd = b(i)/a;
                    bndtype = 0; %vertical (x)
                end
                if(bndtype == 0) %vertical
                    %%disp('vertical bound')
                    if(in)
                        %%disp('lb')
                        patch([vbnd vbnd xl(1) xl(1)],[yl(1) yl(2) yl(2) yl(1)],'y','FaceAlpha',0.3)
                    else
                        %%disp('ub')
                        patch([vbnd vbnd xl(2) xl(2)],[yl(1) yl(2) yl(2) yl(1)],'y','FaceAlpha',0.3)
                    end
                else %horizontal
                    %%disp('horizontal bound')
                    if(in)
                        %%disp('lb')
                        patch([xl(1) xl(1) xl(2) xl(2)],[vbnd yl(1) yl(1) vbnd],'y','FaceAlpha',0.3)
                    else
                        %%disp('ub')
                        patch([xl(1) xl(1) xl(2) xl(2)],[vbnd yl(2) yl(2) vbnd],'y','FaceAlpha',0.3)
                    end
                end
            end
        end
    end
end

%Plot Equality Constraints
if(~isempty(Aeq))
    Aeq = full(Aeq);
    %Check for 1D, fake second column as zeros
    if(size(Aeq,2) == 1)
        Aeq = [Aeq zeros(size(Aeq,1),1)];
    end

    for i = 1:length(beq)
        if(any(Aeq(i,:)))
            if(all(Aeq(i,:)) && beq(i)) %normal, full A + b
                %disp('eq normal')
                c = beq(i) / Aeq(i,2);
                m = -(beq(i) / Aeq(i,2)) / (beq(i) / Aeq(i,1));
                xi = (yl - c) / m;
            elseif(all(Aeq(i,:)) && ~beq(i)) %b is zero, line through origin
                %disp('eq b is zero, full A')
                c = 0;
                m = -Aeq(i,1)/Aeq(i,2);
                xi = (yl - c) / m;
            elseif(beq(i)) %bounds, b scales
                %disp('eq partially empty A, b scales')
                m = inf;
            else %bounds, b is zero
                %disp('eq partially empty A, empty b');
                m = inf;
            end

            %Normal Line
            if(~isinf(m) && ~isnan(m))
                line(xi,yl);
            %Straight Line (H or V)
            else
                if(Aeq(i,1) == 0)
                    %disp('eq A1 empty')
                    a = Aeq(i,2);
                    vbnd = beq(i)/a;
                    line(xl,[vbnd vbnd]);
                else
                    %disp('eq A2 empty')
                    a = Aeq(i,1);
                    vbnd = beq(i)/a;
                    line([vbnd vbnd],yl);
                end

            end
        end
    end
end

hold off;

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function roi = roiCheckCoords(roi, mr, preserveCoords);
% Check that an ROI's coords match the coordinates
% of an MR object's data, applying transformations if
% needed.
%
% roi = roiCheckCoords(roi, mr, [preserveCoords]);
%
% roi: an ROI struct. See roiCreate.
%
% mr: an mr object, which may be the object on which the ROI is
%     defined, or may be a different object. The ROI can still be used
%     if a common coordinate space can be found (e.g., by scaling the
%     ROI's coords against the MR data, or applying a transformation
%     into a mutual space like Scanner space).
%
% preserveCoords: optional flag; if 1, will have the output roi
% have the same # of columns in roi.coords as the input roi.
% If 0, removes redundant columns. [Default 0]
% NOTE: importantly, setting preserveCoords to 1 risks that some
% functions may not work on the transformed ROI, since it will
% include coordinates which are out-of-range of the mr data.
% If preserveCoords==0, will prune these automatically.
%
%
% ras, 10/2005.
if nargin<2, help(mfilename); error('Not enough args.');    end
if notDefined('preserveCoords'), preserveCoords = 0;        end
if notDefined('roi'), error('Empty ROI.');                  end

mr = mrParse(mr);
roi = roiParse(roi, mr);

% recursively check several ROIs if requested
if length(roi) > 1
    for i = 1:length(roi)
        roi(i) = roiCheckCoords(roi(i), mr, preserveCoords);
    end
    return
end

if isempty(roi.coords)	% no need to check
	return
end

xform = [];

% There are a few ways to check if the ROI matches the MR object.
% First, check if the ROI is defined with respect to one of the
% existing spaces:
if ismember(lower(roi.reference), {'raw data in pixels' lower(mr.name)})
    xform = eye(4);     % defined in same space

elseif ischar(roi.reference)
    % see if the roi and mr share a common space
    % (but ignore "standard" spaces: may have the same name but refer to
    % different coordinates)
    std = {'Raw Data in Pixels' 'Raw Data in mm' 'L/R Flipped'};
    commonSpaces = setdiff({mr.spaces.name}, std);
    if ~isempty( cellfind(commonSpaces, roi.reference) )
        fprintf('Coregistering via common space %s \n', roi.reference);
        I = cellfind({mr.spaces.name}, roi.reference);
        xform = inv(mr.spaces(I).xform);
    end
end

if isempty(xform)
    % still haven't gotten xform: try the next step:
    if isfield(roi,'referenceMR') & ~isempty(roi.referenceMR)
        % ensure the reference MR object is loaded
        roi.referenceMR = mrParse(roi.referenceMR);

        % compare the mr name and the referenceMR name
        if isequal(roi.referenceMR.name, mr.name)
            % great, the ROI was defined on this data! We're done.
            return

        else
            % check the spaces in the reference MR and the mr file
            try
                xform = mrBaseXform(mr, roi.referenceMR);
            catch
                disp(lasterr);
                error('Sorry, can''t use this ROI with that MR data.')
            end

        end

    else
        error('Sorry, can''t use this ROI with that MR data.');
    end
end

% if we got here, we have an xform -- apply it to the
% current ROI coords to get it into the mr data space:
if ~isequal(xform, eye(4)) % don't waste time on a trivial xform
	roi = roiXformCoords(roi, xform, mr.voxelSize(1:3));
end
roi.voxelSize = mr.voxelSize;
roi.dimUnits = mr.dimUnits;

% remove redundant coords, if desired
if ~preserveCoords
    roi.coords = intersectCols(roi.coords, roi.coords);

    % also remove coords out of range of mr data
    ok = find(roi.coords(1,:)>=1 & roi.coords(1,:)<=mr.dims(1) & ...
              roi.coords(2,:)>=1 & roi.coords(2,:)<=mr.dims(2) & ...
              roi.coords(3,:)>=1 & roi.coords(3,:)<=mr.dims(3));
    roi.coords = roi.coords(:,ok);
end

roi.reference = mr.name;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function diftab = dif_basis ( ntab, xtab )

%*****************************************************************************80
%
%% DIF_BASIS computes all Lagrange basis polynomials in divided difference form.
%
%  Discussion:
%
%    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
%    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
%    XTAB(J) for J not equal to I, and 1 when J is equal to I.
%
%    The Lagrange basis polynomials have the property that the interpolating
%    polynomial through a set of NTAB data points (XTAB,YTAB) may be
%    represented as
%
%      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
%
%    Higher order interpolation at selected points may be accomplished
%    using repeated X values, and scaled derivative values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer NTAB, the number of X data points XTAB, and the number of
%    basis polynomials to compute.
%
%    Input, real XTAB(NTAB), the X values upon which the
%    Lagrange basis polynomials are to be based.
%
%    Output, real DIFTAB(NTAB,NTAB), the set of divided
%    difference tables.  Column I of DIFTAB contains the table for
%    the I-th Lagrange basis polynomial.
%

%
%  Initialize DIFTAB to the identity matrix.
%
  diftab(1:ntab,1:ntab) = 0.0;
  for i = 1 : ntab
    diftab(i,i) = 1.0;
  end
%
%  Compute each Lagrange basis polynomial.
%
  for i = 1 : ntab
    diftab(1:ntab,i) = ( data_to_dif ( ntab, xtab, diftab(1:ntab,i) ) )';
  end

  return
end


text:
function circle_forward ( a )

%*****************************************************************************80
%
%% CIRCLE_FORWARD plots the image of the unit circle under the 2x2 matrix A.
%
%  Discussion:
%
%    A typical 2x2 matrix will map the unit circle to some kind of tilted ellipse.
%    The aspect ratio of the ellipse (ratio of major to minor axes) is a measure
%    of the conditioning of the matrix.
%
%    A singular matrix maps the unit circle to a line.
%
%    A diagonal matrix maps the unit circle to an ellipse with no tilting.
%    (There is no rotation.)
%
%    An orthogonal matrix maps the unit circle to the unit circle (but points
%    may have rotated.)
%
%    The identity matrix maps the unit circle to itself.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(2,2), the matrix whose mapping of the unit circle
%    is to be studied.
%

%
%  This call produces matrices U, S and V.
%
% [ u, s, v ] = svd ( a );
%
%  This call produces just the singular values, and in a vector, not a matrix.
%
  s = svd ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Singular values of A are %f and %f\n', s(1), s(2) );
%
%  Select N+1 evenly spaced points X on the unit circle.
%
  n = 32;

  i = [ 0 : n ];

  angle = 2 * pi * i / n;

  r = 1;

  x = [ r * cos(angle); r * sin(angle) ];
%
%  AX contains the image of each point under the mapping X -> A*X.
%
  ax =  a * x;
%
%  Plot the corners of the region.  This is a trick to force MATLAB
%  to plot in a square with equal axes.  I don't know why it seems to
%  hard to get this to happen in general!
%
  s = [];

  p_min = min ( min ( min ( ax ) ), min ( min ( x ) ) );
  p_max = max ( max ( max ( ax ) ), max ( max ( x ) ) );

  scatter ( [ p_min, p_max, p_max, p_min ], [ p_min, p_min, p_max, p_max ], s, 'k' );
  axis ( [ p_min, p_max, p_min, p_max ] );

  axis equal
  axis tight
  axis square

  hold on
%
%  Plot the points X on the unit circle.
%
  scatter ( x(1,:), x(2,:), s, 'b', 'filled' );
%
%  Plot the points AX on the image of the unit circle.
%
  scatter ( ax(1,:), ax(2,:), s, 'r', 'filled' );
%
%  Plot vectors for "3-oclock" and "12-oclock".
%
  line ( [ 0, x(1,1)], [ 0, x(2,1)], 'color', 'b' );
  line ( [ 0, x(1,n/4+1)], [ 0, x(2,n/4+1)], 'color', 'b' );
%
%  Plot the images of "3-oclock" and "12-oclock".
%
  line ( [ 0, ax(1,1)], [ 0, ax(2,1)], 'color', 'r' );
  line ( [ 0, ax(1,n/4+1)], [ 0, ax(2,n/4+1)], 'color', 'r' );

  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
  title ( 'X in blue, A*X in red' );

  hold off

  return
end


text:
function varargout = compileAndRun(fname,varargin)
% Compile a Matlab function via emlmex and then run the compiled version.
% The function must reside in its own file and be eml compliant. See the
% help sections on emlmex for details.
%
% fname     - the string name of the function
% varargin  - arguments to the compiled function
%
% Example:
%
% [a b c] = compileAndRun('myfunction',X,Y,Z,3); % where X,Y,Z,3 are inputs
%           to 'myfunction', which returns three outputs.
%
%
% Matthew Dunham
% www.cs.ubc.ca/~mdunham

    if(~exist('emlmex','file'))
        error('emlmex could not be found');
    end

    clear mex;
    if(~exist(fname,'file'))
       error('%s could not be found',fname);
    end

    prevDirectory = pwd;
    cd(fileparts(which(fname)));
    try
        fprintf('\ncompiling %s...\n',fname);
        tic;
        emlmex(fname,'-eg',varargin);
        t = toc;
        fprintf('\ncompiled successfully in %g seconds\n',t);
    catch
        warning('compileAndRun:mex','\n\ncould not compile, performance might be slow.\n\n');
    end
    fprintf('\nexecuting function...\n');
    tic
    nout = nargout([fname,'.m']);
    varargout = {};
    if(nout == 0)
        feval(fname,varargin{:});
    elseif(nout < 0)
        varargout{1} = feval(fname,varargin{:});
    else
        out = '[';
        for i=1:nout
            out = [out,sprintf(' varargout{%d} ',i)]; %#ok
        end
        out = [out,']'];
        eval([out,' = ','feval(fname,varargin{:})']);
    end

    clear mex;
    delete([fname,'.',mexext]);
    cd(prevDirectory);
    t = toc;
    fprintf('\n%s ran in %g seconds.\n',fname,t);
end


text:
function y = givensrot(gjc,g,x,K) %#ok
% y = givensrot(gjc,g,x,K)
%
% GIVENSROT
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Stores the result of a regression algorithm
classdef RegressionResult < Data
    properties (Access = public)
        predictedResults;
        truthResults;
        timestamps;
    end

    methods (Access = public)
        function obj = RegressionResult(predictedResults,truthResults,timestamps)
            if nargin > 0
                obj.predictedResults = predictedResults;
                if nargin > 1
                    obj.truthResults = truthResults;
                    if nargin > 2
                        obj.timestamps = timestamps;
                    end
                end
                if isempty(obj.timestamps)
                    obj.timestamps = 1 : length(predictedResults);
                end
            end
            obj.type = DataType.kRegressionResult;
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [X,freq] = absoluteValFFT(x,Fs)
N = length(x); %get the number of points
k = 0:N-1;     %create a vector from 0 to N-1
T = N/Fs;      %get the frequency interval
freq = k/T;    %create the frequency range
X = fft(x)/N*2; % normalize the data

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2);

%take only the first half of the spectrum
X = X(1:cutOff,:,:);
X = abs(X/N);
freq = freq(1:cutOff);

author: PatternRecognition

repo: OpenBMI


text:
function tests = test_ft_connectivity_pdc

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivity_pdc

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  func = localfunctions;
  for i=1:numel(func)
    fprintf('evaluating %s\n', func2str(func{i}));
    feval(func{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_rpt_chan_chan_freq(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nrpt   = 10;
nchan  = 3;
nfreq  = 4;
H      = randn(nrpt, nchan, nchan, nfreq);

dat = randn(nchan,1000);
noisecov = cov(dat');

result = {};
result{end+1} = ft_connectivity_pdc(H, 'noisecov', []);
result{end+1} = ft_connectivity_pdc(H, 'noisecov', noisecov);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_rpt_chan_chan_freq_time(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nrpt   = 10;
nchan  = 3;
nfreq  = 4;
ntime  = 5;
H      = randn(nrpt, nchan, nchan, nfreq, ntime);

dat = randn(nchan,1000);
noisecov = cov(dat');

result = {};
result{end+1} = ft_connectivity_pdc(H, 'noisecov', []);
result{end+1} = ft_connectivity_pdc(H, 'noisecov', noisecov);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% demo script of getting tightly cropped faces on IJB-B test set.

% data_dir is the root directory of your IJBB data
% code_dir is the root directory of your detector code, e.g. '/mtcnn_face',
% including standard caffe, the scripts and dependent libraries supporting MTCNN
% (configuration details can be found in https://github.com/kpzhang93/MTCNN_face_detection_alignment)
% e.g. /mtcnn_face/caffe, /mtcnn_face/mtcnn, /mtcnn_face/pdollar_toolbox.
% det_tight_data_dir is the root directory for saving crops. Please create the directory "$data_dir/debug"
% in case of converting image format.

% The coordinates of target facial regions are provided by IJB-B.
% For convenience, we concatenate the files including meta information of faces
% (e.g. ijbb_1N_gallery_S1.csv, ijbb_1N_gallery_S2.csv) and adding the term "FACEID" for face indexing ([1:face_num]).

data_dir = YOUR_IJBB_ROOT_DIR;
code_dir = YOUR_CODE_DIR;
table_path = DATA_CSV_PATH;
det_tight_data_dir = YOUR_TARGET_DIR;

database = readtable(table_path);
imglist = database.FILENAME;
x = database.FACE_X;
y = database.FACE_Y;
width = database.FACE_WIDTH;
height = database.FACE_HEIGHT;
faceid = database.FACEID; %index the target faces
%path of toolbox
pdollar_toolbox_path = fullfile(code_dir, 'pdollar_toolbox';
code_path = fullfile(code_dir, 'mtcnn/code/codes/MTCNNv2');%code/codes/MTCNNv2 is the path strcture used in MTCNN's github code.
addpath(genpath(pdollar_toolbox_path));
addpath(code_path);
caffe_model_path = fullfile(code_path, '/model');
%use cpu
%caffe.set_mode_cpu();
gpu_id = 0;
caffe.set_mode_gpu();
caffe.set_device(gpu_id);

%three steps's threshold
threshold=[0.6 0.7 0.7];
%scale factor
factor=0.709;
det_tight_ratio = 0.3;
%load caffe models
prototxt_dir =fullfile(caffe_model_path,'/det1.prototxt');
model_dir = fullfile(caffe_model_path,'/det1.caffemodel');
PNet=caffe.Net(prototxt_dir, model_dir,'test');
prototxt_dir = fullfile(caffe_model_path,'/det2.prototxt');
model_dir = fullfile(caffe_model_path,'/det2.caffemodel');
RNet=caffe.Net(prototxt_dir, model_dir,'test');
prototxt_dir = fullfile(caffe_model_path,'/det3.prototxt');
model_dir = fullfile(caffe_model_path,'/det3.caffemodel');
ONet=caffe.Net(prototxt_dir, model_dir,'test');
prototxt_dir =  fullfile(caffe_model_path,'/det4.prototxt');
model_dir =  fullfile(caffe_model_path,'/det4.caffemodel');
LNet=caffe.Net(prototxt_dir, model_dir,'test');

% images detector failed on
non_match_gt = zeros(length(imglist),1);
for i = 1: length(imglist)
    i
    img_name = imglist{i};
    img_path = fullfile(data_dir, img_name);
    try
 	img = imread(img_path);
    catch
        f1 = fullfile(data_dir, ['debug/', img_name]);
        f2 = fullfile(data_dir, ['debug/1_', img_name]);
        copyfile(img_path, f1);
        system(sprintf('convert -colorspace RGB %s %s', f1, f2));
        img = imread(f2);
    end
    gt_box = [x(i), y(i), x(i) + width(i), y(i) + height(i)];
    % crop a loose candidate region
    center_gt = [x(i) + 0.5 * width(i), y(i) + 0.5 * height(i)];
    extent_gt = 1.4;
    candidate_region = img(max(1, round(center_gt(2) - height(i)* extent_gt /2)): min(size(img,1), ...
        round(center_gt(2) + height(i) * extent_gt /2)), max(1, round(center_gt(1) - width(i) * extent_gt /2)): ...
        min(size(img,2), round(center_gt(1) + width(i) * extent_gt /2)),:);
    img = candidate_region;
    %we recommend you to set minsize as x * short side
    if size(img,3) == 1
        t_img = repmat(img, 1,1,3);
        img = t_img;
    end
    minl = min([size(img,1) size(img,2)]);
    minsize = fix(minl*0.1);
    tic
    [boundingboxes points]=detect_face(img,minsize,PNet,RNet,ONet,LNet,threshold,false,factor);
    toc
    numbox=size(boundingboxes,1);
    scal_ratio_min = 112/128;
    if numbox == 0
        non_match_gt(i) = 1;
        fprintf('%d image cannot match gt\n',i);
    else
      % select the bb close to the center
      bb_center = [(boundingboxes(:,1) + boundingboxes(:,3))/2, (boundingboxes(:,2) + boundingboxes(:,4))/2];
      bb_gt_dist = abs(bb_center - repmat([size(img,2)/2, size(img,1)/2], size(bb_center,1),1));
      [min_dist, region_id]= min(sum(bb_gt_dist, 2));
      org_coord = boundingboxes(region_id, 1:2);
      center_point = [(boundingboxes(region_id,1) + boundingboxes(region_id, 3))/2, ...,
              (boundingboxes(region_id,2) + boundingboxes(region_id,4))/2];
      bb_size = boundingboxes(region_id, 3:4) - org_coord;
      scal_ratio = max(scal_ratio_min, bb_size(1)/bb_size(2));
      tight_bb_w_h = bb_size(2) * [scal_ratio, 1] * (1 + det_tight_ratio);
      org_tight_region = img(max(1, round(center_point(2) - tight_bb_w_h(2)/2)): ...
              min(size(img,1), round(center_point(2) + tight_bb_w_h(2)/2)), ...
              max(1, round(center_point(1) - tight_bb_w_h(1)/2)): ...
              min(size(img,2), round(center_point(1) + tight_bb_w_h(1)/2)),:);
      imwrite(org_tight_region, fullfile(det_tight_data_dir, [num2str(faceid(i)), '.jpg']));
    end
end

% dealing with empty files
extend_ratio = 0.1;
non_ind = find(non_match_gt);
for i = 1:numel(non_ind)
    i
    img_name = imglist{non_ind(i)};
    img_path = fullfile(data_dir, img_name);
    img = imread(img_path);
    %cropping with GT_0.1
    bb_width = width(non_ind(i)) * (1 + extend_ratio);
    bb_height = height(non_ind(i)) * (1 + extend_ratio);
    xy = [x(non_ind(i)), y(non_ind(i))];
    centrepoint = [round(xy(1) + bb_width/2) round(xy(2) + bb_height/2)];
    x1 = centrepoint(1) - round(bb_width / 2);
    y1 = centrepoint(2) - round(bb_height / 2);
    x2 = centrepoint(1) + round(bb_width / 2);
    y2 = centrepoint(2) + round(bb_height / 2);

    x1 = max(1, x1);
    y1 = max(1, y1);
    x2 = min(x2, size(img,2));
    y2 = min(y2, size(img,1));

    crop = img(y1:y2, x1:x2, :);
    imwrite(crop, fullfile(det_tight_data_dir, [num2str(faceid(non_ind(i))), '.jpg']));
end

author: ox-vgg


text:
%Corrects the navigation states using the Kalman generated error values
%if size(att)=[4,1], it assumes att=qbn else it assumes att=Cbn
%vel=Vn
%dx=corrections where position errors are defined as dr^n (order:[pos, vel, att])
%delta_x(k)=K(delta_y(k)-h(delta_x(k|k-1)))=K*delta_y(k)

function [att_new, Ve_new, ecef_new]=correctnav_eframe_v000(att, Ve, ecef, dx)

%Correct the position
ecef_new=ecef-dx(1:3);

%velocity correction
Ve_new=Ve-dx(4:6);

%attitude correction
if (size(att,1)==3) %attitude is represented as DCM
    Cet=rot2dcm_v000(dx(7:9));   %Erroneous to true navigation frame:DCM=(I+S(ang))
    att_new=Cet*att;
else %attitude is quaternion
    %Method I
    qet=rvec2quat_v000(dx(7:9));
    att_new=quatmult_v000(qet, att);

%     %Method II:
%     %Method I involves some trigonometric functions. Here is a
%     %algebraic correction method (See: Chung 1996-Eq:16) (Special thanks to Kaygisiz)
%     R=[-q(2) -q(3) -q(4);q(1) q(4) -q(3);q(-4) q(1) q(2);q(3) q(-2) q(1)];
%     dq=-0.5*R*dx(7:9);
%     att_new=att_new-dq;
%     att_new=att_new/norm(att_new);  %optional
end

repo: ekfmonoslam


text:
function ROCout=roc(varargin)
% ROC - Receiver Operating Characteristics.
% The ROC graphs are a useful tecnique for organizing classifiers and
% visualizing their performance. ROC graphs are commonly used in medical
% decision making.
%
% Syntax: ROCout=roc(x,thresholds,alpha,verbose)
%
% Input: x - This is a Nx2 data matrix. The first column is the column of the data value;
%            The second column is the column of the tag: unhealthy (1) and
%            healthy (0).
%        Thresholds - If you want to use all unique values in x(:,1)
%            then set this variable to 0 or leave it empty;
%            else set how many unique values you want to use (min=3);
%        alpha - significance level (default 0.05)
%        verbose - if you want to see all reports and plots (0-no; 1-yes by
%        default);
%
% Output: if verbose = 1
%         the ROCplots, the sensitivity and specificity at thresholds; the Area
%         under the curve with Standard error and Confidence interval and
%         comment.
%         if ROCout is declared, you will have a struct:
%         ROCout.AUC=Area under the curve (AUC);
%         ROCout.SE=Standard error of the area;
%         ROCout.ci=Confidence interval of the AUC
%         ROCout.co=Cut off points
%         ROCdata.xr and ROCdata.yr points for ROC plot
%
% USING roc WITHOUT ANY DATA, IT WILL RUN A DEMO
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2008) ROC curve: compute a Receiver Operating Characteristics curve.
% http://www.mathworks.com/matlabcentral/fileexchange/19950

%Input Error handling
args=cell(varargin);
nu=numel(args);
if isempty(nu)
    error('Warning: almost the data matrix is required')
elseif nu>4
    error('Warning: Max four input data are required')
end
default.values = {[165 1;140 1;154 1;139 1;134 1;154 1;120 1;133 1;150 1;...
146 1;140 1;114 1;128 1;131 1;116 1;128 1;122 1;129 1;145 1;117 1;140 1;...
149 1;116 1;147 1;125 1;149 1;129 1;157 1;144 1;123 1;107 1;129 1;152 1;...
164 1;134 1;120 1;148 1;151 1;149 1;138 1;159 1;169 1;137 1;151 1;141 1;...
145 1;135 1;135 1;153 1;125 1;159 1;148 1;142 1;130 1;111 1;140 1;136 1;...
142 1;139 1;137 1;187 1;154 1;151 1;149 1;148 1;157 1;159 1;143 1;124 1;...
141 1;114 1;136 1;110 1;129 1;145 1;132 1;125 1;149 1;146 1;138 1;151 1;...
147 1;154 1;147 1;158 1;156 1;156 1;128 1;151 1;138 1;193 1;131 1;127 1;...
129 1;120 1;159 1;147 1;159 1;156 1;143 1;149 1;160 1;126 1;136 1;150 1;...
136 1;151 1;140 1;145 1;140 1;134 1;140 1;138 1;144 1;140 1;140 1;159 0;...
136 0;149 0;156 0;191 0;169 0;194 0;182 0;163 0;152 0;145 0;176 0;122 0;...
141 0;172 0;162 0;165 0;184 0;239 0;178 0;178 0;164 0;185 0;154 0;164 0;...
140 0;207 0;214 0;165 0;183 0;218 0;142 0;161 0;168 0;181 0;162 0;166 0;...
150 0;205 0;163 0;166 0;176 0;],0,0.05,1};
default.values(1:nu) = args;
[x threshold alpha verbose] = deal(default.values{:});
if isvector(x)
    error('Warning: X must be a matrix')
end
if ~all(isfinite(x(:))) || ~all(isnumeric(x(:)))
    error('Warning: all X values must be numeric and finite')
end
x(:,2)=logical(x(:,2));
if all(x(:,2)==0)
    error('Warning: there are only healthy subjects!')
end
if all(x(:,2)==1)
    error('Warning: there are only unhealthy subjects!')
end
if nu>=2
    if isempty(threshold)
        threshold=0;
    else
        if ~isscalar(threshold) || ~isnumeric(threshold) || ~isfinite(threshold)
            error('Warning: it is required a numeric, finite and scalar THRESHOLD value.');
        end
        if threshold ~= 0 && threshold <3
            error('Warning: Threshold must be 0 if you want to use all unique points or >=2.')
        end
    end
    if nu>=3
        if isempty(alpha)
            alpha=0.05;
        else3
            if ~isscalar(alpha) || ~isnumeric(alpha) || ~isfinite(alpha)
                error('Warning: it is required a numeric, finite and scalar ALPHA value.');
            end
            if alpha <= 0 || alpha >= 1 %check if alpha is between 0 and 1
                error('Warning: ALPHA must be comprised between 0 and 1.')
            end
        end
    end
    if nu==4
        verbose=logical(verbose);
    end
end
clear args default nu

tr=repmat('-',1,100);
lu=length(x(x(:,2)==1)); %number of unhealthy subjects
lh=length(x(x(:,2)==0)); %number of healthy subjects
z=sortrows(x,1);
if threshold==0
    labels=unique(z(:,1));%find unique values in z
else
    K=linspace(0,1,threshold+1); K(1)=[];
    labels=quantile(unique(z(:,1)),K)';
end
labels(end+1)=labels(end)+1;
ll=length(labels); %count unique value
a=zeros(ll,2); b=a; c=zeros(ll,1);%array preallocation
ubar=mean(x(x(:,2)==1),1); %unhealthy mean value
hbar=mean(x(x(:,2)==0),1); %healthy mean value
for K=1:ll
    if hbar<ubar
        TP=length(x(x(:,2)==1 & x(:,1)>labels(K)));
        FP=length(x(x(:,2)==0 & x(:,1)>labels(K)));
        FN=length(x(x(:,2)==1 & x(:,1)<=labels(K)));
        TN=length(x(x(:,2)==0 & x(:,1)<=labels(K)));
    else
        TP=length(x(x(:,2)==1 & x(:,1)<labels(K)));
        FP=length(x(x(:,2)==0 & x(:,1)<labels(K)));
        FN=length(x(x(:,2)==1 & x(:,1)>=labels(K)));
        TN=length(x(x(:,2)==0 & x(:,1)>=labels(K)));
    end
    M=[TP FP;FN TN];
    a(K,:)=diag(M)'./sum(M); %Sensitivity and Specificity
    b(K,:)=[a(K,1)/(1-a(K,2)) (1-a(K,1))/a(K,2)]; %Positive and Negative likelihood ratio
    c(K)=trace(M)/sum(M(:)); %Efficiency
end
b(isnan(b))=Inf;

if hbar<ubar
    xroc=flipud(1-a(:,2)); yroc=flipud(a(:,1)); %ROC points
    labels=flipud(labels);
else
    xroc=1-a(:,2); yroc=a(:,1); %ROC points
end

st=[1 mean(xroc) 1]; L=[0 0 0]; U=[Inf 1 Inf];
fo_ = fitoptions('method','NonlinearLeastSquares','Lower',L,'Upper',U,'Startpoint',st);
ft_ = fittype('1-1/((1+(x/C)^B)^E)',...
     'dependent',{'y'},'independent',{'x'},...
     'coefficients',{'B', 'C', 'E'});
cfit = fit(xroc,yroc,ft_,fo_);

xfit=linspace(0,1,500);
yfit=feval(cfit,xfit);
Area=trapz(xfit,yfit); %estimate the area under the curve
%standard error of area
Area2=Area^2; Q1=Area/(2-Area); Q2=2*Area2/(1+Area);
V=(Area*(1-Area)+(lu-1)*(Q1-Area2)+(lh-1)*(Q2-Area2))/(lu*lh);
Serror=realsqrt(V);
%confidence interval
ci=Area+[-1 1].*(realsqrt(2)*erfcinv(alpha)*Serror);
if ci(1)<0; ci(1)=0; end
if ci(2)>1; ci(2)=1; end
m=zeros(1,4);
%z-test
SAUC=(Area-0.5)/Serror; %standardized area
p=1-0.5*erfc(-SAUC/realsqrt(2)); %p-value

if verbose
    %Performance of the classifier
    if Area==1
        str='Perfect test';
    elseif Area>=0.90 && Area<1
        str='Excellent test';
    elseif Area>=0.80 && Area<0.90
        str='Good test';
    elseif Area>=0.70 && Area<0.80
        str='Fair test';
    elseif Area>=0.60 && Area<0.70
        str='Poor test';
    elseif Area>=0.50 && Area<0.60
        str='Fail test';
    else
        str='Failed test - less than chance';
    end
    %display results
    disp('ROC CURVE ANALYSIS')
    disp(' ')
    disp(tr)
    str2=['AUC\t\t\tS.E.\t\t\t\t' num2str((1-alpha)*100) '%% C.I.\t\t\tComment\n'];
    fprintf(str2)
    disp(tr)
    fprintf('%0.5f\t\t\t%0.5f\t\t\t%0.5f\t\t%0.5f\t\t\t%s\n',Area,Serror,ci,str)
    disp(tr)
    fprintf('Standardized AUC\t\t1-tail p-value\n')
    if p<1e-4
        fprintf('%0.4f\t\t\t\t%0.4e',SAUC,p)
    else
        fprintf('%0.4f\t\t\t\t%0.4f',SAUC,p)
    end
    if p<=alpha
        fprintf('\t\tThe area is statistically greater than 0.5\n')
    else
        fprintf('\t\tThe area is not statistically greater than 0.5\n')
    end
    disp(' ')
    %display graph
    H=figure;
    set(H,'Position',[4 402 560 420])
    hold on
    plot([0 1],[0 1],'k');
    plot(xfit,yfit,'marker','none','linestyle','-','color','r','linewidth',2);
    H1=plot(xroc,yroc,'bo');
    set(H1,'markersize',6,'markeredgecolor','b','markerfacecolor','b')
    hold off
    xlabel('False positive rate (1-Specificity)')
    ylabel('True positive rate (Sensitivity)')
    title(sprintf('ROC curve (AUC=%0.4f)',Area))
    axis square
end

if p<=alpha
    ButtonName = questdlg('Do you want to input the true prevalence?', 'Prevalence Question', 'Yes', 'No', 'Yes');
    if strcmp(ButtonName,'Yes')
        ButtonName = questdlg('Do you want to input the true prevalence as:', 'Prevalence Question', 'Ratio', 'Probability', 'Ratio');
        switch ButtonName
            case 'Ratio'
                prompt={'Enter the Numerator or the prevalence ratio:','Enter the denominator or the prevalence ratio:'};
                name='Input for Ratio prevalence';
                Ratio=str2double(inputdlg(prompt,name));
                POD=Ratio(1)/diff(Ratio); %prior odds
            case 'Probability'
                prompt={'Enter the prevalence probability comprised between 0 and 1:'};
                name='Input for prevalence';
                pr=str2double(inputdlg(prompt,name));
                POD=pr/(1-pr); %prior odds
        end
        d=[1./(1+1./(b(:,1).*POD)) 1./(1+(b(:,2).*POD))];
        d((a(:,1)==0 & a(:,2)==1),1)=NaN;
        d((a(:,1)==1 & a(:,2)==0),2)=NaN;
        table=[labels'; a(:,1)'; a(:,2)';c';d(:,1)'.*100; d(:,2)'.*100;]';
        if verbose
            disp('ROC CURVE DATA')
            disp(tr)
            fprintf('Cut-off \tSensitivity\tSpecificity\tEfficiency\tPos.Pred.\tNeg.Pred.\n')
            fprintf('%0.2f\t\t%0.4f\t\t%0.4f\t\t%0.4f\t\t%0.2f\t\t%0.2f\n',table')
            disp(tr)
            disp(' ')
        end
    else
        table=[labels'; a(:,1)'; a(:,2)';c']';
        if verbose
            disp('ROC CURVE DATA')
            disp(tr)
            fprintf('Cut-off \tSensitivity\tSpecificity\tEfficiency\n')
            fprintf('%0.2f\t\t%0.4f\t\t%0.4f\t\t%0.4f\n',table')
            disp(tr)
            disp(' ')
        end
    end
    CSe=find(table(:,2)==max(table(:,2)),1,'first'); %Max sensitivity cut-off
    CSp=find(table(:,3)==max(table(:,3)),1,'last'); %Max specificity cut-off
    CEff=find(table(:,4)==max(table(:,4)),1,'first'); %Max efficiency cut-off
    d=realsqrt(xroc.^2+(1-yroc).^2); %apply the Pitagora's theorem
    [~,CE]=min(d); %Cost-effective cut-off
    xg=linspace(0,max(table(:,1)),500);
    st=[1 mean(table(:,1)) 1]; U=[Inf max(table(:,1)) Inf];
    fo_ = fitoptions('method','NonlinearLeastSquares','Lower',L,'Upper',U,'Startpoint',st);
    fitSe = fit(table(:,1),table(:,2),ft_,fo_);
    st=[-1 mean(table(:,1)) 1]; L=[-Inf 0 0]; U=[0 max(table(:,1)) Inf];
    fo_ = fitoptions('method','NonlinearLeastSquares','Lower',L,'Upper',U,'Startpoint',st);
    fitSp=fit(table(:,1),table(:,3),ft_,fo_);
    st=[min(table(:,4)) 1 mean(table(:,1)) max(table(:,4)) 1]; L=[0 0 0 0 0]; U=[1 Inf max(table(:,1)) 1 Inf];
    ft_ = fittype('D+(A-D)/((1+(x/C)^B)^E)','dependent',{'y'},'independent',{'x'},'coefficients',{'A', 'B', 'C', 'D', 'E'});
    fo_ = fitoptions('method','NonlinearLeastSquares','Lower',L,'Upper',U,'Startpoint',st);
    fitEff=fit(table(:,1),table(:,4),ft_,fo_);
    if verbose
        H2=figure;
        set(H2,'Position',[570 402 868 420])
        hold on
        HSE = plot(xg,feval(fitSe,xg),'marker','none','linestyle','-','color','r','linewidth',2);
        HCSe=plot([table(CSe,1) table(CSe,1)],[0 1],'marker','none','linestyle','--','color','r','linewidth',2);
        HSP = plot(xg,feval(fitSp,xg),'marker','none','linestyle','-','color','g','linewidth',2);
        HCSp=plot([table(CSp,1) table(CSp,1)],[0 1],'marker','none','linestyle','--','color','g','linewidth',2);
        HEFF = plot(xg,feval(fitEff,xg),'marker','none','linestyle','-','color','b','linewidth',2);
        HCEff=plot([table(CEff,1) table(CEff,1)],[0 1],'marker','none','linestyle','--','color','b','linewidth',2);
        HCO=plot([table(CE,1) table(CE,1)],[0 1],'marker','none','linestyle','--','color','m','linewidth',2);
        hold off
        legend([HSE HCSe HSP HCSp HCO HEFF HCEff],...
            'Sensitivity',sprintf('Max Sensitivity cutoff: %0.4f',table(CSe,1)),...
            'Specificity',sprintf('Max Specificity cutoff: %0.4f',table(CSp,1)),...
            sprintf('Cost effective cutoff: %0.4f',table(CE,1)),...
            'Efficiency',sprintf('Max Efficiency cutoff: %0.4f',table(CEff,1)),...
            'Location','BestOutside')
        axis([xg(1) xg(end) 0 1.1])

        fprintf('1) Max Sensitivity Cut-off point= %0.2f\n',table(CSe,1))
        fprintf('2) Max Specificity Cut-off point= %0.2f\n',table(CSp,1))
        fprintf('3) Cost effective Cut-off point= %0.2f\n',table(CE,1)),
        fprintf('4) Max Efficiency Cut-off point= %0.2f\n',table(CEff,1))
        m=table([CSe CSp CE CEff],1);
    end
else
    table=NaN;
end

if nargout
    ROCout.AUC=Area; %Area under the curve
    ROCout.SE=Serror; %standard error of the area
    ROCout.ci=ci; % 95% Confidence interval
    ROCout.co=m; % cut off points
    ROCout.xr=xroc; %graphic x points
    ROCout.yr=yroc; %graphic y points
    ROCout.table=table;
end


text:
function [imref, im, mask] = histology_preprocessing(imref, im, opts)
% HISTOLOGY_PREPROCESSING  Prepare slices for intra-histology registration.
%
% HISTOLOGY_PREPROCESSING converts a list of histology images to grayscale,
% inverts and thresholds them (so that the background is black instead of
% white), extends the histograms to cover the dynamic range, and then
% matches the histograms to a reference slice. This prepares the slices for
% successful intra-histology registration.
%
% The function has two syntaxes. The reference can be provided as an image
% or histogram.
%
% -------------------------------------------------------------------------
% Reference image syntax:
% [IMREF2, IM2, MASK] = HISTOLOGY_PREPROCESSING(IMREF, IM)
%
%   IMREF is the slice whose histogram will be used as reference for the
%   others. IM is the stack of RGB colour or grayscale histology images to
%   preprocess. The images can be provided as:
%
%     * A char with a single path and filename. E.g. 'file.mha'. With this
%       format, an output directory must be provided in OPTS.outdir, or
%       else a temporal one with a random name is created (see below).
%
%     * A cell array of filenames. E.g. {'file1.mha', 'file2.mha'}. IM{I}
%       is the path and filename of the I-th slice. See note about
%       OPTS.outdir in previous item.
%
%     * A plain array IM(R, C, I, CH), R=row, C=column, I=slice,
%       CH=channel.
%
%     * A vector of scimat structs. IM(I) is a scimat struct with the I-th
%       slice. See "help scimat" for details.
%
%   IMREF2, IM2 are the preprocessed images in the same format as the
%   inputs.
%
%   MASK is a binary mask for IM2. This mask can be used to speed up
%   registration (when a mask is provided to elastix, only pixels within
%   the mask are used for the registration metric). The mask includes the
%   tissue and a bit of background around it.
%
% ... = HISTOLOGY_PREPROCESSING(..., OPTS)
%
%   OPTS is a struct with options for the algorithm.
%
%     'outdir': If IM is a list of filenames, this is the path where the
%               preprocessed images are saved to.
%
%     'Grayscale': (def false) true/false whether image will be converted
%              to grayscale.
%
% -------------------------------------------------------------------------
% Reference histogram syntax:
% [HREF2, IM2, MASK] = HISTOLOGY_PREPROCESSING(HREF, IM)
%
%   HREF is a (CH, L)-matrix with the reference histogram, where CH is the
%   number of channels, and L is the number of intensity levels. For
%   example, for an RGB image, CH=3, and for a uint8 class, L=256.
%
%   IM is the image to preprocess, as in the syntax above.
%
%   HREF2 is the histogram with the same preprocessing as the image.
%
%   IM2, MASK are as in the syntax above.
%
% ... = HISTOLOGY_PREPROCESSING(..., OPTS)
%
%   OPTS is as in the syntax above.
%
%
% See also: regmatchedfilt, transfdiffreg.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014-2016 University of Oxford
% Version: 0.5.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% 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 arguments
narginchk(2, 3);
nargoutchk(0, 3);

% defaults
if (nargin < 3 || isempty(opts))
    opts = struct;
end
if (~isfield(opts, 'Grayscale') || isempty(opts.Grayscale))
    opts.Grayscale = false;
end

% check the type and number of the input images
if (ischar(im)) % single filename

    imType = 'filename';
    im = {im};
    N = 1;

elseif (iscell(im)) % list of filenames

    imType = 'filename';
    N = length(im);

elseif (isstruct(im)) % list of scimat structs

    imType = 'scimat';
    N = length(im);

else % plain array

    imType = 'array';
    N = size(im, 3);

end

% convert IMREF to a plain array if it's not already, or keep as histogram
% if it's a histogram
if (ischar(imref)) % filename

    imrefMat = scimat_load(imref);
    imrefMat = imrefMat.data;
    imrefType = 'filename';

elseif (isstruct(imref)) % scimat struct

    imrefMat = imref.data;
    imrefType = 'scimat';

elseif (isa(imref, 'double') && size(imref, 1)==3) % histogram as (3,256) matrix

    imrefType = 'histogram';

else % array

    % move channels to 5th index, so that it's the same format as IM.data
    imrefMat = reshape(imref, size(imref, 1), size(imref, 2), 1, 1, size(imref, 3));
    imrefType = 'array';

end

switch (imrefType)
    case 'histogram'

        % is the histogram has more than one channel, we don't try to
        % convert to grayscale

        % preprocessing of histogram equivalent to preprocessing the
        % reference image (invert, threshold, etc.)
        imref = histogram_preprocessing(imref);

    otherwise

        % grayscale conversion of IMREF, if requested by user
        if (opts.Grayscale)

            imrefMat = rgb2gray(squeeze(imrefMat));

        end

        % preprocessing of the reference image (invert, threshold, etc.)
        imrefMat = individual_image_preprocessing(imrefMat);

end

% if images are provided as filenames, either the user must provide an
% output directory, or we generate a temp one with a random name
if (strcmp(imType, 'filename') ...
        && (~isfield(opts, 'outdir') || isempty(opts.outdir)))
    opts.outdir = tempname;
end

% allocate memory for masks
switch (imType)

    case  'filename' % filename

        mask = cell(1, N);

    case 'scimat' % scimat struct

        mask = im;

    case 'array' % plain array

        mask = zeros(size(im, 1), size(im, 2), size(im, 3), 'uint8');

end

% loop all images
for I = 1:N

    % convert IM to a scimat struct if it's not already
    switch (imType)

        case  'filename' % filename

            im0 = scimat_load(im{I});

        case 'scimat' % scimat struct

            im0 = im(I);

        case 'array' % plain array

            im0 = scimat_im2scimat(reshape(im(:, :, I, :, :), ...
                [size(im, 1) size(im, 2) 1 1 size(im, 4)]));

    end

    switch (imrefType)
        case 'histogram'

            % preprocessing of IM
            mask0 = im0;
            [im0.data, mask0.data] = individual_image_preprocessing(im0.data);

            im0.data = match_histograms_hist(imref, im0.data);

            % grayscale conversion of IM, if requested by user
            if (opts.Grayscale)

                im0.data = rgb2gray(squeeze(im0.data));

            end

        otherwise

            % grayscale conversion of IM, if requested by user
            if (opts.Grayscale)

                im0.data = rgb2gray(squeeze(im0.data));

            end

            % preprocessing of IM
            mask0 = im0;
            [im0.data, mask0.data] = individual_image_preprocessing(im0.data);

            % match IM histogram to IMREF
            im0.data = match_histograms_im(imrefMat, im0.data);

    end

    % convert processed images to output format
    switch (imType)

        case  'filename' % filename

            [~, name, ext] = fileparts(im{I});
            im{I} = [opts.outdir filesep name ext];
            scimat_save(im{I}, im0);
            mask{I} = [opts.outdir filesep 'mask-' name ext];
            scimat_save(mask{I}, mask0);

        case 'scimat' % scimat struct

            im(I) = im0;
            mask(I) = mask0;

        case 'array' % plain array

            im(:, :, I, :) = reshape(im0.data, ...
                [size(im, 1) size(im, 2) 1 size(im, 5)]);
            mask(:, :, I, :) = mask0.data;

    end

end

% convert processed images to output format
switch (imrefType)

    case  'filename' % filename

        error('Not implemented yet')

    case 'scimat' % scimat struct

        imref.data = imrefMat;

    case 'array' % plain array

        imref = reshape(imrefMat, ...
            [size(imrefMat, 1) size(imrefMat, 2) 1 size(imrefMat, 5)]);

end

end

%% nested functions

% individual_image_preprocessing:
%
%   Preprocessing of an individual image. Invert, threshold background,
%   extend histogram to whole dynamic range. Compute a rough binary mask of
%   the object.
function [im, mask] = individual_image_preprocessing(im)

% loop image channels
for I = 1:size(im, 5)

    % select one channel
    ch = im(:, :, :, :, I);

    % the image should have two types of blackground pixels: white-ish and
    % pure black. The white-ish ones are from the original image, and the
    % black ones come from the B-spline fill-in. When we invert the image,
    % we don't want to invert the black ones, so we select non-black pixels
    idx = ch ~= 0;

    % invert image and make lowest intensity = 0
    ch(idx) = max(ch(idx)) - ch(idx);

    % remove background by zeroing anything <= mode (the background forms
    % the largest peak). We are quite conservative here in terms of keeping
    % all the tissue we can, to avoid losing detail, even if that means
    % that we are going to have a bit of background noise
    ch(ch <= mode(double(ch(idx)))) = 0;

    % extend histogram to cover whole dynamic range
    minh = double(min(ch(ch > 0)));
    maxh = double(max(ch(ch > 0)));
    ch = uint8(255 * (double(ch) - minh) / (maxh - minh));

    % replace image channel with the processed one
    im(:, :, :, :, I) = ch;

end

% compute a mask of the foreground
if (size(im, 5) == 1)

    mask = uint8(im > 2);

else

    mask = uint8(rgb2gray(squeeze(im)) > 2);

end

% remove some background noise in the mask
se = strel('disk', 1);
mask = imerode(mask, se);

% dilate mask to cover some background
se = strel('disk', 10);
mask = imdilate(mask, se);

end

% histogram_preprocessing:
%
%   Transform histogram to replicate what individual_image_preprocessing()
%   does to an image
function hgram = histogram_preprocessing(hgram)

% invert histogram, except for the pixels == 0, that continue being 0
hgram(:, 2:end) = hgram(:, end:-1:2);

% find the modes of the histogram
[~, imode] = max(hgram, [], 2);

% number of intensity levels
L = size(hgram, 2);

% loop histogram channels
for CH = 1:size(hgram, 1)

    % crop the histogram <= mode
    hgram(CH, 1) = sum(hgram(CH, 1:imode(CH)));
    hgram(CH, 2:imode(CH)) = 0;

    % extend dynamic range of histogram
    xbin = linspace(1, L-1, L-imode(CH));
    hgram(CH, 2:end) = round(interp1(xbin, hgram(CH, (imode(CH)+1):end), 1:L-1));

end

end

% match_histograms_im:
%
%   Match the histogram of each channel in IM to IMREF, ignoring background
%   voxels. IMREF is a reference image.
%
%   IMREF: reference image, array (R, C, CH)
function im = match_histograms_im(imref, im)

% match histograms to reference slice
for I = 1:size(im, 5)

    % select channels
    chref = imref(:, :, :, :, I);
    ch = im(:, :, :, :, I);

    % match histograms ignoring the background
    idxref = chref > 0;
    idx = ch > 0;
    ch(idx) = imhistmatch(ch(idx), chref(idxref));

    % replace image channel with processed ones
    %imref(:, :, I) = chref;
    im(:, :, :, :, I) = ch;

end

end

% match_histograms_hist:
%
%   Match the histogram of each channel in IM to IMREF, ignoring background
%   voxels. IMREF is a histogram with one or more channels.
%
%   IMREF: histogram, array (CH, intensity)
function im = match_histograms_hist(hgram, im)

% match image histograms to reference histogram
for I = 1:size(im, 5)

    % select image channel
    ch = im(:, :, :, :, I);

    % match histograms ignoring the background
    idx = ch > 0;
    ch(idx) = histeq(ch(idx), hgram(I, 2:end));

    % replace image channel with processed ones
    %imref(:, :, I) = chref;
    im(:, :, :, :, I) = ch;

end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus