text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal, Marc Bakry (c) 2017-2019.                 |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             marc.bakry@polytechnique.edu                               |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtHmxVibroSlab.m                             |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Marc Bakry                  |
%|  ( # )  |   CREATION   : 14.03.2019                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
tol = 1e-3
X0  = [0 0 1]
e   = 0.5;
f   = 900;

% Exterior domain (water)
rho0 = 1000;         % density      (kg.m3)
c0   = 1500;         % celerity     (m.s-1)
k0   = 2*pi/c0*f;    % wave-number  (m-1)
lam0 = c0/f;         % wave-length  (m)

% Interior domain (different from water)
rhoS = 2*rho0;                                 % density                             (kg.m3)
cL   = 2*c0;                                   % celerity of longitudinal waves      (m.s-1)
kL   = 2*pi*f/cL;                              % wave-number of longitudinal waves   (m-1)
lamL = real(cL)/f;                             % wavelength of longitudinal waves    (m)
cT   = 0;                                      % celerity of transverse waves        (m.s-1)
kT   = 2*pi.*f/cT;                             % wave-number of transverse waves     (m-1)
lamT = real(cT)/f;                             % wavelength of transverse waves

% Minimum wavelength
tmp  = [lam0,lamL,lamT];
lmin = min(tmp(tmp>0));

% Slab mesh
L    = 20 * lmin;            % 20 wavelength to simulate infinite slab
nx   = ceil(L/lmin * 6)+1;   % 6 node per wavelength for L
ny   = ceil(L/lmin * 6)+1;   % 6 node per wavelength for L
nz   = ceil(e/lmin * 12)+1;  % 12 node per wavelength for e
N    = nx * ny * nz;         % Total number of nodes
mesh = mshCube(N,[L L e])

% Boundary
bound = mesh.bnd

% Admissible frequency
stp  = mesh.stp;
kmax = 1/stp(2);
if (k0 > kmax) || (kL > kmax/2) || ((kT>kmax/2) && ~isinf(kT))
%     error('frequency is too high for mesh discretization')
end

% Radiative mesh (fixed number of nodes)
radiatx     = mshSquare(1e3,[L L]);
radiaty     = radiatx;
radiatx.vtx = radiatx.vtx(:,[1 3 2]);
radiaty.vtx = radiaty.vtx(:,[3 1 2]);
radiat      = union(radiatx,radiaty);

% Measurement points for trans and refl coeff (1 wavelenth from bound)
Xmes = [0 0 -e/2-lmin ; 0 0 e/2+lmin];

% Cut-off function (50% full, 10% decrease)
cutoffx = vibsCutoff(1,L/5,L/10);
cutoffy = vibsCutoff(2,L/5,L/10);

% Green kernel function
Gxy         = @(X,Y) femGreenKernel(X,Y,'[exp(ikr)/r]',k0);
gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',k0);
gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',k0);
gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',k0);
G0          = '[1/r]';
gradyG0     = 'grady[1/r]';
cteGxy      =  1/(4*pi);
cteG0       =  1/(4*pi);

% Plane wave function
PW         = @(X) exp(1i*k0*X*X0') .* cutoffx(X) .* cutoffy(X);
gradxPW{1} = @(X) 1i*k0*X0(1) .* PW(X);
gradxPW{2} = @(X) 1i*k0*X0(2) .* PW(X);
gradxPW{3} = @(X) 1i*k0*X0(3) .* PW(X);

% Coupling coeff for Brackage-Werner simulation
beta = 1i*k0;

% Graphical representation
figure
plot(bound,'r')
hold on
plot(radiat,log(abs(PW(radiat.vtx))))
% plotNrm(bound,'w')
plot(msh(Xmes),'k')
axis equal
colorbar
view(45,10)

% Quadrature and finite elements (volumn)
omega = dom(mesh,4);
U     = fem(mesh,'P1');

% Quadrature and finite elements (boundary)
sigma = dom(bound,3);
u     = fem(bound,'P1');

% Left-hand side
tic
[A,B,C,D] = vibsHmxBlockOperator(omega,U,sigma,u,cL,cT,rhoS,c0,rho0,f,tol);
toc

% Add dirichlet condition to x and y unknows (penalization)
A(sub2ind(size(A),1:2*length(U),1:2*length(U))) = 1e15;

% Right-hand side
V    = cell(4,1);
V{1} = - integral(sigma,ntimes(U,1),PW);
V{2} = - integral(sigma,ntimes(U,2),PW);
V{3} = - integral(sigma,ntimes(U,3),PW);
V{4} = integral(sigma,ntimes(u),gradxPW);

% Resolution with Schur complement
Fa     = decomposition(A);
LHS    = @(V) D*V - C*(Fa \ (B.Ml*(B.Mr*V)) );
RHS    = V{end} - C*(Fa \ cell2mat(V(1:end-1)) );
mu     = mgcr(LHS,RHS,[],tol,100);
lambda = beta*mu;

% Measure of refexive and transmitted coeff
tic
Pmes = cteGxy .* integral(Xmes,sigma,Gxy,u)*lambda - ...
    cteGxy .* integral(Xmes,sigma,gradyGxy,ntimes(u))*mu;
Pmes(2) = Pmes(2) + PW(Xmes(2,:));
toc

% Finite element radiative operator --> \int_Sy G(x,y) psi(y) dy
tic
Srad = cteGxy .* integral(radiat.vtx,sigma,Gxy,u,tol);
Srad = Srad + cteG0 .* regularize(radiat.vtx,sigma,G0,u);

% Finite element radiative operator --> \int_Sy dnyG(x,y) psi(y) dy
Drad = cteGxy .* integral(radiat.vtx,sigma,gradyGxy,ntimes(u),tol);
Drad = Drad + cteG0 .* regularize(radiat.vtx,sigma,gradyG0,ntimes(u));
toc

% Domain solution
Psca = Srad*lambda - Drad*mu;
Pinc = PW(radiat.vtx);
Ptot = Psca + Pinc;

% Graphical repesentation
figure
subplot(1,2,1)
plot(bound)
hold on
plot(radiat,abs(Psca))
plot(msh(Xmes),abs(Pmes))
title('Scattered')
axis equal
colorbar
view(45,10)

subplot(1,2,2)
plot(bound)
hold on
plot(radiat,abs(Ptot))
plot(msh(Xmes),abs(Pmes))
title('Total')
axis equal
colorbar
view(45,10)

% Analytical solution
tic
c     = ones(length(f),1) * [c0 cL c0];
rho   = [rho0 rhoS rho0];
[R,T] = slabVibro(f,rho,c,e);
toc

% Comparison (db)
ref = 20*log10(abs([R ; T]));
sol = 20*log10(abs(Pmes));

norm(ref-sol)/norm(ref)

disp('~~> Michto gypsilab !')


text:
clf
set(gcf, 'Position', [90   410   651   545]);
goal = [4 8];
start = [7 2];
opt.metric = 'euclidean'
opt.show = 1
opt.movie = [];%'dxform2.mp4'

% make a simple map
occgrid = zeros(10,10);
occgrid(4:6,3:7) = 1;
%occgrid(7:8,7) = 1;  % extra bit

cost0 = occgrid;
cost0(cost0==1) = NaN;

cost = cost0;
cost(cost0==0) = Inf;
cost(goal(2), goal(1)) = 0;

if ~isempty(opt.movie)
    anim = Animate(opt.movie);
end

showpixels(cost, 'contrast', 6, 'fmt', '%.2g', 'cscale', [0 12], 'fontsize', 20, 'infsymbol', 'nancolor', 'nohidenan', 'nohideinf', 'infcolor')

if ~isempty(opt.movie)
    anim.add();
end

switch opt.metric
    case 'cityblock'
        m = [inf 1 inf
              1  0  1
             inf 1 inf];
    case 'euclidean'
        r2 = sqrt(2);
        m = [r2 1 r2
              1 0  1
             r2 1 r2];
    otherwise
        error('unknown distance metric');
end

count = 0;
ninf = 0;

while 1

    cost = imorph(cost, m, 'plusmin');
    count = count+1;
    if opt.show
        % transfer over the finite values
%         k = ~isinf(cost0);
%         cost(k) = cost0(k);

        showpixels(cost, 'contrast', 6, 'fmt', '%.2g', 'cscale', [0 14], 'fontsize', 20, 'nancolor', 'nohideinf', 'infsymbol', 'nohidenan', 'infcolor', 'here')
        xlabel('x', 'FontSize', 20); ylabel('y', 'FontSize', 20)
        %         cmap = [1 0 0; gray(count)];
        %         colormap(cmap)
        %         image(occgrid+1, 'CDataMapping', 'direct');
        %         set(gca, 'Ydir', 'normal');
        %         xlabel('x');
        %         ylabel('y');
        if ~isempty(opt.movie)
            anim.add();
        end
        pause(opt.show);
    end

    ninfnow = sum( isinf(cost(:)) ); % current number of Infs
    if ninfnow == ninf
        % stop if the number of Infs left in the map had stopped reducing
        % it may never get to zero if there are unreachable cells in the map
        break;
    end
    ninf = ninfnow;
end
count
if ~isempty(opt.movie)
    anim.close();
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [fBValue, fAValue, fStdDev, fMc] = plot_FMD2(mCatalog, hAxes, rCalculateMC, rDisplay)
% function [fBValue, fAValue, fStdDev, fMc]
%   = plot_FMD2(mCatalog, bCumulative, hAxes, sSymbol, sColor, bPlotB, nCalculateMC, fBinning)
% --------------------------------------------------------------------------------------------
% Creates and plots a frequency magnitude distribution including the b-value
%
% plot_FMD2(mCatalog) opens a figure and plots the frequency magnitude distribution
%   with standard parameters
%
% Input parameters:
%   mCatalog        Earthquake catalog
%   hAxes           Handle of axes to plot the frequency magnitude distribution
%   rCalculateMC    Record with information about magnitude of completeness computing
%     nMethod       Method to determine the magnitude of completeness (see calc_Mc) (default 1)
%     bConstrain    Constrain Mc to a given range [fMinMc; fMaxMc] (default 0)
%     fMinMc        Minimum Mc when constrained
%     fMaxMc        Maximum Mc when constrained
%     fBinning      Magnitude binning of the catalog (default 0.1)
%   rDisplay        Record with display preferences
%     sSymbol       Symbol for FMD (default 's')
%     sColor        Color of FMD (default 'k')
%     bPlotB        Plot b-value line (=1) or not (=0) (default 1)
%     bCumulative   Plot cumulative frequency magnitude distribution (=1) or non-cumulative (=0)
%
% Output parameters:
%   fBValue         Calculated b-value
%   fAValue         Calculateda-value
%   fStdDev         Standard deviation of b-value
%   fMc             Magnitude of completeness
%
% Danijel Schorlemmer
% November 5, 2003

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Create the frequency magnitude distribution vector
[vFMD, vNonCFMD] = calc_FMD(mCatalog);

% Define missing input parameters
if ~exist('hAxes', 'var')
  figure;
  hAxes = newplot;
end
if ~exist('rCalculateMC', 'var')
  rCalculateMC.nMethod = 1;
  rCalculateMC.bConstrain = 0;
  rCalculateMC.fMinMc = -10;
  rCalculateMC.fMaxMc = 10;
  rCalculateMC.fBinning = 0.1;
end
if ~exist('rDisplay', 'var')
  rDisplay.sSymbol = 's';
  rDisplay.sColor = 'k';
  rDisplay.bPlotB = 1;
  rDisplay.bCumulative = 1;
end

% Activate given axes
axes(hAxes);
set(hAxes, 'NextPlot', 'add');

% Plot the frequency magnitude distribution
if rDisplay.bCumulative
  hPlot = semilogy(vFMD(1,:), vFMD(2,:), [rDisplay.sSymbol rDisplay.sColor]);
else
  hPlot = semilogy(vNonCFMD(1,:), vNonCFMD(2,:), [rDisplay.sSymbol rDisplay.sColor]);
end
set(hAxes, 'YScale', 'log'); % Bug in Matlab

if rDisplay.bPlotB
  % Calculate magnitude of completeness
  [fBValue, vDummy, fMc, fAValue] = calc_BandMc(mCatalog, 0, rCalculateMC.nMethod, rCalculateMC.fBinning, rCalculateMC.bConstrain, rCalculateMC.fMinMc, rCalculateMC.fMaxMc);

  % Determine the positions of 'x'-markers
  nIndexLo = find((vFMD(1,:) < fMc + (rCalculateMC.fBinning/2)) & (vFMD(1,:) > fMc - (rCalculateMC.fBinning/2)));
  fMagHi = vFMD(1,1);
  vSel = vFMD(1,:) <= fMagHi & vFMD(1,:) >= fMc-.0001;
  vMagnitudes = vFMD(1,vSel);

  % Plot the 'x'-marker
  hPlot = semilogy(vFMD(1,nIndexLo), vFMD(2,nIndexLo), ['x' rDisplay.sColor]);
  set(hPlot, 'LineWidth', [2], 'MarkerSize', 12);
  hPlot = semilogy(vFMD(1,1), vFMD(2,1), ['x' rDisplay.sColor]);
  set(hPlot, 'LineWidth', [2], 'MarkerSize', 12)

  % Plot the line representing the b-value
  vPoly = [-1*fBValue fAValue];
  fBFunc = 10.^(polyval(vPoly, vMagnitudes));
  hPlot = semilogy(vMagnitudes, fBFunc, rDisplay.sColor);
  set(hPlot, 'LineWidth', [2]);
end

author: CelsoReyes


text:
function loglik = RicianLogLik(meas, signals, sig)
% Computes the log likelihood of the measurements given the model signals
% for the Rician noise model.
%
% loglik = RicianLogLik(meas, signals, sig) returns the likelihood of
% measuring meas given the signals and the noise standard deviation sig.
%
% meas are the measurements
%
% signals are computed from a model
%
% sig is the standard deviation of the Gaussian distributions underlying
% the Rician distribution.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

sumsqsc = (signals.^2 + meas.^2)./(2*sig.^2);
scp = meas.*signals./(sig.^2);
lb0 = logbesseli0(scp);
logliks = - 2*log(sig) - sumsqsc + log(signals) + lb0;
loglik = sum(logliks);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function allen_ccf_2pi(tv,av,st)
% allen_ccf_2pi(tv,av,st)
% written by Samuel Picard (samuel.picard@ucl.ac.uk)
% based on original allen_ccf_npx tool written by Andy Peters (peters.andrew.j@gmail.com)
%
% GUI for planning 2pi chronic window implant with the Allen CCF
% Part of repository: https://github.com/cortex-lab/allenCCF
% - directions for installing atlas in repository readme
% - some dependent functions from that repository
%
% (optional inputs - if CCF path written in Line 22, loaded automatically)
% tv, av, st = CCF template volume, annotated volume, structure tree
% tv = readNPY('template_volume_10um.npy');
% av = readNPY('annotation_volume_10um_by_index.npy');
% st = loadStructureTree('structure_tree_safe_2017.csv');

% Initialize gui_data structure
gui_data = struct;

% Allen CCF-bregma transform (estimated from eyeballing Paxinos->CCF)
% [AP,DV,ML]
bregma = [540,0,570];

% If not already loaded in, load in atlas
if nargin < 3
    allen_atlas_path = 'C:\Users\Samuel\Documents\GitHub\allenCCF'; %put path to atlas files here
    if isempty(allen_atlas_path)
        error('Enter path where Allen CCF is stored at Line 26');
    end
    tv = readNPY([allen_atlas_path filesep 'template_volume_10um.npy']); % grey-scale "background signal intensity"
    av = readNPY([allen_atlas_path filesep 'annotation_volume_10um_by_index.npy']); % the number at each pixel labels the area, see note below
    st = loadStructureTree([allen_atlas_path filesep 'structure_tree_safe_2017.csv']); % a table of what all the labels mean
end

% Load the colormap (located in the repository, find by associated fcn)
allenCCF_path = fileparts(which('allenCCFbregma'));
cmap_filename = [allenCCF_path filesep 'allen_ccf_colormap_2017.mat'];
load(cmap_filename);

% Set up the gui
window_atlas_gui = figure('Toolbar','none','Menubar','none','color','w', ...
    'Name','Atlas-window viewer','Units','normalized','Position',[0.2,0.2,0.7,0.7]);

% Set up the atlas axes
axes_atlas = subplot(2,2,[1 3]);
[~, brain_outline] = plotBrainGrid([],axes_atlas);
hold(axes_atlas,'on');
axis vis3d equal off manual
view([-30,25]);
caxis([0 300]);
[ap_max,dv_max,ml_max] = size(tv);
xlim([-10,ap_max+10])
ylim([-10,ml_max+10])
zlim([-10,dv_max+10])

% Set up the window area axes
%top of window
axes_window_areas = subplot(2,2,2);
axes_window_areas.ActivePositionProperty = 'position';
window_areas_plot = image(nan(166,166)); %initialized for 5mm window at 30um resolution
axis off auto
window_areas_text = text(axes_window_areas,zeros(1,100),zeros(1,100),repmat({''}, 1, 100));
set(axes_window_areas,'FontSize',11,'XLimSpec', 'Tight','YLimSpec', 'Tight');
colormap(axes_window_areas,cmap);
caxis([1,size(cmap,1)])
%bottom of window
axes_window_areas_bottom = subplot(2,2,4);
axes_window_areas_bottom.ActivePositionProperty = 'position';
window_areas_bottom_plot = image(nan(166,166)); %initialized for 5mm window at 30um resolution
axis off auto
window_areas_bottom_text = text(axes_window_areas_bottom,zeros(1,100),zeros(1,100),repmat({''}, 1, 100));
set(axes_window_areas_bottom,'FontSize',11,'XLimSpec', 'Tight','YLimSpec', 'Tight');
colormap(axes_window_areas_bottom,cmap);
caxis([1,size(cmap,1)])

% Position the axes
set(axes_atlas,'Position',[-0.15,-0.1,1,1.2]);
set(axes_window_areas,'Position',[0.7,0.55,0.2,0.34]);
set(axes_window_areas_bottom,'Position',[0.7,0.1,0.2,0.34]);

% Set the current axes to the atlas (dirty, but some gca requirements)
axes(axes_atlas);

% Set up the window reference/actual
window_centre = [bregma(1),bregma(3),0];
window_depth = 60; % effective z-range of 2p microscope (in 10 ums)
window_diameter = 500; % diameter in 10 ums
window_vector = [window_centre',[window_centre(1),window_centre(2),window_depth]'];
window_vector_line = line(window_vector(1,:),window_vector(2,:),window_vector(3,:),'linewidth',3,'color','b','linestyle','-');
window_rads = linspace(-pi,pi,100);
window_circle = line(0.5*window_diameter*cos(window_rads)+window_centre(1),...
    0.5*window_diameter*sin(window_rads)+window_centre(2),...
    zeros(size(window_rads))+window_centre(3),...
    'linewidth',1.5,'color','r','linestyle','--');

% add titles for window plots
title(axes_window_areas,'Top of window');
title(axes_window_areas_bottom,[sprintf('Bottom (%i ',10*window_depth),'\mu','m)']);

% Set up the text to display coordinates
window_coordinates_text = uicontrol('Style','text','String','', ...
    'Units','normalized','Position',[0,0.95,1,0.05], ...
    'BackgroundColor','w','HorizontalAlignment','left','FontSize',12);

% Store data
gui_data.tv = tv; % Intensity atlas
gui_data.av = av; % Annotated atlas
gui_data.st = st; % Labels table
gui_data.cmap = cmap; % Atlas colormap
gui_data.bregma = bregma; % Bregma for external referencing
gui_data.window_depth = window_depth; % Effective z-depth under window
gui_data.structure_plot_idx = []; % Plotted structures
gui_data.window_angle = [0;0]; % window angles in ML/DV
gui_data.window_centre = window_centre; % Window reference centre on 3D atlas
gui_data.window_diameter = window_diameter;
gui_data.ref_centre = window_centre;

%Store handles
gui_data.handles.cortex_outline = brain_outline;
gui_data.handles.structure_patch = []; % Plotted structures
gui_data.handles.axes_atlas = axes_atlas; % Axes with 3D atlas
gui_data.handles.axes_window_areas = axes_window_areas; % Axes with window areas
gui_data.handles.axes_window_areas_bottom = axes_window_areas_bottom; % Axes with window areas
gui_data.handles.slice_plot = surface('EdgeColor','none'); % Slice on 3D atlas
gui_data.handles.slice_volume = 'tv'; % The volume shown in the slice
gui_data.handles.window_vector = window_vector_line; % Window centre vector on 3D atlas
gui_data.handles.window_circle = window_circle; % Window circle on 3D atlas
gui_data.handles.window_areas_plot = window_areas_plot; % Color-coded window regions
gui_data.handles.window_areas_bottom_plot = window_areas_bottom_plot; % Color-coded window regions
gui_data.handles.window_areas_text = window_areas_text; % Labels for window regions
gui_data.handles.window_areas_bottom_text = window_areas_bottom_text; % Labels for window regions
gui_data.window_coordinates_text = window_coordinates_text; % Window coordinates text

% Make 3D rotation the default state (toggle on/off with 'r')
h = rotate3d(axes_atlas);
h.Enable = 'on';
% Update the slice whenever a rotation is completed
%h.ActionPostCallback = @update_slice;

% Set functions for key presses
hManager = uigetmodemanager(window_atlas_gui);
[hManager.WindowListenerHandles.Enabled] = deal(false);
set(window_atlas_gui,'KeyPressFcn',@key_press);
set(window_atlas_gui,'KeyReleaseFcn',@key_release);

% Upload gui_data
guidata(window_atlas_gui, gui_data);

% Display the first slice and update the window position
update_slice(window_atlas_gui);
update_window_coordinates(window_atlas_gui);

% Display controls
display_controls;

end

function key_press(window_atlas_gui,eventdata)

% Get guidata
gui_data = guidata(window_atlas_gui);

switch eventdata.Key

    case 'uparrow'
        if isempty(eventdata.Modifier)
            % Up: move window anterior
            ap_offset = -10;
            set(gui_data.handles.window_circle,'XData',get(gui_data.handles.window_circle,'XData') + ap_offset);
            set(gui_data.handles.window_vector,'XData',get(gui_data.handles.window_vector,'XData') + ap_offset);
            gui_data.window_centre = gui_data.window_centre + [ap_offset,0,0];
        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-up: increase DV angle
            angle_change = [1;0];
            gui_data = update_window_angle(window_atlas_gui,angle_change);
        elseif any(strcmp(eventdata.Modifier,'alt'))
            % Alt-up: raise window
            dv_offset = -10;
            old_window_vector = cell2mat(get(gui_data.handles.window_vector,{'XData','YData','ZData'})');
            old_window_circle = cell2mat(get(gui_data.handles.window_circle,{'XData','YData','ZData'})');
            move_vector = diff(old_window_vector,[],2)./ ...
                norm(diff(old_window_vector,[],2))*dv_offset;
            new_window_vector = bsxfun(@plus,old_window_vector,move_vector);
            new_window_circle = bsxfun(@plus,old_window_circle,move_vector);
            set(gui_data.handles.window_vector,'XData',new_window_vector(1,:), ...
                'YData',new_window_vector(2,:),'ZData',new_window_vector(3,:));
            set(gui_data.handles.window_circle,'XData',new_window_circle(1,:), ...
                'YData',new_window_circle(2,:),'ZData',new_window_circle(3,:));
            gui_data.window_centre = gui_data.window_centre + move_vector';

        end

    case 'downarrow'
        if isempty(eventdata.Modifier)
            % Down: move window posterior
            ap_offset = 10;
            set(gui_data.handles.window_circle,'XData',get(gui_data.handles.window_circle,'XData') + ap_offset);
            set(gui_data.handles.window_vector,'XData',get(gui_data.handles.window_vector,'XData') + ap_offset);
            gui_data.window_centre = gui_data.window_centre + [ap_offset,0,0];
        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-down: decrease DV angle
            angle_change = [-1;0];
            gui_data = update_window_angle(window_atlas_gui,angle_change);
        elseif any(strcmp(eventdata.Modifier,'alt'))
            % Alt-down: lower window
            dv_offset = 10;
            old_window_vector = cell2mat(get(gui_data.handles.window_vector,{'XData','YData','ZData'})');
            old_window_circle = cell2mat(get(gui_data.handles.window_circle,{'XData','YData','ZData'})');
            move_vector = diff(old_window_vector,[],2)./ ...
                norm(diff(old_window_vector,[],2))*dv_offset;
            new_window_vector = bsxfun(@plus,old_window_vector,move_vector);
            new_window_circle = bsxfun(@plus,old_window_circle,move_vector);
            set(gui_data.handles.window_vector,'XData',new_window_vector(1,:), ...
                'YData',new_window_vector(2,:),'ZData',new_window_vector(3,:));
            set(gui_data.handles.window_circle,'XData',new_window_circle(1,:), ...
                'YData',new_window_circle(2,:),'ZData',new_window_circle(3,:));
            gui_data.window_centre = gui_data.window_centre + move_vector';

        end

    case 'rightarrow'
        if isempty(eventdata.Modifier)
            % Right: move window right
            ml_offset = 10;
            set(gui_data.handles.window_circle,'YData',get(gui_data.handles.window_circle,'YData') + ml_offset);
            set(gui_data.handles.window_vector,'YData',get(gui_data.handles.window_vector,'YData') + ml_offset);
            gui_data.window_centre = gui_data.window_centre + [0,ml_offset,0];
        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-right: increase vertical angle
            angle_change = [0;1];
            gui_data = update_window_angle(window_atlas_gui,angle_change);
        end

    case 'leftarrow'
        if isempty(eventdata.Modifier)
            % Left: move window left
            ml_offset = -10;
            set(gui_data.handles.window_circle,'YData',get(gui_data.handles.window_circle,'YData') + ml_offset);
            set(gui_data.handles.window_vector,'YData',get(gui_data.handles.window_vector,'YData') + ml_offset);
            gui_data.window_centre = gui_data.window_centre + [0,ml_offset,0];
        elseif any(strcmp(eventdata.Modifier,'shift'))
            % Ctrl-left: decrease vertical angle
            angle_change = [0;-1];
            gui_data = update_window_angle(window_atlas_gui,angle_change);
        end

    case 'c'
        % Bring up controls again
        display_controls;

    case 'b'
        % Toggle brain outline visibility
        current_visibility = gui_data.handles.cortex_outline.Visible;
        switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
        set(gui_data.handles.cortex_outline,'Visible',new_visibility);

    case 'a'
        % Toggle plotted structure visibility
        if ~isempty(gui_data.structure_plot_idx)
            current_alpha = get(gui_data.handles.structure_patch(1),'FaceAlpha');
            switch current_alpha
                case 0
                    new_alpha = 0.2;
                case 0.2
                    new_alpha = 1;
                case 1
                    new_alpha = 0;
            end
            set(gui_data.handles.structure_patch,'FaceAlpha',new_alpha);
        end

    case 's'
        % Toggle slice volume/visibility
        slice_volumes = {'tv','av','none'};
        new_slice_volume = slice_volumes{circshift( ...
            strcmp(gui_data.handles.slice_volume,slice_volumes),[0,1])};

        if strcmp(new_slice_volume,'none')
            set(gui_data.handles.slice_plot,'Visible','off');
        else
            set(gui_data.handles.slice_plot,'Visible','on');
        end

        gui_data.handles.slice_volume = new_slice_volume;
        guidata(window_atlas_gui, gui_data);

        update_slice(window_atlas_gui);

    case 'w'
        % Toggle window visibility
        current_visibility = gui_data.handles.window_circle.Visible;
        switch current_visibility; case 'on'; new_visibility = 'off'; case 'off'; new_visibility = 'on'; end;
        set(gui_data.handles.window_circle,'Visible',new_visibility);

    case 'm'
        % Set window position manually and find window angle automatically
        set_window_position(window_atlas_gui);
        % Get updated guidata
        gui_data = guidata(window_atlas_gui);

    case {'equal','add'}
        % Add structure(s) to display
        slice_spacing = 10;

        % Prompt for which structures to show (only structures which are
        % labelled in the slice-spacing downsampled annotated volume)

        if ~any(strcmp(eventdata.Modifier,'shift'))
            % (no shift: list in native CCF order)

            parsed_structures = unique(reshape(gui_data.av(1:slice_spacing:end, ...
                1:slice_spacing:end,1:slice_spacing:end),[],1));
            if ~any(strcmp(eventdata.Modifier,'alt'))
                % (no alt: list all)
                plot_structures_parsed = listdlg('PromptString','Select a structure to plot:', ...
                    'ListString',gui_data.st.safe_name(parsed_structures),'ListSize',[520,500]);
                plot_structures = parsed_structures(plot_structures_parsed);
            else
                % (alt: search list)
                structure_search = lower(inputdlg('Search structures'));
                structure_match = find(contains(lower(gui_data.st.safe_name),structure_search));
                list_structures = intersect(parsed_structures,structure_match);
                if isempty(list_structures)
                    error('No structure search results')
                end

                plot_structures_parsed = listdlg('PromptString','Select a structure to plot:', ...
                    'ListString',gui_data.st.safe_name(list_structures),'ListSize',[520,500]);
                plot_structures = list_structures(plot_structures_parsed);
            end

            if ~isempty(plot_structures)
                for curr_plot_structure = reshape(plot_structures,1,[])
                    % If this label isn't used, don't plot
                    if ~any(reshape(gui_data.av( ...
                            1:slice_spacing:end,1:slice_spacing:end,1:slice_spacing:end),[],1) == curr_plot_structure)
                        disp(['"' gui_data.st.safe_name{curr_plot_structure} '" is not parsed in the atlas'])
                        continue
                    end

                    gui_data.structure_plot_idx(end+1) = curr_plot_structure;

                    plot_structure_color = hex2dec(reshape(gui_data.st.color_hex_triplet{curr_plot_structure},2,[])')./255;
                    structure_3d = isosurface(permute(gui_data.av(1:slice_spacing:end, ...
                        1:slice_spacing:end,1:slice_spacing:end) == curr_plot_structure,[3,1,2]),0);

                    if isempty(gui_data.handles.structure_patch)
                        structure_alpha = 0.2;
                    else
                        structure_alpha = get(gui_data.handles.structure_patch(1),'FaceAlpha');
                    end
                    gui_data.handles.structure_patch(end+1) = patch('Vertices',structure_3d.vertices*slice_spacing, ...
                        'Faces',structure_3d.faces, ...
                        'FaceColor',plot_structure_color,'EdgeColor','none','FaceAlpha',structure_alpha);
                end
            end

        elseif any(strcmp(eventdata.Modifier,'shift'))
            % (shift: use hierarchy search)
            if ~any(strcmp(eventdata.Modifier,'control'))
                % (no ctrl, choose between all structures)
                plot_structures = hierarchicalSelect(gui_data.st);
            else
                % (ctrl: take pre-defined regions: isocortex, cerebellum, SC & IC)
                plot_structures = [7,13,31,108,115,122,158,221,239,246,253,279,298,340,361,368,374,381,809,813,1016,1017];
            end

            if ~isempty(plot_structures) % will be empty if dialog was cancelled
                % get all children of this one
                for plot_structure = plot_structures
                    thisID = gui_data.st.id(plot_structure);
                    idStr = sprintf('/%d/', thisID);
                    theseCh = find(cellfun(@(x)contains(x,idStr), gui_data.st.structure_id_path));

                    % plot the structure
                    slice_spacing = 5;
                    plot_structure_color = hex2dec(reshape(gui_data.st.color_hex_triplet{plot_structure},3,[]))./255;
                    structure_3d = isosurface(permute(ismember(gui_data.av(1:slice_spacing:end, ...
                        1:slice_spacing:end,1:slice_spacing:end),theseCh),[3,1,2]),0);

                    if isempty(gui_data.handles.structure_patch)
                        structure_alpha = 0.2;
                    else
                        structure_alpha = get(gui_data.handles.structure_patch(1),'FaceAlpha');
                    end
                    gui_data.structure_plot_idx(end+1) = plot_structure;
                    gui_data.handles.structure_patch(end+1) = patch('Vertices',structure_3d.vertices*slice_spacing, ...
                        'Faces',structure_3d.faces, ...
                        'FaceColor',plot_structure_color,'EdgeColor','none','FaceAlpha',structure_alpha);

                end

            end

        end

    case {'hyphen','subtract'}
        % Remove structure(s) already plotted
        if ~isempty(gui_data.structure_plot_idx)
            remove_structures = listdlg('PromptString','Select a structure to remove:', ...
                'ListString',gui_data.st.safe_name(gui_data.structure_plot_idx));
            delete(gui_data.handles.structure_patch(remove_structures))
            gui_data.structure_plot_idx(remove_structures) = [];
            gui_data.handles.structure_patch(remove_structures) = [];
        end

    case 'x'
        % Export the window vector coordinates in Allen CCF to the workspace
        window_vector = cell2mat(get(gui_data.handles.window_vector_line,{'XData','YData','ZData'})');
        window_vector_ccf = round(window_vector([1,3,2],:))';
        assignin('base','window_vector_ccf',window_vector_ccf)
        disp('Copied window vector coordinates to workspace');

end

% Upload gui_data
guidata(window_atlas_gui, gui_data);

end

function key_release(window_atlas_gui,eventdata)

% Get guidata
gui_data = guidata(window_atlas_gui);

switch eventdata.Key
    case {'rightarrow','leftarrow','uparrow','downarrow'}
        % Update the window info/slice on arrow release
        update_window_coordinates(window_atlas_gui);
        update_slice(window_atlas_gui);
end

% Upload gui_data
guidata(window_atlas_gui, gui_data);

end

function update_slice(window_atlas_gui,varargin)

% Get guidata
gui_data = guidata(window_atlas_gui);

% Only update the slice if it's visible
if strcmp(gui_data.handles.slice_plot(1).Visible,'on')

    offsets = [0,gui_data.window_depth];
    plot_handles = {'window_areas','window_areas_bottom'};

    for islice = 1:length(offsets)

        % Get current window outline coordinate
        curr_window_circle = cell2mat(get(gui_data.handles.window_circle,{'XData','YData','ZData'})');
        window_vector = cell2mat(get(gui_data.handles.window_vector,{'XData','YData','ZData'})');
        move_vector = diff(window_vector,[],2)./ ...
            norm(diff(window_vector,[],2))*offsets(islice);
        window_circle = bsxfun(@plus,curr_window_circle,move_vector);
        window_centre = gui_data.window_centre + move_vector';

        % Get two vectors on the window plane
        vlength = size(window_circle,2);
        window_vector_1 = window_circle(:,1)' - window_circle(:,round(0.4*vlength))';
        window_vector_2 = window_circle(:,1)' - window_circle(:,round(0.7*vlength))';

        %get the normal vector of the plane
        normal_vector = cross(window_vector_1,window_vector_2);

        % Get the plane offset through the window centre
        window_top = window_centre;
        plane_offset = -(normal_vector*window_top');

        % Define a plane of points to index
        % (the plane grid is defined based on the which cardinal plan is most
        % orthogonal to the plotted plane. this is janky but it works)
        slice_px_space = 3;
        [~,cam_plane] = max(abs(normal_vector./norm(normal_vector)));

        switch cam_plane

            case 1
                [plane_y,plane_z] = meshgrid(1:slice_px_space:size(gui_data.tv,3),1:slice_px_space:size(gui_data.tv,2));
                plane_x = ...
                    (normal_vector(2)*plane_y+normal_vector(3)*plane_z + plane_offset)/ ...
                    -normal_vector(1);

            case 2
                [plane_x,plane_z] = meshgrid(1:slice_px_space:size(gui_data.tv,1),1:slice_px_space:size(gui_data.tv,2));
                plane_y = ...
                    (normal_vector(1)*plane_x+normal_vector(3)*plane_z + plane_offset)/ ...
                    -normal_vector(2);

            case 3
                [plane_x,plane_y] = meshgrid(1:slice_px_space:size(gui_data.tv,1),1:slice_px_space:size(gui_data.tv,3));
                plane_z = ...
                    (normal_vector(1)*plane_x+normal_vector(2)*plane_y + plane_offset)/ ...
                    -normal_vector(3);

        end

        % Get the coordinates on the plane
        x_idx = round(plane_x);
        y_idx = round(plane_y);
        z_idx = round(plane_z);

        % make a mask for the window (in the horizontal plane)
        window_mask = uint8(((plane_x-window_centre(1)).^2 + (plane_y-window_centre(2)).^2 + (plane_z-window_centre(3)).^2)<(gui_data.window_diameter/2)^2);

        % Find plane coordinates in bounds with the volume and the window
        use_xd = x_idx > 0 & x_idx < size(gui_data.tv,1);
        use_yd = y_idx > 0 & y_idx < size(gui_data.tv,3);
        use_zd = z_idx > 0 & z_idx < size(gui_data.tv,2);
        use_idx = use_xd & use_yd & use_zd & window_mask;

        curr_slice_idx = sub2ind(size(gui_data.tv),x_idx(use_idx),z_idx(use_idx),y_idx(use_idx));

        % Find plane coordinates that contain brain
        curr_slice_isbrain = false(size(use_idx));
        curr_slice_isbrain(use_idx) = gui_data.av(curr_slice_idx) > 1;

        % Index coordinates in bounds + with brain
        grab_pix_idx = sub2ind(size(gui_data.tv),x_idx(curr_slice_isbrain),z_idx(curr_slice_isbrain),y_idx(curr_slice_isbrain));

        % Grab pixels from (selected) volume
        curr_slice = nan(size(use_idx));
        switch gui_data.handles.slice_volume
            case 'tv'
                curr_slice(curr_slice_isbrain) = gui_data.tv(grab_pix_idx);
                colormap(gui_data.handles.axes_atlas,'gray');
                caxis([0,255]);
            case 'av'
                curr_slice(curr_slice_isbrain) = gui_data.av(grab_pix_idx);
                colormap(gui_data.handles.axes_atlas,gui_data.cmap);
                caxis([1,size(gui_data.cmap,1)]);
        end

        % Update the slice display
        if islice==2 %1 = top, 2 = bottom
            set(gui_data.handles.slice_plot,'XData',plane_x,'YData',plane_y,'ZData',plane_z,'CData',curr_slice);
        end

        % Update the title with new depth
        if islice==2
            set(eval(['gui_data.handles.axes_',plot_handles{islice},'.Title']),'String',[sprintf('Bottom (%i ',10*offsets(islice)),'\mu','m)']);
        end

        % Update the window overview
        curr_slice_av = nan(size(use_idx));
        curr_slice_av(curr_slice_isbrain) = gui_data.av(grab_pix_idx);
        x_edges = [find(sum(window_mask,1),1,'first'),find(sum(window_mask,1),1,'last')];
        y_edges = [find(sum(window_mask,2),1,'first'),find(sum(window_mask,2),1,'last')];
        cropped_window = curr_slice_av(y_edges(1):y_edges(2),x_edges(1):x_edges(2))';
        set(eval(['gui_data.handles.',plot_handles{islice},'_plot']),'CData',cropped_window);

        %update area labels based on which areas are contained in window
        window_areas = cropped_window;
        unique_areas = unique(window_areas(~isnan(window_areas(:))));
        unique_parent_struct = unique(gui_data.st.parent_structure_id(unique_areas));
        for i_struct = 1:100
            if i_struct<=length(unique_parent_struct)
                areas_in_parent_struct = find(gui_data.st.parent_structure_id==unique_parent_struct(i_struct));
                area_bm = ismember(window_areas,areas_in_parent_struct);
                area_bm_largest = zeros(size(area_bm));
                cc = bwconncomp(area_bm);
                [~,largest_blob] = max(cellfun('length',cc.PixelIdxList));
                area_bm_largest(cc.PixelIdxList{largest_blob}) = 1;
                [area_x,area_y] = find(area_bm_largest==1);
                area_centre = [mean(area_y),mean(area_x)];
                area_string = gui_data.st.acronym{find(gui_data.st.id==unique_parent_struct(i_struct))};
            else
                area_centre = [0,0];
                area_string = '';
            end
            set(eval(['gui_data.handles.',plot_handles{islice},'_text(',num2str(i_struct),')']),...
                'Position',area_centre,'String',area_string,...
                'VerticalAlignment','middle','HorizontalAlignment','center');
        end
    end

    % Upload gui_data
    guidata(window_atlas_gui, gui_data);

end

end

function set_window_position(window_atlas_gui,varargin) %this function should place window at tangent plane

% Get guidata
gui_data = guidata(window_atlas_gui);

prompt_text = { ...
    'AP position (mm from bregma)', ...
    'ML position (mm from bregma)', ...
    'window diameter (mm)', ...
    'maximum depth (microns)'};

new_window_position = cellfun(@str2num,inputdlg(prompt_text,'Set window position',1));

% Convert centre position: mm->CCF
window_centre_ccf_coords = round(gui_data.bregma([1,3])' - 100*[new_window_position(1);-new_window_position(2)]);%not sure why we need this sign change...
window_diameter = new_window_position(3)*100;
window_depth = new_window_position(4)/10;

%load surface data or generate it from scratch
if ~exist('surf_coords','var')
    try
        load('dorsalsurface_nomidline_ccf_coords.mat')
    catch
        %CODE TO GENERATE SURFACE MAP
        surf_coords = nan(size(gui_data.av,1),size(gui_data.av,3));
        for i=1:size(gui_data.av,1)
            for j=1:size(gui_data.av,3)
                idxs = find(squeeze(gui_data.av(i,:,j))>1,1,'first');
                if ~isempty(idxs)
                    surf_coords(i,j) = idxs;
                end
            end
            idx_c = round(size(gui_data.av,3)/2);
            [m,idx_m] = min(surf_coords(i,:));
            idx_m_c = idx_m-idx_c;
            surf_coords(i,idx_c-abs(idx_m_c):idx_c+abs(idx_m_c)) = m;
        end
        save('dorsalsurface_nomidline_ccf_coords.mat',surf_coords);
    end
end

%figure out angle of tangent plane
surf_depth_coords = surf_ccf_coords(window_centre_ccf_coords(1),window_centre_ccf_coords(2));
centre_surf_ccf_coords = [window_centre_ccf_coords(1); surf_depth_coords; window_centre_ccf_coords(2)];
centre_patch = surf_ccf_coords(window_centre_ccf_coords(1)-50:window_centre_ccf_coords(1)+50,window_centre_ccf_coords(2)-50:window_centre_ccf_coords(2)+50);
centre_patch_smooth = imgaussfilt(centre_patch,'FilterSize',9);
[fx,fy] = gradient(centre_patch_smooth,1);
avg_gradient = [-nanmean(fy(:));nanmean(fx(:))];
window_angle_rad = atan(avg_gradient);
window_angle_deg = (window_angle_rad/(2*pi))*360;

%compute new window vector and circle (horizontal)
%window_diameter = gui_data.window_diameter;
window_rads = linspace(-pi,pi,100);
window_centre = [window_centre_ccf_coords', surf_depth_coords-gui_data.ref_centre(3)];
window_circle = [0.5*window_diameter*cos(window_rads)+window_centre(1);...
    0.5*window_diameter*sin(window_rads)+window_centre(2);...
    zeros(size(window_rads))+window_centre(3)];
window_vector = [window_centre',(window_centre+[0,0,window_depth])'];

%rotate these by the tangent angle
eul = [0;window_angle_rad]; %[yaw; pitch; roll]
R = eul2rotm(eul'); %find rotation matrix
window_circle_centred = window_circle - window_centre'; %translate to origin by subtracting its centre coords
window_vector_centred = window_vector - window_centre';
new_window_circle_centred = R*window_circle_centred; %rotate around origin
new_window_vector_centred = R*window_vector_centred;
new_window_circle = new_window_circle_centred + window_centre'; %translate back to its position
new_window_vector = new_window_vector_centred + window_centre';

% update gui_data
set(gui_data.handles.window_circle,'XData',new_window_circle(1,:), ...
    'YData',new_window_circle(2,:), ...
    'ZData',new_window_circle(3,:));
set(gui_data.handles.window_vector,'XData',new_window_vector(1,:), ...
    'YData',new_window_vector(2,:), ...
    'ZData',new_window_vector(3,:));
gui_data.window_angle = window_angle_deg;
gui_data.window_centre = [window_centre_ccf_coords', surf_depth_coords];
gui_data.window_diameter = window_diameter;
gui_data.window_depth = window_depth;

% Upload gui_data
guidata(window_atlas_gui, gui_data);

% Update the slice and window coordinates
update_slice(window_atlas_gui);
update_window_coordinates(window_atlas_gui);

end

function gui_data = update_window_angle(window_atlas_gui,angle_change)

% Get guidata
gui_data = guidata(window_atlas_gui);

% Set new angle
new_angle = gui_data.window_angle + angle_change;
gui_data.window_angle = new_angle;

% Get the positions of the window circle and vector
window_circle = cell2mat(get(gui_data.handles.window_circle,{'XData','YData','ZData'})');
window_vector = cell2mat(get(gui_data.handles.window_vector,{'XData','YData','ZData'})');

% Compute the rotation matrix
eul = [0;(angle_change/360)*2*pi];
R = eul2rotm(eul');

% temporarily translate window position to origin
window_circle_centred = window_circle - gui_data.window_centre';
window_vector_centred = window_vector - gui_data.window_centre';

% rotate the vector around the origin using rotation matrix
new_window_circle_centred = R*window_circle_centred;
new_window_vector_centred = R*window_vector_centred;

% translate rotated window back to the original centre point
new_window_circle = new_window_circle_centred + gui_data.window_centre';
new_window_vector = new_window_vector_centred + gui_data.window_centre';

% update gui_data
set(gui_data.handles.window_circle,'XData',new_window_circle(1,:), ...
    'YData',new_window_circle(2,:), ...
    'ZData',new_window_circle(3,:));
set(gui_data.handles.window_vector,'XData',new_window_vector(1,:), ...
    'YData',new_window_vector(2,:), ...
    'ZData',new_window_vector(3,:));

% Upload gui_data
guidata(window_atlas_gui, gui_data);

end

function update_window_coordinates(window_atlas_gui,varargin)

% Get guidata
gui_data = guidata(window_atlas_gui);

window_bregma_coordinate = round((gui_data.bregma([1,3,2]) - gui_data.window_centre)*10);

% Update the text
window_text = ['Window position: ' ....
    num2str(window_bregma_coordinate(1)) ' AP, ', ...
    num2str(-window_bregma_coordinate(2)) ' ML, ', ...
    num2str(-window_bregma_coordinate(3)) ' DV, ', ...
    num2str(round(gui_data.window_angle(1))) char(176) ' pitch, ' ...
    num2str(round(gui_data.window_angle(2))) char(176) ' roll'];
set(gui_data.window_coordinates_text,'String',window_text);

% Upload gui_data
guidata(window_atlas_gui, gui_data);

end

function display_controls

% Print controls
CreateStruct.Interpreter = 'tex';
CreateStruct.WindowStyle = 'non-modal';
msgbox( ...
    {'\fontsize{12}' ...
    '\bf Window: \rm' ...
    'Arrow keys : translate window' ...
    'Alt/Option up/down : raise/lower window' ...
    'Shift arrow keys : change window angle' ...
    'm : set location manually (tangent to surface)', ...
    '\bf 3D brain areas: \rm' ...
    ' =/+ : add (list selector)' ...
    ' Alt/Option =/+ : add (search)' ...
    ' Shift =/+ : add (hierarchy selector)' ...
    ' Ctrl Shift =/+ : add all regions on surface' ...
    ' - : remove', ...
    '\bf Visibility: \rm' ...
    's : atlas slice (toggle tv/av/off)' ...
    'b : brain outline' ...
    'w : window outline' ...
    'a : 3D brain areas (toggle transparency)' ...
    '\bf Other: \rm' ...
    'x : export window coordinates to workspace' ...
    'c : bring up controls box'}, ...
    'Controls',CreateStruct);

end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
% This program obtains the analytical step response of a parallel RLC
% circuit. In addition, graphs of the voltage across the parallel
% combination and branch currents are obtained.
% H. Saadat, Copyright 1999

function PresponseIR
warning('off','MATLAB:dispatcher:InexactMatch')
Rhandle=findobj('Tag','Rtext');
R=eval(get(Rhandle,'String'));
KOhandle=findobj('Tag','KOhm');
KO=get(KOhandle,'Value');
if KO==1
   R=1000*R;
else, end
Lhandle=findobj('Tag','Ltext');
L=eval(get(Lhandle,'String'));
mHhandle=findobj('Tag','mHenry');
ml=get(mHhandle,'Value');
if ml==1
   L=L/1000;
else, end
Chandle=findobj('Tag','Ctext');
C=eval(get(Chandle,'String'));
micFhandle=findobj('Tag','micFarad');
micF=get(micFhandle,'Value');
if micF==1
   C=C/1000000;
else, end
Rhandle=findobj('Tag','IStext');
Is=eval(get(Rhandle,'String'));
Rhandle=findobj('Tag','I0text');
I0=eval(get(Rhandle,'String'));
Rhandle=findobj('Tag','V0text');
V0=eval(get(Rhandle,'String'));
% Step Response of a parallel RLC circuit
if L==inf | C==inf
   plot(0, 0)
   text(-0.8, 0.80, 'Parameters must lie in the following ranges:', 'color', [0 0 0.6275])
   text(-0.8, 0.65, '0 < R \leq inf,   0 < L < inf,   &    0 < C < inf', 'color', [0 0 0.6275])
   text(-0.8, 0.5, 'Enter the correct value and press Solve', 'color', [0 0 0.6275])
   axis([-1  1  -1  1]);return,  else
end
if R==0 | L==0 | C == 0
plot(0, 0)
title('Short-circuit across the parallel circuit', 'color', [ 1 0 0])
text(-0.8, 0.80, 'Parameters must lie in the following ranges:', 'color', [0 0 0.6275])
text(-0.8, 0.65, '0 < R \leq inf,   0 < L < inf,   &    0 < C < inf', 'color', [0 0 0.6275])
text(-0.8, 0.5, 'Enter the correct value and press Solve', 'color', [0 0 0.6275])
axis([-1  1  -1  1]);return,  else
alpha =1/(2*R*C);
end
w02 = 1/(L*C); w0 = sqrt(w02);
if alpha~= 0
   err=(alpha-w0)/alpha;
   if abs(err) <1e-4
   w0=alpha;
   else, end
else, end
dv = Is/C-I0/C-V0/(R*C);
di = V0/L;
if alpha > w0
 % Overdamped response
    s(1) = -alpha + sqrt(alpha^2 - w02);
    s(2) = -alpha - sqrt(alpha^2 - w02);
    A1 = (dv-s(2)*V0)/(s(1)-s(2));
    A2 = (dv-s(1)*V0)/(s(2)-s(1));
    A1i =(di-s(2)*(I0-Is))/(s(1)-s(2));
    A2i =(di-s(1)*(I0-Is))/(s(2)-s(1));
    tf = 6*max(abs(1/s(1)), abs(1/s(2)));
    t=0:tf/100:tf;
    v = A1*exp(s(1)*t)+A2*exp(s(2)*t);
    iL=Is + A1i*exp(s(1)*t)+A2i*exp(s(2)*t);
    iR = v/R;
    iC= Is -iR-iL;
    it=iR+iL+iC;
    plot(t,iR,'erasemode','none','color',[0.9 0 0.8]), grid
    title(['i_R(t) =(', num2str(A1/R), ') e^{(', num2str(s(1)), 't)} + (', num2str(A2/R), ') e^{(', num2str(s(2)), 't)}'],'color',[0.9 0 0.8])
    xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (O.D.)          t, sec'])
    ylabel('i_R(t), Amps')
    elseif alpha < w0
     if alpha~=0
      % Underdamped response
      tf = 6/alpha;
      t=0:tf/200:tf;
      wd= sqrt(w02 - alpha^2);
      s(1) = -alpha +j*wd;
      s(2) = -alpha -j*wd;
      B1 = V0; B2 = (dv+alpha*B1)/wd;
      B1i = I0 - Is; B2i = (di+alpha*B1i)/wd;
      v = exp(-alpha*t).*(B1*cos(wd*t)+B2*sin(wd*t));
      iL = Is + exp(-alpha*t).*(B1i*cos(wd*t)+B2i*sin(wd*t));
      iR = v/R;
      iC= Is -iR-iL;
      it=iR+iL+iC;
      plot(t,iR,'erasemode','none','color',[.9 0 0.8]), grid
      title(['i_R(t) = e^{(-', num2str(alpha), 't)}[(', num2str(B1/R), ')cos',num2str(wd),'t + (' , num2str(B2/R), ')sin',num2str(wd),'t]'],'color',[0.9 0 0.8])
      xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (U.D.)          t, sec'])
      ylabel('i_R(t), Amps')
      else
      % Undamped Response
      tf = 8*pi/w0;
      t=0:tf/200:tf;
      wd= w0;
      B1 = V0; B2 = (dv)/wd;
      B1i = I0 - Is; B2i = (di)/wd;
      iL = Is + B1i*cos(wd*t)+B2i*sin(wd*t);
      v = B1*cos(wd*t)+B2*sin(wd*t);
      iR = v/R;
      iC= Is -iR-iL;
      it=iR+iL+iC;
      plot(t,iR,'erasemode','none','color',[0.9 0 0.80]), grid
      title('i_R(t) =  0')
      xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (Undamped)          t, sec'])
      ylabel('i_R(t), Amps')
      end
   elseif alpha ==w0
    % Critically damped response
    s(1)=-alpha;  s(2) = s(1);
    D2 = V0;
    D1 = dv+alpha*D2;
    D2i = I0 - Is; D1i = di+alpha*D2i;
    tf = 8/alpha;
    t=0:tf/100:tf;
    iL=Is + exp(-alpha*t).*(D1i*t+D2i);
    v = exp(-alpha*t).*(D1*t+D2);
    iR = v/R;
    iC= Is -iR-iL;
    it=iR+iL+iC;
    plot(t,iR,'erasemode','none','color',[0.9 0 0.8]), grid
    title(['i_R(t) = e^{(-', num2str(alpha), 't)}[(', num2str(D1/R), ')t + (' , num2str(D2/R), ')]'],'color',[0.9 0 0.8])
    xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (C.D.)          t, sec'])
    ylabel('i_R(t), Amps')
    else, end

author: Sable

repo: mcbench-benchmarks


text:
function [alphahat epshat etahat] = fastsmo_int(n, y, mis, anymis, allmis, Hdyn, Zdyn, Tdyn, Rdyn, Qdyn, cdyn, Hmat, Zmat, Tmat, Rmat, Qmat, cmat, a1, P1, tol)

%% Kalman filter %%
[d Fns v invF K L L1 QRt] = kalman_int(5, n, y, mis, anymis, allmis, Hdyn, Zdyn, Tdyn, Rdyn, Qdyn, cdyn, Hmat, Zmat, Tmat, Rmat, Qmat, cmat, a1, P1, tol);

%% Initialization %%
D       = (P1 == Inf);
P1(D)   = 0;
P1_inf  = double(D);
if ~Hdyn, H = Hmat; end
if ~Zdyn, Z = Zmat; end
if ~Tdyn, T = Tmat; end
if ~Rdyn, R = Rmat; end
if ~cdyn, c = cmat; end

%% Disturbance smoothing backwards recursion %%
m   = size(a1, 1);
r   = zeros(m, 1);
r1  = zeros(m, 1);
epshat  = zeros(size(y, 1), n);
etahat  = zeros(size(QRt{1}, 1), n);
for t = n : -1 : 1
    if Tdyn, T = Tmat{t}; end
    etahat(:, t) = QRt{t}*r;
    if allmis(t)
        %% Disturbance smoothing when all observations are missing %%
        epshat(:, t) = 0;
        r = T'*r;
        if t <= d, r1 = T'*r1; end
    else
        if Hdyn, H = Hmat{t}; end
        if Zdyn, Z = Zmat{t}; end
        if anymis(t), Z(mis(:, t),:)=[]; H(mis(:, t),:)=[]; H(:,mis(:, t))=[]; end
        if t > d || ~Fns(t)
            %% Normal disturbance smoothing or when F_inf is zero %%
            epshat(~mis(:, t), t) = H*(invF{t}*v{t} - K{t}'*r);
            r = Z'*invF{t}*v{t} + L{t}'*r;
            if t <= d, r1 = T'*r1; end
        else
            %% Exact initial disturbance smoothing when F_inf is nonsingular %%
            epshat(~mis(:, t), t) = -H*K{t}'*r;
            r1 = Z'*invF{t}*v{t} + L{t}'*r1 + L1{t}'*r;
            r = L{t}'*r;
        end
        if anymis(t), if ~Zdyn, Z = Zmat; end, if ~Hdyn, H = Hmat; end, end
    end
end

%% Fast state smoothing %%
alphahat        = zeros(m, n);
alphahat(:, 1)  = a1 + P1*r + P1_inf*r1;
for t = 1 : n-1
    if Tdyn, T = Tmat{t}; end
    if Rdyn, R = Rmat{t}; end
    if cdyn, c = cmat{t}; end
    alphahat(:, t+1) = c + T*alphahat(:, t) + R*etahat(:, t);
end

repo: REST


text:
function out = innerProduct(f, g)
%INNERPRODUCT Compute the inner product of two DELTAFUN objects.
%   INNERPRODUCT(F, G) is the inner-product of two DELTAFUN object F and G.

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

%% Trivial cases:
if ( isempty(f) || isempty(g) )
    out = [];
    return
end

% Make sure both arguments are DELTAFUNS: (At least one must already be)
if ( ~isa(g, 'deltafun') )
    g = deltafun(g, [], []);
elseif ( ~isa(f, 'deltafun') )
    f = deltafun(f, [], []);
end

%%
% The innerProduct has four contributions:
% < (f + df) , (g + dg ) > = <f, g> + <df, g> + <dg, f> + <df, dg>

% <df, dg>: If two delta functions are at the same location, compute the
% appropriate infinity and return:
[fcomIdx, gcomIdx] = sameDeltaLocs(f, g);
if ( any(fcomIdx) )
    % If there is an overlap of locations, extract the columns corresponding to
    % the first overlap:
    fIdx = find(fcomIdx, 1);
    gIdx = find(gcomIdx, 1);
    df = f.deltaMag(:,fIdx);
    dg = g.deltaMag(:,gIdx);

    % Extract, the first non-zero product term:
    pref = chebfunpref();
    tol = pref.deltaPrefs.deltaTol;
    df = df(find(abs(df) > tol, 1));
    dg = dg(find(abs(dg) > tol, 1));
    % Assign a signed infinity:
    out = inf * sign(df*dg);
    return
end

% Delta functions don't overlap: compute the components of the inner product:
funIP = innerProduct(f.funPart, g.funPart);                   % <f, g>
dfIP = deltaInnerProduct(g.funPart, f.deltaMag, f.deltaLoc);  % <g, df>
dgIP = deltaInnerProduct(f.funPart, g.deltaMag, g.deltaLoc);  % <f, dg>
out = funIP + dfIP + dgIP;

end

function deltaIP = deltaInnerProduct(g, deltaMag, deltaLoc)
%DELTAINNERPRODUCT   Inner product of a DELTAFUN G with delta functions.

% Handle the empty case:
if ( isempty(deltaLoc) )
    deltaIP = 0;
    return
end

% Compute the derivatives needed:
m = size(deltaMag, 1);
maxDiffOrder = m - 1;
G = zeros(m, length(deltaLoc));
G(1,:) = feval(g, deltaLoc);
for k = 1:maxDiffOrder
    g = diff(g);
    G(k+1,:) = feval(g, deltaLoc);
end

% The output is always a scalar double:
deltaIP = 0;
v = ones(maxDiffOrder+1, 1);
v(2:2:end) = -1;
for k = 1:length(deltaLoc)
    % Apply the definition of inner product with delta functions:
    % <sum(ai dirac^(i)(x-xk)), f(x)> = sum ai*(-1)^i f^(i)(xk)
    ipk = (v.*deltaMag(:,k)).' * G(:,k) ;
    deltaIP = deltaIP + ipk;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = repext(varargin)
	varargout{1:max(nargout,1)} = extrep(varargin{:});
end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function [outputlabel, outputindex] = inputlabel2outputlabel(cfg, freq)

% INPUTLABEL2OUTPUTLABEL is a subfunction which outputs the cell-arrays
% outputlabel and the corresponding outputindex, and defines how the
% channels in the original data have to be combined, to provide the
% wished for combination of the channels, as defined in cfg.combinechan
%
% Configuration-options are:
%   cfg.combinechan = 'planar' combines the horizontal and vertical planar-gradients
%                     'pseudomeg' one gradiometer versus the rest
%   TODO: more flexible way of combining, e.g. by providing a cell-array

% 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$

if ~isfield(cfg, 'combinechan'), cfg.combinechan = 'no'; end

if strcmp(cfg.combinechan, 'no')
  % the output labels are similar to the input labels
  outputlabel = {};
  outputindex = {};
  for i=1:length(freq.label)
    outputindex{i} = i;
    outputlabel(i) = freq.label(i);
  end
elseif strcmp(cfg.combinechan, 'planar')
  % find the combination of horizontal and vertical channels that should be combined
  planar      = planarchannelset(freq);
  [sel_dH, H] = match_str(freq.label, planar(:,1));  % indices of the horizontal channels
  [sel_dV, V] = match_str(freq.label, planar(:,2));  % indices of the vertical   channels
  lab_dH      = freq.label(sel_dH);
  lab_dV      = freq.label(sel_dV);

  if length(sel_dH)~=length(sel_dV)
    [int,ih,iv] = intersect(H,V);
    sel_dH = sel_dH(ih);
    sel_dV = sel_dV(iv);
    %error('not all planar channel combinations are complete')
  else
    ih  = 1:length(sel_dH);
    iv  = 1:length(sel_dV);
    int = 1:length(ih);
  end

  % find the other channels that are present in the data
  sel_other = setdiff(1:length(freq.label), [sel_dH(:)' sel_dV(:)']);
  lab_other = freq.label(sel_other);

  % define the channel names after combining the planar combinations
  % they should be sorted according to the order of the horizontal planar channels in the data
  [dum, sel_planar] = match_str(freq.label, planar(:,1));
  if exist('ih', 'var'), sel_planar = sel_planar(int); end
  lab_comb          = planar(sel_planar,3);
  [dum, sel_comb]   = match_str(planar(H(ih),3),planar(V(iv),3));

  outputlabel = {};
  outputindex = {};

  outputlabel = [outputlabel; lab_other];
  outputindex = [outputindex; mat2cell(sel_other', ones(length(sel_other),1), 1)];

  outputlabel = [outputlabel; lab_comb];
  outputindex = [outputindex; mat2cell([sel_dH sel_dV(sel_comb)], ones(length(sel_dH),1), 2)];

elseif strcmp(cfg.combinechan, 'pseudomeg'),
  outputlabel = {};
  outputindex = {};
  for i=1:length(freq.label)
    outputlabel{i} = ['all-',freq.label{i}];
    outputindex{i} = [setdiff(1:length(freq.label),i)];
  end
  for i=1:length(freq.label)
    outputlabel{end+1} = freq.label{i};
    outputindex{end+1} = i;
  end
elseif iscell(cfg.combinechan(1)),
  for i=1:length(cfg.combinechan)
    outputindex{i} = match_str(freq.label, cfg.combinechan{i});
    outputindex{i} = outputindex{i}(:)';
    outputlabel{i} = cell2mat(freq.label(outputindex{i})');
  end
else
  ft_error('unknown combination method');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function iou = evaluateRoomObject3D(sequenceName,gtRoom3D,gtObject3D,outRoom3D,outObject3D,visualization,cutoff,delta)

%% set parameters

% set cutoff distance
if ~exist('cutoff','var')
    cutoff = 5.5;
end

% granuality
if ~exist('delta','var')
    delta  = 0.1;
end

if ~exist('visualization','var')
    visualization = false;
end

%% evaluate the accuracy
rangeX = min(min(outRoom3D(1,:)),min(gtRoom3D(1,:)))-delta:delta:max(max(outRoom3D(1,:)),max(gtRoom3D(1,:)))+delta;
rangeY = min(min(outRoom3D(2,:)),min(gtRoom3D(2,:)))-delta:delta:max(max(outRoom3D(2,:)),max(gtRoom3D(2,:)))+delta;
rangeZ = min(min(outRoom3D(3,:)),min(gtRoom3D(3,:)))-delta:delta:max(max(outRoom3D(3,:)),max(gtRoom3D(3,:)))+delta;
[gridX, gridY, gridZ] = ndgrid(rangeX,rangeY,rangeZ);
gridXYZ = [gridX(:)'; gridY(:)'; gridZ(:)'];
clear gridX gridY gridZ
%plot3(gridXYZ(1,:),gridXYZ(2,:),gridXYZ(3,:),'.');

data=readframe(fullfile('/n/fs/sun3d/data/',sequenceName));
imgSize= size(imread(data.depthpath));

extrinsicsFiles = dirSmart(fullfile('http://sun3d.cs.princeton.edu/data/',sequenceName,'extrinsics/'),'txt');
extrinsicsC2W = permute(reshape(readValuesFromTxt(fullfile('http://sun3d.cs.princeton.edu/data/',sequenceName,'extrinsics',extrinsicsFiles(end).name)),4,3,[]),[2 1 3]);

proj2D = data.K * extrinsicsC2W(1:3,1:3)' * gridXYZ;
proj2D(1:2,:) = proj2D(1:2,:) ./ proj2D([3 3],:);
inCamera = proj2D(3,:) > 1.0 & proj2D(3,:)< cutoff & 0<proj2D(1,:) & proj2D(1,:)<imgSize(2) & 0<proj2D(2,:) & proj2D(2,:)<imgSize(1);
gridXYZ = gridXYZ(:,inCamera);

if visualization

    figure

    %alpha(0.6);
    %title('Prediction Result');

    %figure
    numCorners = size(gtRoom3D,2)/2;
    patch(gtRoom3D(1,1:numCorners),gtRoom3D(2,1:numCorners),gtRoom3D(3,1:numCorners),'g');
    hold on
    patch(gtRoom3D(1,numCorners+1:end),gtRoom3D(2,numCorners+1:end),gtRoom3D(3,numCorners+1:end),'g');
    hold on
    for i=1:numCorners-1
        patch(gtRoom3D(1,[i i+1 numCorners+i+1 numCorners+i]),gtRoom3D(2,[i i+1 numCorners+i+1 numCorners+i]),gtRoom3D(3,[i i+1 numCorners+i+1 numCorners+i]),'g');
        hold on
    end
    patch(gtRoom3D(1,[numCorners 1  numCorners+1 numCorners*2]),gtRoom3D(2,[numCorners 1  numCorners+1 numCorners*2]),gtRoom3D(3,[numCorners 1  numCorners+1 numCorners*2]),'g');
    axis equal

    for j=1:length(gtObject3D.objects)
        if ~isempty(gtObject3D.objects{j})
            plot3(gtObject3D.objects{j}.polygon{1}.X([1:end 1]), repmat(gtObject3D.objects{j}.polygon{1}.Ymin,1,5), gtObject3D.objects{j}.polygon{1}.Z([1:end 1]),'-k','LineWidth',2);
            plot3(gtObject3D.objects{j}.polygon{1}.X([1:end 1]), repmat(gtObject3D.objects{j}.polygon{1}.Ymax,1,5), gtObject3D.objects{j}.polygon{1}.Z([1:end 1]),'-k','LineWidth',2);
            for k=1:4
                plot3(gtObject3D.objects{j}.polygon{1}.X([k k]),[gtObject3D.objects{j}.polygon{1}.Ymin,gtObject3D.objects{j}.polygon{1}.Ymax], gtObject3D.objects{j}.polygon{1}.Z([k k]),'-k','LineWidth',2);
            end
        end
    end
    xlabel('x');
    ylabel('y');
    zlabel('z');

    numCorners = size(outRoom3D,2)/2;
    patch(outRoom3D(1,1:numCorners),outRoom3D(2,1:numCorners),outRoom3D(3,1:numCorners),'k');
    hold on
    patch(outRoom3D(1,numCorners+1:end),outRoom3D(2,numCorners+1:end),outRoom3D(3,numCorners+1:end),'b');
    hold on
    for i=1:numCorners-1
        patch(outRoom3D(1,[i i+1 numCorners+i+1 numCorners+i]),outRoom3D(2,[i i+1 numCorners+i+1 numCorners+i]),outRoom3D(3,[i i+1 numCorners+i+1 numCorners+i]),'r');
        hold on
    end
    patch(outRoom3D(1,[numCorners 1  numCorners+1 numCorners*2]),outRoom3D(2,[numCorners 1  numCorners+1 numCorners*2]),outRoom3D(3,[numCorners 1  numCorners+1 numCorners*2]),'r');
    axis equal

    alpha(0.01);

    for j=1:length(outObject3D.objects)
        if ~isempty(outObject3D.objects{j})
            plot3(outObject3D.objects{j}.polygon{1}.X([1:end 1]), repmat(outObject3D.objects{j}.polygon{1}.Ymin,1,5), outObject3D.objects{j}.polygon{1}.Z([1:end 1]),'-b','LineWidth',2);
            plot3(outObject3D.objects{j}.polygon{1}.X([1:end 1]), repmat(outObject3D.objects{j}.polygon{1}.Ymax,1,5), outObject3D.objects{j}.polygon{1}.Z([1:end 1]),'-b','LineWidth',2);
            for k=1:4
                plot3(outObject3D.objects{j}.polygon{1}.X([k k]),[outObject3D.objects{j}.polygon{1}.Ymin,outObject3D.objects{j}.polygon{1}.Ymax], outObject3D.objects{j}.polygon{1}.Z([k k]),'-b','LineWidth',2);
            end
        end
    end
    xlabel('x');
    ylabel('y');
    zlabel('z');

    %hold on;
    %plot3(gridXYZ(1,:),gridXYZ(2,:),gridXYZ(3,:),'.'); xlabel('x'); ylabel('y'); zlabel('z');
end

numCorners = size(gtRoom3D,2)/2;

%% inside the ground truth
[inon,on] = inpolygon(gridXYZ(1,:),gridXYZ(3,:),gtRoom3D(1,1:numCorners),gtRoom3D(3,1:numCorners));
in = inon&~on;
inGroundTruth = in & min(gtRoom3D(2,:)) < gridXYZ(2,:) & gridXYZ(2,:) < max(gtRoom3D(2,:));

intersectCount = zeros(1,length(inGroundTruth));
for i=find(inGroundTruth)
    [xi,yi] = polyxpoly([0 gridXYZ(1,i)],[0 gridXYZ(3,i)],gtRoom3D(1,1:numCorners),gtRoom3D(3,1:numCorners));
    if ~isempty(xi)
        count = 0;
        for k=1:length(xi)
            disSq = xi(k)^2+yi(k)^2;
            if disSq > 0.5^2
                count=count+1;
            end
        end
        intersectCount(i) = count;
    end
end
inGroundTruth = inGroundTruth & (intersectCount==0);

for j=1:length(gtObject3D.objects)
    if ~isempty(gtObject3D.objects{j})
        [inon, on] = inpolygon(gridXYZ(1,:),gridXYZ(3,:), gtObject3D.objects{j}.polygon{1}.X, gtObject3D.objects{j}.polygon{1}.Z);
        in = inon&~on;
        in = in & gtObject3D.objects{j}.polygon{1}.Ymin < gridXYZ(2,:) & gridXYZ(2,:) < gtObject3D.objects{j}.polygon{1}.Ymax;
        inGroundTruth = inGroundTruth & ~in;
    end
end

if visualization
    hold on;
    plot3(gridXYZ(1,inGroundTruth),gridXYZ(2,inGroundTruth),gridXYZ(3,inGroundTruth),'.b');
end

%% inside the prediction results
[inon,on] = inpolygon(gridXYZ(1,:),gridXYZ(3,:),outRoom3D(1,1:end/2),outRoom3D(3,1:end/2));
in = inon&~on;
inResult = in & min(outRoom3D(2,:)) < gridXYZ(2,:) & gridXYZ(2,:) < max(outRoom3D(2,:));

intersectCount = zeros(1,length(inResult));
for i=find(inResult)
    [xi,yi] = polyxpoly([0 gridXYZ(1,i)],[0 gridXYZ(3,i)],outRoom3D(1,1:end/2),outRoom3D(3,1:end/2));
    if ~isempty(xi)
        count = 0;
        for k=1:length(xi)
            disSq = xi(k)^2+yi(k)^2;
            if disSq > 0.5^2
                count=count+1;
            end
        end
        intersectCount(i) = count;
    end
end
inResult = inResult & (intersectCount==0);

for j=1:length(outObject3D.objects)
    if ~isempty(outObject3D.objects{j})
        [inon, on] = inpolygon(gridXYZ(1,:),gridXYZ(3,:), outObject3D.objects{j}.polygon{1}.X, outObject3D.objects{j}.polygon{1}.Z);
        in = inon&~on;
        in = in & outObject3D.objects{j}.polygon{1}.Ymin < gridXYZ(2,:) & gridXYZ(2,:) < outObject3D.objects{j}.polygon{1}.Ymax;
        inResult = inResult & ~in;
    end
end

if visualization
    hold on;
    plot3(gridXYZ(1,inResult),gridXYZ(2,inResult),gridXYZ(3,inResult),'.k');
end

%% evaluate
intersectionGR = inGroundTruth & inResult;
unionGR = inGroundTruth | inResult;
iou = sum(intersectionGR) / sum(unionGR);

author: thusiyuan


text:
%% A Unit Test Class for confidence interval and fit statistics
classdef confidence_tests < matlab.unittest.TestCase

    properties
        absTol          = 1e-10;
        covAbsTol       = 1e-6;
        confIntAbsTol   = 1e-6;
        confBndsRelTol  = 1e-6;
        pValAbsTol      = 1e-8;
        stdErrRelTol    = 1e-5;
        tStatRelTol     = 1e-5;
        data   = load('Utilities/UnitTests/UnitTestData/confidenceData.mat');
    end

    % Unit Tests
    methods (Test)

        %-- Himmelblau - no weights --%
        function himmelblau(testCase)
            mdl = testCase.data.himmel;
            mdl.fcn = @confidence_tests.himmelblauFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- Himmelblau - with weights --%
        function himmelblauWeights(testCase)
            mdl = testCase.data.himmelw;
            mdl.fcn = @confidence_tests.himmelblauFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- Himmelblau - with lower bound active --%
        function himmelblauLb(testCase)
            mdl = testCase.data.himmellb;
            mdl.fcn = @confidence_tests.himmelblauFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- Himmelblau - with upper bound active --%
        function himmelblauUb(testCase)
            mdl = testCase.data.himmelub;
            mdl.fcn = @confidence_tests.himmelblauFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- Himmelblau - with both bounds active --%
        function himmelblauBothBounds(testCase)
            mdl = testCase.data.himmelbnd;
            mdl.fcn = @confidence_tests.himmelblauFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- Himmelblau - with both bounds and weights --%
        function himmelblauBothBoundsWeights(testCase)
            mdl = testCase.data.himmelbndw;
            mdl.fcn = @confidence_tests.himmelblauFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- SAS - no weights --%
        function sas(testCase)
            mdl = testCase.data.sas;
            mdl.fcn = @confidence_tests.sasFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end

        %-- SAS - with weights --%
        function sasWeights(testCase)
            mdl = testCase.data.sasw;
            mdl.fcn = @confidence_tests.sasFcn;
            % Check fit statistics
            checkFitStats(testCase,mdl);
        end
    end

    methods
       function checkFitStats(testCase, mdl)
            % Build OPTI object
            optiObj = opti('fun',mdl.fcn, 'data', mdl.xdata, mdl.ydata, 'weights', mdl.wts, 'x0', mdl.sol, 'bounds', mdl.lb, mdl.ub);
            % Manually set the solution
            optiObj.setSolution(mdl.SSE, mdl.sol);
            % Calculate the fit statistics;
            stats  = calcStatistics(optiObj, 0.95, false);

            % Check statistcs
            testCase.verifyEqual(mdl.Rsquare, stats.Rsquare, 'AbsTol', testCase.absTol);
            testCase.verifyEqual(mdl.AdjRsquare, stats.AdjRsquare, 'AbsTol', testCase.absTol);
            testCase.verifyEqual(mdl.RMSE, stats.RMSE, 'AbsTol', testCase.absTol);
            testCase.verifyEqual(mdl.DFE, stats.DFE, 'AbsTol');
            if (isfield(mdl, 'cov'))
                testCase.verifyEqual(mdl.cov, stats.Cov, 'AbsTol', testCase.covAbsTol);
            end
            testCase.verifyEqual(mdl.confInt, stats.ConfInt, 'AbsTol', testCase.confIntAbsTol);
            testCase.verifyEqual(mdl.confBnds, stats.ConfBnds.bnds, 'RelTol', testCase.confBndsRelTol);
            if (isfield(mdl, 'param'))
                testCase.verifyEqual(mdl.param(:,2), stats.Param.StdError, 'RelTol', testCase.stdErrRelTol);
                testCase.verifyEqual(mdl.param(:,3), stats.Param.tStat, 'RelTol', testCase.tStatRelTol);
                testCase.verifyEqual(mdl.param(:,4), stats.Param.pValues, 'AbsTol', testCase.pValAbsTol);
            end
        end
    end

    methods (Static)
        function r = himmelblauFcn(theta,p)
            % r = f(x,theta) (x = pressure, theta unknowns)
            r = theta(1)*p./(1+theta(2)*p);
        end

        function k = sasFcn(theta,n)
            k = theta(1)*n./(theta(2) + n);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function data = rle(x)
% data = rle(x) (de)compresses the data with the RLE-Algorithm
%   Compression:
%      if x is a numbervector data{1} contains the values
%      and data{2} contains the run lenths
%
%   Decompression:
%      if x is a cell array, data contains the uncompressed values
%
%      Version 1.0 by Stefan Eireiner (stefan-e@web.de)
%      based on Code by Peter J. Acklam
%      last change 14.05.2004

if iscell(x) % decoding
	i = cumsum([ 1 x{2} ]);
	j = zeros(1, i(end)-1);
	j(i(1:end-1)) = 1;
	data = x{1}(cumsum(j));
else % encoding
	if size(x,1) > size(x,2), x = x'; end % if x is a column vector, tronspose
    i = [ find(x(1:end-1) ~= x(2:end)) length(x) ];
	data{2} = diff([ 0 i ]);
	data{1} = x(i);
end


text:
function [F,G,rho,u,v,p] = EulerFluxes2D(Q, gamma)

% function [F,G,rho,u,v,p] = EulerFluxes2D(Q, gamma)
% Purpose: evaluate primitive variables and Euler flux functions

% extract conserved variables
rho = Q(:,:,1); rhou = Q(:,:,2); rhov = Q(:,:,3); Ener = Q(:,:,4);

% compute primitive variables
u = rhou./rho; v = rhov./rho; p = (gamma-1)*(Ener - 0.5*(rhou.*u + rhov.*v));

% compute flux functions
F = zeros(size(Q));
F(:,:,1) = rhou;
F(:,:,2) = rhou.*u + p;
F(:,:,3) = rhov.*u;
F(:,:,4) = u.*(Ener+p);

G = zeros(size(Q));
G(:,:,1) = rhov;
G(:,:,2) = rhou.*v;
G(:,:,3) = rhov.*v + p;
G(:,:,4) = v.*(Ener+p);
return;

repo: Aero-matlab


text:
function G = compute_periodic_poisson(d, symmetrize)

% compute_periodic_poisson - solve poisson equation
%
%   G = compute_periodic_poisson(d,symmetrize);
%
%   Solve
%       Delta(G) = d
%   with periodic boundary condition.
%   G has zero mean.
%
%   Set symmetrize=1 (default 0) if the data divergence d is not periodic.
%   This will double the size and consider symmetric extension.
%
%   Copyright (c) 2007 Gabriel Peyre

n = size(d,1);

if nargin==2 && symmetrize
    d = perform_size_doubling(d);
    G = compute_periodic_poisson(d);
    G = G(1:n,1:n);
    return;
end

% solve for laplacian
[Y,X] = meshgrid(0:n-1,0:n-1);
mu = sin(X*pi/n).^2; mu = -4*( mu+mu' );
mu(1) = 1; % avoid division by 0

G = fft2(d) ./ mu; G(1) = 0;
G = real( ifft2( G ) );

%%
function g = perform_size_doubling(g)

g = [g;g(end:-1:1,:,:)];
g = [g,g(:,end:-1:1,:)];

author: gpeyre

repo: matlab-toolboxes


text:
%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym laplacian (@var{f})
%% @defmethodx @@sym laplacian (@var{f}, @var{x})
%% Symbolic Laplacian of symbolic expression.
%%
%% The Laplacian of a scalar expression @var{f} is
%% the scalar expression:
%% @example
%% @group
%% syms f(x, y, z)
%% laplacian(f)
%%   @result{} (sym)
%%         2                 2                 2
%%        ∂                 ∂                 ∂
%%       ───(f(x, y, z)) + ───(f(x, y, z)) + ───(f(x, y, z))
%%         2                 2                 2
%%       ∂x                ∂y                ∂z
%% @end group
%% @end example
%%
%% @var{x} can be a scalar, vector or cell list.  If omitted,
%% it is determined using @code{symvar}.
%%
%% Example:
%% @example
%% @group
%% syms x y
%% laplacian(x^3 + 5*y^2)
%%   @result{} (sym) 6⋅x + 10
%% @end group
%% @end example
%%
%% Note: assumes @var{x} is a Cartesian coordinate system.
%%
%% @seealso{@@sym/divergence, @@sym/gradient, @@sym/curl, @@sym/jacobian,
%%          @@sym/hessian}
%% @end defmethod

function g = laplacian(f,x)

  assert (isscalar(f), 'laplacian: only scalar functions supported')

  if (nargin == 1)
    x = symvar(f);
    if (isempty(x))
      x = sym('x');
    end
  elseif (nargin == 2)
    % no-op
  else
    print_usage ();
  end

  if (~iscell(x) && isscalar(x))
    x = {x};
  end

  cmd = { '(f, x) = _ins'
          'g = 0'
          'for y in x:'
          '    g = g + f.diff(y, 2)'
          'return g,' };

  g = pycall_sympy__ (cmd, sym(f), x);

end

%!shared x,y,z
%! syms x y z

%!test
%! % 1D
%! f = x^2;
%! g = diff(f,x,x);
%! assert (isequal (laplacian(f), g))
%! assert (isequal (laplacian(f,{x}), g))
%! assert (isequal (laplacian(f,[x]), g))
%! assert (isequal (laplacian(f,x), g))

%!test
%! % const
%! f = sym(1);
%! g = sym(0);
%! assert (isequal (laplacian(f), g))
%! assert (isequal (laplacian(f,x), g))
%! f = sym('c');
%! assert (isequal (laplacian(f,x), g))

%!test
%! % double const
%! f = 1;
%! g = sym(0);
%! assert (isequal (laplacian(f,x), g))

%!test
%! % 1D fcn in 2d/3d
%! f = sin(2*y);
%! g = -4*f;
%! assert (isequal (laplacian(f), g))
%! assert (isequal (laplacian(f, {x,y}), g))
%! assert (isequal (laplacian(f, {x,y,z}), g))

%!test
%! % 2d fcn in 2d/3d
%! f = sin(exp(x)*y);
%! g = diff(f,x,x) + diff(f,y,y);
%! assert (isequal (laplacian(f), g))
%! assert (isequal (laplacian(f, {x,y}), g))

%!test
%! % 2d fcn in 2d/3d
%! f = sin(exp(x)*y+sinh(z));
%! gr2 = gradient(f, {x,y});
%! divgr2 = divergence(gr2, {x,y});
%! l2 = laplacian(f,{x,y});
%! gr3 = gradient(f, {x,y,z});
%! divgr3 = divergence(gr3, {x,y,z});
%! l3 = laplacian(f,{x,y,z});
%! assert (isAlways (l2 == divgr2))
%! assert (isAlways (l3 == divgr3))

%!error laplacian(sym('x'), sym('x'), 42)
%!error <only scalar> laplacian([sym('x'), sym('x')])

repo: octsympy


text:
function drawlayers(varargin)
%*************************************************
%DATE: 11/8/2007 (created); 	11/8/2007 (modified)
%FUNCTION:	draw layers according to Layerstructure
%INPUTS:	layer structure matrix
%OUTPUTS:	NULL
%muwn.gu@gmail.com
%*************************************************
global Layerprofile;
global layerMatrix;
global layerColor;
layerNrm =  max(max(layerMatrix));
matrixSize =  size(layerMatrix);
lineNrm   =  matrixSize(1,1);
lineLength=  matrixSize(1,2);

for lineIndex = 1: lineNrm	% for layers

	for realLayerIndex = 1 : layerNrm    % for different layers

		if (layerMatrix(lineIndex, 1) == realLayerIndex)
			seg.from = 0;
		end

		for pointIndex = 1 : (lineLength-1)	% for point in one layer

			if (~(layerMatrix(lineIndex, pointIndex) == realLayerIndex)) & (layerMatrix(lineIndex, pointIndex + 1) == realLayerIndex)
			seg.from = (pointIndex)/(lineLength);
			end

			if (layerMatrix(lineIndex, pointIndex) == realLayerIndex) & (~(layerMatrix(lineIndex, pointIndex + 1) == realLayerIndex))
			seg.end = (pointIndex)/(lineLength);
			fillColor	= getlayercolor (realLayerIndex);
			drawsegment (seg, fillColor);
			end
		end

		if (layerMatrix(lineIndex, lineLength) == realLayerIndex)
			seg.end = 1;
			fillColor	= getlayercolor (realLayerIndex);
			drawsegment (seg, fillColor);
		end
	end
			Layerprofile.currentY = Layerprofile.currentY + Layerprofile.lineWidth;
end

Layerprofile.currentY =  Layerprofile.currentY + Layerprofile.stepSpace;

function drawsegment(seg, fillcolor)
%*************************************************
%DATE: 11/8/2007 (created); 	11/8/2007 (modified)
%FUNCTION:
%INPUTS:
%OUTPUTS:
%*************************************************
global Layerprofile;
			xfrom 	= Layerprofile.diameter * seg.from;
			xrange 	= Layerprofile.diameter * seg.end - xfrom;
			rectangle('Position',[Layerprofile.origin.x + xfrom, Layerprofile.currentY, xrange , Layerprofile.lineWidth], 'FaceColor', fillcolor, 'LineStyle', 'none')

function fillColor = getlayercolor(layerIndex)
%*************************************************
%FUNCTION:
%INPUTS:
%OUTPUTS:
%*************************************************
global 	Layerprofile;
switch	layerIndex
	case 1
		fillColor =  Layerprofile.color.silicon;
	case 2
		fillColor =  Layerprofile.color.sidioxi;
	case 3
		fillColor =  Layerprofile.color.metal;
	case 4
		fillColor =  Layerprofile.color.sini;
	case 5
		fillColor =  Layerprofile.color.posi;
	case 9
		fillColor =  Layerprofile.color.phore;
	otherwise
		fillColor = [1, 1, 1];
end


text:
function demo_blockproc_slidingerblets(source,varargin) %RUNASSCRIPT
%DEMO_BLOCKPROC_SLIDINGERBLETS Basic real-time rolling erblet-spectrogram visualization
%   Usage: demo_blockproc_slidingerblets('gspi.wav')
%
%   For additional help call |demo_blockproc_slidingerblets| without arguments.
%
%   This demo shows a simple rolling erblet-spectrogram of whatever is specified in
%   source.

if demo_blockproc_header(mfilename,nargin)
   return;
end

% Control pannel (Java object)
% Each entry determines one parameter to be changed during the main loop
% execution.
p = blockpanel({
               {'GdB','Gain',-20,20,0,21},...
               {'cMult','C mult',0,80,20,81}
               });

fobj = blockfigure();

% Setup blocktream
try
    fs=block(source,varargin{:},'loadind',p);
catch
    % Close the windows if initialization fails
    blockdone(p,fobj);
    err = lasterror;
    error(err.message);
end

% Buffer length (30 ms)
bufLen = floor(30e-3*fs);
zpad = floor(bufLen/2);

% Number of filters
M = 200;
F = frame('erbletfb',fs,2*bufLen+2*zpad,'fractionaluniform','M',M);
Fa = blockframeaccel(F,bufLen,'sliced','zpad',zpad);

flag = 1;
cola = [];
%Loop until end of the stream (flag) and until panel is opened
while flag && p.flag
  % Get parameters
  [gain,mult] = blockpanelget(p,'GdB','cMult');
  gain = 10^(gain/20);
  mult = 10^(mult/20);

  % Read block of length bufLen
  [f,flag] = blockread(bufLen);
  f = f*gain;
  % Apply analysis frame
  c = blockana(Fa, f);
  % Plot
  cola = blockplot(fobj,Fa,mult*c(:,1),cola);

  blockplay(f);
end
blockdone(p,fobj,Fa);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [Mc, Mc95, Mc90] = calc_McBest(magnitudes, binInterval)
    % CALC_MCBEST calculate best Magnitude of Completion
    % [fMc, fMc95, fMc90] = calc_McBest(magnitudes, binInterval)
    % each row of magitudes gets its own calculation

    % report_this_filefun();

    % Magnitude increment
    if ~exist('binInterval','var')
        binInterval = 0.1;
    end
    magnitudes = sort(magnitudes);

    half_bin   = binInterval / 2;
    % First estimation of magnitude of completeness (maximum curvature)
    McStarts = maxCurvature(magnitudes, [-2 : binInterval : 6] , half_bin); % vectorized

    % cheat to allow this to handle multiple rows of magnitudes (independently)
    for jj = 1:size(magnitudes,2)
        valid = ~isnan(magnitudes(:,jj));
        [Mc(jj), Mc95(jj), Mc90(jj)] = calc_McBest_unvectorized(McStarts(jj), magnitudes(valid,jj), binInterval, half_bin);
    end
end

function  [fMc, fMc95, fMc90] = calc_McBest_unvectorized(McStart, magnitudes, binInterval, half_fBin)
    magCenters  = (McStart - 0.9) : binInterval : (McStart + 1.5);
    eachedge    = [magCenters - half_fBin , magCenters(end) + half_fBin];

    magnitudes = flipud(magnitudes); % from biggest to smallest
    nGtEdge = sum(magnitudes > eachedge((end-1) : -1 : 1)); % magnitudes(nGtEdge) gives last event above threshhold
    nGtEdge = fliplr(nGtEdge);
    too_few = nGtEdge < 25;

    results=nan(numel(magCenters),1);
    for idx = 1:numel(magCenters)
        if too_few(idx)
            continue
        end
        hypotheticalMc  = magCenters(idx);
        results(idx) = doCalculation(magnitudes(1:nGtEdge(idx)), binInterval, hypotheticalMc);
    end

    % Evaluation of results

    % Is fMc90 available
    nSel = find(results < 10, 1, 'first' );
    if isempty(nSel)
        fMc90 = NaN;
    else
        fMc90 = magCenters(nSel);
    end

    % Is fMc95 available
    nSel = find(results < 5, 1 );
    if isempty(nSel)
        fMc95 = NaN;
    else
        fMc95 = magCenters(nSel);
    end

    % take results from bins. (I tested against discretize, and this was faster -CGR)

    j =  find(results < 10 , 1, 'first');
    if isempty(j)
        j =  find(results < 15 , 1, 'first' );
    end
    if isempty(j)
        j =  find(results < 20 , 1, 'first' );
    end
    if isempty(j)
        j =  find(results < 25 , 1, 'first' );
    end
    fMc = magCenters(j);
    if isempty(fMc)
        fMc = NaN;
    end
end

function [result] = doCalculation(theseMags, binInterval, hypotheticalMc)
        fBValue   = calc_bmemag(theseMags, binInterval);
        half_fBin = binInterval ./ 2;

        % log10(N)=A-B*M
        vMag    = hypotheticalMc:binInterval:15; % Ending magnitude must be sufficiently high (???what is "sufficently high")
        vNumber = 10.^(log10(numel(theseMags)) - fBValue*(vMag - hypotheticalMc));
        vNumber = round(vNumber);

        % PM=vMag(1:ct);
        PMedges     = [vMag-half_fBin , vMag(end)+half_fBin];
        [bval, ~]   = histcounts(theseMags, PMedges);
        b3          = cumsum(bval,'reverse');    % N for M >= (counted backwards)
        result      = sum(abs(b3 - vNumber)) / sum(b3)*100; %res2
end

function Mc = maxCurvature(m, centers, halfBinwidth)
    % MAXCURVATURE First estimation of magnitude of completeness (maximum curvature)
    % Mc = MAXCURVATURE(m, min_max, binwidth)
    %
    % vectorized
    edges           = [centers - halfBinwidth , centers(end) + halfBinwidth];
    [vEvents, ~]    = histc(m, edges, 1);
    [~,idx] = max(flipud(vEvents));
    nSel = size(vEvents,1) - idx + 1;
    % nSel    = find(vEvents == max(vEvents), 1, 'last' );
    Mc      = centers(nSel);
end


text:
function DatsMetadata = bz_DatFileMetadata(basepath)
% Store basic info about the original dat files: names and bytes
% Brendon Watson 2016-8

%% Input and directory handling
if ~exist('basepath','var')
    basepath = cd;
elseif isempty(basepath)
    basepath = cd;
end

basename = bz_BasenameFromBasepath(basepath);

%% First find out what system was used to record this... based on files present
% check if amplipex (.metas present) or intan (rhd.infos present)... or neither
recsys = '';

%Amplipex case
dmeta = dir(fullfile(basepath,'*.meta'));
if ~isempty(dmeta)
    recsys = 'Amplipex';
end

%Intan case
drhd = dir(basepath);
for didx = length(drhd):-1:1
   ok = 0;
   if drhd(didx).isdir
       t = dir(fullfile(basepath,drhd(didx).name,'info.rhd'));
       if ~isempty(t)%only use folders with info.rhd inside
            recsys = 'Intan';
            ok = 1;
       end
   end
   if ~ok
        drhd(didx) = [];
   end
end

%% Now gather info
switch recsys
    case 'Amplipex'%if from an amplipex
        d = dmeta;
        DatsMetadata.RecordingSystem = 'Amplipex';
        for idx = 1:length(d);
            DatsMetadata.Recordings.Names{idx} = d(idx).name(1:end-5);
            tdat = strcat(d(idx).name(1:end-5),'.dat');
            t = dir(fullfile(basepath,tdat));
            if ~isempty(t) %if original .dat still around
                DatsMetadata.Recordings.Bytes(idx) = t.bytes;
            else%if no .dat around (ie anymore), get bytes from .meta
                disp([tdat ' not found']);
                DatsMetadata.Recordings.Bytes{idx} = str2num(bz_ReadAmplipexMetafileAspects(fullfile(basepath,d(idx).name),'filebytes'));
            end
            NChannels(idx) = str2num(bz_ReadAmplipexMetafileAspects(fullfile(basepath,d(idx).name),'NumChans'));
            SamplingRate(idx) = str2num(bz_ReadAmplipexMetafileAspects(fullfile(basepath,d(idx).name),'SamplingRate'));
            Amplification(idx) = str2num(bz_ReadAmplipexMetafileAspects(fullfile(basepath,d(idx).name),'Gain'));
        end
        if ~isempty(NChannels) %if some files were found
            %Channels
            if sum(abs(diff(NChannels)))%if different numbers of channels differed across recordings in session, error out
                warning('Rercordings contain different numbers of channels - per .meta files.  Will use count from first recording.')
            end
            DatsMetadata.Parameters.NumberOfChannels = NChannels(1);

            %Sampling Rate
            if sum(abs(diff(SamplingRate)))%if different numbers of channels differed across recordings in session, error out
                warning('Rercordings contain different SamplingRates - per .meta files.  Will use number from first recording.')
            end
            DatsMetadata.Parameters.SampleRate = SamplingRate(1);

            %Amplification
            if sum(abs(diff(Amplification)))%if different numbers of channels differed across recordings in session, error out
                warning('Rercordings contain different Amplifications - per .meta files.  Will use number from first recording.')
            end
            DatsMetadata.Parameters.Amplification = Amplification(1);

            DatsMetadata.Parameters.VoltsPerUnit = VoltsPerUnit_Amplirec(basename,basepath);%calculate from .metas/.inis
            DatsMetadata.Parameters.BitsPerSample = 16;%amplirec default
            DatsMetadata.Parameters.VoltageRange = 10;%not used except to make xml
        end
    case 'Intan' %if not amplipex, look for intan
        d = drhd;
        NAmpChannels = [];
        NAuxChannels = [];
        NDigitalChannels = [];
        SamplingRate = [];
        thisidx = 0;

        for didx = 1:length(d)
            if d(didx).isdir
                t = dir(fullfile(basepath,d(didx).name,'info.rhd'));
                if ~isempty(t)%only use folders with info.rhd inside

                   DatsMetadata.RecordingSystem = 'Intan';
                   rhd = fullfile(basepath,d(didx).name,'info.rhd');
    %                ampdats{end+1} = fullfile(basepath,d(didx).name,'amplifier.dat');
    %                recordingdiridxs = didx;

                   thisidx = thisidx + 1;

                   DatsMetadata.Recordings.Names{thisidx} = d(didx).name;

                   % read the info.rhd file for each individual recorded
                   % file/folder
                   try
                      [amplifier_channels, notes, aux_input_channels, spike_triggers,...
                          board_dig_in_channels, supply_voltage_channels, frequency_parameters ] =...
                          read_Intan_RHD2000_file(fullfile(basepath,d(didx).name),'info.rhd');%function from Intan
                      DatsMetadata.Recordings.IntanRHDInfoRaw{thisidx} = v2struct(amplifier_channels, notes, aux_input_channels, spike_triggers,...
                          board_dig_in_channels, supply_voltage_channels, frequency_parameters);
                      NAmpChannels(thisidx) = length(amplifier_channels);
                      NAuxChannels(thisidx) = length(aux_input_channels);
                      NDigitalChannels(thisidx) = length(board_dig_in_channels);
                      SamplingRate(thisidx) = frequency_parameters.amplifier_sample_rate;

                   catch%if read error, set values based on prior
                       if ~isempty(NAmpChannels)
%                            if ~isnan(NAmpChannels(1));
                          NAmpChannels(thisidx) = NAmpChannels(1);
                       else
                          sess = bz_getSessionInfo(basepath,'noPrompts',true);
                          NAmpChannels(thisidx) = sess.nChannels;
                       end
                       if ~isempty(NAuxChannels)
                           NAuxChannels(thisidx) = NAuxChannels(1);
                       else
                           NAuxChannels(thisidx) = nan;
                       end
                       if ~isempty(NDigitalChannels)
                           NDigitalChannels(thisidx) = NDigitalChannels(1);
                       else
                           NDigitalChannels(thisidx) = nan;
                       end
                       if ~isempty(SamplingRate)
                           SamplingRate(thisidx) = SamplingRate(1);
                       else
                           SamplingRate(thisidx) = nan;
                       end
                   end

                   t = dir(fullfile(basepath,d(didx).name,'amplifier.dat'));
                   t2 = dir(fullfile(basepath,d(didx).name,'amplifier_analogin_auxiliary_int16.dat'));
                   if ~isempty(t2) %if original .dat still around.  Handling Luke's RHD software to put all 16bit into same file
                       DatsMetadata.Recordings.Bytes(thisidx) = t2.bytes;
                   elseif ~isempty(t) %if original .dat still around
                       DatsMetadata.Recordings.Bytes(thisidx) = t.bytes;
                   else%if no .dat around
                       disp([DatsMetadata.Recordings.Names{end} ' not found']);
                       timedat = dir(fullfile(basepath,d(didx).name,'time.dat'));
                       if ~isnan(NAmpChannels(end));
                          DatsMetadata.Recordings.Bytes(thisidx) = timedat(1).bytes/2*NAmpChannels(end);
                       else
                           sess = bz_getSessionInfo(basepath);
                           DatsMetadata.Recordings.Bytes(thisidx) = timedat(1).bytes/2*sess.nChannels;
                       end
                   end

               end
            end
        end
        %handle errors created by improper read of rhd.info files
        if isnan(sum(abs(diff(NAmpChannels))))
            warning('At least one read error on rhd.info files.  Using values from other files to populate.')
            ok = find(~isnan(NAmpChannels),1,'first');
            if isempty(ok)
                warning('No good reads on NumberofAmpChannels. Assuming based on probe maps')
                sess = bz_getSessionInfo(basepath);
                NAmpChannels(isnan(NAmpChannels)) = sess.nChannels;
            else
                NAmpChannels(isnan(NAmpChannels)) = NAmpChannels(ok);
            end
        end
        if isnan(sum(abs(diff(NAuxChannels))))
            warning('At least one read error on rhd.info files.  Using values from other files to populate.')
            ok = find(~isnan(NAuxChannels),1,'first');
            if isempty(ok)
                error('No good reads on NAuxChannels, assuming 0')
                NAuxChannels(isnan(NAuxChannels)) = 0;
            else
                NAuxChannels(isnan(NAuxChannels)) = NAuxChannels(ok);
            end
        end
        if isnan(sum(abs(diff(NDigitalChannels))))
            warning('At least one read error on rhd.info files.  Using values from other files to populate.')
            ok = find(~isnan(NDigitalChannels),1,'first');
            if isempty(ok)
                error('No good reads on NDigitalChannels, assuming 0')
                NDigitalChannels(isnan(NDigitalChannels)) = 0;
            else
                NDigitalChannels(isnan(NDigitalChannels)) = NDigitalChannels(ok);
            end
        end
        if isnan(sum(abs(diff(SamplingRate))))
            warning('At least one read error on rhd.info files.  Using values from other files to populate.')
            ok = find(~isnan(SamplingRate),1,'first');
            if isempty(ok)
                error('No good reads on SamplingRate, taking from AnimalMetadata')
                sess = bz_getSessionInfo(basepath);
                SamplingRate(isnan(SamplingRate)) = sess.rates.wideband;
            else
                SamplingRate(isnan(SamplingRate)) = SamplingRate(ok);
            end
        end

        %handle errors where different parameters read for different files (rather than uniform)
        if sum(abs(diff(NAmpChannels)))%if different numbers of channels differed across recordings in session, error out
            warning('Rercordings contain different numbers of channels - per info.rhd files.  Will use count from first recording.')
        end
        DatsMetadata.Parameters.NumberOfChannels = NAmpChannels(1);

        if sum(abs(diff(NAuxChannels)))%if different numbers of channels differed across recordings in session, error out
            warning('Rercordings contain different numbers auxiliary of channels - per info.rhd files.  Will use count from first recording.')
        end
        DatsMetadata.NumAuxiliaryChannels = NAuxChannels(1);

        if sum(abs(diff(NDigitalChannels)))%if different numbers of channels differed across recordings in session, error out
            warning('Rercordings contain different numbers digital of channels - per info.rhd files.  Will use count from first recording.')
        end
        DatsMetadata.NumDigitalChannels = NAuxChannels(1);

        %sampling rate
        if sum(abs(diff(SamplingRate)))%if different numbers of channels differed across recordings in session, error out
            warning('Rercordings contain different SamplingRates - per info.rhd files.  Will use number from first recording.')
        end
        DatsMetadata.Parameters.SampleRate = SamplingRate(1);

        % Not gather-able from intan info.rhd files
        DatsMetadata.Parameters.Amplification = 1;%digitized on chip, let's say 1
        DatsMetadata.Parameters.VoltsPerUnit = 0.0000002;%intan default
        DatsMetadata.Parameters.BitsPerSample = 16;%intan default
        DatsMetadata.Parameters.VoltageRange = 10;%not used except to make xml

        %Now that everything is in order, calculate seconds from bytes
        DatsMetadata.Recordings.Seconds = DatsMetadata.Recordings.Bytes/DatsMetadata.Parameters.NumberOfChannels/2/DatsMetadata.Parameters.SampleRate;
    case '' %if no known recording system used
        warning('No recording system found, using defaults.  May manually set values in guts of this .m file')
        DatsMetadata.RecordingSystem = 'Unknown';
        DatsMetadata.Recordings.Names = {};
        DatsMetadata.Recordings.Bytes = [];
        DatsMetadata.Recordings.Seconds = [];%can manually write in here if desired

        try
            sess = bz_getSessionInfo(basename);
        catch
            disp('No metadata or session info detected, quitting bz_DatFileMetadata')
            return
        end

        DatsMetadata.Parameters.NumberOfChannels = sess.nChannels;
        DatsMetadata.Parameters.SampleRate = sess.rates.wideband;%Lab default
        DatsMetadata.Parameters.Amplification = 1;%Intan digitized on chip, let's say 1
        DatsMetadata.Parameters.VoltsPerUnit = 0.0000002;%Intan default
        DatsMetadata.Parameters.BitsPerSample = sess.nBits;%Intan default
        DatsMetadata.Parameters.VoltageRange = 10;%not used except to make xml
end

save(fullfile(basepath,[basename '_DatsMetadata.mat']),'DatsMetadata')

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode