function price = PROJ_Parisian(N, call, down, S_0, W, H, M, r, rnCHF, T, Gamm, resetting, alph)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for Parisian-style barrier options using PROJ method
% Models Supported: Levy Processes, including jump diffusions and Black-Scholes model
% Returns: price of contract
% Author: Justin Lars Kirkby
%
% ----------------------
% Contract/Model Params
% ----------------------
% S_0 = initial stock price (e.g. 100)
% W   = strike  (e.g. 100)
% r   = interest rate (e.g. 0.05)
% T   = time remaining until maturity (in years, e.g. T=1)
% M   = number of subintervals of [0,T] (total of M+1 monitoring points in time grid, including S_0)
% call = 1 for call (else put)
% down = 1 for down and out (otherwise it's up and out)
% H    = barrier
% Gamm = maximum number of discretely monitored excursions into the knockout region allowed (more than this results in knockout)
% resetting = 1 if a reseting type parisian option (otherwise its cumulative, ie never resets)
% rnCHF = risk netural characteristic function (function handle with single argument)
%
% ----------------------
% Numerical (PROJ) Params
% ----------------------
% alph  = grid with is 2*alph
% N     = number of grid/basis points (power of 2, e.g. 2^12), resolution = 2*alph/(N-1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~(down == 1 && call ~= 1) && ~(down ~=1 && call == 1)
    fprintf('Sorry, currently only Up and out calls, and down and out puts have been coded \n')
    return
end

gamm0 = 1; %HARDCODED: this param would allow us to specify an inital consumed budget, in case of in-progress valuation

dt   = T/M;
nrdt = -r*dt;
h = log(H/S_0);
lws   = log(W/S_0);

%%%%%%  Gaussian Quad Constants
q_plus = (1 + sqrt(3/5))/2;  q_minus = (1 - sqrt(3/5))/2;
b3  = sqrt(15); b4 = b3/10;

dx = 2*alph/(N-1);
xmin = -alph/2;

n_h = floor((h-xmin)/dx +1);
xmin = h - (n_h -1)*dx;    %realign so that h is on the grid (this is important for the case where h = 0)

if h~= 0 %Realign so that h and 0 are both members of grid (if possible)
    nnot =  floor(1-xmin/dx);
    if abs(h) > dx  %so that n_h ~= nnot
        dx = (h - 0)/(n_h - nnot);
        xmin = dx*(1-nnot);  %hence nnot should remain on the grid
        %n_h = floor((h-xmin)/dx +1);  %NOT Numerically Stable
        n_h = floor(nnot + h/dx);  %Numerically Stable
    end
else
    nnot = n_h;
end

a    = 1/dx;
a2   = a^2;
zmin = (1 - N/2)*dx;  %Kbar corresponds to zero

%Cons = 24*a2/N;
Cons2 = 24*a2*exp(nrdt)/N;
dw    = 2*pi*a/N;
grand = dw*(1:N-1);
grand = exp(-1i*zmin*grand).*rnCHF(grand).*(sin(grand/(2*a))./grand).^2./(2+cos(grand/a));
beta  = Cons2*real(fft([1/(24*a2) grand]));   %%%%  NOTE: all toep matrices incorporate exp(-r*dt)

interp_Atend = 0;
if 0 < abs(h) && abs(h)<dx
    interp_Atend = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   DETERMINE COMMON Params
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K     = N/2;
nbar  = floor(a*(lws - xmin)+1);
rho   = lws - (xmin+(nbar - 1)*dx);
zeta  = a*rho;

toepM = [beta(K:-1:1)'; 0 ; beta(2*K-1:-1:K +1)'];   toepM = fft(toepM);

%%%% PAYOFF CONSTANTS-----------------------------------
varthet_01 = exp(.5*dx)*(5*cosh(b4*dx) - b3*sinh(b4*dx) + 4)/18;
varthet_m10 = exp(-.5*dx)*(5*cosh(b4*dx) + b3*sinh(b4*dx) + 4)/18;
varthet_star = varthet_01 + varthet_m10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if down == 1 && call ~= 1  %DOP
    %l = log(H/S_0);

    n_l = n_h;

    zeta_plus = zeta*q_plus; zeta_minus = zeta*q_minus;
    rho_plus = rho*q_plus; rho_minus = rho*q_minus;

    ed1 = exp(rho_minus); ed2 = exp(rho/2); ed3 = exp(rho_plus);

    dbar_1 = zeta^2/2;
    dbar_0 = zeta - dbar_1;         %  dbar_1 = zeta + .5*((zeta - 1)^2 - 1);
    d_0    = zeta*(5*( (1-zeta_minus)*ed1 + (1-zeta_plus)*ed3 ) + 4*(2-zeta)*ed2)/18;
    d_1    = zeta*( 5*(zeta_minus*ed1 + zeta_plus*ed3) + 4*zeta*ed2 )/18;

    %%%%Thet(1)        =  W/2 - H*varthet_01;
    Thet = zeros(K,1);
    Thet(1:nbar-1) =  W - exp(xmin +dx*(0:nbar-2))*S_0*varthet_star;
    Thet(nbar)     =  W*(.5 + dbar_0 - exp(-rho)*(varthet_m10 + d_0));
    Thet(nbar + 1) =  W*(dbar_1 - exp(- rho)*d_1);

    %%%%%%% Initialize Val
    Val = zeros(K,Gamm+1);
    p   = ifft(toepM.*fft([Thet;zeros(K,1)]));
    for j=1:Gamm
       Val(:,j) = p(1:K);
    end
    Thet(1:n_l-1) = 0;  %up to the barrier l, the coefficients are zero
    Thet(n_l) = W/2 - H*varthet_01;  %because this is the case where you can knock out at termination
    p   = ifft(toepM.*fft([Thet;zeros(K,1)]));
    Val(:,Gamm+1) = p(1:K);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% RESETTING PARISIAN
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if resetting == 1
        for m=M-2:-1:0
            %These Thet must be kept outside of loop, else they will be inadvertently altered
            Thet(n_l+1:K-1) = (Val(n_l:K-2,1)+10*Val(n_l+1:K-1,1)+Val(n_l+2:K,1))/12;
            Thet(K)         = (13*Val(K,1)+15*Val(K-1,1)-5*Val(K-2,1)+Val(K-3,1))/48;
            ThetPartial     = (13*Val(n_l,1)+15*Val(n_l+1,1)-5*Val(n_l+2,1)+Val(n_l+3,1))/48;
            for j=1:Gamm

                Thet(1)      = (13*Val(1,j+1)+15*Val(2,j+1)-5*Val(3,j+1)+Val(4,j+1))/48;
                Thet(2:n_l-1) = (Val(1:n_l-2,j+1)+10*Val(2:n_l-1,j+1)+Val(3:n_l,j+1))/12;
                Thet(n_l)     = (13*Val(n_l,j+1)+15*Val(n_l-1,j+1)-5*Val(n_l-2,j+1)+Val(n_l-3,j+1))/48 ...
                                 + ThetPartial;

                p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
                Val(:,j)  = p(1:K);
            end

            %%% Now to Gamm+1
            j = Gamm+1;
            Thet(1:n_l-1)   = 0;
            Thet(n_l)       = ThetPartial;

            p    = ifft(toepM.*fft([Thet; zeros(K,1)]));
            Val(:,j)  = p(1:K);
        end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% CUMULATIVE PARISIAN
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    else
        for m=M-2:-1:0

            for j=1:Gamm
                Thet(1)      = (13*Val(1,j+1)+15*Val(2,j+1)-5*Val(3,j+1)+Val(4,j+1))/48;
                Thet(2:n_l-1) = (Val(1:n_l-2,j+1)+10*Val(2:n_l-1,j+1)+Val(3:n_l,j+1))/12;

                Thet(n_l)     = (13*Val(n_l,j+1)+15*Val(n_l-1,j+1)-5*Val(n_l-2,j+1)+Val(n_l-3,j+1))/48 ...
                              + (13*Val(n_l,j)+15*Val(n_l+1,j)-5*Val(n_l+2,j)+Val(n_l+3,j))/48;

                Thet(n_l+1:K-1) = (Val(n_l:K-2,j)+10*Val(n_l+1:K-1,j)+Val(n_l+2:K,j))/12;
                Thet(K)        = (13*Val(K,j)+15*Val(K-1,j)-5*Val(K-2,j)+Val(K-3,j))/48;

                p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
                Val(:,j)  = p(1:K);
            end

            %%% Now to Gamm+1
            j = Gamm+1;
            Thet(1:n_l-1)   = 0;
            Thet(n_l)       = (13*Val(n_l,j)+15*Val(n_l+1,j)-5*Val(n_l+2,j)+Val(n_l+3,j))/48;
            Thet(n_l+1:K-1) = (Val(n_l:K-2,j)+10*Val(n_l+1:K-1,j)+Val(n_l+2:K,j))/12;

            p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
            Val(:,j)  = p(1:K);
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif down ~=1 && call == 1  %UOC
    %u    = log(H/S_0);
    n_u = n_h;

    sigma = 1 - zeta; sigma_plus = (q_plus-.5)*sigma; sigma_minus = (q_minus-.5)*sigma;

    es1 = exp(dx*sigma_plus); es2 = exp(dx*sigma_minus);

    dbar_0 = .5 + zeta*(.5*zeta-1);
    dbar_1 = sigma*(1 - .5*sigma);

    d_0 = exp((rho+dx)*.5)*sigma^2/18*(5*((1-q_minus)*es2 +(1-q_plus)*es1) + 4);
    d_1 = exp((rho-dx)*.5)*sigma/18*(5*( (.5*(zeta+1) +sigma_minus)*es2 + (.5*(zeta+1) +sigma_plus)*es1 ) + 4*(zeta+1) );

    %%%%Thet(1)        =  W/2 - H*varthet_01;
    Thet = zeros(K,1);
    Thet(nbar)     =  W*(exp(-rho)*d_0 - dbar_0);
    Thet(nbar + 1) =  W*(exp(dx-rho)*(varthet_01 +d_1) -(.5 + dbar_1) );
    Thet(nbar + 2:K)  = exp(xmin +dx*(nbar+1:K-1))*S_0*varthet_star - W;

    Val = zeros(K,Gamm+1);
    p   = ifft(toepM.*fft([Thet;zeros(K,1)]));
    for j=1:Gamm
       Val(:,j) = p(1:K);
    end

    Thet(n_u+1:K) = 0;  %after the barrier u, the coefficients are zero
    Thet(n_u) = H*varthet_m10 - .5*W;  %because this is the case where you can knock out at termination
    p   = ifft(toepM.*fft([Thet;zeros(K,1)]));
    Val(:,Gamm+1) = p(1:K);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% RESETTING PARISIAN
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% For the resetting, its more efficient to pull the
    %%% Thet(1:n_u-1) ouside of the loop through j = 1:Gamm
    if resetting == 1
        for m=M-2:-1:0
            %NOTE: these must be defined outside of loop, else we change them and then the next function requires the unchanged value!

            Thet(1)      = (13*Val(1,1)+15*Val(2,1)-5*Val(3,1)+Val(4,1))/48;
            Thet(2:n_u-1) = (Val(1:n_u-2,1)+10*Val(2:n_u-1,1)+Val(3:n_u,1))/12;
            ThetPartial = (13*Val(n_u,1)+15*Val(n_u-1,1)-5*Val(n_u-2,1)+Val(n_u-3,1))/48 ;

            for j=1:Gamm
                Thet(n_u)     = ThetPartial ...
                    + (13*Val(n_u,j+1)+15*Val(n_u+1,j+1)-5*Val(n_u+2,j+1)+Val(n_u+3,j+1))/48;
                Thet(n_u+1:K-1) = (Val(n_u:K-2,j+1)+10*Val(n_u+1:K-1,j+1)+Val(n_u+2:K,j+1))/12;
                Thet(K)        = (13*Val(K,j+1)+15*Val(K-1,j+1)-5*Val(K-2,j+1)+Val(K-3,j+1))/48;

                p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
                Val(:,j)  = p(1:K);
            end

            %%% Now to Gamm+1
            j = Gamm+1;
            Thet(n_u)     = ThetPartial;
            Thet(n_u+1:K) = 0;

            p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
            Val(:,j)  = p(1:K);
        end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% CUMULATIVE PARISIAN
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    else
        for m=M-2:-1:0
            for j=1:Gamm
                Thet(1)      = (13*Val(1,j)+15*Val(2,j)-5*Val(3,j)+Val(4,j))/48;
                Thet(2:n_u-1) = (Val(1:n_u-2,j)+10*Val(2:n_u-1,j)+Val(3:n_u,j))/12;

                Thet(n_u)     = (13*Val(n_u,j)+15*Val(n_u-1,j)-5*Val(n_u-2,j)+Val(n_u-3,j))/48 ...
                              + (13*Val(n_u,j+1)+15*Val(n_u+1,j+1)-5*Val(n_u+2,j+1)+Val(n_u+3,j+1))/48;

                Thet(n_u+1:K-1) = (Val(n_u:K-2,j+1)+10*Val(n_u+1:K-1,j+1)+Val(n_u+2:K,j+1))/12;
                Thet(K)        = (13*Val(K,j+1)+15*Val(K-1,j+1)-5*Val(K-2,j+1)+Val(K-3,j+1))/48;

                p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
                Val(:,j)  = p(1:K);
            end

            %%% Now to Gamm+1
            j = Gamm+1;
            Thet(1)       = (13*Val(1,j)+15*Val(2,j)-5*Val(3,j)+Val(4,j))/48;
            Thet(2:n_u-1) = (Val(1:n_u-2,j)+10*Val(2:n_u-1,j)+Val(3:n_u,j))/12;
            Thet(n_u)     = (13*Val(n_u,j)+15*Val(n_u-1,j)-5*Val(n_u-2,j)+Val(n_u-3,j))/48;
            Thet(n_u+1:K)  =0;

            p    = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));
            Val(:,j)  = p(1:K);
        end
    end
end

if interp_Atend ~= 1
    price = Val(nnot,gamm0);

else  %%% INTERPOLATION
%     Use 5 Point Cubic Interpolation
%     xnot = xmin +(nnot-1)*dx;
%     xs = [xnot-2*dx,xnot-dx,xnot,xnot+dx,xnot+2*dx];
%     ys = [Val(nnot-2,1),Val(nnot-1,1),Val(nnot,1),Val(nnot+1,1),Val(nnot+2,1)];
%     price = spline(xs,ys,0);

    dd = 0 - (xmin+ (nnot -1)*dx);
    price = Val(nnot,gamm0) + (Val(nnot+1,gamm0) - Val(nnot,gamm0))*dd/dx;
end

end


text:
function [rotDetector,angles]=coord2angles(geo,det_origin,up_vec,norm_vec)
%COORD2ANGLES converts real world coordinates of the detector to TIGRE angles
%
%
% The detector origin is assumed to have no detector offset!
%
% There are 2 relevant rotations here. The ones around the detector center
% and the ones around the object center. From the origin and 2 vectors we
% can obtain both.

% Input checking:

assert(((geo.DSD-geo.DSO)-sqrt(sum(det_origin(:).^2))<1e-4),'Detector is not at geo.DOD distance from origin.');

assert(dot(up_vec,norm_vec)<1e-4,'Up vector and normal vector are not orthonormal');

% Lets normalize the vectors
up_vec=up_vec/sqrt(sum(up_vec(:).^2));
norm_vec=norm_vec/sqrt(sum(norm_vec(:).^2));

%% first lets obtain local rotation (around detector origin)
% Lets define the vectors we want

non_rotated_norm_vec=-det_origin/sqrt(sum(det_origin(:).^2));

% The up vector is harder, as there are infinte posibilities. We chose the
% one that is angularly closer to the real one. This is done by projecting
% the vector onto the plane
% https://math.stackexchange.com/questions/633181/formula-to-project-a-vector-onto-a-plane

non_rotated_up_vec=up_vec-(dot(up_vec,non_rotated_norm_vec)*non_rotated_norm_vec);

% solve! (Wahba's problem using SVD)

B=up_vec*non_rotated_up_vec.'+norm_vec*non_rotated_norm_vec.'; % this should be 3x3
assert(isequal(size(B),[3 3]),'input vectors are not 3x1, (they are likely 1x3)')

[U,~,V]=svd(B);
M=eye(3); M(3,3)=det(U)*det(V);
R=U*M*V';

% This rotation is roll-pitch-yaw.
% Computing Euler angles from a rotation matrix -Greg Slabaugh

if (abs(R(3,1))-1) < 1e-10
    pitch1=-asin(R(3,1));
    % pitch2=pi-pitch1;
    yaw1=atan2(R(3,2)/cos(pitch1),R(3,3)/cos(pitch1));
    % yaw2=atan2(R(3,2)/cos(pitch2),R(3,3)/cos(pitch2));
    roll1=atan2(R(2,1)/cos(pitch1),R(1,1)/cos(pitch1));
    % roll2=atan2(R(2,1)/cos(pitch2),R(1,1)/cos(pitch2));
else
    roll1=0;
    if sign(R(3,1))==1
        pitch1=-pi/2;
        yaw1=-roll1+atan2(-R(1,2),-R(1,3));
    else
        pitch1=pi/2;
        yaw1=roll1+atan2(R(1,2),R(1,3));
    end
end

rotDetector=[roll1;pitch1;yaw1];

%% Now lets get global rotation around the origin. We need new vectors

origin_norm_vec=det_origin+non_rotated_norm_vec;
origin_up_vec=det_origin+non_rotated_up_vec;

% now we can reuse the variables
non_rotated_norm_vec=[geo.DSO-1;0;0];
non_rotated_up_vec=[geo.DSO;0;1];

B=origin_up_vec*non_rotated_up_vec.'+origin_norm_vec*non_rotated_norm_vec.'; % this should be 3x3
[U,~,V]=svd(B);
M=eye(3); M(3,3)=det(U)*det(V);
R=U*M*V';

if (abs(R(3,3))-1) < 1e-10

    angle21=acos(R(3,3));
    angle22=angle21+pi;
    angle31=atan2(R(3,2)/sin(angle21),-R(3,1)/sin(angle21));
    angle32=atan2(R(3,2)/sin(angle22),-R(3,1)/sin(angle22));
    angle11=atan2(R(2,3)/sin(angle21),R(1,3)/sin(angle21));
    angle12=atan2(R(2,3)/sin(angle22),R(1,3)/sin(angle22));
else
    angle31=0;
    if sign(R(3,3))==1
        angle21=pi/2;
        angle11=-angle31+atan2(R(1,1),R(2,1));
    else
        angle21=-pi/2;
        angle21=pi/2;
        angle11=angle31+atan2(-R(1,1),R(2,1));
    end
end
angles=[angle11;angle21;angle31];

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function Pe = bit_channel_upgrading_procedure(W, z, miu)
N = length(z);
if N == 1
    Pe = 0.5 * sum(min(W));
    disp(['Bit index = ' num2str(z) '  ML detection Bit Error rate = ' num2str(Pe)])
else
    W_up = get_W_up(W);
    W_up = LR_sort(W_up);
    W_up_after_erasure_symbol_merge = erasure_symbol_merge(W_up);
    W_up_after_merge = upgrading_merge(W_up_after_erasure_symbol_merge, miu);
    Pe1 = bit_channel_upgrading_procedure(W_up_after_merge, z(1 : N/2), miu);

    W_down= get_W_down(W);
    W_down = LR_sort(W_down);
    W_down_after_erasure_symbol_merge = erasure_symbol_merge(W_down);
    W_down_after_merge = upgrading_merge(W_down_after_erasure_symbol_merge, miu);
    Pe2 = bit_channel_upgrading_procedure(W_down_after_merge, z(N/2 + 1 : end), miu);

    Pe = [Pe1 Pe2];
end
end

% function Pe = bit_channel_upgrading_procedure(W, z, miu)
% N = length(z);
% m = round(log2(N));
% Pe = zeros(N, 1);
% for k = N - 7 : N - 1
%     char_bin_expansion = dec2bin(k, m);
%     W_tmp = W;
%     for i_level = 1 : m
%         if char_bin_expansion(i_level) == '0'
%             W_up = get_W_up(W_tmp);
%             W_up = LR_sort(W_up);
%             W_up_after_erasure_symbol_merge = erasure_symbol_merge(W_up);
%             W_tmp = upgrading_merge(W_up_after_erasure_symbol_merge, miu);
%         else
%             W_down = get_W_down(W_tmp);
%             W_down = LR_sort(W_down);
%             W_down_after_erasure_symbol_merge = erasure_symbol_merge(W_down);
%             W_tmp = upgrading_merge(W_down_after_erasure_symbol_merge, miu);
%         end
%     end
%     Pe(k + 1) = 0.5 * sum(min(W_tmp));
%     disp(['Bit index = ' num2str(k + 1) ' ML detection Bit Error rate = ' num2str(Pe(k + 1))])
% end
% end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
% STE_WLP Short-time-energy-weighted linear prediction
%   [A,w] = ste_wlp(s,p,m,k)
%
% Description
%   This function fits linear prediction coefficients to the analysis
%   frame using the basic form of weighted linear prediction (WLP), i.e.,
%   by weighting each value of the squared prediction error by the short-time
%   energy (STE) of the previous samples.
%
% Inputs
%   s      : Speech signal frame [samples]
%   p      : Order of WLP analysis
%   m      : Length of the STE window (default to p)
%   k      : Delay of the STE window (default to 1)
%
% Outputs
%   A      : Linear prediction inverse filter coefficients
%   w      : the STE weighting function
%
% Notes
%   Even though coefficients are solved using a criterion similar to the
%   autocorrelation method of linear prediction, because of the weighting
%   the filter is not guaranteed to be stable. If a stable synthesis filter
%   is required, use STE_SWLP instead.
%
% Example
%   A = ste_wlp(s,p) gives the linear predictive inverse filter
%   coefficients optimized using STE-WLP
%
% References
%  [1] C. Ma, Y. Kamp and L. F. Willems, "Robust signal selection for
%  linear prediction analysis of voiced speech", Speech Communication, vol.
%  12, no. 1, pp. 6981, 1993.
%
% Copyright (c) 2013 Aalto University
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Common Speech Processing Repository
% http://covarep.github.io/covarep/
%
% Octave Compatible
%
% Author
%  Jouni Pohjalainen jouni.pohjalainen@aalto.fi
%
% $Id <info set by the versioning system> $

function [A,w] = ste_wlp(s,p,m,k)

% handle the special case of all-zero frames
if all(s==0)
  s = eps*randn(size(s));
end

N = length(s);

% default value for the STE window length
if nargin<3
    m = p;
end

% default value for the STE window lag
if nargin<4
    k = 1;
end

% compute STE weighting function
w = stew(s,m,p,k)+eps;

% apply the square root of the weighting function to the delayed versions
% of the signal from lag 0 to lag p
wsr = sqrt(w(1:(N+p)));
Y = zeros(N+p,p+1);
for i1=0:p
    Y(:,i1+1) = [zeros(i1,1);s;zeros(p-i1,1)].*wsr;
end

% compute weighted autocorrelations
R = (Y'*Y)/N;

% solve the p predictor coefficients
A = R(2:end,2:end)\R(2:end,1);

% convert to inverse filter form
A = [1;-A]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function w = stew(x,m,p,k)
% STE weighting for an order p predictor using window of length m delayed
% by k samples
w = conv([zeros(k,1);x(:)].^2,ones(m,1));
w = [w;zeros(p-m,1)];

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
% NSA v1 (Aybat et al. 2011)
% process_video('RPCA', 'NSA1', 'dataset/demo.avi', 'output/demo_NSA1.avi');
stdev = 1;
tol = 5e-6; % optimality tolerance for stopping_type 1
L = nsa_v1(M,stdev,tol,1);
S = M - L;

author: andrewssobral

repo: lrslibrary


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)

   if ~exist('hdr','var') || ~exist('filetype','var') || ~exist('fileprefix','var') || ~exist('machine','var')
      error('Usage: [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,[img_idx],[dim5_idx],[dim6_idx],[dim7_idx],[old_RGB]);');
   end

   if ~exist('img_idx','var') || isempty(img_idx) || hdr.dime.dim(5)<1
      img_idx = [];
   end

   if ~exist('dim5_idx','var') || isempty(dim5_idx) || hdr.dime.dim(6)<1
      dim5_idx = [];
   end

   if ~exist('dim6_idx','var') || isempty(dim6_idx) || hdr.dime.dim(7)<1
      dim6_idx = [];
   end

   if ~exist('dim7_idx','var') || isempty(dim7_idx) || hdr.dime.dim(8)<1
      dim7_idx = [];
   end

   if ~exist('old_RGB','var') || isempty(old_RGB)
      old_RGB = 0;
   end

   %  check img_idx
   %
   if ~isempty(img_idx) && ~isnumeric(img_idx)
      error('"img_idx" should be a numerical array.');
   end

   if length(unique(img_idx)) ~= length(img_idx)
      error('Duplicate image index in "img_idx"');
   end

   if ~isempty(img_idx) && (min(img_idx) < 1 || max(img_idx) > hdr.dime.dim(5))
      max_range = hdr.dime.dim(5);

      if max_range == 1
         error(['"img_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"img_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim5_idx
   %
   if ~isempty(dim5_idx) && ~isnumeric(dim5_idx)
      error('"dim5_idx" should be a numerical array.');
   end

   if length(unique(dim5_idx)) ~= length(dim5_idx)
      error('Duplicate index in "dim5_idx"');
   end

   if ~isempty(dim5_idx) && (min(dim5_idx) < 1 || max(dim5_idx) > hdr.dime.dim(6))
      max_range = hdr.dime.dim(6);

      if max_range == 1
         error(['"dim5_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim5_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim6_idx
   %
   if ~isempty(dim6_idx) && ~isnumeric(dim6_idx)
      error('"dim6_idx" should be a numerical array.');
   end

   if length(unique(dim6_idx)) ~= length(dim6_idx)
      error('Duplicate index in "dim6_idx"');
   end

   if ~isempty(dim6_idx) && (min(dim6_idx) < 1 || max(dim6_idx) > hdr.dime.dim(7))
      max_range = hdr.dime.dim(7);

      if max_range == 1
         error(['"dim6_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim6_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim7_idx
   %
   if ~isempty(dim7_idx) && ~isnumeric(dim7_idx)
      error('"dim7_idx" should be a numerical array.');
   end

   if length(unique(dim7_idx)) ~= length(dim7_idx)
      error('Duplicate index in "dim7_idx"');
   end

   if ~isempty(dim7_idx) && (min(dim7_idx) < 1 || max(dim7_idx) > hdr.dime.dim(8))
      max_range = hdr.dime.dim(8);

      if max_range == 1
         error(['"dim7_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim7_idx" should be an integer within the range of [' range '].']);
      end
   end

   [img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB);

   return					% load_nii_img

%---------------------------------------------------------------------
function [img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)

   switch filetype
   case {0, 1}
      fn = [fileprefix '.img'];
   case 2
      fn = [fileprefix '.nii'];
   end

   fid = fopen(fn,'r',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   end

   %  Set bitpix according to datatype
   %
   %  /*Acceptable values for datatype are*/
   %
   %     0 None                     (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN
   %     1 Binary                         (ubit1, bitpix=1) % DT_BINARY
   %     2 Unsigned char         (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8
   %     4 Signed short                  (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16
   %     8 Signed integer                (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32
   %    16 Floating point    (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32
   %    32 Complex, 2 float32      (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
   %    64 Double precision  (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64
   %   128 uint8 RGB                 (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24
   %   256 Signed char            (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8
   %   511 Single RGB              (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
   %   512 Unsigned short               (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16
   %   768 Unsigned integer             (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32
   %  1024 Signed long long              (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
   %  1280 Unsigned long long           (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64
   %  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
   %  1792 Complex128, 2 float64  (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
   %  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
   %
   switch hdr.dime.datatype
   case   1,
      hdr.dime.bitpix = 1;  precision = 'ubit1';
   case   2,
      hdr.dime.bitpix = 8;  precision = 'uint8';
   case   4,
      hdr.dime.bitpix = 16; precision = 'int16';
   case   8,
      hdr.dime.bitpix = 32; precision = 'int32';
   case  16,
      hdr.dime.bitpix = 32; precision = 'float32';
   case  32,
      hdr.dime.bitpix = 64; precision = 'float32';
   case  64,
      hdr.dime.bitpix = 64; precision = 'float64';
   case 128,
      hdr.dime.bitpix = 24; precision = 'uint8';
   case 256
      hdr.dime.bitpix = 8;  precision = 'int8';
   case 511
      hdr.dime.bitpix = 96; precision = 'float32';
   case 512
      hdr.dime.bitpix = 16; precision = 'uint16';
   case 768
      hdr.dime.bitpix = 32; precision = 'uint32';
   case 1024
      hdr.dime.bitpix = 64; precision = 'int64';
   case 1280
      hdr.dime.bitpix = 64; precision = 'uint64';
   case 1792,
      hdr.dime.bitpix = 128; precision = 'float64';
   otherwise
      error('This datatype is not supported');
   end

   hdr.dime.dim(find(hdr.dime.dim < 1)) = 1;

   %  move pointer to the start of image block
   %
   switch filetype
   case {0, 1}
      fseek(fid, 0, 'bof');
   case 2
      fseek(fid, hdr.dime.vox_offset, 'bof');
   end

   %  Load whole image block for old Analyze format or binary image;
   %  otherwise, load images that are specified in img_idx, dim5_idx,
   %  dim6_idx, and dim7_idx
   %
   %  For binary image, we have to read all because pos can not be
   %  seeked in bit and can not be calculated the way below.
   %
   if hdr.dime.datatype == 1 || isequal(hdr.dime.dim(5:8),ones(1,4)) || ...
	(isempty(img_idx) && isempty(dim5_idx) && isempty(dim6_idx) && isempty(dim7_idx))

      %  For each frame, precision of value will be read
      %  in img_siz times, where img_siz is only the
      %  dimension size of an image, not the byte storage
      %  size of an image.
      %
      img_siz = prod(hdr.dime.dim(2:8));

      %  For complex float32 or complex float64, voxel values
      %  include [real, imag]
      %
      if hdr.dime.datatype == 32 || hdr.dime.datatype == 1792
         img_siz = img_siz * 2;
      end

      %MPH: For RGB24, voxel values include 3 separate color planes
      %
      if hdr.dime.datatype == 128 || hdr.dime.datatype == 511
	 img_siz = img_siz * 3;
      end

      img = fread(fid, img_siz, sprintf('*%s',precision));

      d1 = hdr.dime.dim(2);
      d2 = hdr.dime.dim(3);
      d3 = hdr.dime.dim(4);
      d4 = hdr.dime.dim(5);
      d5 = hdr.dime.dim(6);
      d6 = hdr.dime.dim(7);
      d7 = hdr.dime.dim(8);

      if isempty(img_idx)
         img_idx = 1:d4;
      end

      if isempty(dim5_idx)
         dim5_idx = 1:d5;
      end

      if isempty(dim6_idx)
         dim6_idx = 1:d6;
      end

      if isempty(dim7_idx)
         dim7_idx = 1:d7;
      end
   else

      d1 = hdr.dime.dim(2);
      d2 = hdr.dime.dim(3);
      d3 = hdr.dime.dim(4);
      d4 = hdr.dime.dim(5);
      d5 = hdr.dime.dim(6);
      d6 = hdr.dime.dim(7);
      d7 = hdr.dime.dim(8);

      if isempty(img_idx)
         img_idx = 1:d4;
      end

      if isempty(dim5_idx)
         dim5_idx = 1:d5;
      end

      if isempty(dim6_idx)
         dim6_idx = 1:d6;
      end

      if isempty(dim7_idx)
         dim7_idx = 1:d7;
      end

      %  compute size of one image
      %
      img_siz = prod(hdr.dime.dim(2:4));

      %  For complex float32 or complex float64, voxel values
      %  include [real, imag]
      %
      if hdr.dime.datatype == 32 || hdr.dime.datatype == 1792
         img_siz = img_siz * 2;
      end

      %MPH: For RGB24, voxel values include 3 separate color planes
      %
      if hdr.dime.datatype == 128 || hdr.dime.datatype == 511
         img_siz = img_siz * 3;
      end

      % preallocate img
      img = zeros(img_siz, length(img_idx)*length(dim5_idx)*length(dim6_idx)*length(dim7_idx) );
      currentIndex = 1;

      for i7=1:length(dim7_idx)
         for i6=1:length(dim6_idx)
            for i5=1:length(dim5_idx)
               for t=1:length(img_idx)

                  %  Position is seeked in bytes. To convert dimension size
                  %  to byte storage size, hdr.dime.bitpix/8 will be
                  %  applied.
                  %
                  pos = sub2ind([d1 d2 d3 d4 d5 d6 d7], 1, 1, 1, ...
			img_idx(t), dim5_idx(i5),dim6_idx(i6),dim7_idx(i7)) -1;
                  pos = pos * hdr.dime.bitpix/8;

                  if filetype == 2
                     fseek(fid, pos + hdr.dime.vox_offset, 'bof');
                  else
                     fseek(fid, pos, 'bof');
                  end

                  %  For each frame, fread will read precision of value
                  %  in img_siz times
                  %
                  img(:,currentIndex) = fread(fid, img_siz, sprintf('*%s',precision));
                  currentIndex = currentIndex +1;

               end
            end
         end
      end
   end

   %  For complex float32 or complex float64, voxel values
   %  include [real, imag]
   %
   if hdr.dime.datatype == 32 || hdr.dime.datatype == 1792
      img = reshape(img, [2, length(img)/2]);
      img = complex(img(1,:)', img(2,:)');
   end

   fclose(fid);

   %  Update the global min and max values
   %
   hdr.dime.glmax = double(max(img(:)));
   hdr.dime.glmin = double(min(img(:)));

   %  old_RGB treat RGB slice by slice, now it is treated voxel by voxel
   %
   if old_RGB && hdr.dime.datatype == 128 && hdr.dime.bitpix == 24
      % remove squeeze
      img = (reshape(img, [hdr.dime.dim(2:3) 3 hdr.dime.dim(4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
      img = permute(img, [1 2 4 3 5 6 7 8]);
   elseif hdr.dime.datatype == 128 && hdr.dime.bitpix == 24
      % remove squeeze
      img = (reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
      img = permute(img, [2 3 4 1 5 6 7 8]);
   elseif hdr.dime.datatype == 511 && hdr.dime.bitpix == 96
      img = double(img(:));
      img = single((img - min(img))/(max(img) - min(img)));
      % remove squeeze
      img = (reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
      img = permute(img, [2 3 4 1 5 6 7 8]);
   else
      % remove squeeze
      img = (reshape(img, [hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
   end

   if ~isempty(img_idx)
      hdr.dime.dim(5) = length(img_idx);
   end

   if ~isempty(dim5_idx)
      hdr.dime.dim(6) = length(dim5_idx);
   end

   if ~isempty(dim6_idx)
      hdr.dime.dim(7) = length(dim6_idx);
   end

   if ~isempty(dim7_idx)
      hdr.dime.dim(8) = length(dim7_idx);
   end

   return						% read_image

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function initials = tc_condInitials(condNames);
%
% initials = tc_condInitials(condNames);
%
% Provides brief initials for a set of condition names.
%
% PROBLEM: when creating bar graphs or other images where conditions
% are plotted against one another, oftentimes there isn't space for the
% full condition label. (And often there's a legend with the full label
% next to the colors). So, just as a reminder, we want to have some
% abbreviated form of the condition names.
%
% SOLUTION: This code figures out 'good' initials, meaning initials that
% are the smallest possible to be unambiguous. If each condition name starts
% with a different letter, it will just be the first letters. If some are
% the same, it will use heuristics like finding capital letters or spacers
% (e.g. '_'), and making each initial unique.
%
% The returned initials are all capital, in a cell array of strings.
%
% ras, 02/06.
nConds = length(condNames);

% check that there aren't duplicate condition names: if so, the
% code below may get stuck in a loop (because there are always non-unique
% identifiers), so fix them by appending numbers
condNames = lower(condNames);
for i = 1:nConds
    test = cellfind(condNames, condNames{i});
    if length(test) > 1
        msg = sprintf(['Condition name %s is not unique! ' ...
                       'Appears in entries %s. ' ...
                       'Am numbering these duplicate entries. '], ...
                       condNames{i}, num2str(test));
        warning(msg);
        for j = test
            condNames{j} = [condNames{j} num2str(j)];
        end
    end
end

% remove any leading and trailing spaces from condNames:
for i = 1:nConds
    while isequal(condNames{i}(1), ' ')
        condNames{i} = condNames{i}(2:end);
    end
end

% initialize the initials ( :-) ) as the first letter of each cond name
for i = 1:nConds
    initials{i} = upper(condNames{i}(1));
end

% if each letter is unique, great! We're done.
if length(unique(initials))==length(initials)
    return;
end

% if we got here, there are some non-unique letters. Find them and fix
% them:

% First, go through each non-unique cond name, and see if we can find
% a 'second word', meaning a string after a space or an underscore,
% or a capital letter (or number). If we find one for a condition name,
% reshape that condition name to be 2D (words by characters in word).
for i = 1:nConds
    str = condNames{i};

    % blank spacers
    spacers = [findstr('_', str) findstr(' ', str)];

    % upper-case letters (ASCII range 65-90): find the preceding letter as
    % the spacer:
    spacers = [spacers find(ismember(str, 65:90))-1];

    % do the same operation for numbers as for upper-case letters, only
	% the ASCII range for number strings is 48-57:
    spacers = [spacers find(ismember(str, 48:57))-1];

    % sort into ascending indices in the string; remove spacers at end
    % points:
    spacers = spacers(spacers>1 & spacers<length(str));
    spacers = unique( sort(spacers) );
    if ~isempty(spacers)
        for j = spacers
            str = strvcat(str(1:j-1), str(j+1:end));
        end
        condNames{i} = str;
    end
end

notUnique = ones(size(initials)); % initialize index of remaining non-unique entries
while any(notUnique)
    % find which entries are repeated
    for i = 1:length(initials)
        if length(cellfind(initials, initials{i})) > 1
            notUnique(i) = 1;
        else
            notUnique(i) = 0;
        end
    end

    % for each non-unique entry, add one more character
    for j = find(notUnique)
        N = min( length(initials{j})+1, prod(size(condNames{j})) );
        initials{j} = upper(condNames{j}(1:N));

		% if the initials are as long as the condition name, stop.
		if N == numel(condNames{j})
			notUnique(j) = 0;
		end
    end
end

initials = upper(initials);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
fprintf('All results should be zero\n');
pochhammer(5.6,0)
pochhammer(1e4,3) - 27.63132109093154778328
pochhammer(1e4,1e3) - 9258.70469397627827724045
pochhammer(1e5,3) - 34.538806394660688260
pochhammer(1e6,3) - 41.44653467389032232
pochhammer(0.1,100) - 357.3415591371226056543663

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function obj = setAxis(obj, axis)
    % set the joint rotation axis vector in the joint coordinate
    %
    % Parameters:
    % axis: the axis vector @type rowvec

    if iscolumn(axis)
        axis = axis';
    end
    % validate if it is a numeric 1x3 vector
    validateattributes(axis, {'numeric'},{'size',[1,3]});

    tol = 1e-1;
    assert(abs(norm(axis) - 1) < tol,...
        'The axis must be a unit vector.');
    obj.Axis = axis;

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function params = gpSubspaceExtractParam(model)

% GPSUBSPACEEXTRACTPARAM
%
% COPYRIGHT : Carl Henrik Ek, 2008

% GP

params = gpExtractParam(model);

return;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function data = reshapeData(~, data, nVars)
%RESHAPEDATA   Extract half of the values since the values have been doubled-up
%with the DFS method.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Note: This is used at the end of the time-stepping in SPINOPERATOR/SOLVEPDE
% when constructing the SHPEREFUN output from the discrete values.

% Number of grid points in each direction:
N = length(data{end})/nVars;

% Loop over the entries of the CELL-ARRAY DATA:
for k = 1:length(data)
    % Extract half of the values and the real part since complex-valued
    % SPHEREFUN are not supported; do it for each variable:
    vals = [];
    for i = 1:nVars
        idx = (i-1)*N + 1;
        temp = data{k}(idx:idx+N-1, :);
        vals = [vals; real(temp([floor(N/2)+1:N 1], :))]; %#ok<*AGROW>
    end
    data{k} = vals;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function  varargout = drawPlatform(plane, siz, varargin)
%DRAWPLATFORM Draw a rectangular platform with a given size.
%
%   drawPlatform(PLANE, SIZ) draws a rectangular platform with the
%   dimensions specified by SIZ. If SIZ contains only one value instead of
%   two the platform will be quadratic.
%
%   drawPlatform(...,'PropertyName',PropertyValue,...) sets the value of
%   the specified patch property. Multiple property values can be set with
%   a single statement. See function patch for details.
%
%   drawPlane3d(AX,...) plots into AX instead of GCA.
%
%   H = drawPlatform(...) returns a handle H to the patch object.
%
%   Example
%
%     p0 = [1 2 3];
%     v1 = [1 0 1];
%     v2 = [0 -1 1];
%     plane = [p0 v1 v2];
%     axis([-10 10 -10 10 -10 10]);
%     drawPlatform(plane, [7,3])
%     set(gcf, 'renderer', 'zbuffer');
%
%   See also
%   planes3d, createPlane, patch

% ------
% Author: oqilipo
% Created: 2018-08-09
% Copyright 2018

%% Parse inputs

% extract axis handle
if numel(plane) == 1 && ishandle(plane)
    hAx = plane;
    plane = siz;
    siz = varargin{1};
    varargin(1) = [];
else
    hAx = gca;
end

% parse optional arguments
p = inputParser;
addRequired(p, 'plane', @(x) size(x,1)==1 && isPlane(x))
addRequired(p, 'siz', @(x)validateattributes(x,{'numeric'},...
    {'size',[1, nan],'positive','nonnan','real','finite'}))
parse(p, plane, siz)

if ~isempty(varargin)
    if length(varargin) == 1
        if isstruct(varargin{1})
            % if options are specified as struct, need to convert to
            % parameter name-value pairs
            varargin = [fieldnames(varargin{1}) struct2cell(varargin{1})]';
            varargin = varargin(:)';
        else
            % if option is a single argument, assume it corresponds to
            % plane color
            varargin = {'FaceColor', varargin{1}};
        end
    end
else
    % default face color
    varargin = {'FaceColor', 'm'};
end

if numel(siz) == 1
    siz(2) = siz(1);
end

%% Algorithm
% Calculate vertex points of the platform
pts(1,:) = planePoint(plane, [1,1]*0.5.*siz);
pts(2,:) = planePoint(plane, [1,-1]*0.5.*siz);
pts(3,:) = planePoint(plane, [-1,-1]*0.5.*siz);
pts(4,:) = planePoint(plane, [-1,1]*0.5.*siz);

pf.vertices = pts;
pf.faces = [1 2 3 4];

% Draw the patch
h = patch(hAx, pf, varargin{:});

%% Parse outputs
% Return handle to plane if needed
if nargout > 0
    varargout{1} = h;
end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [logp, yhat, res] = tapas_condhalluc_obs(r, infStates, ptrans)
% Calculates the log-probability of response y=1 under the unit-square sigmoid model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2015-2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (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/>.

% Transform beta to its native space
be = exp(ptrans(1));

% Initialize returned log-probabilities as NaNs so that NaN is
% returned for all irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Check input format
if size(r.u,2) ~= 2
    error('tapas:hgf:CondHalluc:InputsIncompatible', 'Inputs incompatible with condhalluc_obs observation model. See tapas_condhalluc_obs_config.m.')
end

% Get true-positive rate corresponding to stimuli
tp = r.u(:,2);

% Weed irregular trials out
mu1hat = infStates(:,1,1);
mu1hat(r.irr) = [];
y = r.y(:,1);
y(r.irr) = [];
tp(r.irr) = [];

% Calculate belief x using Bayes' theorem
x = tp.*mu1hat./(tp.*mu1hat + (1-mu1hat).^2);

% Belief is mu1hat in trials where there is no tone
x(find(tp==0)) = mu1hat(find(tp==0));

% Calculate log-probabilities for non-irregular trials
reg = ~ismember(1:n,r.irr);
logp(reg) = -log(1+exp(-be.*(2.*x-1).*(2.*y-1)));
yhat(reg) = x;
res(reg) = (y-x)./sqrt(x.*(1-x));

return;

author: translationalneuromodeling

repo: tapas


text:
function [data_mean,data_diff,md,sd] = bland_altman_x_gt(data1,data2)
% Function to generate Bland Altman plots. Barry Greene, September 2008
% Bland, J.M., Altman, D.G. 'Statistical methods for assessing agreement ...
% between two methods of clinical measurement'(1986) Lancet, 1 (8476), pp. 307-310.
%
% Inputs: data1: ground truth
%         data2: estimated
%
% Produces Bland Altman plot with mean difference and mean difference +/-
% 2*SD difference lines.

[m,n] = size(data1);
if(n>m)
    data1 = data1';
end

if(size(data1)~=size(data2))
    error('Data matrices must be the same size')
end

data_mean = data1;  % Mean of values from each instrument
data_diff = data2 - data1;              % Difference between data from each instrument
md = mean(data_diff);               % Mean of difference between instruments
sd = std(data_diff);                % Std dev of difference between instruments

% figure;
plot(data_mean,data_diff,'ok','MarkerSize',3,'LineWidth',1);   % Bland Altman plot
hold on; x = [min(data_mean),max(data_mean)];
y = md*ones(1,2); plot(x,y,'-k'); % Mean difference line
y = md+2*sd*ones(1,2);
plot(x,y,'--k'); % Mean plus 2*SD line
% text(x(2),y(2),'+2 SD');
y = md-2*sd*ones(1,2);
plot(x,y,'--k'); % Mean minus 2*SD line

% text(x(2),y(2),'-2 SD');
% grid on
% title('Bland Altman plot','FontSize',9)
% xlabel('Mean of two measures','FontSize',8)
% ylabel('Difference between two measures','FontSize',8)

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [surfhandle, T2] = ft_plot_slice(dat, varargin)

% FT_PLOT_SLICE plots a single slice that cuts through a 3-D volume and interpolates
% the data if needed.
%
% Use as
%   ft_plot_slice(dat, ...)
% or
%   ft_plot_slice(dat, mask, ...)
% where dat and mask are equal-sized 3-D arrays.
%
% Additional options should be specified in key-value pairs and can be
%   'transform'    = 4x4 homogeneous transformation matrix specifying the mapping from
%                    voxel coordinates to the coordinate system in which the data are plotted.
%   'location'     = 1x3 vector specifying a point on the plane which will be plotted
%                    the coordinates are expressed in the coordinate system in which the
%                    data will be plotted. location defines the origin of the plane
%   'orientation'  = 1x3 vector specifying the direction orthogonal through the plane
%                    which will be plotted (default = [0 0 1])
%   'unit'         = string, can be 'm', 'cm' or 'mm (default is automatic)
%   'coordsys'     = string, assume the data to be in the specified coordinate system (default = 'unknown')
%   'resolution'   = number (default = 1 mm)
%   'datmask'      = 3D-matrix with the same size as the data matrix, serving as opacitymap
%                    If the second input argument to the function contains a matrix, this
%                    will be used as the mask
%   'maskstyle'    = string, 'opacity' or 'colormix', defines the rendering
%   'background'   = needed when maskstyle is 'colormix', 3D-matrix with
%                    the same size as the data matrix, serving as
%                    grayscale image that provides the background
%   'opacitylim'   = 1x2 vector specifying the limits for opacity masking
%   'interpmethod' = string specifying the method for the interpolation, see INTERPN (default = 'nearest')
%   'colormap'     = string, see COLORMAP
%   'clim'         = 1x2 vector specifying the min and max for the colorscale
%   'facealpha'    = transparency when no mask is specified, between 0 and 1 (default = 1)
%
% You can plot the slices from the volume together with an intersection of the slices
% with a triangulated surface mesh (e.g. a cortical sheet) using
%   'intersectmesh'       = triangulated mesh, see FT_PREPARE_MESH
%   'intersectcolor'      = string, color specification
%   'intersectlinestyle'  = string, line specification
%   'intersectlinewidth'  = number
%
% See also FT_PLOT_ORTHO, FT_PLOT_MONTAGE, FT_SOURCEPLOT

% Undocumented options
%   'plotmarker'     = Nx3 matrix with points to be plotted as markers, e.g. dipole positions
%   'markersize'
%   'markercolor'

% Copyrights (C) 2010-2014, Jan-Mathijs Schoffelen
% Copyrights (C) 2014-2022, Robert Oostenveld and Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

persistent dim X Y Z

if isstruct(dat) && isfield(dat, 'anatomy') && isfield(dat, 'transform')
  % the input is an MRI structure, call this function recursively
  varargin = ft_setopt(varargin, 'transform', dat.transform);
  if isfield(dat, 'coordsys')
    varargin = ft_setopt(varargin, 'coordsys', dat.coordsys);
  end
  if isfield(dat, 'unit')
    varargin = ft_setopt(varargin, 'unit', dat.unit);
  end
  dat = dat.anatomy;
  [surfhandle, T2] = ft_plot_slice(dat, varargin{:});
  return
end

% the data can have up to 5 dimensions, including time and/or frequency
% the size function in MATLAB 2022b can take multiple input arguments, but not in 2018b
datsize = [size(dat, 1) size(dat, 2) size(dat, 3) size(dat, 4) size(dat, 5)];

if isequal(dim, datsize(1:3))
  % reuse the persistent variables to speed up subsequent calls with the same input
else
  % construct the persistent variables
  dim = datsize(1:3);
  [X, Y, Z] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));
end

if any(dim==1)
  ft_error('it is not possible to plot a volume that consists of a single slice');
end

% parse first input argument(s), it is either
% (dat, varargin)
% (dat, msk, varargin)
% (dat, [], varargin)

if numel(varargin)>0 && (isempty(varargin{1}) || isnumeric(varargin{1}) || islogical(varargin{1}))
  msk      = varargin{1};
  varargin = varargin(2:end);
end

% get the optional input arguments
transform           = ft_getopt(varargin, 'transform', eye(4));
loc                 = ft_getopt(varargin, 'location');
ori                 = ft_getopt(varargin, 'orientation', [0 0 1]);
coordsys            = ft_getopt(varargin, 'coordsys');
unit                = ft_getopt(varargin, 'unit');       % the default will be determined further down
resolution          = ft_getopt(varargin, 'resolution'); % the default depends on the units and will be determined further down
datmask             = ft_getopt(varargin, 'datmask');
maskstyle           = ft_getopt(varargin, 'maskstyle', 'opacity');
background          = ft_getopt(varargin, 'background');
opacitylim          = ft_getopt(varargin, 'opacitylim');
interpmethod        = ft_getopt(varargin, 'interpmethod', 'nearest');
cmap                = ft_getopt(varargin, 'colormap');
clim                = ft_getopt(varargin, 'clim');
doscale             = ft_getopt(varargin, 'doscale', true); % only scale when necessary (time consuming), i.e. when plotting as grayscale image & when the values are not between 0 and 1
surfhandle          = ft_getopt(varargin, 'surfhandle', []);
patchhandle         = ft_getopt(varargin, 'patchhandle', []);
tag                 = ft_getopt(varargin, 'tag', ''); % this is used to keep the thee intersecting slices or subplots apart

mesh                = ft_getopt(varargin, 'intersectmesh');
intersectcolor      = ft_getopt(varargin, 'intersectcolor', 'yrgbmyrgbm');
intersectlinewidth  = ft_getopt(varargin, 'intersectlinewidth', 2);
intersectlinestyle  = ft_getopt(varargin, 'intersectlinestyle');

plotmarker          = ft_getopt(varargin, 'plotmarker');
markersize          = ft_getopt(varargin, 'markersize', 'auto');
markercolor         = ft_getopt(varargin, 'markercolor', 'w');
facealpha           = ft_getopt(varargin, 'facealpha', 1); % only applies when no mask is specified

% convert from yes/no/true/false/0/1 into a proper boolean
doscale = istrue(doscale);

if ~isa(dat, 'double')
  dat = cast(dat, 'double');
end

if exist('msk', 'var') && isempty(datmask)
  ft_warning('using the second input argument as mask rather than the one from the varargin list');
  datmask = msk; clear msk;
end

% normalise the orientation vector to one
ori = ori./sqrt(sum(ori.^2));

% set the default location
if isempty(loc) && (isempty(transform) || isequal(transform, eye(4)))
  loc = (dim+1)./2;
elseif isempty(loc)
  loc = [0 0 0];
end

% shift the location to be along the orientation vector
loc = ori*dot(loc,ori);

dointersect = ~isempty(mesh);

% the mesh should be a cell-array
if isstruct(mesh)
  tmp = mesh;
  mesh = cell(size(tmp));
  for i=1:numel(tmp)
    mesh{i} = tmp(i);
  end
elseif isempty(mesh)
  mesh = {};
end

% replace pnt by pos
for k = 1:numel(mesh)
  mesh{k} = fixpos(mesh{k});
end

% the mesh should be a structure with both pos and tri
for k = 1:numel(mesh)
  if ~isfield(mesh{k}, 'pos') || ~isfield(mesh{k}, 'tri')
    mesh{k}.pos = [];
    mesh{k}.tri = [];
  end
end

domask = ~isempty(datmask);
if domask
  % check whether the mask is ok
  if ~isequal(size(dat), size(datmask)) && ~isequal(cmap, 'rgb')
    % the exception is when the functional data is to be interpreted as rgb
    ft_error('the mask data should have the same dimensions as the functional data');
  end
end

dobackground = ~isempty(background);
if dobackground
  % check whether the background is ok
  if ~isequal(size(dat), size(background))
    error('the background data should have the same dimensions as the functional data');
  end
end

% determine the voxel center
% voxel_center_vc = [X(:) Y(:) Z(:)];
% voxel_center_hc = ft_warp_apply(transform, voxel_center_vc);

% determine the edges, i.e. the corner points of each voxel
% [Xe, Ye, Ze] = ndgrid(0:dim(1), 0:dim(2), 0:dim(3));
% Xe = Xe+0.5;
% Ye = Ye+0.5;
% Ze = Ze+0.5;
% voxel_edge_vc = [Xe(:) Ye(:) Ze(:)];
% voxel_edge_hc = ft_warp_apply(transform, voxel_edge_vc);

% determine the corner points of the box that encompasses the whole data
% extend the box with half a voxel in all directions to get the outer edge
corner_vc = [
  0.5        0.5        0.5
  0.5+dim(1) 0.5        0.5
  0.5+dim(1) 0.5+dim(2) 0.5
  0.5        0.5+dim(2) 0.5
  0.5        0.5        0.5+dim(3)
  0.5+dim(1) 0.5        0.5+dim(3)
  0.5+dim(1) 0.5+dim(2) 0.5+dim(3)
  0.5        0.5+dim(2) 0.5+dim(3)
  ];
corner_hc = ft_warp_apply(transform, corner_vc);

if isempty(unit)
  if ~isequal(transform, eye(4))
    % estimate the geometrical units we are dealing with
    unit = ft_estimate_units(norm(range(corner_hc)));
  else
    % units are in voxels, these are assumed to be close to mm
    unit = 'mm';
  end
end

if isempty(resolution)
  % the default resolution is 1 mm
  resolution = ft_scalingfactor('mm', unit);
end

% determine whether interpolation is needed
dointerp = false;
dointerp = dointerp || sum(sum(transform-eye(4)))~=0;
dointerp = dointerp || ~all(round(loc)==loc);
dointerp = dointerp || sum(ori)~=1;
dointerp = dointerp || ~(resolution==round(resolution));
% determine the caller function and toggle dointerp to true, if ft_plot_slice has been called from ft_plot_montage
% this is necessary for the correct allocation of the persistent variables
st = dbstack;
if ~dointerp && numel(st)>1 && strcmp(st(2).name, 'ft_plot_montage'), dointerp = true; end

% define 'x' and 'y' axis in projection plane, the definition of x and y is more or less arbitrary
[x, y] = projplane(ori);
% z = ori;

% project the corner points onto the projection plane
corner_pc = zeros(size(corner_hc));
for i=1:8
  corner   = corner_hc(i, :) - loc(:)';
  corner_pc(i,1) = dot(corner, x);
  corner_pc(i,2) = dot(corner, y);
  corner_pc(i,3) = 0;
end

% get the transformation matrix from the projection plane to head coordinates
T2 = [x(:) y(:) ori(:) loc(:); 0 0 0 1];

% get the transformation matrix from projection plane to voxel coordinates
T3 = transform\T2;

min_corner_pc = min(corner_pc, [], 1);
max_corner_pc = max(corner_pc, [], 1);
% round the bounding box limits to the nearest mm
switch unit
  case 'm'
    min_corner_pc = ceil(min_corner_pc*100)/100;
    max_corner_pc = floor(max_corner_pc*100)/100;
  case 'cm'
    min_corner_pc = ceil(min_corner_pc*10)/10;
    max_corner_pc = floor(max_corner_pc*10)/10;
  case 'mm'
    min_corner_pc = ceil(min_corner_pc);
    max_corner_pc = floor(max_corner_pc);
end

% determine a grid of points in the projection plane
xplane = min_corner_pc(1):resolution:max_corner_pc(1);
yplane = min_corner_pc(2):resolution:max_corner_pc(2);
zplane = 0;
[Xi, Yi, Zi]      = ndgrid(xplane, yplane, zplane);
siz               = [size(squeeze(Xi)) size(dat,4)];
interp_center_pc  = [Xi(:) Yi(:) Zi(:)];
% interp_center_hc = ft_warp_apply(T2, interp_center_pc);

% get the positions of the points in the projection plane in voxel coordinates
interp_center_vc = ft_warp_apply(T3, interp_center_pc);

Xi = reshape(interp_center_vc(:, 1), siz(1:2));
Yi = reshape(interp_center_vc(:, 2), siz(1:2));
Zi = reshape(interp_center_vc(:, 3), siz(1:2));

% check whether the values in the axes are close enough to integer
tol = nanmean([diff(unique(Xi(:)));diff(unique(Yi(:)))])./100;
isintegerXi = issufficientlyinteger(Xi(:),tol);
isintegerYi = issufficientlyinteger(Yi(:),tol);
isintegerZi = issufficientlyinteger(Zi(:),tol);

% check whether it's possible to select an orthogonal plane
[islineXi, lineXi] = isline(Xi);
[islineYi, lineYi] = isline(Yi);
[islineZi, lineZi] = isline(Zi);

use_interpn = ~isequal(transform, eye(4)) || ~isequal(interpmethod, 'nearest') || ~all([isintegerXi isintegerYi isintegerZi]);
get_slice   = ~use_interpn && all([islineXi islineYi islineZi]) && all([isintegerXi isintegerYi isintegerZi]);
if use_interpn
  V  = interpn(X, Y, Z, dat, Xi, Yi, Zi, interpmethod);
  if domask,       Vmask = interpn(X, Y, Z, datmask,    Xi, Yi, Zi, interpmethod); end
  if dobackground, Vback = interpn(X, Y, Z, background, Xi, Yi, Zi, interpmethod); end
elseif get_slice
  % something more efficient than an interpolation can be done:
  % just select the appropriate plane, and permute to get the orientation
  % right in the plots, this has something to do with ndgrid vs meshgrid I think
  permutevec = [2 1];
  if ndims(dat)>3
    permutevec = [permutevec 3:ndims(dat)];
  end
  if numel(unique(lineXi(:)))==1
    lineXi = lineXi(1);
  elseif numel(unique(lineYi(:)))==1
    lineYi = lineYi(1);
  elseif numel(unique(lineZi(:)))==1
    lineZi = lineZi(1);
  end
  V = permute(reshape(dat(lineXi,lineYi,lineZi,:), siz(permutevec(1:ndims(dat)-1))), permutevec);
  if domask,       Vmask = permute(reshape(datmask(lineXi,lineYi,lineZi,:),       siz(permutevec(1:2))), [2 1]); end
  if dobackground, Vback = permute(reshape(background(lineXi,lineYi,lineZi,:), siz(permutevec(1:2))), [2 1]); end
else
  % use sub2ind in the unlikely case that it's an oblique plane, parallel
  % to one of the axes with only integer indices
  % this fails for rgb data
  V = dat(sub2ind(dim, Xi(:), Yi(:), Zi(:)));
  V = reshape(V, siz);
end

if all(isnan(V(:)))
  % the projection plane lies completely outside the box spanned by the data
else
  % trim the edges of the projection plane
  [sel1, sel2] = tight(V(:,:,1));
  V  = V (sel1,sel2,:);
  Xi = Xi(sel1,sel2);
  Yi = Yi(sel1,sel2);
  Zi = Zi(sel1,sel2);
  if domask
    Vmask = Vmask(sel1,sel2);
  end
  if dobackground
    Vback = Vback(sel1,sel2);
  end
end

if dobackground
  % convert the background plane to a grayscale image
  bmin  = nanmin(background(:));
  bmax  = nanmax(background(:));
  Vback = (Vback-bmin)./(bmax-bmin);
  Vback(~isfinite(Vback)) = 0;
  Vback = cat(3, Vback, Vback, Vback);
end

interp_center_vc = [Xi(:) Yi(:) Zi(:)]; clear Xi Yi Zi
interp_center_pc = ft_warp_apply(inv(T3), interp_center_vc);

% determine a grid of points in the projection plane
% this reconstruction is needed since the edges may have been trimmed off
xplane = min(interp_center_pc(:, 1)):resolution:max(interp_center_pc(:, 1));
yplane = min(interp_center_pc(:, 2)):resolution:max(interp_center_pc(:, 2));
zplane = 0;

[Xi, Yi, Zi] = ndgrid(xplane, yplane, zplane); % 2D cartesian grid of projection plane in plane voxels
siz          = size(squeeze(Xi));

% extend with one voxel along dim 1
Xi = cat(1, Xi, Xi(end,:)+mean(diff(Xi,[],1),1));
Yi = cat(1, Yi, Yi(end,:)+mean(diff(Yi,[],1),1));
Zi = cat(1, Zi, Zi(end,:)+mean(diff(Zi,[],1),1));
% extend with one voxel along dim 2
Xi = cat(2, Xi, Xi(:,end)+mean(diff(Xi,[],2),2));
Yi = cat(2, Yi, Yi(:,end)+mean(diff(Yi,[],2),2));
Zi = cat(2, Zi, Zi(:,end)+mean(diff(Zi,[],2),2));
% shift with half a voxel along dim 1 and 2
Xi = Xi-0.5*resolution;
Yi = Yi-0.5*resolution;
% Zi = Zi; % do not shift along this direction

interp_edge_pc = [Xi(:) Yi(:) Zi(:)]; clear Xi Yi Zi
interp_edge_hc = ft_warp_apply(T2, interp_edge_pc);

if false
  % plot all objects in head coordinates
  ft_plot_mesh(voxel_center_hc,   'vertexmarker', 'o')
  ft_plot_mesh(voxel_edge_hc,     'vertexmarker', '+')
  ft_plot_mesh(corner_hc,         'vertexmarker', '*')
  ft_plot_mesh(interp_center_hc,  'vertexmarker', 'o', 'vertexcolor', 'r')
  ft_plot_mesh(interp_edge_hc,    'vertexmarker', '+', 'vertexcolor', 'r')
  axis on
  grid on
  xlabel('x')
  ylabel('y')
  zlabel('z')
end

if isempty(cmap)
  % treat as gray value: scale and convert to rgb
  if doscale
    dmin = min(dat(:));
    dmax = max(dat(:));
    V    = (V-dmin)./(dmax-dmin);
    clear dmin dmax
  end
  V(~isfinite(V)) = 0;

  % deal with clim for RGB data here, where the purpose is to increase the
  % contrast range, rather than shift the average grey value
  if ~isempty(clim)
    V = (V-clim(1))./clim(2);
    V(V>1)=1;
  end

  % convert into RGB values, e.g. for the plotting of anatomy
  V = cat(3, V, V, V);
end

% get positions of the voxels in the interpolation plane in head coordinates
Xh = reshape(interp_edge_hc(:,1), siz+1);
Yh = reshape(interp_edge_hc(:,2), siz+1);
Zh = reshape(interp_edge_hc(:,3), siz+1);

% do the actual plotting of the slice
if ~domask
  % no masked slice to be plotted
  if isempty(surfhandle)
    % create surface object
    surfhandle = surface(Xh, Yh, Zh, V);
    set(surfhandle, 'linestyle', 'none');
    set(surfhandle, 'FaceAlpha', facealpha);
  else
    % update the colordata in the surface object
    set(surfhandle, 'Cdata', V);
    set(surfhandle, 'Xdata', Xh);
    set(surfhandle, 'Ydata', Yh);
    set(surfhandle, 'Zdata', Zh);
    set(surfhandle, 'FaceAlpha', facealpha);
  end
elseif domask
  % what should be done depends on the maskstyle
  switch maskstyle
    case 'opacity'
      if dobackground
        ft_warning('specifying maskstyle = ''opacity'' causes the supplied background image not to be used');
      end
      if isempty(surfhandle)
        % create surface object
        surfhandle = surface(Xh, Yh, Zh, V);
        set(surfhandle, 'linestyle', 'none');
      else
        % update the colordata in the surface object
        set(surfhandle, 'Cdata', V);
        set(surfhandle, 'Xdata', Xh);
        set(surfhandle, 'Ydata', Yh);
        set(surfhandle, 'Zdata', Zh);
      end
      if islogical(Vmask), Vmask = double(Vmask); end
      set(surfhandle, 'FaceColor', 'texture');
      set(surfhandle, 'FaceAlpha', 'texturemap'); %flat
      set(surfhandle, 'AlphaDataMapping', 'scaled');
      set(surfhandle, 'AlphaData', Vmask);
      if ~isempty(opacitylim)
        alim(opacitylim)
      end
    case 'colormix'
      if isempty(cmap), error('using ''colormix'' as maskstyle requires an explicitly defined colormap'); end
      if ischar(cmap),  cmap = strrep(cmap, 'default', 'parula'); cmap = ft_colormap(cmap); end
      V = bg_rgba2rgb(Vback,V,cmap,clim,Vmask,'rampup',opacitylim);
      if isempty(surfhandle)
        % create surface object
        surfhandle = surface(Xh, Yh, Zh, V);
        set(surfhandle, 'linestyle', 'none');
      else
        % update the colordata in the surface object
        set(surfhandle, 'Cdata', V);
        set(surfhandle, 'Xdata', Xh);
        set(surfhandle, 'Ydata', Yh);
        set(surfhandle, 'Zdata', Zh);
      end
    otherwise
      error('unsupported maskstyle');
  end
end

if ~isempty(coordsys) && ~strcmp(coordsys, 'unknown')
  % convert 'neuromag' to 'ras', etc.
  coordsys = generic(coordsys);

  % determine the center of each of the 6 faces of the box that encompasses the whole data
  % shift all of them inward by about 5%
  center_vc = [
    dim(1)*0.05 dim(2)*0.50 dim(3)*0.50
    dim(1)*0.95 dim(2)*0.50 dim(3)*0.50
    dim(1)*0.50 dim(2)*0.05 dim(3)*0.50
    dim(1)*0.50 dim(2)*0.95 dim(3)*0.50
    dim(1)*0.50 dim(2)*0.50 dim(3)*0.05
    dim(1)*0.50 dim(2)*0.50 dim(3)*0.95
    ];
  center_hc = ft_warp_apply(transform, center_vc);

  % the order of the center points is arbitrary and does not match coordsys/xyz
  % hence we have to determine them for each x/y/z direction
  minx = min(center_hc(:,1));
  maxx = max(center_hc(:,1));
  miny = min(center_hc(:,2));
  maxy = max(center_hc(:,2));
  minz = min(center_hc(:,3));
  maxz = max(center_hc(:,3));
  % also compute the midpoints, i.e. the center of the cube itself
  midx = mean(center_hc(:,1));
  midy = mean(center_hc(:,2));
  midz = mean(center_hc(:,3));

  if isequal(ori, [1 0 0])
    % the slice is perpendicular to the x-axis
    delete(findall(gca, 'Tag', ['coordsyslabel_' tag])) % remove the labels from the previous call
    text(loc(1), miny, midz, upper(flipletter(coordsys(2))), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(loc(1), maxy, midz, upper(           coordsys(2) ), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(loc(1), midy, minz, upper(flipletter(coordsys(3))), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(loc(1), midy, maxz, upper(           coordsys(3) ), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
  elseif isequal(ori, [0 1 0])
    % the slice is perpendicular to the y-axis
    delete(findall(gca, 'Tag', ['coordsyslabel_' tag])) % remove the labels from the previous call
    text(minx, loc(2), midz, upper(flipletter(coordsys(1))), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(maxx, loc(2), midz, upper(           coordsys(1) ), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(midx, loc(2), minz, upper(flipletter(coordsys(3))), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(midx, loc(2), maxz, upper(           coordsys(3) ), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
  elseif isequal(ori, [0 0 1])
    % the slice is perpendicular to the z-axis
    delete(findall(gca, 'Tag', ['coordsyslabel_' tag])) % remove the labels from the previous call
    text(minx, midy, loc(3), upper(flipletter(coordsys(1))), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(maxx, midy, loc(3), upper(           coordsys(1) ), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(midx, miny, loc(3), upper(flipletter(coordsys(2))), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
    text(midx, maxy, loc(3), upper(           coordsys(2) ), 'Color', 'y', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Tag', ['coordsyslabel_' tag]);
  end
end % if coordsys

% plot the intersection with a mesh
if dointersect
  % determine three points on the plane
  inplane = eye(3) - (eye(3) * ori') * ori;
  v1 = loc + inplane(1,:);
  v2 = loc + inplane(2,:);
  v3 = loc + inplane(3,:);

  if isempty(patchhandle) || length(patchhandle)~=length(mesh)
    % try to find the handles of all patches
    patchhandle = findall(gca, 'tag', ['intersectmesh_' tag]);
  end

  if length(patchhandle)~=length(mesh)
    % the patch handles do not make sense, start from scratch
    delete(findall(gca, 'tag', ['intersectmesh_' tag]))
    patchhandle = nan(size(mesh));
  end

  for k = 1:numel(mesh)
    [xmesh, ymesh, zmesh] = intersect_plane(mesh{k}.pos, mesh{k}.tri, v1, v2, v3);

    % draw each individual line segment of the intersection
    if ~isempty(xmesh)
      if ~ishandle(patchhandle(k))
        patchhandle(k) = patch(xmesh', ymesh', zmesh', nan(1, size(xmesh, 1)));
        set(patchhandle(k), 'tag', ['intersectmesh_' tag]);
        if ~isempty(intersectcolor),     set(patchhandle(k), 'EdgeColor', intersectcolor(k));  end
        if ~isempty(intersectlinewidth), set(patchhandle(k), 'LineWidth', intersectlinewidth); end
        if ~isempty(intersectlinestyle), set(patchhandle(k), 'LineStyle', intersectlinestyle); end
      else
        set(patchhandle(k), 'XData', xmesh', 'YData', ymesh', 'ZData', zmesh', 'FaceVertexCdata', nan(size(xmesh,1),1));
      end
    end
  end
end

if ~isempty(cmap) && ~isequal(cmap, 'rgb') && ~isequal(maskstyle, 'colormix')
  ft_colormap(cmap);
  if ~isempty(clim)
    caxis(clim);
  end
end

if ~isempty(plotmarker)
  % determine three points on the plane
  inplane = eye(3) - (eye(3) * ori') * ori;
  v1 = loc + inplane(1,:);
  v2 = loc + inplane(2,:);
  v3 = loc + inplane(3,:);
  pr = nan(size(plotmarker,1), 3);
  d  = nan(size(plotmarker,1), 1);
  for k = 1:size(plotmarker,1)
    [pr(k,:), d(k,:)] = ptriprojn(v1, v2, v3, plotmarker(k,:));
  end
  sel = d<eps*1e8;
  if sum(sel)>0
    ft_plot_dipole(pr(sel,:), repmat([0;0;1], 1, size(pr,1)), 'length', 0, 'color', markercolor, 'diameter', markersize);
  end
end

% update the axes to ensure that the whole volume fits
ax = [min(corner_hc) max(corner_hc)];
axis(ax([1 4 2 5 3 6])); % reorder into [xmin xmax ymin ymaz zmin zmax]

st = dbstack;
if numel(st)>1
  % ft_plot_slice has been called from another function, probably ft_plot_ortho
  % assume the remainder of the axis settings to be handled there
else
  set(gca,'xlim',[min(Xh(:))-0.5*resolution max(Xh(:))+0.5*resolution]);
  set(gca,'ylim',[min(Yh(:))-0.5*resolution max(Yh(:))+0.5*resolution]);
  set(gca,'zlim',[min(Zh(:))-0.5*resolution max(Zh(:))+0.5*resolution]);

  set(gca,'dataaspectratio',[1 1 1]);
  % axis equal; % this for some reason does not work robustly when drawing intersections, replaced by the above
  axis vis3d
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y] = projplane(z)
[u, s, v] = svd([eye(3) z(:)]);
x = u(:, 2)';
y = u(:, 3)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [sel1, sel2] = tight(V)
% make a selection to cut off the nans at the edges
sel1 = sum(~isfinite(V), 2)<size(V, 2);
sel2 = sum(~isfinite(V), 1)<size(V, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bool = issufficientlyinteger(X, tolerance)
%isinteger only checks for integer class, so will always return false with
%double integers
bool = all(abs(X-round(X))<tolerance);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bool, lineX] = isline(X)
%isline returns an array of the values in X are columnwise or rowwise the
%same, otherwise returns false
if isequal(X(ones(1,size(X,1)),:),X)
  lineX = X(1,:);
  bool = true;
elseif isequal(X(:,ones(1,size(X,2))),X)
  lineX = X(:,1)';
  bool = true;
else
  lineX = [];
  bool = false;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function letter = flipletter(letter)
switch letter
  case 'a'
    letter = 'p';
  case 'p'
    letter = 'a';
  case 'l'
    letter = 'r';
  case 'r'
    letter = 'l';
  case 'i'
    letter = 's';
  case 's'
    letter = 'i';
  otherwise
    ft_error('incorrect letter')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function coordsys = generic(coordsys)
mapping = {
  'ctf',       'als'
  'bti',       'als'
  '4d',        'als'
  'yokogawa',  'als'
  'eeglab',    'als'
  'neuromag',  'ras'
  'itab',      'ras'
  'acpc',      'ras'
  'spm',       'ras'
  'mni',       'ras'
  'fsaverage', 'ras'
  'tal',       'ras'
  'scanras',   'ras'
  'scanlps',   'lps'
  'dicom',     'lps'
  };

sel = find(strcmp(mapping(:,1), coordsys));
if length(sel)==1
  coordsys = mapping{sel,2};
end

if ~all(ismember(coordsys, 'lrapis'))
  ft_error('cannot convert "%s" to a generic coordinate system label', coordsys);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [keypointsAll,bboxAll] = backprojectKeypoints(expidx)

fprintf('backprojectKeypoints()\n');

% load annolist with full size images
[annolistOrig, imgidxs] = getAnnolist(expidx);
for imgidx = 1:length(annolistOrig)
    annolistOrig(imgidx).annorect = [];
end

p = rcnn_exp_params(expidx);
load(p.testGT);
if (~exist('annolist','var'))
    annolist = single_person_annolist;
end
annolist_crops = annolist;

[~,parts] = util_get_parts24();

% load bboxes
bboxAll = getBBoxAll(expidx,parts,length(annolist_crops));
imgidxs_orig = zeros(length(bboxAll),1);
jidxs = [1:6 9:length(bboxAll(1).det)];

fprintf('backproject bboxes()\n');
for imgidx = 1:length(bboxAll)
    fprintf('.');
    [~,name] = fileparts(annolist_crops(imgidx).image.name);
    load([fileparts(p.testGT) '/T_' name(3:end)],'T');
    imgidxs_orig(imgidx) = str2double(name(3:7));

    det = bboxAll(imgidx).det;
    for jidx = jidxs
        detNew1 = ([det{jidx}(:,1:2) ones(size(det{jidx},1),1)]*T');
        detNew2 = ([det{jidx}(:,3:4) ones(size(det{jidx},1),1)]*T');
        det{jidx}(:,1:4) = [detNew1(:,1:2) detNew2(:,1:2)];
        %             figure(101); clf; imagesc(imread(annolist(imgidx_orig).image.name)); axis equal; hold on;
        %             plot(det{jidx}(:,1),det{jidx}(:,2),'ro','MarkerFaceColor','r','MarkerEdgeColor','k','MarkerSize',5);
    end

    bboxAll(imgidx).det = det;
    bboxAll(imgidx).imgname = annolistOrig(imgidxs_orig(imgidx)).image.name;
    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(bboxAll));
    end
end
fprintf(' done\n');

% merge detections
imgidxs_merged = imgidxs_orig(1);
imgidxs_merged_rel = 1;
for imgidx = 2:length(bboxAll)
    if (imgidxs_orig(imgidx) == imgidxs_merged(end))
        for jidx = jidxs
            bboxAll(imgidxs_merged_rel(end)).det{jidx} = [bboxAll(imgidxs_merged_rel(end)).det{jidx}; bboxAll(imgidx).det{jidx}];
        end
    else
        imgidxs_merged = [imgidxs_merged; imgidxs_orig(imgidx)];
        imgidxs_merged_rel = [imgidxs_merged_rel; imgidx];
    end
end
bboxAll = bboxAll(imgidxs_merged_rel);

fprintf('nms()\n');
for imgidx = 1:length(bboxAll)
    fprintf('.');
    for jidx = jidxs
        keep = nms(bboxAll(imgidx).det{jidx}, 0.5);
        bboxAll(imgidx).det{jidx} = bboxAll(imgidx).det{jidx}(keep, :);
    end
    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(bboxAll));
    end
end
fprintf(' done\n');

for imgidx = 1:length(bboxAll)
    keypointsAll(imgidx).imgname = bboxAll(imgidx).imgname;
    keypointsAll(imgidx).det = cell(16,1);
end

% convert to keypoints
for imgidx = 1:length(bboxAll)
    points = nan(16,2);
    for jidx = jidxs
        det = bboxAll(imgidx).det{jidx};
        x = mean(det(:,[1 3]),2);
        y = mean(det(:,[2 4]),2);
        keypointsAll(imgidx).det{jidx,:} = [[x y] det(:,5)];
        [val,idx] = max(det(:,5));
        points(jidx,:) = [x(idx) y(idx)];
    end
%     figure(101); clf; imagesc(imread(keypointsAll(imgidx).imgname)); axis equal; hold on;
%     plot(points(:,1),points(:,2),'ro','MarkerFaceColor','r','MarkerEdgeColor','k','MarkerSize',5);
%     set(gca,'Ydir','reverse');
end

path = fileparts(p.evalTest);
fnameKeypoints = [path '/keypointsAll'];
save(fnameKeypoints, 'keypointsAll');

fnameBBox = [path '/bboxAll'];
save(fnameBBox, 'bboxAll');
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut