text:
% Originally Coded By: Niall Mangan
% Copyright 2017, All Rights Reserved
% Code by Niall Mangan for paper "Inferring biological networks by sparse
% identification of nonlinear dynamics"
% by N. M. Mangan S. L. Brunton, J. L. Proctor, and J. N. Kutz
% Original Code: https://github.com/niallmm/iSINDy
%%
% This file is modified to compare the performance of i-SINDy on
% insufficient data.
% Modified by:K, 2019/07/16
% Find equations for yeast glycolysis state variable 6
%%
clc;clear all;close all;

addpath('./utils');
addpath('./bioutils');
FolderName='Results2';

% define libarary parameters
laurentorder = 0;
polyorder = 6;
usesine = 0;
dyorder = 1;

% clear variables from other solutions.
clear Theta Thetastring Xi indTheta lambdavec numterms errorv indopt
clear indTheta1 Xi1 numterms1 nT

% Here we load the previously simulated data
load('TrainingData.mat')

%% Define some parameters for the implicit SINDy
% Initial lambda value, which is the value used for soft thresholding in ADM
% To make implicit SINDy work, you need to carefully tune these two
% parameters. The below values are the best one that I found.
tol = 1e-4;
lambda = 3e-4;

% counter
jj = 1;

% initialize the number of nonzero terms found for the lambda
num= 1;
errorvec= 0;

% Define maximum iteration you need
MaxIter = 1000;

% Plot option, we set it to zero
plottag=0;

% Define the for loop
percent_start=0.3;
percent_end=1;
d_percent=0.1;

% Define how many results you want to get for each percentage
N_Iter=20;

for iter=1:N_Iter
    for percent=percent_start:d_percent:percent_end
        fprintf('\n\n\t Using %i percent of the data...\n',percent*100)

        fprintf('\n\t\t Calculating for the %i time...\n',iter)
        % Get the data length
        new_length=round(percent*length(xt));

        % Shuffel the original data
        Sequence=randperm(size(xt,1));
        xt_dum=xt(Sequence,:);
        dxt_dum=dxt(Sequence,:);

        % Assign the value to the new variables
        Data=xt_dum(1:new_length,:);dData=dxt_dum(1:new_length,:);

        % Define the number of states
        n=size(Data,2);

        % pool Data  (i.e., build library of nonlinear time series)
        [Theta, Thetastring] = poolDatady(Data,n,polyorder,usesine, laurentorder, dData(:,6), dyorder);

        % for now calculate null space using null function
        nT = null(Theta);

        % Get the sparse vector
        tic
        [indTheta1, Xi1, numterms1] = ADMinitvary(nT,lambda,MaxIter,tol, plottag);
        fprintf('\t\t\t Calculation finished! Used %d seconds.\n',toc)

        % Save the calculation result of current iteration
        fprintf('\n\t\t\t\t Saving the result...\n')
        tic
        cc=clock;
        ResultName=strcat(FolderName,'/implicit_SINDY_Data_Length_',num2str(percent*100),'_','Iter',num2str(iter),'_',num2str(cc(3)),'_',num2str(cc(4)),'_',num2str(cc(5)),'_',num2str(round(cc(6))),'.mat');
        save(ResultName,'Xi1','indTheta1','percent')
        fprintf('\n\t\t\t\t Saving finished! Using %i seconds...\n',toc)

    end
end

author: dynamicslab

repo: SINDy-PI


text:
function [x,y] = intline(x1,x2,y1,y2)
%INTLINE Integer-coordinate line drawing algorithm.
%   [X,Y] = INTLINE(X1,X2,Y1,Y2) computes an approximation to the line
%   segment joining (X1,Y1) and (X2,Y2) with integer coordinates. X1,
%   X2, Y1, and Y2 should be integers. INTLINE is reversible; that is,
%   INTLINE(X1,X2,Y1,Y2) produces the same results as
%   FLIPUD(INTLINE(X2,X1,Y2,Y1)).

%   Copyright 1993-2019 The MathWorks, Inc.
%   License: https://github.com/mathworks/matlab-color-tools/blob/master/license.txt

dx = abs(x2 - x1);
dy = abs(y2 - y1);

% Check for degenerate case.
if ((dx == 0) && (dy == 0))
   x = x1;
   y = y1;
  return;
end

flip = 0;
if (dx >= dy)
   if (x1 > x2)
      % Always "draw" from left to right.
      t = x1; x1 = x2; x2 = t;
      t = y1; y1 = y2; y2 = t;
      flip = 1;
   end
   m = (y2 - y1)/(x2 - x1);
   x = (x1:x2).';
   y = round(y1 + m*(x - x1));
else
   if (y1 > y2)
      % Always "draw" from bottom to top.
      t = x1; x1 = x2; x2 = t;
      t = y1; y1 = y2; y2 = t;
      flip = 1;
   end
   m = (x2 - x1)/(y2 - y1);
   y = (y1:y2).';
   x = round(x1 + m*(y - y1));
end

if (flip)
   x = flipud(x);
   y = flipud(y);
end

author: dipum

repo: dipum-toolbox

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


text:
function fname_out = savechgtab(fname, chgtab, warnings)
%SAVECHGTAB  Save a change table to a file.
%   SAVECHGTAB(FNAME, CHGTAB)
%   SAVECHGTAB(FNAME, CHGTAB, WARNINGS)
%   FNAME = SAVECHGTAB(FNAME, ...)
%
%   Writes a CHGTAB, suitable for use with APPLY_CHANGES to a file
%   specified by FNAME. If FNAME ends with '.mat' it saves CHGTAB
%   and WARNINGS to a MAT-file as the variables 'chgtab' and 'warnings',
%   respectively. Otherwise, it saves an M-file function that returns
%   the CHGTAB, with the optional WARNINGS in comments.
%
%   Optionally returns the filename, with extension added if necessary.
%
%   Input:
%       FNAME :  name of the file to be saved
%       CHGTAB :  change table suitable for use with APPLY_CHANGES
%       WARNINGS : optional cell array of warning messages (to be
%                   included in comments), such as those returned by
%                   PSSECON2CHGTAB
%
%   Output(s):
%       FNAME :  name of the file, with extention added if necessary

%   MATPOWER
%   Copyright (c) 2017, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% handle input args
if nargin < 3
    warnings = {};
end

%% define named indices into data matrices
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% constant names
bus_cols = {'BUS_I', 'BUS_TYPE', 'PD', 'QD', 'GS', 'BS', 'BUS_AREA', ...
    'VM', 'VA', 'BASE_KV', 'ZONE', 'VMAX', 'VMIN'};
br_cols = {'F_BUS', 'T_BUS', 'BR_R', 'BR_X', 'BR_B', 'RATE_A', ...
    'RATE_B', 'RATE_C', 'TAP', 'SHIFT', 'BR_STATUS', 'ANGMIN', 'ANGMAX'};
gen_cols = {'GEN_BUS', 'PG', 'QG', 'QMAX', 'QMIN', 'VG', 'MBASE', ...
    'GEN_STATUS', 'PMAX', 'PMIN', 'PC1', 'PC2', 'QC1MIN', 'QC1MAX', ...
    'QC2MIN', 'QC2MAX', 'RAMP_AGC', 'RAMP_10', 'RAMP_30', 'RAMP_Q', 'APF'};
gc_cols = {'CT_MODCOST_X', 'CT_MODCOST_F', '0', ...  %% use val+3
    'MODEL', 'STARTUP', 'SHUTDOWN', 'NCOST', 'COST'};
load_cols = {'CT_LOAD_ALL_PQ', 'CT_LOAD_FIX_PQ', 'CT_LOAD_DIS_PQ', ...
    'CT_LOAD_ALL_P', 'CT_LOAD_FIX_P', 'CT_LOAD_DIS_P'};
chgtypes = {'CT_REP', 'CT_REL', 'CT_ADD'};

%% verify valid filename
[pathstr, fcn_name, extension] = fileparts(fname);
if isempty(extension)
    extension = '.m';
end
if regexp(fcn_name, '\W')
    old_fcn_name = fcn_name;
    fcn_name = regexprep(fcn_name, '\W', '_');
    fprintf('WARNING: ''%s'' is not a valid function name, changed to ''%s''\n', old_fcn_name, fcn_name);
end
fname = fullfile(pathstr, [fcn_name extension]);

%% open and write the file
if strcmp(upper(extension), '.MAT')     %% MAT-file
    save(fname, 'chgtab', 'warnings');
else                                %% M-file
    %% open file
    [fd, msg] = fopen(fname, 'wt');     %% print it to an M-file
    if fd == -1
        error(['savechgtab: ', msg]);
    end

    %% function header, etc.
    fprintf(fd, 'function chgtab = %s\n', fcn_name);
    fprintf(fd, '%%%s  Returns a change table suitable for use with APPLY_CHANGES\n', upper(fcn_name));

    if ~isempty(warnings)
        fprintf(fd, '%%\n');
        fprintf(fd, '%% Warnings:\n');
        for k = 1:length(warnings)
            fprintf(fd, '%%     %s\n', warnings{k});
        end
    end

    fprintf(fd, '\ndefine_constants;\n');

    fprintf(fd, '\n%%%% Change Table\n');
    fprintf(fd, '%%\tlabel\tprob\ttable\trow\tcol\tchgtype\tnewval\n');
    fprintf(fd, 'chgtab = [\n');
    for c = 1:size(chgtab, 1)
        fprintf(fd, '\t%d\t%g\t', ...       %% CT_LABEL, CT_PROB
            chgtab(c, CT_LABEL), chgtab(c, CT_PROB)');
        switch chgtab(c, CT_TABLE)          %% CT_TABLE
            case CT_TBUS
                fprintf(fd, '%s\t', 'CT_TBUS');
                col = bus_cols{chgtab(c, CT_COL)};
            case CT_TGEN
                fprintf(fd, '%s\t', 'CT_TGEN');
                col = gen_cols{chgtab(c, CT_COL)};
            case CT_TGENCOST
                fprintf(fd, '%s\t', 'CT_TGENCOST');
                if chgtab(c, CT_COL)+3 > 8
                    col = sprintf('%s+%d', gc_cols{5+3}, chgtab(c, CT_COL)-5);
                else
                    col = gc_cols{chgtab(c, CT_COL)+3};
                end
            case CT_TBRCH
                fprintf(fd, '%s\t', 'CT_TBRCH');
                col = br_cols{chgtab(c, CT_COL)};
            case CT_TLOAD
                fprintf(fd, '%s\t', 'CT_TLOAD');
                col = load_cols{abs(chgtab(c, CT_COL))};
                if chgtab(c, CT_COL) < 0
                    col = ['-' col];
                end
            case CT_TAREABUS
                fprintf(fd, '%s\t', 'CT_TAREABUS');
                col = bus_cols{chgtab(c, CT_COL)};
            case CT_TAREAGEN
                fprintf(fd, '%s\t', 'CT_TAREAGEN');
                col = gen_cols{chgtab(c, CT_COL)};
            case CT_TAREAGENCOST
                fprintf(fd, '%s\t', 'CT_TAREAGENCOST');
                if chgtab(c, CT_COL)+3 > 8
                    col = sprintf('%s+%d', gc_cols{5+3}, chgtab(c, CT_COL)-5);
                else
                    col = gc_cols{chgtab(c, CT_COL)+3};
                end
            case CT_TAREABRCH
                fprintf(fd, '%s\t', 'CT_TAREABRCH');
                col = br_cols{chgtab(c, CT_COL)};
            case CT_TAREALOAD
                fprintf(fd, '%s\t', 'CT_TAREALOAD');
                col = load_cols{abs(chgtab(c, CT_COL))};
                if chgtab(c, CT_COL) < 0
                    col = ['-' col];
                end
        end
        fprintf(fd, '%g\t%s\t%s\t%g;\n', ...%% CT_ROW, %% CT_COL, CT_CHGTYPE, CT_NEWVAL
            chgtab(c, CT_ROW), col, chgtypes{chgtab(c, CT_CHGTYPE)}, chgtab(c, CT_NEWVAL));
    end
    fprintf(fd, '];\n');

    %% close file
    if fd ~= 1
        fclose(fd);
    end
end

if nargout > 0
    fname_out = fname;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function closeAllCERRFigures(varargin)
%Close all figures related to CERR, with the exception of the transverse viewer, which closes itself.
%
% 15 Aug 02, V H Clark
% LM: 30 Dec 02, JOD
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS planC

if nargin == 0

	fieldsC = fieldnames(stateS.handle);

	for i = 1 : length(fieldsC)
      if ~strcmp(fieldsC{i},'CERRSliceViewer')
        hV = getfield(stateS.handle,fieldsC{i});
        for j = 1 : length(hV)
          try
            delete(hV(j))
          end
        end
      end
	end

	if ~stateS.workspacePlan & stateS.planLoaded
        clear global planC
	end
	clear global stateS

%If a selective close, only close certain figures.
elseif nargin == 1 & strcmpi(varargin{1},'Selective');
    try, delete(stateS.handle.doseManagementFig);,end;
    try, delete(stateS.handle.doseSubtractionMenuFig);,end;
    try, delete(stateS.handle.IMRTMenuFig);,end;
    try, delete(stateS.handle.DVHMenuFigure);,end;
    try,
        for i=1:length(stateS.handle.DVHPlot)
            try, delete(stateS.handle.DVHPlot),end;
        end
    end
    try, delete(stateS.handle.metricSelectionFig);,end;
    try, delete(stateS.handle.doseShadowFig);,end;
    try, delete(stateS.handle.graphicalComparisonFig);,end;
    try, delete(stateS.handle.aboutCERRFig);,end;
    try, delete(stateS.handle.structureFusionFig);,end;
    try, delete(stateS.handle.navigationMontage);,end;
    try, delete(stateS.handle.deepLearnSegFig);,end;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pde = TorusTimeInitial4(mum,mup,sigm,sigp,epsm,epsp,omega,x0,y0,z0,r1,r2,a,b,intPt)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'Mu',@Mu,'Mum',@Mum,'Mup',@Mup,'one',@one,...
    'Epslon',@Epslon,'Epslonm',@Epslonm,'Epslonp',@Epslonp,...
    'Sig',@Sig,'Sigm',@Sigm,'Sigp',@Sigp,...
    'E1',@E1,'E2',@E2,'E3',@E3,...
    'Et1',@Et1,'Et2',@Et2,'Et3',@Et3);

pde.mum = mum;
pde.mup = mup;
pde.sigm = sigm;
pde.sigp = sigp;
pde.epsm = epsm;
pde.epsp = epsp;
%% interface function
    function u = intf(x,y,z)
        u = (sqrt((x-x0).^2+(y-y0).^2)-r2).^2 + (z-z0).^2 - r1^2;
    end

%% exact solution
    function u = exactu1(x,y,z,t)
        u = um1(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id),t);
        %u = t*ones(size(x));
    end
    function u = exactu2(x,y,z,t)
        u = um2(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id),t);
        %u = t*ones(size(x));
    end
    function u = exactu3(x,y,z,t)
        u = um3(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id),t);
        %u = t*ones(size(x));
    end
    function u = um1(x,y,z,t)
        u = zeros(size(x));
    end
    function u = um2(x,y,z,t)
        u = exp(-b*(a*(x-intPt)-omega*t).^2);
    end
    function u = um3(x,y,z,t)
        u = zeros(size(x));
    end
    function u = up1(x,y,z,t)
        u = zeros(size(x));
        %u = t*ones(size(x));
    end
    function u = up2(x,y,z,t)
        u = exp(-b*(a*(x-intPt)-omega*t).^2);
    end
    function u = up3(x,y,z,t)
        u = zeros(size(x));
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z,t)
        u = Dxum(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id),t);
    end
    function u = Dyu(x,y,z,t)
        u = Dyum(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id),t);
    end
    function u = Dzu(x,y,z,t)
        u = Dzum(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id),t);
    end
    function u = Dxum(x,y,z,t)
        u = zeros(size(x));
    end
    function u = Dyum(x,y,z,t)
        u = zeros(size(x));
    end
    function u = Dzum(x,y,z,t)
        u = -exp(-b*(a*(x-intPt)-omega*t).^2)*2*a*b.*(a*(x-intPt)-omega*t);
    end
    function u = Dxup(x,y,z,t)
        u = zeros(size(x));
    end
    function u = Dyup(x,y,z,t)
        u = zeros(size(x));
    end
    function u = Dzup(x,y,z,t)
        u = -exp(-b*(a*(x-intPt)-omega*t).^2)*2*a*b.*(a*(x-intPt)-omega*t);
    end

%% right hand side function
    function u = f1(x,y,z,t)
        u = fm1(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id),t);
    end
    function u = f2(x,y,z,t)
        u = fm2(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id),t);
    end
    function u = f3(x,y,z,t)
        u = fm3(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id),t);
    end

    function u = fm1(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fm2(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fm3(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp1(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp2(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp3(x,y,z,t)
        u = zeros(size(x));
    end

%% Diffusion coefficient function
    function u = Mu(x,y,z)
        u = Mum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Mup(x(id),y(id),z(id));
    end
    function u = Mum(x,y,z)
        u = mum^(-1)*ones(size(x));
    end
    function u = Mup(x,y,z)
        u = mup^(-1)*ones(size(x));
    end

%% Mass coefficient function
    function u = Epslon(x,y,z)
        u = Epslonm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Epslonp(x(id),y(id),z(id));
    end
    function u = Epslonm(x,y,z)
        u = epsm*ones(size(x));
    end
    function u = Epslonp(x,y,z)
        u = epsp*ones(size(x));
    end

    function u = Sig(x,y,z)
        u = Sigm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Sigp(x(id),y(id),z(id));
    end
    function u = Sigm(x,y,z)
        u = sigm*ones(size(x));
    end
    function u = Sigp(x,y,z)
        u = sigp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function I1 = Normalize_Img(I)

if ndims(I) == 3
    I = rgb2gray(I);
end
I1 = imresize(I, [24 24], 'bicubic');

author: decouples


text:
function [ t_stop, y_stop ] = p35_stop ( neqn )

%*****************************************************************************80
%
%% P35_STOP returns the stopping point for problem p35.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
   -0.584253e-04; ...
    0.359969e-03 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = sldrawgraph(G, X, Xt, ppedges, ppsn, pptn)
%SLDRAWGRAPH Draws a graph
%
% $ Syntax $
%   - sldrawgraph(G, X)
%   - sldrawgraph(G, X, Xt)
%   - sldrawgraph(G, X, [], ppedges)
%   - sldrawgraph(G, X, [], ppedges, ppsn)
%   - sldrawgraph(G, X, Xt, ppedges)
%   - sldrawgraph(G, X, Xt, ppedges, pptn)
%   - h = sldrawgraph(...)
%
% $ Arguments $
%   - G:        The graph in any acceptable form
%   - X:        The sample matrix of (source) nodes (2xn or 3xn)
%   - Xt:       The sample matrix of (target) nodes (2xn or 3xn)
%   - ppedges:  The cell array of parameters for plotting graph edges
%   - ppsn:     The cell array of parameters for plotting source nodes
%   - pptn:     The cell array of parameters for plotting target nodes
%   - h:        The vector of handles to all plotted objects
%
% $ Description $
%   - sldrawgraph(G, X) draws a graph with the coordinates of the
%     graph nodes given in X (a 2 x n matrix for 2D points or 3 x n
%     matrix for 3D points), using default plotting parameters.
%
%   - sldrawgraph(G, X, Xt) draws a bigraph using default plotting
%     parameters.
%
%   - sldrawgraph(G, X, [], ppedges) draws a graph using the specified
%     parameters to draw edges.
%
%   - sldrawgraph(G, X, [], ppedges, ppsn) draws a graph using the
%     specified parameters to draw edges and nodes.
%
%   - sldrawgraph(G, X, Xt, ppedges) draws a bigraph using the specified
%     parameters to draw edges.
%
%   - sldrawgraph(G, X, Xt, ppedges, ppsn, pptn) draws a bigraph using the
%     specified parameters to draw edges and nodes.
%
%   - h = sldrawgraph(...) returns the handles to the plotted objects.
%
% $ History $
%   - Created by Dahua Lin, on Sep 11st, 2006
%

%% parse and verify input arguments

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

EG = sledgeset(G, 'off');
gi = slgraphinfo(EG);
n = gi.n;
nt = gi.nt;

d = size(X, 1);
if ~isnumeric(X) || ndims(X) ~= 2 || (d ~= 2 && d ~= 3)
    error('sltoolbox:invalidarg', ...
        'The X should be a 2D numeric matrix, and should have 2 or 3 rows.');
end

if nargin < 3 || isempty(Xt)
    Xt = X;
    bi = false;
else
    if ~isnumeric(Xt) || ndims(Xt) ~= 2
        error('sltoolbox:invalidarg', ...
            'The X should be a 2D numeric matrix, and should have 2 or 3 rows.');
    end
    if size(Xt, 1) ~= d
        error('sltoolbox:sizmismatch', ...
            'The sample dimensions in Xt is not the same as X');
    end
    bi = true;
end

if size(X, 2) ~= n || size(Xt, 2) ~= nt
    error('sltoolbox:sizmismatch', ...
        'The size of graph is not consistent with the sample numbers');
end

if nargin < 4 || isempty(ppedges)
    ppedges = {};
end

if nargin < 5 || isempty(ppsn)
    ppsn = {};
end

if nargin < 6 || isempty(pptn)
    pptn = ppsn;
end

if nargout >= 1
    output_h = true;
else
    output_h = false;
end

%% Prepare Data for plotting

% produce pruned NaN-separated coordinate list

edges = EG.edges;
edges = slpruneedgeset(n, nt, edges);

I = edges(:,1);
J = edges(:,2);
nedges = length(I);

xc = [X(1, I); Xt(1, J); NaN(1, nedges)];
xc = xc(:);
yc = [X(2, I); Xt(2, J); NaN(1, nedges)];
yc = yc(:);
if d == 3
    zc = [X(3, I); Xt(3, J); NaN(1, nedges)];
    zc = zc(:);
end

%% plot edges

% plot the edges
if d == 2
    ch = plot(xc, yc, ppedges{:});
else
    ch = plot3(xc, yc, zc, ppedges{:});
end
if output_h
    h = ch;
end
clear xc yc zc;

%% plot source nodes

if ~isempty(ppsn)
    ppsn = [ppsn, {'LineStyle', 'none'}];
    hold on;
    if d == 2
        ch = plot(X(1,:), X(2,:), ppsn{:});
    else
        ch = plot3(X(1,:), X(2,:), X(3,:), ppsn{:});
    end
    if output_h
        h = [h; ch];
    end
end

%% plot target nodes

if bi && ~isempty(pptn)
    pptn = [pptn, {'LineStyle', 'none'}];
    hold on;
    if d == 2
        ch = plot(Xt(1,:), Xt(2,:), pptn{:});
    else
        ch = plot3(Xt(1,:), Xt(2,:), Xt(3,:), pptn{:});
    end
    if output_h
        h = [h; ch];
    end
end

author: lmthang


text:
%======================================================================
%SLIC demo
% Copyright (C) 2015 Ecole Polytechnique Federale de Lausanne
% File created by Radhakrishna Achanta
% Please also read the copyright notice in the file slicmex.c
%======================================================================
%Input parameters are:
%[1] 8 bit images (color or grayscale)
%[2] Number of required superpixels (optional, default is 200)
%[3] Compactness factor (optional, default is 10)
%
%Ouputs are:
%[1] labels (in raster scan order)
%[2] number of labels in the image (same as the number of returned
%superpixels
%
%NOTES:
%[1] number of returned superpixels may be different from the input
%number of superpixels.
%[2] you must compile the C file using mex slicmex.c before using the code
%below
%======================================================================
%img = imread('someimage.jpg');
% img = imread('bee.jpg');
% [labels, numlabels] = slicmex(img,500,20);%numlabels is the same as number of superpixels
img = imread('/srv/glusterfs/daid/data/cityscape/leftImg8bit/train/aachen/aachen_000001_000019_leftImg8bit.png');
[labels, numlabels] = slicmex(img,2000,10);%numlabels is the same as number of superpixels
imagesc(labels);

author: sakaridis


text:
%SVC Trainable classifier: Support Vector Machine
%
% 	[W,J] = SVC(A,KERNEL,C)
%   [W,J] = A*SVC([],KERNEL,C)
%   [W,J] = A*SVC(KERNEL,C)
%
% INPUT
%   A	      Dataset
%   KERNEL  - Untrained mapping to compute kernel by A*(A*KERNEL) during
%             training, or B*(A*KERNEL) during testing with dataset B.
%           - String to compute kernel matrices by FEVAL(KERNEL,B,A)
%           Default: linear kernel (PROXM('p',1))
%   C       Regularization parameter (optional; default: 1)
%
% OUTPUT
%   W       Mapping: Support Vector Classifier
%   J       Object indices of support objects
%
% DESCRIPTION
% Optimizes a support vector classifier for the dataset A by quadratic
% programming. The non-linearity is determined by the kernel.
% If KERNEL = 0 it is assumed that A is already the kernelmatrix (square).
% In this case also a kernel matrix B should be supplied at evaluation by
% B*W or PRMAP(B,W).
%
% There are several ways to define KERNEL, e.g. PROXM('r',1) for a
% radial basis kernel or by USERKERNEL for a user defined kernel.
%
% If C is NaN this regularisation parameter is optimised by REGOPTC.
%
% See for more possibilties SVCINFO
%
% EXAMPLE
% a = gendatb;                     % generate banana classes
% [w,J] = a*svc(proxm('p',3));  % compute svm with 3rd order polynomial
% a*w*testc                        % show error on train set
% scatterd(a)                      % show scatterplot
% plotc(w)                         % plot classifier
% hold on;
% scatterd(a(J,:),'o')             % show support objcts
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, PROXM, USERKERNEL, NUSVC, RBSVC, LIBSVC, REGOPTC

% Copyright: D. de Ridder, D. Tax, S. Verzakov, R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function [W,J,C,nu,alginf] = svc(varargin)

  mapname = 'SVC';
  argin = shiftargin(varargin,{'prmapping','char','cell'});
  argin = setdefaults(argin,[],proxm([],'p',1),1,1);

  if mapping_task(argin,'definition')

    W = define_mapping(argin,'untrained',mapname);

	elseif mapping_task(argin,'training')			% Train a mapping.

    [a,kernel,C,Options] = check_for_old_call(argin);
    if isequal(Options,1), Options = []; end

    DefOptions.mean_centering   = 1;
    DefOptions.pd_check         = 1;
    DefOptions.bias_in_admreg   = 1;
    DefOptions.pf_on_failure    = 1;
    DefOptions.multiclass_mode  = 'single';

    Options = updstruct(DefOptions,Options,1);

    islabtype(a,'crisp');
    isvaldfile(a,1,2); % at least 1 object per class, 2 classes
    a = testdatasize(a,'objects');
    [m,k,c] = getsize(a);
    nlab = getnlab(a);

    % The SVC is basically a 2-class classifier. More classes are
    % handled by mclassc.

    if c == 2   % two-class classifier

      if (isnan(C))|length(C)>1     % optimize trade-off parameter
        defs = {proxm([],'p',1),1,[]};
        if isnan(C) % use a default range of C values
          parmin_max = [0,0;1e-2,1e2;0,0];   % kernel and Options can not be optimised
        else % use the range that is supplied by the user
          parmin_max = [0 0; min(C) max(C); 0 0];
          C = NaN;
        end
        [W,J,C,nu,alginf] = regoptc(a,mfilename,{kernel,C,Options},defs,[2],parmin_max,testc([],'soft'));
      else
        % Compute the parameters for the optimization:
        y = 3 - 2*nlab;
        if isequal(kernel,0)
          s = [];
          u = [];
          in_size = 0; % to allow old and new style calls
        else
          if Options.mean_centering
            u = mean(a);         % shift origin for better accuracy
            a = a - repmat(u,[m,1]);
          else
            u = [];
          end
          in_size = k;
        end
        K = compute_kernel(a,a,kernel);
        K = min(K,K');   % make sure kernel is symmetric
        [v,J,C,nu] = svo(+K,y,C,Options);

        % Store the results:
        if ~isequal(kernel,0)
          s = a(J,:);
        end
        W = prmapping(mfilename,'trained',{u,s,v,kernel,J},getlablist(a),in_size,2);
        % Note: even in case kernel is a mapping, we store the untrained
        % mapping, as training is just administration
        W = cnormc(W,a);
        W = setcost(W,a);

        alginf.svc_type = 'C-SVM';
        alginf.kernel = kernel;
        alginf.C   = C;
        alginf.nu  = nu;
        alginf.nSV = length(J);
        alginf.classsizes = [nnz(y==1), nnz(y==-1)];
        alginf.pf = isnan(nu);
        W = setname(W,'SVM');

      end

    else   % multi-class classifier:

      [W,J,C,nu,alginf] = mclassc(a,prmapping(mfilename,{kernel,C,Options}),Options.multiclass_mode);
      W = W*classc;

    end
    W = setname(W,mapname);

  else % Evaluation

    [a,v] = deal(argin{1:2});
    w = +v;
    m = size(a,1);

    if isempty(w{2})
      K = a; % user supplied testkernel
      J = w{5};
      if size(K,2) > length(J) & size(K,2) >= max(J)
        K = K(:,J); % probably full test kernel
      elseif size(K,2) ~= length(J)
        error('Wrong test kernel supplied')
      end
    else
      % The first parameter w{1} stores the mean of the dataset. When it
      % is supplied, remove it from the dataset to improve the numerical
      % precision. Then compute the kernel matrix using proxm:
      if isempty(w{1})
        K = a*(w{2}*w{4});
      else
        K = (a-ones(m,1)*w{1})*(w{2}*w{4});
      end
    end
    % Data is mapped by the kernel, now we just have a linear
    % classifier  w*x+b:
    d = [K ones(m,1)] * w{3};
    W = setdat(a,[d -d],v);

  end

return

function K = compute_kernel(a,s,kernel)

	% compute a kernel matrix for the objects a w.r.t. the support objects s
	% given a kernel description

	if  ischar(kernel) % routine supplied to compute kernel
		K = feval(kernel,a,s);
	elseif iscell(kernel)
		K = feval(kernel{1},a,s,kernel{2:end});
	elseif ismapping(kernel)
		K = a*prmap(s,kernel);
	elseif kernel == 0 % we have already a kernel
		K = a;
	else
		error('Do not know how to compute kernel matrix')
	end

	K = +K;

return

function [a,kernel,C,par] = check_for_old_call(argin)

  [a,kernel,C,par] = deal(argin{:});
  if ischar(kernel) && exist(kernel,'file') ~= 2
    kernel = proxm(kernel,C);
    C = par;
    par = [];
  end

return

author: marianux

repo: ecg-kit

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


text:
function [xo]=groupthresh(xi,lambda,varargin)
%GROUPTHRESH   Group thresholding
%   Usage:  xo=groupthresh(xi,lambda);
%
%   `groupthresh(x,lambda)` performs group thresholding on *x*, with
%   threshold *lambda*.  *x* must be a two-dimensional array, the first
%   dimension labelling groups, and the second one labelling members. This
%   means that the groups are the row vectors of the input (the vectors
%   along the 2nd dimension).
%
%   Several types of grouping behaviour are available:
%
%   * `groupthresh(x,lambda,'group')` shrinks all coefficients within a given
%     group according to the value of the $l^2$ norm of the group in
%     comparison to the threshold *lambda*. This is the default.
%
%   * `groupthresh(x,lambda,'elite')` shrinks all coefficients within a
%     given group according to the value of the $l^1$ norm of the
%     group in comparison to the threshold value *lambda*.
%
%   `groupthresh(x,lambda,dim)` chooses groups along dimension
%   *dim*. The default value is $dim=2$.
%
%   `groupthresh` accepts all the flags of |thresh| to choose the
%   thresholding type within each group and the output type (full / sparse
%   matrix). Please see the help of |thresh| for the available
%   options. Default is to use soft thresholding and full matrix output.
%
%   See also:  thresh
%
%   Demos:  demo_audioshrink
%
%   References: Kowalski08sparsity kowalski2009mixed yu2008audio

%   AUTHOR : Kai Siedenburg, Bruno Torresani.
%   REFERENCE: OK

if nargin<2
  error('Too few input parameters.');k
end;

if (prod(size(lambda))~=1 || ~isnumeric(lambda))
  error('lambda must be a scalar.');
end;

% Define initial value for flags and key/value pairs.
definput.import={'thresh','groupthresh'};
definput.importdefaults={'soft'};
definput.keyvals.dim=2;

[flags,keyvals,dim]=ltfatarghelper({'dim'},definput,varargin);

% kv.dim (the time or frequency selector) is handled by assert_sigreshape_pre
[xi,L,NbMembers,NbGroups,dim,permutedsize,order]=assert_sigreshape_pre(xi,[],dim,'GROUPTHRESH');

if flags.do_sparse
  xo = sparse(size(xi));
else
  xo = zeros(size(xi));
end;

if flags.do_group

  groupnorm = sqrt(sum(abs(xi).^2));
  w = thresh(groupnorm, lambda, flags.iofun,flags.outclass)./groupnorm;

  % Clean w for NaN. NaN appears if the input has a group with norm
  % exactly 0.
  w(isnan(w)) = 0;

  xo = bsxfun(@times,xi,w);

end

if flags.do_elite
  for ii=1:NbGroups,
    y = sort(abs(xi(:,ii)),'descend');
    rhs = cumsum(y);
    rhs = rhs .* lambda ./ (1 + lambda * (1:NbMembers)');
    M_ii = find(diff(sign(y-rhs)));
    if (M_ii~=0)
      tau_ii = lambda * norm(y(1:M_ii),1)/(1+lambda*M_ii);
    else
      tau_ii = 0;
    end

    % FIXME: The following line does not work for sparse matrices.
    xo(:,ii) = thresh(xi(:,ii),tau_ii,flags.iofun,flags.outclass);
  end
end;

xo=assert_sigreshape_post(xo,dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ y, m ] = month_carry_islamic ( y, m )

%*****************************************************************************80
%
%% MONTH_CARRY_ISLAMIC carries a year of months on the Islamic calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the year and month.
%    On output, M is no greater than 12.
%
  while ( 1 )

    months = year_length_months_islamic ( y );

    if ( m <= months )
      break
    end

    m = m - months;
    y = y + 1;

  end

  return
end


text:
function Value = GetArg(S, Name, Default)

try
    Value = S.(Name);
catch
    error('GetArg: unknown parameter ''%s''.',Name);
end

if isempty(Value)
    if nargin < 3
        error('Parameter ''%s'' must be assigned a value.',Name);
    end
    Value = Default;
end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function latin_cover_test ( )

%*****************************************************************************80
%
%% LATIN_COVER_TEST tests the LATIN_COVER library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 August 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATIN_COVER_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LATIN_COVER library.\n' );

  latin_cover_test01 ( );
  latin_cover_test02 ( );
  latin_cover_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATIN_COVER_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Simulation code for [1], presented in ICIP'09 (http://www.icip2009.org/)
%
% [1] Paul Rodriguez and Brendt Wohlberg, "A Generalized Vector-Valued Total
%     Variation Algorithm", in Proceedings of IEEE International Conference
%     on Image Processing (ICIP), (Cairo, Egypt), doi:10.1109/ICIP.2009.5413587 ,
%     pp. 1309--1312, Nov 2009
%
% Legal:
%   irnIcip09.m is based on NUMIPAD (http://numipad.sf.net). NUMIPAD is free
%   software, you can redistribute it and/or modify it under the terms of
%   the GNU General Public License (version 2).
%
%   The NUMIPAD library is being developed under U.S. Government contract
%   W-7405-ENG-36 for Los Alamos National Laboratory.
%
% Authors
%   Paul Rodriguez    prodrig@pucp.edu.pe
%   Brendt Wohlberg   brendt@tmail.lanl.gov

clear all;
%close all;

SHOW_IMGS = false;
%  SHOW_IMGS = true;

%  example = 'icip09'
%  example = 'l1deconv';
%  example = 'l2deconv';
example = 'l1denoise';
%  example = 'l2denoise';

extreme = 0;

nmppath;

BKS_CODE = exist('MainRestoration');
if( (BKS_CODE == 0) && ( strcmp(example, 'l1deconv') || strcmp(example, 'l2deconv') || strcmp(example, 'icip09') ) )
  disp('NOTE:');
  disp('  The function MainRestoration (code for [BKS]) is not in your path...');
  disp(sprintf('  disabling BKS simulations for %s.\n',example));
  disp('  [BKS] L. Bar, A. Brook, N. Sochen and N. Kiryati');
  disp('        "Deblurring of Color Images Corrupted by Impulsive Noise" ');
  disp('        IEEE Transactions on Image Processing, 16 (1101-1111), 2007');
end

BnG_CODE = exist('tvdenoise');
if( (BnG_CODE == 0) && ( strcmp(example, 'l2denoise') || strcmp(example, 'icip09') ) )
  disp('NOTE:');
  disp('  The function tvdenoise (code for [BnG]) is not in your path...');
  disp(sprintf('  disabling BnG simulations for %s.\n',example));
  disp('  [BnG] an implementation of the fast dual minimization of VTV [1].');
  disp('        Code may be downloaded from:');
  disp('        http://www.mathworks.fr/matlabcentral/fileexchange/16236');
  disp('  [1] X. Bresson and T. Chan');
  disp('      "Fast dual minimization of the vectorial total variation norm and');
  disp('       applications to color image processing"');
  disp('      Journal of Inverse Problems and Imaging, 2:4(455--484), 2008');
end

if( exist('color_imgs/peppers_color.png') && exist('color_imgs/mandrill_color.png') ...
    && exist('color_imgs/lena_color_256.png') )
  disp(sprintf('\nRunning %s simulation... \n', example));
else
  disp(' ');
  disp('One or more test images are not in the current directory...');
  disp('You may download them from:');
  disp('http://sites.google.com/a/istec.net/prodrig/Home/en/pubs');
  disp('look for "test images" under "A Generalized Vector-Valued Total ');
  disp('Variation Algorithm"');
  disp(' ');
  disp('Exiting simulation code...');
  return;
end

if( strcmp(example, 'icip09') )

  str_all = [];
  ICIP_FLAG = 1;

else

  ICIP_FLAG = 0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Input images        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pepImg = imread('peppers_color.png');
pepImg = double(pepImg)/255.0;

lenImg = imread('lena_color_256.png');
lenImg = double(lenImg)/255.0;

mdrilImg = imread('mandrill_color.png');
mdrilImg = double(mdrilImg)/255.0;

%%%%%%%%%%%%%%%%%%%%%%%%
%     Normalize        %
%%%%%%%%%%%%%%%%%%%%%%%%

Normalize = @(x) (x - min(x(:)))/(max(x(:)) - min(x(:)));

%%%%%%%%%%%%%%%%%%%%%%%%
%       kernels        %
%%%%%%%%%%%%%%%%%%%%%%%%

kernel_BSK = fspecial('disk',3.2);

K_BSK = @(x) imfilter(x, kernel_BSK, 'symmetric','conv');
KT_BSK = @(x) K_BSK(x);
KC_BSK = {K_BSK, KT_BSK};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Blurred & noisy images        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lenImgBlur = K_BSK(lenImg);

lenImgBlur_01L1 = imnoise(lenImgBlur, 'salt & pepper', 0.1);
lenImgBlur_03L1 = imnoise(lenImgBlur, 'salt & pepper', 0.3);

% ---

lenImgBlur_01L2 = imnoise(lenImgBlur, 'gaussian', 0, ...
                          0.01*max(lenImgBlur(:)) ); %NOTE sigma^2 in imnoise

lenImgBlur_005L2 = imnoise(lenImgBlur,'gaussian', 0, ...
                           0.0025*max(lenImgBlur(:)) );

lenImgBlur_001L2 = imnoise(lenImgBlur,'gaussian', 0, ...
                           0.0001*max(lenImgBlur(:)) );
% -- BSK example ---:
lenImgBlur_bskL2 = imnoise(lenImgBlur,'gaussian', 0, ...
                           0.00001*max(lenImgBlur(:)) );

% ---

mdrilImg_01L1 = imnoise(mdrilImg, 'salt & pepper', 0.1);
mdrilImg_03L1 = imnoise(mdrilImg, 'salt & pepper', 0.3);

mdrilImg_01L2 = imnoise(mdrilImg,  'gaussian', 0, ...
                        0.01*max(mdrilImg(:)) ); %NOTE sigma^2 in imnoise
mdrilImg_005L2 = imnoise(mdrilImg, 'gaussian', 0, ...
                         0.0025*max(mdrilImg(:)) );

% ---

pepImg_01L1 = imnoise(pepImg, 'salt & pepper', 0.1);
pepImg_03L1 = imnoise(pepImg, 'salt & pepper', 0.3);

pepImg_01L2 = imnoise(pepImg, 'gaussian', 0, 0.01*max(pepImg(:)) );
pepImg_005L2 = imnoise(pepImg, 'gaussian', 0, 0.0025*max(pepImg(:)) );

%-----------------------------------------------------------------------------

if( strcmp(example,'l1deconv') || strcmp(example,'icip09') )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L1 Deconvolved        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%noise level: 0.1

%-- IRN
lambda = 0.035;
pars = irntvInputPars('l1tv');

pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 4;
%  pars.pcgtol_ini   = 1e-4;
pars.U0           = lenImgBlur_01L1;

tic;
IRN_lenImgBlur_01L1 = irntv(lenImgBlur_01L1, KC_BSK, lambda, pars);
tirn_01L1n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_lenImgBlur_01L1) );
  axis image; axis off;
  title(sprintf('Deconvolved Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, IRN_lenImgBlur_01L1), tirn_01L1n0));
end

if(BKS_CODE)
%-- BSK MSTV (Mumford-Shah TV)

Params = SetParams;
Params.beta=0.5;              % check - same value as bar-2007-deblurring
Params.alpha=0.1;             % check
Params.epsilon=0.1;           % check
Params.gamma = 2*10^(-3);     % as "\mu in bar-2007-deblurring table IV

%  [uh_mstv,V] = MainRestoration(z, kernel, 'L1', 'MSTV', Params);
tic;
[MSTV_lenImgBlur_01L1, V_MSTV_01L1] = MainRestoration(lenImgBlur_01L1, ...
                                            kernel_BSK, 'L1', 'MSTV', Params);
tbsk_01L1n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( MSTV_lenImgBlur_01L1 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, MSTV_lenImgBlur_01L1), tbsk_01L1n0));
end

%-- BSK 1^1-TV
Params = SetParams;
Params.beta=0.1;              % check - same value as bar-2007-deblurring

tic;
[BKSL1_lenImgBlur_01L1, V_BKSL1_01L1] = MainRestoration(lenImgBlur_01L1, ...
                                              kernel_BSK, 'L1', 'L1', Params);
tbsk_01L1n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( BKSL1_lenImgBlur_01L1 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - BKSL1. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, BKSL1_lenImgBlur_01L1), tbsk_01L1n1));
end

end % _END_ if(BKS_CODE)
%------------------------------------------------------------------------------

%noise level: 0.3

lambda = 0.070;
pars = irntvInputPars('l1tv');

pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 6;
%  pars.pcgtol_ini   = 1e-4;
pars.U0           = lenImgBlur_03L1;

tic;
IRN_lenImgBlur_03L1 = irntv(lenImgBlur_03L1, KC_BSK, lambda, pars);
tirn_03L1n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_lenImgBlur_03L1) );
  axis image; axis off;
  title(sprintf('Deconvolved Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, IRN_lenImgBlur_03L1), tirn_03L1n0));
end

if(BKS_CODE)
%-- BSK MSTV (Mumford-Shah TV)

Params = SetParams;
Params.beta=1.1;            % as describe in bar-2007-deblurring, table IV
Params.alpha=0.5;           % as describe in bar-2007-deblurring, table IV
Params.epsilon=0.1;         % check
Params.gamma = 2*10^(-3);

%  [uh_mstv,V] = MainRestoration(z, kernel, 'L1', 'MSTV', Params);
tic;
[MSTV_lenImgBlur_03L1, V_MSTV_03L1] = MainRestoration(lenImgBlur_03L1, ...
                                            kernel_BSK, 'L1', 'MSTV', Params);
tbsk_03L1n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( MSTV_lenImgBlur_03L1 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, MSTV_lenImgBlur_03L1), tbsk_03L1n0));
end

%-- BSK L1

Params = SetParams;
Params.beta=0.2;            % as describe in bar-2007-deblurring, table IV

tic;
[BKSL1_lenImgBlur_03L1, V_MSTV_03L1] = MainRestoration(lenImgBlur_03L1, ...
                                              kernel_BSK, 'L1', 'L1', Params);
tbsk_03L1n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( BKSL1_lenImgBlur_03L1 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, BKSL1_lenImgBlur_03L1), tbsk_03L1n1));
end

end % _END_  if(BKS_CODE)

if(ICIP_FLAG == 0)

  disp('L1 Deconvolve Vector TV');
  disp(' ');
  disp('                          SNR (db)                         Time (s)');
  disp(' Img      Noise    VTV-IRN    MSTV    BKSL1          VTV IRN    MSTV    BKSL1');
  disp(' ');
  if(BKS_CODE)
  disp(sprintf('Lena       0.1      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_01L1), snr(lenImg, MSTV_lenImgBlur_01L1), snr(lenImg, BKSL1_lenImgBlur_01L1), ...
     tirn_01L1n0, tbsk_01L1n0, tbsk_01L1n1));
  disp(sprintf('           0.3      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_03L1), snr(lenImg, MSTV_lenImgBlur_03L1), snr(lenImg, BKSL1_lenImgBlur_03L1), ...
     tirn_03L1n0, tbsk_03L1n0, tbsk_03L1n1));
  else
  disp(sprintf('Lena       0.1      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_01L1), NaN, NaN, ...
     tirn_01L1n0, NaN, NaN));
  disp(sprintf('           0.3      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_03L1), NaN, NaN, ...
     tirn_03L1n0, NaN, NaN));

  end % _END_ if(BKS_CODE)

else % IF(icip)

  str = sprintf('L1 Deconvolve Vector TV\n\n');
  str_all = [str_all str];
  str = sprintf('                          SNR (db)                         Time (s)\n');
  str_all = [str_all str];
  str = sprintf(' Img      Noise    VTV-IRN    MSTV    BKSL1          VTV IRN    MSTV    BKSL1\n\n');
  str_all = [str_all str];

  if(BKS_CODE)
  str = sprintf('Lena       0.1      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f\n', ...
     snr(lenImg, IRN_lenImgBlur_01L1), snr(lenImg, MSTV_lenImgBlur_01L1), snr(lenImg, BKSL1_lenImgBlur_01L1), ...
     tirn_01L1n0, tbsk_01L1n0, tbsk_01L1n1);
  str_all = [str_all str];
  str = sprintf('           0.3      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f\n\n\n', ...
     snr(lenImg, IRN_lenImgBlur_03L1), snr(lenImg, MSTV_lenImgBlur_03L1), snr(lenImg, BKSL1_lenImgBlur_03L1), ...
     tirn_03L1n0, tbsk_03L1n0, tbsk_03L1n1);
  str_all = [str_all str];
  else
  str = sprintf('Lena       0.1      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f\n', ...
     snr(lenImg, IRN_lenImgBlur_01L1), NaN, NaN, ...
     tirn_01L1n0, NaN, NaN);
  str_all = [str_all str];
  str = sprintf('           0.3      %5.1f     %5.1f    %5.1f          %5.1f     %5.1f    %5.1f\n\n\n', ...
     snr(lenImg, IRN_lenImgBlur_03L1), NaN, NaN, ...
     tirn_03L1n0, NaN, NaN);
  str_all = [str_all str];

  end % _END_ if(BKS_CODE)

end

end % _END_ if( strcmp(example,'l1deconv') || strcmp(example,'icip09') )

%-----------------------------------------------------------------------------

if( strcmp(example,'l2deconv') || strcmp(example,'icip09') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L2 Deconvolved        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(extreme == 1)

%noise level: 0.1

%-- IRN
lambda  = 0.04;
pars = irntvInputPars('l2tv');

pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 3;
pars.pcgtol_ini   = 1e-4;
%  pars.U0           = lenImgBlur_01L2;

tic;
IRN_lenImgBlur_01L2 = irntv(lenImgBlur_01L2, KC_BSK, lambda, pars);
tirn_01L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_lenImgBlur_01L2) );
  axis image; axis off;
  title(sprintf('Deconvolved Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, IRN_lenImgBlur_01L2), tirn_01L2n0));
end

if(BKS_CODE)
%-- BSK MSTV (Mumford-Shah TV)

Params      = SetParams;
Params.beta = 0.01*max(lenImgBlur(:));

tic;
[MSTV_lenImgBlur_01L2, V_MSTV_001L2] = MainRestoration(lenImgBlur_01L2, ...
                                              kernel_BSK, 'L2', 'MS', Params);
tbsk_01L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( MSTV_lenImgBlur_01L2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, MSTV_lenImgBlur_01L2), tbsk_01L2n0));
end

end % _END_ if(BKS_CODE)
end % _END_ if(extreme == 1)
%-----------------------------------------------------------------------------

%noise level: 0.05

%-- IRN
lambda  = 0.01;
pars = irntvInputPars('l2tv');

pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 3;
%  pars.pcgtol_ini   = 1e-4;
pars.U0           = lenImgBlur_005L2;

tic;
IRN_lenImgBlur_005L2 = irntv(lenImgBlur_005L2, KC_BSK, lambda, pars);
tirn_005L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_lenImgBlur_005L2) );
  axis image; axis off;
  title(sprintf('Deconvolved Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, IRN_lenImgBlur_005L2), tirn_005L2n0));
end

if(BKS_CODE)
%-- BSK

Params      = SetParams;
Params.beta = 0.0025*max(lenImgBlur(:));

tic;
[MSTV_lenImgBlur_005L2, V_MSTV_05L2] = MainRestoration(lenImgBlur_005L2, ...
                                              kernel_BSK, 'L2', 'MS', Params);
tbsk_005L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( MSTV_lenImgBlur_005L2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, MSTV_lenImgBlur_005L2), tbsk_005L2n0));
end

%-- BSK l2-VTV

Params      = SetParams;
Params.beta = 0.0025*max(lenImgBlur(:));

tic;
[BSKL2_lenImgBlur_005L2, V_MSTV_05L2] = MainRestoration(lenImgBlur_005L2, ...
                                              kernel_BSK, 'L2', 'L1', Params);
tbsk_005L2n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( BSKL2_lenImgBlur_005L2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, BSKL2_lenImgBlur_005L2), tbsk_005L2n1));
end

end % _END_ if(BKS_CODE)

%-----------------------------------------------------------------------------

%noise level: 0.01

%-- IRN
lambda  = 0.0005;

pars = irntvInputPars('l2tv');

pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 3;
%  pars.pcgtol_ini   = 1e-4;
pars.U0           = lenImgBlur_001L2;

tic;
IRN_lenImgBlur_001L2 = irntv(lenImgBlur_001L2, KC_BSK, lambda, pars);
tirn_001L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_lenImgBlur_001L2) );
  axis image; axis off;
  title(sprintf('Deconvolved Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, IRN_lenImgBlur_001L2), tirn_001L2n0));
end

if(BKS_CODE)
%-- BSK MSTV (Mumford-Shah TV)

Params      = SetParams;
Params.beta = 0.0001;

tic;
[MSTV_lenImgBlur_001L2, V_MSTV_001L2] = ...
  MainRestoration(Normalize(lenImgBlur_001L2), kernel_BSK, 'L2', 'MS', Params);
tbsk_001L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( MSTV_lenImgBlur_001L2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, MSTV_lenImgBlur_001L2), tbsk_001L2n0));
end

%-- BSK L2-VTV

Params      = SetParams;
Params.beta = 0.0001;

tic;
[BSKL2_lenImgBlur_001L2, V_MSTV_001L2] = ...
  MainRestoration(Normalize(lenImgBlur_001L2), kernel_BSK, 'L2', 'L1', Params);
tbsk_001L2n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( BSKL2_lenImgBlur_001L2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, BSKL2_lenImgBlur_001L2), tbsk_001L2n1));
end

end % _END_ if(BKS_CODE)

%-----------------------------------------------------------------------------

%noise level: sqrt(1e-5)   BSK example

%-- IRN
lambda  = 0.0001;
pars = irntvInputPars('l2tv');

pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 3;
%  pars.pcgtol_ini   = 1e-4;
pars.U0           = lenImgBlur_bskL2;

tic;
IRN_lenImgBlur_bskL2 = irntv(lenImgBlur_bskL2, KC_BSK, lambda, pars);
tirn_bskL2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_lenImgBlur_bskL2) );
  axis image; axis off;
  title(sprintf('Deconvolved Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, IRN_lenImgBlur_bskL2), tirn_bskL2n0));
end

if(BKS_CODE)

%-- BSK MSTV (Mumford-Shah TV)

Params      = SetParams;
Params.beta = 0.00001;

tic;
[MSTV_lenImgBlur_bskL2, V_MSTV_bskL2] = ...
  MainRestoration(Normalize(lenImgBlur_bskL2), kernel_BSK, 'L2', 'MS', Params);
tbsk_bskL2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( MSTV_lenImgBlur_bskL2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, MSTV_lenImgBlur_bskL2), tbsk_bskL2n0));
end

%-- BSK l2-VTV

Params      = SetParams;
Params.beta = 0.00001;

tic;
[BSKL2_lenImgBlur_bskL2, V_MSTV_bskL2] = ...
  MainRestoration(Normalize(lenImgBlur_bskL2), kernel_BSK, 'L2', 'L1', Params);
tbsk_bskL2n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( BSKL2_lenImgBlur_bskL2 );
  axis image; axis off;
  title(sprintf('Deconvolved Image - MSTV. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(lenImg, BSKL2_lenImgBlur_bskL2), tbsk_bskL2n1));
end

end % _END_ if(BKS_CODE)

%-----------------------------------------------------------------------------

if(ICIP_FLAG == 0)

  disp('L2 Deconvolve Vector TV');
  disp(' ');
  disp('                        SNR (db)                    Time (s)');
  disp(' Img      Noise    VTV-IRN  L2-MS(BSK)   BSK-L2          VTV IRN   L2-MS(BKS)  BSK-L2');
  disp(' ');
  if(BKS_CODE)
  disp(sprintf('Lenna      0.05      %5.1f     %5.1f     %5.1f           %5.1f      %5.1f    %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_005L2), snr(lenImg, MSTV_lenImgBlur_005L2), snr(lenImg, BSKL2_lenImgBlur_005L2), ...
     tirn_005L2n0, tbsk_005L2n0, tbsk_005L2n1));
  disp(sprintf('Lenna      0.01      %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_001L2), snr(lenImg, MSTV_lenImgBlur_001L2), snr(lenImg, BSKL2_lenImgBlur_001L2), ...
     tirn_001L2n0, tbsk_001L2n0, tbsk_001L2n1));
  disp(sprintf('Lenna   sqrt(1e-5)   %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_bskL2), snr(lenImg, MSTV_lenImgBlur_bskL2), snr(lenImg, BSKL2_lenImgBlur_bskL2), ...
     tirn_bskL2n0, tbsk_bskL2n0, tbsk_bskL2n1));
  else
  disp(sprintf('Lenna      0.05      %5.1f     %5.1f     %5.1f           %5.1f      %5.1f    %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_005L2), NaN, NaN, ...
     tirn_005L2n0, NaN, NaN));
  disp(sprintf('Lenna      0.01      %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_001L2), NaN, NaN, ...
     tirn_001L2n0, NaN, NaN));
  disp(sprintf('Lenna   sqrt(1e-5)   %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f', ...
     snr(lenImg, IRN_lenImgBlur_bskL2), NaN, NaN, ...
     tirn_bskL2n0, NaN, NaN));
  end % _END_ if(BKS_CODE)

else

  str = sprintf('L2 Deconvolve Vector TV\n\n');
  str_all = [str_all str];
  str = sprintf('                        SNR (db)                    Time (s)\n');
  str_all = [str_all str];
  str = sprintf(' Img      Noise    VTV-IRN  L2-MS(BSK)   BSK-L2          VTV IRN   L2-MS(BKS)  BSK-L2\n\n');
  str_all = [str_all str];

  if(BKS_CODE)
  str = sprintf('Lena       0.05      %5.1f     %5.1f     %5.1f           %5.1f      %5.1f    %5.1f\n', ...
     snr(lenImg, IRN_lenImgBlur_005L2), snr(lenImg, MSTV_lenImgBlur_005L2), snr(lenImg, BSKL2_lenImgBlur_005L2), ...
     tirn_005L2n0, tbsk_005L2n0, tbsk_005L2n1);
  str_all = [str_all str];
  str = sprintf('Lena       0.01      %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f\n', ...
     snr(lenImg, IRN_lenImgBlur_001L2), snr(lenImg, MSTV_lenImgBlur_001L2), snr(lenImg, BSKL2_lenImgBlur_001L2), ...
     tirn_001L2n0, tbsk_001L2n0, tbsk_001L2n1);
  str_all = [str_all str];
  str = sprintf('Lena    sqrt(1e-5)   %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f\n\n\n', ...
     snr(lenImg, IRN_lenImgBlur_bskL2), snr(lenImg, MSTV_lenImgBlur_bskL2), snr(lenImg, BSKL2_lenImgBlur_bskL2), ...
     tirn_bskL2n0, tbsk_bskL2n0, tbsk_bskL2n1);
  str_all = [str_all str];
  else
  str = sprintf('Lena       0.05      %5.1f     %5.1f     %5.1f           %5.1f      %5.1f    %5.1f\n', ...
     snr(lenImg, IRN_lenImgBlur_005L2), NaN, NaN, ...
     tirn_005L2n0, NaN, NaN);
  str_all = [str_all str];
  str = sprintf('Lena       0.01      %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f\n', ...
     snr(lenImg, IRN_lenImgBlur_001L2), NaN, NaN, ...
     tirn_001L2n0, NaN, NaN);
  str_all = [str_all str];
  str = sprintf('Lena    sqrt(1e-5)   %5.1f     %5.1f      %5.1f          %5.1f      %5.1f     %5.1f\n\n\n', ...
     snr(lenImg, IRN_lenImgBlur_bskL2), NaN, NaN, ...
     tirn_bskL2n0, NaN, NaN);
  str_all = [str_all str];
  end % _END_ if(BKS_CODE)

end

end % _END_ if( strcmp(example,'l2deconv') || strcmp(example,'icp09') )

%-----------------------------------------------------------------------------

if( strcmp(example,'l1denoise') || strcmp(example,'icip09') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L1 Denoise        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%noise level: 0.1

%-- IRN
lambda  = 1.1;

pars = irntvInputPars('l1tv');

pars.pcgtol_ini = 1e-4;
pars.epsF       = 1e-2;
pars.epsR       = 1e-4;
pars.loops      = 2;

%-- peppers
tic;
IRN_pepImg_01L1 = irntv(pepImg_01L1, [], lambda, pars);
tirn_01L1n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_pepImg_01L1) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(pepImg, IRN_pepImg_01L1), tirn_01L1n0));
end

%-- mandrill
tic;
IRN_mdrilImg_01L1 = irntv(mdrilImg_01L1, [], lambda, pars);
tirn_01L1n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_mdrilImg_01L1) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(mdrilImg, IRN_mdrilImg_01L1), tirn_01L1n1));
end

%-----------------------------------------------------------------------------

%noise level: 0.3

%-- IRN
lambda  = 1.2;

pars = irntvInputPars('l1tv');

pars.pcgtol_ini = 1e-4;
pars.epsF       = 1e-2;
pars.epsR       = 1e-4;
pars.loops      = 2;

%-- peppers
tic;
IRN_pepImg_03L1 = irntv(pepImg_03L1, [], lambda, pars);
tirn_03L1n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_pepImg_03L1) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(pepImg, IRN_pepImg_03L1), tirn_03L1n0));
end

%-- mandrill

tic;
IRN_mdrilImg_03L1 = irntv(mdrilImg_03L1, [], lambda, pars);
tirn_03L1n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_mdrilImg_03L1) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(mdrilImg, IRN_mdrilImg_03L1), tirn_03L1n1));
end

%-----------------------------------------------------------------------------

if(ICIP_FLAG == 0)

  disp('L1 Denoise Vector TV (Vector IRN algorithm)');
  disp(' ');
  disp('                   SNR (db)          Time (s)');
  disp(' Img      Noise    VTV IRN           VTV IRN');
  disp(' ');
  disp(sprintf('Peppers    0.1       %5.1f            %5.1f  ', ...
      snr(pepImg, IRN_pepImg_01L1), tirn_01L1n0));
  disp(sprintf('           0.3       %5.1f            %5.1f  ', ...
      snr(pepImg, IRN_pepImg_03L1), tirn_03L1n0));
  disp(sprintf('Mandrill   0.1       %5.1f            %5.1f  ', ...
      snr(mdrilImg, IRN_mdrilImg_01L1), tirn_01L1n1));
  disp(sprintf('           0.3       %5.1f            %5.1f  ', ...
      snr(mdrilImg, IRN_mdrilImg_03L1), tirn_03L1n1));

else

  str = sprintf('L1 Denoise Vector TV (Vector IRN algorithm)\n\n');
  str_all = [str_all str];
  str = sprintf('                   SNR (db)          Time (s)\n');
  str_all = [str_all str];
  str = sprintf(' Img      Noise    VTV IRN           VTV IRN\n\n');
  str_all = [str_all str];

  str = sprintf('Peppers    0.1       %5.1f            %5.1f  \n', ...
      snr(pepImg, IRN_pepImg_01L1), tirn_01L1n0);
  str_all = [str_all str];
  str = sprintf('           0.3       %5.1f            %5.1f  \n', ...
      snr(pepImg, IRN_pepImg_03L1), tirn_03L1n0);
  str_all = [str_all str];
  str = sprintf('Mandrill   0.1       %5.1f            %5.1f  \n', ...
      snr(mdrilImg, IRN_mdrilImg_01L1), tirn_01L1n1);
  str_all = [str_all str];
  str = sprintf('           0.3       %5.1f            %5.1f  \n\n\n', ...
      snr(mdrilImg, IRN_mdrilImg_03L1), tirn_03L1n1);
  str_all = [str_all str];

end

end % _END_ if( strcmp(example,'l1denoise') || strcmp(example,'icip09') )

%-----------------------------------------------------------------------------

if( strcmp(example,'l2denoise') || strcmp(example,'icip09') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L2 Denoise        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%noise level: 0.1

%-- IRN
lambda  = 0.3;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini = 1e-4;
%  pars.epsF       = 1e-1;
%  pars.epsR       = 1e-2;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops      = 1;

tic;
IRN_pepImg_01L2 = irntv(pepImg_01L2, [], lambda, pars);
tirn_01L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_pepImg_01L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(pepImg, IRN_pepImg_01L2), tirn_01L2n0));
end

if(BnG_CODE)
%% http://www.mathworks.fr/matlabcentral/fileexchange/16236
%%  Pascal Getreuer based on X. Bresson and T.F. Chan,
%%  "Fast Minimization of the Vectorial Total Variation Norm and Applications
%%  to Color Image Processing", CAM Report 07-25.

lambda  = 0.15;

tic;
BnG_pepImg_01L2 = tvdenoise(pepImg_01L2, 1/lambda);
tirn_01L2n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(BnG_pepImg_01L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - BnG. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(pepImg, BnG_pepImg_01L2), tirn_01L2n1));
end

end % _END_ if(BnG_CODE)

%------------

%-- IRN
lambda  = 0.3;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini = 1e-4;
%  pars.epsF       = 1e-1;
%  pars.epsR       = 1e-2;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops      = 1;

tic;
IRN_mdrilImg_01L2 = irntv(mdrilImg_01L2, [], lambda, pars);
tirn_01L2n2 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_mdrilImg_01L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(mdrilImg, IRN_mdrilImg_01L2), tirn_01L2n2));
end

if(BnG_CODE)
%% http://www.mathworks.fr/matlabcentral/fileexchange/16236
%%  Pascal Getreuer based on X. Bresson and T.F. Chan,
%%  "Fast Minimization of the Vectorial Total Variation Norm and Applications
%%  to Color Image Processing", CAM Report 07-25.

lambda  = 0.15;

tic;
BnG_mdrilImg_01L2 = tvdenoise(mdrilImg_01L2, 1/lambda);
tirn_01L2n3 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(BnG_mdrilImg_01L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - BnG. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(mdrilImg, BnG_mdrilImg_01L2), tirn_01L2n3));
end

end % _END_ if(BnG_CODE)

%-----------------------------------------------------------------------------

%noise level: 0.05

%-- IRN
lambda  = 0.1;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini = 1e-4;
%  pars.epsF       = 1e-1;
%  pars.epsR       = 1e-2;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops      = 1;

tic;
IRN_pepImg_005L2 = irntv(pepImg_005L2, [], lambda, pars);
tirn_005L2n0 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_pepImg_005L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(pepImg, IRN_pepImg_005L2), tirn_005L2n0));
end

if(BnG_CODE)
%% http://www.mathworks.fr/matlabcentral/fileexchange/16236
%%  Pascal Getreuer based on X. Bresson and T.F. Chan,
%%  "Fast Minimization of the Vectorial Total Variation Norm and Applications
%%  to Color Image Processing", CAM Report 07-25.

lambda  = 0.05;

tic;
BnG_pepImg_005L2 = tvdenoise(pepImg_005L2, 1/lambda);
tirn_005L2n1 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(BnG_pepImg_005L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - BnG. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(pepImg, BnG_pepImg_005L2), tirn_005L2n1));
end

end % _END_ if(BnG_CODE)

%------------

%-- IRN
lambda  = 0.1;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini = 1e-4;
%  pars.epsF       = 1e-1;
%  pars.epsR       = 1e-2;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops      = 1;

tic;
IRN_mdrilImg_005L2 = irntv(mdrilImg_005L2, [], lambda, pars);
tirn_005L2n2 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(IRN_mdrilImg_005L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - Vector IRN. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(mdrilImg, IRN_mdrilImg_005L2), tirn_005L2n2));
end

if(BnG_CODE)
%% http://www.mathworks.fr/matlabcentral/fileexchange/16236
%%  Pascal Getreuer based on X. Bresson and T.F. Chan,
%%  "Fast Minimization of the Vectorial Total Variation Norm and Applications
%%  to Color Image Processing", CAM Report 07-25.

lambda  = 0.05;

tic;
BnG_mdrilImg_005L2 = tvdenoise(mdrilImg_005L2, 1/lambda);
tirn_005L2n3 = toc;

if(SHOW_IMGS)
  figure; imagesc( Normalize(BnG_mdrilImg_005L2) );
  axis image; axis off;
  title(sprintf('Denoised Image - BnG. SNR: %4.1fdB.\n Time %4.1f sec', ...
               snr(mdrilImg, BnG_mdrilImg_005L2), tirn_005L2n3));
end

end % _END_ if(BnG_CODE)

%-----------------------------------------------------------------------------

if(ICIP_FLAG == 0)

  disp('L2 Denoise Vector TV (Vector IRN algorithm)');
  disp(' ');
  disp('                         SNR (db)                    Time (s)');
  disp(' Img      Noise   VTV-IRN   bresson-2008-fast     VTV IRN   bresson-2008-fast');
  disp(' ');
  if(BnG_CODE)
  disp(sprintf('Peppers,   0.1      %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(pepImg, IRN_pepImg_01L2), snr(pepImg, BnG_pepImg_01L2), ...
      tirn_01L2n0, tirn_01L2n1));
  disp(sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(pepImg, IRN_pepImg_005L2), snr(pepImg, BnG_pepImg_005L2), ...
      tirn_005L2n0, tirn_005L2n1));
  disp(sprintf('Mandrill,  0.1      %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(mdrilImg, IRN_mdrilImg_01L2), snr(mdrilImg, BnG_mdrilImg_01L2), ...
      tirn_01L2n2, tirn_01L2n3));
  disp(sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(mdrilImg, IRN_mdrilImg_005L2), snr(mdrilImg, BnG_mdrilImg_005L2), ...
      tirn_005L2n2, tirn_005L2n3));
  else
  disp(sprintf('Peppers,   0.1      %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(pepImg, IRN_pepImg_01L2), NaN, ...
      tirn_01L2n0, NaN));
  disp(sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(pepImg, IRN_pepImg_005L2), NaN, ...
      tirn_005L2n0, NaN));
  disp(sprintf('Mandrill,  0.1      %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(mdrilImg, IRN_mdrilImg_01L2), NaN, ...
      tirn_01L2n2, NaN));
  disp(sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f', ...
      snr(mdrilImg, IRN_mdrilImg_005L2), NaN, ...
      tirn_005L2n2, NaN));
  end % _END_ if(BnG_CODE)

else

  str = sprintf('L2 Denoise Vector TV (Vector IRN algorithm)\n\n');
  str_all = [str_all str];
  str = sprintf('                         SNR (db)                    Time (s)\n');
  str_all = [str_all str];
  str = sprintf(' Img      Noise   VTV-IRN   bresson-2008-fast     VTV IRN   bresson-2008-fast\n\n');
  str_all = [str_all str];

  if(BnG_CODE)
  str = sprintf('Peppers,   0.1      %5.1f      %5.1f                %5.1f         %5.1f\n', ...
      snr(pepImg, IRN_pepImg_01L2), snr(pepImg, BnG_pepImg_01L2), ...
      tirn_01L2n0, tirn_01L2n1);
  str_all = [str_all str];
  str = sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f\n', ...
      snr(pepImg, IRN_pepImg_005L2), snr(pepImg, BnG_pepImg_005L2), ...
      tirn_005L2n0, tirn_005L2n1);
  str_all = [str_all str];
  str = sprintf('Mandrill,  0.1      %5.1f      %5.1f                %5.1f         %5.1f\n', ...
      snr(mdrilImg, IRN_mdrilImg_01L2), snr(mdrilImg, BnG_mdrilImg_01L2), ...
      tirn_01L2n2, tirn_01L2n3);
  str_all = [str_all str];
  str = sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f\n\n\n', ...
      snr(mdrilImg, IRN_mdrilImg_005L2), snr(mdrilImg, BnG_mdrilImg_005L2), ...
      tirn_005L2n2, tirn_005L2n3);
  str_all = [str_all str];
  else
  str = sprintf('Peppers,   0.1      %5.1f      %5.1f                %5.1f         %5.1f\n', ...
      snr(pepImg, IRN_pepImg_01L2), NaN, ...
      tirn_01L2n0, NaN);
  str_all = [str_all str];
  str = sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f\n', ...
      snr(pepImg, IRN_pepImg_005L2), NaN, ...
      tirn_005L2n0, NaN);
  str_all = [str_all str];
  str = sprintf('Mandrill,  0.1      %5.1f      %5.1f                %5.1f         %5.1f\n', ...
      snr(mdrilImg, IRN_mdrilImg_01L2), NaN, ...
      tirn_01L2n2, NaN);
  str_all = [str_all str];
  str = sprintf('           0.05     %5.1f      %5.1f                %5.1f         %5.1f\n\n\n', ...
      snr(mdrilImg, IRN_mdrilImg_005L2), NaN, ...
      tirn_005L2n2, NaN);
  str_all = [str_all str];
  end % _END_ if(BnG_CODE)

end

end % _END_ if( strcmp(example,'l2denoise') || strcmp(example,'icp09') )

%-----------------------------------------------------------------------------

if(ICIP_FLAG == 1)

  disp(str_all);
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [HDR]=sseek(HDR,offset,origin)
% SSEEK repositions file position indicator
%  HDR = sseek(HDR,offset,origin)
%
%  offset       number of samples relative to origin
%  origin
%        -1,'bof':   begin of file
%         0,'cof':   current file position
%        +1,'eof':   end of file
%
% See also: SOPEN, SREAD, SWRITE, SCLOSE, SSEEK, SREWIND, STELL, SEOF

%	$Id: sseek.m 2205 2009-10-27 12:18:15Z schloegl $
%	(C) 1997-2005,2006,2007 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

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

if strcmp(origin,'bof')
	origin = -1;
elseif strcmp(origin,'cof')
	origin = 0;
elseif strcmp(origin,'eof')
	origin = 1;
end;

if origin == -1,
        HDR.FILE.POS = offset;
        if strmatch(HDR.TYPE,{'ACQ','BDF','EDF','GDF','EPL','MFER','SCP','native','TMS32','WG1'}),
	elseif HDR.FILE.FID>2,
                POS = HDR.HeadLen+HDR.AS.bpb*offset;
                if POS~=ceil(POS),  % for alpha format
                        fprintf(HDR.FILE.stderr,'Error STELL (alpha): starting position is non-integer\n');
                        return;
                end
                HDR.FILE.status = fseek(HDR.FILE.FID,POS,-1);
        end

elseif origin == 0,
        HDR.FILE.POS = HDR.FILE.POS + offset;
        if strmatch(HDR.TYPE,{'ACQ','BDF','EDF','GDF','EPL','MFER','SCP','native','TMS32','WG1'}),
	elseif HDR.FILE.FID>2,
                POS = HDR.AS.bpb*offset;
                if POS~=ceil(POS),  % for alpha format
                        fprintf(HDR.FILE.stderr,'Error STELL (alpha): starting position is non-integer\n');
                        return;
                end
                HDR.FILE.status = fseek(HDR.FILE.FID,POS,0);
        end

elseif origin == 1,
	if 0, %strmatch(HDR.TYPE,{}),
		HDR.FILE.POS = HDR.NRec+offset;
		HDR.FILE.status = fseek(HDR.FILE.FID,HDR.AS.bpb*offset,1);
	elseif strmatch(HDR.TYPE,{'ACQ','BDF','EDF','GDF','EPL','Sigma'}),
		HDR.FILE.POS = HDR.NRec*HDR.SPR+offset;
                if HDR.FILE.POS < 0,
                        HDR.FILE.POS = 0;
                elseif HDR.FILE.POS > HDR.NRec*HDR.SPR,
                        HDR.FILE.POS = HDR.NRec*HDR.SPR;
                end;
		HDR.FILE.status = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.FILE.POS*HDR.AS.bpb,-1);
	elseif strmatch(HDR.TYPE,{'CTF','Nicolet'}),
		POS = HDR.AS.endpos+offset*HDR.AS.bpb;
		HDR.FILE.status = fseek(HDR.FILE.FID,POS,-1);
		HDR.FILE.POS = (POS-HDR.HeadLen)/HDR.AS.bpb;
	elseif strmatch(HDR.TYPE,{'AINF','BKR','ET-MEG','ISHNE','RG64','MIT','LABVIEW','SMA','BVbinmul','BCI2000'}),
		HDR.FILE.POS = HDR.AS.endpos+offset;
		HDR.FILE.status = fseek(HDR.FILE.FID,HDR.AS.bpb*offset,1);
	elseif strmatch(HDR.TYPE,{'CNT','EEG','AVG','EGI','SND','WAV','AIF','CFWB','DEMG'}),
		HDR.FILE.POS = HDR.AS.endpos+offset;
		HDR.FILE.status = fseek(HDR.FILE.FID,HDR.HeadLen+(HDR.AS.endpos+offset)*HDR.AS.bpb,-1);
        elseif strmatch(HDR.TYPE,{'alpha'}),
                POS = HDR.HeadLen+(HDR.AS.endpos+offset)*HDR.AS.bpb;
                if POS~=ceil(POS),  % for alpha format
                        fprintf(HDR.FILE.stderr,'Error STELL (alpha): starting position is non-integer\n');
                        return;
                end
		HDR.FILE.POS = HDR.AS.endpos+offset;
		HDR.FILE.status = fseek(HDR.FILE.FID,POS,-1);
        elseif strmatch(HDR.TYPE,{'RDF','SIGIF'}),
		HDR.FILE.POS = length(HDR.Block.Pos)+offset;
        elseif strmatch(HDR.TYPE,{'BVascii','BVbinvec','EEProbe-CNT','EEProbe-AVR','FIF','native','MFER','SCP','TMS32','WG1'}),
		HDR.FILE.POS = HDR.AS.endpos+offset;
	else
		fprintf(HDR.FILE.stderr,'Warning SSEEK: format %s not supported.\n',HDR.TYPE);
	end;
else
        fprintf(2,'error SSEEK: 3rd argument "%s" invalid\n',origin);
        return;
end;

author: PatternRecognition

repo: OpenBMI


text:
% collection of m-files to extend the statistics toolbox with the
% NIG distribution
%
% INVGRND       random number generator for the inverse Gaussian
%               distribution
% NIGCDF        cumulative distribution function for NIG distribution
% NIGINV        inverse cdf, computes the quantiles of the NIG distribution
% NIGPAR        parameters for NIG distribution given mean, standard
%               deviation, skewness and kurtosis
% NIGPDF        probability density function for NIG distribution
% NIGRND        random number generator for the normal inverse Gaussian
%               distribution
% NIGSTAT       mean, standard deviation, skewness and kurtosis for a
%               given parameter set

author: Sable

repo: mcbench-benchmarks


text:
function proposals = genProposalsByOrientation(img, region, regionPerim, estimated_orientation, region_comp_infos)
    global param;

    if(isnan(estimated_orientation))
        proposals = zeros(0, 9);
        return
    end

    [comp_cluster_ind, ~] = compCluster(img, region_comp_infos, ...
                                                    estimated_orientation, ...
                                                    param.minCompHeightSimilarity, ...
                                                    param.maxCompOrientationDiff, ...
                                                    param.minIoUDiff);

    clusterIdx = unique(comp_cluster_ind);
    clusterCount = length(clusterIdx);
    proposals = zeros(clusterCount, 9);
    confident_idx = false(clusterCount, 1);

    for c = 1 : clusterCount
        clusterCompInfos = region_comp_infos(comp_cluster_ind == clusterIdx(c));
        if(length(clusterCompInfos) > 0)
            confident_idx(c) = true;
            proposals(c, 1 : 8) = getProposal(img, region, regionPerim, clusterCompInfos, estimated_orientation);
            proposals(c, 9) = estimated_orientation;
        end
    end

    proposals = proposals(confident_idx, :);

    %% show proposals
    if true && param.debug
        imshow(img);
        hold on;
        for i = 1 : clusterCount
            x_arr = [proposals(i, 1 : 2 : 8), proposals(i, 1)];
            y_arr = [proposals(i, 2 : 2 : 8), proposals(i, 2)];

            plot(x_arr, y_arr, 'color', rand(3,1));
        end
        hold off;
    end

    %% show clusters
    if false && param.debug
        debug_bbox_gathered = zeros(0, 4);
        color_gathered = zeros(0, 3);
        for c = 1: clusterCount
            clusterCompInfos = region_comp_infos(comp_cluster_ind == clusterIdx(c));
            debug_bbox = zeros(length(clusterCompInfos), 4);
            color = repmat(rand(1,3), length(clusterCompInfos), 1);

            for d_i = 1 : length(clusterCompInfos)
                debug_bbox(d_i,:) = clusterCompInfos{d_i}.box;
            end

            debug_bbox_gathered = cat(1, debug_bbox_gathered, debug_bbox);
            color_gathered = cat(1, color_gathered, color);
        end
        show_bbox(img, debug_bbox_gathered, color_gathered);
    end
end

function [comp_cluster_ind, clusterCount] = compCluster(img, region_comp_infos, estimated_orientation, minCompHeightSimilarity, maxCompOrientationDiff, minIoUDiff)
    global param;

    compCount = length(region_comp_infos);

    if(compCount == 0)
        comp_cluster_ind = zeros(0, 1);
        clusterCount = 0;
        return;
    end

    if(compCount == 1)
        comp_cluster_ind = ones(1,1);
        clusterCount = 1;
        return;
    end

    comp_cluster_ind = zeros(compCount, 1, 'uint16');
    clusterCount = 0;

    box4d = zeros(compCount, 4);
    heights = zeros(compCount, 1);

    boxCenter = zeros(compCount, 2);
    boxCenter_U = zeros(compCount, 2);
    boxCenter_D = zeros(compCount, 2);
    boxCenter_L = zeros(compCount, 2);
    boxCenter_R = zeros(compCount, 2);

    for i = 1 : compCount
        box4d(i, :) = region_comp_infos{i}.box;
        [rectX, rectY, ~, ~, sideLenght] = computeMergedBox(box4d(i, :));
        [~, heights(i)] = findWidthAndHeight(rectX, rectY, sideLenght, estimated_orientation);

        boxCenter(i, 2) = round(box4d(i, 1) + 0.5*box4d(i, 3));
        boxCenter(i, 1) = round(box4d(i, 2) + 0.5*box4d(i, 4));

        boxCenter_U(i, 2) = round(box4d(i, 1) + 0.5*box4d(i, 3));
        boxCenter_U(i, 1) = box4d(i, 2);

        boxCenter_D(i, 2) = round(box4d(i, 1) + 0.5*box4d(i, 3));
        boxCenter_D(i, 1) = box4d(i,2) + box4d(i,4) - 1;

        boxCenter_L(i, 2) = box4d(i,1);
        boxCenter_L(i, 1) = round(box4d(i, 2) + 0.5*box4d(i, 4));

        boxCenter_R(i, 2) = box4d(i,1) + box4d(i,3) - 1;
        boxCenter_R(i, 1) = round(box4d(i, 2) + 0.5*box4d(i, 4));
    end

    heightSimilarity = zeros(compCount);
    IoUDiff = zeros(compCount);
    orientationDiffOfBox = zeros(compCount);
    dists = zeros(compCount);

    for i = 1 : compCount

        %% Compute height similarity
        heightSimilarity(i,:) = min( ...
                                    min(box4d(i, 4), box4d(:, 4)) ./ max(box4d(i, 4), box4d(:, 4)), ...
                                    min(box4d(i, 3), box4d(:, 3)) ./ max(box4d(i, 3), box4d(:, 3)));
       %% Compute IoU diff
        intArea = rectint(box4d(i,:), box4d);
        unionArea = box4d(i,3) * box4d(i,4) + box4d(:, 3).* box4d(:, 4);
        IoUDiff(i, :) = intArea ./ unionArea';

        %% Compute orientation diff of box
        orientationDiff_U = computeOrientationDiff(boxCenter_U, boxCenter_U(i,:), estimated_orientation);
        orientationDiff_D = computeOrientationDiff(boxCenter_D, boxCenter_D(i,:), estimated_orientation);
        orientationDiff_L = computeOrientationDiff(boxCenter_L, boxCenter_L(i,:), estimated_orientation);
        orientationDiff_R = computeOrientationDiff(boxCenter_R, boxCenter_R(i,:), estimated_orientation);
        orientationDiffOfBox(i, :) = min(min(orientationDiff_U, orientationDiff_D), ...
                                    min(orientationDiff_L, orientationDiff_R));

        %% Compute dist
        dists(i, :) = ((boxCenter(:, 1) - boxCenter(i, 1)).^2 + (boxCenter(:, 2) - boxCenter(i,2)).^2).^0.5;
    end

    if true && param.debug
        imshow(img);
        for i = 1 : compCount
            rectangle('position', box4d(i, :), 'edgecolor','y');
        end
    end

    for i = 1 : compCount
        if comp_cluster_ind(i) == 0
            clusterCount = clusterCount + 1;
            current_cluster_comp_ind = false(compCount, 1);
            current_cluster_comp_ind(i) = 1;

            isUpdate = true;
            while isUpdate
                isUpdate = false;

                %% Compute boxCenter of merged bbox
                mergedBoxCenter_U = [mean(box4d(current_cluster_comp_ind, 2)), ...
                                        mean(box4d(current_cluster_comp_ind, 1) + 0.5*box4d(current_cluster_comp_ind, 3))];

                mergedBoxCenter_D = [mean(box4d(current_cluster_comp_ind, 2) + box4d(current_cluster_comp_ind, 4) - 1), ...
                                        mean(box4d(current_cluster_comp_ind, 1) + 0.5*box4d(current_cluster_comp_ind, 3))];

                mergedBoxCenter_L = [mean(box4d(current_cluster_comp_ind, 2) + 0.5*box4d(current_cluster_comp_ind, 4)), ...
                                        mean(box4d(current_cluster_comp_ind, 1))];

                mergedBoxCenter_R = [mean(box4d(current_cluster_comp_ind, 2) + 0.5*box4d(current_cluster_comp_ind, 4)), ...
                                        mean(box4d(current_cluster_comp_ind, 1) + box4d(current_cluster_comp_ind, 3) - 1)];

                %% Compute orientation diff
                orientationDiff_U = computeOrientationDiff(boxCenter_U, mergedBoxCenter_U, estimated_orientation);
                orientationDiff_D = computeOrientationDiff(boxCenter_D, mergedBoxCenter_D, estimated_orientation);
                orientationDiff_L = computeOrientationDiff(boxCenter_L, mergedBoxCenter_L, estimated_orientation);
                orientationDiff_R = computeOrientationDiff(boxCenter_R, mergedBoxCenter_R, estimated_orientation);

                orientationDiff = min(min(orientationDiff_U, orientationDiff_D), min(orientationDiff_L, orientationDiff_R))';

                heightSimilarity_rule = heightSimilarity(current_cluster_comp_ind,:) > minCompHeightSimilarity;
                orientationDiff_rule = orientationDiff < maxCompOrientationDiff;
                IoUDiff_rule = IoUDiff(current_cluster_comp_ind, :) > minIoUDiff;
                orientationDiffOfBox_rule = heightSimilarity_rule ...
                                            & (orientationDiffOfBox(current_cluster_comp_ind, :) < maxCompOrientationDiff) ...
                                            & (dists(current_cluster_comp_ind, :) < param.maxDistRatio * mean(box4d(current_cluster_comp_ind, 4)));

                spatial_rule = ((sum(heightSimilarity_rule, 1) > 0) & sum(orientationDiffOfBox_rule, 1) > 0) ...
                                | (sum(IoUDiff_rule, 1) > 0) ...
                                | (orientationDiff_rule & (sum(heightSimilarity_rule, 1) > 0));
                %spatial_rule = sum((heightSimilarity_rule & orientationDiff_rule) | IoUDiff_rule, 1) > 0;

                satisfiedComp = spatial_rule' ...
                            & (~current_cluster_comp_ind);

                satisfiedIdx = find(satisfiedComp);
                if(~isempty(satisfiedIdx))
                    isUpdate = true;
                    current_cluster_comp_ind(satisfiedIdx) = 1;
%                     x1_tmp = min(min(box4d(satisfiedIdx, 1)), ...
%                                  mergedBox(1));
%                     x2_tmp = max(max(box4d(satisfiedIdx, 1) + box4d(satisfiedIdx, 3) - 1), ...
%                                  mergedBox(1) + mergedBox(3) - 1);
%                     y1_tmp = min(min(box4d(satisfiedIdx, 2)), ...
%                                  mergedBox(2));
%                     y2_tmp = max(max(box4d(satisfiedIdx, 2) + box4d(satisfiedIdx, 4) - 1), ...
%                                  mergedBox(2) + mergedBox(4) - 1);
%                     mergedBox(1) = x1_tmp;
%                     mergedBox(2) = y1_tmp;
%                     mergedBox(3) = x2_tmp - x1_tmp + 1;
%                     mergedBox(4) = y2_tmp - y1_tmp + 1;

                    if true &&  param.debug
                        clusterCompInfos = region_comp_infos(current_cluster_comp_ind);
                        color = zeros(length(clusterCompInfos), 3);
                        color(:, 2) = 1;
                        color(1, :) = [1, 0, 0];
                        imshow(img);
                        for d_i = 1 : length(clusterCompInfos)
                           rectangle('position', clusterCompInfos{d_i}.box, 'edgecolor', color(d_i,:));
                        end
                    end
                end
            end
            comp_cluster_ind(current_cluster_comp_ind) = clusterCount;
        end
    end
end

author: stupidZZ


text:
function M = spheresymmetricfactory(n)
% Returns a manifold struct to optimize over unit-norm symmetric matrices.
%
% function M = spheresymmetricfactory(n)
%
% Manifold of n-by-n real symmetric matrices of unit Frobenius norm.
% The metric is such that the sphere is a Riemannian submanifold of the
% space of nxn symmetric matrices with the usual trace inner product, i.e.,
% the usual metric <A, B> = trace(A'*B).
%
% See also: spherefactory obliquefactory spherecomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 17, 2015.
% Contributors:
% Change log:

    M.name = @() sprintf('Sphere of symmetric matrices of size %d', n);

    M.dim = @() n*(n+1)/2 - 1;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) real(acos(x(:).'*y(:)));

    M.typicaldist = @() pi;

    M.proj = @proj;
    function xdot = proj(x, d)
        d = (d+d.')/2;
        xdot = d - x*(x(:).'*d(:));
    end

    M.tangent = @proj;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
	M.egrad2rgrad = @proj;

	M.ehess2rhess = @ehess2rhess;
	function rhess = ehess2rhess(x, egrad, ehess, u)
        % these are not explicitly required, given the use.
        % egrad = (egrad + egrad.')/2;
        % ehess = (ehess + ehess.')/2;
        rhess = proj(x, ehess) - (x(:)'*egrad(:))*u;
	end

    M.exp = @exponential;

    M.retr = @retraction;

    M.log = @logarithm;
    function v = logarithm(x1, x2)
        v = proj(x1, x2 - x1);
        di = M.dist(x1, x2);
        % If the two points are "far apart", correct the norm.
        if di > 1e-6
            nv = norm(v, 'fro');
            v = v * (di / nv);
        end
    end

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() random(n);

    M.randvec = @(x) randomvec(n, x);

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(n);

    M.transp = @(x1, x2, d) proj(x2, d);

    M.pairmean = @pairmean;
    function y = pairmean(x1, x2)
        y = x1+x2;
        y = y / norm(y, 'fro');
    end

    % TODO : check isometry and fix.
    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
    M.vecmatareisometries = @() false;

end

% Exponential on the sphere
function y = exponential(x, d, t)

    if nargin == 2
        t = 1;
    end

    td = t*d;

    nrm_td = norm(td, 'fro');

    if nrm_td > 4.5e-8
        y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
    else
        % If the step is too small to accurately evaluate sin(x)/x,
        % then sin(x)/x is almost indistinguishable from 1.
        y = x + td;
        y = y / norm(y, 'fro');
    end

end

% Retraction on the sphere
function y = retraction(x, d, t)

    if nargin == 2
        t = 1;
    end

    y = x + t*d;
    y = y / norm(y, 'fro');

end

% Uniform random sampling on the sphere.
function x = random(n)

    x = randn(n);
    x = (x + x.')/2;
    x = x/norm(x, 'fro');

end

% Random normalized tangent vector at x.
function d = randomvec(n, x)

    d = randn(n);
    d = (d + d.')/2;
    d = d - x*(x(:).'*d(:));
    d = d / norm(d, 'fro');

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function result = op_selective_search_boxes(first_el, last_el, imdb)
% result = op_selective_search_boxes(first_el, last_el, imdb)
%
% Loop slice operation for distributed computation of selective search
% boxes.
%
% This function depends on simple-cluster-lib, which is specific to
% the Berkeley cluster and not useful to the general public (and
% hence not available). This file exists because it's convenient for
% me to keep it in the repository.
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

fast_mode = true;
mean_num = 0;
mean_time = 0;

result = cell(last_el-first_el+1, 1);
for i = first_el:last_el
  i_ = i-first_el+1;

  fprintf('%d/%d (%s) ...', i, last_el, imdb.image_ids{i});
  im = imread(imdb.image_at(i));
  th = tic();
  result{i_} = selective_search_boxes(im, fast_mode);
  t = toc(th);

  mean_num = (mean_num * (i_-1) + size(result{i_}, 1))/i_;
  mean_time = (mean_time * (i_-1) + t)/i_;
  fprintf('%.2fs...%d boxes (means: %.2fs %.1f boxes)\n', t, ...
      size(result{i_}, 1), mean_time, mean_num);
end

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function value = r8_lgmc ( x )

%*****************************************************************************80
%
%% R8_LGMC evaluates the log gamma correction factor for an R8 argument.
%
%  Discussion:
%
%    For 10 <= X, compute the log gamma correction factor so that
%
%      log ( gamma ( x ) ) = log ( sqrt ( 2 * pi ) )
%                          + ( x - 0.5 ) * log ( x ) - x
%                          + r8_lgmc ( x )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the correction factor.
%
  persistent algmcs
  persistent nalgm
  persistent xbig
  persistent xmax

  if ( isempty ( nalgm ) )

    algmcs = [ ...
      +0.1666389480451863247205729650822, ...
      -0.1384948176067563840732986059135E-04, ...
      +0.9810825646924729426157171547487E-08, ...
      -0.1809129475572494194263306266719E-10, ...
      +0.6221098041892605227126015543416E-13, ...
      -0.3399615005417721944303330599666E-15, ...
      +0.2683181998482698748957538846666E-17, ...
      -0.2868042435334643284144622399999E-19, ...
      +0.3962837061046434803679306666666E-21, ...
      -0.6831888753985766870111999999999E-23, ...
      +0.1429227355942498147573333333333E-24, ...
      -0.3547598158101070547199999999999E-26, ...
      +0.1025680058010470912000000000000E-27, ...
      -0.3401102254316748799999999999999E-29, ...
      +0.1276642195630062933333333333333E-30 ]';

    nalgm = r8_inits ( algmcs, 15, r8_mach ( 3 ) );
    xbig = 1.0 / sqrt ( r8_mach ( 3 ) );
    xmax = exp ( min ( log ( r8_mach ( 2 ) / 12.0 ), ...
      - log ( 12.0 * r8_mach ( 1 ) ) ) );
  end

  if ( x < 10.0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_LGMC - Fatal error!\n' );
    fprintf ( 1, '  X must be at least 10.\n' );
    error ( 'R8_LGMC - Fatal error!' )

  elseif ( x < xbig )

    value = r8_csevl ( 2.0 * ( 10.0 / x ) ...
      * ( 10.0 / x ) - 1.0, algmcs, nalgm ) / x;

  elseif ( x < xmax )

    value = 1.0 / ( 12.0 * x );

  else

    value = 0.0;

  end

  return
end


text:
function result = cross_entropy(head_embedding, tail_embedding, head, tail, weights, a, b, same_embedding)
%CROSS_ENTROPY Given a distance for each 1-simplex in low-dimensional space
% and the original weights of the 1-simplices in high-dimensional space,
% compute the approximation to the cross-entropy between the two simplicial
% complexes. This calculation uses the modified smooth formula Phi for
% low-dimensional weight that is used in the stochastic gradient descent.
%
% result = cross_entropy(head_embedding, tail_embedding, head, tail, weights, a, b, same_embedding)
%
% Parameters
% ----------
% dists: array of size (n_1_simplices, 1)
%     The current distance between the two endpoints of the 1-simplex in
%     low-dimensional Euclidean space.
%
% weights: array of size (n_1_simplices, 1)
%     The original weights assigned to the 1-simplices in high-dimensional
%     space.
%
% a: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% b: double
%     Parameter of differentiable approximation of right adjoint functor.
%
% Returns
% -------
% result: double
%     The total approximated cross entropy between the two simplicial complexes.
%
% See also: NEG_SAMPLING_OBJECTIVE
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
    n1 = size(head_embedding, 1);
    n2 = size(tail_embedding, 1);

    if n1*n2 > 1e8
        error('HALTED: MATLAB usually freezes for embeddings this large.');
    end

    full_dists = pdist2(head_embedding, tail_embedding);
    full_weights = full(sparse(head, tail, weights, n1, n2));
    if same_embedding
        full_weights = full_weights + eye(n1);
    end
    Phi = ones(size(full_weights))./(1 + a*(full_dists.^(2*b)));
    fw0 = full_weights == 0;
    fw1 = full_weights == 1;
    other = ~fw0 & ~fw1;
    Phi_summands = zeros(size(full_weights));
    Phi_summands(fw0) = log(1-Phi(fw0));
    Phi_summands(fw1) = log(Phi(fw1));
    Phi_summands(other) = full_weights(other).*log(Phi(other)) + (1-full_weights(other)).*log(1-Phi(other));

    result = -sum(sum(Phi_summands));
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function eventLog = ma_executeSetState(setStateEvent, initialState, eventNum, celBodyData)
%ma_executeSetState Summary of this function goes here
%   Detailed explanation goes here

    event = setStateEvent;

    if(strcmpi(event.subType,'setState'))
        eventLog = [event.epoch, ...
                    event.rVect(1), event.rVect(2), event.rVect(3), ...
                    event.vVect(1), event.vVect(2), event.vVect(3), ...
                    event.centralBody.id, ...
                    event.dryMass, event.fuelOxMass, event.monoMass, event.xenonMass, ...
                    eventNum];
    elseif(strcmpi(event.subType,'estLaunch'))
        eventLog = ma_executeLaunch(initialState, event.launch, celBodyData);

        eventLog(:,8) = event.centralBody.id;
        eventLog(:,9) = event.dryMass;
        eventLog(:,10) = event.fuelOxMass;
        eventLog(:,11) = event.monoMass;
        eventLog(:,12) = event.xenonMass;
        eventLog(:,13) = eventNum;

%         if(isempty(initialState))
%             eventLog(:,8) = event.centralBody.id;
%             eventLog(:,9) = event.dryMass;
%             eventLog(:,10) = event.fuelOxMass;
%             eventLog(:,11) = event.monoMass;
%             eventLog(:,12) = event.xenonMass;
%             eventLog(:,13) = eventNum;
%         else
%             eventLog(:,8) = initialState(8);
%             eventLog(:,9) = initialState(9);
%             eventLog(:,10) = initialState(10);
%             eventLog(:,11) = initialState(11);
%             eventLog(:,12) = initialState(12);
%             eventLog(:,13) = eventNum;
%         end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%This folder contains the so-called ACQUIRE functions for bbci_apply.
%These functions acquire small bolcks of signals (and maybe event markers
%from a specific measurement device, and provide it for online processing.
%If no new data is available, they return an empty structure.
%
%These functions have the following format:
%
% --- --- --- ---
%BBCI_ACQUIRE_XYZ - Online data acquisition from device XYZ
%
%Synopsis:
%  STATE= bbci_acquire_XYZ('init', <PARAM>)
%  [CNTX, MRKTIME, MRKDESC, STATE]= bbci_acquire_XYZ(STATE)
%  bbci_acquire_XYZ('close')
%  bbci_acquire_XYZ('close', STATE)
%
%Arguments:
%  PARAM - Optional arguments, specific to XYZ.
%
%Output:
%  STATE - Structure characterizing the incoming signals; fields:
%     'fs', 'clab', and intern stuff
%  CNTX - 'acquired' signals [Time x Channels]
%  The following variables hold the markers that have been 'acquired' within
%  the current block (if any).
%  MRKTIME - DOUBLE: [1 nMarkers] position [msec] within data block.
%      A marker occurrence within the first sample would give
%      MARTIME= 1/STATE.fs.
%  MRKDESC - CELL {1 nMarkers} descriptors like 'S 52'
% --- --- --- ---
%
%List of ACQUIRE functions (prefix bbci_acquire_ is left out)
% bv:      Acquire data from BV Recorder (option 'Remote Data Access' must
%          be enabled in the BV Recorder settings!)
% nirx:    Acquire data from a NIRx system
% offline: Simulate online acquisition by returning small chunks of signals
%          from an initially given data file.
% randomSignals: Generate random signals
% lsl:     Acquire data from Lab Streaming Layer(LSL). LSL itself can hold
%           data from most common EEG and other recording devices.

author: bbci