text:
function [nav_e] = ins(imu, gnss, att_mode)
% ins: inertial navigation system.
%
% INPUT
%   imu, IMU data structure.
%         t: Ix1 time vector (seconds).
%        fb: Ix3 accelerations vector in body frame XYZ (m/s^2).
%        wb: Ix3 turn rates vector in body frame XYZ (radians/s).
%       arw: 1x3 angle random walks (rad/s/root-Hz).
%       vrw: 1x3 velocity random walks (m/s^2/root-Hz).
%     g_std: 1x3 gyros standard deviations (radians/s).
%     a_std: 1x3 accrs standard deviations (m/s^2).
%    gb_sta: 1x3 gyros static biases or turn-on biases (radians/s).
%    ab_sta: 1x3 accrs static biases or turn-on biases (m/s^2).
%    gb_dyn: 1x3 gyros dynamic biases or bias instabilities (radians/s).
%    ab_dyn: 1x3 accrs dynamic biases or bias instabilities (m/s^2).
%   gb_corr: 1x3 gyros correlation times (seconds).
%   ab_corr: 1x3 accrs correlation times (seconds).
%    gb_psd: 1x3 gyros dynamic biases root-PSD (rad/s/root-Hz).
%    ab_psd: 1x3 accrs dynamic biases root-PSD (m/s^2/root-Hz);
%      freq: 1x1 sampling frequency (Hz).
% ini_align: 1x3 initial attitude at t(1).
% ini_align_err: 1x3 initial attitude errors at t(1).
%
%	gnss, GNSS data structure.
%         t: Gx1 time vector (seconds).
%       lat: Gx1 latitude (radians).
%       lon: Gx1 longitude (radians).
%         h: Gx1 altitude (m).
%       vel: Gx3 NED velocities (m/s).
%       std: 1x3 position standard deviations (rad, rad, m).
%      stdm: 1x3 position standard deviations (m, m, m).
%      stdv: 1x3 velocity standard deviations (m/s).
%      larm: 3x1 lever arm from IMU to GNSS antenna (x-fwd, y-right, z-down) (m).
%      freq: 1x1 sampling frequency (Hz).
%   zupt_th: 1x1 ZUPT threshold (m/s).
%  zupt_win: 1x1 ZUPT time window (seconds).
%       eps: 1x1 time interval to compare current IMU time to current GNSS time vector (s).
%
%  att_mode: attitude mode string.
%      'quaternion': attitude updated in quaternion format. Default value.
%             'dcm': attitude updated in Direct Cosine Matrix format.
%
% OUTPUT
%   nav_e, INS/GNSS navigation estimates data structure.
%         t: Ix1 INS time vector (seconds).
%        tg: Gx1 GNSS time vector, when Kalman filter was executed (seconds).
%      roll: Ix1 roll (radians).
%     pitch: Ix1 pitch (radians).
%       yaw: Ix1 yaw (radians).
%       vel: Ix3 NED velocities (m/s).
%       lat: Ix1 latitude (radians).
%       lon: Ix1 longitude (radians).
%         h: Ix1 altitude (m).
%        xi: Gxn Kalman filter a priori states.
%        xp: Gxn Kalman filter a posteriori states.
%         z: Gxr INS/GNSS measurements
%         v: Gxr Kalman filter innovations.
%         b: Gxr Kalman filter biases compensations, [gb_dyn ab_dyn].
%         A: Gxn^2 Kalman filter transition-state matrices, one matrix per
%            row ordered by columns.
%        Pp: Gxn^2 Kalman filter a posteriori covariance matrices, one
%         matrix per row ordered by columns.
%        Pi: Gxn^2 Kalman filter a priori covariance matrices, one matrix
%            per row ordered by columns.
%         K: Gx(n*r) Kalman gain matrices
%         S: Gxr^2 Innovation matrices
%        ob: Gx1 Number of observable states after each GNSS data arriving
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics, vol. 17,
% issue 2, pp. 110-120, 2015. Alg. 2.
%
%   Groves, P.D. (2013), Principles of GNSS, Inertial, and
% Multisensor Integrated Navigation Systems (2nd Ed.). Artech House.
%
%   Crassidis, J.L. and Junkins, J.L. (2011). Optimal Esti-
% mation of Dynamic Systems, 2nd Ed. Chapman and Hall/CRC, USA.
%
%   ZUPT algothim based on Groves, Chapter 15, "INS Alignment, Zero Updates,
% and Motion Constraints".
%
%   ins_gps.m, ins_gnss function is based on that previous NaveGo function.
%
% Version: 012
% Date:    2022/07/19
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

if nargin < 3, att_mode  = 'quaternion'; end

%% ZUPT ALGORITHM

% zupt_flag = false;

%% PREALLOCATION

% Kalman filter dimensions
n = 15; % number of states
r = 6;  % number of sensors
% q = 12; % number of inputs

% Constant matrices
% I = eye(3);
% O = zeros(3);

% Length of INS time vector
LI = length(imu.t);

% Length of GNSS time vector
LG = length(gnss.t);

% Preallocation of attitude vectors
roll_e  = zeros (LI, 1);
pitch_e = zeros (LI, 1);
yaw_e   = zeros (LI, 1);

% Preallocation of velocity vector
vel_e   = zeros (LI, 3);

% Preallocation of gravity vector
gn_e   = zeros (LI, 3);

% Preallocation of position vectors
lat_e    = zeros (LI, 1);
lon_e    = zeros (LI, 1);
h_e      = zeros (LI, 1);

% Preallocation of Kalman filter matrices for later performance analysis
xi = zeros(LG, n);      % Evolution of Kalman filter a priori states
xp = zeros(LG, n);      % Evolution of Kalman filter a posteriori states
z = zeros(LG, r);       % INS/GNSS measurements
v = zeros(LG, r);       % Kalman filter innovations

A  = zeros(LG, n^2);    % Transition-state matrices
Pi = zeros(LG, n^2);    % A priori covariance matrices
Pp = zeros(LG, n^2);    % A posteriori covariance matrices
K  = zeros(LG, n*r);    % Kalman gain matrices
S  = zeros(LG, r^2);    % Innovation matrices
ob = zeros(LG, 1);      % Number of observable states at each GNSS data arriving

b = zeros(LG, r);       % Biases compensantions after Kalman filter correction

%% INITIAL VALUES AT INS TIME = 1

% Initial attitude
roll_e(1)  = imu.ini_align(1);
pitch_e(1) = imu.ini_align(2);
yaw_e(1)   = imu.ini_align(3);
DCMnb = euler2dcm([roll_e(1); pitch_e(1); yaw_e(1);]);
DCMbn = DCMnb';
qua   = euler2qua([roll_e(1) pitch_e(1) yaw_e(1)]);

% Initial velocity
vel_e(1,:) = gnss.vel(1,:);

% Initial position
lat_e(1) = gnss.lat(1);
lon_e(1) = gnss.lon(1);
h_e(1)   = gnss.h(1);

% Initial dynamic biases
gb_dyn = imu.gb_dyn';
ab_dyn = imu.ab_dyn';

% Turn-rates update with both updated velocity and position
omega_ie_n = earth_rate(lat_e(1));
omega_en_n = transport_rate(lat_e(1), vel_e(1,1), vel_e(1,2), h_e(1));

% Gravity update
gn_e(1,:) = gravity(lat_e(1), h_e(1));

%% INITIALIZATION OF KALMAN FILTER MATRICES

% Prior estimates
% kf.xi = [ zeros(1,9), imu.gb_dyn, imu.ab_dyn ]';  % Error vector state
% kf.Pi = diag([imu.ini_align_err, gnss.stdv, gnss.std, imu.gb_dyn, imu.ab_dyn].^2);
%
% kf.Q  = diag([imu.arw, imu.vrw, imu.gb_psd, imu.ab_psd].^2);
%
% fn = DCMbn * (imu.fb(1,:)' - ab_dyn - imu.ab_sta');
% wn = DCMbn * (imu.wb(1,:)' - gb_dyn - imu.gb_sta');
%
% % Vector to update matrix F
% upd = [gnss.vel(1,:) gnss.lat(1) gnss.h(1) fn' wn'];
%
% % Update matrices F and G
% [kf.F, kf.G] = F_update(upd, DCMbn, imu);
%
% [RM,RN] = radius(gnss.lat(1));
% Tpr = diag([(RM + gnss.h(1)), (RN + gnss.h(1)) * cos(gnss.lat(1)), -1]);  % radians-to-meters
%
% % Update matrix H
% kf.H = [ O I O O O ;
%     O O Tpr O O ; ];
% kf.R = diag([gnss.stdv gnss.stdm]).^2;
% kf.z = [ gnss.stdv, gnss.stdm ]';
%
% % Propagate prior estimates to get xp(1) and Pp(1)
% kf = kf_update( kf );
%
% % Initial matrices for Kalman filter performance analysis
% xi(1,:) = kf.xi';
% xp(1,:) = kf.xp';
% Pi(1,:) = reshape(kf.Pi, 1, n^2);
% Pp(1,:) = reshape(kf.Pp, 1, n^2);
% K(1,:)  = reshape(kf.K, 1, n*r);
% S(1,:)  = reshape(kf.S, 1, r^2);
% v(1,:)  = kf.v';
% z(1,:)  = kf.z';
% b(1,:) = [gb_dyn', ab_dyn'];

%% INS (IMU) TIME IS THE MASTER CLOCK
for i = 2:LI

    %% INERTIAL NAVIGATION SYSTEM (INS)

    % Print a dot on console every 10,000 INS executions
    if (mod(i,10000) == 0), fprintf('. ');  end
    % Print a return on console every 200,000 INS executions
    if (mod(i,200000) == 0), fprintf('\n'); end

    % IMU sampling interval
    dti = imu.t(i) - imu.t(i-1);

    % Inertial sensors corrected with a posteriori KF biases estimation and
    % deterministic static biases
    wb_corrected = imu.wb(i,:)' - gb_dyn - imu.gb_sta';
    fb_corrected = imu.fb(i,:)' - ab_dyn - imu.ab_sta';
    fn = DCMbn * fb_corrected;
%     wn = DCMbn * wb_corrected;

    % Velocity update
    vel = vel_update(fn, vel_e(i-1,:), omega_ie_n, omega_en_n, gn_e(i-1,:)', dti);
    vel_e (i,:) = vel;

    % Position update
    pos = pos_update([lat_e(i-1) lon_e(i-1) h_e(i-1)], vel_e(i,:), dti);
    lat_e(i) = pos(1);
    lon_e(i) = pos(2);
    h_e(i)   = pos(3);

    if ( h_e(i) < 0)
        h_e(i) = 10;
    end

    % Turn-rates update with both updated velocity and position
    omega_ie_n = earth_rate(lat_e(i));
    omega_en_n = transport_rate(lat_e(i), vel_e(i,1), vel_e(i,2), h_e(i));

    % Gravity update
    gn_e(i,:) = gravity(lat_e(i), h_e(i));

    % Attitude update
    [qua, DCMbn, euler] = att_update(wb_corrected, DCMbn, qua, ...
        omega_ie_n, omega_en_n, dti, att_mode);
    roll_e(i) = euler(1);
    pitch_e(i)= euler(2);
    yaw_e(i)  = euler(3);

    %% ZUPT DETECTION ALGORITHM
    idz = floor( gnss.zupt_win / dti ); % Index to set ZUPT window time

    if ( i > idz )

        % Mean velocity value for the ZUPT window time
        vel_m = mean (vel_e(i-idz:i , :));

        % If mean velocity value is under the ZUPT threshold velocity...
        if (abs(vel_m) < gnss.zupt_th)

            % Current attitude is equal to the mean of previous attitudes
            % inside the ZUPT window time
            roll_e(i)  = mean (roll_e(i-idz:i , :));
            pitch_e(i) = mean (pitch_e(i-idz:i , :));
            yaw_e(i)   = mean (yaw_e(i-idz:i , :));

            % Current position is equal to the mean of previous positions
            % inside the ZUPT window time
            lat_e(i) = mean (lat_e(i-idz:i , :));
            lon_e(i) = mean (lon_e(i-idz:i , :));
            h_e(i)   = mean (h_e(i-idz:i , :));

            % Alternative attitude ZUPT correction
            % roll_e(i)  = (roll_e(i-idz , :));
            % pitch_e(i) = (pitch_e(i-idz , :));
            % yaw_e(i)   = (yaw_e(i-idz, :));
            % lat_e(i) = (lat_e(i-idz:i , :));
            % lon_e(i) = (lon_e(i-idz:i , :));
            % h_e(i)   = (h_e(i-idz:i , :));

%             zupt_flag = true;

            % fprintf(' z\n')       % DEBUG
        end
    end

%     %% KALMAN FILTER UPDATE
%
%     % Check if there is a new GNSS measurement to process at current INS time
%     gdx =  find (gnss.t >= (imu.t(i) - gnss.eps) & gnss.t < (imu.t(i) + gnss.eps));
%
%     if ( ~isempty(gdx) && gdx > 1)
%
%         %  gdx       % DEBUG
%
%         %% MEASUREMENTS
%
%         % Meridian and normal radii of curvature update
%         [RM,RN] = radius(lat_e(i));
%
%         % Radians-to-meters matrix
%         Tpr = diag([(RM + h_e(i)), (RN + h_e(i)) * cos(lat_e(i)), -1]);
%
%         % Position innovations in meters with lever arm correction
%         zp = Tpr * ([lat_e(i); lon_e(i); h_e(i);] - [gnss.lat(gdx); gnss.lon(gdx); gnss.h(gdx);]) ...
%             + (DCMbn * gnss.larm);
%
%         % Velocity innovations with lever arm correction
%         zv = (vel_e(i,:) - gnss.vel(gdx,:) - ((omega_ie_n + omega_en_n) * (DCMbn * gnss.larm ))' ...
%             + (DCMbn * skewm(wb_corrected) * gnss.larm )' )';
%
%         %% KALMAN FILTER
%
%         % GNSS sampling interval
%         dtg = gnss.t(gdx) - gnss.t(gdx-1);
%
%         % Vector to update matrix F
%         upd = [vel_e(i,:) lat_e(i) h_e(i) fn' wn'];
%
%         % Matrices F and G update
%         [kf.F, kf.G] = F_update(upd, DCMbn, imu);
%
%         % Matrix H update
%         if(zupt_flag == false)
%             kf.H = [ O I O O O ;
%                 O O Tpr O O ; ];
%             kf.R = diag([gnss.stdv gnss.stdm]).^2;
%             kf.z = [ zv' zp' ]';
%         else
%             kf.H = [ O I O O O ; ];
%             kf.R = diag([gnss.stdv]).^2;
%             kf.z = zv;
%         end
%
%         % a posteriori states are forced to be zero (error-state approach)
%         kf.xp = zeros(n , 1);
%         % Execution of the extended Kalman filter
%         kf = kalman(kf, dtg);
%
%         %% OBSERVABILITY
%
%         % Number the observable states at current GNSS time
%         ob(gdx) = rank(obsv(kf.F, kf.H));
%
%         %% INS/GNSS CORRECTIONS
%
%         % Quaternion correction
%         qua_skew = -skewm(qua(1:3));    % According to Crassidis, qua_skew should be
%                                         % positive, but if positive NaveGo diverges.
%         % Crassidis, Eq. A.174a
%         Xi = [qua(4)*eye(3) + qua_skew; -qua(1:3)'];
%
%         % Crassidis, Eq. 7.34
%         qua = qua + 0.5 .* Xi * kf.xp(1:3);
%         qua = qua / norm(qua);          % Brute-force normalization
%
%         % DCM correction
%         DCMbn = qua2dcm(qua);
%
%         % Attitude correction, method 1
%         %         euler = qua2euler(qua);
%         %         roll_e(i) = euler(1);
%         %         pitch_e(i)= euler(2);
%         %         yaw_e(i)  = euler(3);
%
%         % Attitude correction, method 2
%         roll_e(i)  = roll_e(i)  - kf.xp(1);
%         pitch_e(i) = pitch_e(i) - kf.xp(2);
%         yaw_e(i)   = yaw_e(i)   - kf.xp(3);
%
%         % Velocity correction
%         vel_e(i,1) = vel_e(i,1) - kf.xp(4);
%         vel_e(i,2) = vel_e(i,2) - kf.xp(5);
%         vel_e(i,3) = vel_e(i,3) - kf.xp(6);
%
%         % Position correction
%         lat_e(i) = lat_e(i) - kf.xp(7);
%         lon_e(i) = lon_e(i) - kf.xp(8);
%         h_e(i)   = h_e(i)   - kf.xp(9);
%
%         % Biases estimation
%         gb_dyn   = -kf.xp(10:12);
%         ab_dyn   = -kf.xp(13:15);
%
%         % Matrices for later Kalman filter performance analysis
%         xi(gdx,:) = kf.xi';
%         xp(gdx,:) = kf.xp';
%         b(gdx,:) = [gb_dyn', ab_dyn'];
%         A(gdx,:)  = reshape(kf.A,  1, n^2);
%         Pi(gdx,:) = reshape(kf.Pi, 1, n^2);
%         Pp(gdx,:) = reshape(kf.Pp, 1, n^2);
%
%         if(zupt_flag == false)
%             v(gdx,:)  = kf.v';
%             z(gdx,:)  = kf.z';
%             K(gdx,:)  = reshape(kf.K, 1, n*r);
%             S(gdx,:)  = reshape(kf.S, 1, r^2);
%         else
%             zupt_flag = false;
%             z(gdx,:)  = [ kf.z' 0 0 0 ]';
%             v(gdx,:)  = [ kf.v' 0 0 0 ]';
%             K(gdx,1:n*3) = reshape(kf.K, 1, n*3);
%             S(gdx,1:9)  = reshape(kf.S, 1, 3^2);
%         end
%     end
end

%% Summary from INS/GNSS integration

nav_e.t     = imu.t(1:i, :);    % INS time vector
nav_e.tg    = gnss.t;           % GNSS time vector, which is the time vector when the Kalman filter was executed
nav_e.roll  = roll_e(1:i, :);   % Roll
nav_e.pitch = pitch_e(1:i, :);  % Pitch
nav_e.yaw   = yaw_e(1:i, :);    % Yaw
nav_e.vel   = vel_e(1:i, :);    % NED velocities
nav_e.lat   = lat_e(1:i, :);    % Latitude
nav_e.lon   = lon_e(1:i, :);    % Longitude
nav_e.h     = h_e(1:i, :);      % Altitude
nav_e.gn    = gn_e(1:i, :);     % Gravity estimation in the nav-frame.

nav_e.xi    = xi;       % A priori states
nav_e.xp    = xp;       % A posteriori states
nav_e.z     = z;        % INS/GNSS measurements
nav_e.v     = v;        % Kalman filter innovations
nav_e.b     = b;        % Biases compensations

nav_e.A     = A;        % Transition matrices
nav_e.Pi    = Pi;       % A priori covariance matrices
nav_e.Pp    = Pp;       % A posteriori covariance matrices
nav_e.K     = K;        % Kalman gain matrices
nav_e.S     = S;        % Innovation matrices
nav_e.ob    = ob;       % Number of observable states after each GNSS data arriving

fprintf('\n');

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function dydt=DouPenODE(t,y,m1,m2,a1,a2,L1,I1,I2,k1,k2,g)
dydt=[
    y(3,:);
    y(4,:);
    (-(I2*k1*y(3,:) + I2*k2*y(3,:) - I2*k2*y(4,:) + a2^2*k1*m2*y(3,:) + a2^2*k2*m2*y(3,:) - a2^2*k2*m2*y(4,:) + L1*a2^3*m2^2*y(4,:).^2.*sin(y(1,:) - y(2,:)) - (L1*a2^2*g*m2^2*sin(y(1,:)))/2 - I2*L1*g*m2*sin(y(1,:)) - (L1*a2^2*g*m2^2*sin(y(1,:) - 2*y(2,:)))/2 - I2*a1*g*m1*sin(y(1,:)) + (L1^2*a2^2*m2^2*y(3,:).^2.*sin(2*y(1,:) - 2*y(2,:)))/2 + L1*a2*k2*m2*y(3,:).*cos(y(1,:) - y(2,:)) - L1*a2*k2*m2*y(4,:).*cos(y(1,:) - y(2,:)) - a1*a2^2*g*m1*m2*sin(y(1,:)) + I2*L1*a2*m2*y(4,:).^2.*sin(y(1,:) - y(2,:)))./(I1*I2 + L1^2*a2^2*m2^2 + I2*L1^2*m2 + I2*a1^2*m1 + I1*a2^2*m2 - L1^2*a2^2*m2^2*cos(y(1,:) - y(2,:)).^2 + a1^2*a2^2*m1*m2));
    ((I1*k2*y(3,:) - I1*k2*y(4,:) + L1^2*k2*m2*y(3,:) - L1^2*k2*m2*y(4,:) + a1^2*k2*m1*y(3,:) - a1^2*k2*m1*y(4,:) + L1^3*a2*m2^2*y(3,:).^2.*sin(y(1,:) - y(2,:)) + L1^2*a2*g*m2^2*sin(y(2,:)) + I1*a2*g*m2*sin(y(2,:)) + (L1^2*a2^2*m2^2*y(4,:).^2.*sin(2*y(1,:) - 2*y(2,:)))/2 + L1*a2*k1*m2*y(3,:).*cos(y(1,:) - y(2,:)) + L1*a2*k2*m2*y(3,:).*cos(y(1,:) - y(2,:)) - L1*a2*k2*m2*y(4,:).*cos(y(1,:) - y(2,:)) - L1^2*a2*g*m2^2*cos(y(1,:) - y(2,:)).*sin(y(1,:)) + a1^2*a2*g*m1*m2*sin(y(2,:)) + I1*L1*a2*m2*y(3,:).^2.*sin(y(1,:) - y(2,:)) + L1*a1^2*a2*m1*m2*y(3,:).^2.*sin(y(1,:) - y(2,:)) - L1*a1*a2*g*m1*m2*cos(y(1,:) - y(2,:)).*sin(y(1,:)))./(I1*I2 + L1^2*a2^2*m2^2 + I2*L1^2*m2 + I2*a1^2*m1 + I1*a2^2*m2 - L1^2*a2^2*m2^2*cos(y(1,:) - y(2,:)).^2 + a1^2*a2^2*m1*m2))];

author: dynamicslab

repo: SINDy-PI


text:
%% alternating minimization scheme for low-rank+sparse decomposition

Z_norm=norm(Z(:));
stat.ener=[]; stat.ener(1)=vm.f(A,B);
if input_sig==1
A_norm=norm(A_0(:));
B_norm=norm(B_0(:));
stat.err_a=[]; stat.err_a(1)=norm(A(:)-A_0(:))/A_norm;
stat.err_b=[]; stat.err_b(1)=norm(B(:)-B_0(:))/B_norm;
end
stat.f_a=[]; stat.f_a(1)=norm(proj_fr(A+B-Z,U,V),'fro');
% stat.f_b=[]; %stat.f_b(1)=norm(proj_l0(A+B-Z,B),'fro');
% stat.ts1=[]; stat.ts2=[];
% stat.vi1=[]; stat.vi2=[];
stat.ang1=[]; stat.ang2=[];
stat.ls_count=[];
stat.cg_count=[];
stat.n3=[]; stat.n3(1)=n3;
stat.n4=[]; stat.n4(1)=n4;
stat.sfg_a=[]; stat.sfg_b=[];
stat.cpu=[]; %stat.cpu=time_main;
if ~isempty(A_ref), stat.diff_a=[]; stat.diff_a=norm(A(:)-A_ref(:))/A_norm; end
if ~isempty(B_ref), stat.diff_b=[]; stat.diff_b=norm(B(:)-B_ref(:))/B_norm; end
A_old=A; B_old=B;

% I1=eye(n1); I2=eye(n2);
waiting=waitbar(0,'iterating...');
%time_main=tic;
for k=1:slv.kmax
    waiting=waitbar(k/slv.kmax);

    % A-subproblem
    if sg.sig~=1
    % first shot via partial SVD
    [U,S,V]=lansvd(Z-B,n3);
    A=U*S*V';

    % safeguard by projected dogleg method
    v1=A_old(:)-A(:);
    stat.ang1(end+1)=(vm.f(A_old,B)-vm.f(A,B))/norm(v1(:))^2;

    if sg.sig==2 && stat.ang1(end)<sg.thr
        A=A_old;
        stat.sfg_a(end+1)=1;
        disp(['A safeguarded by Riemannian optim at iter ',num2str(k)])
        riem_optim
    else
        stat.sfg_a(end+1)=0;
    end

    else    % directly via projected dogleg method
        riem_optim
    end

    % B-subproblem by sorting
    B=rtr_l0(Z-A,n4);

    % safeguard
    v1=B_old(:)-B(:);
    stat.ang2(end+1)=(vm.f(A,B_old)-vm.f(A,B))/norm(v1(:))^2;

    if sg.sig~=0 && stat.ang2(end)<sg.thr
        B=Z-A; B(B_old==0)=0;
        stat.ang2(end)=.5;
        stat.sfg_b(end+1)=1;
        disp(['B safeguarded at iter ',num2str(k)])
    else
        stat.sfg_b(end+1)=0;
    end

%     fprintf('|A^k-A^{k-1}|=%1.4e\n',norm(A(:)-A_old(:)))

    % trimming
    trimming

    % stats
    %stat.cpu(end+1)=stat.cpu(1)+toc(time_main);
    stat.ener(end+1)=vm.f(A,B);
    if input_sig==1
    stat.err_a(end+1)=norm(A(:)-A_0(:))/A_norm;
    stat.err_b(end+1)=norm(B(:)-B_0(:))/B_norm;
    end
    stat.f_a(end+1)=norm(proj_fr(vm.g(A,B_old),U,V),'fro');
%     stat.f_b(end+1)=norm(proj_l0(A+B-Z,B),'fro');
    if ~isempty(A_ref), stat.diff_a(end+1)=norm(A(:)-A_ref(:))/A_norm; end
    if ~isempty(B_ref), stat.diff_b(end+1)=norm(B(:)-B_ref(:))/B_norm; end
    stat.n3(end+1)=n3; stat.n4(end+1)=n4;

%     if sg.sig==1, cg.tol=1e-2*(stat.f_a(end)/stat.f_a(1))^.2; end

    % stopping rule
    if stat.f_a(end)<slv.ktol*stat.f_a(1)
%     if stat.err_a(end)<slv.ktol
% 	if stat.ener(end-1)-stat.ener(end)<slv.ktol*stat.ener(end-1)
        break
    end

    A_old=A;
    B_old=B;
end
%toc(time_main)
close(waiting)

% [sum(stat.ls_count),sum(stat.cg_count)]

author: andrewssobral

repo: lrslibrary


text:
% OtTableFields Get a list of fields in a string array.
%
% Usage:
%    list = OtTableFields(htable)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * 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
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% 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 OWNER 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.

% 2006-02-28 Created (Johan Meijdam, VORtech)

function list = OtTableFields(htable)
   list = Ot('table.fields', htable);
end

author: Sable

repo: mcbench-benchmarks


text:
 function pot = potential_fun(ptype, delta, param, varargin)
%function pot = potential_fun(ptype, delta, param, varargin)
%|
%| Define roughness penalty potential functions (as strum).
%|
%| The penalty will have the form
%|	R(x) = sum_k w_k * potential_k([Cx]_k, delta_k)
%| where w_k is provided elsewhere, not here!
%|
%| in
%|	ptype		quad broken huber hyper2 hyper3 cauchy qgg2 genhub
%|			lange1 lange3 (fair) geman&mcclure gf1 li98cfs
%|			Recommended: 'hyper3'
%|			To list all possible choices, use:
%|			potential_fun('list') for "smooth" options
%|			potential_fun('list1') for "non-smooth" options
%|	delta		scalar, or image-sized array;
%|			"cutoff" parameter for edge-preserving regularization
%|	param		optional additional parameter(s) for some choices:
%|				'gf1' (generalized Fair) : [a b]
%|				'qgg2' : q
%|				'genhub' & 'stevenson94dpr' : [p q]
%|				'table1' : [dz, dpot([1:K] * dz)]
%|
%| option
%|	'dummy'		include extra dummy argument for backward compatibility
%|			e.g., pot.potk([], C*x). default: 0
%|
%| out
%|	pot		strum object, with data: delta and param
%|	methods:
%|		pot.potk(C*x)	potential function value
%|		pot.wpot(C*x)	potential 'weights' (aka half-quad. curvatures)
%|		pot.dpot(C*x)	potential derivative
%|		pot.shrink(b, reg)	proximal (shrinkage) operation:
%|					argmin_z 1/2 |z - b|^2 + reg * pot(z)
%|		pot.plot()	plot all of the above functions
%|
%| trick: for ptype 'gf1-fit', the param argument should be:
%|	{'name of potential to fit', points, param}
%| and this function returns the [a b] parameters needed for a subsequent
%| call with ptype 'gf1'
%|
%| Copyright 2004-5-18, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(ptype, 'test') ir_potential_fun_test, return, end
if streq(ptype, 'test_lpp') ir_potential_fun_test_lpp, return, end

if streq(ptype, 'list') % return list of all "smooth" choices
	pot = {'quad', 'huber', 'huber2', ...
	'hyper2', 'hyper3', 'cauchy', ...
	'qgg2', 'genhub', 'lange1', 'lange3', 'geman&mcclure', 'gf1', ...
	'table0', 'table1', 'li98cfs'};
%	'stevenson94dpr' omit this one; use 'genhub' instead
return
end

if streq(ptype, 'list1') % return list of 'l1' type choices (mostly non-smooth)
	pot = {'l0', 'l1', 'lpp', 'tav', 'broken', 'fair-l1'};
return
end

if ~isvar('delta'), delta = []; end
if ~isvar('param'), param = []; end

arg.dummy = false;
arg.scale = 1;
arg = vararg_pair(arg, varargin);

if streq(ptype, 'gf1-fit') % trick: just return parameters for this case
	pot = ir_potential_fun_gf1_fit(param{:});
return
end

[pot.type, pot.delta, pot.param, potk, wpot, dpot, shrink] = ...
	ir_potential_fun_parse(ptype, delta(:), param(:));

if arg.scale ~= 1
	potk = @(pot, z) arg.scale * potk(pot, z);
	dpot = @(pot, z) arg.scale * dpot(pot, z);
	wpot = @(pot, z) arg.scale * wpot(pot, z);
end

if arg.dummy
	potk = @(pot, dum, z) potk(pot, z);
	dpot = @(pot, dum, z) dpot(pot, z);
	wpot = @(pot, dum, z) wpot(pot, z);
end

meth = {'potk', potk, 'wpot', wpot, 'dpot', dpot, 'shrink', shrink, ...
	'plot', @ir_potential_fun_plot};
pot = strum(pot, meth);

end % potential_fun()

% ir_potential_fun_gf1_fit()
% trick: gf1-fit means match gf1 to named potential at given points [0, s1, s2]
% where the given point values are s = z / delta, i.e., in "delta" units
% param: {'name of potential to fit', points, param}
function ab = ir_potential_fun_gf1_fit(ptype, sv, param)
sv = sv(:);
pot = potential_fun(ptype, 1, param); % trick: delta=1 here
pt = pot.wpot(sv);

% ab = [1 1; -pt']' \ ((pt-1) ./ sv);

s1 = sv(1);
s2 = sv(2);
w1 = pt(1);
w2 = pt(2);
ab(1) = (w2 * s2 * (1 - w1) - w1 * s1 * (1 - w2)) ...
		/ ((w1 - w2) * s1 * s2);
ab(2) = (s2 * (1 - w1) - s1 * (1 - w2)) ...
		/ ((w1 - w2) * s1 * s2);
end % ir_potential_fun_gf1_fit()

% ir_potential_fun_parse()
function [ptype, delta, param, potk, wpot, dpot, shrink] = ...
	ir_potential_fun_parse(ptype, delta, param)

dpot = [];
shrink = [];

% trick: huber2 is just a huber with delta / 2
% so that weighting function drops to 1/2 at delta, like hyper3 etc.
if streq(ptype, 'huber2')
	ptype = 'huber';
	delta = delta / 2;
end

% trick: hyper3 is just a hyperbola with delta scaled by sqrt(3)
% to approximately "match" the properties of 'cauchy' (old erroneous 'hyper')
if streq(ptype, 'hyper3')
	ptype = 'hyper2';
	delta = delta / sqrt(3);
end

if streq(ptype, 'gf1') && param(1) == 0
	if param(2) ~= 1
		fail 'only b=1 makes sense for gf1 with a=0'
	end
	ptype = 'lange3'; param = []; % trick: gf1 with a=0 and b=1 is lange3
end

switch ptype

% quadratic potential function
case 'quad'
	potk = @(pot, z) (abs(z).^2) / 2;
	wpot = @(pot, z) ones(size(z), class(z));
	dpot = @(pot, z) z;
	shrink = @(pot, b, reg) b ./ (1 + reg);

% broken parabola
case 'broken'
	potk = @(pot, z) min(z.^2, pot.delta.^2)/2;
	wpot = @(pot, z) ir_tonumeric(abs(z) < pot.delta, z);
	dpot = @(pot, z) z .* (abs(z) < pot.delta);
	shrink = @(pot, b, reg) ir_broken_shrink(b, reg, pot.delta);

% l0 (hard threshold)
case 'l0'
	potk = @(pot, z) ir_tonumeric(z ~= 0, z);
	wpot = @(pot, z) nan(size(z), class(z));
	dpot = @(pot, z) nan(size(z), class(z));
	shrink = @(pot, b, reg) b .* (abs(b) > sqrt(2*reg));

% absolute value (i.e., l1, for soft threshold)
case {'abs', 'l1'}
	potk = @(pot, z) abs(z);
	wpot = @(pot, z) 1 ./ abs(z); % invalid at z=0
	dpot = @(pot, z) sign(z);
	shrink = @(pot, b, reg) sign(b) .* max(abs(b) - reg, 0);

% l_p^p aka generalized gaussian (gg)
case 'lpp'
	potk = @(pot, z) abs(z).^param;
	dpot = @(pot, z) param * sign(z) .* abs(z).^(param-1); % invalid at z=0
	wpot = @(pot, z) param .* abs(z).^(param-2); % invalid at z=0
	switch param
	case 0
		potk = @(pot, z) ir_tonumeric(z ~= 0, z);
		dpot = @(pot, z) zeros(size(z), class(z)); % meaningless
		shrink = @(pot, b, reg) b .* (abs(b) > sqrt(2*reg));
	case 1
		potk = @(pot, z) abs(z);
		dpot = @(pot, z) sign(z); % invalid at t=0
		shrink = @(pot, b, reg) sign(b) .* max(abs(b) - reg, 0);
	case 0.5 % p=1/2
		shrink = @(pot, b, reg) ir_lpp_1_2_shrink(b, reg);
	case 4/3. % p=4/3
		shrink = @(pot, b, reg) ir_lpp_4_3_shrink(b, reg);
	case 1.5 % p=3/2
		shrink = @(pot, b, reg) sign(b) .* ...
			( sqrt((3/4*reg).^2 + abs(b)) - 3/4*reg ).^2;
		% from eqn. (4.5) in chaux:07:avf
	%	z + 9 * reg.^2 * sign(z) .* (1 - sqrt(1 + 16*abs(z)/(9*reg^2))) / 8;

	otherwise
		warn('no shrink for p=%g', param)
	end

% truncated absolute value
case 'tav'
	potk = @(pot, z) min(abs(z), pot.delta);
	wpot = @(pot, z) (1 ./ abs(z)) .* (abs(z) < pot.delta); % bad at t=0
	dpot = @(pot, z) sign(z) .* (abs(z) < pot.delta);
	shrink = @(pot, b, reg) ir_tav_shrink(b, reg, pot.delta);

% huber potential function
case 'huber'
	potk = @(pot, z) huber_pot(z, pot.delta);
	wpot = @(pot, z) huber_wpot(z, pot.delta);
	dpot = @(pot, z) huber_dpot(z, pot.delta);
	shrink = @(pot, b, reg) ir_huber_shrink(b, reg, pot.delta);

% cauchy penalty: d^2 / 2 * log(1 + (t/d)^2) (not convex!)
case 'cauchy'
	potk = @(pot, z) pot.delta.^2 / 2 .* log(1 + abs(z ./ pot.delta).^2);
	wpot = @(pot, z) 1 ./ (1 + abs(z ./ pot.delta).^2);
	dpot = @(pot, z) z ./ (1 + abs(z ./ pot.delta).^2);
	shrink = @(pot, b, reg) cauchy_shrink(b, reg, pot.delta);

% Geman&McClure penalty: d^2 / 2 * |z/d|^2 / (1 + |z/d|^2)
% Not convex!
case 'geman&mcclure'
	potk = @(pot, z) pot.delta.^2 / 2 .* abs(z/pot.delta).^2 ./ (1 + abs(z ./ pot.delta).^2);
	wpot = @(pot, z) 1 ./ (1 + abs(z ./ pot.delta).^2).^2;
	dpot = @(pot, z) z ./ (1 + abs(z ./ pot.delta).^2).^2;

% gf1: Generalized Fair 1st-order
% wpot(z) = (1 + a * |z/d|) / (1 + b * |z/d|)
case 'gf1'
	potk = @(pot, z) gf1_potk(z, pot.delta, pot.param(1), pot.param(2));
	wpot = @(pot, z) (1 + pot.param(1) .* abs(z ./ pot.delta)) ...
		./ (1 + pot.param(2) .* abs(z ./ pot.delta));
	shrink = @(pot, b, reg) ...
		ir_gf1_shrink(b, reg, pot.delta, pot.param(1), pot.param(2));

% hyperbola penalty: d^2 * [ sqrt(1 + (z/d)^2) - 1 ]
case 'hyper2'
	potk = @(pot, z) pot.delta.^2 .* (sqrt(1 + abs(z ./ pot.delta).^2) - 1);
	wpot = @(pot, z) 1 ./ sqrt(1 + abs(z ./ pot.delta).^2);
	dpot = @(pot, z) z ./ sqrt(1 + abs(z ./ pot.delta).^2);

case 'hyper'
	error 'use "cauchy" or "hyper3" not "hyper" now'

% Lange1 penalty
case 'lange1'
	potk = @(pot, z) abs(z).^2 / 2 ./ (1+abs(z./pot.delta));
	wpot = @(pot, z) (1 + abs(z ./ pot.delta) / 2) ./ (1 + abs(z ./ pot.delta)).^2;

% Lange3 penalty
case {'lange3', 'fair'}
	potk = @(pot, z) pot.delta.^2 .* (abs(z./pot.delta) - log(1+abs(z./pot.delta)));
	wpot = @(pot, z) 1 ./ (1 + abs(z ./ pot.delta));
	dpot = @(pot, z) z ./ (1 + abs(z ./ pot.delta));
	% caution: no built-in shink, use 'fair-l1' if you want shrink!

% Fair potential "rounded corner" approximation to l1
case 'fair-l1'
	potk = @(pot, z) abs(z) - pot.delta .* log(1+abs(z./pot.delta));
	wpot = @(pot, z) 1 ./ (pot.delta + abs(z));
	dpot = @(pot, z) z ./ (pot.delta + abs(z));
	shrink = @(pot, b, reg) fair_l1_shrink(b, reg, pot.delta);

% li98cfs
case 'li98cfs'
	% f = @(x) atan(x) / x - 0.5; fsolve(f, 2.3)
	delta = delta / 2.3311;
	potk = @(pot, z) ir_li98cfs_potk(z, pot.delta);
	wpot = @(pot, z) ir_li98cfs_wpot(z, pot.delta);

% qgg2: q-generalized gaussian for p=2, due to Thibault, Sauer, Bouman
% q = "param", same as lange1 when q=1
case 'qgg2'
	potk = @(pot, z) z.^2 / 2 ./ (1+abs(z./pot.delta).^(2-pot.param));
	wpot = @(pot, z) (1 + abs(z ./ pot.delta).^(2-pot.param) * pot.param ...
		 / 2) ./ (1 + abs(z ./ pot.delta).^(2-pot.param)).^2;

% genhub : generalized Huber (switch between two generalized gaussians)
% same as Huber when p=2 and q=1
% p is power near 0, q is asymptotic power
case 'genhub'
	potk = @(pot, z) ir_genhub_potk(z, pot.delta, pot.param(1), pot.param(2));
	wpot = @(pot, z) ir_genhub_wpot(z, pot.delta, pot.param(1), pot.param(2));

% from stevenson:94:dpr
% p = param(1), q = param(2), same as Huber when p=2 and q=1 ???
case 'stevenson94dpr'
	warn('Potential %s does not work well; use "genhub" instead', ptype)
	potk = @(pot, z) ir_stevenson94dpr_potk(z, pot.delta, pot.param(1), pot.param(2));
	wpot = @(pot, z) ones(size(z), class(z)); % bogus attempt at upper bound

% tabulate derivative of pot at z_k = dz * k for k=1,...,K
% param is samples of derivative dpot(z_k)
% dpot(0) = 0, and use sample-and-hold interpolation of dpot()
% except use linear interpolation over [0 dz]
case 'table0'
	if ~isnumeric(param) || numel(param) < 10
%	if ~iscell(param) || numel(param) ~= 2
%		fail 'for table0 param must be {dz, dpot([1:K]*dz}'
		fail 'for table0 param must be [dz, dpot([1:K]*dz)]'
	end
%	param = ir_table0_setup(param{1}, col(param{2}));
	param = ir_table0_setup(param(1), col(param(2:end)));
	potk = @(pot, z) ir_table0_potk(z, pot.param);
	wpot = @(pot, z) ir_table0_wpot(z, pot.param);
	shrink = @(pot, b, reg) ir_table0_shrink(b, reg, pot.param);

% tabulate derivative of pot at z_k = dz * k for k=1,...,K
% param is samples of derivative dpot(z_k)
% dpot(0 = 0, and use linear interpolation of dpot(t)
case 'table1'
	if ~isnumeric(param) || numel(param) < 10
%	if ~iscell(param) || numel(param) ~= 2
%		fail 'for table1 param must be {dz, dpot([1:K]*dz}'
		fail 'for table1 param must be [dz, dpot([1:K]*dz)]'
	end
%	param = ir_table1_setup(param{1}, col(param{2}));
	param = ir_table1_setup(param(1), col(param(2:end)));
	potk = @(pot, z) ir_table1_potk(z, pot.param);
	wpot = @(pot, z) ir_table1_wpot(z, pot.param);
	shrink = @(pot, b, reg) ir_table1_shrink(b, reg, pot.param);

otherwise
	fail('Unknown potential "%s"', ptype)
end

if isempty(dpot) % default is z * wpot(z)
	dpot = @(pot, z) z .* wpot(pot, z);
end
if isempty(shrink)
	shrink = @ir_potential_fun_shrink;
end

end % ir_potential_fun_parse()

% ir_potential_fun_shrink()
% find argmin_z 1/2 |z - b|^2 + reg * pot(z)
% default method uses fzero() which will be slow!
function out = ir_potential_fun_shrink(pot, b, reg)
out = zeros(size(b), class(b));
opt = optimset('tolx', 1e-7);
for ii=1:numel(b)
	a = abs(b(ii));
	s = sign(b(ii));
	cost = @(z) 0.5 * (z - a).^2 + reg * pot.potk(z);
%	dfun = @(z) z - a + reg * pot.dpot(z);
	try
	%	out(ii) = s * fzero(dfun, a);
	%	out(ii) = s * fminsearch(cost, a, opt);
		z0 = s * fminsearch(cost, 0, opt);
		z1 = s * fminsearch(cost, a, opt);
		if cost(z0) < cost(z1)
			out(ii) = z0;
		else
			out(ii) = z1;
		end
	catch
		whos
		keyboard
	end
end % for
end % ir_potential_fun_shrink

% ir_table0_setup()
% dz	[1]	z spacing
% dk	[K]	dpot([1:K] * dz)
function out = ir_table0_setup(dz, dk)
sk = dk(1)*dz/2 + dz * [0; cumsum(dk(1:end-1))]; % [K]
out.dz = dz;
out.dk = dk; % [1:K] samples of dpot
out.sk = sk; % [K] cumulative sums for pot(z)
K = numel(dk);
out.K = K;
end % ir_table0_setup

% ir_table0_potk()
function out = ir_table0_potk(z, param)
dz = param.dz;
dk = param.dk;
sk = param.sk;
z = abs(z);
k = floor(z / dz);
k = min(k, param.K); % at most K
big = z > dz;
out = zeros(size(z), class(z));
out(~big) = 0.5 * dk(1) / dz * z(~big).^2;
k = k(big);
out(big) = sk(k) + dk(k) .* (z(big) - k * dz);
end % ir_table0_potk

% ir_table0_wpot()
function out = ir_table0_wpot(z, param)
dz = param.dz;
dk = param.dk; % [K]
z = abs(z);
k = floor(z / dz);
k = min(k, param.K); % at most K
out = repmat(dk(1) / dz, size(z));
big = k > 0;
k = k(big);
z = z(big);
out(big) = dk(k) ./ z;
end % ir_table0_wpot

% ir_table0_shrink()
function out = ir_table0_shrink(b, reg, param)
K = param.K;
dk = param.dk;
dz = param.dz;
zk = (1:K)' * dz;
bk = zk + reg * dk; % must be done here because depends on reg
ck = bk + dz;
tmp = [0; col([bk ck]')];
out = [zk zk+dz];
out = [0; col([zk zk+dz]')];
out = sign(b) .* interp1(tmp, out, abs(b), 'linear', 'extrap');
end % ir_table0_shrink

% ir_table1_setup()
% dz	[1]	z spacing
% dk	[K]	dpot([1:K] * dz)
function out = ir_table1_setup(dz, dk)
K = numel(dk);
k = [0:K]'; % [K+1]
zk = dz * k; % [K+1]
dk0 = [0; dk]; % [K+1] prepend sample at zero
ck = [diff(dk0); 0] / dz; % [K+1] curvatures 0:K
tmp = (dk0 - zk .* ck) * dz ...
	+ dz^2/2 * ck .* ((k+1).^2 - k.^2); % [K+1]
sk = cumsum( [0; tmp] );
out.dz = dz;
out.dk = dk; % [1:K] samples of dpot
out.ck = ck; % [K+1] curvatures 0:K
out.sk = sk; % [K] cumulative sums for pot(t)
out.K = K;
end % ir_table1_setup

% ir_table1_potk()
function out = ir_table1_potk(z, param)
sk = param.sk;
dz = param.dz;
ck = param.ck;
dk0 = [0; param.dk]; % [K+1] prepend sample at zero
z = abs(z);
k = floor(z / dz);
k = min(k, param.K); % at most K
sk = sk(1+k); % matlab indexing
ck = ck(1+k); % matlab indexing
dk = dk0(1+k);
zk = k * dz;
ck = reshape(ck, size(z));
dk = reshape(dk, size(z));
sk = reshape(sk, size(z));
out = sk + (dk - zk .* ck) .* (z - zk) + ck / 2 .* (z.^2 - zk.^2);
end % ir_table1_potk

% ir_table1_wpot()
function out = ir_table1_wpot(z, param)
dz = param.dz;
ck = param.ck;
dk = param.dk;
z = abs(z);
k = floor(z / dz);
k = min(k, param.K); % at most K
out = repmat(dk(1) / dz, size(z));
big = k > 0;
k = k(big);
z = z(big);
dk = reshape(dk(k), size(z));
ck = reshape(ck(k+1), size(z));
out(big) = (dk + (z - k * dz) .* ck) ./ z;
end % ir_table1_wpot

% ir_table1_shrink()
% unfortunately this routine works only for a single scalar reg value.
function out = ir_table1_shrink(z, reg, param)
K = param.K;
zk = (1:K)' * param.dz;
bk = zk + reg * param.dk; % must be done here because depends on reg
out = sign(z) .* interp1([0; bk], [0; zk], abs(z), 'linear', 'extrap');
end % ir_table1_shrink

% gf1_potk()
% gf1: generalized fair 1st-order potential
function pot = gf1_potk(z, delta, a, b)
atd = abs(z ./ delta);
pot = delta.^2 ./ (2 * b.^3) * ...
	(2 * b.^2 .* atd + a .* b.^2 .* atd.^2 ...
	- 2 * a .* b .* atd + 2 * (a-b) .* log(1 + b .* atd));

if 0 % symbolic check
	syms x
	syms a positive
	syms b positive
	syms z positive
	int(x*(1+a*x) / (1+b*x), x, 0, z)
end
end % gf1_potk

% ir_broken_shrink()
function out = ir_broken_shrink(z, reg, delta)
out = z ./ (1 + reg);
big = delta * (1 + reg) < abs(z);
out(big) = z(big);
end % ir_broken_shrink

% cauchy_shrink()
function out = cauchy_shrink(z, reg, delta)
z = z(:);
coef = ones(numel(z),1);
coef = [1/delta^2*coef -z/delta^2 (1+reg)*coef -z];
out = zeros(size(z));
for ii=1:numel(z)
	tmp = roots(coef(ii,:));
	pick = tmp == real(tmp); % empirically, 3rd root is often real
	if sum(pick) ~= 1 % if multiple real roots, empirically pick largest
		pick = imax(abs(tmp));
	end
	out(ii) = tmp(pick);
end
end % cauchy_shrink

% ir_huber_shrink()
function out = ir_huber_shrink(z, reg, delta)
out = z ./ (1 + reg);
big = delta .* (1 + reg) < abs(z);
if numel(reg) > 1, reg = reg(big); end
if numel(delta) > 1, delta = delta(big); end
z = z(big);
out(big) = z .* (1 - reg .* delta ./ abs(z));
end % ir_huber_shrink

% fair_l1_shrink()
function out = fair_l1_shrink(z, reg, delta)
out = sign(z) .* (abs(z) - (delta + reg) ...
	+ sqrt( (delta + reg - abs(z)).^2 + 4 * delta .* abs(z) )) ./ 2;
end % fair_l1_shrink

% ir_gf1_shrink()
function out = ir_gf1_shrink(z, reg, delta, a, b)
u = a / delta;
v = b / delta;
out = sign(z) .* (v .* abs(z) - (1+reg) ...
	+ sqrt( (1 + reg - v.*abs(z)).^2 + 4 * (v + reg .* u) .* abs(z) )) ...
	./ (2 * (v + reg .* u));
end % ir_gf1_shrink

% ir_lpp_1_2_shrink()
% for l_p with p=1/2
% for a > 0 and z = t^2 the minimizer solves t^3 - a t + reg/2 = 0
function out = ir_lpp_1_2_shrink(b, reg)
sb = sign(b); a = abs(b);
%reg = 1;
%a = 50; roots([1 0 -a reg/2])
%a = linspace(0, 4, 401);
out = zeros(size(b), class(b));
%big = true(size(b));
%{
% loop method
for ib=1:numel(b)
	tmp = roots([1 0 -a(iz) reg/2])'
	out(ib) = max(tmp).^2;
end
%}
%{
% hyperbolic method for one real root:
big = 27 * (reg/2)^2 > 4 * a.^3;
t0 = -2 * sqrt(a/3) .* cosh(1/3 * acosh(3/2 * (reg/2) ./ a .* sqrt(3./a)));
%}
%{
% https://en.wikipedia.org/wiki/Cubic_function#Vieta.27s_substitution
% vieta solves t^3 + p t + q = 0 using t = w - p / (3w), i.e. p=-a and q=reg/2
% for which w^6 + q w^3 - p^3/27 = 0
% i.e. (w^3)^2 + reg/2 (w^3) + a^3/27 = 0
w3 = (-reg/2 + sqrt((reg/2).^2 - 4 * a.^3/27)) / 2; % w^3 solution by quad form
w1 = w3 .^ (1/3);
t0 = w1 + a ./ (3*w1);
out(big) = t0(big).^2;
%}
% trigonometric method for three real roots:
% https://en.wikipedia.org/wiki/Cubic_function#Trigonometric_method_for_three_real_roots
t0 = 2 * sqrt(a/3) .* cos(1/3 * acos(3*(reg/2)/2./(-a) .* sqrt(3./a)));
a_1_2 = 3/2 * reg^(2/3);
big = a > a_1_2;
out(big) = t0(big).^2;
% plot(z, out, '-o', z, z, ':', z, z-reg, '--')
out = out .* sb;
end % ir_lpp_1_2_shrink

% ir_lpp_4_3_shrink()
% for l_p with p=4/3
% from eqn. (4.5) in chaux:07:avf
function out = ir_lpp_4_3_shrink(b, reg)
sb = sign(b); a = abs(b);
x = sqrt(a.^2 + 256 * reg^3 / 729);
out = sb .* (a + 4 * reg / (3 * 2^1/3) * ((x - a).^1/3 - (x + a).^1/3));
end % ir_lpp_4_3_shrink

% ir_tav_shrink()
function out = ir_tav_shrink(z, reg, delta)
out = zeros(size(z), class(z));
big = reg < abs(z) & abs(z) < reg + delta;
out(big) = z(big) .* (1 - reg ./ abs(z(big)));
big = reg + delta <= abs(z);
out(big) = z(big);
end % ir_tav_shrink

% ir_li98cfs_potk()
function pot = ir_li98cfs_potk(z, d)
pot = d.^2 .* ((z ./ d) .* atan(z ./ d) - 0.5 * log(1 + (z ./ d).^2));
end % ir_li98cfs_potk

% ir_genhub_potk()
function pot = ir_genhub_potk(z, d, p, q)
pot = 0.5 * abs(z) .^ p .* (abs(z) <= d) ...
 + 0.5 * (p ./ q .* d .^ (p-q) .* abs(z) .^ q ...
 + (1 - p ./ q) .* d .^ p) .* (abs(z) > d);
end % ir_genhub_potk

% ir_genhub_wpot()
function pot = ir_genhub_wpot(z, d, p, q)
pot = p / 2 .* (d .^ (p-q)) .* (abs(z) .^ (q-2)) .* (abs(z) > d);
ii = abs(z) <= d;
pot(ii) = p / 2 .* (abs(z(ii)) .^ (p-2));
%pot = p / 2 .* abs(t) .^ (p-2) .* (abs(t) <= d) ...
% + p / 2 .* d .^ (p-q) .* abs(t) .^ (q-2) .* (abs(t) > d);
end % ir_genhub_wpot

% ir_stevenson94dpr_potk()
function pot = ir_stevenson94dpr_potk(z, d, p, q)
% pr [d p q]
tmp1 = (0.5 * abs(z) .^ p) .* (abs(z) <= d);
tmp2 = 0.5 * ( (p .* (d .^ (p-1)) .* abs(z) - p .* (d .^ p) ...
 + (1 ./ q) .^ (1 ./ (q-1)) ) .^ q ...
 + d .^ p - (1 ./ q) .^ (q ./ (q-1)) ) .* (abs(z) > d);
pot = tmp1 + tmp2;
end % ir_stevenson94dpr_potk

% ir_tonumeric()
% convert x to type of y
function out = ir_tonumeric(x, y)
switch class(y)
case 'double'
	out = double(x);
case 'single'
	out = single(x);
otherwise
	fail('unknown type %s', class(y))
end
end % ir_tonumeric

% ir_potential_fun_plot()
function dummy = ir_potential_fun_plot(pot)
z = linspace(-1,1,101)*2;
if isvar('pot.delta')
	z = z * pot.delta;
end
if ~im, return, end
im plc 2 2
im subplot 1
plot(z, pot.potk(z), '-', z, z.^2/2, ':', z, abs(z), ':')
axis([min(z) max(z) minmax(pot.potk(z))'])

im subplot 2
plot(z, pot.dpot(z), '-', z, z, ':')
axis([min(z) max(z) minmax(pot.dpot(z))'])

im subplot 3
plot(z, pot.wpot(z), '-', z, 1+0*z, ':')
axis([min(z) max(z) 0 1])

im subplot 4
reg = 1;
tmp = pot.shrink(z, reg);
plot(z, tmp, '-', z, z, ':')
axis([min(z) max(z) minmax(tmp)'])

dummy = [];
end % ir_potential_fun_plot

% ir_potential_fun_test()
% test routine
% examine potential functions after rescaling.
function ir_potential_fun_test

delta = 10; zmax = 4 * delta; reg = 1.5 * delta;
plist = potential_fun('list');
%plist = {'l1', 'fair-l1'}; delta = 0.5;
%plist = {'quad', 'li98cfs', 'hyper3', 'huber2'}; % show li98cfs roughly hyper3
%plist = {'quad', 'genhub', 'huber', 'stevenson94dpr'};
%plist = {'genhub'}
%plist = {'hyper3', 'qgg2', 'huber'}; delta = 10; zmax = 50;
%plist = {'lange3', 'qgg2'}; zmax = 200;
%plist = {'qgg2', 'gf1-fit'};
%plist = potential_fun('list1'); delta = 0.5;
%plist = {plist{:}, 'quad', 'gf1', 'huber', 'broken'}; % todo: more!
%plist = {'qgg2', 'table1', 'table0'};
%plist = {'qgg2', 'table1'}; fname = 'fig_reg_pot_table1_qgg2';
%plist = {'qgg2', 'table0'}; fname = 'fig_reg_pot_table0_qgg2';
%plist = {'qgg2', 'table1', 'gf1'}; fname = 'fig_reg_pot_table0_gf1_qgg2';
%plist = {'cauchy', 'l1'};
zz = zmax * linspace(-1, 1, 2001)';
bb = 3 * max(delta + reg, delta * (1+reg)) * linspace(-1, 1, 301)';
ps = [];
lshrink = {};
for ii=1:numel(plist)
	ptype = plist{ii}; % pr ptype
	if streq(ptype, 'quad')
		leg{ii} = ptype;
	else
		leg{ii} = [ptype sprintf(', $\\delta = %g$', delta)];
	end

	switch ptype
	case {'gf1', 'gf1-fit'}
		param = potential_fun('gf1-fit', nan, {'qgg2', [1 10], 1.2});
		ptype = 'gf1'; % trick
		leg{ii} = [leg{ii} sprintf(' %.3g %.4f', param(1), param(2))];
	case 'qgg2'
		param = 1.2;
		leg{ii} = [leg{ii} sprintf(', $q = %g$', param)];
	case 'genhub'
		param = [2.0 1.2];
		leg{ii} = [leg{ii} sprintf(', $p=%g, q=%g$', param(1), param(2))];
	case 'stevenson94dpr'
		param = [2 1.2];
		leg{ii} = [leg{ii} sprintf(', $p=%g, q=%g$', param(1), param(2))];
	case {'table0', 'table1'}
		dz = delta / 20;
		tmp = potential_fun('qgg2', delta, 1.2);
%		param = {dz, tmp.dpot([1:1e4] * dz)};
		param = [dz, tmp.dpot([1:1e4] * dz)];
%		tmp = sprintf(' K = %d dz = %g', numel(param{2}), param{1});
		tmp = sprintf(', $K = %d$, $\\Delta t = %g$', numel(param)-1, param(1));
		leg{ii} = [leg{ii} tmp];
	otherwise
		param = [];
	end

	pot = potential_fun(ptype, delta, param);
	pp(:,ii) = pot.potk(zz);
	pw(:,ii) = pot.wpot(zz);
	pd(:,ii) = pot.dpot(zz);

	% replace 0 with 1 to make (slow) figure showing fzero-based shrinkage
	if 0 || ~streq(func2str(pot.meth.shrink), 'ir_potential_fun_shrink')
		tmp = pot.shrink(bb, reg);
		if any(tmp ~= 0)
			ps(:,end+1) = pot.shrink(bb, reg);
			lshrink{end+1} = leg{ii};
		end
	end

	if 0 % test vs old
		try
			opot = potential_func(ptype, delta, param);
			opp = opot.potk(opot, zz);
			opw = opot.wpot(opot, zz);
			opd = opot.dpot(opot, zz);
		catch
			printm('%s no old version', ptype)
		%	keyboard % problem with private scope of ir_li98cfs_wpot
			continue
		end
		if ~isequal(opp, pp(:,ii)), 'p bug', ptype, keyboard, end
		if ~isequal(opw, pw(:,ii)), 'w bug', ptype, keyboard, end
		if ~isequal(opd, pd(:,ii)), 'd bug', ptype, keyboard, end
	end
end

if im
%	set(0,'DefaultAxesLineStyleOrder', '-|:')
	clf, pl = @(i) subplot(410 + i);
	pl(1), plot(zz, pp), title 'potk'
	axis tight, axisy([-0.0 2.5] * delta^2)
	ir_legend(leg, 'location', 'north')
	pl(2), plot(zz, pw), title 'wpot'
	axis tight, axisy(0, 1.1), ytick([0 1])
	pl(3), plot(zz, pd), title 'dpot'
	axis tight, axisy([-1 1] * 1.3 * delta)
	xlabelf '$z$'
%	ir_savefig eps_c fname
	% check derivatives
	pl(4)
	plot(zz, pd)
	title 'dpot check'
	hold on
	d = diffc(pp) / (zz(2)-zz(1));
	plot(zz(1:end-1), d(1:end-1,:), '--', ...
		zz(1:end-1), d(1:end-1,:)-pd(1:end-1,:), ':')
	hold off
	axis tight, axisy([-1 1] * 1.3 * delta)
end

if 1 && im
	if ~isempty(lshrink)
		prompt
		clf
		plot(bb, ps, '-', bb, bb, '-')
		axis equal, axis square
		axis([-1 1 -1 1]*400)
		xtick([-1 0 1] * 400)
		ytick([-1 0 1] * 400)
		ir_legend(lshrink, 'location', 'southeast')
		xlabelf '$c$'
%		ylabelf 'xhat(z)'
		ylabelf '$\hat{z}(c)$'
		grid
%		ir_savefig cw fig_reg_pot_table0_qgg2_shrink

		if 0 % figure for book
			plot(bb, ps(:,[3 2]) - repmat(ps(:,1), [1 2]), 'o')
			plot(	bb, ps(:,3) - ps(:,1), 'bo', ...
				bb, ps(:,2) - ps(:,1), 'gx')
			axis([0 500 -0.05 0.5]), ytick([0 0.05 0.1 0.5])
			xtick([0 500])
			ir_fontsize label 18
			ir_fontsize text 18
			xlabelf 'c', ylabelf 'shrinkage error for QGG2'
			ir_legend(lshrink([3 2]), 1)
%			ir_savefig cw fig_reg_pot_table01_qgg2_shrinker
			keyboard
		end
	end
end

% check 'scale' option
pot1 = potential_fun('lange3', 10, []);
pot2 = potential_fun('lange3', 10, [], 'scale', 2);
jf_equal(2 * pot1.potk(zz), pot2.potk(zz))
jf_equal(2 * pot1.dpot(zz), pot2.dpot(zz))
jf_equal(2 * pot1.wpot(zz), pot2.wpot(zz))

end % ir_potential_fun_test

% ir_potential_fun_test_lpp()
function ir_potential_fun_test_lpp
params = [0 0.01 0.1 1/2 3/4 1 4/3];
params = [0 1/2 1];
reg = 4*2^4; zmax = 2 * reg;
zz = zmax * linspace(-1, 1, 2001)';
bb = 2 * reg * linspace(-1, 1, 201);
%a_1_2 = 3 * (reg/4)^(2/3) % saddle point for p=1/2
a_1_2 = 3/2 * reg^(2/3); % breakpoint for p=1/2
if 0
	a = a_1_2;
	pot = potential_fun('lpp', nan, 1/2);
	z0 = ir_potential_fun_shrink(pot, a+1e-9, reg)
	cost = @(z) 1/2 * (z - a).^2 + reg * pot.potk(z);
	t0 = sqrt(z0)
	t0^3 - 2*a*t0 + 2*reg % should be 0
	t0^3 - a*t0 + reg/2 % should be 0
	3/2*reg - a * t0 % should be 0
	(3^(3/2) * reg)/(4*a^(3/2)) % should be 1/sqrt(2)
	cost(0)
	cost(z0)
%	syms x; solve(3 * acos(x) - acos(-x), x) % 1/sqrt(2)
return
end
tmp = [sqrt(2*reg) reg a_1_2];
tmp = outer_sum(tmp, [-1 +1]*1e-5);
bb = sort([bb col(tmp)' -col(tmp)'])'; % trick
ps = [];
lshrink = {};
if 0
	pot = potential_fun('lpp', nan, 1/2);
	potk = @(z) pot.potk(z);

%	y = sqrt(2*reg); % l0
%	y = reg; % l1
	y = a_1_2 + 1e-9; % p=1/2
	if 0
		ir_potential_fun_shrink(pot, y, reg)
		ir_potential_fun_shrink(pot, y, reg)
	return
	end

	if 0
		z1 = pot.shrink(y, reg)
		d1 = z1 - y + reg * pot.dpot(z1) % 1st derivative
		d2 = 1 + reg * (1/2)*(-1/2) * z1^(-3/2) % 2nd derivative > 0?

		dfit = 1/2 * (y - zz).^2;
		cost = dfit + reg * potk(zz);
		plot(zz, dfit, '--', zz, reg*potk(zz), '-', zz, cost, '-')
		tmp = minmax(zz(cost <= 3*min(cost)));
		tmp(1) = min(tmp(1), -10);
		tmp(2) = max(tmp(2), reg/3);
		xlim(tmp), ylim([0.0 2*1.2] * min(cost)), grid
	%	ytick([0 round(min(cost))])
		ytick([0 min(cost)])
		xtick([0 z1 a_1_2]), xlabelf '$z$', ylabelf '$\Psi(z)$'
	return
	end

	zs = ir_potential_fun_shrink(pot, bb, reg);
	zp = pot.shrink(bb, reg);

	plot(bb, zs, '-o', bb, zp, '.-', bb, bb, ':')
	legend({'fmin', 'trig'}, 'location', 'northwest')
	xlabelf '$b$'
	ylabelf '$\hat{z}(b)$'
	axis equal, axis square
	axis([-1 1 -1 1]*1.1*reg)
	tmp = [sqrt(2*reg) a_1_2 reg]; % l0 and l1 breakpoints
	xtick([-tmp 0 tmp])
	ytick([-1 0 1] * reg)
	grid
%	ir_savefig cw fig_?
return
end

for ip=1:length(params)
	pot = potential_fun('lpp', 0, params(ip));

	% replace 0 with 1 to make (slow) figure showing fzero-based shrinkage
	if 0 || ~streq(func2str(pot.meth.shrink), 'ir_potential_fun_shrink')
		tmp = pot.shrink(bb, reg);
		if any(tmp ~= 0)
			ps(:,end+1) = pot.shrink(bb, reg);
			lshrink{end+1} = sprintf('p=%g', params(ip));
		end
	end

	clf
	plot(bb, ps, '-', bb, bb, ':')
	legend(lshrink{:}, 'location', 'southeast')
	xlabelf '$b$'
	ylabelf '$\hat{z}(b)$'
	axis equal, axis square
	axis([-1 1 -1 1]*1.5*reg)
	tmp = [sqrt(2*reg) a_1_2 reg]; % l0 and l1 breakpoints
	xtick([-tmp 0 tmp])
	ytick([-1 0 1] * reg)
	grid
%	ir_savefig cw fig_?
end

end % ir_potential_fun_test_lpp()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function xyz_ls=plot_ROI_3D(CONST,idlist)
% plot yx,yx image
% Note imagesc yx axis are mirrored to scatter3d
% change y to -y to match the two

marker_size=64;

z_ratio=19.7; % adjust scale between z and xy

Nid=length(idlist);
xyz_ls=zeros(3,Nid);

% for i=1:Nid
%     center_i=CONST.CIF(idlist(i));
%     xyz_ls(3,i)=center_i.slice*z_ratio;
%     xyz_ls(2,i)=center_i.center(1);
%     xyz_ls(1,i)=center_i.center(2);
% end;
CIFs=CONST.CIF(idlist);
xyz_ls(3,:)=[CIFs.slice]*z_ratio;
xyz_ls([2,1],:)=reshape([CIFs.center],2,[]);

marker_option={marker_size,'b','filled'};

scatter3(xyz_ls(1,:),-xyz_ls(2,:),xyz_ls(3,:),marker_option{:});
axis equal;

end

author: xiuyechen

repo: FishExplorer


text:
function suborder = triangle_nco_suborder ( rule, suborder_num )

%*****************************************************************************80
%
%% TRIANGLE_NCO_SUBORDER returns the suborders for an NCO rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Input, integer SUBORDER_NUM, the number of suborders of the rule.
%
%    Output, integer SUBORDER(SUBORDER_NUM), the suborders of the rule.
%
  if ( rule == 1 )
    suborder(1:suborder_num) = [ ...
      1 ];
  elseif ( rule == 2 )
    suborder(1:suborder_num) = [ ...
      3 ];
  elseif ( rule == 3 )
    suborder(1:suborder_num) = [ ...
      3, 3 ];
  elseif ( rule == 4 )
    suborder(1:suborder_num) = [ ...
      3, 6, 1 ];
  elseif ( rule == 5 )
    suborder(1:suborder_num) = [ ...
      3, 6, 3, 3 ];
  elseif ( rule == 6 )
    suborder(1:suborder_num) = [ ...
      3, 6, 6, 3, 3 ];
  elseif ( rule == 7 )
    suborder(1:suborder_num) = [ ...
      3, 6, 6, 3, 3, 6, 1 ];
  elseif ( rule == 8 )
    suborder(1:suborder_num) = [ ...
      3, 6, 6, 3, 6, 6, 3, 3 ];
  elseif ( rule == 9 )
    suborder(1:suborder_num) = [ ...
      3, 6, 6, 3, 6, 6, 3, 6, 3, 3  ];

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_NCO_SUBORDER - Fatal error!\n' );
    fprintf ( 1, '  Illegal RULE = %d\n', rule );
    error ( 'TRIANGLE_NCO_SUBORDER - Fatal error!' )

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hh = quiver(varargin)
%QUIVER Quiver plot.
%   QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
%   at the points (x,y).  The matrices X,Y,U,V must all be the same size
%   and contain corresponding position and velocity components (X and Y
%   can also be vectors to specify a uniform grid).  QUIVER automatically
%   scales the arrows to fit within the grid.
%
%   QUIVER(U,V) plots velocity vectors at equally spaced points in
%   the x-y plane.
%
%   QUIVER(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the
%   arrows to fit within the grid and then stretches them by S.  Use
%   S=0 to plot the arrows without the automatic scaling.
%
%   QUIVER(...,LINESPEC) uses the plot linestyle specified for
%   the velocity vectors.  Any marker in LINESPEC is drawn at the base
%   instead of an arrow on the tip.  Use a marker of '.' to specify
%   no marker at all.  See PLOT for other possibilities.
%
%   QUIVER(...,'filled') fills any markers specified.
%
%   H = QUIVER(...) returns a vector of line handles.
%
%   Example:
%      [x,y] = meshgrid(-2:.2:2,-1:.15:1);
%      z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
%      contour(x,y,z), hold on
%      quiver(x,y,px,py), hold off, axis image
%
%   See also FEATHER, QUIVER3, PLOT.

%   Clay M. Thompson 3-3-94
%   Copyright 1984-2002 The MathWorks, Inc.
%   $Revision: 2.1 $  $Date: 2005/05/23 16:20:12 $

% Arrow head parameters
alpha = 0.33; % Size of arrow head relative to the length of the vector
beta = 0.33;  % Width of the base of the arrow head relative to the length
autoscale = 1; % Autoscale if ~= 0 then scale by this.
plotarrows = 1; % Plot arrows
sym = '';

filled = 0;
ls = '-';
ms = '';
col = '';

nin = nargin;
% Parse the string inputs
while isstr(varargin{nin}),
  vv = varargin{nin};
  if ~isempty(vv) & strcmp(lower(vv(1)),'f')
    filled = 1;
    nin = nin-1;
  else
    [l,c,m,msg] = colstyle(vv);
    if ~isempty(msg),
      error(sprintf('Unknown option "%s".',vv));
    end
    if ~isempty(l), ls = l; end
    if ~isempty(c), col = c; end
    if ~isempty(m), ms = m; plotarrows = 0; end
    if isequal(m,'.'), ms = ''; end % Don't plot '.'
    nin = nin-1;
  end
end

error(nargchk(2,5,nin));

% Check numeric input arguments
if nin<4, % quiver(u,v) or quiver(u,v,s)
  [msg,x,y,u,v] = xyzchk(varargin{1:2});
else
  [msg,x,y,u,v] = xyzchk(varargin{1:4});
end
if ~isempty(msg), error(msg); end

if nin==3 | nin==5, % quiver(u,v,s) or quiver(x,y,u,v,s)
  autoscale = varargin{nin};
end

% Scalar expand u,v
if prod(size(u))==1, u = u(ones(size(x))); end
if prod(size(v))==1, v = v(ones(size(u))); end

if autoscale,
  % Base autoscale value on average spacing in the x and y
  % directions.  Estimate number of points in each direction as
  % either the size of the input arrays or the effective square
  % spacing if x and y are vectors.
  if min(size(x))==1, n=sqrt(prod(size(x))); m=n; else [m,n]=size(x); end
  delx = diff([min(x(:)) max(x(:))])/n;
  dely = diff([min(y(:)) max(y(:))])/m;
  del = delx.^2 + dely.^2;
  if del>0
    len = sqrt((u.^2 + v.^2)/del);
    maxlen = max(len(:));
  else
    maxlen = 0;
  end

  if maxlen>0
    autoscale = autoscale*0.9 / maxlen;
  else
    autoscale = autoscale*0.9;
  end
  u = u*autoscale; v = v*autoscale;
end

ax = newplot;
next = lower(get(ax,'NextPlot'));
hold_state = ishold;

% Make velocity vectors
x = x(:).'; y = y(:).';
u = u(:).'; v = v(:).';
uu = [x;x+u;repmat(NaN,size(u))];
vv = [y;y+v;repmat(NaN,size(u))];

h1 = plot(uu(:),vv(:),[col ls]);

if plotarrows,
  % Make arrow heads and plot them
  hu = [x+u-alpha*(u+beta*(v+eps));x+u; ...
        x+u-alpha*(u-beta*(v+eps));repmat(NaN,size(u))];
  hv = [y+v-alpha*(v-beta*(u+eps));y+v; ...
        y+v-alpha*(v+beta*(u+eps));repmat(NaN,size(v))];
  hold on
  h2 = plot(hu(:),hv(:),[col ls]);
else
  h2 = [];
end

if ~isempty(ms), % Plot marker on base
  hu = x; hv = y;
  hold on
  h3 = plot(hu(:),hv(:),[col ms]);
  if filled, set(h3,'markerfacecolor',get(h1,'color')); end
else
  h3 = [];
end

if ~hold_state, hold off, view(2); set(ax,'NextPlot',next); end

if nargout>0, hh = [h1;h2;h3]; end

author: strawlab

repo: MultiCamSelfCal


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% periodic signals

t=1:5;
x=sin(t);
T=2*pi;
for k=1:10
xk(k,:)=sin(t+k*T);
end
xk

%Sum of periodic continuous time signals
 t=0:.1:6*pi;
 x=cos(t)+sin(3*t);
 plot(t,x)

 %	Construction of periodic signals
 figure
 [s,t]=gensig('square',3,20,0.01);
 plot(t,s)
 ylim([-.3 1.3])

 figure
 [s,t]=gensig('pulse',2,10);
 plot(t,s)
 ylim([-.3 1.3])

 figure
 t=0:0.1:10;
 s=square(t);
 plot(t,s);
 ylim([-1.3 1.3])

 figure
 s2=square(2*pi*t);
 plot(t,s2);
 ylim([-1.3 1.3]) ;

 figure
 t=0:0.1:20;
 s=sawtooth(t);
 plot(t,s);
 ylim([-1.3 1.3])

 figure
 t=0:.1:10;
 x=t.*exp(-t);
 xp=repmat(x,1,8);
 tp=linspace(0,80,length(xp));
 plot(tp,xp)

 figure
 N=10;
 x=[1 1 -1 -1];
 xp=repmat(x,1,N);
 n=0:length(xp)-1;
 stem(n,xp);


text:
% Test size bounds with varying mesh gradation rates.
clearvars; clc;

addpath('..')
addpath(genpath('../utilities/'))
addpath(genpath('../datasets/'))
addpath(genpath('../m_map/'))

RESO_TOL = 95; %percentage of resolution in bounds

bbox = [166 176;		% lon_min lon_max
    -48 -40]; 		% lat_min lat_max
min_el    = 1e3;  		% minimum resolution in meters.
max_el    = 100e3; 		% maximum resolution in meters.
max_el_ns = 5e3;        % maximum resolution nearshore in meters.
grade     = [0.15; 0.25; 0.35]; 		% mesh grade in decimal percent.
R         = 3;    		% number of elements to resolve feature width.
coastline = 'GSHHS_f_L1';
gdat = geodata('shp',coastline,'bbox',bbox,'h0',min_el);

for i = 1 : 3 % for each grade
    fh = edgefx('geodata',gdat,...
        'fs',R,'max_el_ns',max_el_ns,...
        'max_el',max_el,'g',grade(i));
    mshopts = meshgen('ef',fh,'bou',gdat,'plot_on',0,'nscreen',5,'proj','trans');
    mshopts = mshopts.build;
    m1 = mshopts.grd;

    [bars,barlen] = GetBarLengths(m1,0);
    % sort bar lengths in descending order
    [barlen,IA] = sort(barlen,'descend');
    bars = bars(IA,:);
    % get the minimum bar length for each node
    [B1,IB] = unique(bars(:,1),'last');
    [B2,IC] = unique(bars(:,2),'last');
    d1 = NaN*m1.p(:,1); d2 = NaN*m1.p(:,1);
    d1(B1) = barlen(IB); d2(B2) = barlen(IC);
    reso = min(d1,d2);

    reso_in_bounds = 100*sum(reso > min_el & reso < max_el)/length(reso);
    if reso_in_bounds < RESO_TOL
        error(['Resolution bounds does not match for grade ' num2str(grade(i)) ...
               '. Got ' num2str(reso_in_bounds) '% of vertices with resolution in bounds']);
        exit(1)
    end
    disp(['Passed for ' num2str(grade(i)) '. ' ...
          num2str(reso_in_bounds) '% of vertices have resolution in bounds']);
end

author: CHLNDDEV


text:
function r8vec_split_test ( )

%*****************************************************************************80
%
%% R8VEC_SPLIT_TEST tests R8VEC_SPLIT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 25;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_SPLIT_TEST\n' );
  fprintf ( 1, '  R8VEC_SPLIT splits a vector into\n' );
  fprintf ( 1, '  entries less than and greater than a\n' );
  fprintf ( 1, '  splitting value.\n' );

  b = 0.0;
  c = 10.0;
  seed = 123456789;

  [ a, seed ] = r8vec_uniform_ab ( n, b, c, seed );

  a(1:n) = round ( a(1:n) ) / 2.0;

  split = 0.5 * ( a(1) + a(n) );

  r8vec_print ( n, a, '  The array:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Splitting value is %f\n', split );

  [ a, isplit ] = r8vec_split ( n, a, split );

  r8vec_print ( n, a, '  The split array:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Array entries <= SPLIT up to index %d\n', isplit );

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: finds coefficients of cubic spline interpolant
%
% INPUTs: (p1,p2): interpolating mediary points
%
% Author: Nick Battista
% Date: February 23, 2018
% Institution: TCNJ
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function interp_Function_Coeffs(p1,p2)

%
% Set up MATRIX for linear system of cubic interpolant
%
mat1 = [1 0 0 0 0 0 0 0 0 0 0 0];
mat2 = [0 1 0 0 0 0 0 0 0 0 0 0];
mat3 = [0 0 1 0 0 0 0 0 0 0 0 0];
mat4 = [1 p1 p1^2 p1^3 -1 -p1 -p1^2 -p1^3 0 0 0 0];

mat5 = [0 1 2*p1 3*p1^2 0 -1 -2*p1 -3*p1^2 0 0 0 0];
mat6 = [0 0 2 6*p1 0 0 -2 -6*p1 0 0 0 0];
mat7 = [0 0 0 0 1 p2 p2^2 p2^3 -1 -p2 -p2^2 -p2^3];
mat8 = [0 0 0 0 0 1 2*p2 3*p2^2 0 -1 -2*p2 -3*p2^2];

mat9 = [0 0 0 0 0 0 2 6*p2 0 0 -2 -6*p2];
mat10= [0 0 0 0 0 0 0 0 1 1 1 1];
mat11= [0 0 0 0 0 0 0 0 0 1 2 3];
mat12= [0 0 0 0 0 0 0 0 0 0 2 6];

mat = [mat1; mat2; mat3; mat4; mat5; mat6; mat7; mat8; mat9; mat10; mat11; mat12];

%
% RHS of Linear System for cubic interpolant
%
rhs = [0 0 0 0 0 0 0 0 0 1 0 0]';

%
% Give coefficients
%
coeffs = mat\rhs

%------------------------------------
% Plot Interpolation Polynomial
%               (and its derivatives)
%------------------------------------

%Plotting t-Values
ds = 0.005;
t = 0:ds:1;
t2=t/2;

%Plot Attributes
ms = 20;

%
% PLOT: interpolating function, g(t)
figure(1)
for i=1:length(t)
    plot(t2(i),g(coeffs,t(i),p1,p2),'.','MarkerSize',ms); hold on;
end
xlabel('t');
ylabel('g(t)');
set(gca,'FontSize',18)

%
% PLOT: interpolating function's derivative, dg(t)/dt
figure(2)
for i=1:length(t)
    plot(t2(i),gP(coeffs,t(i),p1,p2),'.','MarkerSize',ms); hold on;
end
xlabel('t');
ylabel('dg/dt');
set(gca,'FontSize',18)

%
% PLOT: interpolating function's 2nd derivative, d^2g(t)/dt^2
figure(3)
for i=1:length(t)
    plot(t2(i),gPP(coeffs,t(i),p1,p2),'.','MarkerSize',ms); hold on;
end
xlabel('t');
ylabel('d^2g/dt^2');
set(gca,'FontSize',18)

%
% Store coefficients for plotting
%
% vec = [];
% for i=1:12
%     vec = [vec coeffs(i)];
% end

%strName = 'coeffs_Vec';
%print_Matrix_To_Txt_File(coeffs,strName);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION prints matrix to file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Matrix_To_Txt_File(a,strName)

nameTxt = [strName '.txt'];

fid = fopen(nameTxt, 'wt'); % Open for writing
for i=1:size(a,1)
   fprintf(fid, '%.14f ', a(i,:));
   fprintf(fid, '\n');
end
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: evaluate interpolating polynomial
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = g(coeffs,x,p1,p2)

if x<=p1
    val = coeffs(1) + coeffs(2)*x + coeffs(3)*x^2 + coeffs(4)*x^3;
elseif x<=p2
    val = coeffs(5) + coeffs(6)*x + coeffs(7)*x^2 + coeffs(8)*x^3;
else
    val = coeffs(9) + coeffs(10)*x + coeffs(11)*x^2 + coeffs(12)*x^3;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: evaluate interpolating polynomial 1st derivative
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = gP(coeffs,x,p1,p2)

if x<=p1
    val = coeffs(2) + 2*coeffs(3)*x + 3*coeffs(4)*x^2;
elseif x<=p2
    val = coeffs(6) + 2*coeffs(7)*x + 3*coeffs(8)*x^2;
else
    val = coeffs(10) + 2*coeffs(11)*x + 3*coeffs(12)*x^2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: evaluate interpolating polynomial 2nd derivative
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = gPP(coeffs,x,p1,p2)

if x<=p1
    val = 2*coeffs(3) + 6*coeffs(4)*x;
elseif x<=p2
    val = 2*coeffs(7) + 6*coeffs(8)*x;
else
    val = 2*coeffs(11) + 6*coeffs(12)*x;
end

author: nickabattista


text:
function Iout= ImRegular(Iin)

         m = size(Iin, 1);
         n = size(Iin, 2);
         Imin = min(Iin(:));
         Iout = Iin;

         if Imin<0
             a = -Imin;
             for j = 1:m
                 for i = 1:n
                     if Iin(j,i) < a
                         Iout(j,i) = (Iin(j,i)-Imin)^2/(4*a);
                     end
                 end
             end
         end

         Imax = max(Iout(:));
         if Imax>255
             b = 510-Imax;
             for j = 1:m
                 for i = 1:n
                     if Iout(j,i) > b
                         Iout(j,i) = (Iout(j,i)-Imax)^2/(4*(b-255)) + 255;
                     end
                 end
             end
         end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function s = s_low ( s )

%*****************************************************************************80
%
%% S_LOW replaces all uppercase letters by lowercase ones.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be transformed.
%
%    Output, string S, a lowercase copy of the input string.
%
  nchar = s_len_trim ( s );

  for i = 1 : nchar
    s(i) = ch_low ( s(i) );
  end

  return
end


text:
function triangulation_test17 ( )

%*****************************************************************************80
%
%% TEST17 tests TRIANGULATION_ORDER3_EXAMPLE1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST17\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_EXAMPLE1_SIZE gives the sizes\n' );
  fprintf ( 1, '    for an example triangulation;\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_EXAMPLE1 returns the information\n' );
  fprintf ( 1, '    for an example triangulation;\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_PRINT prints a triangulation.\n' );
%
%  Get the sizes.
%
  [ node_num, triangle_num, hole_num ] = ...
    triangulation_order3_example1_size ( );
%
%  Get the data.
%
  [ node_xy, triangle_node, triangle_neighbor ] = ...
    triangulation_order3_example1 ( node_num, triangle_num );

  triangulation_order3_print ( node_num, triangle_num, node_xy, ...
    triangle_node, triangle_neighbor );

  return
end


text:
function mean = half_normal_mean ( a, b )

%*****************************************************************************80
%
%% HALF_NORMAL_MEAN returns the mean of the Half Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a + b * sqrt ( 2.0 / pi );

  return
end


text:
function sc = zncc(I,J,SI,SII,SJ,SJJ)

% ZNCC  Zero mean, Normalized Correlation Coefficient
%   ZNCC(I,J) computes the zncc score of matrices I and J:
%
%               1     sum( sum ( (I-mean(I)) * (J-mean(J)) ) )
%       ZNCC = --- * ------------------------------------------
%               N               std(I) * std(J)
%
%   where N = prod(size(I))
%   and   size(I) = size(J)
%
%   ZNCC(I,J,SI,SII,SJ,SJJ) accepts useful intermediate results
%   that permit to speed up the calculations. These are:
%
%     SI  = sum(sum(I))
%     SII = sum(sum(I.*I))
%     SJ  = sum(sum(J))
%     SJJ = sum(sum(J.*J))
%
%   See also PATCHCORR, SSD, CENSUS

switch nargin
    case 6
        SIJ = sum(sum(I.*J));
    case 5
        SJJ = sum(sum(I.*I));
        SIJ = sum(sum(I.*J));
    case 4
        SJ  = sum(sum(J));
        SJJ = sum(sum(J.*J));
        SIJ = sum(sum(I.*J));
    case 3
        SII = sum(sum(I.*I));
        SJ  = sum(sum(J));
        SJJ = sum(sum(J.*J));
        SIJ = sum(sum(I.*J));
    case {1,2}
        SI  = sum(sum(I));
        SII = sum(sum(I.*I));
        SJ  = sum(sum(J));
        SJJ = sum(sum(J.*J));
        SIJ = sum(sum(I.*J));
end

N = numel(I);

% This is the equivalent, faster formula than that given
% in the help section:
%
%   sc = (N*SIJ - SI*SJ)/sqrt((N*SII - SI^2)*(N*SJJ - SJ^2)+eps);
%
% we make a test on the variance of I and J to avoid degenerated cases

den = (N*SII - SI^2)*(N*SJJ - SJ^2);

if den < N
%     figure(5);image(J);colormap(gray(255));axis image
    sc = 0;
else
    sc = (N*SIJ - SI*SJ)/sqrt(den);
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function selectionIndexArrayCell = tapas_uniqc_convert_selection_range_to_array(selectionIndexRangeCell)
% Converts selection name/range pairs to cells of single selections
%
%   selectionIndexArrayCell = tapas_uniqc_convert_selection_range_to_array(selectionIndexRangeCell)
%
% IN
%   selectionIndexRangeCell     cell(1,2*dimLabels) of dimLabel /
%                               dimValueRange pairs,
%                               e.g., {'coils', 1:8, 'echo', 1:3}
%
% OUT
%   selectionIndexArrayCell     cell(nValuesDim1,...,nValuesDim1) of
%                               dimLabel / dimValue pairs as used in
%                               MrDimInfo.split (selectionArray)
%                               e.g.,
%                               {'coils', 1, 'echo', 1}, ..., {'coils', 1, 'echo', 3}
%                               ...
%                               {'coils', 8, 'echo', 1}, ..., {'coils', 8, 'echo', 3}
%
% EXAMPLE
%   selectionIndexArrayCell = tapas_uniqc_convert_selection_range_to_array(...
%       {'coils', 1:8, 'echo', 1:3});
%   selectionIndexRangeCell = tapas_uniqc_convert_selection_array_to_range(...
%       selectionIndexArrayCell); % should return {'coils', 1:8, 'echo', 1:3}
%
%   See also tapas_uniqc_convert_selection_array_to_range MrDimInfo.split

% Author:   Lars Kasper
% Created:  2018-05-04
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

indexRanges = selectionIndexRangeCell(2:2:end);

nElementsPerDim = cellfun(@numel, indexRanges);
nDims = numel(nElementsPerDim);

selectionIndexArrayCell = cell(nElementsPerDim);

indexGrid = cell(1,nDims);
[indexGrid{:}] = ndgrid(indexRanges{:});

nElements = prod(nElementsPerDim);
for iElement = 1:nElements
    selectionIndexArrayCell{iElement} = cell(1,2*nDims);
    for iDim = 1:nDims
        selectionIndexArrayCell{iElement}{2*iDim-1} = ...
            selectionIndexRangeCell{2*iDim-1};
        selectionIndexArrayCell{iElement}{2*iDim} = ...
            indexGrid{iDim}(iElement);
    end
end

author: translationalneuromodeling

repo: tapas


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
Q3 = -pi:0.1:pi;
for j=1:numcols(Q3);
    M = p560.inertia([0 0 Q3(j) 0 0 0]);
    M22(j) = M(2,2);
end
plot(Q3, M22)
xlabel('q_3 (rad)');
ylabel('M_{22}');

author: Allopart

repo: rbpf-gmapping


text:
% pop_compareerps() - Compare the (ERP) averages of two datasets.
%
% Usage:
%       >> pop_compareerps( ALLEEG, datasetlist, chansubset, title);
% Inputs:
%   ALLEEG      - array of datasets
%   datasetlist - list of datasets
%   chansubset  - vector of channel subset
%   title       - plot title
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eeglab(), plottopo()

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

% 01-25-02 reformated help & license -ad
% 03-11-02 added empty ALLEEG check -ad
% 03-18-02 added channel subset -ad
% 03-18-02 added title -ad & sm
function com = pop_compareerps( ALLEEG, setlist, chansubset, plottitle);

com = '';
if nargin < 1
   help pop_compareerps;
   return;
end;
if isempty(ALLEEG)
    error('pop_compareerps: cannot process empty sets of data');
end;
if exist('plottitle') ~= 1
	plottitle = '';
end;
if nargin < 2
   % which set to save
	% -----------------
	promptstr    = { 'List of datasets to compare (ex: 1 3 4):' ...
	                 'Channels subset to consider ([]=all):' ...
	                 'Plot title ([]=automatic):' };
	inistr       = { '1' '' '' };
	result       = inputdlg2( promptstr, 'Compare dataset ERPs -- pop_compareerps()', 1,  inistr, 'pop_compareerps');
	if length(result) == 0 return; end;
	setlist   	 = eval( [ '[' result{1} ']' ] );
	chansubset   = eval( [ '[' result{2} ']' ] );
	if isempty( chansubset ), chansubset = 1:ALLEEG(setlist(1)).nbchan; end;
	plottitle    = result{3};
    if isempty(plottitle)
        plottitle = [ 'Compare datasets number (blue=first; red=sec.)' int2str(setlist) ];
    end;
    figure;
end;
try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
tracing = [];
for setindex = setlist
	tracing  = [ tracing squeeze(mean(ALLEEG(setindex).data,3))];
end;

% save channel names
% ------------------
plottopo( tracing, ALLEEG(setlist(1)).chanlocs, ALLEEG(setlist(1)).pnts, [ALLEEG(setlist(1)).xmin ALLEEG(setlist(1)).xmax 0 0]*1000, plottitle, chansubset );

com = sprintf('figure; pop_compareerps( %s, [%s], [%s], ''%s'');', inputname(1), num2str(setlist), num2str(chansubset), plottitle);
return;

author: PatternRecognition

repo: OpenBMI


text:
% This file will generate the necessary LHS library to identify the PDE.
% Please modify this file to code the library you want.
%
% Last Updated: 2019/06/19
% Coded By: K

function [Data,Sym_Struct]=LHS_Guess_PDE_Ex1(U,Ut,Utt,Ux,Uxx,Uxxx)
%% First get the size of the u vector.
[Data_Length,~]=size(U);

%Also create the symbolic variable
syms u ut utt ux uxx uxxx

Data=[];
Index=1;

Org_Data=[U Ut Utt Ux Uxx Uxxx];
Org_Sym=[u ut utt ux uxx uxxx];

for i=2:size(Org_Data,2)
    Data(:,Index)=Org_Data(:,1).*Org_Data(:,i);
    Sym_Struct{1,Index}=Org_Sym(1,1)*Org_Sym(1,i);
    Index=Index+1;
end


text:
function u_est = CG_MultiScale_LAP3D_Fast(I1, I2, FilterSizes, PreFilt, MedFilt, uin)
% The function implements a fast multi-scale framework for the 3D version
% of the LAP optical flow algorithm. Instead of downsampling the images, the
% framework changes the size of the all-pass filters used in the LAP
% algorithm. The filter basis used in the LAP algorithm spans the
% derivatives of a Gaussian filter.
%
% Function uses faster interpolation (but less accurate).
%
% Note that this implementation is for greyscale images only.
%
% Input parameters:
%       I1 and I2   -> Input images of size M by N (Greyscale)
%       FilterSizes -> Controls the size of the filters (vector with
%                      values for R
%                      i.e. filter size is [2R + 1] by [2R + 1]
%       PreFilt     -> Optional parameter to decide whether to highpass
%                       filter the images. (Default = 1 (Yes), 0 = (No))
%       MedFilt     -> Optional parameter to decide whether to median
%                       filter the flow at fine scales. (Default = 1 (Yes), 0 = (No))
%
% Outputs:
%       u_est       -> Estimate of the optical flow
%
%
% Date: 08/07/2015, Author: Chris Gilliam
%
% modifications for MRI application:
% - different input parameter: Filter size instead of upper parameter R
% - initialization in single precision (except filter kernel)
% 26.08.2015, Verena Neumann/Thomas Kuestner

if nargin <= 3,
    PreFilt = 1;
    MedFilt = 1;
elseif nargin <= 4,
    MedFilt = 1;
end

% Obtain the dimensions of the images:
[M, N, P] = size(I1);

% Initialise filter functions:
funs = Filter_Functions;

% Local Gaussian filter (used in high pass filtering)
h = exp(-(-2:2).^2/2);
h = h./sum(h(:));

% Initial optical flow estimate
if(nargin < 6)
    u_holder = repmat({zeros(M,N,P,'single')}, 1, 3);
else
    u_holder = uin;
end

% define half support of filters (i.e. R)
amp_array = FilterSizes;

% Initialise local counter:
num_level = 0;

% Local I1 variable:
if PreFilt == 1,
    im1 = I1 - funs.Filter_General(I1, h); % High pass filtering using gaussian filter h
%       im1 = funs.Filter_Laplacian(im1);  % high-pass using generalised Laplacian filter
else
    im1 = I1;
end

% Start estimating the optical flow
for l = 1:length(amp_array),
    num_level = num_level + 1;
    disp(['Level ', int2str(num_level), '/', int2str(size(amp_array,2))]);

    % Define filter parameter R at each iteration
    amp_size = amp_array(l);

    % Load Filter Basis:
    Basis_Set = loadbasis(3,amp_size);

    tic;
    if l == 1,
        % No warping for first iteration:
        if PreFilt == 1,
            I2_shift = I2 - funs.Filter_General(I2, h); % High pass filtering using gaussian filter h
%             I2_shift = funs.Filter_Laplacian(I2);       % high-pass using generalised Laplacian filter
        else
            I2_shift = I2;
        end
    else
        % Warp I2 closer to I1 using current optical flow estimate
%         I2_shift = imshift_3D(I2,{-u_holder{1}, -u_holder{2}, -u_holder{3}}, 'shiftedlinear');
        I2_shift = ShiftedLinear_Interp_3D(-u_holder{1},-u_holder{2},-u_holder{3},I2);
        if PreFilt == 1,
            I2_shift = I2_shift - funs.Filter_General(I2_shift, h); % High pass filtering using gaussian filter h
%           I2_shift = funs.Filter_Laplacian(I2_shift); % high-pass using generalised Laplacian filter
        end
    end

    disp(['Filter size = ', num2str(2*amp_size+1), ', Image Warping time = ', num2str(toc,3)]);

    % Using basis functions estimate optical flow on a local scale:
    tic;
    [uest_Orig, ~] = optiflowFilter3D(im1, I2_shift, Basis_Set);
    disp(['Filter size = ', num2str(2*amp_size+1), ', LAP Algorithm time = ', num2str(toc,3)]);

    % Clean optical flow
    % Stage 1: Remove nan's in the optical flow using inpainting:
    if sum(isnan(uest_Orig{1}(:))) >= M*N*P,
        error('All NaN. Suggest reduce Level_Num');
    end
    tic;
    [uest_Clean, ~] = cleanOF3D(uest_Orig);

    % Stage 2: Remove flow elements that corresponding to large warping
    % errors.
    R = round(2*amp_size);
    k1 = -R:R;
    uest_Clean = cellfun(@funs.Filter_Gauss, uest_Clean, {R,R,R}, {k1,k1,k1}, 'uni',false);

    disp(['Filter size = ', num2str(2*amp_size+1), ', Post-Processing time = ', num2str(toc,3)]);

    % Rescale optical flow and add to estimate
    u_holder = cellfun(@plus,u_holder,uest_Clean,'uniformoutput',false);

    % Refinement of OF at highest level based on errors in shifted
    %           images
    if amp_size <= 2 && MedFilt == 1,
        u_holder = Cleaning_Procedure(u_holder);
    end

end

u_est = u_holder;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Embedded functions %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u_out = Cleaning_Procedure(u_in)
% function cleans the estimate of the optical flow. The cleaning process
% comprises a robust smoothing stage using two Median filters (of different
% sizes)
%

% Define size of median filters
B1 = 11;
B2 = 3;

% Two part median filtering:
% fine scale
u_out = cellfun(@medfilt3, u_in, {B2,B2,B2}, {'symmetric','symmetric','symmetric'},'uni',false);
% coarse scale
u_out = cellfun(@medfilt3, u_out, {B1,B1,B1}, {'symmetric','symmetric','symmetric'},'uni',false);
end

author: thomaskuestner


text:
function sphere_integrals_test ( )

%*****************************************************************************80
%
%% SPHERE_INTEGRALS_TEST tests the SPHERE_INTEGRALS library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_INTEGRALS_TEST:\n' )
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SPHERE_INTEGRALS library.\n' );

  sphere_integrals_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_INTEGRALS_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function l = gaussian2PriorLogProb(prior, x)

% GAUSSIAN2PRIORLOGPROB Log probability of Gaussian prior.
% COPYRIGHT: Neil D. Lawrence 2004
% MODIFICATIONS: Andreas C. Damianou, 2013
% SHEFFIELDML

% Compute log prior
x_c = x - prior.mean;
l = length(prior.mean) * (-.5*sum(sum(prior.precision*x_c.*x_c + log(2*pi) - log(prior.precision))));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function rx = mrRx(vol,ref,varargin)
%
% rx = mrRx([vol],[referenceVol],[options]):
%
% A GUI to create, edit, and save 4 x 4 transformations to apply
% to volume matrices. Can be used to coregister a volume to a
% reference volume, evaluate motion and motion correction, and
% perform motion correction.
%
% vol: a 2D or 3D matrix (4D support may come later).
%
% referenceVol: Another volume, which is not transformed,
% to which the first volume may be compared.
%
% The reference volume need not be the same size as the
% transformed volume. But, it will be assumed they're at
% the same resolution, unless the 'volRes' or 'refRes' options
% are applied.
%
% Both vol and referenceVol can be supplied as strings specifying
% paths to volume files. Can load volumes in any of the formats
% supported by loadVolume (type help loadVolume for a list). This
% includes vAnatomy .dat, analyze .img, I-files, DICOM, and others.
%
% If invoked without input arguments, just opens the control
% window, from which transformed and reference volumes
% can be loaded with dialogs.
%
% Returns rx,a struct with info about the UI state, similar to a view.
% The key field in rx is rx.xform, which specifies a 4x4 affine
% transformation matrix to transform coordinates in the rx to coordinates
% in the volume. That is:
%   volCoords = rx.xform * rxCoords;
% though these transformations can be properly accessed with the functions
% rx2vol and vol2rx.
%
% Options:
% rxRes, [1 x 3 matrix]: specify the size of voxels in the prescription, in
%                 mm, as a 1x3 vector.
% volRes, [1 x 3 matrix]: specify the size of voxels in the prescribed volume,
%                  in mm, as a 1x3 vector.
% refRes, [1 x 3 matrix]: specify the size of voxels in the reference volume,
%                  in mm, as a 1x3 vector.
% rxDims, [1 x 3 matrix]: specify the dims (# of voxels -- rows cols slices)
%                  for the prescription.
% rxSizeMM, [1 x 3 matrix]: specify the total size of the prescription in mm.
%                   This is overridden if rxDims is specified.
%
% Note these parameters can be used to specify vol
% and ref as being at different resolutions, and the
% prescription at a third resolution.
%
% 02/17/05 ras.
if notDefined('vol')
    vol = [];
elseif ischar(vol)
    [vol volVoxelSize] = loadVolume(vol); %,'reorient');
	varargin = [varargin {'volRes' volVoxelSize}];
end

if notDefined('ref')
    ref = [];
elseif ischar(ref)
    ref = loadVolume(ref); %,'reorient');
end

% initialize the rx struct
rx = rxInit(vol, ref, varargin);

% open the main control figure:
rx = rxOpenControlFig(rx);

if ~isempty(vol)
    rx = rxOpenInterpFig(rx);
end

if ~isempty(ref)
    rx = rxOpenRefFig(rx);
end

% add rx as user data in control fig
set(rx.ui.controlFig,'UserData',rx);

rxRefresh(rx,1);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% =========================================================================
% FUNCTION
% scd_scheme_display.m
%
% INPUT
% gradient_vectors			nx3
% OR scheme             nx9
%
% OUTPUT
% (-)
%
% EXAMPLE
%   scd_bvecs_display(scheme)
%
% EXAMPLE 2
%   scd_bvecs_display(scd_schemefile_read('qspace.scheme'))
%
% COMMENTS
% Julien Cohen-Adad 2009-10-02
% =========================================================================
%
% See also scd_schemefile_read
function scd_scheme_display(gradient_vectors,display_voronoi)

%% Check if gradient_vectors is a schemefile
if isstr(gradient_vectors)
    error('Input should be a matrix not a filename. Use scd_schemefile_read')
end

if size(gradient_vectors,2)>3
    gradient_vectors = gradient_vectors(:,1:3).*repmat(gradient_vectors(:,4),[1 3]);
end

% gradient_vectors=gradient_vectors/max(max(gradient_vectors));
gradient_norm = sqrt(gradient_vectors(:,1).^2+gradient_vectors(:,2).^2+gradient_vectors(:,3).^2);
display_3d = 0;
colorindex = jet(1000);
colorprct = max(1,round(gradient_norm/max(gradient_norm)*1000));
color = colorindex(colorprct,:);
% display gradients
if display_3d
	for i = 1:size(gradient_vectors,1)
		plot3(gradient_vectors(i,1),gradient_vectors(i,2),gradient_vectors(i,3),'.','MarkerSize',10,'Color',color(i,:))
		hold on
	end
	xlabel('X')
	ylabel('Y')
	zlabel('Z')
	axis vis3d;
	view(3), axis equal
	axis on, grid
	rotate3d on;
end

% display gradients
subplot(2,2,1)
for i = 1:size(gradient_vectors,1)
    plot3(gradient_vectors(i,1),gradient_vectors(i,2),gradient_vectors(i,3),'.','MarkerSize',10,'Color',color(i,:))
    hold on
end
xlabel('X')
ylabel('Y')
zlabel('Z')
axis vis3d;
view(3), axis equal
axis on, grid
rotate3d on;
view(0,0)
lim = max(max(abs([xlim; ylim; zlim])));
xlim([-lim lim]); ylim([-lim lim]); zlim([-lim lim]);

subplot(2,2,2)
for i = 1:size(gradient_vectors,1)
    plot3(gradient_vectors(i,1),gradient_vectors(i,2),gradient_vectors(i,3),'.','MarkerSize',10,'Color',color(i,:))
    hold on
end
xlabel('X')
ylabel('Y')
zlabel('Z')
axis vis3d;
view(3), axis equal
axis on, grid
rotate3d on;
view(90,0)
lim = max(max(abs([xlim; ylim; zlim])));
xlim([-lim lim]); ylim([-lim lim]); zlim([-lim lim]);

subplot(2,2,3)
for i = 1:size(gradient_vectors,1)
    plot3(gradient_vectors(i,1),gradient_vectors(i,2),gradient_vectors(i,3),'.','MarkerSize',10,'Color',color(i,:))
    hold on
end
xlabel('X')
ylabel('Y')
zlabel('Z')
axis vis3d;
view(3), axis equal
axis on, grid
rotate3d on;
view(0,90)
lim = max(max(abs([xlim; ylim; zlim])));
xlim([-lim lim]); ylim([-lim lim]); zlim([-lim lim]);

% % Voronoi transformation
if exist('display_voronoi','var') && display_voronoi
	X = gradient_vectors;
	h_fig = figure('name','Voronoi');
	[V,C] = voronoin(X);
	K = convhulln(X);
	d = [1 2 3 1];       % Index into K
	for i = 1:size(K,1)
	   j = K(i,d);
	   h(i) = patch(X(j,1),X(j,2),X(j,3),i,'FaceColor','white','FaceLighting','phong','EdgeColor','black');
	end
	hold off
	view(2)
	axis off
	axis equal
	colormap(gray);
	% title('One cell of a Voronoi diagram')
	axis vis3d;
	rotate3d on;
	print(h_fig,'-dpng',strcat(['fig_voronoi.png']));
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function check = f_check ( m, n )

%*****************************************************************************80
%
%% F_CHECK checks the parameters of the F PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the parameters of the PDF.
%    0 < M
%    0 < N
%
%    Output, logical CHECK, is TRUE if the parameters are legal.
%
  if ( m ~= round ( m ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'F_CHECK - Fatal error!\n' );
    fprintf ( 1, '  M is not an integer.\n' );
    check = 0;
  end

  if ( m <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'F_CHECK - Fatal error!\n' );
    fprintf ( 1, '  M <= 0.\n' );
    check = 0;
  end

  if ( n ~= round ( n ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'F_CHECK - Fatal error!\n' );
    fprintf ( 1, '  N is not an integer.\n' );
    check = 0;
  end

  if ( n <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'F_CHECK - Fatal error!\n' );
    fprintf ( 1, '  N <= 0.\n' );
    check = 0;
  end

  check = 1;

  return
end


text:
function [tracks_interp,trk_mean_length] = trk_interp(tracks,nPoints_new,spacing,tie_at_center)
%TRK_INTERP - Interpolate tracks with cubic B-splines
%Streamlines will be resampled to have a new number of vertices. May be useful
%to groom your tracks first. Can be multithreaded if Parallel Computing Toolbox
%is installed.
%
% Syntax: [tracks_interp,trk_mean_length] = trk_interp(tracks,nPoints_new,spacing,tie_at_center)
%
% Inputs:
%    tracks      - Struc array output of TRK_READ [1 x nTracks]
%    nPoints_new - Constant # mode: Number of vertices for each streamline
%                  (spacing between vertices will vary between streamlines)
%                  (Default: 100)
%    spacing     - Constant spacing mode: Spacing between each vertex (# of
%                  vertices will vary between streamlines). Note: Only supply
%                  nPoints_new *OR* spacing!
%    tie_at_center - (Optional) Use with nPoints_new to add an additional
%                  "tie-down" point at the midpoint of the tract. Recommended.
%                  http://github.com/johncolby/along-tract-stats/wiki/tie-at-center
%
% Outputs:
%    tracks_interp   - Interpolated tracks in matrix format.
%                      [nPoints_new x 3 x nTracks]
%    trk_mean_length - The length of the mean tract geometry if using the
%                      spacing parameter, above. Useful to then normalize track
%                      lengths between subjects.
%
% Example:
%    exDir                   = '/path/to/along-tract-stats/example';
%    subDir                  = fullfile(exDir, 'subject1');
%    trkPath                 = fullfile(subDir, 'CST_L.trk');
%    volPath                 = fullfile(subDir, 'dti_fa.nii.gz');
%    volume                  = read_avw(volPath);
%    [header tracks]         = trk_read(trkPath);
%    tracks_interp           = trk_interp(tracks, 100);
%    tracks_interp           = trk_flip(header, tracks_interp, [97 110 4]);
%    tracks_interp_str       = trk_restruc(tracks_interp);
%    [header_sc tracks_sc]   = trk_add_sc(header, tracks_interp_str, volume, 'FA');
%    [scalar_mean scalar_sd] = trk_mean_sc(header_sc, tracks_sc);
%
% Other m-files required: Curve Fitting Toolbox (aka Spline Toolbox)
% Subfunctions: none
% MAT-files required: none
%
% See also: TRK_READ, SPLINE

% Author: John Colby (johncolby@ucla.edu)
% UCLA Developmental Cognitive Neuroimaging Group (Sowell Lab)
% Apr 2010

if nargin<4, tie_at_center = []; end
if nargin<3, spacing = []; end
if nargin<2 || isempty(nPoints_new), nPoints_new = 100; end

tracks_interp   = zeros(nPoints_new, 3, length(tracks));
trk_mean_length = [];
pp = repmat({[]},1,length(tracks));

% Interpolate streamlines so that each has the same number of vertices, spread
% evenly along its length (i.e. vertex spacing will vary between streamlines)
for iTrk=1:length(tracks)
    tracks_tmp = tracks(iTrk);

    % Martin: bugfix for single point fibers
    if tracks_tmp.nPoints == 1
        tracks_interp(:,:,iTrk) = repmat(tracks_tmp.matrix, nPoints_new, 1);
        continue
    end

    % Determine streamline segment lengths
    segs = sqrt(sum((tracks_tmp.matrix(2:end,1:3) - tracks_tmp.matrix(1:(end-1),1:3)).^2, 2));
    dist = [0; cumsum(segs)];

    % Remove duplicates
    [dist I J]= unique(dist);
    % Martin: bugfix for fibers with same points
    if length(dist) == 1
        dist = 0:1;
        I = 1:2;
    end

    % Fit spline
    % Martin: bugfix for fibers with more than 3 dimensions
    pp{iTrk} = spline(dist, tracks_tmp.matrix(I,1:3)');

    % Resample streamline along the spline
    tracks_interp(:,:,iTrk) = ppval(pp{iTrk}, linspace(0, max(dist), nPoints_new))';
end

% Interpolate streamlines so that the vertices have equal spacing for a central
% "tie-down" origin. This means streamlines will have varying #s of vertices
if ~isempty(spacing)
    % Calculate streamline lengths
    lengths = trk_length(tracks_interp);

    % Determine the mean tract geometry and grab the middle vertex
    track_mean      = mean(tracks_interp, 3);
    trk_mean_length = trk_length(track_mean);
    middle          = track_mean(round(length(track_mean)/2),:);

    % Interpolate streamlines again, but this time sample with constant vertex
    % spacing for all streamlines. This means that the longer streamlines will now
    % have more vertices.
    tracks_interp = repmat(struct('nPoints', 0, 'matrix', [], 'tiePoint', 0), 1, length(tracks));
    for iTrk=1:length(tracks)
        tracks_interp(iTrk).matrix  = ppval(pp{iTrk}, 0:spacing:lengths(iTrk))';
        tracks_interp(iTrk).nPoints = size(tracks_interp(iTrk).matrix, 1);

        % Also determine which vertex is the "tie down" point by finding the one
        % closest to the middle point of the mean tract geometry
        dists = sqrt(sum(bsxfun(@minus, tracks_interp(iTrk).matrix, middle).^2,2));
        [tmp, ind] = min(dists);
        tracks_interp(iTrk).tiePoint = ind;
    end
end

% Streamlines will all have the same # of vertices, but now they will be spread
% out so that an equal proportion lies on either side of a central origin.
if ~isempty(nPoints_new) && ~isempty(tie_at_center)
    % Make nPoints_new odd
    nPoints_new_odd = floor(nPoints_new/2)*2+1;

    % Calculate streamline lengths
    lengths = trk_length(tracks);

    % Determine the mean tract geometry and grab the middle vertex
    track_mean      = mean(tracks_interp, 3);
    trk_mean_length = trk_length(track_mean);
    middle          = track_mean(round(length(track_mean)/2),:);

    tracks_interp_tmp = zeros(nPoints_new_odd, 3, length(tracks));

    for iTrk=1:length(tracks)
        dists = sqrt(sum(bsxfun(@minus, tracks_interp(:,:,iTrk), middle).^2,2));
        [tmp, ind] = min(dists);

        first_half  = ppval(pp{iTrk}, linspace(0, lengths(iTrk)*(ind/nPoints_new), ceil(nPoints_new_odd/2)))';
        second_half = ppval(pp{iTrk}, linspace(lengths(iTrk)*(ind/nPoints_new), lengths(iTrk), ceil(nPoints_new_odd/2)))';
        tracks_interp_tmp(:,:,iTrk) = [first_half; second_half(2:end,:)];
    end

    tracks_interp = tracks_interp_tmp;
end

author: brainstorm-tools


text:
function [ ] = make_dir( folderPath )
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

folder = fileparts(folderPath);

if(~isdir(folder))
    mkdir(folder);
end

author: adikhosla

repo: feature-extraction