text:
clear; clc; close all;

addpath('../libs/exportFig');
addpath('../libs/fun4MeanShift');
addpath('../libs/layerExt');
addpath('../libs/myFunctions/');
path_to_matconvnet = '../libs/matconvnet-1.0-beta23_modifiedDagnn';
run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));

% set GPU
gpuId = 3;
gpuDevice(gpuId);
%% read matconvnet model
load('imdb_complete_on_server.mat');
imdb.meta.meanvalue = reshape(imdb.meta.meanvalue,[1 1 3]);
imdb.path_to_image = '/home/skong2/dataset/';
imdb.path_to_annotation = '/home/skong2/dataset/';

path_to_model = 'main001_ft4COCO_baseline_v1';
modelName = 'pairMMAbsReg_net-epoch-1.mat';

netbasemodel = load( fullfile('./exp', path_to_model, modelName) );
netbasemodel = netbasemodel.net;
netbasemodel = dagnn.DagNN.loadobj(netbasemodel);
netbasemodel.move('gpu') ;
netbasemodel.mode = 'normal'; % {test, normal} !! normal --> per-instance normalization
netbasemodel.conserveMemory = 1;

netbasemodel.removeLayer('obj_MM');
netbasemodel.removeLayer('obj_reg');
netbasemodel.removeLayer('output_cosSim');
netbasemodel.removeLayer('output_l2norm');

sName = 'res9_conv';
res9_conv = netbasemodel.layers(netbasemodel.getLayerIndex(sName)).outputIndexes;
netbasemodel.vars(res9_conv).precious = 1;

lName = 'output_l2norm';
netbasemodel.addLayer(lName, L2normalization(), sName, lName) ;
sName = lName;
outputIdx_l2norm = netbasemodel.layers(netbasemodel.getLayerIndex('output_l2norm')).outputIndexes;
netbasemodel.vars(outputIdx_l2norm).precious = 1;

loopNum = 16;
predInstanceMaskMat_loop0 = cell(1,3);
predInstanceMaskMat_loops = cell(loopNum,3);

% alpha = 0.5;
% GaussianBandwidth = 1-alpha;
GaussianBandwidth = 0.2;

sName_l2norm = lName;
for loopIdx = 1:loopNum
    [netbasemodel, sName, sName_l2norm] = addOneLoop_forMeanShiftGrouping(netbasemodel, sName_l2norm, loopIdx, GaussianBandwidth);

    keepLayerName = sprintf('loop%d_meanshift_S_is_XX', loopIdx);
    netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.analyzeGradient = false;

    keepLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
    netbasemodel.vars(netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).outputIndexes).precious = 1;

    rmLayerName = sprintf('loop%d_meanshift_cosSim', loopIdx);
    netbasemodel.removeLayer(rmLayerName); % remove layer
end
%% visualize single image
idxList = 1:length(imdb.val_image);
for idx = idxList
    %% read image and annotation
    cur_path_to_image = fullfile(imdb.path_to_image, imdb.val_image{idx});
    cur_path_to_annot = fullfile(imdb.path_to_annot, imdb.val_annot{idx});
    cur_path_to_image = strrep(cur_path_to_image,'\','/');
    cur_path_to_annot = strrep(cur_path_to_annot,'\','/');

    imgOrg = imread(cur_path_to_image);
    cur_segMap = load(cur_path_to_annot);
    cur_instMap = cur_segMap.gtMat.instMap;
    cur_segMap = cur_segMap.gtMat.segMap;
    gtOrg = cur_instMap;

    sz = size(gtOrg);
    reSZ = round(sz/8)*8;
    img = imresize(imgOrg, reSZ);

    instMap_div8 = imresize(cur_instMap, reSZ/8, 'nearest');
    img = bsxfun(@minus, single(img), imdb.meta.meanvalue) ;
    inputs = {'data', gpuArray(single(img))};
    netbasemodel.eval(inputs) ;
    %% gather result
    keepLayerName = 'output_l2norm';
    keepLayerName = netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).outputIndexes;
    feaMap = gather(netbasemodel.vars(keepLayerName).value);

    feaMapSize = size(feaMap);
    randRGB4Pred1 = reshape(feaMap, [], size(feaMap,3));
    rng(777); randProj = randn(feaMapSize(3), 3);
    randRGB4Pred1 = randRGB4Pred1*randProj;
    randRGB4Pred1 = reshape( randRGB4Pred1, [feaMapSize(1), feaMapSize(2), 3] );
    randRGB4Pred1 = rescaleFeaMap(randRGB4Pred1);

    rng(77); randProj = randn(feaMapSize(3), 3);
    randRGB4Pred2 = reshape(feaMap, [], size(feaMap,3));
    randRGB4Pred2 = randRGB4Pred2*randProj;
    randRGB4Pred2 = reshape( randRGB4Pred2, [feaMapSize(1), feaMapSize(2), 3] );
    randRGB4Pred2 = rescaleFeaMap(randRGB4Pred2);

    rng(7); randProj = randn(feaMapSize(3), 3);
    randRGB4Pred3 = reshape(feaMap, [], size(feaMap,3));
    randRGB4Pred3 = randRGB4Pred3*randProj;
    randRGB4Pred3 = reshape( randRGB4Pred3, [feaMapSize(1), feaMapSize(2), 3] );
    randRGB4Pred3 = rescaleFeaMap(randRGB4Pred3);

    randRGB4Pred1 = imresize(randRGB4Pred1, [size(imgOrg,1),size(imgOrg,2)]);
    randRGB4Pred2 = imresize(randRGB4Pred2, [size(imgOrg,1),size(imgOrg,2)]);
    randRGB4Pred3 = imresize(randRGB4Pred3, [size(imgOrg,1),size(imgOrg,2)]);

    %     imwrite(uint8(255*randRGB4Pred1), sprintf('%s/valId%04d/randRGB4Pred1.png', ...
    %         prefix, idx));
    %     imwrite(uint8(255*randRGB4Pred2), sprintf('%s/valId%04d/randRGB4Pred2.png', ...
    %         prefix, idx));
    %     imwrite(uint8(255*randRGB4Pred3), sprintf('%s/valId%04d/randRGB4Pred3.png', ...
    %         prefix, idx));
    %% after some loops of mean shift
    keepLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopNum);
    keepLayerName = netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).outputIndexes;
    feaMap = gather(netbasemodel.vars(keepLayerName).value);

    feaMapSize = size(feaMap);
    randRGB4Pred1_afterMShift = reshape(feaMap, [], size(feaMap,3));
    rng(777); randProj = randn(feaMapSize(3), 3);
    randRGB4Pred1_afterMShift = randRGB4Pred1_afterMShift*randProj;
    randRGB4Pred1_afterMShift = reshape( randRGB4Pred1_afterMShift, [feaMapSize(1), feaMapSize(2), 3] );
    randRGB4Pred1_afterMShift = rescaleFeaMap(randRGB4Pred1_afterMShift);

    rng(77); randProj = randn(feaMapSize(3), 3);
    randRGB4Pred2_afterMShift = reshape(feaMap, [], size(feaMap,3));
    randRGB4Pred2_afterMShift = randRGB4Pred2_afterMShift*randProj;
    randRGB4Pred2_afterMShift = reshape( randRGB4Pred2_afterMShift, [feaMapSize(1), feaMapSize(2), 3] );
    randRGB4Pred2_afterMShift = rescaleFeaMap(randRGB4Pred2_afterMShift);

    rng(7); randProj = randn(feaMapSize(3), 3);
    randRGB4Pred3_afterMShift = reshape(feaMap, [], size(feaMap,3));
    randRGB4Pred3_afterMShift = randRGB4Pred3_afterMShift*randProj;
    randRGB4Pred3_afterMShift = reshape( randRGB4Pred3_afterMShift, [feaMapSize(1), feaMapSize(2), 3] );
    randRGB4Pred3_afterMShift = rescaleFeaMap(randRGB4Pred3_afterMShift);

    randRGB4Pred1_afterMShift = imresize(randRGB4Pred1_afterMShift, [size(imgOrg,1),size(imgOrg,2)]);
    randRGB4Pred2_afterMShift = imresize(randRGB4Pred2_afterMShift, [size(imgOrg,1),size(imgOrg,2)]);
    randRGB4Pred3_afterMShift = imresize(randRGB4Pred3_afterMShift, [size(imgOrg,1),size(imgOrg,2)]);

    figHandler = figure(1);
    set(figHandler, 'Position', [100 100 1500 800]) % [1 1 width height]
    subplot(3,3,1); imshow(imgOrg); title('imgOrg');
    subplot(3,3,2); imagesc(cur_segMap); axis off image; title('segMap');
    subplot(3,3,3); imagesc(cur_instMap); axis off image; title('instMap');

    subplot(3,3,4); imagesc(randRGB4Pred1); axis off image; title('randRGB4Pred1');
    subplot(3,3,5); imagesc(randRGB4Pred2); axis off image; title('randRGB4Pred2');
    subplot(3,3,6); imagesc(randRGB4Pred3); axis off image; title('randRGB4Pred3');
    subplot(3,3,7); imagesc(randRGB4Pred1_afterMShift); axis off image; title('randRGB4Pred1_afterMShift','Interpreter','none');
    subplot(3,3,8); imagesc(randRGB4Pred2_afterMShift); axis off image; title('randRGB4Pred2_afterMShift','Interpreter','none');
    subplot(3,3,9); imagesc(randRGB4Pred3_afterMShift); axis off image; title('randRGB4Pred3_afterMShift','Interpreter','none');
    drawnow;
end
%% leaving blank
%{
%}

author: aimerykong


text:
% Chapter 8 - Planar Systems.
% Program_8c - Phase Portrait of a Nonlinear System (Fig. 8.12).
% Copyright Birkhauser 2013. Stephen Lynch.

% Phase portrait of a nonlinear system of ODE's.
% IMPORTANT - Program_8a is vectorfield.m.
clear
% sys=inline('[x(2);x(1)*(1-(x(1))^2)+x(2)]','t', 'x');
sys = @(t,x) [x(2);x(1)*(1-(x(1))^2)+x(2)];
vectorfield(sys,-3:.5:3,-3:.5:3);
     hold on
     sep=1;
     for x0=-3:sep:3
         for y0=-3:sep:3
            [ts,xs] = ode45(sys,[0 6],[x0 y0]);
            plot(xs(:,1),xs(:,2))
         end
     end
     for x0=-3:sep:3
         for y0=-3:sep:3
            [ts,xs] = ode45(sys,[0 -6],[x0 y0]);
            plot(xs(:,1),xs(:,2))
         end
     end
     hold off
axis([-3 3 -3 3])
fsize=15;
set(gca,'XTick',-3:1:3,'FontSize',fsize)
set(gca,'YTick',-3:1:3,'FontSize',fsize)
xlabel('x(t)','FontSize',fsize)
ylabel('y(t)','FontSize',fsize)
hold off

% End of Program_8c.

author: Sable

repo: mcbench-benchmarks


text:
function [normal,edgeLength,unitNormal] = edgenormal(node,edge)
% compute normal vector of a triangular face

edgeVec = node(edge(:,2),:)-node(edge(:,1),:);
normal = [edgeVec(:,2), -edgeVec(:,1)];
if nargout > 1
    edgeLength = sqrt(sum(normal.^2,2));
    if nargout > 2
        unitNormal = normal./repmat(edgeLength,[1,2]);
    end
end

repo: ifem


text:
function obj = updateHistogram(obj,histIndex)

% x:...[DONE]
% y:...[DONE]
% histnorm:...[DONE]
% name:...[DONE]
% autobinx:...[DONE]
% nbinsx:...[DONE]
% xbins:...[DONE]
% autobiny:...[DONE]
% nbinsy:...[DONE]
% ybins:...[DONE]
% text:...[NOT SUPPORTED IN MATLAB]
% error_y:...[HANDLED BY ERRORBARSERIES]
% error_x:...[HANDLED BY ERRORBARSERIES]
% opacity: --- [TODO]
% xaxis:...[DONE]
% yaxis:...[DONE]
% showlegend:...[DONE]
% stream:...[HANDLED BY PLOTLYSTREAM]
% visible:...[DONE]
% type:...[DONE]
% orientation:...[DONE]

% MARKER:
% color: ...[DONE]
% size: ...[NA]
% symbol: ...[NA]
% opacity: ...[TODO]
% sizeref: ...[NA]
% sizemode: ...[NA]
% colorscale: ...[NA]
% cauto: ...[NA]
% cmin: ...[NA]
% cmax: ...[NA]
% outliercolor: ...[NA]
% maxdisplayed: ...[NA]

% MARKER LINE:
% color: ...[DONE]
% width: ...[DONE]
% dash: ...[NA]
% opacity: ...[TODO]
% shape: ...[NA]
% smoothing: ...[NA]
% outliercolor: ...[NA]
% outlierwidth: ...[NA]

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

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(histIndex).AssociatedAxis);

%-HIST DATA STRUCTURE- %
hist_data = get(obj.State.Plot(histIndex).Handle);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

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

%-hist xaxis-%
obj.data{histIndex}.xaxis = ['x' num2str(xsource)];

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

%-hist yaxis-%
obj.data{histIndex}.yaxis = ['y' num2str(ysource)];

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

%-hist type-%
obj.data{histIndex}.type = 'histogram';

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

%-HIST XAXIS-%
obj.data{histIndex}.histfunc= 'count';

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

orientation = histogramOrientation(hist_data);

switch orientation
    case 'v'

        %-hist x data-%
        xdata = mean(hist_data.XData(2:3,:));

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

        %-hist y data-%
        xlength = 0;
        for d = 1:length(xdata)
            obj.data{histIndex}.x(xlength + 1: xlength + hist_data.YData(2,d)) = repmat(xdata(d),1,hist_data.YData(2,d));
            xlength = length(obj.data{histIndex}.x);
        end

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

        %-hist autobinx-%
        obj.data{histIndex}.autobinx = false;

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

        %-hist xbins-%
        xbins.start = hist_data.XData(2,1);
        xbins.end = hist_data.XData(3,end);
        xbins.size = diff(hist_data.XData(2:3,1));
        obj.data{histIndex}.xbins = xbins;

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

    case 'h'

        %-hist y data-%
        ydata = mean(hist_data.YData(2:3,:));

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

        %-hist y data-%
        ylength = 0;
        for d = 1:length(ydata)
            obj.data{histIndex}.y(ylength + 1: ylength + hist_data.XData(2,d)) = repmat(ydata(d),1,hist_data.XData(2,d));
            ylength = length(obj.data{histIndex}.y);
        end

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

        %-hist autobiny-%
        obj.data{histIndex}.autobiny = false;

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

        %-hist ybins-%
        ybins.start = hist_data.YData(2,1);
        ybins.end = hist_data.YData(3,end);
        ybins.size = diff(hist_data.YData(2:3,1));
        obj.data{histIndex}.ybins = ybins;

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

end

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

%-hist name-%
obj.data{histIndex}.name = hist_data.DisplayName;

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

%-layout barmode-%
obj.layout.barmode = 'group';

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

%-layout bargap-%
obj.layout.bargap = (hist_data.XData(3,1)-hist_data.XData(2,2))/(hist_data.XData(3,1)-hist_data.XData(2,1));

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

%-hist line width-%
obj.data{histIndex}.marker.line.width = hist_data.LineWidth;

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

%-hist opacity-%
if ~ischar(hist_data.FaceAlpha)
    obj.data{histIndex}.opacity = hist_data.FaceAlpha;
end

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

%-hist marker-%
obj.data{histIndex}.marker = extractPatchFace(hist_data);

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

%-hist visible-%
obj.data{histIndex}.visible = strcmp(hist_data.Visible,'on');

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

%-hist showlegend-%
leg = get(hist_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{histIndex}.showlegend = showleg;

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

end

author: goGPS-Project


text:
load('net_weights.mat');
load('out_0_7layer_out.mat');
original_img = squeeze(original_img);
conv1_1 = squeeze(conv1_1);
conv1_2 = squeeze(conv1_2);
conv2_1 = squeeze(conv2_1);
conv2_2 = squeeze(conv2_2);
conv3_1 = squeeze(conv3_1);
conv3_2 = squeeze(conv3_2);
conv3_3 = squeeze(conv3_3);
conv4_1 = squeeze(conv4_1);
conv4_2 = squeeze(conv4_2);
conv4_3 = squeeze(conv4_3);
conv5_1 = squeeze(conv5_1);
conv5_2 = squeeze(conv5_2);
conv5_3 = squeeze(conv5_3);
img_gradient = squeeze(img_gradient);
pool1 = squeeze(pool1);
pool2 = squeeze(pool2);
pool3 = squeeze(pool3);
pool4 = squeeze(pool4);
pool_5 = squeeze(pool_5);
rpn_bbox_pred = squeeze(rpn_bbox_pred);
rpn_cls_prob = squeeze(rpn_cls_prob);
rpn_cls_prob_reshape = squeeze(rpn_cls_prob_reshape);
rpn_cls_score = squeeze(rpn_cls_score);
rpn_cls_score_reshape = squeeze(rpn_cls_score_reshape);
rpn_conv_3x3 = squeeze(rpn_conv_3x3);

%setup matconvnet
run /home/spc-public/Yixuan/matconvnet/matconvnet-1.0-beta24/matlab/vl_setupnn;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Function to load the stored optical flow based on some given method
%
%    Copyright (C) 2013  Anestis Papazoglou
%
%    You can redistribute and/or modify this software for non-commercial use
%    under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%    For commercial use, contact the author for licensing options.
%
%    Contact: a.papazoglou@sms.ed.ac.uk

function result = loadFlow( options)

    file = fullfile( options.outfolder, 'flow', 'flow.mat' );
    if( exist( file, 'file' ) )
        flow = load( file );
        if( isfield( flow, 'input' ) )
            result = flow.input;
        elseif( isfield( flow, 'flow' ) )
            result = flow.flow;
        else
            warning( '%s: no known field found\n', file );
            result = [];
        end
    else
        warning( '%s not found\n', file );
        result = [];
    end

end


text:
function [circle, labels] = imEnclosingCircle(img, varargin)
% Minimal enclosing circle of a region.
%
%   CIRC = imEnclosingCircle(IMG)
%   Computes the minimal enclosing circle around a binary region, or
%   around each labeled region in the input image.
%
%
%   CIRC = imEnclosingCircle(IMG, SPACING);
%   CIRC = imEnclosingCircle(IMG, SPACING, ORIGIN);
%   Specifies the spatial calibration of image. Both SPACING and ORIGIN are
%   1-by-2 row vectors. SPACING = [SX SY] contains the size of a pixel.
%   ORIGIN = [OX OY] contains the center position of the top-left pixel of
%   image.
%   If no calibration is specified, spacing = [1 1] and origin = [1 1] are
%   used. If only the sapcing is specified, the origin is set to [0 0].
%
%   CIRC = imEnclosingCircle(..., LABELS)
%   Specify the labels for which the enclosing circle needs to be computed.
%   The result is a N-by-3 array with as many rows as the number of labels.
%
%
%   Example
%   % Draw a commplex region together with its enclosing circle
%     img = imread('circles.png');
%     imshow(img); hold on;
%     circ = imEnclosingCircle(img);
%     drawCircle(circ)
%
%   % Compute and display the equivalent ellipses of several regions
%     img = imread('rice.png');
%     img2 = img - imopen(img, ones(30, 30));
%     lbl = bwlabel(img2 > 50, 4);
%     circles = imEnclosingCircle(lbl);
%     imshow(img); hold on;
%     drawCircle(circles, 'linewidth', 2, 'color', 'g');
%
%   See also
%     drawCircle, enclosingCircle
%     imInscribedCircle, imInertiaEllipse, imInertiaEllipse
%
%   References
%   Based on a file from Yazan Ahed (yash78@gmail.com)
%   which was rewritten from a Java applet by Shripad Thite:
%   http://heyoka.cs.uiuc.edu/~thite/mincircle/
%

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

%% Process input arguments

% default values
spacing = [1 1];
origin  = [1 1];
calib   = false;

% extract spacing
if ~isempty(varargin) && sum(size(varargin{1}) == [1 2]) == 2
    spacing = varargin{1};
    varargin(1) = [];
    calib = true;
    origin = [0 0];
end

% extract origin
if ~isempty(varargin) && sum(size(varargin{1}) == [1 2]) == 2
    origin = varargin{1};
end

% check if labels are specified
labels = [];
if ~isempty(varargin) && size(varargin{1}, 2) == 1
    labels = varargin{1};
end

%% Initialisations

% extract the set of labels, without the background
if isempty(labels)
    labels = imFindLabels(img);
end
nLabels = length(labels);

% allocate memory for result
circle = zeros(nLabels, 3);

%% Iterate over labels

for i = 1:nLabels
    % extract points of the current region
    [y, x] = find(img==labels(i));

    % works on convex hull (faster), or on original points if the hull
    % could not be computed
    try
        inds = convhull(x, y);
        pts = [x(inds) y(inds)];
    catch %#ok<CTCH>
        pts = [x y];
    end

    % works on pixel corners rather than on pixel centers
    x = pts(:,1);
    y = pts(:,2);
    pts = [x-0.5 y-0.5 ; x+0.5 y-0.5 ; x-0.5 y+0.5 ; x+0.5 y+0.5];

    % convex hull again, to remove duplicates
    inds = convhull(pts(:,1), pts(:,2));
    pts = pts(inds, :);

    % transform to physical space if needed
    if calib
        pts(:,1) = (pts(:,1)-1) * spacing(1) + origin(1);
        pts(:,2) = (pts(:,2)-1) * spacing(2) + origin(2);
    end

    circle(i,:) = recurseCircle(size(pts, 1), pts, 1, zeros(3, 2));
end

function circ = recurseCircle(n, p, m, b)
%    n: number of points given
%    m: an argument used by the function. Always use 1 for m.
%    bnry: an argument (3x2 array) used by the function to set the points that
%          determines the circle boundary. You have to be careful when choosing this
%          array's values. I think the values should be somewhere outside your points
%          boundary. For my case, for example, I know the (x,y) I have will be something
%          in between (-5,-5) and (5,5), so I use bnry as:
%                       [-10 -10
%                        -10 -10
%                        -10 -10]

if m == 4
    circ = createCircle(b(1,:), b(2,:), b(3,:));
    return;
end

circ = [Inf Inf 0];

if m == 2
    circ = [b(1,1:2) 0];
elseif m == 3
    c = (b(1,:) + b(2,:))/2;
    circ = [c distancePoints(b(1,:), c)];
end

for i = 1:n
    if distancePoints(p(i,:), circ(1:2)) > circ(3)
        if sum(b(:,1)==p(i,1) & b(:,2)==p(i,2)) == 0
            b(m,:) = p(i,:);
            circ = recurseCircle(i, p, m+1, b);
        end
    end
end

function dist = distancePoints(p1, p2)

dist = hypot(p2(:,1) - p1(:,1), p2(:,2) - p1(:,2));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function f = compose(f, op, g, data, pref)
%COMPOSE   Compose with a DELTAFUN is not allowed.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

error('CHEBFUN:DELTAFUN:compose:notSupported', ...
    'Composition with DELTAFUN is not supported.');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef Adam < ALGORITHM
% <single> <real> <large/none>
% Adaptive moment estimation
% alpha ---     1 --- Learning rate
% beta1 ---   0.9 --- A parameter within [0 1)
% beta2 --- 0.999 --- A parameter within [0 1)

%------------------------------- Reference --------------------------------
% D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv
% preprint arXiv:1412.6980, 2014.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [alpha,beta1,beta2] = Algorithm.ParameterSet(1,0.9,0.999);

            %% Generate random solution
            X  = Problem.Initialization(1);
            m0 = zeros(1,Problem.D);
            v0 = zeros(1,Problem.D);

            %% Optimization
            k = 1;
            while Algorithm.NotTerminated(X)
                gk = Problem.CalObjGrad(X.dec);
                m  = beta1*m0 + (1-beta1)*gk;
                v  = beta2*v0 + (1-beta2)*gk.^2;
                X  = Problem.Evaluation(X.dec-alpha*(m/(1-beta1.^k))./(sqrt(v/(1-beta2.^k))+1e-8));
                k  = k + 1;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function obj_io_test01 ( input_file_name )

%*****************************************************************************80
%
%% TEST01 tests OBJ_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  OBJ_SIZE determines the size of various objects\n' );
  fprintf ( 1, '  in an OBJ file.\n' );

  [ node_num, face_num, normal_num, order_max ] = obj_size ( input_file_name );

  obj_size_print ( input_file_name, node_num, face_num, normal_num, order_max );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [results, unmatched] = parse_params(args, defaults, varargin)
%PARSE_PARAMS Parses a set of name-value pairs.
% Usage:
%   parsed = parse_params(args, defaults)
%   parsed = parse_params(args, defaults, 'Name', Val, ...)
%   [parsed, unmatched] = parse_params(...)
%
% Args:
%   args: a cell array or structure of name-value pairs
%   defaults: a cell array or structure of name-value pairs
%
% Params:
%   This function takes any public properties of inputParser as parameters,
%   e.g., 'StructExpand'.
%
% Returns:
%   results: a structure containing the parsed results
%   unmatched: a structure containing any args that did not have a default
%
% See also: inputParser, struct2nameval, nameval2struct

% Make sure args are name-val cells
if isstruct(args); args = struct2nameval(args); end

% Make sure defaults are structs
if iscell(defaults); defaults = nameval2struct(defaults); end

% Default inputParser parameters
caller = get_caller_name();
parser_defaults.CaseSensitive = false;
parser_defaults.FunctionName = caller;
parser_defaults.KeepUnmatched = true;
parser_defaults.PartialMatching = true;
parser_defaults.StructExpand = true;

% Parse inputParser parameters
p = inputParser();
p.FunctionName = mfilename;
p = add_params(p, parser_defaults);
p.parse(varargin{:});
parser_params = p.Results;

% Create new inputParser instance
p = inputParser;

% Set inputParser properties from parsed parameters
param_names = fieldnames(parser_params);
for i = 1:numel(param_names)
    p.(param_names{i}) = parser_params.(param_names{i});
end

% Add parameters
p = add_params(p, defaults);

% Parse
p.parse(args{:});

% Return results of parsing
results = p.Results;
unmatched = p.Unmatched;

end

function p = add_params(p, defaults)

% Add defaults as parameters
names = fieldnames(defaults);
for i = 1:numel(names)
    p.addParameter(names{i}, defaults.(names{i}));
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function varargout = mlog(command,directory,varargin)
%function varargout = mlog(command,directory,varargin)
%
%Matlab-based log of files and folders
%
%This function currently supports two options:
%
%1> mlog('LOG',directory,save_flag,logName), %or
% [CERRVer,VerDate] = CERRCurrentVersion;
% mlog('LOG',getCERRPath,1,[CERRVer,'_', VerDate])
% directory: is the absolute path of the directory to be logged.
% save_flag: is a flag to tell the routine whether or not to save the log (1=save).
% logName: is the fileName to save log.
%2> [filesChangedC,filesRemovedC,filesAddedC,dirsChangedC,dirsRemovedC,dirsAddedC] = mlog('DIFF',getCERRPath);
%
%APA 05/27/2008
%
% 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/>.

switch upper(command)

    case 'LOG'
        clear logFilesAndDir
        logS = logFilesAndDir(directory);
        logDir = fullfile(directory,'Log','versions');
        if ~exist(logDir,'dir')
            mkdir(directory,'versions')
        end
        if length(varargin)>1 && varargin{1}==1 %(save_flag = varargin{1};)
            logName = varargin{2};
            save(fullfile(logDir,[logName,'.mat']),'logS')
        end
        varargout{1} = logS;

    case 'DIFF'
        hWait = waitbar(0,'Checking integrity of CERR m-files. Please wait...');
        logPath = fullfile(getCERRPath,'Log','versions');
        if ~exist(logPath,'dir')
            return;
        else
            %Get a fresh log
            clear logFilesAndDir
            logNewS = mlog('LOG',getCERRPath, 1, 'logCurrent');
            dirS = dir(logPath);
            dirS = dirS(3:end); %Assume 1st two directories are . and ..
        end
        [CERRVer, VerDate] = CERRCurrentVersion;
        %Find the log-file associated with the current version of CERR
        indFile = strmatch([CERRVer,'_', VerDate,'.mat'],{dirS.name},'exact');
        if ~isempty(indFile) %compare
            load(fullfile(logPath,dirS(indFile).name)) %loads logS variable
            [filesChangedC,filesRemovedC,filesAddedC,dirsChangedC,dirsRemovedC,dirsAddedC] = compareLogs(logS,logNewS,hWait);
            if isempty(filesChangedC) && isempty(filesRemovedC) && isempty(filesAddedC) && isempty(dirsChangedC) && isempty(dirsRemovedC) && isempty(dirsAddedC)
                disp('Nothing changed in CERR on your machine.')
            end
        else
            warning('Log not found...')
            close(hWait)
            varargout{1} = [];
            varargout{2} = [];
            varargout{3} = [];
            varargout{4} = [];
            varargout{5} = [];
            varargout{6} = [];
            return
        end

        varargout{1} = filesChangedC;
        varargout{2} = filesRemovedC;
        varargout{3} = filesAddedC;
        varargout{4} = dirsChangedC;
        varargout{5} = dirsRemovedC;
        varargout{6} = dirsAddedC;

end

return;
%% end of main

%---------------------- SUPPORTING FUNCTIONS

function [filesChangedC,filesRemovedC,filesAddedC,dirsChangedC,dirsRemovedC,dirsAddedC] = compareLogs(logOldS,logNewS,hWait)
%function compareLoggedDir(logOld,logNew)
%This function lists the files and directories in logOld which do not exist
%in logNew
%
%APA, 11/10/2006

%Check files
fNamesC = fieldnames(logOldS.fileLog);
numFiles = length(fNamesC);
dNamesC = fieldnames(logOldS.dirLog);
numDirs = length(dNamesC);
fileCheckV = zeros(numFiles,1);
indRemoved = [];

for i = numFiles:-1:1
    waitbar((numFiles-i)*1/(numFiles+numDirs),hWait)
    try
        %wy fileCheckV(i) = abs(logOldS.fileLog.(fNamesC{i}) - logNewS.fileLog.(fNamesC{i}));
        fileCheckV(i) = strcmpi(logOldS.fileLog.(fNamesC{i}), logNewS.fileLog.(fNamesC{i}));
    catch
        indRemoved = [indRemoved i];
    end
end

%%
%wy indChanged = find(fileCheckV>eps);
indChanged = find(fileCheckV==0);
filesChangedC = fNamesC(indChanged);
filesRemovedC = fNamesC(indRemoved);
fnamesNewC = fieldnames(logNewS.fileLog);
indAdded = ~ismember(fnamesNewC,fNamesC);
filesAddedC = fnamesNewC(indAdded);

%Check directories
dirCheckV = zeros(numDirs,1);
indRemoved = [];
for i = numDirs:-1:1
    waitbar(numFiles/(numFiles+numDirs) + (numDirs-i)*1/(numFiles+numDirs),hWait)
    try
        dirCheckV(i) = logOldS.dirLog.(dNamesC{i}) - logNewS.dirLog.(dNamesC{i});
    catch
        indRemoved = [indRemoved i];
    end
end
indChanged = find(dirCheckV>eps);
dirsChangedC = dNamesC(indChanged);
dirsRemovedC = dNamesC(indRemoved);
dnamesNewC = fieldnames(logNewS.dirLog);
indAdded = ~ismember(dnamesNewC,dNamesC);
dirsAddedC = dnamesNewC(indAdded);

close(hWait)

return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function result=cosmo_parcellfun(nproc,func,arg_cell,varargin)
% applies a function to elements in a cell in parallel
%
% result=cosmo_parcellfun(nproc,func,arg_cell,...)
%
% Inputs:
%   nproc                   Maximum number of processes to run in parallel
%   func                    Function handle that takes a single input
%                           argument and gives a single output
%   arg_cell                Cell with arguments to be given to func
%   'UniformOutput',o_u     If false, then the output is converted to a
%                           numeric or logical array. Default: true
%
% Output:
%   result                  Cell with the same size as arg_cell, with
%                               result{i}=func(arg_cell{i})
%                           If o_u is true, then result is converted to a
%                           numeric or logical array, assuming that each
%                           output is a scalar. If any output is not a
%                           scalar while o_u is true, an exception is
%                           thrown.
%

    default=struct();
    default.UniformOutput=true;

    opt=cosmo_structjoin(default,varargin{:});

    check_input(nproc,func,arg_cell,opt);

    % see how many processes to use
    narg_cell=numel(arg_cell);

    if narg_cell==1
        nproc_to_use=1;
    else
        nproc_available=cosmo_parallel_get_nproc_available(opt);
        nproc_to_use=min([nproc,narg_cell,nproc_available]);
    end

    if nproc_to_use<=1
        helper_func=@run_single_thread;
    else
        is_matlab=cosmo_wtf('is_matlab');
        if is_matlab
            helper_func=@run_parallel_matlab;
        else
            helper_func=@run_parallel_octave;
        end
    end

    result=helper_func(nproc_to_use,func,arg_cell,opt);

function args=get_extra_builtin_cellfun_args(opt)
    if opt.UniformOutput
        args={};
    else
        args={'UniformOutput',false};
    end

function result=run_single_thread(nproc,func,arg_cell,opt)
% single thread, Matlab of Octave --- redirect to cellfun
    result_cell=cellfun(func,arg_cell,'UniformOutput',false);
    result=convert_to_uniform_output_if_necessary(result_cell,opt);

function result=run_parallel_matlab(nproc,func,arg_cell,opt)
% multi-thread, Matlab
    narg_cell=numel(arg_cell);
    result_cell=cell(size(arg_cell));

    parfor (i=1:narg_cell, nproc)
        result_cell{i}=func(arg_cell{i});
    end

    result=convert_to_uniform_output_if_necessary(result_cell,opt);

function result=convert_to_uniform_output_if_necessary(result_cell,opt)
    if opt.UniformOutput
        is_uniform_func=@(x)isequal(size(x),[1 1]);
        is_uniform_mask=cellfun(is_uniform_func, result_cell);
        if ~all(is_uniform_mask)
            error(['Not all outputs are scalar. Use\n'...
                    '    ''UniformOutput'',false\n'...
                    'to return cell output']);
        end

        % concatenate and put in original shape
        result=reshape(cat(1,result_cell{:}),size(result_cell));
    else
        result=result_cell;
    end

function result=run_parallel_octave(nproc,func,arg_cell,opt)
% multi-thread, Octave
    extra_octave_args={'VerboseLevel',0};

    cellfun_args=get_extra_builtin_cellfun_args(opt);
    result=parcellfun(nproc,func,arg_cell,...
                                cellfun_args{:},...
                                extra_octave_args{:});

function check_input(nproc,func,arg_cell,opt)
    if ~(isnumeric(nproc) && ...
            isscalar(nproc) && ...
            round(nproc)==nproc && ...
            nproc>0)
        error(['nproc must be a positive integer. Use nproc=inf to use '...
                'as many processes as there are cores available']);
    end

    if ~isa(func,'function_handle')
        error('second argument must be a function handle');
    end

    if ~iscell(arg_cell)
        error('third argument must be a cell');
    end

    if ~(islogical(opt.UniformOutput) && ...
                isscalar(opt.UniformOutput))
        error('option ''UniformOutput'' must be scalar boolean');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
close all; clear all; clc

addpath('apm')

s = 'http://byu.apmonitor.com';
a = 'parameter_regression';

% clear application and load new model and data files
apm(s,a,'clear all');
apm_load(s,a,'model.apm');
csv_load(s,a,'data.csv');

% estimation mode
apm_option(s,a,'apm.imode',5);
apm_option(s,a,'apm.solver',3);

% classify variables
apm_info(s,a,'FV','K1');
apm_info(s,a,'FV','K2');
apm_info(s,a,'FV','K3');
apm_info(s,a,'FV','tau12');
apm_info(s,a,'FV','tau3');
apm_info(s,a,'MV','Q1');
apm_info(s,a,'MV','Q2');
apm_info(s,a,'CV','TC1');
apm_info(s,a,'CV','TC2');

% set status on
apm_option(s,a,'K1.status',1);
apm_option(s,a,'K2.status',1);
apm_option(s,a,'K3.status',1);
apm_option(s,a,'tau12.status',1);
apm_option(s,a,'tau3.status',0);

% set feedback status on
apm_option(s,a,'Q1.fstatus',1);
apm_option(s,a,'Q2.fstatus',1);
apm_option(s,a,'TC1.fstatus',1);
apm_option(s,a,'TC2.fstatus',1);

% optimize parameters
output = apm(s,a,'solve');
disp(output)

% retrieve solution
y = apm_sol(s,a);
z = y.x;

% optimized parameter values
K1 = apm_tag(s,a,'K1.newval');
K2 = apm_tag(s,a,'K2.newval');
K3 = apm_tag(s,a,'K3.newval');
tau12 = apm_tag(s,a,'tau12.newval');
tau3 = apm_tag(s,a,'tau3.newval');

% display values
disp(['K1:   ' num2str(K1)])
disp(['K2:   ' num2str(K2)])
disp(['K3:   ' num2str(K3)])
disp(['tau12:   ' num2str(tau12)])
disp(['tau3:   ' num2str(tau3)])

% read data.csv file for plotting
data = csvread('data.csv',1);
t = data(:,1);
T1meas = data(:,4);
T2meas = data(:,5);

% Plot results
figure(1)
subplot(3,1,1)
plot(t/60,T1meas,'b-','LineWidth',2)
hold on
plot(z.time/60,z.tc1,'g:','LineWidth',2)
ylabel('Temperature (degC)')
legend('T_1 measured','T_1 optimized')

subplot(3,1,2)
plot(t/60,T2meas,'k-','LineWidth',2)
hold on
plot(z.time/60,z.tc2,'r:','LineWidth',2)
ylabel('Temperature (degC)')
legend('T_2 measured','T_2 optimized')

subplot(3,1,3)
plot(z.time/60,z.q1,'g-','LineWidth',2)
hold on
plot(z.time/60,z.q2,'k--','LineWidth',2)
ylabel('Heater Output')
legend('Q_1','Q_2')

xlabel('Time (min)')

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function index = subv2ind(siz,sub)
%SUBV2IND   Linear index from subscript vector.
% SUBV2IND(SIZ,SUB) returns an equivalent single index corresponding to a
% subscript vector for an array of size SIZ.
% If SUB is a matrix, with subscript vectors as rows, then the result is a
% column vector.
%
% This is the opposite of IND2SUBV, so that
%   SUBV2IND(SIZ,IND2SUBV(SIZ,IND)) == IND.
%
% See also IND2SUBV, SUB2IND.

% Written by Tom Minka
% Part of Tom Minka's lightspeed package.
% (c) Microsoft Corporation. All rights reserved.

prev_cum_size = [1 cumprod(siz(1:end-1))];
%index = (sub-1)*prev_cum_size' + 1;
index = sub*prev_cum_size' - sum(prev_cum_size) + 1;

author: andrewssobral

repo: lrslibrary


text:
function h = quadsurf(x, y, z, c, varargin)
% QUADSURF  Quadrangular mesh surface plot
%
% quadsurf(X, Y, Z)
%
%   X, Y, Z are matrices with the coordinates of the mesh vertices. Note
%   that with this function it is possible to plot meshes that bend on
%   themselves or are rolled up. This is not possible with Matlab's
%   function surf(), that requires a unique value of Z for each (X,Y) pair.
%
%   This is an example of matrix X for a mesh with 10 rows and 7 columns.
%
%   X(1,1)-----X(1,2)---...---X(1,7)
%    |          |        |      |
%    |          |        |      |
%    |          |        |      |
%   X(2,1)-----X(2,2)---...---X(2,7)
%    |          |        |      |
%   ...        ...      ...    ...
%    |          |        |      |
%   X(10,1)----X(20,2)--...---X(20,7)
%
% quadsurf(X, Y, Z, C, <parameter/value pairs>)
%
%   C is a matrix of the same size as X, Y and Z, and specifies the colour
%   of the vertices by indexing into the colormap. By default, C=Z.
%
%   The X,Y,Z,C quad can be followed by parameter/value pairs to specify
%   additional properties of the Patch, e.g.
%
%     quadsurf(x, y, z, c, 'EdgeColor', 'white');
%
% H = quadsurf(...)
%
%   H returns a patch handle. Patches are children of AXES objects.
%
% See also: trisurf, trimesh, surf, mesh.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(3, Inf);
nargoutchk(0, 1);

% defaults
if (nargin < 4)
    c = z;
end

% size of the inputs
[R, C] = size(x);
if (any(size(y)~=[R, C]) || any(size(z)~=[R, C]))
    error('X, Y, Z must have the same size')
end

% patch takes quadrangles as rows of 4 elements. Here we create matrices
% with the same size as X,Y,Z, for each of the 4 vertices of a quadrangle.
% Note that we have a lot of redundancy, but this is the patch convention
%
% vx1       vx2
%  ----------
%  |        |
%  |        |
%  |        |
%  |        |
%  ----------
% vx4       vx3
vx1 = x(1:R-1, 1:C-1);
vx2 = x(1:R-1, 2:C);
vx3 = x(2:R, 2:C);
vx4 = x(2:R, 1:C-1);
vy1 = y(1:R-1, 1:C-1);
vy2 = y(1:R-1, 2:C);
vy3 = y(2:R, 2:C);
vy4 = y(2:R, 1:C-1);
vz1 = z(1:R-1, 1:C-1);
vz2 = z(1:R-1, 2:C);
vz3 = z(2:R, 2:C);
vz4 = z(2:R, 1:C-1);
c1 = c(1:R-1, 1:C-1);
c2 = c(1:R-1, 2:C);
c3 = c(2:R, 2:C);
c4 = c(2:R, 1:C-1);

% format the vertices matrices, and create the patch image
h = patch([vx1(:)'; vx2(:)'; vx3(:)'; vx4(:)'], ...
    [vy1(:)'; vy2(:)'; vy3(:)'; vy4(:)'], ...
    [vz1(:)'; vz2(:)'; vz3(:)'; vz4(:)'], ...
    [c1(:)'; c2(:)'; c3(:)'; c4(:)'], varargin{:});

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus