text:
function [anal] = er_chopTSeriesMulti3D(view,coords,scans,varargin);
% [anal] = er_chopTSeriesMulti3D(view,[roi],[scans],[options]);
%
% Chop an event-related tSeries according to specified parfiles,
% separately for each voxel in an ROI.
%
% Returns two output structs. results is a struct with
% certain summary statistics across voxels, such as the mean
% amplitudes for each condition and voxel, sems, etc. anal
% is a struct array containing one entry for each voxel. Each
% entry contains fields from er_chopTSeries (type
% help on that function for a list).
%
% Note: this does a minimal analysis only on each voxel, for
% speed purposes, unless the 'fullanal' option is passed in.
% In this case, calculations of SNR, fMRI relative amplitudes,
% and t-tests for significant activation are performed for
% each voxel, which would not otherwise be performed.
% Be warned, however, that the analysis could take twice to four
% times as long.
%
% 08/04 ras.
global dataTYPES;

if ieNotDefined('coords')
    rois = viewGet(view,'rois');
    selRoi = viewGet(view,'selectedroi');
    coords = rois(selRoi).coords;
end

dt = viewGet(view,'curdt');

if ieNotDefined('scans')
    [scans dt] = er_getScanGroup(view);
    view = viewSet(view,'curdt',dt);
end

%%%%% params/defaults %%%%%
barebones = 0;          % if 0, do full analysis; if 1, do minimal analysis
normBsl = 1;            % flag to zero baseline or not
alpha = 0.05;           % threshold for significant activations
bslPeriod = -6:0;       % period to use as baseline in t-tests, in seconds
peakPeriod = 6:12;       % period to look for peaks in t-tests, in seconds
timeWindow = -6:22;     % seconds relative to trial onset to take for each trial
onsetDelta = 0;         % # secs to shift onsets in parfiles, relative to time course
snrConds = [];          % For calculating SNR, which conditions to use (if empty, use all)
waitbarFlag = 0;        % flag to show a graphical mrvWaitbar to show load progress

%%%%% get the parfile info
trials = er_concatParfiles(view,scans);
condNums = unique(trials.cond(trials.cond > 0));
whichconds=condNums;

%%%%% parse the options %%%%%
varargin = unNestCell(varargin);
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch lower(varargin{i})
        case 'barebones', barebones = 1;
        case 'normbsl', normBsl = varargin{i+1};
        case 'alpha', alpha = varargin{i+1};
        case 'peakperiod', peakPeriod = varargin{i+1};
        case 'bslperiod', bslPeriod = varargin{i+1};
        case 'timewindow', timeWindow = varargin{i+1};
        case 'onsetdelta', onsetDelta = varargin{i+1};
        case 'snrconds', snrConds = varargin{i+1};
        case 'whichconds', whichconds = varargin{i+1};
        case 'mrvWaitbar', waitbarFlag = 1;
        otherwise, % ignore
        end
    end
end

%%%%% check if the full analysis option is set
chk = 0;
for i = 1:length(varargin)
    if isequal(lower(varargin{i}),'fullanal')
        chk = 1;
        break;
    end
end
if chk==0
    varargin{end+1} = 'barebones';
end

%%%%% get the tSeries for each voxel / concat across scans
dt = viewGet(view,'curdt');
textstring=sprintf('Loading tSeries from %s datatype scans %d-%d...',dataTYPES(dt).name,min(scans),max(scans));
h = mrvWaitbar(0,textstring);

tSeries = [];
for s = 1:length(scans)
     raw = ~(detrendFlag(view,s));
    subt = getTseriesOneROI(view,coords,scans(s),raw);
    tSeries = [tSeries; subt{1}];
    mrvWaitbar(s/length(scans),h);
end

close(h)

%%%%% run the standard choptSeries analysis on each voxel

nVoxels = size(tSeries,2);
fprintf(1,'num voxels in view %s is %5d \n',view.viewType,nVoxels);

% jl101104 peak peropd 4:12 onset delta 0
% jv062004 onset delta=-4 peakPeriod 6:16
% ras062304 onset delta=0 peakPeriod 6:14-
% kgs080304  peak period 4:14 delta=0
% kgs061504 peak period 4:14 delta  onset=0
% skipped 3 frames adjusted par files;  use mc data
% ras 041702 6:14 offset -6 blocks of 16s
% js061504  4:12 onset 0 very noisy scan did not cover FFA

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (3): Settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
defaults{1} = num2str(timeWindow);
defaults{2} = num2str(bslPeriod);
defaults{3} = num2str(peakPeriod);
defaults{4} = num2str(alpha);
defaults{5} = num2str(onsetDelta);
defaults{6} = num2str(whichconds);
defaults{7} = num2str(normBsl);
defaults{8}='y';
defaults{9}='y';
defaults{10}='n';

prompt{1} = 'Time Window, in seconds (incl. pre-stimulus onset period):';
prompt{2} = 'Baseline Period, in seconds:';
prompt{3} = 'Peak Period, in seconds:';
prompt{4} = 'Alpha for significant activation (peak vs baseline):';
prompt{5} = 'Shift onsets relative to time course, in seconds:';
prompt{6} = 'Which conditions to calculate reliability?';
prompt{7} = 'Normalize all trials during the baseline period? (0 for no, 1 for yes)';
prompt{8}=  'Plot matrix of all voxels amps [y/n]?';
prompt{9}=  'Conduct reliability analysis [y/n]?';
prompt{10}= 'Load tcUI for ROI data [y/n]?';
AddOpts.Resize = 'on';
AddOpts.Interpreter = 'tex';
AddOpts.WindowStyle = 'Normal';
answers = inputdlg(prompt,'Chop tSeries-Multi Settings...',1,defaults,AddOpts);

% exit if cancel is selected
if isempty(answers)
    return;
end

% parse the user responses / defaults
timeWindow = str2num(answers{1});
bslPeriod = str2num(answers{2});
peakPeriod = str2num(answers{3});
alpha = str2num(answers{4});
onsetDelta = str2num(answers{5});
whichconds = str2num(answers{6});
normBsl = str2num(answers{7});

if answers{8}=='y' | answers{8}=='Y'
    plotampsmat=1;
else
    plotampsmat=0;
end

if answers{9}=='y' | answers{9}=='Y'
    reliability=1;
else
    reliability=0;
end
if answers{10}=='y' | answers{10}=='Y'
    tcUI=1;
else
    tcUI=0;
end
roiName=view.ROIs(selRoi).name;

h = mrvWaitbar(0,'Chopping tSeries...');
anal = er_chopTSeries3D(tSeries,trials,roiName,...
              'peakPeriod',peakPeriod,...
              'bslPeriod',bslPeriod,...
              'timeWindow',timeWindow,...
              'alpha',alpha,...
              'onsetDelta',onsetDelta,...
              'whichconds',whichconds,...
              'normBsl',normBsl,...
              'plotampsmat',plotampsmat);
close(h)
% save data
newdir=view.viewType;
cd (newdir)
dataTypedir=getDataTypeName(view);
cd (dataTypedir)

if ~exist('ROI4D')
        mkdir 'ROI4D'
end
cd ..    % newdir
cd ..    % data dir

anal.roiName=view.ROIs(selRoi).name;
whichconds
if reliability
    [rc,pc,rcov, pcov, odd, even]=pattern_analysis(anal.allamps(:,whichconds,:),anal.allampsminusDC(:,whichconds,:),anal.labels(whichconds),anal.roiName);
    anal.pc=pc;
    anal.rc=rc;
    anal.rcov=rcov;
    anal.pcov=pcov;
end

% too long but makes it platform insensitive
cd (newdir)
cd (dataTypedir)
cd ('ROI4D')
filename=view.ROIs(selRoi).name;
eval (['save ' filename ' anal']);
fprintf(1, 'saved anal to %s/%s/ROI4D/%s.mat \n ', newdir,dataTypedir,filename);
cd ..
cd ..
cd ..
if tcUI
    fprintf(1, 'Loading mean ROI timecourse UI..........')
    tc = timeCourseUI(view,view.ROIs(selRoi).name,scans,dt);
end
return


text:
classdef TimeMagnitudePlotter < TimeSomethingPlotter
    % Used to create time-mag plots
    %
    % plotting multiple times will simply replace the existing plot
    %
    % Example:
    %   TimeMagnitudePlotter.plot(catalog,bigCat)
    %   ax = TimeMagnitudePlotter.axes;
    %   TimeMagnitudePlotter.clear();

    %TODO add ability to plot "big" events

    methods
        function obj=TimeMagnitudePlotter()
            obj@TimeSomethingPlotter('Magnitude','');
        end

        function pl=plot(obj,ax,catalog,bigcat)
            % plot plot a time-mag series for this catalog, with symbol sizes representing
            % event size
            % pl = plot(catalog)
            %
            if ~exist('ax','var') || isempty(ax)
                % ax = findobj('Tag','time_mag_axis');
                f=figure('Name','Time-Magnitude Plot','NumberTitle','off', ...
                    'Tag',obj.Tags.Figure);
                addAboutMenuItem();
                obj.ax=axes(f);
            else
                obj.ax=ax;
            end
            if isempty(obj.ax.Tag)
                obj.ax.Tag=obj.Tags.Axes;
            end

            obj.ax.Visible = 'off';

            % plotting from ZERO magnitude is arbitrary, and stemplot becomes unusable if all magnitudes
            % are below zero. Therefore, make sure stems are always going up
            minMag= min(catalog.Magnitude);
            baseValue = min([0 , floor(minMag)]);

            if obj.hasLotsOfEvents(catalog)
                pl=obj.scatter(catalog);
            else
                pl=obj.stem(catalog, 'BaseValue', baseValue);
            end

            obj.prepare_axes(catalog);

            if exist('bigcat','var')
                obj.overlayBigEvents(bigcat);
            end
            obj.ax.Visible = 'on';
        end

    end
end

author: CelsoReyes


text:
function cS = chromite

cs_Chr = crystalSymmetry('m3m');
N = Miller({1,1,0},cs_Chr);
dist = [1];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
function X = imag(X)
%IMAG (overloaded)

X.basis = imag(X.basis);
X = clean(X);
if isa(X,'sdpvar')
    X.conicinfo = [0 0];

    if length(X.midfactors)>0
        newleftfactors = {};
        newmidfactors = {};
        newrightfactors = {};
        for i = 1:length(X.midfactors)
            ar = real(X.leftfactors{i});
            ai = imag(X.leftfactors{i});
            br = real(X.rightfactors{i});
            bi = imag(X.rightfactors{i});

            if nnz(ai)>0 & nnz(br)>0
                newleftfactors{end+1} = ai;
                newmidfactors{end+1} = X.midfactors{i};
                newrightfactors{end+1} = br;
            end

            if nnz(ar)>0 & nnz(bi)>0
                newleftfactors{end+1} = ar;
                newmidfactors{end+1} = X.midfactors{i};
                newrightfactors{end+1} = bi;
            end

        end
        X.leftfactors = newleftfactors;
        X.midfactors = newmidfactors;
        X.rightfactors = newrightfactors;
    end
end

%
% Y = X;
% x_lmi_variables = X.lmi_variables;
% lmi_variables = [];
% n = X.n;
% m = X.m;
% imagX = imag(X.basis(:,1));
% Y.basis = imagX(:);
%
% j = 1;
% for i = 1:length(x_lmi_variables)
%     imagX = imag(X.basis(:,i+1));
%     if (norm(imagX,inf)>0)
%         Y.basis(:,j+1) = imagX(:);
%         lmi_variables = [lmi_variables x_lmi_variables(i)];
%         j = j+1;
%     end
% end
% if isempty(lmi_variables)
%     Y = full(reshape(Y.basis,n,m));
% else
%     Y.lmi_variables = lmi_variables;
%     % Reset info about conic terms
%     Y.conicinfo = [0 0];
% end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% AK 20160415 Clustering Play

%clearvars -except M_0_fluo;

% Export M0_fluo all cells, around 40k
f_cells = M_0_fluo;
num_cells = length(f_cells);

num_select = 2000;
clear idx_select;
idx_select = randperm(num_cells, num_select);
f_cells_select = f_cells(idx_select,:);

%% Correlation Distance
tic
Dpairs = pdist(f_cells_select,'correlation');
toc

figure;
h = histogram(1-Dpairs);%Plot corr, not corr dist
title('Correlation Distance');

%% Euclidean Distance

% tic
% Dpairs = pdist(f_cells_select,'euclidean');
% toc
%
% figure;
% histogram(Dpairs);
% title('Euclidean Distance');

