text:
function calpak_test75 ( )

%*****************************************************************************80
%
%% CALPAK_TEST75 tests YMD_INC_YMD_COMMON and YMDF_DIF_YMDF_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST75\n' );
  fprintf ( 1, '  For the Common calendar,\n' );
  fprintf ( 1, '  YMD_INC_YMD_COMMON increments a YMDF date by YMDF;\n' );
  fprintf ( 1, '  YMDF_DIF_YMDF_COMMON finds the YMDF difference.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Date1        increment    Date2          difference\n' );
  fprintf ( 1, '\n' );

  y1 = 1900;
  m1 = 5;
  d1 = 27;
  f1 = 0.0;

  yn1 = 50;
  mn1 = 9;
  dn1 = 10;
  fn2 = 0.0;

  [ y2, m2, d2 ] = ymd_inc_ymd_common ( y1, m1, d1, yn1, mn1, dn1 );
  f2 = 0.0;

  [ yn2, mn2, dn2, fn2, ierror ] = ymdf_dif_ymdf_common ( y1, m1, d1, f1, ...
  y2, m2, d2, f2 );

  s1 = ymdf_to_s_common ( y1, m1, d1, f1 );

  s2 = ymdf_to_s_common ( y2, m2, d2, f2 );

  fprintf ( 1, '  %15s  %3d  %3d  %3d  %15s  %3d  %3d  %3d\n', ...
    s1, yn1, mn1, dn1, s2, yn2, mn2, dn2 );

  return
end


text:
function [model, pca_model] = project_model(model, coeff, k)
% [model, pcamodel] = project_model(model, coeff, k)
%
% Project a model's filters onto the top k PCA eigenvectors
% stored in the columns of the matrix coeff.  The output
% variable 'model' holds the original model augmented to
% hold the PCA filters as extra data.  The output variable
% 'pcamodel' has its filters replaced with the PCA filters.

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% take the top k eigenvectors from coeff as the projection matrix
coeff = coeff(:, 1:k);
% augment the projection matrix by adding a vector with all zeros ...
coeff = padarray(coeff, [1 1], 0, 'post');
% ... except in the last position to preserve the occlusion feature
coeff(end,end) = 1;
% save the projection matrix in the model
model.pca_coeff = coeff;
% Make a new model with projected filters
pca_model = model;
for i = 1:model.numfilters
  bl = model.filters(i).blocklabel;
  % w is reshaped and appropriately flipped
  w = model_get_block(model, model.filters(i));
  w_pca = project(w, coeff);
  if model.filters(i).flip
    pca_model.blocks(bl).type      = block_types.PCAFilter;
    pca_model.blocks(bl).shape(3)  = k+1;
    pca_model.blocks(bl).w_flipped = w_pca;

    model.blocks(bl).w_pca_flipped = w_pca;
    model.blocks(bl).w_flipped     = w;
  else
    pca_model.blocks(bl).type      = block_types.PCAFilter;
    pca_model.blocks(bl).shape(3)  = k+1;
    pca_model.blocks(bl).w         = w_pca;

    model.blocks(bl).w_pca         = w_pca;
    model.blocks(bl).w             = w;
  end
  pca_model.blocks(bl).dim = numel(w_pca);
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function cost = computeBoundariesCost(opts, X, theta_x, theta_y, theta_z, T)
    H = eye(4);
    H(1:3,1:3) = rotx(theta_x) * roty(theta_y) * rotz(theta_z);
    H(1:3,4) = T';
    X_transformed = H * X;
    cost_x_pos = sum(X_transformed(1, X_transformed(1,:)>0));
    cost_x_neg = sum(X_transformed(1, X_transformed(1,:)<0));

    cost_y_pos = sum(X_transformed(2, X_transformed(2,:)>0));
    cost_y_neg = sum(X_transformed(2, X_transformed(2,:)<0));

    cost_z_pos = sum(X_transformed(3, X_transformed(3,:)>0));
    cost_z_neg = sum(X_transformed(3, X_transformed(3,:)<0));
    figure(1)
    hold on
    scatter3(X_transformed(1,:), X_transformed(2,:), X_transformed(3,:), 'r.')
end

author: UMich-BipedLab


text:
function str = cfg_textfill(obj, left, right, tflag)

% function str = cfg_textfill(obj, left, right)
% Fill a text object, so that the left part is left justified and the
% right part right justified. If tflag is set, try to fit text in widget
% by truncating right until at least 5 characters are displayed.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_textfill.m 7473 2018-11-06 10:26:44Z guillaume $

rev = '$Rev: 7473 $'; %#ok

if ~ishandle(obj)
    cfg_message('matlabbatch:usage',...
          'First input must be a graphics handle.');
else
    try
        get(obj,'string');
    catch
        cfg_message('matlabbatch:usage',...
              'Input object must have a ''string'' property.');
    end;
end;
if ~iscellstr(left)
    if ischar(left)
        left = cellstr(left);
    else
        cfg_message('matlabbatch:usage',...
              'Second input must be a string array or cellstr.');
    end;
    % add one space as delimiter
    left = strcat(left, {' '});
end;
if ~iscellstr(right)
    if ischar(right)
        right = cellstr(right);
    else
        cfg_message('matlabbatch:usage',...
              'Third input must be a string array or cellstr.');
    end;
end;
if numel(left) ~= numel(right)
    cfg_message('matlabbatch:usage',...
          'Second and third input must have the same number of lines.');
end;

if exist('OCTAVE_VERSION', 'builtin'), set(0,'currentfigure',ancestor(obj,'figure')); end
TempObj=copyobj(obj,get(obj,'Parent'));
set(TempObj,'Visible','off','Max',100);

% Find max extent of left string
lext = cfg_maxextent(TempObj, left);
mlext = max(lext);

% Find max extent of right string
rext = cfg_maxextent(TempObj, right);
mrext = max(rext);

% Find extent of single space
% Work around MATLAB inaccuracy by measuring 100 spaces and dividing by 100
spext = cfg_maxextent(TempObj, {repmat(' ',1,100)})/100;

% try to work out slider size
pos = get(TempObj,'Position');
oldun = get(TempObj,'units');
set(TempObj,'units','points');
ppos = get(TempObj,'Position');
set(TempObj,'units',oldun);
sc = pos(3)/ppos(3);
% assume slider width of 15 points
swidth=15*sc;

pos = get(obj,'Position');
width = max(mlext+mrext,pos(3)-swidth);
if tflag && mlext+mrext > pos(3)-swidth
    newrext = max(5*spext, pos(3)-swidth-mlext);
    trind = find(rext > newrext);
    for k = 1:numel(trind)
        tright = right{trind(k)};
        set(TempObj,'String',['...' tright]);
        ext = get(TempObj,'Extent');
        while ext(3) > newrext
            tright = tright(2:end);
            set(TempObj,'String',['...' tright]);
            ext = get(TempObj,'Extent');
        end;
        right{trind(k)} = ['...' tright];
        rext(trind(k)) = ext(3);
    end;
    % Re-estimate max extent of trimmed right string
    rext = cfg_maxextent(TempObj, right);
    mrext = max(rext);
    width = mlext+mrext;
end;

fillstr = cell(size(left));
for k = 1:numel(left)
    fillstr{k} = repmat(' ',1,floor((width-(lext(k)+rext(k)))/spext));
end;
str = strcat(left, fillstr, right);
delete(TempObj);

author: spm


text:
function [J, Ybus, Yf, Yt] = makeJac(baseMVA, bus, branch, gen, fullJac)
%MAKEJAC  Forms the power flow Jacobian.
%   J = MAKEJAC(MPC)
%   J = MAKEJAC(MPC, FULLJAC)
%   J = MAKEJAC(BASEMVA, BUS, BRANCH, GEN)
%   J = MAKEJAC(BASEMVA, BUS, BRANCH, GEN, FULLJAC)
%   [J, YBUS, YF, YT] = MAKEJAC(MPC)
%
%   Returns the power flow Jacobian and, optionally, the system admittance
%   matrices. Inputs can be a MATPOWER case struct or individual BASEMVA,
%   BUS, BRANCH and GEN values. Bus numbers must be consecutive beginning
%   at 1 (i.e. internal ordering). If the FULLJAC argument is present and
%   true, it returns the full Jacobian (sensitivities of all bus injections
%   w.r.t all voltage angles/magnitudes) as opposed to the reduced version
%   used in the Newton power flow updates. The units for all quantities are
%   in per unit with radians for voltage angles.
%
%   Note: This function builds the Jacobian from scratch, rebuilding the
%         YBUS matrix in the process. You probably don't want to use this
%         in performance critical code.
%
%   See also MAKEYBUS, EXT2INT

%   MATPOWER
%   Copyright (c) 1996-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.

if nargin < 4
    mpc     = baseMVA;
    if nargin > 1
        fullJac = bus;
    else
        fullJac = 0;
    end
    baseMVA = mpc.baseMVA;
    bus     = mpc.bus;
    branch  = mpc.branch;
    gen     = mpc.gen;
elseif nargin < 5
    fullJac = 0;
end

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%% build Ybus
[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);

%% extract voltage
V = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));

%% make sure we use generator setpoint voltage for PV and slack buses
on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
gbus = gen(on, GEN_BUS);                %% what buses are they at?
k = find(bus(gbus, BUS_TYPE) == PV | bus(gbus, BUS_TYPE) == REF);
V(gbus(k)) = gen(on(k), VG) ./ abs(V(gbus(k))).* V(gbus(k));

%% build Jacobian
[dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V);
if fullJac
    j11 = real(dSbus_dVa);
    j12 = real(dSbus_dVm);
    j21 = imag(dSbus_dVa);
    j22 = imag(dSbus_dVm);
else
    %% get bus index lists of each type of bus
    [ref, pv, pq] = bustypes(bus, gen);

    j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
    j12 = real(dSbus_dVm([pv; pq], pq));
    j21 = imag(dSbus_dVa(pq, [pv; pq]));
    j22 = imag(dSbus_dVm(pq, pq));
end

J = [   j11 j12;
        j21 j22;    ];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%Extracts the grid for the grid defined on
%http://www.vision.ee.ethz.ch/~cwengert/calibration_toolbox.php
%
%Input:         im      The image to process
%               show    0 to not display any information graphically
%               minArea The minimum area a grid dot occupies
%               maxArea The minimum area a grid dot occupies
%Output:        success 1 if everything went ok
%               stats   The processed blobs of this image
%               id0     The index of the main bar
%               id1     The index of the second main bar
%               searchPoints    The potential grid points
%
%
%Christian Wengert
%Computer Vision Laboratory
%ETH Zurich
%Sternwartstrasse 7
%CH-8092 Zurich
%www.vision.ee.ethz.ch/cwengert
%wengert@vision.ee.ethz.ch

function [success, stats, cnt, id0, id1, searchPoints] = gridextractor(im, show, minArea, maxArea)
    %crop image
    if(nargin==1)
        show=0;
    end
    success = 0;
    stats = [];
    cnt = 0;
    searchPoints = [];
    %Indices for main bars
    idx = [];
    ids = [];
    id0 = -1;
    id1 = -1;
    x_ = [];
    X_ = [];
    %Show the image
    if(show)
        if(isfloat(im))
             imshow(im/255.0),hold on
        else
            imshow(im),hold on
        end
    end
    %Get the size of the image
    [h,w,bpp] = size(im);
    if(h<=0 & w<=0)
        return
    end
    if(bpp==3)
        img = rgb2gray(im);
    else
        img = im;
    end

    %used to filter for ellipticity
    minEllipticity = 1.95;
    maxEllipticity = 8;
    minSolidity = 0.80;
     %Used to filter by closeness to border->to close to border cannot be
     %correct
    borderRect = [0+50, 0+50, w-50, h-50];

    x0 = w/2;
    y0 = h/2;
    imc = im;
    for i=1:w
        for j=1:h
            r = sqrt((x0-i)*(x0-i) + (y0-j)*(y0-j));
            if(r>160)
                imc(j,i) = 0;
            end
        end
    end

    %Rather median filter
    imb = medfilt2(img,[3 3]);

    %Threshold
    h1 = fspecial('gaussian',100, 15);
    imf= (imfilter(im, h1));
    imt = ( imclose(medfilt2(im<imf,[3 3]), ones(5)))    ;
    %Label
    iml = bwlabel(imt,8);
    %Extract features
    stats=regionprops(iml,'Centroid','Area','Orientation','MajorAxisLength','MinorAxisLength','ConvexHull','BoundingBox','Solidity');

    %Extract values and put them into vectors, just makes it easier
    areas = [stats.Area];
    pts = [stats.Centroid];
    ptsx = pts(1:2:end);
    ptsy = pts(2:2:end);
    majors = [stats.MajorAxisLength];
    minors = [stats.MinorAxisLength];
    solidity = [stats.Solidity];

    %Now get the two major things
    cnt = 0;
    for i=1:length(stats)
        %Filter on closeness to border
        if(stats(i).Centroid(1)>borderRect(1) & stats(i).Centroid(1)<borderRect(3) & stats(i).Centroid(2) > borderRect(2) & stats(i).Centroid(2) <borderRect(4))
            %Filter area
            if(areas(i)>minArea & areas(i)<maxArea & solidity(i) > minSolidity)
                idx = [idx;i];
                %Filter on ellipticity
                if((majors(i)/minors(i))>minEllipticity & (majors(i)/minors(i))<maxEllipticity )
                    if(show)
                        plot(stats(i).Centroid(1),stats(i).Centroid(2),'r+')
                    end
                    ids = [ids;i];
                    cnt = cnt+1;
                    %Delete it from searchPoints, serves as an ellipticity
                    %filter!
                    idx(end) = [];
                end
            end
        end
    end
    %Show area filtered blobs
    if(show)
          for i=1:length(idx)
            plot(stats(idx(i)).Centroid(1),stats(idx(i)).Centroid(2),'go')
        end
    end
    try
        %Extract the main bars
        if(cnt==2)
            if(areas(ids(1))>areas(ids(2)))
                id0 = ids(1); id1 = ids(2);
            else
                id0 = ids(2); id1 = ids(1);
            end
            %Test angle between main lines
            dir0 = (stats(id1).Centroid-stats(id0).Centroid);
            dir0 = dir0/norm(dir0);
            %get the orientation of main bar
            orientation = stats(id0).Orientation;
            y = tand(orientation);
            dir90 = [-1,y];
            dir90 = dir90/norm(dir90);
            %Check angle
            theta = acos((dir0*dir90')/(norm(dir0)*norm(dir90)))/pi*180;
            if(theta<70 | theta>110)
                 disp('gridextractor:: Could not find both main bars, bar-count = 1. ')
                 stats = [];
                 searchPoints = [];
                 success = 0;
                 cnt = 0;
                 id0 = 0;
                 id1 = 0;
            else
                searchPoints = [areas(idx);...
                                ptsx(idx);...
                                ptsy(idx)];
                success = 1;
                if(show)
                    text(stats(id0).Centroid(1)+5,stats(id0).Centroid(2)+5,'0','Color','m');
                    text(stats(id1).Centroid(1)+5,stats(id1).Centroid(2)+5,'1','Color','m');
                end
            end
        elseif(cnt>2)
            disp('gridextractor:: Found too many possible candidates, trying to choose the correct ones')
            %get the two blobs that are closest to each other!
            a = [];
            xb = [];
            %Sort by Area
            for i=1:length(ids)
                xb = [xb, [stats(ids(i)).Centroid(1);stats(ids(i)).Centroid(2)]];
                if(xb(1,i) > 50 | xb(1,i) < (w-50) | xb(2,i) > 50 | xb(2,i) < (h-50))
                    a = [a;stats(ids(i)).Area];
                else
                    a = [a;-1];
                end
            end

            [a,ixxx] = sort(a);
            %Take the two biggest ones
            id0 = ids(ixxx(end));
            id1 = ids(ixxx(end-1));

             %Directions
            dir0 = (stats(id1).Centroid-stats(id0).Centroid);
            dir0 = dir0/norm(dir0);
            %get the orientation of main bar
            orientation = stats(id0).Orientation;
            y = tand(orientation);
            dir90 = [-1,y];
            dir90 = dir90/norm(dir90);
            %Check angle
            theta = acos((dir0*dir90')/(norm(dir0)*norm(dir90)))/pi*180;
            if(theta<80 | theta>100)
                id1 = ids(ixxx(end-2));
                %Test again
                dir0 = (stats(id1).Centroid-stats(id0).Centroid);
                dir0 = dir0/norm(dir0);
                %get the orientation of main bar
                orientation = stats(id0).Orientation;
                y = tand(orientation);
                dir90 = [-1,y];
                dir90 = dir90/norm(dir90);
                %Check angle
                theta = acos((dir0*dir90')/(norm(dir0)*norm(dir90)))/pi*180;
                if(theta<80 | theta>100)
                     disp('gridextractor:: Could not find both main bars, bar-count = 1. ')
                     stats = [];
                     searchPoints = [];
                     success = 0;
                     cnt = 0;
                     id0 = 0;
                     id1 = 0;
                end
            end

            if(show)
                text(stats(id0).Centroid(1)+5,stats(id0).Centroid(2)+5,'0','Color','y');
                text(stats(id1).Centroid(1)+5,stats(id1).Centroid(2)+5,'1','Color','y');
            end
            cnt = 2;
            success = 1;
            searchPoints = [areas(idx);...
                            ptsx(idx);...
                            ptsy(idx)];
        else
            disp('gridextractor:: Could not find both main bars, bar-count = 1. ')
            stats = [];
            searchPoints = [];
            success = 0;
            cnt = 0;
            id0 = 0;
            id1 = 0;
        end
    catch
        disp('gridextractor:: Problems extracting the main bars. Please verify whether it is visible in this image.')
        s = lasterror;
        disp(['gridextractor:: ' s.message ' in function ' s.stack.name ' in ' s.stack.file ' @line ' num2str(s.stack.line)])
        stats = [];
        searchPoints = [];
        success = 0;
        cnt = 0;
        id0 = 0;
        id1 = 0;
    end

author: christianwengert


text:
function [data,units] = compute_darea_inmost_wing(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  dareal = diff(trx(fly).wing_areal_mm);
  darear = diff(trx(fly).wing_arear_mm);
  data{i} = dareal;
  idx = trx(fly).wing_angler(1:end-1) <= -trx(fly).wing_anglel(1:end-1);
  data{i}(idx) = darear(idx);

  data{i} = data{i} ./ trx(fly).dt;
end
units = parseunits('mm^2/s');

author: kristinbranson

repo: JAABA


text:
function [scores, rects, pts, pts_valid] = detect_img( param, caffe_solver, thres )
    num = length(list);
    test_time = ceil(num / length(param.gpu_id));
    for i = 1 : test_time
        target_scale_now = param.max_img *
        for j = 1 : length(param.gpu_id)
            [ img, scale ] = scale_img_to_target( img, target )
        end
    end

end

author: liuyuisanai

repo: RSA-for-object-detection


text:
%UNDISTORT  Transforms an image to compensate for lens distortion
%
%     dst = cv.undistort(src, cameraMatrix, distCoeffs)
%     dst = cv.undistort(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ Input (distorted) image.
% * __cameraMatrix__ Input camera matrix `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
% * __distCoeffs__ Input vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy]` of 4, 5, 8, 12 or 14
%   elements. If the vector is empty, the zero distortion coefficients are
%   assumed.
%
% ## Output
% * __dst__ Output (corrected) image that has the same size and type as `src`.
%
% ## Options
% * __NewCameraMatrix__ Camera matrix of the distorted image. By default, it
%   is the same as `cameraMatrix` but you may additionally scale and shift the
%   result by using a different matrix. Not set by default
%
% The function transforms an image to compensate radial and tangential lens
% distortion.
%
% The function is simply a combination of cv.initUndistortRectifyMap (with
% unity `R`) and cv.remap (with bilinear interpolation). See the former
% function for details of the transformation being performed.
%
% Those pixels in the destination image, for which there is no correspondent
% pixels in the source image, are filled with zeros (black color).
%
% A particular subset of the source image that will be visible in the
% corrected image can be regulated by `NewCameraMatrix`. You can use
% cv.getOptimalNewCameraMatrix to compute the appropriate `NewCameraMatrix`
% depending on your requirements.
%
% The camera matrix and the distortion parameters can be determined using
% cv.calibrateCamera. If the resolution of images is different from the
% resolution used at the calibration stage, `fx`, `fy`, `cx` and `cy` need to
% be scaled accordingly, while the distortion coefficients remain the same.
%
% See also: cv.undistortPoints, cv.calibrateCamera,
%  cv.initUndistortRectifyMap, cv.remap, undistortImage
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function normal_01_cdf_values_test ( )

%*****************************************************************************80
%
%% NORMAL_01_CDF_VALUES_TEST demonstrates the use of NORMAL_01_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NORMAL_01_CDF_VALUES_TEST:\n' );
  fprintf ( 1, '  NORMAL_01_CDF_VALUES stores values of\n' );
  fprintf ( 1, '  the Normal 01 Cumulative Density Function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X                  CDF(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = normal_01_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %24.16f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function M = fixedrankfactory_3factors_preconditioned(m, n, k)
% Manifold of m-by-n matrices of rank k with three factor quotient geometry.
%
% function M = fixedrankfactory_3factors_preconditioned(m, n, k)
%
% This geometry is tuned to least squares problems such as low-rank matrix
% completion with ell-2 loss.
%
% A point X on the manifold is represented as a structure with three
% fields: L, S and R. The matrices L (mxk) and R (nxk) are orthonormal,
% while the matrix S (kxk) is a full rank matrix such that X = L*S*R'.
%
% Tangent vectors are represented as a structure with three fields: L, S
% and R.
%
% Please cite the Manopt paper as well as the research paper:
%     @InProceedings{mishra2014r3mc,
%       Title        = {{R3MC}: A {R}iemannian three-factor algorithm for low-rank matrix completion},
%       Author       = {Mishra, B. and Sepulchre, R.},
%       Booktitle    = {{53rd IEEE Conference on Decision and Control}},
%       Year         = {2014},
%       Organization = {{IEEE CDC}}
%     }
%
%
% See also: fixedrankfactory_3factors fixedrankfactory_2factors_preconditioned

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:
%
%	April 04, 2015 (BM):
%       Cosmetic changes including avoiding storing the inverse of a kxk matrix.

    M.name = @() sprintf('LSR'' (tuned for least square problems) quotient manifold of %dx%d matrices of rank %d', m, n, k);

    M.dim = @() (m+n-k)*k;

    % Some precomputations at the point X that are to be used in the inner product (and
    % pretty much everywhere else).
    function X = prepare(X)
        if ~all(isfield(X,{'StS','SSt'}) == 1)
            X.SSt = X.S*X.S';
            X.StS = X.S'*X.S;
        end
    end

    % The choice of metric is motivated by symmetry and tuned to least square
    % objective function.
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        X = prepare(X);

        ip = trace(X.SSt*(eta.L'*zeta.L)) + trace(X.StS*(eta.R'*zeta.R)) ...
            + trace(eta.S'*zeta.S);
    end

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(x, y) error('fixedrankfactory_3factors_preconditioned.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    skew = @(X) .5*(X-X');
    symm = @(X) .5*(X+X');

    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)
        X = prepare(X);

        SSL = X.SSt;
        ASL = 2*symm(SSL*(egrad.S*X.S'));

        SSR = X.StS;
        ASR = 2*symm(SSR*(egrad.S'*X.S));

        [BL, BR] = tangent_space_lyap(X.S, ASL, ASR); % It computes the solution without calling Matlab's Lyap.

        rgrad.L = (egrad.L - X.L*BL)/X.SSt;
        rgrad.R = (egrad.R - X.R*BR)/X.StS;
        rgrad.S = egrad.S;

        % Debug
        %         BL1 = lyap(SSL, -ASL); % Alternate way
        %         BR1 = lyap(SSR, -ASR);
        %         norm(skew(X.SSt*(rgrad.L'*X.L) + rgrad.S*X.S'), 'fro')
        %         norm(skew(X.StS*(rgrad.R'*X.R) - X.S'*rgrad.S), 'fro')

    end

    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)
        X = prepare(X);

        % Riemannian gradient.
        SSL = X.SSt;
        ASL = 2*symm(SSL*(egrad.S*X.S'));
        SSR = X.StS;
        ASR = 2*symm(SSR*(egrad.S'*X.S));
        [BL, BR] = tangent_space_lyap(X.S, ASL, ASR);

        rgrad.L = (egrad.L - X.L*BL)/X.SSt;
        rgrad.R = (egrad.R - X.R*BR)/X.StS;
        rgrad.S = egrad.S;

        % Directional derivative of the Riemannian gradient.
        ASLdot = 2*symm((2*symm(X.S*eta.S')*(egrad.S*X.S')) + X.SSt*(ehess.S*X.S' + egrad.S*eta.S')) - 4*symm(symm(eta.S*X.S')*BL);
        ASRdot = 2*symm((2*symm(X.S'*eta.S)*(egrad.S'*X.S)) + X.StS*(ehess.S'*X.S + egrad.S'*eta.S)) - 4*symm(symm(eta.S'*X.S)*BR);

        %         SSLdot = X.SSt;
        %         SSRdot = X.StS;
        %         BLdot = lyap(SSLdot, -ASLdot);
        %         BRdot = lyap(SSRdot, -ASRdot);

        [BLdot, BRdot] = tangent_space_lyap(X.S, ASLdot, ASRdot);

        Hess.L = (ehess.L - eta.L*BL - X.L*BLdot - 2*rgrad.L*symm(eta.S*X.S'))/X.SSt;
        Hess.R = (ehess.R - eta.R*BR - X.R*BRdot - 2*rgrad.R*symm(eta.S'*X.S))/X.StS;
        Hess.S = ehess.S;

        % BM: Till this, everything seems correct.
        % We still need a correction factor for the non-constant metric
        % that is imposed.
        % The computation of the correction factor owes itself to the Koszul formula.
        % This corresponds to the Riemannian connection in the Euclidean space with the
        % scaled metric.
        Hess.L = Hess.L + (eta.L*symm(rgrad.S*X.S') + rgrad.L*symm(eta.S*X.S'))/X.SSt;
        Hess.R = Hess.R + (eta.R*symm(rgrad.S'*X.S) + rgrad.R*symm(eta.S'*X.S))/X.StS;
        Hess.S = Hess.S - symm(rgrad.L'*eta.L)*X.S - X.S*symm(rgrad.R'*eta.R);

        % The Riemannian connection on the quotient space is the
        % projection of the Riemmian connection in the ambient space onto the tangent space of the total space and
        % then onto the horizontal space.
        % This is accomplished by the following operation.
        Hess = M.proj(X, Hess);

        % Debug
        %         norm(skew(X.SSt*(Hess.L'*X.L) + Hess.S*X.S'))
        %         norm(skew(X.StS*(Hess.R'*X.R) - X.S'*Hess.S))

    end

    M.proj = @projection;
    function etaproj = projection(X, eta)
        X = prepare(X);

        % First, projection onto the tangent space of the total space.
        SSL = X.SSt;
        ASL = 2*symm(X.SSt*(X.L'*eta.L)*X.SSt);
        BL = lyap(SSL, -ASL);
        eta.L = eta.L - X.L*(BL/X.SSt);

        SSR = X.StS;
        ASR = 2*symm(X.StS*(X.R'*eta.R)*X.StS);
        BR = lyap(SSR, -ASR);
        eta.R = eta.R - X.R*(BR/X.StS);

        % Project onto the horizontal space
        PU = skew((X.L'*eta.L)*X.SSt) + skew(X.S*eta.S');
        PV = skew((X.R'*eta.R)*X.StS)  + skew(X.S'*eta.S);
        [Omega1, Omega2] = coupled_lyap(X.S, PU, PV);
        %         norm(2*skew(Omega1*X.SSt) - PU -(X.S*Omega2*X.S'),'fro' )
        %         norm(2*skew(Omega2*X.StS) - PV -(X.S'*Omega1*X.S),'fro' )
        %

        etaproj.L = eta.L - (X.L*Omega1);
        etaproj.S = eta.S - (X.S*Omega2 - Omega1*X.S) ;
        etaproj.R = eta.R - (X.R*Omega2);

        % Debug
        %         norm(skew(X.SSt*(etaproj.L'*X.L) + etaproj.S*X.S'))
        %         norm(skew(X.StS*(etaproj.R'*X.R) - X.S'*etaproj.S))
        %
        %         norm(skew(X.SSt*(etaproj.L'*X.L) - X.S*etaproj.S'))
        %         norm(skew(X.StS*(etaproj.R'*X.R) + etaproj.S'*X.S))

    end

    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;

    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        Y.S = (X.S + t*eta.S);
        Y.L = uf((X.L + t*eta.L));
        Y.R = uf((X.R + t*eta.R));

        Y = prepare(Y);
    end

    M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y = retraction(X, eta, t);
        warning('manopt:fixedrankfactory_3factors_preconditioned:exp', ...
            ['Exponential for fixed rank ' ...
            'manifold not implemented yet. Used retraction instead.']);
    end

    M.hash = @(X) ['z' hashmd5([X.L(:) ; X.S(:) ; X.R(:)])];

    M.rand = @random;
    % Factors L and R live on Stiefel manifolds, hence we will reuse
    % their random generator.
    stiefelm = stiefelfactory(m, k);
    stiefeln = stiefelfactory(n, k);
    function X = random()
        X.L = stiefelm.rand();
        X.R = stiefeln.rand();
        X.S = diag(1+rand(k, 1));

        X = prepare(X);
    end

    M.randvec = @randomvec;
    function eta = randomvec(X)
        % A random vector on the horizontal space
        eta.L = randn(m, k);
        eta.R = randn(n, k);
        eta.S = randn(k, k);
        eta = projection(X, eta);
        nrm = M.norm(X, eta);
        eta.L = eta.L / nrm;
        eta.R = eta.R / nrm;
        eta.S = eta.S / nrm;
    end

    M.lincomb = @lincomb;

    M.zerovec = @(X) struct('L', zeros(m, k), 'S', zeros(k, k), ...
        'R', zeros(n, k));

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

    % vec and mat are not isometries, because of the unusual inner metric.
    M.vec = @(X, U) [U.L(:) ; U.S(:); U.R(:)];
    M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
        'S', reshape(u((m*k+1): m*k + k*k), k, k), ...
        'R', reshape(u((m*k+ k*k + 1):end), n, k));
    M.vecmatareisometries = @() false;

end

% Linear combination of tangent vectors
function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>

    if nargin == 3
        d.L = a1*d1.L;
        d.R = a1*d1.R;
        d.S = a1*d1.S;
    elseif nargin == 5
        d.L = a1*d1.L + a2*d2.L;
        d.R = a1*d1.R + a2*d2.R;
        d.S = a1*d1.S + a2*d2.S;
    else
        error('Bad use of fixedrankfactory_3factors_preconditioned.lincomb.');
    end

end

function A = uf(A)
    [L, unused, R] = svd(A, 0); %#ok
    A = L*R';
end

function[BU, BV] = tangent_space_lyap(R, E, F)
    % We intent to solve a linear system    RR^T  BU + BU RR^T  = E
    %                                       R^T R BV + BV R^T R = F
    % for BU and BV.
    %
    % This can be solved using two calls to the Matlab's lyap.
    % However, we can still have a more efficient implementation
    % that does not require the full functionaliyt of Matlab's lyap.

    [U, Sigma, V] = svd(R);
    E_mod = U'*E*U;
    F_mod = V'*F*V;
    b1 = E_mod(:);
    b2 = F_mod(:);

    r = size(Sigma, 1);
    sig = diag(Sigma); % all the singular values in a vector
    sig1 = sig*ones(1, r); % columns repeat
    sig1t = sig1'; % rows repeat
    s1 = sig1(:);
    s2 = sig1t(:);

    % The block elements
    a =  s1.^2 + s2.^2; % a column vector

    % Solve the linear system of equations
    cu = b1./a; %a.\b1;
    cv = b2./a; %a.\b2;

    % Matricize
    CU = reshape(cu, r, r);
    CV = reshape(cv, r, r);

    % Do the similarity transforms
    BU = U*CU*U';
    BV = V*CV*V';

    % %% Debug
    %
    % norm(R*R'*BU + BU*R*R' - E, 'fro');
    % norm((Sigma.^2)*CU + CU*(Sigma.^2) - E_mod, 'fro');
    % norm(a.*cu - b1, 'fro');
    %
    % norm(R'*R*BV + BV*R'*R - F, 'fro');
    %
    % BU1 = lyap(R*R', - E);
    % norm(R*R'*BU1 + BU1*R*R' - E, 'fro');
    %
    % BV1 = lyap(R'*R, - F);
    % norm(R'*R*BV1 + BV1*R'*R - F, 'fro');
    %
    % % as accurate as the lyap
    % norm(BU - BU1, 'fro')
    % norm(BV - BV1, 'fro')
end

function[Omega1, Omega2] = coupled_lyap(R, E, F)
    % We intent to solve the coupled system of Lyapunov equations
    %
    % RR^T Omega1 + Omega1 RR^T  - R Omega2 R^T = E
    % R^T R Omega2 + Omega1 R^T R  - R^T Omega2 R = F,
    %
    % for Omega1 and Omega2, both are skew symmetric matrices.
    %
    % Below is an efficient implementation

    [U, Sigma, V] = svd(R);
    E_mod = U'*E*U;
    F_mod = V'*F*V;
    b1 = E_mod(:);
    b2 = F_mod(:);

    r = size(Sigma, 1);
    sig = diag(Sigma); % All the singular values in a vector
    sig1 = sig*ones(1, r); % Columns repeat
    sig1t = sig1'; % Rows repeat
    s1 = sig1(:);
    s2 = sig1t(:);

    % The block elements
    a =  s1.^2 + s2.^2; % A column vector
    c = s1.*s2;

    % Solve directly using the formula
    % A = diag(a);
    % C = diag(c);
    % Y1_sol = (A*(C\A) - C) \ (b2 + A*(C\b1));
    % Y2_sol = A\(b2 + C*Y1_sol);

    Y1_sol = (b2 + (a./c).*b1) ./ ((a.^2)./c - c);
    Y2_sol = (b2 + c.*Y1_sol)./a;

    % Matricize
    Omega1 = reshape(Y1_sol, r, r);
    Omega2 = reshape(Y2_sol, r, r);

    % Do the similarity transforms
    Omega1 = U*Omega1*U';
    Omega2 = V*Omega2*V';

    % %% Debug: whether we have the right solution.
    % norm(R*R'*Omega1 + Omega1*R*R'  - R*Omega2*R' - E, 'fro')
    % norm(R'*R*Omega2 + Omega2*R'*R  - R'*Omega1*R - F, 'fro')
end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [f, fval] = test(varargin)

% Images
if nargin == 0
    A = imread('pout.tif');
    B = imread('cameraman.tif');
elseif nargin == 2
    A = varargin{1};
    B = varargin{2};
    if ischar(A)
        A = imread(A);
    end;
    if ischar(B)
        B = imread(B);
    end;
end;

% Histograms
nbins = 10;
[ca ha] = imhist(A, nbins);
[cb hb] = imhist(B, nbins);

% Features
f1 = ha;
f2 = hb;

% Weights
w1 = ca / sum(ca);
w2 = cb / sum(cb);

% Earth Mover's Distance
[f, fval] = emd(f1, f2, w1, w2, @gdf);

% Results
wtext = sprintf('fval = %f', fval);
figure('Name', wtext);
subplot(121);imshow(A);title('first image');
subplot(122);imshow(B);title('second image');

end

author: Sable

repo: mcbench-benchmarks


text:
function [CX, sse] = vgg_kmeans(X, nclus, varargin)

% VGG_KMEANS    initialize K-means clustering
%               [CX, sse] = vgg_kmeans(X, nclus, optname, optval, ...)
%
%               - X: input points (one per column)
%               - nclus: number of clusters
%               - opts (defaults):
%                    maxiters (inf): maxmimum number of iterations
%                    mindelta (eps): minimum change in SSE per iteration
%                       verbose (1): 1=print progress
%
%               - CX: cluster centers
%               - sse: SSE

% Author: Mark Everingham <me@robots.ox.ac.uk>
% Date: 13 Jan 03

opts = struct('maxiters', inf, 'mindelta', eps, 'verbose', 1);
if nargin > 2
    opts=vgg_argparse(opts,varargin);
end

perm=randperm(size(X,2));
CX=X(:,perm(1:nclus));

sse0 = inf;
iter = 0;
while iter < opts.maxiters

    tic;
    [CX, sse] = vgg_kmiter(X, CX);
    t=toc;

    if opts.verbose
        fprintf('iter %d: sse = %g (%g secs)\n', iter, sse, t)
    end

    if sse0-sse < opts.mindelta
        break
    end

    sse0=sse;
    iter=iter+1;

end

author: rksltnl


text:
function [ x, seed ] = sphere_unit_sample_2d ( seed )

%*****************************************************************************80
%
%% SPHERE_UNIT_SAMPLE_2D picks a random point on the unit sphere (circle) in 2D.
%
%  Discussion:
%
%    The unit sphere in 2D satisfies:
%
%      X * X + Y * Y = 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(2), a random point on the unit circle.
%
%    Output, integer SEED, a seed for the random number generator.
%
  dim_num = 2;

  [ u, seed ] = r8_uniform_01 ( seed );

  x(1) = cos ( 2.0 * pi * u );
  x(2) = sin ( 2.0 * pi * u );

  return
end


text:
function v = get_path(arg)

if strcmp(arg, 'tri')
	v = get_path_tri();
end

if strcmp(arg, 'sep')
	v = get_path_sep();
end

end

repo: toflow


text:
% demo code for 3D human pose estimation from a monocular RGB image
% this demo code uses CPM (Convolutional Pose Machines) for 2D pose
% estimation, however, any 2D pose estimator from image can be applied.
% prediction: The 2D pose, a 14X2 matrix, 14 corresponding position as
% listed in the following order, 3 columns with [x y] coordinates.
% (head)
% (neck)
% (right shoulder)
% (right elbow)
% (right wrist)
% (left shoulder)
% (left elbow)
% (left wrist)
% (right hip)
% (right knee)
% (right ankle)
% (left hip)
% (left knee)
% (left ankle)

%% setup
% param = model_config();
load('3D_library.mat');

%% input image path
test_image = ['imgs/im0371.jpg'];
im = imread(test_image);

% load the precomputed 2D pose derived by Convolutional Pose Machines,
% other 2D pose estimation methods is applicable, format the 2D pose in the
% definition described above
load('0371.mat');

%{
% get image size
[H,W,~] = size(im);

%% detect human using CPM
rectangle(1) = 1;
rectangle(2) = 1;
rectangle(3) = W-1;
rectangle(4) = H-1;
[heatMaps, human_c] = applyModel_human3d(test_image, param, rectangle, 2);

%% 2D pose estimation using CPM with the bounding box specified by the
% detected human postition
rectangle(1) = human_c(1) - 0.5*W;
rectangle(2) = human_c(2) - 0.5*H;
rectangle(3) = 1*W;
rectangle(4) = 1*H;
[heatMaps, prediction] = applyModel_human3d(test_image, param, rectangle, 1);
%}
Prediction{1}=prediction;

% manually adjust two hips' position due to the variation across MPII and
% H36M
prediction(9,2) = prediction(9,2) - 20;
prediction(12,2) = prediction(12,2) - 20;

%% extract the nearest neighbor
[j_p] = NN_pose(s1_s9_2d_n,s1_s9_3d,prediction);
% [j_p] = kNN_pose_procrus(s1_s9_2d_n,s1_s9_3d,prediction, 10);

% compute the scale between pixel and real world to recover real size of
% prediction
scale = (max(j_p(:,2))-min(j_p(:,2)))/(max(prediction(:,2))-min(prediction(:,2)));

% predict the depth of each joint by the exemplar
prediction(:,3) = j_p(:,3)/scale;

%% visualization
H=figure();
subplot(1,3,[1 2]);imshow(im);
% self-defined ground plane
prediction(:,3) = prediction(:,3) - min(prediction(:,3)) + 1300;
an_x_m = min(prediction(:,1));
an_x_M = max(prediction(:,1));
an_y_M = max(prediction(11,2),prediction(14,2));
x = [an_x_m-100 an_x_M+100 an_x_m-100 an_x_M+100];
y = [an_y_M-10 an_y_M-10 an_y_M+40 an_y_M+40];
z = [min(prediction(:,3)-120) min(prediction(:,3)-120) max(prediction(:,3)+120) max(prediction(:,3)+120)];
A = [x(:) y(:) z(:)];
% solve for the ground plane
[n,v,m,aved]=plane_fit(A);

% draw the prediction on the input image
vis_2d(prediction);

% draw the 3D pose
subplot(1,3,3);
vis_3d(prediction);

% draw the ground plane
planeplot(A,n,m)

% draw camera
scale = 50;
P = scale*[0 0 0;0.5 0.5 0.8; 0.5 -0.5 0.8; -0.5 0.5 0.8;-0.5 -0.5 0.8];
cen = mean(prediction);P1=(P+repmat([cen(1:2), 800],[5,1]));
line([P1(1,1) P1(2,1)],[P1(1,2) P1(2,2)],[P1(1,3) P1(2,3)],'color','k')
line([P1(1,1) P1(3,1)],[P1(1,2) P1(3,2)],[P1(1,3) P1(3,3)],'color','k')
line([P1(1,1) P1(4,1)],[P1(1,2) P1(4,2)],[P1(1,3) P1(4,3)],'color','k')
line([P1(1,1) P1(5,1)],[P1(1,2) P1(5,2)],[P1(1,3) P1(5,3)],'color','k')

line([P1(2,1) P1(3,1)],[P1(2,2) P1(3,2)],[P1(2,3) P1(3,3)],'color','k')
line([P1(3,1) P1(5,1)],[P1(3,2) P1(5,2)],[P1(3,3) P1(5,3)],'color','k')
line([P1(5,1) P1(4,1)],[P1(5,2) P1(4,2)],[P1(5,3) P1(4,3)],'color','k')
line([P1(4,1) P1(2,1)],[P1(4,2) P1(2,2)],[P1(4,3) P1(2,3)],'color','k')

% adjust the viewing angle
view(26,-56);
axis equal
axis off

author: flyawaychase


text:
% Estimation of the standard deviation of an image assuming that the noise
% is Gaussian white additive. The file estimates the standard deviation
% using the median filter on the fine scale subband of the wavelet
% decomposition.
 function sd_estimate=sdest(x)
 [ca,ch,cv,cd] = dwt2(x,'sym4','mode','sym');
 sd_estimate = mad(cd(:),1)/0.6745

author: thomaskuestner


text:
%IRECTIFY Rectify stereo image pair
%
% [OUT1,OUT2] = IRECTIFY(F, M, IM1, IM2) is a rectified pair of images
% corresponding to IM1 and IM2.  F (3x3) is the fundamental matrix relating
% the two views and M is a FeatureMatch object containing point correspondences
% between the images.
%
% [OUT1,OUT2,H1,H2] = IRECTIFY(F, M, IM1, IM2) as above but also returns
% the homographies H1 and H2 that warp IM1 to OUT1 and IM2 to OUT2 respectively.
%
% Notes::
% - The resulting image pair are epipolar aligned, equivalent to the view
%   if the two original camera axes were parallel.
% - Rectified images are required for dense stereo matching.
% - The effect of lense distortion is not removed, use the camera calibration
%   toolbox to unwarp each image prior to rectification.
% - The resulting images may have negative disparity.
% - Some output pixels may have no corresponding input pixels and will be
%   set to NaN.
%
% See also FeatureMatch, ISTEREO, HOMWARP, CentralCamera.

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

function [Img1_new, Img2_new, H12,H21] = irectify(F, m, Img1, Img2)
% http://se.cs.ait.ac.th/cvwiki/matlab:tutorial:rectification

F12 = F';

[rows,cols,depth] = size(Img1);

% Get homographies.

x1 = e2h( m.p1 );
x2 = e2h( m.p2 );

[H12,H21,bSwap] = rectify_homographies( F12, x1, x2, rows, cols );

[w1,off1] = homwarp(H12, Img1, 'full');
[w2,off2] = homwarp(H21, Img2, 'full');

% fix the vertical alignment of the images by padding
dy = off1(2) - off2(2);
if dy < 0
    w1 = ipad(w1, 'b', -dy);
    w2 = ipad(w2, 't', -dy);
else
    w1 = ipad(w1, 't', dy);
    w2 = ipad(w2, 'b', dy);
end

[w1,w2] = itrim(w1, w2);

if nargout == 0
    stdisp(w1, w2)
else
    Img1_new = w1;
    Img2_new = w2;
end

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

function [H1,H2,bSwap] = rectify_homographies( F, x1, x2, rows, cols )

  % F: a fundamental matrix

  % x1 and x2: corresponding points such that x1_i' * F * x2_i = 0

  % Initialize

  H1 = [];
  H2 = [];
  bSwap = 0;

  % Center of image

  cy = round( rows/2 );
  cx = round( cols/2 );

  % Fix F to be rank 2 to numerical accuracy

  [U,D,V] = svd( F );
  D(3,3) = 0;
  F = U*D*V';

  % Get epipole.  e12 is the epipole in image 1 for camera 2.

  e12 = null( F' );             % Epipole in image 1 for camera 2
  e21 = null( F );              % Epipole in image 2 for camera 1

  % Put epipoles in front of camera

  if e12 < 0, e12 = -e12; end;
  if e21 < 0, e21 = -e21; end;

  % Make sure the epipoles are inside the images

  check_epipoles_in_image( e12, e21, rows, cols );

  % Check that image 1 is to the left of image 2

%   if e12(1)/e12(3) < cx
%     fprintf( 1, 'Swapping left and right images...\n' );
%     tmp = e12;
%     e12 = e21;
%     e21 = tmp;
%     F = F';
%     bSwap = 1;
%   end;

  % Now we have
  % F' * e12 = 0,
  % F  * e21 = 0,

  % Let's get the rectifying homography Hprime for image 1 first

  Hprime = map_to_infinity( e12, cx, cy );
  e12_new = Hprime * e12;
  % Normalize Hprime so that Hprime*eprime = (1,0,0)'
  Hprime = Hprime / e12_new(1);
  e12_new = Hprime * e12;
  fprintf( 1, 'Epipole 1/2 mapped to infinity: (%g, %g, %g)\n', e12_new );

  % Get canonical camera matrices for F12 and compute H0, one possible
  % rectification homography for image 2

  [P,Pprime] = get_canonical_cameras( F );
  M = Pprime(:,1:3);
  H0 = Hprime * M;

  % Test that F12 is a valid F for P,Pprime

  test_p_f( P, Pprime, F );

  % Now we need to find H so that the epipolar lines match
  % each other, i.e., inv(H)' * l = inv(Hprime)' * lprime
  % and the disparity is minimized, i.e.,
  % min \sum_i d(H x_i, Hprime xprime_i)^2

  % Transform data initially according to Hprime (img 1) and H0 (img 2)

  x1hat = Hprime * x1;
  x1hat = x1hat ./ repmat( x1hat(3,:), 3, 1 );
  x2hat = H0 * x2;
  x2hat = x2hat ./ repmat( x2hat(3,:), 3, 1 );
  rmse_x = sqrt( mean( (x1hat(1,:) - x2hat(1,:) ).^2 ));
  rmse_y = sqrt( mean( (x1hat(2,:) - x2hat(2,:) ).^2 ));
  fprintf( 1, 'Before Ha, RMSE for corresponding points in Y: %g X: %g\n', ...
           rmse_y, rmse_x );

  % Estimate [ a b c ; 0 1 0 ; 0 0 1 ] aligning H, Hprime

  n = size(x1,2);
  A = [ x2hat(1,:)', x2hat(2,:)', ones(n,1) ];
  b = x1hat(1,:)';
  abc = A\b;
  HA = [ abc' ; 0 1 0 ; 0 0 1 ];
  H = HA*H0;
  x2hat = H * x2;
  x2hat = x2hat ./ repmat( x2hat(3,:), 3, 1 );
  rmse_x = sqrt( mean(( x1hat(1,:) - x2hat(1,:) ).^2 ));
  rmse_y = sqrt( mean(( x1hat(2,:) - x2hat(2,:) ).^2 ));
  fprintf( 1, 'After Ha, RMSE for corresponding points in Y: %g X: %g\n', ...
           rmse_y, rmse_x );

  % Return the homographies as appropriate

  if bSwap
    H1 = H;
    H2 = Hprime;
  else
    H1 = Hprime;
    H2 = H;
  end;

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

function check_epipoles_in_image( e1, e2, rows, cols )

  % Check whether given epipoles are in the image or not

  if abs( e1(3) ) < 1e-6 & abs( e2(3) ) < 1e-6, return; end;

  e1 = e1 / e1(3);
  e2 = e2 / e2(3);
  if ( e1(1) <= cols & e1(1) >= 1 & e1(2) <= rows & e1(2) >= 1 ) | ...
     ( e2(1) <= cols & e2(1) >= 1 & e2(2) <= rows & e2(2) >= 1 )
    err_msg = sprintf( 'epipole (%g,%g) or (%g,%g) is inside image', ...
                       e1(1:2), e2(1:2) );
    error( [ err_msg, ' -- homography does not work in this case!' ] );
  end;

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

function [P,Pprime] = get_canonical_cameras( F )

  % Get the "canonical" cameras for given fundamental matrix
  % according to Hartley and Zisserman (2004), p256, Result 9.14

  % But ensure that the left 3x3 submatrix of Pprime is nonsingular
  % using Result 9.15, that the general form is
  % [ skewsym( e12 ) * F + e12 * v', k * e12 ] where v is an arbitrary
  % 3-vector and k is an arbitrary scalar

  P = [ 1 0 0 0
        0 1 0 0
        0 0 1 0 ];

  e12 = null( F' );
  M = skew( e12 ) * F + e12 * [1 1 1];
  Pprime = [ M, e12 ];

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

function test_p_f( P, Pprime, F )

  % Test that camera matrices Pprime and P are consistent with
  % fundamental matrix F
  % Meaning  (Pprime*X)' * F * (P*X) = 0,  for all X in 3space

  % Get the epipole in camera 1 for camera 2

  C2 = null( P );
  eprime = Pprime * C2;

  % Construct F from Pprime, P, and eprime

  Fhat = skew( eprime ) * Pprime * pinv( P );

  % Check that it's close to F

  alpha = Fhat(:)\F(:);
  if norm( alpha*Fhat-F ) > 1e-10
    fprintf( 1, 'Warning: supplied camera matrices are inconsistent with F\n' );
  else
    fprintf( 1, 'Supplied camera matrices OK\n' );
  end;

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

function H = map_to_infinity( x, cx, cy )

  % Given a point and the desired origin (point of minimum projective
  % distortion), compute a homograph H = G*R*T taking the point to the
  % origin, rotating it to align with the X axis, then mapping it to
  % infinity.

  % First map cx,cy to the origin

  T = [ 1 0 -cx
        0 1 -cy
        0 0 1 ];
  x = T * x;

  % Now rotate the translated x to align with the X axis.

  cur_angle = atan2( x(2), x(1) );
  R = [ cos( -cur_angle ), -sin( -cur_angle ), 0
        sin( -cur_angle ),  cos( -cur_angle ), 0
                        0,                  0, 1 ];
  x = R * x;

  % Now the transformation G mapping x to infinity

  if abs( x(3)/norm(x) ) < 1e-6
      % It's already at infinity
      G = eye(3)
  else
      f = x(1)/x(3);
      G = [    1   0  0
               0   1  0
             -1/f  0  1 ];
  end;

  H = G*R*T;

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [x, timing] = blendenpik_under_alternative(A, b, params)
% [x, timing] = blendenpik_under(A, b, params)
%
% Alternative function for solving the equation
% x = arg min norm(A * x - b, 2) using Blendenpik, where A has
% more columns than rows. Less stable than the regular version.
%
% "params" - parameters governning the method.
%    params.type - type of mixing transform. Optional values: 'DCT', 'DHT', 'WHT'.
%                  Default is DHT.
%    params.gamma - gamma * m columns will be sampled (A is m-by-n).
%                   Default is 4.
%    params.preprocess_steps - number of mixing steps to do in advance.
%                               Default is 1.
%    params.maxcond - maximum condition number of the preconditioner.
%                     Default is 1 / (5 * epsilon_machine).
%    params.tol - convergence thershold for LSQR.
%                 Default is 1e-14.
%    params.maxit - maximum number of LSQR iterations.
%                   Default is 1000.
%    params.lsvec - whether to output in "timing" the LSQR residuals.
%                   Default is false.
%    params.use_full_lsqr - whether to use LSQR with full
%                           orthogonalization. Useful for
%                           preprocess_steps=0.
%                           Default is false.
%
% Output:
%   x - the solution.
%   timing - statistics on the time spent on various phases.
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.

if (nargin < 3)
    params = struct;
end

if (~isfield(params, 'type'));
    params.type = 'DHT';
end

if (~isfield(params, 'gamma'))
    params.gamma = 4;
end

if (~isfield(params, 'preprocess_steps'))
    params.preprocess_steps = 1;
end

if (~isfield(params, 'tol'))
    params.tol = 1e-14;
end

if (~isfield(params, 'maxit'))
    params.maxit = 1000;
end

if (~isfield(params, 'maxcond'))
    params.maxcond = 1 / (5 * eps);
end

if (~isfield(params, 'lsvec'))
    params.lsvec = false;
end

if (~isfield(params, 'slight_coherence'))
    params.slight_coherence = 0;
end

if (~isfield(params, 'use_full_lsqr'))
    params.use_full_lsqr = false;
end

tstart = wtime;

%% Build preconditioner
t1 = wtime;

[m, n] = size(A);
AT = A';
t0 = wtime;
frut_D = sign(randn(n, 1));
B = fast_unitary_transform(AT, frut_D, params.type);
nn = size(B, 1);
frut_D(n+1:nn) = 0;
frut_time = wtime - t0;
disp(sprintf('\t\tRandom unit diagonal + unitary transformation time: %.2f sec', frut_time));
timing.precond_timing.timing.frut_time = frut_time;

t0 = wtime;
t = params.gamma * m / nn;
s = rand(nn, 1);
SB = B(find(s < t), :);
sample_time = wtime - t0;
disp(sprintf('\t\tRandom sampling time: %.2f sec', sample_time));
timing.precond_timing.sample_time = sample_time;

t0 = wtime;
[Y, tau] = mex_dgeqrf(SB);
qr_time = wtime - t0;
disp(sprintf('\t\tQR on random sample time: %.2f sec', qr_time));
timing.qr_time = qr_time;

timing.preprocess_total_time = wtime - t1;
disp(sprintf('\tPreprocessing time: %.2f sec', timing.preprocess_total_time));

%% One solution
t1 = wtime;
y = mex_dtrsm(Y, b, 1.0, 'L', 'U', 'T', 'N');
ssz = size(Y, 1);
p0 = mex_dormqr(Y, tau, [y; zeros(ssz-m, 1)], 'L', 'N');
p1 = zeros(nn, 1); p1(s < t) = p0;
p2 = fast_unitary_transform(p1, frut_D, ['I' params.type]);
p = p2(1:n);
timing.find_non_min_time = wtime - t1;
disp(sprintf('\tFind a non-minimal solution: %.2f sec', timing.find_non_min_time));

%% LSQR
t1 = wtime;
R = triu(Y(1:m, 1:m));
% clear Y;
% clear SB
% clear tau;
if (~params.lsvec)
    if (~params.use_full_lsqr)
        [x0, timing.lsqr_its] = dense_lsqr(AT, p, R, params.tol, params.maxit);
    else
        [x0, timing.lsqr_its] = dense_full_lsqr(AT, p, R, params.tol, params.maxit);
    end
else
    [x0, timing.lsqr_its, timing.lsvec, timing.resvec] = dense_lsqr(AT, p, R, params.tol, params.maxit);
end
timing.lsqr_time = wtime - t1;
disp(sprintf('\tLSQR time: %.2f sec', timing.lsqr_time));

%% Mult by AT to get final solution
t1 = wtime;
x = AT * x0;
timing.multAT_time = wtime - t1;
disp(sprintf('\tMultiply by A'' time: %.2f sec', timing.multAT_time));
timing.total_time = wtime - tstart;
disp(sprintf('Total time: %.2f sec', timing.total_time));


text:
function labels_plot(expdirs,jabs,nframesplot,varargin)
% Labelgram

% AL 20150128: This is looking more like ethogram_plot
% AL 20150202: even worse, xlsExport

if ~exist('expdirs','var')
  expdirs = {};
end
assert(iscellstr(expdirs));
if ~exist('jabs','var')
  jabs = [];
end
if isempty(jabs)
  [tfsuccess,jabs] = ExpPP.uiGetJabFiles;
  if ~tfsuccess
    return;
  end
end
if ~exist('nframesplot','var')
  nframesplot = [];
elseif ~isempty(nframesplot)
  assert(isnumeric(nframesplot));
  assert(isscalar(nframesplot) || numel(nframesplot)==2,'Invalid nframesplot');
end

[xlsexport,xlsfile,ethcoreargexport] = myparse(varargin,...
  'xlsexport',false,...
  'xlsfile',[],...
  'ethcoreargexport',false);

assert(isscalar(jabs),'Currently expect single jab/macguffin.');
if iscellstr(jabs)
  jabfiles = jabs;
  jabs = loadAnonymous(jabs{1});
elseif isa(jabs,'Macguffin');
  jabfiles = {[]};
else
  assert(false,'''jabs'' argument must be either a jab filename or object.');
end

if isempty(expdirs)
  expdirs = jabs.expDirNames;
end
Nexp = numel(expdirs);

jab = jabs;
jab.modernize();

% Special case for AH exps. Call ExpPP.loadexps on all experiments now.
% Any experiments that do not load are excluded.
behs = Labels.verifyBehaviorNames(jab.behaviors.names);
doautomarks = ExpPP.doNamesSpanAllBasicBehaviors(behs);
if doautomarks
  % guess in this case require jabfiles to be nonempty, ie jabs cannot have
  % been obj
  dPP = ExpPP.loadexps(expdirs,jabfiles,'includelabels',true);
  if isempty(dPP)
    error('labels_plot:noExps','ExpPP failed to load any experiments.');
  end
  expsLoaded = {dPP.expfull}';
  tfloaded = ismember(expdirs,expsLoaded);
  if ~all(tfloaded)
    tmpstr = expdirs(~tfloaded);
    tmpstr = sprintf('%s\n',tmpstr{:});
    warningNoTrace('labels_plot:expsExcluded','Excluding %d incomplete experiments from labelgram: %s.',...
      nnz(~tfloaded),tmpstr);
  end
  expdirs = expsLoaded;
  Nexp = numel(expdirs);
end

% get the labels for expdirs
assert(numel(jab.expDirNames)==numel(jab.labels));
[tf,loc] = ismember(expdirs,jab.expDirNames);
lbls = Labels.labels(Nexp);
lbls(tf) = jab.labels(loc(tf));

% Trx/T0/T1/nframesplot
trxfname = jab.file.trxfilename;
trxstarts = nan(Nexp,1);
trxends = nan(Nexp,1);
for i = 1:Nexp
  trxfile = fullfile(expdirs{i},trxfname);
  tmp = load(trxfile,'trx');
  trxstarts(i) = min([tmp.trx.firstframe]);
  trxends(i) = max([tmp.trx.endframe]);
end
T0 = min(trxstarts);
T1 = max(trxends)-1;
if ~isempty(nframesplot)
  if numel(nframesplot)==1,
    T1 = min(T1,T0+nframesplot-1);
  else
    assert(numel(nframesplot)==2);
    T0 = nframesplot(1);
    T1 = min(T1,nframesplot(2));
  end
end

if doautomarks
  % in this case:
  % - use exptags from jab if possible

  dPP = ExpPP.ensureCustomGroupInit(dPP);
  assert(isequal({dPP.expfull}',expdirs));
  [tfExpInJab,loc] = ismember(expdirs,jab.expDirNames);
  if all(tfExpInJab)
    if isstruct(jab) && isfield(jab,'expDirTags') || ...
       isobject(jab) && isprop(jab,'expDirTags')
      exptags = jab.expDirTags(loc);
      exptags = ExperimentTags.cleanLegacyNotTags(exptags);
    else
      exptags = [];
    end
    if isempty(exptags)
      exptags = ExperimentTags.expTags(expdirs);
    end

    allUniqueTags = ExperimentTags.allUniqueTags(exptags);
    allUniqueTagsValid = matlab.lang.makeValidName(allUniqueTags);
    if numel(unique(allUniqueTagsValid))~=numel(allUniqueTagsValid)
      error('ethogram_plot:tag','Repeated tags.'); % very rare
    end
    Ntag = numel(allUniqueTags);
    for iTag = 1:Ntag
      tagRaw = allUniqueTags{iTag};
      cgrpname = allUniqueTagsValid{iTag};
      tfTag = ExperimentTags.findTag(exptags,tagRaw);
      dPP = ExpPP.createCustomGroup(dPP,tfTag,cgrpname);
    end
  end

  ethplotcoreargs = {'doautomarks',true,'automarkdata',dPP,'expppstatprefix','labl','expppcontrolargs',{'labelgram' true}};
else
  ethplotcoreargs = cell(1,0);
end

% go
m = Labels.labelMatrix(lbls,T0,T1,behs);
ethplotcoreargs = [{figure,m,T0,T1,expdirs,behs,repmat(jabfiles,numel(behs),1)} ethplotcoreargs];
if ethcoreargexport
  fname = sprintf('ethcoreargs.%s.mat',datestr(now,'yyyymmddTHHMMSS'));
  save(fname,'ethplotcoreargs');
end
[boutmat,line_names] = ethogram_plot_core(ethplotcoreargs{:});

if xlsexport
  if isempty(xlsfile)
    [xlsfile,xlspath] = uiputfile('bouts.txt','Export to file');
    if isequal(xlsfile,0)
      xlsfile = [];
    else
      xlsfile = fullfile(xlspath,xlsfile);
    end
  elseif exist(xlsfile,'file')
    warning('ethogram_plot:fileExists','Overwriting file ''%s''.',xlsfile);
  end

  if ~isempty(xlsfile)
    Ethogram.boutExport(xlsfile,boutmat,line_names,behs);
    fprintf(1,'Exported to tab-delimited file ''%s''.\n',xlsfile);
  end
end


text:
% STD_ERPIMAGEPLOT - Commandline function to plot cluster ERPimage or channel erpimage.
%
% Usage:
%   >> [STUDY] = std_erpimageplot(STUDY, ALLEEG, key1, val1, key2, val2);
%   >> [STUDY data times freqs pgroup pcond pinter] = ...
%                std_erpimageplot(STUDY, ALLEEG ...);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the datasets in the STUDY.
%                Note: ALLEEG for a STUDY set is typically created using LOAD_ALLEEG.
%
% Additional help:
% Inputs and output of this function are strictly identical to the STD_ERSPPLOT.
% See the help message of this function for more information.
%
% See also: STD_ERSPPLOT
%
% Authors: Arnaud Delorme, UCSD/CERCO, August, 2011-

% Copyright (C) Arnaud Delorme, arno@ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [STUDY, allitc, alltimes, allfreqs, pgroup, pcond, pinter, events] = std_erpimageplot(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_erpimageplot;
    return;
end

events = [];
STUDY = pop_erpimparams(STUDY, 'default');
if strcmpi(STUDY.etc.erpimparams.concatenate, 'off')

    % use std_erspplot for stats when concatenate is off
    [STUDY allitc alltimes allfreqs pgroup pcond pinter events] = std_erspplot(STUDY, ALLEEG, 'datatype', 'erpim', varargin{:});

else
    params = STUDY.etc.erpimparams;
    [ opt moreparams ] = finputcheck( varargin, { ...
        'design'      'integer' [] STUDY.currentdesign;
        'topotime'    'real'    [] params.topotime;
        'topotrial'  'real'     [] params.topotrial;
        'timerange'   'real'    [] params.timerange;
        'trialrange'  'real'    [] params.trialrange;
        'colorlimits' 'real'    [] params.colorlimits; % ERPimage
        'statistics'  'string'  [] params.statistics;
        'groupstats'  'string'  [] params.groupstats;
        'condstats'   'string'  [] params.condstats;
        'threshold'   'real'    [] params.threshold;
        'naccu'       'integer' [] params.naccu;
        'mcorrect'    'string'  [] params.mcorrect;
        'erpimageopt' 'cell'    [] params.erpimageopt;
        'concatenate'  'string' { 'on','off' }  params.concatenate;
        'channels'    'cell'    []              {};
        'clusters'    'integer' []              [];
        'comps'       {'integer','string'}  []              []; % for backward compatibility
        'plotsubjects' 'string' { 'on','off' }  'off';
        'plotmode'    'string' { 'normal','condensed','none' }  'normal';
        'subject'     'string'  []              '' }, 'std_erpimageplot', 'ignore');

    if ~isempty(opt.topotime) && ~isempty(opt.topotrial)
        error('Cannot plot topography when ERP-image is in trial concatenation mode');
    end
    if ~isempty(opt.trialrange)
        error('Cannot select trial range when ERP-image is in trial concatenation mode');
    end
    if strcmpi(opt.groupstats, 'on') || strcmpi(opt.condstats, 'on')
        disp('Warning: cannot perform statistics when ERP-image is in trial concatenation mode');
    end

    % options
    if ~isempty(opt.colorlimits), options = { 'caxis' opt.colorlimits opt.erpimageopt{:} };
    else                          options = { 'cbar' 'on' opt.erpimageopt{:}  };
    end

    if ~isempty(opt.channels)
        [STUDY allerpimage alltimes alltrials tmp events] = std_readersp(STUDY, ALLEEG, 'channels', opt.channels, 'infotype', 'erpim', 'subject', opt.subject, ...
            'concatenate', 'on', 'timerange', opt.timerange, 'design', opt.design);

        % get figure title
        % ----------------
        locs = eeg_mergelocs(ALLEEG.chanlocs);
        locs = locs(std_chaninds(STUDY, opt.channels));
        allconditions = STUDY.design(opt.design).variable(1).value;
        allgroups     = STUDY.design(opt.design).variable(2).value;
        alltitles = std_figtitle('condnames', allconditions, 'cond2names', allgroups, 'chanlabels', { locs.labels }, ...
            'subject', opt.subject, 'valsunit', 'ms', 'datatype', 'ERPIM');

        figure;
        for iCond = 1:length(allconditions)
            for iGroup = 1:length(allgroups)
                tmpevents = events{iCond, iGroup};
                if isempty(tmpevents), tmpevents = zeros(1, size(allerpimage{iCond, iGroup},2)); end
                subplot(length(allconditions), length(allgroups), (iCond-1)*length(allgroups) + iGroup);

                % use color scale for last plot
                if ~isempty(opt.colorlimits) && iCond == length(allconditions) && iGroup == length(allgroups)
                    options = { options{:} 'cbar' 'on' };
                end
                erpimage(allerpimage{iCond, iGroup}, tmpevents, alltimes, alltitles{iCond, iGroup}, params.smoothing, params.nlines, options{:});
            end
        end

    else
        for cInd = 1:length(opt.clusters)
            [STUDY allerpimage alltimes alltrials tmp events] = std_readersp(STUDY, ALLEEG, 'clusters', opt.clusters(cInd), 'infotype', 'erpim', 'subject', opt.subject, ...
                'concatenate', 'on', 'timerange', opt.timerange, 'design', opt.design);
            % get figure title
            % ----------------
            locs = eeg_mergelocs(ALLEEG.chanlocs);
            locs = locs(std_chaninds(STUDY, opt.channels));
            allconditions = STUDY.design(opt.design).variable(1).value;
            allgroups     = STUDY.design(opt.design).variable(2).value;
            alltitles = std_figtitle('condnames', allconditions, 'cond2names', allgroups, 'clustname', STUDY.cluster(opt.clusters(cInd)).name, ...
                'subject', opt.subject, 'valsunit', 'ms', 'datatype', 'ERPIM');

            figure;
            for iCond = 1:length(allconditions)
                for iGroup = 1:length(allgroups)
                    tmpevents = events{iCond, iGroup};
                    if isempty(tmpevents), tmpevents = zeros(1, size(allerpimage{iCond, iGroup},2)); end
                    subplot(length(allconditions), length(allgroups), (iCond-1)*length(allgroups) + iGroup);

                    % use color scale for last plot
                    if ~isempty(opt.colorlimits) && iCond == length(allconditions) && iGroup == length(allgroups)
                        options = { options{:} 'cbar' 'on' };
                    end
                    erpimage(allerpimage{iCond, iGroup}, tmpevents, alltimes, alltitles{iCond, iGroup}, params.smoothing, params.nlines, options{:});
                end
            end
        end
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function beta = volumeCompressibility(S)
% computes the volume compressibility of an elasticity tensor
%
% Syntax
%   beta = volumeCompressibility(S)
%
% Input
%  C - elastic stiffness @tensor
%
% Output
%  beta - volume compressibility
%
% Description
%
% $$\beta(x) = S_{iikk}$$
%

% compute tensor product
beta = EinsteinSum(S,[-1 -1 -2 -2]);

% beta = 1/S.bulkModulus(S)

author: mtex-toolbox

repo: mtex

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


text:
function a=getNextMPartition(param1,m)
%%GETNEXTMPARTITION Get the next partition of the integer n into m
%                   parts. This is m nonzero items such that the sum is n.
%                   The partitions are visited in reverse lexicographic
%                   order (colex order). The function can be called with
%                   two parameters to get the first n,m partition, or with
%                   1 parameter (that being the previous partition) to
%                   get the next partition. The total number of
%                   partitions can be found using the function
%                   numMPartitions(n,m). As combinations are to
%                   permutations, partitions are to compositions: The order
%                   of the terms in an m-partitions do not matter, whereas
%                   they matter in a composition.
%
%INPUTS: param1 If two parameters are passed to the function, then param1
%            is n, the integer to be partitioned. Otherwise, param1 is a,
%            the current partition of the chosen integer n into m parts.
%            This means that a must be length m and its elements must sum
%            to n. The elements must be ordered in decreasing order. The
%            first partition in the series must be a(1)=n-m+1; a(2:m)=1;
%            Note that n>=m>=2 if a is passed, because a size 1 partition
%            only has one value (the initial value).
%          m This parameter is not passed if the next m-partition in the
%            sequence is desired. However, if one wishes to obtain the
%            first value in the sequence, then m is the positive integer
%            number of parts into which the integer n is split, n>=m;
%
%OUTPUTS: a The next partition in the series or the empty matrix if the
%           parition given was the last partition in the series. If the
%           function was called getNextMPartition(n,m), then a is the first
%           parition.
%
%The algorithm is Algorithm H in Chapter 7.2.1.4 of [1].
%
%The function can either be called as
%a=getNextMPartition(n,m);
%to get the first parition of the integer n into m parts or as
%a=getNextMPartition(a);
%to get the next partition.
%
%Often, one might prefer all partitions of m integers that sum to a
%constant value n INCLUDING zeros. This function can be used to produce
%such a modified set of partitions by initializing it as
%a=getNextMPartition(n+m,m);
%and for every composition generated, including the first, the modified
%composition starting from zero is given by a-1.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Return the first partition in the sequence if two parameters are passed.
if(nargin>1)
    n=param1;
    a(1)=n-m+1;
    a(2:m)=1;
    return;
end

a=param1;

m=length(a);

%Deal with the special case of a 1D partition. In this instance, there is
%only 1 possible solution, so if one is calling this function again, it
%will be past the end of all possible solutions.
if(m==1)
    a=[];
    return;
end

while(a(2)<a(1)-1)
    %Step H3 in Knuth
    a(1)=a(1)-1;
    a(2)=a(2)+1;
    return;
end

%Step H4 in Knuth, with an added check to avoid reading past the end of the
%array if m=2.
j=3;
if(j>m)
    a=[];
    return;
end
s=a(1)+a(2)-1;
while(a(j)>=a(1)-1)
    s=s+a(j);
    j=j+1;
    %If we have passed the last partition
    %(Termination part of step H5)
    if(j>m)
       a=[];
       return;
    end
end

%Step H5 in Knuth.
x=a(j)+1;
a(j)=x;
j=j-1;

%Step H6 in Knuth
while(j>1)
    a(j)=x;
    s=s-x;
    j=j-1;
    a(1)=s;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% read a g2o data file describing a 2D SLAM instance
function graph = read_graph(filename)

fid = fopen(filename, 'r');

graph = struct (
  'x', [],
  'edges', [],
  'idLookup', struct
);

disp('Parsing File');
while true
  ln = fgetl(fid);
  if (ln == -1)
    break;
  end
  tokens = strsplit(ln, ' ', true);
  double_tokens = str2double(tokens);

  tk = 2;
  if (strcmp(tokens(1), 'VERTEX_SE2') != 0)
    id = int32(double_tokens(tk++));
    values = double_tokens(tk:tk+2)'; tk += 3;
    graph.idLookup = setfield(graph.idLookup, num2str(id), struct('offset', length(graph.x), 'dimension', length(values)));
    graph.x = [graph.x; values];
  elseif (strcmp(tokens(1), 'VERTEX_XY') != 0)
    id = int32(double_tokens(tk++));
    values = double_tokens(tk:tk+1)'; tk += 2;
    graph.idLookup = setfield(graph.idLookup, num2str(id), struct('offset', length(graph.x), 'dimension', length(values)));
    graph.x = [graph.x; values];
  elseif (strcmp(tokens(1), 'EDGE_SE2') != 0)
    fromId = int32(double_tokens(tk++));
    toId = int32(double_tokens(tk++));
    measurement = double_tokens(tk:tk+2)'; tk += 3;
    uppertri = double_tokens(tk:tk+5)'; tk += 6;
    information = [uppertri(1), uppertri(2), uppertri(3);
                   uppertri(2), uppertri(4), uppertri(5);
                   uppertri(3), uppertri(5), uppertri(6)];
    graph.edges = [graph.edges; struct(
      'type', 'P',
      'from', fromId,
      'to', toId,
      'measurement', measurement,
      'information', information)];
  elseif (strcmp(tokens(1), 'EDGE_SE2_XY') != 0)
    fromId = int32(double_tokens(tk++));
    toId = int32(double_tokens(tk++));
    measurement = double_tokens(tk:tk+1)'; tk += 2;
    uppertri = double_tokens(tk:tk+2)'; tk += 3;
    information = [uppertri(1), uppertri(2); uppertri(2), uppertri(3)];
    graph.edges = [graph.edges; struct(
      'type', 'L',
      'from', fromId,
      'to', toId,
      'measurement', measurement,
      'information', information)];
  end

end

% setup the index into the state vector
disp('Preparing helper structs');
for eid = 1:length(graph.edges)
  graph.edges(eid).fromIdx = getfield(graph.idLookup, num2str(graph.edges(eid).from)).offset + 1;
  graph.edges(eid).toIdx = getfield(graph.idLookup, num2str(graph.edges(eid).to)).offset + 1;
end

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function [Y,skel, channels] = loadMocapData()

%YA = vargplvmLoadData('hierarchical/demHighFiveHgplvm1',[],[],'YA');
%{
curDir = pwd;
cd ../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/02/
fileNameAsf='02.asf';
fileNameAmc='02_05.amc';
skel = acclaimReadSkel(fileNameAsf);
[channels, skel] = acclaimLoadChannels(fileNameAmc, skel);
% Remove root node?
channels(:, [1 3]) = zeros(size(channels, 1), 2);
skelPlayData(skel, channels, 1/40);
cd(curDir)
%}
%%
clear
close all
[Y, lbls, Ytest, lblstest,skel] = lvmLoadData2('cmuXNoRoot', '02', {'01','05','10'});
Yorig = Y;
% REmove motion in xz?

seq = cumsum(sum(lbls)) - [1:31];
% Ywalk = Y(1:3:70,:); % Orig: Y(1:85,:);
% Ypunch = Y(86:13:548,:); % Orig: Y(86:548,:);
% Ywash = Y(549:20:1209,:); % Orig:  Y(549:1209,:);
Ywalk = Y(1:2:70,:); % Orig: Y(1:85,:);
Ypunch = Y(86:13:548,:); % Orig: Y(86:548,:);
Ywash = Y(830:6:1140,:); % Orig:  Y(549:1209,:);

[channels xyzDiffIndices] = skelGetChannels(Ywalk);
Ywalk(:, xyzDiffIndices) = zeros(size(Ywalk(:, xyzDiffIndices) ));
[channels xyzDiffIndices] = skelGetChannels(Ypunch);
Ypunch(:, xyzDiffIndices) = zeros(size(Ypunch(:, xyzDiffIndices) ));
[channels xyzDiffIndices] = skelGetChannels(Ywash);
Ywash(:, xyzDiffIndices) = zeros(size(Ywash(:, xyzDiffIndices) ));

Y = [Ywalk; Ywash];
[channels] = skelGetChannels(Y);
close; skelPlayData(skel, channels, 1/5);

%%
%{
[channelsWalk] = skelGetChannels(Ywalk);
[channelsPunch] = skelGetChannels(Ypunch);
[channelsWash] = skelGetChannels(Ywash);
close; skelPlayData(skel, channelsWalk, 1/20);
close; skelPlayData(skel, channelsPunch, 1/20);
close; skelPlayData(skel, channelsWash, 1/20);
%}

%{
try
    load '../cmu13Data.mat'
catch
    [Y, lbls] = lvmLoadData2('cmu13');
    seq = cumsum(sum(lbls)) - [1:31];

    % load data
    [Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13NoRoot');
    skel = acclaimReadSkel('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13.asf');
    % (I think) any motion of the specific subject would do here, just to get
    % the channels
    [tmpchan, skel] = acclaimLoadChannels('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13_13.amc', skel);
end
%}

%{
See cmu49 for more coherent motions
try
    load '../cmu13Data.mat'
catch
    [Y, lbls] = lvmLoadData2('cmu13');
    seq = cumsum(sum(lbls)) - [1:31];

    % load data
    [Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13NoRoot');
    skel = acclaimReadSkel('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13.asf');
    % (I think) any motion of the specific subject would do here, just to get
    % the channels
    [tmpchan, skel] = acclaimLoadChannels('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13_13.amc', skel);
end
Yjump = Y(35:90,:);
Yjacks = Y(677:end,:);
channels = skelGetChannels(Yjump);
close; skelPlayData(skel, channels, 1/20);
%}

%%
function createCmuData(subject, motions)
baseDir = datasetsDirectory;
dirSep = filesep;
[Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13');
skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '13' dirSep '13.asf']);
seq = cumsum(sum(lbls)) - [1:31];
[tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '13' dirSep '13_16.amc'], skel);
[Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13NoRoot');

try
    load '../cmu13Data.mat'
catch
    [Y, lbls] = lvmLoadData2('cmu13');
    seq = cumsum(sum(lbls)) - [1:31];

    % load data
    [Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13NoRoot');
    skel = acclaimReadSkel('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13.asf');
    % (I think) any motion of the specific subject would do here, just to get
    % the channels
    [tmpchan, skel] = acclaimLoadChannels('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13_13.amc', skel);
end

%%
%{
%%%%%%%%%%%%%%%%%
%baseDir = datasetsDirectory;
%dirSep = filesep;
%[Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13');
%skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '13' dirSep '13.asf']);
%seq = cumsum(sum(lbls)) - [1:31];
%[tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '13' dirSep '13_16.amc'], skel);
%[Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13NoRoot');

try
    load '../cmu13Data.mat'
catch
    [Y, lbls] = lvmLoadData2('cmu13');
    seq = cumsum(sum(lbls)) - [1:31];

    % load data
    [Y, lbls, Ytest, lblstest] = lvmLoadData2('cmu13NoRoot');
    skel = acclaimReadSkel('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13.asf');
    % (I think) any motion of the specific subject would do here, just to get
    % the channels
    [tmpchan, skel] = acclaimLoadChannels('../../../vargplvmDEPENDENCIES/DATASETS0p1371/mocap/cmu/13/13_13.amc', skel);
end

Y = Y(1:12:1430,:);
Y = Y([1:45 72:end],:);
Y = Y([1:50 70:end],:);

[Ywalk, lbls, Ytest, lblstest] = lvmLoadData('cmu35gplvm');
Y = [Y; Ywalk(1:12:100,:)];

channels = skelGetChannels(Y);

%skelPlayData(skel, channels, 1/20);
%}

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function output = F_multi_softmax(input_layer, TaskVocabSizes)
input = input_layer.a;
[D,T,N] = size(input);

if N>1; [mask, variableLength] = CheckTrajectoryLength(input); end

if length(TaskVocabSizes)==1    % this is for when you have the same task, e.g. skip gram
    if D>TaskVocabSizes
        input2 = reshape(input, TaskVocabSizes, D/TaskVocabSizes, T, N);
        fakeLayer.a = input2;
        output = F_softmax(fakeLayer);
        output = reshape(output, D, T, N);
    else        % this is normal softmax
        output = F_softmax(input_layer);
    end

else    % this is for when you have different tasks
    % to be implemented
end

if N>1 && variableLength; output = PadShortTrajectory(output, mask, -1e10); end

end

author: singaxiong

repo: SignalGraph


text:
% function sem2legend

%legendArray = get(hndl,'DisplayName')
fnddata = findobj('parent',hndl);
legendArray = get(fnddata,'DisplayName');
legendArray = legendArray';
legendArray = legendArray(~cellfun(@isempty, legendArray));
legendArray2 = {''};
nnlg = 2*length(legendArray);
gg=length(legendArray); hh=1;
for bb=1:nnlg
        if mod(bb,2)~=0
                legendArray2{bb} = legendArray{gg};
                gg=gg-1;
        else
                legendArray2{bb} = sprintf('s.e.m. for Bin %g', hh);
                hh=hh+1;
        end
end
legendArray2 = fliplr(legendArray2);
legendArray2 = [{''} legendArray2];
set(fnddata,{'DisplayName'}, legendArray2')

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function ctr = ctrLoad(filename,xform)
% Load Contrack (ctr) structure to file.
%
%  ctr = ctrLoad(filename)
%
% Loads the Contrack structure to a file from ASCII format.
% I think there is a problem in how the xform is specified here.  By
% default, it seems to me that the transform is only for AC-PC space.
%
% Probably we should be more careful or helpful or something.

if (notDefined('filename'))
    filename = mrvSelectFile('r','txt',[],'Save Contrack parameters');
    if isempty(filename), disp('User canceled'); return; end
end
if notDefined('xform'), xform = eye(4); disp('Identity xform'); end

% Open  the file, checking that the user sent in a good file path
fileDir = fileparts(filename);
if ~isdir(fileDir), error('No directory %s',fileDir);
else fid = fopen(filename,'rt');
    if(fid == -1)
        error('Can not read ConTrack parameters file.');
    end
end

% What we will return.
ctr = ctrCreate;

ctr = ctrSet(ctr, 'tensors_filename', getStringFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'fa_filename', getStringFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'pdf_filename', getStringFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'mask_filename', getStringFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'Xmask_filename', getStringFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'require_way', getStringFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'desired_samples', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'max_nodes', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'min_nodes', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'step_size', getScalarFromLine(getLine(fid)));
coords = getBoxCoordsFromFile(fid,xform);
ctr = ctrSet(ctr,'roi', coords, 1, 'coords');
ctr = ctrSet(ctr,'roi', getStringFromLine(getLine(fid)), 1, 'seed_region');

%ctr = ctrSet(ctr,'roi', getStringFromLine(getLine(fid)), 1, 'valid_cortex');
coords = getBoxCoordsFromFile(fid,xform);
ctr = ctrSet(ctr,'roi', coords, 2, 'coords');
ctr = ctrSet(ctr,'roi', getStringFromLine(getLine(fid)), 2, 'seed_region');
%ctr = ctrSet(ctr,'roi', getStringFromLine(getLine(fid)), 2, 'valid_cortex');
ctr = ctrSet(ctr, 'save_out_spacing', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'fa_absorb', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'abs_normal', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'abs_penalty', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'smooth_std', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'angle_cutoff', getScalarFromLine(getLine(fid)));
ctr = ctrSet(ctr, 'shape_params_vec', getVectorFromLine(getLine(fid)));

fclose(fid);

return;

%--------------------
function line = getLine(fid)
line = fgetl(fid); if ~ischar(line),   error('Bad file.'),   end
return;

%--------------------
function string = getStringFromLine(line)
[foo string filepart] = strread(line,'%s %s %s','delimiter', ':');
%Reconstruct filename if necessary
if( ~isempty(filepart) ), string = strcat(string{:}, ':', filepart{:});
else                      string = string{:};
end
return;

%--------------------
function scalar  = getScalarFromLine(line)
[foo str_scalar] = strread(line,'%s %s','delimiter', ':');
scalar           = strread(str_scalar{:},'%f');
return;

%--------------------
function vector = get3VectorFromLine(line)
[foo str_vector] = strread(line,'%s %s','delimiter', ':');
vector = strread(str_vector{:},'%f','delimiter',',');
return;

%--------------------
function vector = getVectorFromLine(line)
[foo str_vector] = strread(line,'%s %s','delimiter', ':');
str_vector = strread(str_vector{:},'%*c %s','delimiter',']');
vector = strread(str_vector{:},'%f','delimiter',',');
return;

%--------------------
function coords = getBoxCoordsFromFile(fid,xform)
center = get3VectorFromLine(getLine(fid));
center = center(:)';
length = get3VectorFromLine(getLine(fid));
length = length(:)';
% Just put the corners of the box into coords
% XXX HACKY
% Assuming that we are going to add to the length to convert from points to
% voxels
coords = [center - (length-2)/2; center + (length-2)/2];
% Transform coords into the given space
coords = mrAnatXformCoords(xform, coords);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%*****************************************************************************
%% HSDsqlpcheckconvg: check convergence.
%%
%% ZpATynorm, AX, normX, normZ are with respect to the
%% original variables, not the HSD variables.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************************

function [param,breakyes,use_olditer,msg] = HSDsqlpcheckconvg(param,runhist)

termcode    = param.termcode;
iter        = param.iter;
obj         = param.obj;
relgap      = param.relgap;
gap         = param.gap;
prim_infeas = param.prim_infeas;
dual_infeas = param.dual_infeas;
mu          = param.mu;
prim_infeas_bad = param.prim_infeas_bad;
dual_infeas_bad = param.dual_infeas_bad;
printlevel  = param.printlevel;
stoplevel   = param.stoplevel;
inftol      = param.inftol;
gaptol      = param.gaptol;
kap         = param.kap;
tau         = param.tau;
theta       = param.theta;
breakyes    = 0;
use_olditer = 0;
msg         = [];
infeas      = max(prim_infeas,dual_infeas);
prim_infeas_min = min(param.prim_infeas_min, max(prim_infeas,1e-10));
dual_infeas_min = min(param.dual_infeas_min, max(dual_infeas,1e-10));
%%
err = max(infeas,relgap);
if (obj(2) > 0); homRd = param.ZpATynorm/obj(2); else homRd = inf; end
if (obj(1) < 0); homrp = norm(param.AX)/(-obj(1)); else homrp = inf; end
if (param.normX > 1e15*param.normX0 || param.normZ > 1e15*param.normZ0)
    termcode = 3;
    breakyes = 1;
end
if (homRd < min(1e-6,1e-2*sqrt(err*inftol)) && tau < 1e-4 ...
        && prim_infeas > 0.5*runhist.pinfeas(iter)) ...
        || (homRd < 10*tau && tau < 1e-7)
    termcode = 1;
    breakyes = 1;
end
if (homrp < min(1e-6,1e-2*sqrt(err*inftol)) && tau < 1e-4 ...
        && dual_infeas > 0.5*runhist.dinfeas(iter)) ...
        || (homrp < 10*tau && tau < 1e-7)
    termcode = 2;
    breakyes = 1;
end
if (err < gaptol)
    msg = sprintf('Stop: max(relative gap,infeasibilities) < %3.2e',gaptol);
    if (printlevel); fprintf('\n  %s',msg); end
    termcode = 0;
    breakyes = 1;
end
min_prim_infeas = min(runhist.pinfeas(1:iter));
prim_infeas_bad = prim_infeas_bad + (prim_infeas > ...
    max(1e-10,5*min_prim_infeas) && (min_prim_infeas < 1e-2));
if (mu < 1e-6)
    idx = max(1,iter-1): iter;
elseif (mu < 1e-3);
    idx = max(1,iter-2): iter;
else
    idx = max(1,iter-3): iter;
end
idx2 = max(1,iter-4): iter;
gap_ratio2 = runhist.gap(idx2+1)./runhist.gap(idx2);
gap_slowrate = min(0.8,max(0.6,2*mean(gap_ratio2)));
gap_ratio = runhist.gap(idx+1)./runhist.gap(idx);
pstep = runhist.step(iter+1);
if (infeas < 1e-4 || prim_infeas_bad) && (relgap < 1e-3) ...
        && (iter > 5) && (prim_infeas > (1-pstep/2)*runhist.pinfeas(iter))
    gap_slow = all(gap_ratio > gap_slowrate) && (relgap < 1e-3);
    min_pinfeas = min(runhist.pinfeas);
    if (relgap < 0.1*infeas) ...
            && ((runhist.step(iter+1) < 0.5) || (min_pinfeas < min(1e-6,0.1*prim_infeas))) ...
            && (dual_infeas > 0.9*runhist.dinfeas(iter) || (dual_infeas < 1e-2*gaptol))
        msg = 'Stop: relative gap < infeasibility';
        if (printlevel); fprintf('\n  %s',msg); end
        termcode = -1;
        breakyes = 1;
    elseif (gap_slow) && (infeas > 0.9*runhist.infeas(iter)) ...
            && (theta < 1e-8)
        msg = 'Stop: progress is too slow';
        if (printlevel); fprintf('\n  %s',msg); end
        termcode = -5;
        breakyes = 1;
    end
elseif (prim_infeas_bad) && (iter >50) && all(gap_ratio > gap_slowrate)
    msg = 'Stop: progress is bad';
    if (printlevel); fprintf('\n  %s',msg); end
    termcode = -5;
    breakyes = 1;
elseif (infeas < 1e-8) && (gap > 1.2*mean(runhist.gap(idx)))
    msg = 'Stop: progress is bad*';
    if (printlevel); fprintf('\n  %s',msg); end
    termcode = -5;
    breakyes = 1;
end
if (err < 1e-3) && (iter > 10) ...
        && (runhist.pinfeas(iter+1) > 0.9*runhist.pinfeas(max(1,iter-5))) ...
        && (runhist.dinfeas(iter+1) > 0.9*runhist.dinfeas(max(1,iter-5))) ...
        && (runhist.relgap(iter+1)  > 0.1*runhist.relgap(max(1,iter-5)));
    msg = 'Stop: progress is bad**';
    if (printlevel); fprintf('\n  %s',msg); end
    termcode = -5;
    breakyes = 1;
end
if (infeas > 100*max(1e-12,min(runhist.infeas)) && relgap < 1e-4)
    msg = 'Stop: infeas has deteriorated too much';
    if (printlevel); fprintf('\n  %s, %3.1e',msg,infeas); end
    use_olditer = 1;
    termcode = -7;
    breakyes = 1;
end
if (min(runhist.infeas) < 1e-4 || prim_infeas_bad) ...
        && (max(runhist.infeas) > 1e-4) && (iter > 5)
    relgap2 = abs(diff(obj))/(1+mean(abs(obj)));
    if (relgap2 < 1e-3);
        step_short = all(runhist.step(iter:iter+1) < 0.1) ;
    elseif (relgap2 < 1)
        idx = max(1,iter-3): iter+1;
        step_short = all(runhist.step(idx) < 0.05);
    else
        step_short = 0;
    end
    if (step_short)
        msg = 'Stop: steps too short consecutively';
        if (printlevel); fprintf('\n  %s',msg); end
        termcode = -5;
        breakyes = 1;
    end
end
if (iter > 3 && iter < 20) && (max(runhist.step(max(1,iter-3):iter+1)) < 1e-3) ...
        && (infeas > 1) && (min(homrp,homRd) > 1000*inftol)
    if (stoplevel >= 2)
        msg = 'Stop: steps too short consecutively';
        if (printlevel); fprintf('\n  %s',msg); end
        termcode = -5;
        breakyes = 1;
    end
end
if (pstep < 1e-4) && (err > 1.1*max(runhist.relgap(iter),runhist.infeas(iter)))
    msg = 'Stop: steps are too short';
    if (printlevel); fprintf('\n  %s',msg); end
    use_olditer = 1;
    termcode = -5;
    breakyes = 1;
end
if (iter == param.maxit)
    termcode = -6;
    msg = 'Stop: maximum number of iterations reached';
    if (printlevel); fprintf('\n  %s',msg); end
end
if (infeas < 1e-8 && relgap < 1e-10 && kap < 1e-13 && theta < 1e-15)
    msg = 'Stop: obtained accurate solution';
    if (printlevel); fprintf('\n  %s',msg); end
    termcode = 0;
    breakyes = 1;
end
param.prim_infeas_bad = prim_infeas_bad;
param.prim_infeas_min = prim_infeas_min;
param.dual_infeas_bad = dual_infeas_bad;
param.dual_infeas_min = dual_infeas_min;
param.termcode = termcode;
%%*****************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% 3D Hyperelastic registration of NIREP data to compare with LDDMM
% as described in Sec. 4 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
%
% =========================================================================

close all; clc; clear all;

setupNIREPDataNA02;
%setupNIREPDataNA03;
%setupNIREPDataNA04;
%setupNIREPDataNA05;

imgModel('reset','imgModel','linearInterMex');
trafo('reset','trafo','rigid3D');
distance('reset','distance','SSD');
viewImage('reset','viewImage','viewOrthoSlices2D');

% example specific parameters
switch example
    case ['3D-nirep-',dataset]
        alpha(1) = 100;
        alpha(2) = 1;
        alpha(3) = .1;
        alpha(4) = 1;
        pad = 0;
        minLevel = 5;
        maxLevel = 7;
        maxIterNPIR = 50;
    otherwise
        error('MP-LDDMM example %s - nyi');
end

Dshow = @(T,R,omega,m) viewIP(abs(T-R),omega,m,'colormap',gray(256));
FAIRplots('set','Dshow',Dshow);

solve = true;
postproc = true;

% 2) setup regularizer
regularizer('reset','regularizer','mfHyperElastic','alpha',alpha(1),...
    'alphaLength',alpha(2),'alphaArea',alpha(3),'alphaVolume',alpha(4));

filename = ['HyperElastic-',example,'-',regularizer,'-alpha-',num2str(alpha(1)),'-',num2str(alpha(2)),'-',num2str(alpha(3)),'-',num2str(alpha(4))];

if (solve)
    diary(fullfile('results',[filename,'.log']));
    % finally: run the MultiLevel Non-Parametric Image Registration
    [yc,wc,his] = MLIR(ML,'parametric',0,'NPIRLS',@ArmijoBacktrack,'minLevel',minLevel,'maxLevel',maxLevel,'maxIterNPIR',maxIterNPIR);
    diary('off')
    [~,interOpts] = inter;
    [~,distOpts ] = distance;
    [~,regOpts  ] = regularizer;

    save(fullfile('results',[filename,'.mat']),'yc','wc','his','m','minLevel','maxLevel','interOpts','distOpts','regOpts')
end

if postproc
    % generate figures
    close all;
    load(fullfile('results',[filename,'.mat']));

    computeOverlap;

    % generate figures
    figDir = ['results/',filename];
    if ~exist(figDir,'dir'), mkdir(figDir); end;
    dim  = numel(omega)/2;

    viewImage('reset','viewImage','viewOrthoSlices2D','colormap','gray');
    fig = figure(1); clf;
    fig.Name = 'dataR';
    viewImage(dataR,omega,m);
    cax = caxis;

    fig = figure(2); clf;
    fig.Name = 'dataT';
    viewImage(dataT,omega,m);
    caxis(cax);

    fig = figure(3); clf;
    fig.Name = 'res0';
    % viewImage(abs(dataR-dataT),omega,m)
    viewOrthoSlices2D(abs(dataR-dataT),omega,m,'color1',[0 0 0],'color2',[0 0 0]);
    colormap gray;
    colormap(flipud(colormap));
    % caxd = caxis;
    caxd = cax;

    Topt = imgModel(dataT,omega,center(yc,m));
    fig = figure(4); clf;
    fig.Name = 'Topt';
    viewImage(Topt,omega,m);
    % colormap gray
    caxis(cax);

    fig = figure(5); clf;
    fig.Name = 'resOpt';
    % viewImage(abs(dataR(:)-Topt),omega,m)
    viewOrthoSlices2D(abs(dataR(:)-Topt),omega,m,'color1',[0 0 0],'color2',[0 0 0]);
    caxis(caxd);
    colormap gray;
    colormap(flipud(colormap));

    fig = figure(6); clf;
    fig.Name = 'vol';
    Jac = geometry(yc,m,'Jac','matrixFree',1,'omega',omega);
    viewOrthoSlices2D(Jac,omega,m,'color1',[0 0 0],'color2',[0 0 0]);
    colormap jet;
    cb = colorbar('East');
    cpos = cb.Position;
    cb.Position(4) = cpos(4)*.4;
    cb.Position(2) = .56;
    cb.Ticks = [min(Jac) max(Jac)];
    caxis([min(Jac) max(Jac)]);
    cb.TickLabels = {sprintf('%1.2f',min(Jac)), sprintf('%1.2f',max(Jac))};
    cb.FontSize = 50;
    cb.Color = 'k';

    %str = {'dataR','dataT','res0','Topt','resOpt','yOpt','yInv','vol','T+yc'};
    str = {'dataR','dataT','res0','Topt','resOpt','vol'};
    for k=1:6
        figure(k);
        axis off;
        title([]);
        printFigure(gcf,fullfile(figDir,['LDDMM-' str{k} '_' regularizer '.png']),...
            'printOpts','-dbmp','printFormat','.bmp');
    end
end