%% Fit random variance
% Fit Gaussian to part above 1
%h = histogram(1-Dpairs);
idx = find(h.BinEdges < 0.1 & h.BinEdges > -0.1);
allBins = h.BinEdges(1:end-1) + h.BinWidth/2;
allCorrVals = h.Values;
unCorrBins = h.BinEdges(idx) + h.BinWidth/2;
unCorrVals = h.Values(idx+1);

figure;
options = fitoptions('gauss1', 'Lower', [-Inf 1 -Inf], 'Upper', [Inf 1 Inf]);
% Forcing b1 = 1 doesn't work well
%f = fit(unCorrBins',unCorrVals','gauss1',options);
f = fit(unCorrBins',unCorrVals','gauss1');
plot(f,allBins,allCorrVals);
%plot(f,unCorrBins,unCorrVals);
hold on;
grid on;
title('Pairwise Cell Correlations');
xlabel('Correlation')
ylabel('Frequency')

%% Cumulative Significant Correlations

fracCorrSig = (allCorrVals'-f(allBins));
%figure;plot(allBins,fracCorrSig);
cumCorrSig = cumsum(fracCorrSig);
figure;plot(allBins,cumCorrSig/cumCorrSig(end));
grid on;

author: xiuyechen

repo: FishExplorer


text:
function suppData_GatherResultsRF(pathRF,pathSUPP)

startpath = pwd;

% INITIALIZATION
tables = struct;
% nameOutcomes = {'Locoregional','Distant','Death','Death','DeathSign'};
% fSetNames = {'PET','CT','PETCT','CT','CT'};
nameOutcomes = {'Locoregional','Distant','Death'};
fSetNames = {'PET','CT','CT'};
metrics = {'AUC','Sensitivity','Specificity','Accuracy','CI'}; nMetrics = numel(metrics);
strMetrics = [];
for m = 1:nMetrics
    strMetrics = [strMetrics,metrics{m},','];
end
strMetrics(end) = [];
paths = {pathSUPP,pathSUPP,pathRF}; nComp = numel(paths);
prefixRead = {'radiomics','clinical',''};
fields = {'radiomics','clinical','combined'};

% PROCESSING
for c = 1:nComp
    if strcmp(fields{c},'clinical')
        nExp = 3; % 3 clinical outcomes
        rowNames = cell(nExp,1);
        for o = 1:nExp
            rowNames{o} = nameOutcomes{o};
        end
    else
        nExp = numel(nameOutcomes);
        rowNames = cell(nExp,1);
        for o = 1:nExp
            rowNames{o} = [nameOutcomes{o},'_',fSetNames{o}];
        end
    end
    cd(paths{c})
    for m = 1:nMetrics
        eval([metrics{m},' = zeros(nExp,1);']);
    end
    for o = 1:nExp
        if strcmp(fields{c},'clinical') || strcmp(fields{c},'volClinical')
            load(['testResultsRF',prefixRead{c},'_',nameOutcomes{o}]) % results gets out of there
        else
            load(['testResultsRF',prefixRead{c},'_',fSetNames{o},'_',nameOutcomes{o}]) % results gets out of there
        end
        for m = 1:nMetrics
            eval([metrics{m},'(o,1) = results.',metrics{m},';']);
        end
    end
    eval(['tableExp = table(',strMetrics,',''RowNames'',rowNames);']);
    tables.(fields{c}) = tableExp;
end

cd(pathSUPP), save('summaryTables','tables')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [data_train, labels_train, data_devel, labels_devel, raw_devel, PC, means_norm, stds_norm, devel_ids, devel_success] = ...
    Prepare_HOG_AU_data(train_users, devel_users, au_train, rest_aus, Bosphorus_dir, params_data_dir)

%%
addpath(genpath('../data extraction/'));

% First extracting the labels
[ labels_train, valid_ids_train, filenames ] = extract_Bosphorus_labels(Bosphorus_dir, train_users, au_train);

[ labels_other, ~, ~ ] = extract_Bosphorus_labels(Bosphorus_dir, train_users, rest_aus);

% Reading in the HOG data (of only relevant frames)
[train_appearance_data, valid_ids_train_hog, vid_ids_train_string] = Read_HOG_files(filenames, params_data_dir);

[train_geom_data] = Read_geom_files(filenames,  params_data_dir);

% Subsample the data to rebalance it
if(numel(train_users) > 0)
    reduced_inds = false(size(labels_train,1),1);
    reduced_inds(labels_train > 0) = true;

    % make sure the same number of positive and negative samples is taken
    pos_count = sum(labels_train > 0);
    neg_count = sum(labels_train == 0);

    num_other = floor( pos_count / (size(labels_other, 2)));

    inds_all = 1:size(labels_train,1);

    for i=1:size(labels_other, 2)+1

        if(i > size(labels_other, 2))
            % fill the rest with a proportion of neutral
            inds_other = inds_all(sum(labels_other,2)==0 & ~labels_train );
                num_other_i = min(numel(inds_other), pos_count - sum(labels_train(reduced_inds,:)==0));
        else
                % take a proportion of each other AU
            inds_other = inds_all(labels_other(:, i) & ~labels_train );
            num_other_i = min(numel(inds_other), num_other);
        end
        inds_other_to_keep = inds_other(round(linspace(1, numel(inds_other), num_other_i)));
        reduced_inds(inds_other_to_keep) = true;

    end

    % Remove invalid ids based on CLM failing or AU not being labelled
    reduced_inds(~valid_ids_train) = false;
    reduced_inds(~valid_ids_train_hog) = false;

    labels_other = labels_other(reduced_inds, :);
    labels_train = labels_train(reduced_inds,:);
    train_appearance_data = train_appearance_data(reduced_inds,:);
    train_geom_data = train_geom_data(reduced_inds,:);
    vid_ids_train_string = vid_ids_train_string(reduced_inds,:);
end
%% Extract devel data

% First extracting the labels
[ labels_devel, valid_ids_devel, filenames_devel ] = extract_Bosphorus_labels(Bosphorus_dir, devel_users, au_train);

% Reading in the HOG data (of only relevant frames)
[devel_appearance_data, valid_ids_devel_hog, vid_ids_devel_string] = Read_HOG_files(filenames_devel, params_data_dir);
devel_success = valid_ids_devel_hog;
devel_ids = vid_ids_devel_string;

[devel_geom_data] = Read_geom_files(filenames_devel, params_data_dir);

% Peforming zone specific masking
if(au_train < 8 || au_train == 43 || au_train == 45) % upper face AUs ignore bottom face
    % normalise the data
    pca_file = '../../pca_generation/generic_face_upper.mat';
    load(pca_file);
elseif(au_train > 9) % lower face AUs ignore upper face and the sides
    % normalise the data
    pca_file = '../../pca_generation/generic_face_lower.mat';
    load(pca_file);
elseif(au_train == 9) % Central face model
    % normalise the data
    pca_file = '../../pca_generation/generic_face_rigid.mat';
    load(pca_file);
end

% Grab all data for validation as want good params for all the data
raw_devel = cat(2, devel_appearance_data, devel_geom_data);

devel_appearance_data = bsxfun(@times, bsxfun(@plus, devel_appearance_data, -means_norm), 1./stds_norm);

data_devel = devel_appearance_data * PC;

data_devel = cat(2, data_devel, devel_geom_data);

valid_ids_devel = valid_ids_devel & devel_success;
data_devel = data_devel(valid_ids_devel,:);
labels_devel = labels_devel(valid_ids_devel,:);
raw_devel = raw_devel(valid_ids_devel,:);
devel_success = devel_success(valid_ids_devel,:);
devel_ids = devel_ids(valid_ids_devel);

if(numel(train_users) > 0)
    train_appearance_data = bsxfun(@times, bsxfun(@plus, train_appearance_data, -means_norm), 1./stds_norm);

    data_train = train_appearance_data * PC;
    data_train = cat(2, data_train, train_geom_data);
else
    data_train = [];
end

geom_size = max(size(train_geom_data, 2), size(devel_geom_data, 2));

PC_n = zeros(size(PC)+geom_size);
PC_n(1:size(PC,1), 1:size(PC,2)) = PC;
PC_n(size(PC,1)+1:end, size(PC,2)+1:end) = eye(geom_size);
PC = PC_n;

means_norm = cat(2, means_norm, zeros(1, geom_size));
stds_norm = cat(2, stds_norm, ones(1, geom_size));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function tapas_hgf_jget_plotTraj(r)
% Plots the estimated trajectories for the HGF perceptual model for
% the JGET project
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_hgf_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Time axis
if size(r.u,2) > 1 && ~isempty(find(strcmp(fieldnames(r.c_prc),'irregular_intervals'))) && r.c_prc.irregular_intervals
    t = r.u(:,end)';
else
    t = ones(1,size(r.u,1));
end

ts = cumsum(t);
ts = [0, ts];

% Do we know the generative parameters?
if size(r.u,2) > 2
    genpar = true;
    mean   = r.u(:,2);
    sd     = r.u(:,3);
else
    genpar = false;
end

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = length(r.p_prc.p)/8;
end

% Upper levels
for j = 1:l-1

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Left subplot (x)                       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(l+1,2,2*j-1);

    if plotsd == true
        upperprior = r.p_prc.mux_0(l-j+1) +1.96*sqrt(r.p_prc.sax_0(l-j+1));
        lowerprior = r.p_prc.mux_0(l-j+1) -1.96*sqrt(r.p_prc.sax_0(l-j+1));
        upper = [upperprior; r.traj.mux(:,l-j+1)+1.96*sqrt(r.traj.sax(:,l-j+1))];
        lower = [lowerprior; r.traj.mux(:,l-j+1)-1.96*sqrt(r.traj.sax(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mux_0(l-j+1); r.traj.mux(:,l-j+1)], 'b', 'LineWidth', 1.5);
    hold all;
    plot(0, r.p_prc.mux_0(l-j+1), 'ob', 'LineWidth', 1.5); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu x_', num2str(l-j+1)]);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Right subplot (alpha)                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(l+1,2,2*j);

    if plotsd == true
        upperprior = r.p_prc.mua_0(l-j+1) +1.96*sqrt(r.p_prc.saa_0(l-j+1));
        lowerprior = r.p_prc.mua_0(l-j+1) -1.96*sqrt(r.p_prc.saa_0(l-j+1));
        upper = [upperprior; r.traj.mua(:,l-j+1)+1.96*sqrt(r.traj.saa(:,l-j+1))];
        lower = [lowerprior; r.traj.mua(:,l-j+1)-1.96*sqrt(r.traj.saa(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mua_0(l-j+1); r.traj.mua(:,l-j+1)], 'b', 'LineWidth', 1.5);
    hold all;
    plot(0, r.p_prc.mua_0(l-j+1), 'ob', 'LineWidth', 1.5); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of \alpha_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu \alpha_', num2str(l-j+1)]);
end

% Input level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Left subplot (x)                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(l+1,2,2*l-1);

if plotsd == true
    upperprior = r.p_prc.mux_0(1) +1.96*sqrt(r.p_prc.sax_0(1));
    lowerprior = r.p_prc.mux_0(1) -1.96*sqrt(r.p_prc.sax_0(1));
    upper = [upperprior; r.traj.mux(:,1)+1.96*sqrt(r.traj.sax(:,1))];
    lower = [lowerprior; r.traj.mux(:,1)-1.96*sqrt(r.traj.sax(:,1))];

    plot(0, upperprior, 'or', 'LineWidth', 1);
    hold all;
    plot(0, lowerprior, 'or', 'LineWidth', 1);
    fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
plot(ts, [r.p_prc.mux_0(1); r.traj.mux(:,1)], 'r', 'LineWidth', 1.5);
hold all;
plot(0, r.p_prc.mux_0(1), 'or', 'LineWidth', 1.5); % prior
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0.6 0]); % inputs
if genpar
    plot(ts(2:end), mean, '-', 'Color', 'k', 'LineWidth', 1); % mean of input distribution
    plot(ts(2:end), mean +1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
    plot(ts(2:end), mean -1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
end
xlim([0 ts(end)]);
title(['Input u (green) and posterior expectation of x_1 (red) for \kappa_x=', ...
       num2str(r.p_prc.kax), ', \omega_x=', num2str(r.p_prc.omx)], 'FontWeight', 'bold');
ylabel('u, \mu x_1');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Right subplot (alpha)                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(l+1,2,2*l);

if plotsd == true
    upperprior = r.p_prc.mua_0(1) +1.96*sqrt(r.p_prc.saa_0(1));
    lowerprior = r.p_prc.mua_0(1) -1.96*sqrt(r.p_prc.saa_0(1));
    upper = [upperprior; r.traj.mua(:,1)+1.96*sqrt(r.traj.saa(:,1))];
    lower = [lowerprior; r.traj.mua(:,1)-1.96*sqrt(r.traj.saa(:,1))];

    transupperprior = sqrt(exp(r.p_prc.kau *upperprior +r.p_prc.omu));
    translowerprior = sqrt(exp(r.p_prc.kau *lowerprior +r.p_prc.omu));
    transupper = sqrt(exp(r.p_prc.kau *upper +r.p_prc.omu));
    translower = sqrt(exp(r.p_prc.kau *lower +r.p_prc.omu));

    plot(0, transupperprior, 'or', 'LineWidth', 1);
    hold all;
    plot(0, translowerprior, 'or', 'LineWidth', 1);
    fill([ts, fliplr(ts)], [(transupper)', fliplr((translower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
transmuaprior = sqrt(exp(r.p_prc.kau *r.p_prc.mua_0(1) +r.p_prc.omu));
plot(ts, [transmuaprior; sqrt(exp(r.p_prc.kau *r.traj.mua(:,1) +r.p_prc.omu))], 'r', 'LineWidth', 1.5);
hold all;
plot(0, transmuaprior, 'or', 'LineWidth', 1.5); % prior
if genpar
    plot(ts(2:end), sd, '--', 'Color', 'k', 'LineWidth', 1);
end
xlim([0 ts(end)]);
title(['Belief on noise (red) for \kappa_\alpha=', ...
       num2str(r.p_prc.kaa), ', \omega_\alpha=', num2str(r.p_prc.oma)], 'FontWeight', 'bold');
ylabel('\mu \alpha_1');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Decision model                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(l+1,2,2*l+1);

if plotsd == true
    upper = r.traj.muxhat(:,1)+1.96*sqrt(r.p_obs.ze +r.traj.saxhat(:,1));
    lower = r.traj.muxhat(:,1)-1.96*sqrt(r.p_obs.ze +r.traj.saxhat(:,1));

    fill([ts(2:end), fliplr(ts(2:end))], [(upper)', fliplr((lower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    hold all;
end
plot(ts(2:end), r.traj.muxhat(:,1), 'Color', [153/256 17/256 153/256], 'LineWidth', 1.5);
hold all;
plot(ts(2:end), r.y(:,1), '.', 'Color', [1 0.65 0], 'MarkerSize', 15); % responses
if genpar
    plot(ts(2:end), mean, '-', 'Color', 'k', 'LineWidth', 1); % mean of input distribution
    plot(ts(2:end), mean +1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
    plot(ts(2:end), mean -1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1); % 95% interval of input distribution
end
xlim([1 ts(end)]);
title('Decision model: prediction of decision (purple) and decision (orange)', 'FontWeight', 'bold');
ylabel('y, \^{\mu} x_1');
xlabel({'Trial number', ' '}); % A hack to get the relative subplot sizes right
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Learning rate                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if genpar
subplot(l+1,2,2*l+2);
[AX, H1, H2 ] = plotyy(ts(2:end), mean, ts(2:end), r.traj.lrx(:,1));
hold all;
ylim(AX(1), [min(mean -1.96.*sd-3), max(mean +1.96.*sd+3)]);
plot(AX(1), ts(2:end), mean +1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1.1); % 95% interval of input distribution
plot(AX(1), ts(2:end), mean -1.96.*sd, '--', 'Color', 'k', 'LineWidth', 1.1); % 95% interval of input distribution
set(H1, 'Color', 'k', 'LineWidth', 1.1);
set(H2, 'Color', [178/256, 34/256, 34/256], 'LineWidth', 1.5);
set(AX(1), 'YColor', 'k');
set(AX(2), 'YColor', 'k');
xlim(AX(1), [1 ts(end)]);
xlim(AX(2), [1 ts(end)]);
title('Learning rate (bordeaux) and input sampling distribution (black)', 'FontWeight', 'bold');
ylabel(AX(1), 'Input');
ylabel(AX(2), '\sigma_x');
xlabel('Trial number');
hold off;
end

author: translationalneuromodeling

repo: tapas


text:
function [el, lab] = read_elec(fn)

% READ_ELEC reads "la/mu" electrode parameters from a MBF electrode file
% which are used to position them on a triangulated surface
%
% [el, lab] = read_elec(filename)
%
% where el = [tri, la, mu]
% and lab contains the electrode labels (if present)
%
% See also READ_TRI, TRANSFER_ELEC

% Copyright (C) 1998, Robert Oostenveld
%
% 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$

fid = fopen_or_error(fn, 'rt');

% read the number of electrodes
Nel = sscanf(fgetl(fid), '%d');

% read the electrode triangle, lambda and mu
for i=1:Nel
  str = fgetl(fid);
  el(i,:)  = sscanf(str, '%f %f %f')';
  indx = find(str=='!');
  if (indx)
    lab(i,:) = sprintf('%6s', str((indx+1):length(str)));
  end
end
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function sp=orthpoly(n,opfamily,alpha,beta)
% orthpoly: generate an n'th order (symbolic) orthogonal polynomial
% usage: sp=orthpoly(n,opfamily,alpha,beta);
%
% arguments: (input)
%  n        - order of polynomial requested
%
%  opfamily - (OPTIONAL) character string - specifies the
%             type of orthogonal polynomial to be generated.
%             Legal polytypes are: {'legendre', 'hermite',
%             '1cheby', '2cheby', 'laguerre', 'jacoby'}
%
%             Default == 'legendre'
%
%             Note: 1cheby refers to first kind
%             chebychev polynomials (Tn(x)),
%             and 2cheby refers to second kind
%             chebychev polynomials (Un(x))
%
%             The domain of support is [-1,1] for the polynomial
%             families 'legendre', '1cheby', '2cheby', 'jacoby'.
%             It is [0,inf] for 'laguerre', and [-inf,inf] for
%             the 'hermite' family.
%
%             See Abramowitz & Stegun for in-depth information
%             on all of these polynomial families.
%
%  alpha,beta - scalar - orthogonal family parameters -
%             See Abramowitz & Stegun for more information.
%
%             These parameters are ignored for 'legendre',
%             'hermite', '1cheby' and '2cheby' polys.
%
% arguments: (output)
%  p        - sympoly object containing the polynomial

% defaults
if (nargin<2)|isempty(opfamily)
  % default is legendre
  opfamily = 'legendre';
else
  if ~ischar(opfamily)
    error 'Opfamily must be a character string'
  end
  valtypes = {'legendre', 'hermite', 'cheby1', '1chebychev' ...
            'cheby2', '2chebychev', 'laguerre', 'jacobi'};

  ptind = strmatch(opfamily,valtypes);
  if isempty(ptind)
    % no match
    error 'Invalid opfamily'
  elseif length(ptind)>1
    % match more than one
    error 'Ambiguous opfamily'
  else
    % there was exactly one match
    opfamily = valtypes{ptind};
  end
end

if nargin<1
  error 'Polynomial order must be supplied. No default for n.'
end

if (nargin<3) || isempty(alpha)
  % default is 0
  alpha = 0;
end

if (nargin<4) || isempty(beta)
  % default is 0
  beta = 0;
end

% initialize (-1)'th and zero'th order sympolys
% for the three term recurrence relation.
pnm1=sympoly(0);
pn=sympoly(1);
x=sympoly('x');

if n==0
  % p0 is easy to get
  sp=pn;
  return
elseif n<0
  error 'Polynomial order must not be negative'
end

% iterate to get polynomial
for i=0:(n-1)
  switch opfamily
    case 'legendre'
      pnp1=((2*i+1)*x*pn - i*pnm1)./(i+1);
      % normalize so that Tn(1)=1
      pnp1=pnp1./double(subs(pnp1,'x',1));

    case 'hermite'
      % neither alpha nor beta is meaningful
      pnp1=2*pn*x - 2*i*pnm1;

    case 'laguerre'
      % only alpha is meaningful here
      pnp1=((2*i+alpha+1)*pn - pn*x - (i+alpha)*pnm1)./(i+1);

    case {'cheby1' '1chebychev'}
      % first kind chebychev
      % neither alpha nor beta is meaningful
      pnp1=(2*pn*x - pnm1);
      % normalize so that Tn(1)=1
      pnp1=pnp1/double(subs(pnp1,'x',1));

    case {'cheby2' '2chebychev'}
      % second kind chebychev
      % neither alpha nor beta is meaningful
      pnp1=(2*pn*x - pnm1);

    case 'jacobi'
      % both alpha and beta are needed
      if (alpha~=0) | (beta~=0)
        a1n=2*(i+1)*(i+alpha+beta+1)*(2*i+alpha+beta);
        a2n=(2*i+alpha+beta+1)*(alpha^2-beta^2);
        if (2*i+alpha+beta)<=150
          a3n=gamma(2*i+alpha+beta+3)./gamma(2*i+alpha+beta);
        else
          a3n=exp(gammaln(2*i+alpha+beta+3)-gammaln(2*i+alpha+beta));
        end
        a4n=2*(i+alpha)*(i+beta)*(2*i+alpha+beta+2);
        pnp1=(a2n*pn + a3n*pn*x - a4n*pnm1)/a1n;
      else
        pnp1=((2*i+1)*pn*x - i*pnm1)/(i+1);
      end
  end
  pnm1=pn;
  pn=pnp1;

end

% return the generated sympoly
sp=pnp1;

author: Sable

repo: mcbench-benchmarks


text:
function xform = find_xform(c, d)
%Finds the xform that maps from c to d
%in this case the xform is a simple translation and scaling
%xform is a [3 x 3] matrix which maps coordinates from c's frame to
%d's frame.  the transformation is applied to homogeneous
%coordinaes

%convert bounding box to cornners
xs(:,1) = c([1 2])';
xs(:,2) = c([3 2])';
xs(:,3) = c([1 4])';
xs(:,4) = c([3 4])';

ys(:,1) = d([1 2])';
ys(:,2) = d([3 2])';
ys(:,3) = d([1 4])';
ys(:,4) = d([3 4])';

xs(3,:) = 1;
ys(3,:) = 1;

A=ys*pinv(xs);
A(abs(A)<.000001) = 0;

xform = A;

author: quantombone

repo: exemplarsvm


text:
function diceC = ellipseDice(E1,E2,varargin)
% Compute the percentage overlap of a pair of ellipses
%
% Syntax
%   diceC = ellipseDice(E1,E2,varargin)
%
% Input
%   E1, E2 - Two structs that contain the ellipse parameters
%       .center
%       .sigma
%       .theta
%
% Optional key/value pairs
%   spatial samples - vector of the sample points in deg (default -10:0.1:10);
%   show -  Create an image showing the overlap
%
% Outputs
%   diceC - Dice coefficient = 2*areaOfE1E2Intersection / eArea1 + eArea2 ;
%
% Author, Wandell January 17, 2020
%
% See also
%   ellipseInterior, ellipsePoints, ...

% Examples:
%{
   E1.center = [0,0]; E1.sigma = [3,1]; E1.theta = pi;
   E2.center = [-1,0]; E2.sigma = [3,1]; E2.theta = pi/2;
   diceC = ellipseDice(E1,E2,'show',true);
%}
%{
   E1.center = [0,0]; E1.sigma = [1,1]; E1.theta = pi;
   E2.center = [-1,0]; E2.sigma = [3,3]; E2.theta = pi/2;
   diceC = ellipseDice(E1,E2,'show',true);
%}
%{
   E1.center = [0,0]; E1.sigma = [1,1]; E1.theta = pi;
   E2.center = [0,0]; E2.sigma = [1,1]; E2.theta = pi/2;
   diceC = ellipseDice(E1,E2,'show',true);
%}

%% Input parameters

% Remove spaces and force lower case
varargin = mrvParamFormat(varargin);

p = inputParser;

% Validation function
vFunc = @(x)(isstruct(x) && isfield(x,'center') && isfield(x,'sigma') && isfield(x,'theta'));
p.addRequired('E1',vFunc);
p.addRequired('E2',vFunc);

p.addParameter('spatialsamples',(-10:0.05:10),@isvector);
p.addParameter('show',false,@islogical);

p.parse(E1,E2,varargin{:});

samples = p.Results.spatialsamples;
show    = p.Results.show;

%% Compute

% 1 inside, 0 outside.
[img1, nSamples] = ellipseInterior('center',E1.center, ...
    'sigma',E1.sigma, ...
    'theta',E1.theta',...
    'spatial samples',samples);
img2   = ellipseInterior('center',E2.center,...
    'sigma',E2.sigma,...
    'theta',E2.theta,...
    'spatial samples',samples);

%%
overlapArea = sum(img1(:) .* img2(:));
diceC = 2*overlapArea/ (sum(img1(:)) + sum(img2(:)));

if show
    img = img1 + img2;
    mrvNewGraphWin; colormap([0.2 0.3 0.4; 0.6 0.6 0.6; 1 1 1]);
    image(samples,samples,img + 1); axis image;
    xlabel('Deg'); ylabel('Deg'); grid on
end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function b_nm = beamWeightsFromFunction(fHandleArray, order)
%BEAMWEIGHTSFROMFUNCTION Generate spherical beamweights for arbitrary pattern
%
%   This function computes the SH coefficients of an arbitrary target
%   beampattern, by performing the SHT on it. The pattern should be
%   specified as a function handle that expects a Kx1 vector of azimuths as
%   first argument, and a vector of Kx1 elevations as the second argument,
%   and return the functions values for these K directions.
%   Multiple patterns can be passed as an array of function handles.
%   The order specifies the maximum SH order of the beamweights.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BEAMWEIGHTSFROMFUNCTION.M - 10/4/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get a dense uniform design
[~, dirsAziElev] = getTdesign(20);
% convert azi-elev to azi-incl
dirsAziInc = [dirsAziElev(:,1) pi/2-dirsAziElev(:,2)];

Nfunc = length(fHandleArray);
b_nm = zeros((order+1)^2, Nfunc);
for nf =1:Nfunc
    if Nfunc==1, fPattern = fHandleArray;
    else fPattern = fHandleArray{nf};
    end

    % compute the pattern at grid points
    F_dirs = fPattern(dirsAziElev(:,1), dirsAziElev(:,2));

    % get harmonic coefficients
    b_nm(:,nf) = directSHT(order, F_dirs, dirsAziInc, 'real');
end

author: polarch

repo: Spherical-Array-Processing


text:
function destination = BV_CombineWithMask(sourceFile, destinationFile, maskFile, outFile, normalizeImages)
%
%   destination = BV_CombineWithMask(sourceFile, destinationFile, maskFile, outFile)
%
% Combines Analyze images using a mask
%
% Like this:
%    destination(mask) = source(mask);
%
% EXAMPLE:
%       BV_CombineWithMask(homog_brain, orig, brain_mask);
%
% $Author: bob $
% $Date: 2003/12/20 00:01:50 $

if (~exist('normalizeImages','var'))
    normalizeImages=0;
end
if ~exist('outFile','var')
    outFile = [destinationFile,'_Combined']
end

% Load in the images:
[source,mmPerVox]      = loadAnalyze(sourceFile,'',1);
[destination,mmPerVox] = loadAnalyze(destinationFile,'',1);

% Assume that mask is 0s and 1s
mask = loadAnalyze(maskFile,'',1);
mask = mask~=0;

% Do some normalization if required
if (normalizeImages);
    disp('Normalizing...');
    source = double(source);
    source = source-min(source(:));
    source = source/max(source(:))*65000;
    source = uint16(source);
    destination = double(destination);
    destination = destination-min(destination(:));
    destination = destination/max(destination(:))*65000;
    destination = uint16(destination);
end

% Take the locations in the Source that are in the mask positions
% and copy them into the Destination file.
destination(mask) = source(mask);

if nargout == 0
    saveAnalyze(double(destination),outFile, mmPerVox);
    fprintf('Saving combined file:  %s\n',outFile);
end

return;


text:
% I = PATCH_INPAINT(Iorg, Mask, verbose) inpaints an image
%
%   Iorg  -  original image
%   Mask  -  replacement mask; pixels where mask==1 will be inpainted
%
%   I     - inpainted image
function A = patch_inpaint(Aorg,Morg, verbose, sigma)

if nargin < 3
    verbose = true;
end
if nargin < 4
    sigma = 0.25;
end

if ~exist('CSH_nn.m','file')
    error('CSH_nn not found. Please download and add to path:  http://www.eng.tau.ac.il/~simonk/CSH/index.html');
end

width = 8;
csh_iterations = 5;
k = 1;
calcBnn = 0;

diffthresh = 10;

% Also accept RGB masks, but only use first channel
if size(Morg,3) > 1
    Morg = Morg(:,:,1);
end

% Determinte starting scale
[m n rgb] = size(Aorg);
startscale = -ceil(log2(min(m,n))) + 5;
scale = 2^(startscale);

% Resize image to starting scale
A = imresize(Aorg,scale);
M = imresize(Morg,scale);
M(M>0)=1;
M3 = repmat(M,[1 1 3])==1;

% Random starting guess for inpainted image
[m n rgb] = size(A);
Rnd = uint8(255*rand(m,n,rgb));
A(M3) = Rnd(M3);

% Go through all scales
for logscale = startscale:0

    scale = 2^(logscale);

    % Maximum number of iterations on this scale;
    % oscillations are possible
    iterations = 30;

    if verbose
        fprintf('Scale = 2^%d\n',logscale);
    end

    for iter = 1:iterations
        if verbose
            fprintf('  Iteration %2d/%2d',iter,iterations);
            imshow(A);
            pause(0.001)
        end

        B = A;
        B(M3)=0;

        % Compute NN field
        CSH_ann = CSH_nn(A,B,width,csh_iterations,k,calcBnn,M);

        % Now be work in double precision
        A = double(A)./255;

        % Create new image by letting each patch vote
        R = zeros(size(A));
        Rcount = zeros(m,n);
        for i = 1:m
            for j = 1:n
                if M(i,j) == 1
                    if 1<=i && i+width-1<=m && 1<=j && j+width-1<=n
                        patch = A(i:i+width-1,j:j+width-1,:);
                        i2 = CSH_ann(i,j,2);
                        j2 = CSH_ann(i,j,1);
                        patch2 = A(i2:i2+width-1,j2:j2+width-1,:);

                        d = sum( (patch(:)-patch2(:)).^2 );
                        sim = exp( -d / (2*sigma^2) );

                        pi = i:i+width-1;
                        pj = j:j+width-1;
                        R(pi,pj,:) = R(pi,pj,:) + sim*patch2;
                        Rcount(pi,pj) = Rcount(pi,pj) + sim;
                    end
                end
            end
        end

        % Normalize and
        Rcount = repmat(Rcount,[1 1 3]);
        R(Rcount>0) = R(Rcount>0) ./ Rcount(Rcount>0);
        % Keep pixels outside mask
        R(~M3)=A(~M3);
        % Convert back to uint8
        Aprev = 255*A;
        A = uint8(255*R);

        if iter>1
            % Measure how much image has changed
            diff = sum( (double(A(:))-double(Aprev(:))).^2 ) / sum(M(:)>0);
            if verbose
                fprintf(' diff = %f\n',diff);
            end
            % Stop iterating if change is low
            if diff < diffthresh
                break;
            end
        elseif verbose
            fprintf('\n');
        end
    end

    %Upsample A for the next scale
    if logscale < 0
        Adata = imresize(Aorg,2*scale);
        [m n rgb] = size(Adata);
        A = imresize(A,[m n]);

        M = imresize(Morg,[m n]);
        M(M>0)=1;
        M3 = repmat(M,[1 1 3])==1;

        % Outside mask, A is equal to original image
        A(~M3) = Adata(~M3);
    end
end


text: d

author: jianxiongxiao

repo: ProfXkit


text:
function pass = test_biharm(pref)
% Check the CHEBFUN2 BIHARM command.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e6*pref.cheb2Prefs.chebfun2eps;

% Function to be used:
ff = @(x,y) x.^2.*y.^2;

% Bihamrmonic operator applied to ff:
f = chebfun2(ff);
fB = biharm(f);

% Exact solution:
gg = @(x,y) 8;
g = chebfun2(gg);

% Compare:
pass = ( norm(fB - g) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
filename='Cantilever_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function linpack_z_test ( )

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

  linpack_z_test01 ( );
  linpack_z_test02 ( );
  linpack_z_test03 ( );
  linpack_z_test04 ( );
  linpack_z_test05 ( );
  linpack_z_test06 ( );
  linpack_z_test07 ( );
  linpack_z_test08 ( );
  linpack_z_test09 ( );

  linpack_z_test10 ( );
  linpack_z_test11 ( );
  linpack_z_test12 ( );
  linpack_z_test13 ( );
  linpack_z_test14 ( );
  linpack_z_test15 ( );
  linpack_z_test16 ( );
  linpack_z_test17 ( );
  linpack_z_test18 ( );
  linpack_z_test19 ( );

  linpack_z_test20 ( );
  linpack_z_test21 ( );
  linpack_z_test22 ( );
  linpack_z_test23 ( );
  linpack_z_test24 ( );
  linpack_z_test25 ( );
  linpack_z_test26 ( );
  linpack_z_test27 ( );
  linpack_z_test28 ( );
  linpack_z_test29 ( );

  linpack_z_test30 ( );
  linpack_z_test31 ( );
  linpack_z_test32 ( );
  linpack_z_test33 ( );
  linpack_z_test34 ( );
  linpack_z_test345 ( );
  linpack_z_test35 ( );
  linpack_z_test36 ( );
  linpack_z_test37 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINPACK_Z_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ x, w ] = ccoh_rule ( n )

%*****************************************************************************80
%
%% CCOH_RULE computes a Clenshaw Curtis Open Half rule.
%
%  Discussion:
%
%    Our convention is that the abscissas are numbered from left to right.
%
%    The rule is defined on [-1,1].
%
%    The integral to approximate:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  x =  ( 1 : 2 : 2 * n - 1 )' / ( 2 * n );
%
%  Chebyshev transformation from [0,1] to [-1,1].
%
  x = cos ( x * pi );
%
%  Compute the weights.
%
  x_min = -1.0;
  x_max = +1.0;

  w = nc_compute ( n, x_min, x_max, x );

  return
end


text:
classdef MLSSolver < pf2odfSolver
%
% The class MLSSolver implements the modified least squares solver for
% reconstructing an ODF from arbitrarily scattered pole figure intensities.
% The resulting ODF is represented as a weighted sum of unimodal
% components. The shape and the number of component centers can be
% specified. The algorithm is explained in detail in *A novel pole figure
% inversion method: specification of the MTEX algorithm*, Hielscher,
% Schaeben: J. of Appl. Cryst., 41(6), 2008.
%
% Syntax
%
%   solver = MLSSolver(pf,'resolution',5*degree,'halfwidth',7.5*degree);
%   [odf,alpha] = solver.calcODF;
%
% Class Properties
%  psi     - @SO3Kernel describing the shape of the unimodal components
%  S3G     - the centers of the unimodal components as @SO3Grid in orientation space
%  c       - weighting coefficients to the unimodal components
%  weights -
%  zrm     - @zeroRangeMethod
%  ghostCorrection - whether to use ghost correction
%  iterMax - max number of iterations
%  iterMin - min number of iterations
%
% Dependent Class Properties
%  odf - @SO3Fun the reconstructed ODF
%
% See also
% PoleFigureTutorial PoleFigure2ODFAmbiguity PoleFigure2ODFGhostCorrection

  properties
    psi     % SO3Kernel function
    S3G     % SO3Grid
    c       % current coefficients
    c0      % constant portion
    weights % cell
    zrm     % zero range method
    ghostCorrection = 1
    iterMax = 10; % max number of iterations
    iterMin = 5;  % max number of iterations
    lambda = 0;   % regularisation parameter
    RM = [];      % regularization matrix
  end

  properties (Access = private)
    nfft_gh  % list of nfft plans
    nfft_r   % list of nfft plans
    A        % legendre coefficients of the kernel function
    refl     % cell
    u
    a
    alpha
  end

  properties (Dependent = true)
    odf % the reconstructed @SO3Fun
  end

  methods
    function solver = MLSSolver(pf,varargin)
      % constructor

      if nargin == 0, return; end

      % ensure intensities are non negative
      %solver.pf = unique(max(pf,0));
      solver.pf = max(pf,0);

      % normalize very different polefigures
      mm = max(pf.intensities(:));

      for i = 1:pf.numPF
        if mm > 5*max(pf.allI{i}(:))
          pf.allI{i} = pf.allI{i} * mm/5/max(pf.allI{i}(:));
        end
      end

      % generate discretization of orientation space
      solver.S3G = getClass(varargin,'SO3Grid');
      if isempty(solver.S3G)
        if pf.allR{1}.isOption('resolution')
          res = pf.allR{1}.resolution;
        else
          res = 5*degree;
        end
        res = get_option(varargin,'resolution',res);
        solver.S3G = equispacedSO3Grid(solver.pf.CS,solver.pf.SS,'resolution',res);
      end

      % zero range method
      if check_option(varargin,{'ZR','zero_range','zeroRange'})
        solver.zrm = zeroRangeMethod(solver.pf);
      end

      % get kernel
      psi = SO3DeLaValleePoussinKernel('halfwidth',...
        get_option(varargin,{'HALFWIDTH','KERNELWIDTH'},solver.S3G.resolution,'double'));
      solver.psi = getClass(varargin,'SO3Kernel',psi);

      % get other options
      solver.iterMin = get_option(varargin,'iterMin',solver.iterMin);
      solver.iterMax = get_option(varargin,'iterMax',solver.iterMax);

      % start vector
      solver.c = get_option(varargin,'C0',[]);

      % ghost correction
      solver.ghostCorrection = ~check_option(varargin,'noGhostCorrection');

      % compute quadrature weights
      if numProper(solver.SS) == 1
        solver.weights = cellfun(@(r) calcQuadratureWeights(r),solver.pf.allR,'UniformOutput',false);
      else
        solver.weights = num2cell(1./length(pf,[]));
      end

      % regularisation
      solver.lambda = get_option(varargin,'regularisation',0);

    end

    function delete(solver)
      % destructor

      % free all nfft plans
      for i = 1:length(solver.nfft_gh)
        nfsftmex('finalize',solver.nfft_gh(i));
      end

      for i = 1:length(solver.nfft_r)
        nfsftmex('finalize',solver.nfft_r(i));
      end

    end

    function odf = get.odf(solver)
      if solver.c0 > 0.99
        odf = SO3FunHarmonic(1,solver.CS, solver.SS);
      else
        odf = SO3FunRBF(solver.S3G,solver.psi,...
          (1-solver.c0)*solver.c./sum(solver.c), solver.c0);
      end
    end

  end

  methods (Static = true)

    function check

      cs = crystalSymmetry('222');
      odf = unimodalODF(orientation.id(cs));
      r = equispacedS2Grid('upper','antipodal','resolution',5*degree);
      h = Miller({1,0,0},{1,1,1},{1,1,0},cs);
      pf = calcPoleFigure(odf,h,r);
      odf = calcODF(pf)

      plotPDF(odf,h)
      plotFibre(odf,fibre.alpha(cs))

      mtexdata dubna
      tic
      solver = MLSSolver;
      solver.pf = pf;
      solver.S3G = equispacedSO3Grid(pf.CS,'resolution',2.5*degree);
      solver.psi = SO3DeLaValleePoussinKernel('halfwidth',2.5*degree);
      solver.weights = repcell(1,numPF(pf),1);
      solver.c = ones(length(solver.S3G),1) ./ length(solver.S3G);

      solver.init
      odf = solver.calcODF;
      delete(solver);
      toc

    end

  end

end


text:
function p = dirichlet_logProb_fast(a, bar_p)

p = gammaln(sum(a)) - sum(gammaln(a)) + sum((a-1).*bar_p);
K = length(a);
flops(flops + (K+1)*flops_digamma + 3*K);

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function z = power( x, y, op )

%POWER   Internal cvx version.

persistent BP
if isempty( BP ),
    BP.map = cvx_remap( ...
        ... % Invalid combinations: zero ^ negative, posynomial ^ negative
        { { 'zero', 'g_posynomial', 'gg_monomial' }, { 'negative' } }, ...
        ... % Constant
        { { 'constant' } }, ...
        ... % constant ^ convex/concave
        { { 'positive' }, { 'convex', 'concave' } }, ...
        ... % geometric ^ real
        { { 'l_valid_' }, { 'real' } }, ...
        ... % other non-constant ^ real
        { { 'valid' }, { 'real' } }, [0,1,2,3,4] );
    BP.funcs = { @power_c, @power_e, @power_g, @power_p };
    BP.constant = 1;
end
if isnumeric( y ) && numel( y ) == 1,
    switch y,
        case   1, z = x;           return
        case   2, z = square( x ); return
        case 0.5, z = sqrt( x );   return
        case  -1, z = recip( x );  return
    end
end
if nargin < 3, op = '.^'; end
BP.name = op;
z = cvx_binary_op( BP, x, y );

function z = power_c( x, y )
z = builtin( 'power', x, y );

function z = power_e( x, y )
z = exp( log( cvx_constant( x ) ) .* y );

function z = power_g( x, y )
z = exp( log( x ) .* y );

function z = power_p( x, y )
nx = numel( x );
ny = numel( y );
y  = cvx_constant( y );
vx = cvx_classify( x );
nz = max( nx, ny );
if ny == 1
    pu = y;
    multp = false;
else
    pu = sort( y' ); %#ok'
    pu = pu( [ true, diff(pu) ~= 0 ] );
    multp = length( pu ) ~= 1;
end
if multp,
    zp = cvx( nz, [] );
    vp = vx;
    xp = x;
    yp = y;
end
errs = {};
for pk = pu,
    if multp,
        tp = yp == pk;
        if nx > 1,
            x = cvx_fastref( xp, tp );
            x = cvx( [size(x,2),1], x );
            vx = vp( tp );
        end
        if ny > 1,
            y = yp( tp );
        end
    end
    y1 = y(1);
    if y == 2,
        z = square( x );
    elseif y1 > 1,
        z = power_cvx( x, vx, y1 );
    elseif y1 == 1,
        z = x;
    elseif y1 > 0,
        z = power_ccv( x, vx, y1 );
    elseif y1 == 0,
        z = cvx( ones(size(x)) );
    else
        z = power_neg( x, vx, y1 );
    end
    if islogical( z ),
        errs(end+1,:) = { cvx_subsref( x, z ), y }; %#ok
        continue
    end
    if nx ~= nz,
        z = repmat( z, [nz,1] );
    end
    if multp,
        zp = cvx_fastasgn( zp, tp, z );
    end
end
if multp
    z = zp;
end
if ~isempty( errs ),
    cvx_dcp_error( errs, op );
elseif multp
    z = zm;
end

%
% P > 1, integer: X affine, p-convex, n-concave
% P > 1, non-integer: X affine, p-convex
%

function y = power_cvx_fast( x, p ) %#ok
ne = 0;
nx = numel( x );
while rem(p,2) == 0,
    ne = ne + 1;
    p = p * 0.5;
end
cvx_begin
    epigraph variable y(nx) nonnegative_
    if ne > 0,
        variables z2(nx,ne-1)
        if p == 1,
            w = y;
        else
            variable w(nx)
        end
        { [x,z2], 0.5, [z2,w] } == rotated_lorentz( [ nx, ne ], 3 ); %#ok
    else
        w = x;
    end
    if p > 1,
        { [y,ones(nx,1)], w } == geo_mean_cone( [ nx, 2 ], 2, [1/p,1-1/p], 'func' );  %#ok
    end
cvx_end

function y = power_cvx( x, v, p )
persistent remap remap_i remap_e
if isempty( remap_e ),
    remap   = cvx_remap( 'affine', 'p_convex' );
    remap_i = cvx_remap( 'p_convex' ) - cvx_remap( 'n_concave' );
    remap_e = remap_i | remap;
end
if ~isempty( v ),
    isint = rem( p, 1 ) == 0;
    if isint,
        isevn = rem( p, 2 ) == 0;
        if isevn,
            v = remap_e( v );
        else
            v = remap_i( v );
        end
    else
        v = remap( v );
    end
    if ~all( v ),
        y = v == 0;
        return
    end
    v = 1 - 2 * ( v(:) == 3 );
    x = linearize( v .* x );
else
    isint = false;
end
y = power_cvx_fast( x, p );
if isint,
    y = v .* y;
end

%
% 0 < P < 1: X concave
%

function y = power_ccv( x, v, p )
persistent remap
if isempty( remap ),
    remap = ~cvx_remap( 'concave' );
end
v = remap( v );
if any( v ),
    y = v;
    return
end
cvx_begin
    hypograph variable y(numel(x)) nonnegative_
    power_cvx_fast( y, 1.0 / p ) <= x; %#ok
cvx_end

%
% P < 0, integer: p-concave, n-convex
% P < 0, non-integer: p-concave
%

function y = power_neg( x, v, p )
persistent remap remap_i
if isempty( remap_i ),
    remap   = cvx_remap( 'p_concave' );
    remap_i = remap - cvx_remap( 'n_convex' );
end
if rem( p, 1 ) == 0,
    % Integer: positive and negative lobes
    v = remap_i( v );
else
    % Non-integer: positive lobe only
    v = remap( v );
end
if ~all( v ),
    y = v == 0;
    return
end
nx = numel(x);
z  = linearize( v .* x );
cvx_begin
    epigraph variable y(nx)
    if p == 1,
        { 2, z, y } == rotated_lorentz( [nx,1], 2, 0 ); %#ok
    else
        { [z,y], 1 } == geo_mean_cone( [nx,2], 2, [-p,1], 'func' ); %#ok
    end
    cvx_setnneg(y);
cvx_end
y = v .* y;

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function id_searchwindows( ref_VAD, ref_Nsamples, deg_VAD, deg_Nsamples);

global MINUTTLENGTH Downsample MINUTTLENGTH SEARCHBUFFER
global Crude_DelayEst Nutterances UttSearch_Start UttSearch_End

Utt_num = 1;
speech_flag = 0;

VAD_length= floor( ref_Nsamples/ Downsample);
del_deg_start= MINUTTLENGTH- Crude_DelayEst/ Downsample;
del_deg_end= floor((deg_Nsamples- Crude_DelayEst)/ Downsample)-...
    MINUTTLENGTH;

for count= 1: VAD_length
    VAD_value= ref_VAD(count);
    if( (VAD_value> 0) && (speech_flag== 0) )
        speech_flag= 1;
        this_start= count;
        UttSearch_Start(Utt_num)= count- SEARCHBUFFER;
        if( UttSearch_Start(Utt_num)< 0 )
            UttSearch_Start(Utt_num)= 0;
        end
    end

    if( ((VAD_value== 0) || (count == (VAD_length-1))) && ...
            (speech_flag == 1) )
        speech_flag = 0;
        UttSearch_End(Utt_num) = count + SEARCHBUFFER;
        if( UttSearch_End(Utt_num) > VAD_length - 1 )
            UttSearch_End(Utt_num) = VAD_length -1;
        end

        if( ((count - this_start) >= MINUTTLENGTH) &&...
                (this_start < del_deg_end) &&...
                (count > del_deg_start) )
            Utt_num= Utt_num + 1;
        end
    end
end
Utt_num= Utt_num- 1;
Nutterances = Utt_num;

% fprintf( 1, 'Nutterances is %d\n', Nutterances);

% fid= fopen( 'mat_utt.txt', 'wt');
% fprintf( fid, '%d\n', UttSearch_Start( 1: Nutterances));
% fprintf( fid, '\n');
% fprintf( fid, '%d\n', UttSearch_End( 1: Nutterances));
% fclose(fid);

author: vipchengrui

repo: traditional-speech-enhancement


text:
function varargout = transformPoint(varargin)
%TRANSFORMPOINT Apply an affine transform to a point or a point set.
%
%   PT2 = transformPoint(PT1, TRANSFO);
%   Returns the result of the transformation TRANSFO applied to the point
%   PT1. PT1 has the form [xp yp], and TRANSFO is either a 2-by-2, a
%   2-by-3, or a 3-by-3 matrix,
%
%   Format of TRANSFO can be one of :
%   [a b]   ,   [a b c] , or [a b c]
%   [d e]       [d e f]      [d e f]
%                            [0 0 1]
%
%   PT2 = transformPoint(PT1, TRANSFO);
%   Also works when PTA is a N-by-2 array representing point coordinates.
%   In this case, the result PT2 has the same size as PT1.
%
%   [X2, Y2] = transformPoint(X1, Y1, TRANS);
%   Also works when PX1 and PY1 are two arrays the same size. The function
%   transforms each pair (PX1, PY1), and returns the result in (X2, Y2),
%   which has the same size as (PX1 PY1).
%
%
%   See also
%     points2d, transforms2d, translation, rotation
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-04-06
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% parse input arguments
if length(varargin) == 2
    var = varargin{1};
    px = var(:,1);
    py = var(:,2);
    trans = varargin{2};
elseif length(varargin) == 3
    px = varargin{1};
    py = varargin{2};
    trans = varargin{3};
else
    error('wrong number of arguments in "transformPoint"');
end

% apply linear part of the transform
px2 = px * trans(1,1) + py * trans(1,2);
py2 = px * trans(2,1) + py * trans(2,2);

% add translation vector, if exist
if size(trans, 2) > 2
    px2 = px2 + trans(1,3);
    py2 = py2 + trans(2,3);
end

% format output arguments
if nargout < 2
    varargout{1} = [px2 py2];
elseif nargout
    varargout{1} = px2;
    varargout{2} = py2;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% GET THE 3D TRAJECTORY TRAILS AS A GIF  [ UPDATED ]
function [currentFigure,figureHandle] = GetFigure_isometricGif(SIM,objectIndex,DATA,currentFigure)
% This function generates an animated .gif representing the object
% trajectories over the complete timeseries.

% CONFIGURE THE PLOT ATTRIBUTES
figurePath = strcat(SIM.outputPath,'isometric.gif');
figureHandle = figure('Name','OpenMAS isometric timelapse (GIF)');

% MAXIMISE GRAPH SIZE IN WINDOW
setappdata(figureHandle, 'SubplotDefaultAxesLocation', [0.08, 0.08, 0.90, 0.88]);
set(figureHandle,'Position',DATA.figureProperties.windowSettings);        % [x y width height]
set(figureHandle,'Color',DATA.figureProperties.figureColor);               % Background colour
if DATA.figureProperties.publish
    set(figureHandle,'MenuBar','none');
    set(figureHandle,'ToolBar','none');
    set(figureHandle,'Visible','off');
end
% AXES SETTINGS
ax = axes(figureHandle);

% GET THE SIZE OF THE PLOTTED SET
% The data has already been pre-formatted so that the number of states
% aligns with the number of expected frames in the annimation.
% IMPORT THE OBJECT DATA FROM TEMP FILE
load([SIM.outputPath,SIM.systemFile]);
framesToPlot = size(eval(sprintf('objectID%d',SIM.OBJECTS(1).objectID)),2);   % The variable name in the workspace
F(framesToPlot) = struct('cdata',[],'colormap',[]);

% ax.NextPlot = 'replaceChildren';
inclinationAngle = 62;
viewPoint = -10;
viewRate = 0.05;
for frame = 1:framesToPlot
    hold on;
    % Must be done on a step by step base to get the annimations correct.

    % BUILD THE FRAME
    if frame == 1
        % UPDATE OBJECT HANDLES
        for entity = 1:SIM.totalObjects
            % GET THE OBJECT HANDLE
            objectHandle = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entity).objectID};
            % PULL VARIABLE FROM WORKSPACE
            variableLabel = sprintf('objectID%d',SIM.OBJECTS(entity).objectID); % The variable name in the workspace
            globalStates = eval(variableLabel);                              % Evalutate object states

            % TRACE INITIAL POSITIONS
            traceHandle(entity) = plot3(ax,globalStates(1,frame),globalStates(2,frame),globalStates(3,frame));
            if numel(objectHandle.GEOMETRY.vertices) > 0
                % GET THE GLOBAL POSE AT STEP
                [R_frame] = OMAS_geometry.quaternionToRotationMatrix(globalStates(7:end,frame));
                % BUILD MARKER
                markerHandle(entity) = patch(ax,...
                    'Vertices',objectHandle.GEOMETRY.vertices*R_frame + globalStates(1:3,frame)',...
                    'Faces',objectHandle.GEOMETRY.faces,...
                    'FaceColor',SIM.OBJECTS(entity).colour,...
                    'EdgeColor',DATA.figureProperties.edgeColor,...
                    'EdgeAlpha',DATA.figureProperties.edgeAlpha,...
                    'FaceLighting',DATA.figureProperties.faceLighting,...
                    'FaceAlpha',DATA.figureProperties.faceAlpha,...
                    'LineWidth',DATA.figureProperties.patchLineWidth);             % Patch properties
            else
                % INITIAL MARKER PLOT
                markerHandle(entity) = plot3(ax,globalStates(1,frame),globalStates(2,frame),globalStates(3,frame));
                % ASSIGN MARKER PROPERTIES
                markerHandle(entity).Marker = SIM.OBJECTS(entity).symbol;
                markerHandle(entity).MarkerSize = DATA.figureProperties.markerSize;
                markerHandle(entity).MarkerFaceColor = SIM.OBJECTS(entity).colour;
                markerHandle(entity).MarkerEdgeColor = DATA.figureProperties.markerEdgeColor;
                markerHandle(entity).Color = SIM.OBJECTS(entity).colour;
            end
            % ASSIGN TRACE PROPERTIES
            traceHandle(entity).LineStyle = DATA.figureProperties.lineStyle;
            traceHandle(entity).LineWidth = DATA.figureProperties.lineWidth;
            traceHandle(entity).Color = SIM.OBJECTS(entity).colour;
        end
        % ANNOTATION WITH THE CURRENT TIME
        clockHandle = annotation('textbox',[0.025 0.025 0.15 0.06],'String','Time:','FitBoxToText','off');
        set(clockHandle,...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'FontSize',DATA.figureProperties.axisFontSize,...
            'FontWeight',DATA.figureProperties.fontWeight);
        % Title
        title(ax,...
            sprintf('Object global trajectories over a period of %ss',num2str(SIM.TIME.endTime)),...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'FontSize',DATA.figureProperties.titleFontSize,...
            'FontSmoothing','on');
        % X-Label
        xlabel(ax,'x (m)',...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'FontSize',DATA.figureProperties.axisFontSize,...
            'FontSmoothing','on');
        % Y-Label
        ylabel(ax,'y (m)',...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'FontSize',DATA.figureProperties.axisFontSize,...
            'FontSmoothing','on');
        % Z-Label
        zlabel(ax,'z (m)',...
            'Interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'FontSize',DATA.figureProperties.axisFontSize,...
            'FontSmoothing','on');
        % Axes
        set(ax,...
            'TickLabelInterpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName,...
            'Fontweight',DATA.figureProperties.fontWeight,...
            'FontSize',DATA.figureProperties.axisFontSize,...
            'FontSmoothing','on',...
            'Color',DATA.figureProperties.axesColor,...
            'GridLineStyle','--',...
            'GridAlpha',0.25,...
            'GridColor','k');

        axisLimit = (DATA.figureProperties.objectMaximalRadii + DATA.figureProperties.maxAbsPosition);
        xlim(ax,[-axisLimit,axisLimit]);
        ylim(ax,[-axisLimit,axisLimit]);
        zlim(ax,[-axisLimit,axisLimit]);
        axis manual;
        % Legend
        legend(markerHandle,DATA.figureProperties.legendEntries,...
            'Location','northeastoutside',...
            'interpreter',DATA.figureProperties.interpreter,...
            'fontname',DATA.figureProperties.fontName);
        view([viewPoint inclinationAngle]);     % Set initial view angle
        grid on; box on; hold off; grid minor;
    else
        % CONTINUED FRAMES
        % FOR EACH PLOT ENTITY
        for entity = 1:SIM.totalObjects
            % GET THE OBJECT HANDLE
            objectHandle = objectIndex{SIM.globalIDvector == SIM.OBJECTS(entity).objectID};
            % PULL VARIABLE FROM WORKSPACE
            variableLabel = sprintf('objectID%d',SIM.OBJECTS(entity).objectID); % The variable name in the workspace
            globalStates = eval(variableLabel);                                 % Evalutate object states
            % CHECK IF UPDATE IS NECESSARY
            if any(isnan(globalStates(:,frame)))                                % Object is static, freeze its position
                continue
            end

            % HANDLE DIFFERENT REPRESENTATIONS
            if numel(objectHandle.GEOMETRY.vertices) > 0
                % GET THE GLOBAL POSE AT STEP
                [R_frame] = OMAS_geometry.quaternionToRotationMatrix(globalStates(7:end,frame));
                % BUILD MARKER
                set(markerHandle(entity),'Vertices',objectHandle.GEOMETRY.vertices*R_frame + globalStates(1:3,frame)');
            else
                % UPDATE MARKER DATA
                set(markerHandle(entity),'XData',globalStates(1,frame));
                set(markerHandle(entity),'YData',globalStates(2,frame));
                set(markerHandle(entity),'ZData',globalStates(3,frame));
            end

            % UPDATE TRACE DATA
            if frame <= DATA.figureProperties.tailLength/SIM.TIME.dt
                % TRANSITIONING PERIOD
                tailTrace = globalStates(1:3,1:frame);                     % All points upto the step position
            elseif (frame - DATA.figureProperties.tailLength/SIM.TIME.dt) > 0
                % IF THE TRAIL IS NOW A SUBSET
                tailTrace = globalStates(1:3,(frame-(DATA.figureProperties.tailLength/SIM.TIME.dt)):frame);        % All points upto the step position
            else
                tailTrace = globalStates(1:3,frame);
            end
            % UPDATE TRACE DATA
            set(traceHandle(entity),'XData', tailTrace(1,:));
            set(traceHandle(entity),'YData', tailTrace(2,:));
            set(traceHandle(entity),'ZData', tailTrace(3,:));
        end
    end
    % UPDATE TIMESTAMP ANNOTATION
    set(clockHandle,'String',sprintf('Time: %ss',num2str(frameTimeVector(frame))));
    % FORCE IMAGE WRITE
    drawnow();

    % COLLECT FRAMES
    F(frame) = getframe(figureHandle);
    im = frame2im(F(frame));
    [imind,cm] = rgb2ind(im,256);

    % APPEND THE FRAMES TO GIF
    if frame == 1
        imwrite(imind,cm,figurePath,'gif','Loopcount',inf,'DelayTime',(1/DATA.figureProperties.fps));
    else
        imwrite(imind,cm,figurePath,'gif','WriteMode','append','DelayTime',(1/DATA.figureProperties.fps));
    end
%     viewPoint = viewPoint + viewRate;
end
% FIGURE COMPLETE
currentFigure = currentFigure + 1;
% CLEAN UP
clearvars -except currentFigure figureHandle
end

repo: OpenMAS


text:
%  Figure 10.54      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_54.m is a script to generate Fig. 10.54, the step response
% for the fuel/air ratio with a nonlinear sensor.  A (3,3)
% Pade approximant is used to approximate the delay.
clf;
% construct the F/A dynamics with the sensor time-constant too
f =[-50     0     0;
     0    -1     0;
    10    10   -10];

g =[25.0000;
    0.5000;
         0];
h =[0     0     1];
j=0;
[f3,g3,h3,j3 ]=pade(.2,3); %  the delay Pade model
% put the plant and the delay in series
[fp,gp,hp,jp]=series(f,g,h,j,f3,g3,h3,j3);
% construct the controller
nc=.1*[1 .3];
dc=[1 0]; % the PI controller in polynomial form
[fc,gc,hc,jc]= tf2ss(nc, dc); % put the controller in state space form
[fol,gol,hol,jol] = series(fc,gc,hc,jc,fp,gp,hp,jp);
% get a discrete model for Ts = .01
[phi,gam] = c2d(fol,gol,.01);
%
% form the closed-loop difference equation
x=[0 0 0 0 0 0 0]';
yout=[];

for t=0:.01:20;
 y=hol*x ;
 e=.4*sat(150*(.068-y));
 x=phi*x+gam*e;
 yout=[yout [e;y]];
end

t=0:.01:20;
subplot(211);
plot(t,yout(1,:));
axis( [0 20 -.6 .6]);
xlabel('Time (sec)');
ylabel('e');
title('Fig. 10.54 (b) Error plot for nonlinear control of F/A')
grid on;

subplot(212);
plot(t,yout(2,:));
axis([0 20 0 .09]);
xlabel('Time (sec)');
ylabel('F/A');
title('Fig. 10. 54 (c) Output F/A ratio for nonlinear control')
hold off;
%grid
nicegrid


text:
function [val] = filetype_check_header(filename, head, offset)

% FILETYPE_CHECK_HEADER helper function to determine the file type
% by reading the first number of bytes of a file and comparing them
% to a known string (c.f. magic number).

% Copyright (C) 2003-2022 Robert Oostenveld
%
% 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$

% these are for speeding up subsequent calls with the same input arguments
persistent previous_argin previous_argout

if nargin<3
  offset = 0;
end

current_argin = {filename, head, offset};
if ~isempty(previous_argin) && isequal(current_argin, previous_argin)
  % don't do the detection again, but return the previous value from cache
  val = previous_argout;
  return
end

if iscell(filename)
  % compare the header of multiple files
  val = zeros(size(filename));
  for i=1:length(filename)
    val(i) = filetype_check_header(filename{i}, head, offset);
  end
elseif isfolder(filename)
  % a directory cannot have a header
  val = false;
elseif ~exist(filename, 'file')
  val = false;
  cache = false; % the file does not exist now, but can exist later
else
  % read the first few bytes from the file and compare them to the desired header
  fid = fopen(filename, 'rb');
  if fid<0
    ft_warning('could not open %s', filename);
    val = false;
  else
    fseek(fid, offset, 'cof');
    if iscell(head)
      len = zeros(size(head));
      for i=1:length(head)
        len(i) = length(head{i});
      end
      [str, siz] = fread(fid, max(len), 'uint8');
      fclose(fid);
      for i=1:length(head)
        val = all(str(1:len(i))==head{i}(:));
        if val
          break
        end
      end
    else
      [str, siz] = fread(fid, length(head), 'uint8=>char');
      fclose(fid);
      if siz~=length(head)
        ft_warning('could not read the header from %s', filename);
        val = false;
      else
        val = all(str(:)==head(:));
      end
    end
  end
end

% remember the current input and output arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
current_argout  = val;
previous_argin  = current_argin;
previous_argout = current_argout;

return % main()


text:
function f = replaceBoundaryRoots(f)
%REPLACEBOUNDARYROOTS  Replace boundary roots of a SINGFUN.
%  F = REPLACEBOUNDARYROOTS(F) returns a SINGFUN which has both exponents less
%  than 1 by absorbing the integer part of any boundary exponents larger than 1
%  into its smoothPart.
%
% See also EXTRACTBOUNDARYROOTS.

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

% Grab the exponents:
exps = get(f, 'exponents');

% Grab the indice for exponents larger or equal to 1:
ind = ( exps >= 1 );

% Both exponents are less than 1:
if ( ~any( ind ) )
    return
end

% Sort out the new exponents and the order of the boundary roots which need to
% be absorbed into the smoothPart:
newExps = exps;
newExps(ind) = exps(ind) - floor(exps(ind));
pow = exps - newExps;

% Compute the factor from function roots:
mult = singfun.constructSmoothPart(@(x) (x+1).^pow(1).*(1-x).^pow(2), [], []);
f.smoothPart = f.smoothPart.*mult;

% Update the exponents:
f.exponents = newExps;

end


text:
function UNew = elasticPeriodic2D(varargin);
% elasticPeriodic2D: solve elastic registraion in 2D with periodic
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[U,F,mu,lambda,gamma,PixSize,NumPix] = parse_inputs(varargin{:});

% construct filters that implement discretized Navier-Lame equations
d1 = [1;-2;1]/(PixSize(1)^2);
d2 = [1 -2 1]/(PixSize(2)^2);
d12 = [1 0 -1;0 0 0;-1 0 1]/(4*PixSize(1)*PixSize(2));

[A11,A22] = deal(zeros(3,3));
A11(2,2) = gamma;
A11(:,2) = A11(:,2) + (lambda+2*mu)*d1;
A11(2,:) = A11(2,:) + mu*d2;
A22(2,2) = gamma;
A22(:,2) = A22(:,2) + mu*d1;
A22(2,:) = A22(2,:) + (lambda+2*mu)*d2;

A12 = d12*(lambda+mu)/4;
A21 = A12;

% add displacement vectors to multiple of force field
F = gamma*U + F;

% multiply force field by adjoint of Navier-Lame equations
Fnew = zeros(NumPix(1),NumPix(2),2);
Fnew(:,:,1) = imfilter(F(:,:,1),A22,'replicate') - imfilter(F(:,:,2),A12,'replicate');
Fnew(:,:,2) = imfilter(F(:,:,2),A11,'replicate') - imfilter(F(:,:,1),A21,'replicate');

% compute Fourier transform of new force field
FnewF1 = fft2(Fnew(:,:,1));
FnewF2 = fft2(Fnew(:,:,2));

% construct images of coordinates scaled by pi/(N or M)
[alpha,beta] = ndgrid(2*pi*(0:(NumPix(1)-1))/NumPix(1),2*pi*(0:(NumPix(2)-1))/NumPix(2));

% construct LHS factor
T = 2*cos(alpha) + 2*cos(beta) - 4;
LHSfactor = (gamma + (lambda+2*mu).*T).*(gamma + mu.*T);

% if gamma is zero, set origin term to 1, as DC term does not matter
if isequal(gamma,0),
    LHSfactor(1,1) = 1;
end

% solve for FFT of U
UF1 = FnewF1./LHSfactor;
UF2 = FnewF2./LHSfactor;

% perform inverse fft and concatenate
UNew = cat(3,ifft2(UF1,'symmetric'),ifft2(UF2,'symmetric'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,F,mu,lambda,gamma,PixSize,NumPix] = parse_inputs(varargin);

% get displacement field and check size
U = varargin{1};
F = varargin{2};
gamma = varargin{3};
PixSize = varargin{4}(1:2);
NumPix = [varargin{5} varargin{6}];
mu = varargin{8};
lambda = varargin{9};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| 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                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtDomIntSing2D.m                             |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 25.11.2018                                    |
%|  / 0 \  |   LAST MODIF : 25.11.2018                                    |
%| ( === ) |   SYNOPSIS   : Singular integration over a triangle          |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

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

% Edge definition
a = 0;
b = 1;
N = 100;

% Reference mesh (edge)
vtx  = [a 0 0 ; b 0 0];
elt  = [1 2];
mesh = msh(vtx,elt);

% Observation points (100 random & 11 segments points)
X = [-2 + 5*rand(N,3);
    [(-2:0.2:3)' zeros(26,2)] ];
X(:,3) = 0;

% Graphical representation
plot(mesh,1)
hold on
plotNrm(mesh,'b')
plot3(X(:,1),X(:,2),X(:,3),'or')
axis equal
view(0,90)
grid on

% Analytical integration (exact)
tic
S     = mesh.vtx;
n     = mesh.nrm;
tau   = mesh.tgt;
[logRa,rlogRa,gradlogRa] = domSemiAnalyticInt2D(X,S,n,tau);
toc

% 1D simpson integration ("exact")
tic
logRs     = zeros(size(logRa));
rlogRs    = zeros(size(rlogRa));
gradlogRs = zeros(size(gradlogRa));
for i = 1:size(X,1)
    % Function |r| in the 2D plane
    fun = @(x1) sqrt( (x1-X(i,1)).^2 + X(i,2).^2 );

    % Scalar integration log(|r|)
    logRs(i) = integral(@(x1) log(fun(x1)),a,b,'Reltol',1e-12);

    % Vectorial integration rlog(|r|)
    rlogRs(i,1) = integral(@(x1) (x1-X(i,1)).*log(fun(x1)),a,b,'Reltol',1e-12);
    rlogRs(i,2) = -X(i,2) .*  logRs(i);
    rlogRs(i,3) = 0;

    % Vectoral integration gradlog(|r|)
    if (X(i,2) == 0) && (X(i,1)>=a && X(i,1)<=b)
        gradlogRs(i,:) = 0;
    else
        gradlogRs(i,1) = integral(@(x1) (x1-X(i,1))./fun(x1).^2,a,b);
        gradlogRs(i,2) = integral(@(x1) -X(i,2)./fun(x1).^2,a,b);
        gradlogRs(i,3) = 0;
    end
end
toc

% Relative errors log(|r|)
figure
semilogy( abs(logRa-logRs)./abs(logRs) )
title('Relative error (log) : \int log(|r|)')
grid on
I = 1:N;
disp('Relative error (Linf) : \int log(|r|)')
norm(logRa(I)-logRs(I),'inf')./norm(logRa(I),'inf')
norm(logRa-logRs,'inf')./norm(logRa,'inf')

% Relative errors rlog(|r|)
figure
semilogy( abs(rlogRa-rlogRs)./abs(rlogRs) )
title('Relative error (log) : \int rlog(|r|)')
grid on
I = 1:N;
disp('Relative error (Linf) : \int rlog(|r|)')
norm(rlogRa(I)-rlogRs(I),'inf')./norm(rlogRa(I),'inf')
norm(rlogRa-rlogRs,'inf')./norm(rlogRa,'inf')

% Relative errors gradlog(|r|)
figure
semilogy( abs(gradlogRa-gradlogRs)./abs(gradlogRs) )
title('Relative error (log) : \int gradlog(|r|)')
grid on
I = 1:N;
disp('Relative error (Linf) : \int gradlog(|r|)')
norm(gradlogRa(I)-gradlogRs(I),'inf')./norm(gradlogRa(I),'inf')
norm(gradlogRa-gradlogRs,'inf')./norm(gradlogRa,'inf')

disp('~~> Michto gypsilab !')


text:
function im_out = linear_stretch(im_in)
% The top 0.00001 pixels of the image are clipped to 1, and the
% image is linearly stretched.
img1 = max(im_in, [],3);   % maximum between RGB channels
img1 = img1(:);
img1 = sort(img1, 'descend');
idx = round(0.00001*length(img1));
scale = img1(idx);
im_out = im_in./scale;
im_out = min(im_out,1);

author: danaberman

repo: underwater-hl


text:
classdef FreeFemFileIdenticalPrinter < FreeFemFilePrinter

    methods (Access = protected)

       function printingLines(obj)
            nlines = numel(obj.linesToPrint);
            for iline = 1:nlines
                line = obj.linesToPrint{iline};
                fprintf(obj.fileID,'%s\r\n',line);
            end
        end

    end

end


text:
report_this_filefun(mfilename('fullpath'));

l = isnan(tmap);
tmap(l) = 1;

%l = tmap <  0;
%tmap(l) = 0;

%l = Z < 0.12;
%Z(l) = 0.14;

[lat,lon] = meshgrat(tmap,tmapleg);
%[lat,lon] = meshgrat(vlat,vlon);

%[smap,smapleg] = country2mtx('peru',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

%sw = interp2(lon0,lat0,smap,lon,lat);
%l =  isnan(sw) == 1 | sw == 2 ;

ren = interp2(X,Y,Z,lon,lat);

mi = min(min(ren));

ren(l) = mi-0.01;

% ren = interp2(X,Y,Z,lon,lat);

% mi = min(min(ren));
l =  isnan(ren);
ren(l) = 0;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[s4_south s3_north],'MapLonLimit',[s2_west s1_east])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',10);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

load worldlo
h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',2)
% h2 = displaym(PPpoint);
%h = displaym(PPtext); trimcart(h);

%pl = plotm(lima(:,2),lima(:,1),'hw');
%set(pl,'LineWidth',1.5,'MarkerSize',12,...
% 'MarkerFaceColor','y','MarkerEdgeColor','k')

pl = plotm(dam(:,2),dam(:,1),'^w');
set(pl,'LineWidth',1.5,'MarkerSize',12,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')
%pl = plotm(coastline(:,2), coastline(:,1),'w','Linewidth',2);

ci = worldlo('PPpoint');
cx = ci(1).long;
cy = ci(1).lat;
hold on
plotm(cy,cx,'sk','Markersize',12,'Markerfacecolor',[ 1 1 1])

ri = worldlo('DNline');
hold on

rx = [ri(1).long ; ri(2).long];
ry = [ri(1).lat ; ri(2).lat ];
hold on
% plotm(ry,rx,'b','Linewidth',2);

zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface

j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.93 0.9 0.9] ; j];

caxis([0. 0.28])
colormap(j);

axis off; set(gcf,'color','w')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','k','flinewidth',3);

setm(gca,'mlabellocation',3)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',3)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','k','Fontweight','bold','FontSize',12,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.86 0.35 0.01 0.3],'TickDir','out','Ycolor','k','Xcolor','k',...
    'Fontweight','bold','FontSize',12);
%set(gcf,'Inverthardcopy','off');

author: CelsoReyes


text:
%% DEMO_stent_hexahedral_sweeping_02
% Below is a demonstration for:
%
% * Creating a hexahedral mesh for a vascular stent by sweeping allong a
% curve and copying over the segments.

%% Keywords
% * Sweeping, sweepLoft
% * Hexahedral mesh
% * stent, vascular
% * Exporting Abaqus, .inp

%%
clear; close all; clc;

%%
% plot settings

fontSize=25;
markerSize=10;
lineWidth=3;

%% Contol parameters

controlParameterSet.stentRadius=3; %The outer radius of the stent
controlParameterSet.numPeriodsWave=10; %The number of periods to use for a sinusoidal modulation
numStepsPeriod=100; %Number of sweeping steps allong a single period segment for sweeping
controlParameterSet.stentSectionHeight=0.1; %Height of the stent wire
controlParameterSet.stentSectionWidth=0.1; %Width of the stent wire
controlParameterSet.numStepsCircumference=(controlParameterSet.numPeriodsWave*numStepsPeriod)+1; %Number of sweeping steps across curve
controlParameterSet.overSampleFactorCurve=10; %Oversample factor curve
controlParameterSet.numSplitSteps_axial=1;
controlParameterSet.numSplitSteps_inward=1;
controlParameterSet.plotOn=0;
% controlParameterSet.waveAmplitude=0.9; %Amplitude of the sinusoidal modulation

sheetLayerThickness=0.025;
numStepsSheet=1;

%%

numSegments=8;
waveAmplitudes=0.6*ones(1,numSegments);
waveAmplitudes(2)=0.9;
offsetLevels=waveAmplitudes*2;

offSetTotal=0;

cFigure; hold on;
title('Stent hexahedral mesh','fontSize',fontSize);
colormap(gjet(4)); caxis([1 4]); icolorbar;
axisGeom;
camlight headlight;
drawnow;

E_stent_cell=cell(numSegments,1);
V_stent_cell=cell(numSegments,1);
E_sheet_cell=cell(numSegments,1);
V_sheet_cell=cell(numSegments,1);
for q=1:1:numSegments

    controlParameterSet.waveAmplitude=waveAmplitudes(q); %Amplitude of the sinusoidal modulation
    [E,V]=stentSegmentDesign(controlParameterSet);

    offSetTotal=offSetTotal+offsetLevels(q);
    V(:,3)=V(:,3)+offSetTotal;

    %%

    C=hexVol(E,V); %Get hexahedral element volumes

    [F,CF]=element2patch(E,C); %Create face data for plotting

    [indBoundary]=tesBoundary(F,V);
    faceMarker=ones(size(E,1),1)*(1:6); %The 6 face colors for the hexahedral faces
    faceMarker=faceMarker(:); %Force as a column
    Fb=F(indBoundary,:); %Select the boundary faces (which will exclude tops (1) and bottoms (2))
    faceBoundaryMarker=faceMarker(indBoundary,:)-2; %Get boundary colors and subtract 2 so they are 1-4

    %%
    gpatch(Fb,V,faceBoundaryMarker,'k',1);

    %%

    F_inner = Fb(faceBoundaryMarker==2,:);
    [edgesBoundaryInner]=patchBoundary(F_inner,V);

    edgesBottom=F_inner(:,[4 1]);
    edgesTop=F_inner(:,[2 3]);

    edgesBoundaryInnerTop=edgesBoundaryInner(all(ismember(edgesBoundaryInner,edgesTop),2),:);
    edgesBoundaryInnerBottom=edgesBoundaryInner(all(ismember(edgesBoundaryInner,edgesBottom),2),:);

    indCurveTop=edgeListToCurve(edgesBoundaryInnerTop);
    indCurveTop=flip(indCurveTop(1:end-1));
    indCurveBottom=edgeListToCurve(edgesBoundaryInnerBottom);
    indCurveBottom=indCurveBottom(1:end-1);

    plotV(V(indCurveTop(:),:),'r-','LineWidth',lineWidth);
    plotV(V(indCurveBottom(:),:),'b-','LineWidth',lineWidth);

    drawnow;

   %%
   if q==1
       [FQ,VQ]=patchCleanUnused(F_inner,V);
   else
       cPar.closeLoopOpt=1;
       cPar.patchType='quad';
       [Fq,Vq]=polyLoftLinear(V_curveTopPrevious,V(indCurveBottom(:),:),cPar);
%        gpatch(Fq,Vq,'rw','rw',1);
       [F_inner_clean,V_inner_clean]=patchCleanUnused(F_inner,V);
       [FQ,VQ]=joinElementSets({Fq,F_inner_clean},{Vq,V_inner_clean});
       [FQ,VQ]=mergeVertices(FQ,VQ);
%        gpatch(FQ,VQ,'rw','rw',1);
%        patchNormPlot(FQ,VQ);
   end

   [E_sheet,V_sheet,Fq1,Fq2]=quadThick(FQ,VQ,1,sheetLayerThickness,numStepsSheet);

   [F_sheet]=element2patch(E_sheet); %Create face data for plotting
   gpatch(F_sheet,V_sheet,'gw','gw',1);
   %        patchNormPlot(F_sheet,V_sheet);

   E_sheet_cell{q}=E_sheet;
   V_sheet_cell{q}=V_sheet;

   V_curveTopPrevious=V(indCurveTop(:),:);

   E_stent_cell{q}=E;
   V_stent_cell{q}=V;

end

%% Merge components

[E_stent,V_stent,C_stent]=joinElementSets(E_stent_cell,V_stent_cell);
[E_sheet,V_sheet,C_sheet]=joinElementSets(E_sheet_cell,V_sheet_cell);
[E,V,C]=joinElementSets({E_stent,E_sheet},{V_stent,V_sheet},{C_stent,C_sheet+max(C_stent)});
[E,V]=mergeVertices(E,V);

[F,CF]=element2patch(E,C); %Create face data for plotting

%%
cFigure; hold on;
title('Stent hexahedral mesh','fontSize',fontSize);
gpatch(F,V,CF,'none',1);
% patchNormPlot(F,V);
colormap gjet; icolorbar
axisGeom;
camlight headlight;
drawnow;

%% Export inp file
%
% elementStruct.E=E;
% elementStruct.E_ind=(1:size(E,1))';
% elementStruct.E_type='*ELEMENT, TYPE=C3D8, ELSET=PART-STENT';
% nodeStruct.N=V;
% nodeStruct.N_ind=(1:size(V,1))';
%
% pathName = fileparts(fileparts(mfilename('fullpath')));
% fileName=fullfile(pathName,'data','INP','stentMeshSheet.inp');
% export_INP(elementStruct,nodeStruct,fileName);

%% FUNCTIONS

function [E,V]=stentSegmentDesign(controlParameterSet)

%% parse input

stentRadius=controlParameterSet.stentRadius; %The outer radius of the stent
numPeriodsWave=controlParameterSet.numPeriodsWave; %The number of periods to use for a sinusoidal modulation
waveAmplitude=controlParameterSet.waveAmplitude; %Amplitude of the sinusoidal modulation
stentSectionHeight=controlParameterSet.stentSectionHeight; %Height of the stent wire
stentSectionWidth=controlParameterSet.stentSectionWidth; %Width of the stent wire
numStepsCircumference=controlParameterSet.numStepsCircumference; %Number of sweeping steps across curve
overSampleFactorCurve=controlParameterSet.overSampleFactorCurve; %Oversample factor curve
numSplitSteps_axial=controlParameterSet.numSplitSteps_axial;
numSplitSteps_inward=controlParameterSet.numSplitSteps_inward;
plotOn=controlParameterSet.plotOn;

%% plot settings
if plotOn==1
    fontSize=25;
    markerSize=10;
    lineWidth=1;
end

%% Build stent section
% The rectangular stent wire section is created here.

V_section=[-stentSectionWidth/2  stentSectionHeight/2 0; ...
    stentSectionWidth/2  stentSectionHeight/2 0; ...
    stentSectionWidth/2 -stentSectionHeight/2 0; ...
    -stentSectionWidth/2 -stentSectionHeight/2 0; ...
    ];

%%
% V=isualize stent section
if plotOn==1
    cFigure; hold on;
    title('Stent section','fontSize',fontSize);
    plotV(V_section,'b.-','lineWidth',lineWidth,'MarkerSize',markerSize);
    view(2); axis tight; axis equal; grid on; box on;
    set(gca,'fontSize',fontSize);
    drawnow;
end

%% Create guide curve
% The sweepLoft (see |HELP_sweepLoft|) is created here. First and angle
% based parameterization is created. Next this curve is evenly sample
% across the curve length (see |HELP_evenlySampleCurve|).
t=linspace(0,2*pi,numStepsCircumference*overSampleFactorCurve); %Angles
t=t(1:end-1); %Remove last point so it is not closed for resampling
x=stentRadius.*sin(t); %x-coordinates
y=stentRadius.*cos(t); %y-coordinates
z=waveAmplitude.*sin(numPeriodsWave*t); %z-coordinates
V_guide_curve=[x(:) y(:) z(:)]; %Collected curve nodes
[V_guide_curve] = evenlySampleCurve(V_guide_curve,numStepsCircumference-1,'pchip',1); %Resample curve evenly
V_guide_curve(end+1,:)=V_guide_curve(1,:); %Append start to end so it is a closed loop

%%
% Visualize guide curve
if plotOn==1
    cFigure; hold on;
    title('Stent guide curve','fontSize',fontSize);
    plotV(V_guide_curve,'k.-','lineWidth',lineWidth,'MarkerSize',markerSize);
    axisGeom;
    drawnow;
end

%% Position stent section at the start and end of the guide curve
% Next the section is translated and rotated so it is placed at the start
% of the guide curve such that the curve normal points allong the curve.

% Create rotation matrix
n3=vecnormalize(V_guide_curve(2,:)-V_guide_curve(1,:)); %Out of section normal direction z ish direction
[~,indMin]=min(dot(n3(ones(1,2),:),[1 0 0; 0 1 0],2)); %Get index most appropriate initial other axis
switch indMin
    case 1
        n1=[1 0 0]; %Initialized x direction
        n2=vecnormalize(cross(n3,n1)); %y ish direction
        n1=vecnormalize(cross(n2,n3)); %Proper x ish direction
        R=[n1; n2; n3]; %Rotation matrix
    case 2
        n2=[0 1 0]; %Initialized y direction
        n1=vecnormalize(cross(n2,n3)); %x ish direction
        n2=vecnormalize(cross(n3,n1)); %Proper y ish direction
        R=[n1; n2; n3]; %Rotation matrix
end

p1=V_guide_curve(1,:); %The start node
V_section=V_section*R; %Rotate the section
V_section=V_section+p1(ones(size(V_section,1),1),:); % Translate coordinate to start

%%
% Visualize guide curve

if plotOn==1
    cFigure; hold on;
    title('Stent section positioned on guide curve','fontSize',fontSize);
    plotV(V_guide_curve,'k-','lineWidth',1);
    plotV(V_section,'k.-','lineWidth',lineWidth,'MarkerSize',markerSize);
    quiverVec(p1,n1,1,'r');
    quiverVec(p1,n2,1,'g');
    quiverVec(p1,n3,1,'b');
    axisGeom;
    drawnow;
end

%% Sweeping section allong curve
% Normally |sweepLoft| produces patch data as an output (e.g. faces and
% vertices). However these outputs are supressed here and the coordinate
% mesh output is instead used to create a hexahedral mesh. See also |HELP_sweepLoft|

numTwist=0; %Number of additional twists of loft feature around guide curve
numStepsSweep=numStepsCircumference; %Number of steps for loft feature from sketch 1 to sketch 2
[~,~,~,S]=sweepLoft(V_section,V_section,n3,n3,V_guide_curve,numStepsSweep,numTwist,0);

%% Construct hexahedral mesh

X=S.X'; Y=S.Y'; Z=S.Z'; %Coordinate matrices
V=[X(:) Y(:) Z(:)]; %Create node list

F=reshape((1:1:size(V,1)),4,size(V,1)/4)'; %All top and bottom faces
E=[F(2:end,:) F(1:end-1,:)]; %The hexahedral elements
[E,V]=mergeVertices(E,V); %Merge nodes (start and end are not shared yet)

%% Refine mesh
% The swept mesh can be refined through slitting. The splitting can be
% homogeneous or only in a particular direction (see HELP_subHex|)
% Split method explanation:
% 1: Overall splitting in all directions
% 2: Split allong curve direction
% 3: Split axially
% 4: Splint inward

splitMethod=3;
nRefine=numSplitSteps_axial;
[E,V]=subHex(E,V,nRefine,splitMethod);

splitMethod=4;
nRefine=numSplitSteps_inward;
[E,V]=subHex(E,V,nRefine,splitMethod);

%%
% Visualize hexahedral mesh

if plotOn==1

    [F]=element2patch(E); %Create face data for plotting

    cFigure; hold on;
    title('Stent hexahedral mesh','fontSize',fontSize);
    plotV(V_guide_curve,'k-','lineWidth',3);
    gpatch(F,V,'gw','k',1);
    patchNormPlot(F,V);
    axisGeom;
    camlight headlight;
    drawnow;
end

%%

end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON