text:
function [streams,fileheader] = load_xdf(filename,varargin)
% Import an XDF file.
% [Streams,FileHeader] = load_xdf(Filename, Options...)
%
% This is a MATLAB importer for mult-stream XDF (Extensible Data Format) recordings. All
% information covered by the XDF 1.0 specification is imported, plus any additional meta-data
% associated with streams or with the container file itself.
%
% See https://github.com/sccn/xdf/ for more information on XDF.
%
% The function supports several further features, such as compressed XDF archives, robust
% time synchronization, support for breaks in the data, as well as some other defects.
%
% In:
%   Filename : name of the file to import (*.xdf or *.xdfz)
%
%   Options... : A list of optional name-value arguments for special use cases. The allowed names
%                are listed in the following:
%
%                Parameters that control various processing features:
%
%                'Verbose' : Whether to print verbose diagnostics. (default: false)
%
%                'HandleClockSynchronization' : Whether to enable clock synchronization based on
%                                               ClockOffset chunks. (default: true)
%
%                'HandleJitterRemoval' : Whether to perform jitter removal for regularly sampled
%                                        streams. (default: true)
%
%                'OnChunk' : Function that is called for each chunk of data as it
%                            is being retrieved from the file; the function is allowed to modify the
%                            data (for example, sub-sample it). The four input arguments are 1) the
%                            matrix of [#channels x #samples] values (either numeric or 2d cell
%                            array of strings), 2) the vector of unprocessed local time stamps (one
%                            per sample), 3) the info struct for the stream (same as the .info field
%                            in the final output, buth without the .effective_srate sub-field), and
%                            4) the scalar stream number (1-based integers). The three return values
%                            are 1) the (optionally modified) data, 2) the (optionally modified)
%                            time stamps, and 3) the (optionally modified) header (default: []).
%
%                'DisableVendorSpecifics' : Whether to perform certain vendor or system specific
%                                           operations. One example is the "BrainVision RDA" data
%                                           where it is necessary to pass the time stamps of a newly
%                                           introduced marker identifier channel on to the actual markers
%                                           in order to keep them perfectly in sync with the EEG data.
%                                           'DisableVendorSpecifics' takes the base names of certain
%                                           streams derived from the same source as a cell array of
%                                           strings (e.g. 'BrainVision RDA'). It is also possible to
%                                           disable vendor specifics altogether by providing the
%                                           value 'all'.
%
%                Parameters for advanced failure recovery in clock synchronization:
%
%                'HandleClockResets' : Whether the importer should check for potential resets of the
%                                      clock of a stream (e.g. computer restart during recording, or
%                                      hot-swap). Only useful if the recording system supports
%                                      recording under such circumstances. (default: true)
%
%                'ClockResetThresholdStds' : A clock reset must be accompanied by a ClockOffset
%                                            chunk being delayed by at least this many standard
%                                            deviations from the distribution. (default: 5)
%
%                'ClockResetThresholdSeconds' : A clock reset must be accompanied by a ClockOffset
%                                               chunk being delayed by at least this many seconds.
%                                               (default: 5)
%
%                'ClockResetThresholdOffsetStds' : A clock reset must be accompanied by a
%                                                  ClockOffset difference that lies at least this many
%                                                  standard deviations from the distribution. (default: 10)
%
%                'ClockResetThresholdOffsetSeconds' : A clock reset must be accompanied by a
%                                                     ClockOffset difference that is at least this
%                                                     many seconds away from the median. (default: 1)
%
%                'ClockResetMaxJitter' : Maximum tolerable jitter (in seconds of error) for clock
%                                        reset handling. (default: 5)
%
%                'CorrectStreamLags' : Apply lag correction described by timing spec sheets in the file.
%                                      (default: true)
%
%                Parameters for jitter removal in the presence of data breaks:
%
%                'JitterBreakThresholdSeconds' : An interruption in a regularly-sampled stream of at least this
%                                                many seconds will be considered as a potential break (if also
%                                                the BreakThresholdSamples is crossed) and multiple segments
%                                                will be returned. Default: 1
%
%                'JitterBreakThresholdSamples' : An interruption in a regularly-sampled stream of at least this
%                                                many samples will be considered as a potential break (if also
%                                                the BreakThresholdSeconds is crossed) and multiple segments
%                                                will be returned. Default: 500
%
%                Parameters for streams that can drop samples:
%
%                'FrameRateAccuracy' : How accurate the nominal frame rate is.
%                                      Used for can_drop_samples == true frames to determine
%                                      the maximum possible number of dropped frames. Default: 0.05%
%
% Out:
%   Streams : cell array of structs, one for each stream; the structs have the following content:
%             .time_series field: contains the stream's time series [#Channels x #Samples]
%                                 this matrix is of the type declared in .info.channel_format
%             .time_stamps field: contains the time stamps for each sample (synced across streams)
%
%             .info field: contains the meta-data of the stream (all values are strings)
%               .name: name of the stream
%               .type: content-type of the stream ('EEG','Events', ...)
%               .channel_format: value format ('int8','int16','int32','int64','float32','double64','string')
%               .nominal_srate: nominal sampling rate of the stream (as declared by the device);
%                               zero for streams with irregular sampling rate
%               .effective_srate: effective (measured) sampling rate of the stream, if regular
%                                 (otherwise omitted)
%               .desc: struct with any domain-specific meta-data declared for the stream; see
%                      www.xdf.org for the declared specifications
%
%             .segments field: struct array containing segment ranges for regularly sampled
%                              time series with breaks (not present if the stream is irregular)
%               .index_range: 1st and last index of the segment within the .time_series/.time_stamps
%                             arrays
%               .t_begin: time of the 1st sample in the segment, in seconds
%               .t_end: time of the last sample in the segment, in seconds
%               .duration: duration of the segment, in seconds
%               .num_samples: number of samples in the segment
%               .effective_srate: effective (i.e. measured) sampling rate within the segment
%
%   FileHeader : struct with file header contents in the .info field
%
% Examples:
%   % load the streams contained in a given XDF file
%   streams = load_xdf('C:\Recordings\myrecording.xdf')
%
% License:
%     This file is covered by the BSD license.
%
%     Copyright (c) 2012, Christian Kothe
%     Portions Copyright (c) 2010, Wouter Falkena
%     All rights reserved.
%
%     Redistribution and use in source and binary forms, with or without
%     modification, are permitted provided that the following conditions are
%     met:
%
%         * Redistributions of source code must retain the above copyright
%           notice, this list of conditions and the following disclaimer.
%         * Redistributions in binary form must reproduce the above copyright
%           notice, this list of conditions and the following disclaimer in
%           the documentation and/or other materials provided with the distribution
%
%     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%     AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%     IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%     ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%     LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%     CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%     SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%     INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%     CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%     ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%     POSSIBILITY OF SUCH DAMAGE.
%
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-04-22
%
%                                Contains portions of xml2struct Copyright (c) 2010, Wouter Falkena,
%                                ASTI, TUDelft, 21-08-2010
%
%                                Contains frame snapping routine to handle dropped video frames by
%                                Matthew Grivich.
%
%                                version 1.13
LIBVERSION = '1.13';
% check inputs
opts = cell2struct(varargin(2:2:end),varargin(1:2:end),2);
if ~isfield(opts,'OnChunk')
    opts.OnChunk = []; end
if ~isfield(opts,'Verbose')
    opts.Verbose = false; end
if ~isfield(opts,'HandleClockSynchronization')
    opts.HandleClockSynchronization = true; end
if ~isfield(opts,'HandleClockResets')
    opts.HandleClockResets = true; end
if ~isfield(opts,'HandleJitterRemoval')
    opts.HandleJitterRemoval = true; end
if ~isfield(opts,'JitterBreakThresholdSeconds')
    opts.JitterBreakThresholdSeconds = 1; end
if ~isfield(opts,'JitterBreakThresholdSamples')
    opts.JitterBreakThresholdSamples = 500; end
if ~isfield(opts,'ClockResetThresholdSeconds')
    opts.ClockResetThresholdSeconds = 5; end
if ~isfield(opts,'ClockResetThresholdStds')
    opts.ClockResetThresholdStds = 5; end
if ~isfield(opts,'ClockResetThresholdOffsetSeconds')
    opts.ClockResetThresholdOffsetSeconds = 1; end
if ~isfield(opts,'ClockResetThresholdOffsetStds')
    opts.ClockResetThresholdOffsetStds = 10; end
if ~isfield(opts,'WinsorThreshold')
    opts.WinsorThreshold = 0.0001; end
if ~isfield(opts,'ClockResetMaxJitter')
    opts.ClockResetMaxJitter = 5; end
if ~isfield(opts,'DisableVendorSpecifics')
    opts.DisableVendorSpecifics = {}; end
if ~isfield(opts,'CorrectStreamLags')
    opts.CorrectStreamLags = true; end
if ~isfield(opts,'FrameRateAccuracy')
    opts.FrameRateAccuracy = .05; end
if ~exist(filename,'file')
    error(['The file "' filename '" does not exist.']); end

% uncompress if necessary (note: "bonus" feature, not part of the XDF 1.0 spec)
[p,n,x] = fileparts(filename);
if strcmp(x,'.xdfz')
    % idea for this type of approach by Michael Kleder
    import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
    src = java.io.FileInputStream(filename);
    flt = java.util.zip.InflaterInputStream(src);
    filename = [p filesep n '_temp_uncompressed' x];
    dst = java.io.FileOutputStream(filename);
    copier = InterruptibleStreamCopier.getInterruptibleStreamCopier;
    copier.copyStream(flt,dst);
    dst.close();
    src.close();
end

streams = {};                                   % cell array of returned streams (in the order of appearance in the file)
idmap = sparse(2^31-1,1);                       % remaps stream id's onto indices in streams
temp = struct();                                % struct array of temporary per-stream information
fileheader = struct();                          % the file header
f = fopen(filename,'r','ieee-le.l64');          % file handle
closer = onCleanup(@()close_file(f,filename));  % object that closes the file when the function exits

% there is a fast C mex file for the inner loop, but it's
% not necessarily available for every platform
have_mex = exist('load_xdf_innerloop','file');
if ~have_mex
    if opts.Verbose
        disp(['NOTE: apparently you are missing a compiled binary version of the inner loop code.',...
            ' Attempting to download...']);
    end

    fname = ['load_xdf_innerloop.' mexext];
    mex_url = ['https://github.com/sccn/xdf/releases/download/v',...
        LIBVERSION, '/', fname];
    [this_path, this_name, this_ext] = fileparts(mfilename('fullpath'));
    try
        have_mex = true;
        websave(fullfile(this_path, fname), mex_url);
    catch ME
        if opts.Verbose
            disp(['Unable to download the compiled binary version for your platform.',...
                ' Using the slow MATLAB code instead.']);
        end
        have_mex = false;
        %rethrow(ME);
    end
end

% ======================
% === parse the file ===
% ======================

% read [MagicCode]
if ~strcmp(fread(f,4,'*char')','XDF:')
    error(['This is not a valid XDF file (' filename ').']); end

% for each chunk...

if opts.Verbose; fprintf('Now reading from %s ...', filename); end;
while 1
    % read [NumLengthBytes], [Length]
    len = double(read_varlen_int(f));
    if ~len
        break; end
    % read [Tag]
    switch fread(f,1,'uint16')
        case 3 % read [Samples] chunk
            try
                % read [StreamId]
                id = idmap(fread(f,1,'uint32'));
                if have_mex
                    % read the chunk data at once
                    data = fread(f,len-6,'*uint8');
                    % run the mex kernel
                    [values,timestamps] = load_xdf_innerloop(data, temp(id).chns, temp(id).readfmt, temp(id).sampling_interval, temp(id).last_timestamp);
                    temp(id).last_timestamp = timestamps(end);
                else % fallback MATLAB implementation
                    % read [NumSampleBytes], [NumSamples]
                    num = read_varlen_int(f);
                    % allocate space
                    timestamps = zeros(1,num);
                    if strcmp(temp(id).readfmt,'*string')
                        values = cell(temp(id).chns,num);
                    else
                        values = zeros(temp(id).chns,num);
                    end
                    % for each sample...
                    for s=1:num
                        % read or deduce time stamp
                        if fread(f,1,'*uint8')
                            timestamps(s) = fread(f,1,'double');
                        else
                            timestamps(s) = temp(id).last_timestamp + temp(id).sampling_interval;
                        end
                        % read the values
                        if strcmp(temp(id).readfmt,'*string')
                            for v = 1:size(values,1)
                                values{v,s} = fread(f,double(read_varlen_int(f)),'*char')'; end
                        else
                            values(:,s) = fread(f,size(values,1),temp(id).readfmt);
                        end
                        temp(id).last_timestamp = timestamps(s);
                    end
                end
                % optionally send through the OnChunk function
                if ~isempty(opts.OnChunk)
                    [values,timestamps,streams{id}] = opts.OnChunk(values,timestamps,streams{id},id); end %#ok<*AGROW>
                % append to the time series...
                temp(id).time_series{end+1} = values;
                temp(id).time_stamps{end+1} = timestamps;
            catch e
                % an error occurred (perhaps a chopped-off file): emit a warning
                % and return the file up to this point
                warning(e.identifier,e.message);
                break;
            end
        case 2 % read [StreamHeader] chunk
            % read [StreamId]
            streamid = fread(f,1,'uint32');
            id = length(streams)+1;
            idmap(streamid) = id; %#ok<SPRIX>
            % read [Content]
            header = parse_xml_struct(fread(f,len-6,'*char')');
            streams{id} = header;
            if opts.Verbose
                fprintf(['  found stream ' header.info.name '\n']); end
            % generate a few temporary fields
            temp(id).chns = str2num(header.info.channel_count); %#ok<*ST2NM>
            temp(id).srate = str2num(header.info.nominal_srate);
            temp(id).last_timestamp = 0;
            temp(id).time_series = {};
            temp(id).time_stamps = {};
            temp(id).clock_times = [];
            temp(id).clock_values = [];
            if temp(id).srate > 0
                temp(id).sampling_interval = 1/temp(id).srate;
            else
                temp(id).sampling_interval = 0;
            end
            % fread parsing format for data values
            temp(id).readfmt = ['*' header.info.channel_format];
            if strcmp(temp(id).readfmt,'*double64') && ~have_mex
                temp(id).readfmt = '*double'; end % for fread()
        case 6 % read [StreamFooter] chunk
            % read [StreamId]
            id = idmap(fread(f,1,'uint32'));
            % read [Content]
            footer = parse_xml_struct(fread(f,len-6,'*char')');
            streams{id} = hlp_superimposedata(footer,streams{id});
        case 1 % read [FileHeader] chunk
            fileheader = parse_xml_struct(fread(f,len-2,'*char')');
        case 4 % read [ClockOffset] chunk
            % read [StreamId]
            id = idmap(fread(f,1,'uint32'));
            % read [CollectionTime]
            temp(id).clock_times(end+1) = fread(f,1,'double');
            % read [OffsetValue]
            temp(id).clock_values(end+1) = fread(f,1,'double');
        otherwise
            % skip other chunk types (Boundary, ...)
            fread(f,len-2,'*uint8');
    end
end

% concatenate the signal across chunks
for k=1:length(temp)
    try
        temp(k).time_series = [temp(k).time_series{:}];
        temp(k).time_stamps = [temp(k).time_stamps{:}];
    catch e
        disp(['Could not concatenate time series for stream ' streams{k}.info.name '; skipping.']);
        disp(['Reason: ' e.message]);
        temp(k).time_series = [];
        temp(k).time_stamps = [];
    end
end

% ===================================================================
% === perform (fault-tolerant) clock synchronization if requested ===
% ===================================================================

if opts.HandleClockSynchronization
    if opts.Verbose
        disp('  performing clock synchronization...'); end
    for k=1:length(temp)
        if ~isempty(temp(k).time_stamps)
            try
                clock_times = temp(k).clock_times;
                clock_values = temp(k).clock_values;
                if isempty(clock_times)
                    error('No clock offset values present.'); end
            catch
                disp(['No clock offsets were available for stream "' streams{k}.info.name '"']);
                continue;
            end

            % detect clock resets (e.g., computer restarts during recording) if requested
            % this is only for cases where "everything goes wrong" during recording
            % note that this is a fancy feature that is not needed for normal XDF compliance
            if opts.HandleClockResets
                % first detect potential breaks in the synchronization data; this is only necessary when the
                % importer should be able to deal with recordings where the computer that served a stream
                % was restarted or hot-swapped during an ongoing recording, or the clock was reset otherwise
                time_diff = diff(clock_times);
                value_diff = abs(diff(clock_values));
                % points where a glitch in the timing of successive clock measurements happened
                time_glitch = (time_diff < 0 | (((time_diff - median(time_diff)) ./ median(abs(time_diff-median(time_diff)))) > opts.ClockResetThresholdStds & ...
                    ((time_diff - median(time_diff)) > opts.ClockResetThresholdSeconds)));
                % points where a glitch in successive clock value estimates happened
                value_glitch = (value_diff - median(value_diff)) ./ median(abs(value_diff-median(value_diff))) > opts.ClockResetThresholdOffsetStds & ...
                    (value_diff - median(value_diff)) > opts.ClockResetThresholdOffsetSeconds;
                % points where both a time glitch and a value glitch co-occur are treated as resets
                resets_at = time_glitch & value_glitch;
                % determine the [begin,end] index ranges between resets
                if any(resets_at)
                    tmp = find(resets_at)';
                    tmp = [tmp tmp+1]';
                    tmp = [1 tmp(:)' length(resets_at)];
                    ranges = num2cell(reshape(tmp,2,[])',2);
                    if opts.Verbose
                        disp(['  found ' num2str(nnz(resets_at)) ' clock resets in stream ' streams{k}.info.name '.']); end
                else
                    ranges = {[1,length(clock_times)]};
                end
            else
                % otherwise we just assume that there are no clock resets
                ranges = {[1,length(clock_times)]};
            end

            % calculate clock offset mappings for each data range
            mappings = {};
            for r=1:length(ranges)
                idx = ranges{r};
                if idx(1) ~= idx(2)
                    % to accomodate the Winsorizing threshold (in seconds) we rescale the data (robust_fit sets it to 1 unit)
                    mappings{r} = robust_fit([ones(idx(2)-idx(1)+1,1) clock_times(idx(1):idx(2))']/opts.WinsorThreshold, clock_values(idx(1):idx(2))'/opts.WinsorThreshold);
                else
                    mappings{r} = [clock_values(idx(1)) 0]; % just one measurement
                end
            end

            if length(ranges) == 1
                % apply the correction to all time stamps
                temp(k).time_stamps = temp(k).time_stamps + (mappings{1}(1) + mappings{1}(2)*temp(k).time_stamps);
            else
                % if there are data segments measured with different clocks we need to
                % determine, for any time stamp lying between two segments, to which of the segments it belongs
                clock_segments = zeros(size(temp(k).time_stamps));  % the segment index to which each stamp belongs
                begin_of_segment = 1;                               % first index into time stamps that belongs to the current segment
                end_of_segment = NaN; %#ok<NASGU>                   % last index into time stamps that belongs to the current segment
                for r=1:length(ranges)-1
                    cur_end_time = clock_times(ranges{r}(2));       % time at which the current segment ends
                    next_begin_time = clock_times(ranges{r+1}(1));  % time at which the next segment begins
                    % get the data that is not yet processed
                    remaining_indices = begin_of_segment:length(temp(k).time_stamps);
                    if isempty(remaining_indices)
                        break; end
                    remaining_data = temp(k).time_stamps(remaining_indices);
                    if next_begin_time > cur_end_time
                        % clock jumps forward: the end of the segment is where the data time stamps
                        % lie closer to the next segment than the current in time
                        end_of_segment = remaining_indices(min(find([abs(remaining_data-cur_end_time) > abs(remaining_data-next_begin_time),true],1)-1,length(remaining_indices)));
                    else
                        % clock jumps backward: the end of the segment is where the data time stamps
                        % jump back by more than the max conceivable jitter (as any negative delta is jitter)
                        end_of_segment = remaining_indices(min(find([diff(remaining_data) < -opts.ClockResetMaxJitter,true],1),length(remaining_indices)));
                    end
                    % assign the segment of data points to the current range
                    % go to next segment
                    clock_segments(begin_of_segment:end_of_segment) = r;
                    begin_of_segment = end_of_segment+1;
                end
                % assign all remaining time stamps to the last segment
                clock_segments(begin_of_segment:end) = length(ranges);
                % apply corrections on a per-segment basis
                for r=1:length(ranges)
                    temp(k).time_stamps(clock_segments==r) = temp(k).time_stamps(clock_segments==r) + (mappings{r}(1) + mappings{r}(2)*temp(k).time_stamps(clock_segments==r)); end
            end
        end
    end
end

% ===========================================
% === perform jitter removal if requested ===
% ===========================================
if opts.HandleJitterRemoval
    % jitter removal is a bonus feature that yields linearly increasing timestamps from data
    % where samples had been time stamped with some jitter (e.g., due to operating system
    % delays)
    if opts.Verbose
        disp('  performing jitter removal...'); end

    for k=1:length(temp)

        if ~isempty(temp(k).time_stamps) && temp(k).srate

            if isfield(streams{k}.info.desc, 'synchronization') && ...
                isfield(streams{k}.info.desc.synchronization, 'can_drop_samples') && ...
                strcmp(streams{k}.info.desc.synchronization.can_drop_samples, 'true')
                temp(k).time_stamps = droppedFramesCorrection(temp(k).time_stamps,temp(k).srate, opts.FrameRateAccuracy);
           else

                % identify breaks in the data
                diffs = diff(temp(k).time_stamps);
                breaks_at = diffs > max(opts.JitterBreakThresholdSeconds,opts.JitterBreakThresholdSamples*temp(k).sampling_interval);
                if any(breaks_at)
                    % turn the break mask into a cell array of [begin,end] index ranges
                    tmp = find(breaks_at)';
                    tmp = [tmp tmp+1]';
                    tmp = [1 tmp(:)' length(breaks_at)];
                    ranges = num2cell(reshape(tmp,2,[])',2);
                    if opts.Verbose
                        disp(['  found ' num2str(nnz(breaks_at)) ' data breaks in stream ' streams{k}.info.name '.']); end
                else
                    ranges = {[1,length(temp(k).time_stamps)]};
                end

                % process each segment separately
                segments = repmat(struct(),1,length(ranges));
                for r=1:length(ranges)
                    range = ranges{r};
                    segments(r).num_samples = range(2)-range(1)+1;
                    segments(r).index_range = range;
                    if segments(r).num_samples > 0
                        indices = segments(r).index_range(1):segments(r).index_range(2);
                        % regress out the jitter
                        mapping = temp(k).time_stamps(indices) / [ones(1,length(indices)); indices];
                        temp(k).time_stamps(indices) = mapping(1) + mapping(2) * indices;
                    end
                    % calculate some other meta-data about the segments
                    segments(r).t_begin = temp(k).time_stamps(range(1));
                    segments(r).t_end = temp(k).time_stamps(range(2));
                    segments(r).duration = segments(r).t_end - segments(r).t_begin;
                    segments(r).effective_srate = (segments(r).num_samples - 1) / segments(r).duration;
                end

                % calculate the weighted mean sampling rate over all segments
                temp(k).effective_rate = sum(bsxfun(@times,[segments.effective_srate],[segments.num_samples]/sum([segments.num_samples])));

                % transfer the information into the output structs
                streams{k}.info.effective_srate = temp(k).effective_rate;
                streams{k}.segments = segments;
            end
        end
    end
else
    % calculate effective sampling rate
    for k=1:length(temp)
        temp(k).effective_srate = (length(temp(k).time_stamps) - 1) / (temp(k).time_stamps(end) - temp(k).time_stamps(1)); end
end

% copy the information into the output
for k=1:length(temp)
    offset_mean = 0;
    if opts.CorrectStreamLags && ...
        isfield(streams{k}.info.desc, 'synchronization') && ...
        isfield(streams{k}.info.desc.synchronization, 'offset_mean')
                offset_mean = str2num(streams{k}.info.desc.synchronization.offset_mean);
    end
    streams{k}.time_series = temp(k).time_series;
    streams{k}.time_stamps = temp(k).time_stamps - offset_mean;
end

% =========================================
% === peform vendor specific operations ===
% =========================================

if ~any(strcmp('all',opts.DisableVendorSpecifics))

    % BrainVision RDA
    targetName = 'BrainVision RDA';
    if ~any(strcmp(opts.DisableVendorSpecifics,targetName))

        % find a target EEG stream...
        for k=1:length(streams)

            if strcmp(streams{k}.info.name,targetName) % Is a BrainVision RDA stream?
                mkChan = [];
                for iChan = 1:length( streams{ k }.info.desc.channels.channel ) % Find marker index channel (any channel, not necessary last)
                    if strcmp( streams{ k }.info.desc.channels.channel{ iChan }.label, 'MkIdx' ) && strcmp( streams{ k }.info.desc.channels.channel{ iChan }.type, 'Marker' ) && strcmp( streams{ k }.info.desc.channels.channel{ iChan }.unit, 'Counts (decimal)' )
                        mkChan = iChan;
                        break % Only one marker channel expected
                    end
                end
                if ~isempty( mkChan ) % Has a marker index channel?
                    for m = 1:length( streams ) % find a corresponding indexed marker stream...
                        if strcmp( streams{ m }.info.name, [ targetName ' Markers' ] ) && strcmp( streams{ m }.info.hostname, streams{ k }.info.hostname ) && strncmp( streams{ m }.info.source_id, streams{ k }.info.source_id, length( streams{ k }.info.source_id ) )
                            if opts.Verbose
                                disp( [ '  performing ', targetName, ' specific tasks for stream ', num2str( k ), '...' ] );
                            end
                            streams = ProcessBVRDAindexedMarkers( streams, k, m, mkChan );
                        end
                    end

					% Remove marker index channel
					streams{ k }.time_series( mkChan, : ) = [];
					streams{ k }.info.desc.channels.channel( mkChan ) = [];

					% Decrement channel count by 1
					streams{ k }.info.channel_count = num2str( str2num( streams{ k }.info.channel_count ) - 1 );

                end

            end

        end

    end

end

end

% ========================
% === helper functions ===
% ========================

% read a variable-length integer
function num = read_varlen_int(f)
try
    switch fread(f,1,'*uint8')
        case 1
            num = fread(f,1,'*uint8');
        case 4
            num = fread(f,1,'*uint32');
        case 8
            num = fread(f,1,'*uint64');
        otherwise
            error('Invalid variable-length integer encountered.');
    end
catch %#ok<*CTCH>
    num = 0;
end
end

% close the file and delete temporary data
function close_file(f,filename)
fclose(f);
if strfind(filename,'_temp_uncompressed.xdf')
    delete(filename); end
end

% parse a simplified (attribute-free) subset of XML into a MATLAB struct
function result = parse_xml_struct(str)
import org.xml.sax.InputSource
import javax.xml.parsers.*
import java.io.*
tmp = InputSource();
tmp.setCharacterStream(StringReader(str));
result = parseChildNodes(xmlread(tmp));

% this is part of xml2struct (slightly simplified)
    function [children,ptext] = parseChildNodes(theNode)
        % Recurse over node children.
        children = struct;
        ptext = [];
        if theNode.hasChildNodes
            childNodes = theNode.getChildNodes;
            numChildNodes = childNodes.getLength;
            for count = 1:numChildNodes
                theChild = childNodes.item(count-1);
                [text,name,childs] = getNodeData(theChild);
                if (~strcmp(name,'#text') && ~strcmp(name,'#comment'))
                    if (isfield(children,name))
                        if (~iscell(children.(name)))
                            children.(name) = {children.(name)}; end
                        index = length(children.(name))+1;
                        children.(name){index} = childs;
                        if(~isempty(text))
                            children.(name){index} = text; end
                    else
                        children.(name) = childs;
                        if(~isempty(text))
                            children.(name) = text; end
                    end
                elseif (strcmp(name,'#text'))
                    if (~isempty(regexprep(text,'[\s]*','')))
                        if (isempty(ptext))
                            ptext = text;
                        else
                            ptext = [ptext text];
                        end
                    end
                end
            end
        end
    end

% this is part of xml2struct (slightly simplified)
    function [text,name,childs] = getNodeData(theNode)
        % Create structure of node info.
        name = char(theNode.getNodeName);
        if ~isvarname(name)
            name = regexprep(name,'[-]','_dash_');
            name = regexprep(name,'[:]','_colon_');
            name = regexprep(name,'[.]','_dot_');
        end
        [childs,text] = parseChildNodes(theNode);
        if (isempty(fieldnames(childs)))
            try
                text = char(theNode.getData);
            catch
            end
        end
    end
end

function x = robust_fit(A,y,rho,iters)
% Perform a robust linear regression using the Huber loss function.
% x = robust_fit(A,y,rho,iters)
%
% Input:
%   A : design matrix
%   y : target variable
%   rho : augmented Lagrangian variable (default: 1)
%   iters : number of iterations to perform (default: 1000)
%
% Output:
%   x : solution for x
%
% Notes:
%   solves the following problem via ADMM for x:
%     minimize 1/2*sum(huber(A*x - y))
%
% Based on the ADMM Matlab codes also found at:
%   http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-03-04

if ~exist('rho','var')
    rho = 1; end
if ~exist('iters','var')
    iters = 1000; end
Aty = A'*y;
L = sparse(chol(A'*A,'lower')); U = L';
z = zeros(size(y)); u = z;
for k = 1:iters
    x = U \ (L \ (Aty + A'*(z - u)));
    d = A*x - y + u;
    z = rho/(1+rho)*d + 1/(1+rho)*max(0,(1-(1+1/rho)./abs(d))).*d;
    u = d - z;
end
end

function res = hlp_superimposedata(varargin)
% Merge multiple partially populated data structures into one fully populated one.
% Result = hlp_superimposedata(Data1, Data2, Data3, ...)
%
% The function is applicable when you have cell arrays or structs/struct arrays with non-overlapping
% patterns of non-empty entries, where all entries should be merged into a single data structure
% which retains their original positions. If entries exist in multiple data structures at the same
% location, entries of later items will be ignored (i.e. earlier data structures take precedence).
%
% In:
%   DataK : a data structure that should be super-imposed with the others to form a single data
%           structure
%
% Out:
%   Result : the resulting data structure
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-08-19

% first, compactify the data by removing the empty items
compact = varargin(~cellfun('isempty',varargin));
% start with the last data structure, then merge the remaining data structures into it (in reverse
% order as this avoids having to grow arrays incrementally in typical cases)
res = compact{end};
for k=length(compact)-1:-1:1
    res = merge(res,compact{k}); end
end

function A = merge(A,B)
% merge data structures A and B
if iscell(A) && iscell(B)
    % make sure that both have the same number of dimensions
    if ndims(A) > ndims(B)
        B = grow_cell(B,size(A));
    elseif ndims(A) < ndims(B)
        A = grow_cell(A,size(B));
    end
    % make sure that both have the same size
    if all(size(B)==size(A))
        % we're fine
    elseif all(size(B)>=size(A))
        % A is a minor of B: grow A
        A = grow_cell(A,size(B));
    elseif all(size(A)>=size(B))
        % B is a minor of A: grow B
        B = grow_cell(B,size(A));
    else
        % A and B have mixed sizes... grow both as necessary
        M = max(size(A),size(B));
        A = grow_cell(A,M);
        B = grow_cell(B,M);
    end
    % find all non-empty elements in B
    idx = find(~cellfun(@(x)isequal(x,[]),B));
    if ~isempty(idx)
        % check if any of these is occupied in A
        clean = cellfun('isempty',A(idx));
        if ~all(clean)
            % merge all conflicting items recursively
            conflicts = idx(~clean);
            for k=conflicts(:)'
                A{k} = merge(A{k},B{k}); end
            % and transfer the rest
            if any(clean)
                A(idx(clean)) = B(idx(clean)); end
        else
            % transfer all to A
            A(idx) = B(idx);
        end
    end
elseif isstruct(A) && isstruct(B)
    % first make sure that both have the same fields
    fnA = fieldnames(A);
    fnB = fieldnames(B);
    if isequal(fnA,fnB)
        % we're fine
    elseif isequal(sort(fnA),sort(fnB))
        % order doesn't match -- impose A's order on B
        B = orderfields(B,fnA);
    elseif isempty(setdiff(fnA,fnB))
        % B has a superset of A's fields: add the remaining fields to A, and order them according to B
        remaining = setdiff(fnB,fnA);
        for fn = remaining'
            A(1).(fn{1}) = []; end
        A = orderfields(A,fnB);
    elseif isempty(setdiff(fnB,fnA))
        % A has a superset of B's fields: add the remaining fields to B, and order them according to A
        remaining = setdiff(fnA,fnB);
        for fn = remaining'
            B(1).(fn{1}) = []; end
        B = orderfields(B,fnA);
    else
        % A and B have incommensurable fields; add B's fields to A's fields, add A's fields to B's
        % and order according to A's fields
        remainingB = setdiff(fnB,fnA);
        for fn = remainingB'
            A(1).(fn{1}) = []; end
        remainingA = setdiff(fnA,fnB);
        for fn = remainingA'
            B(1).(fn{1}) = []; end
        B = orderfields(B,A);
    end
    % that being established, convert them to cell arrays, merge their cell arrays, and convert back to structs
    merged = merge(struct2cell(A),struct2cell(B));
    A = cell2struct(merged,fieldnames(A),1);
elseif isstruct(A) && ~isstruct(B)
    if ~isempty(B)
        error('One of the sub-items is a struct, and the other one is of a non-struct type.');
    else
        % we retain A
    end
elseif isstruct(B) && ~isstruct(A)
    if ~isempty(A)
        error('One of the sub-items is a struct, and the other one is of a non-struct type.');
    else
        % we retain B
        A = B;
    end
elseif iscell(A) && ~iscell(B)
    if ~isempty(B)
        error('One of the sub-items is a cell array, and the other one is of a non-cell type.');
    else
        % we retain A
    end
elseif iscell(B) && ~iscell(A)
    if ~isempty(A)
        error('One of the sub-items is a cell array, and the other one is of a non-cell type.');
    else
        % we retain B
        A = B;
    end
elseif isempty(A) && ~isempty(B)
    % we retain B
    A = B;
elseif isempty(B) && ~isempty(A)
    % we retain A
elseif ~isequal(A,B)
    % we retain A and warn about dropping B
    disp('Two non-empty (and non-identical) sub-elements occupied the same index; one was dropped. This warning will only be displayed once.');
end
end

function C = grow_cell(C,idx)
% grow a cell array to accomodate a particular index
% (assuming that this index is not contained in the cell array yet)
tmp = sprintf('%i,',idx);
eval(['C{' tmp(1:end-1) '} = [];']);
end

%function written by Matthew Grivich to handle case where frame rate is
%consistent but with dropped events, like a video camera.
function frameTimesModeled = droppedFramesCorrection(frameTimes, nominalFrameRate, frameRateAccuracy)

  %  figure
  %  plot(frameTimes(1:end-1), frameTimes(2:end)-frameTimes(1:end-1));
  %  xlabel('Frame Time (s)');
  %  ylabel('Frame Interval (s)');

    nFrames = length(frameTimes); %nFrames shown, does not included dropped.

    %calculates the maximum conceivable number of dropped frames, given the
    %nominal frame rate and the expected frame rate accuracy.
    maxDropped = round((frameTimes(end)-frameTimes(1))*nominalFrameRate*(1+frameRateAccuracy) - nFrames);
    stds = zeros(0, 1); %initialize array of zero length
    for iteration = 1:maxDropped
        interval = (frameTimes(end)-frameTimes(1))/(nFrames-1+iteration-1);
   %     fprintf('frequency: %1.20f\n',1/interval);
        frameNumbers = zeros(1,length(frameTimes));

        for i=1:length(frameNumbers)

            frameNumbers(i) =  round((frameTimes(i)-frameTimes(1))/interval)+1;

        end
        %remove zero frame intervals
        for i=length(frameNumbers)-1:-1:1
            if(frameNumbers(i) == frameNumbers(i+1))
                frameNumbers(i) = frameNumbers(i) -1;
            end
        end

        pf = polyfit(frameNumbers, frameTimes,1);
        %text(0,.0004, sprintf('interval: %1.20f\n', interval));
    %    fprintf('interval after fit: %1.20f\n',pf(1));

        frameTimesModeled = polyval(pf,frameNumbers);
  %     plot(stds)
      %  plot(frameTimesModeled, frameTimesModeled - frameTimes);
  %      text(0,.02, sprintf('dropped: %d\n', iteration-1));
  %      xlabel('Frame Time (s)')
  %      ylabel('Frame Time Modeled - Frame Time (s)');
  %     pause(0.001);

        stds(iteration) = std(frameTimesModeled-frameTimes);
        %if (mean(stds) - min(stds))/std(stds) > 10
        %    break;
        %end

    end

 %   figure
 %   plot(stds);
    dropped = find(stds==min(stds)) - 1;

    interval = (frameTimes(end)-frameTimes(1))/(nFrames-1+dropped);
   % fprintf('interval: %1.20f\n',interval);
    frameNumbers = zeros(1,length(frameTimes));

    %Find closest frame.
    for i=1:length(frameNumbers)
        frameNumbers(i) =  round((frameTimes(i)-frameTimes(1))/interval)+1;
    end
    %remove zero frame intervals
    for i=length(frameNumbers)-1:-1:1
        if(frameNumbers(i) >= frameNumbers(i+1))
            frameNumbers(i) = frameNumbers(i+1) -1;
        end
    end

    pf = polyfit(frameNumbers, frameTimes,1);
    %text(0,.0004, sprintf('interval: %1.20f\n', interval));
%    fprintf('interval after fit: %1.20f\n',pf(1));

    frameTimesModeled = polyval(pf,frameNumbers);
%figure
%        plot(frameTimesModeled, smooth(frameTimesModeled - frameTimes,21));
%        xlabel('Frame Time (s)')
%        ylabel('Frame Time Modeled - Frame Time (s)');
 %      pause(0.5);

%    fprintf('Interval: %1.20f\n',pf(1));
%    fprintf('Dropped Frames: %d\n', dropped);

   % figure
   % plot(frameTimes(1:end-1), pf(1)*(frameNumbers(2:end)-frameNumbers(1:end-1)));
   % xlabel('Frame Time (s)');
   % ylabel('Frame Interval (s)');

end

function streams = ProcessBVRDAindexedMarkers( streams, dataStream, mkStream, mkChan )

clearMarkers = [];

for iMrk = 1:length( streams{ mkStream }.time_series )

    % Decode marker
    MrkInfo = regexp( streams{ mkStream }.time_series{ iMrk }, 'mk(?<idx>\d+)=(?<str>.*)', 'names' );

    if ~isempty( MrkInfo.idx )

        % Find corresponding sample in marker index channel
        lat = find( streams{ dataStream }.time_series( mkChan, : ) == str2double( MrkInfo.idx ) );
        offset = streams{ mkStream }.time_stamps( iMrk ) - streams{ dataStream }.time_stamps( lat );

        % Is the index unique (overflow)?
        if length( lat ) > 1
            [ minOffset, minOffsetIdx ] = min( abs( offset ) ); %#ok<ASGLU>
            lat = lat( minOffsetIdx );
            offset = offset( minOffsetIdx );
        end

        % Sanity check
        if offset > 10
            warning( 'Time stamp difference between indexed marker %s (%.3f s) and corresponding sample in marker channel (%.3f s) exceeding threshold.\n', MrkInfo.idx, streams{ mkStream }.time_stamps( iMrk ), streams{ dataStream }.time_stamps( lat ) )
        end

        % Copy time stamp and rewrite marker
        if ~isempty( lat )
            streams{ mkStream }.time_stamps( iMrk ) = streams{ dataStream }.time_stamps( lat );
            streams{ mkStream }.time_series{ iMrk } = MrkInfo.str;
        else
            warning( 'No corresponding sample found in marker channel for indexed marker %s. Removing...', MrkInfo.idx )
            clearMarkers = [ clearMarkers iMrk ]; %#ok<AGROW>
        end

    end

end

% Remove markers without corresponding marker channel sample
streams{ mkStream }.time_stamps( clearMarkers ) = [];
streams{ mkStream }.time_series( clearMarkers ) = [];

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [p1Clip, p2Clip]=clipLineSegment2Rect(rectMaxMin,p1,p2)
%%CLIPLINESEGMENT2RECT Given a 2D rectangle (axis-aligned) and given the
%           two endpoints of a line segment clip the line segment to the
%           rectangle. If the segment and the rectangle do not overlap,
%           then an empty matrix is returned.
%
%INPUTS: rectMaxMin A length-4 array specifying the maximum and minimum
%                   bounds in x and y. rectMaxMin=[xMin;xMax;yMin;yMax].
%                   The minimum values must be <=the maximum values.
%                p1 A length-2 array specifying the start of the line
%                   segment.
%                p2 A length-2 array specifying the end of the line
%                   segment.
%
%OUTPUTS: p1Clip, p2Clip These are the same p1 and p2 but have been clipped to the
%                start and end of the line segment. If nothing clips, the
%                empty matrices are returned.
%
%The algorithm of [1] is used, without speed improvements that are possible
%from direct substitution, which is mentioned in Section IX of [1]. The
%performance of the algorithm is discussed in [2].
%
%EXAMPLE:
% rect=[-1;1;-1;1];
% numLines=6;
% p1=zeros(2,numLines);
% %Lines entirely inside the box.
% p1(:,1)=[0;0]; p2(:,1)=[0.8;0.8];
% p1(:,2)=[0.6;0]; p2(:,2)=[-0.8;0.8];
% %Lines leaving via each edge of the box.
% p1(:,3)=[0.6;-0.1]; p2(:,3)=[1.6;-0.2];
% p1(:,4)=[-0.6;0.1]; p2(:,4)=[-1.6;0.2];
% p1(:,5)=[-0.1;0.5]; p2(:,5)=[-0.2;2.8];
% p1(:,6)=[0.1;-0.5]; p2(:,6)=[0.2;-1.8];
%
% figure(1)
% clf
% hold on
% %Draw the bounding box lines, extending them.
% extDist=3;
% plot([-extDist,extDist],[-1,-1],'-k','linewidth',2)
% plot([-extDist,extDist],[1,1],'-k','linewidth',2)
% plot([-1,-1],[-extDist,extDist],'-k','linewidth',2)
% plot([1,1],[-extDist,extDist],'-k','linewidth',2)
% for k=1:numLines
%     %Plot the line segments
%     plot([p1(1,k),p2(1,k)],[p1(2,k),p2(2,k)],'-b','linewidth',4)
%     %Clip the segments
%     [p1Clip,p2Clip]=clipLineSegment2Rect(rect,p1(:,k),p2(:,k));
%     %Plot the clipped segments.
%     if(~isempty(p1Clip))
%         plot([p1Clip(1),p2Clip(1)],[p1Clip(2),p2Clip(2)],'-r','linewidth',2)
%     end
% end
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis",
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.
%[2] F. Dévai, "Analysis of the Nichioll-Lee-Nicholl Algorithm," in
%    Proceedings of the International Conference on Computational Science
%    and its Applications, Signapore, pp. 726-736, 9-12 May, 2005.
%
%August 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xMin=rectMaxMin(1);
xMax=rectMaxMin(2);
yMin=rectMaxMin(3);
yMax=rectMaxMin(4);

x1=p1(1);
y1=p1(2);
x2=p2(1);
y2=p2(2);

if(x1<xMin)
    [x1,y1,x2,y2,isVisible]=leftColumn(xMin,yMax,xMax,yMin,x1,y1,x2,y2);
elseif(x1>xMax)
    [x1,y1]=rotate180c(x1,y1);
    [x2,y2]=rotate180c(x2,y2);

    [x1,y1,x2,y2,isVisible]=leftColumn(-xMax,-yMin,-xMin,-yMax,x1,y1,x2,y2);

    [x1,y1]=rotate180c(x1,y1);
    [x2,y2]=rotate180c(x2,y2);
else
    [x1,y1,x2,y2,isVisible]=centerColumn(xMin,yMax,xMax,yMin,x1,y1,x2,y2);
end

if(isVisible==false)
    %If the line is not in the box at all.
    p1Clip=[];
    p2Clip=[];
else
    p1Clip=[x1;y1];
    p2Clip=[x2;y2];
end
end

function [x1,y1,x2,y2,isVisible]=leftColumn(xL,yT,xR,yB,x1,y1,x2,y2)
%%LEFTCOLUMN The case where (x1,y1) is to the left of the vertical line
%            x=xL.

if(x2<xL)
    isVisible=false;
elseif(y1>yT)
    [x1,y1,x2,y2,isVisible]=topLeftCorner(xL,yT,xR,yB,x1,y1,x2,y2);
elseif(y1<yB)
    y1=reflectxAxis(y1);
    y2=reflectxAxis(y2);
    [x1,y1,x2,y2,isVisible]=topLeftCorner(xL,-yB,xR,-yT,x1,y1,x2,y2);

    y1=reflectxAxis(y1);
    y2=reflectxAxis(y2);
else
    [x1,y1,x2,y2,isVisible]=leftEdge(xL,yT,xR,yB,x1,y1,x2,y2);
end
end

function [x1,y1,x2,y2,isVisible]=topLeftCorner(xL,yT,xR,yB,x1,y1,x2,y2)
%TOPLEFTCORNER The case where (x1,y1) is in the top-left corner.

    if(y2>yT)
        isVisible=false;
    else
        relx2=x2-x1;
        rely2=y2-y1;

        topProd=(yT-y1)*relx2;
        leftProd=(xL-x1)*rely2;
        %If (x2,y2) is below the line connecting (x1,y1) and (xL,yT).
        if(topProd>leftProd)
            [x1,y1,x2,y2,isVisible]=leftBottomRegion(xL,xR,yB,x1,y1,x2,y2,relx2,rely2,leftProd);
        else
            [x1,y1]=reflectxMinusy(x1,y1);
            [x2,y2]=reflectxMinusy(x2,y2);

            [x1,y1,x2,y2,isVisible]=leftBottomRegion(-yT,-yB,-xR,x1,y1,x2,y2,-rely2,-relx2,topProd);

            [x1,y1]=reflectxMinusy(x1,y1);
            [x2,y2]=reflectxMinusy(x2,y2);
        end
    end
end

function [x1,y1,x2,y2,isVisible]=leftBottomRegion(xL,xR,yB,x1,y1,x2,y2,relx2,rely2,leftProd)
%%LEFTBOTTOMREGION The case where (x1,y1) is in the left bottom corner.

if(y2>=yB)
    if(x2>xR)
        %The right-edge intersection.
        y2=y1+(xR-x1)*rely2/relx2;
        x2=xR;
    end

    %The left-edge intersection.
    y1=y1+leftProd/relx2;
    x1=xL;
    isVisible=true;
else
    bottomProd=(yB-y1)*relx2;

    %If (x2,y2) is below the (x1,y1) to (xR,yB) line.
    if(bottomProd>leftProd)
        isVisible=false;
    else
        if(x2>xR)
            rightProd=(xR-x1)*rely2;

            %If (x2,y2) is below the (x1,y1) to (xR,yB) line .
            if(bottomProd>rightProd)
                %The bottom-edge intersection.
                x2=x1+bottomProd/rely2;
                y2=yB;
            else
                %The right-edge intersection.
                y2=y1+rightProd/relx2;
                x2=xR;
            end
        else
            %The bottom-edge intersection.
            x2=x1+bottomProd/rely2;
            y2=yB;
        end

        %The left-edge intersection.
        y1=y1+leftProd/relx2;
        x1=xL;
        isVisible=true;
    end
end
end

function [x1,y1,x2,y2,isVisible]=leftEdge(xL,yT,xR,yB,x1,y1,x2,y2)
%%LEFTEDGE (x1,y1) is in the left edge.

    if(x2<xL)
        isVisible=false;
    elseif(y2<yB)
        [x1,y1,x2,y2,isVisible]=p2Bottom(xL,xR,yB,x1,y1,x2,y2);
    elseif(y2>yT)
        y1=reflectxAxis(y1);
        y2=reflectxAxis(y2);

        [x1,y1,x2,y2,isVisible]=p2Bottom(xL,xR,-yT,x1,y1,x2,y2);

        y1=reflectxAxis(y1);
        y2=reflectxAxis(y2);
    else
        relx2=x2-x1;
        rely2=y2-y1;

        if(x2>xR)%The right-edge intersection.
            y2=y1+(xR-x1)*rely2/relx2;
            x2=xR;
        end

        %The left-edge intersection.
        y1=y1+(xL-x1)*rely2/relx2;
        x1=xL;
        isVisible=true;
    end
end

function [x1,y1,x2,y2,isVisible]=centerColumn(xL,yT,xR,yB,x1,y1,x2,y2)
%CENTERCOLUMN The case where (x1,y1) is between the left and right
%             boundaries.

    if(y1>yT)
        [x1,y1]=rotate270c(x1,y1);
        [x2,y2]=rotate270c(x2,y2);
        [x1,y1,x2,y2,isVisible]=leftEdge(-yT,xR,-yB,xL,x1,y1,x2,y2);
        [x1,y1]=rotate90c(x1,y1);
        [x2,y2]=rotate90c(x2,y2);
    elseif(y1<yB)
        [x1,y1]=rotate90c(x1,y1);
        [x2,y2]=rotate90c(x2,y2);
        [x1,y1,x2,y2,isVisible]=leftEdge(yB,-xL,yT,-xR,x1,y1,x2,y2);
        [x1,y1]=rotate270c(x1,y1);
        [x2,y2]=rotate270c(x2,y2);
    else
        [x1,y1,x2,y2,isVisible]=inside(xL,yT,xR,yB,x1,y1,x2,y2);
    end
end

function [x1,y1,x2,y2,isVisible]=inside(xL,yT,xR,yB,x1,y1,x2,y2)
%%INSIDE The case where (x1,y1) is inside the rectangle.

if(x2<xL)
    [x1,y1,x2,y2]=p2Left(xL,yT,yB,x1,y1,x2,y2);
elseif(x2>xR)
	[x1,y1]=rotate180c(x1,y1);
	[x2,y2]=rotate180c(x2,y2);
    [x1,y1,x2,y2]=p2Left(-xR,-yB,-yT,x1,y1,x2,y2);
    [x1,y1]=rotate180c(x1,y1);
	[x2,y2]=rotate180c(x2,y2);
elseif(y2>yT)
    x2=x1+(x2-x1)*(yT-y1)/(y2-y1);
    y2=yT;
elseif(y2<yB)
    x2=x1+(x2-x1)*(yB-y1)/(y2-y1);
    y2=yB;
end

isVisible=true;
end

function [x1,y1,x2,y2]=p2Left(xL,yT,yB,x1,y1,x2,y2)
%%P2LEFT

    if(y2>yT)
        [x1,y1,x2,y2]=p2LeftTop(xL,yT,x1,y1,x2,y2);
    elseif(y2<yB)
        [x1,y1]=rotate90c(x1,y1);
        [x2,y2]=rotate90c(x2,y2);
        [x1,y1,x2,y2]=p2LeftTop(yB,-xL,x1,y1,x2,y2);
        [x1,y1]=rotate270c(x1,y1);
        [x2,y2]=rotate270c(x2,y2);
    else
        y2=y1+(y2-y1)*(xL-x1)/(x2-x1);
        x2=xL;
    end
end

function [x1,y1,x2,y2]=p2LeftTop(xL,yT,x1,y1,x2,y2)
%%P2LEFTTOP

    relx2=x2-x1;
    rely2=y2-y1;
    leftProduct=rely2*(xL-x1);
    topProduct=relx2*(yT-y1);

    if(topProduct>leftProduct)
        x2=x1+topProduct/rely2;
        y2=yT;
    else
        y2=y1+leftProduct/relx2;
        x2=xL;
    end
end

function [x1,y1,x2,y2,isVisible]=p2Bottom(xL,xR,yB,x1,y1,x2,y2)
%%P2BOTTOM (x1,y1) is in the left edge; (x2,y2) is not beyond the left
%          edge,  and (x2,y2) is beyond the bottom edge.

    relx2=x2-x1;
    rely2=y2-y1;
    leftProd=(xL-x1)*rely2;
    bottomProd=(yB-y1)*relx2;
    if(bottomProd>leftProd)
        isVisible=false;
    else
       if(x2<=xR)
           x2=x1+bottomProd/rely2;
           y2=yB;
       else
           rightProd=(xR-x1)*rely2;
           if(bottomProd>rightProd)
               x2=x1+bottomProd/rely2;
               y2=yB;
           else
               y2=y1+rightProd/relx2;
               x2=xR;
           end
       end
       y1=y1+leftProd/relx2;
       x1=xL;
       isVisible=true;
    end
end

function [x,y]=reflectxMinusy(x,y)
%%REFLECTXMINUSY This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis",
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

t=x;
x=-y;
y=-t;

end

function y=reflectxAxis(y)
%%REFLECTAXIS This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis",
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

y=-y;

end

function [x,y]=rotate90c(x,y)
%%ROTATE90C This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis",
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

t=x;
x=y;
y=-t;

end

function [x,y]=rotate180c(x,y)
%%ROTATE180C This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis",
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

x=-x;
y=-y;
end

function [x,y]=rotate270c(x,y)
%%ROTATE270C This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis",
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

t=x;
x=-y;
y=t;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y = icplxdual2D(w, J, Fsf, sf)

% Inverse Dual-Tree Complex 2D Discrete Wavelet Transform
%
% USAGE:
%   y = icplxdual2D(w, J, Fsf, sf)
% INPUT:
%   w - wavelet coefficients
%   J - number of stages
%   Fsf - synthesis filters for final stage
%   sf - synthesis filters for preceeding stages
% OUTPUT:
%   y - output array
% See cplxdual2D
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

for j = 1:J
    for m = 1:3
        [w{j}{1}{1}{m} w{j}{2}{2}{m}] = pm(w{j}{1}{1}{m},w{j}{2}{2}{m});
        [w{j}{1}{2}{m} w{j}{2}{1}{m}] = pm(w{j}{1}{2}{m},w{j}{2}{1}{m});
    end
end

y = zeros(size(w{1}{1}{1}{1})*2);
for m = 1:2
    for n = 1:2
        lo = w{J+1}{m}{n};
        for j = J:-1:2
            lo = sfb2D(lo, w{j}{m}{n}, sf{m}, sf{n});
        end
        lo = sfb2D(lo, w{1}{m}{n}, Fsf{m}, Fsf{n});
        y = y + lo;
    end
end

% normalization
y = y/2;

author: thomaskuestner


text:
%% Delete path
%
% Delete all subdirectories under pathstr to search path.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

rmpath(genpath(pwd));
savepath;

repo: ifem


text:
classdef Builder_DesignVarMonitor_Abstract < handle

    properties (Access = protected)
        monitor
    end

    methods (Access = protected, Abstract)

        createCamera(obj)

    end

    methods (Access = public)

        function build(obj,monitor)
            obj.init(monitor)
            obj.setCamera();
        end

    end

    methods (Access = private)

        function init(obj,monitor)
            obj.monitor = monitor;
        end

        function setCamera(obj)
            cam = obj.createCamera();
            obj.monitor.setCamera(cam);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = ARMAX_close_dialog(varargin)
% ARMAX_CLOSE_DIALOG M-file for ARMAX_close_dialog.fig
%      ARMAX_CLOSE_DIALOG by itself, creates a new ARMAX_CLOSE_DIALOG or raises the
%      existing singleton*.
%
%      H = ARMAX_CLOSE_DIALOG returns the handle to a new ARMAX_CLOSE_DIALOG or the handle to
%      the existing singleton*.
%
%      ARMAX_CLOSE_DIALOG('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ARMAX_CLOSE_DIALOG.M with the given input arguments.
%
%      ARMAX_CLOSE_DIALOG('Property','Value',...) creates a new ARMAX_CLOSE_DIALOG or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ARMAX_close_dialog_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ARMAX_close_dialog_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ARMAX_close_dialog

% Last Modified by GUIDE v2.5 19-Aug-2006 01:03:58

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ARMAX_close_dialog_OpeningFcn, ...
                   'gui_OutputFcn',  @ARMAX_close_dialog_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ARMAX_close_dialog is made visible.
function ARMAX_close_dialog_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ARMAX_close_dialog (see VARARGIN)

% Choose default command line output for ARMAX_close_dialog
handles.output = 'Yes';

% Update handles structure
guidata(hObject, handles);

% Insert custom Title and Text if specified by the user
% Hint: when choosing keywords, be sure they are not easily confused
% with existing figure properties.  See the output of set(figure) for
% a list of figure properties.
if(nargin > 3)
    for index = 1:2:(nargin-3),
        if nargin-3==index, break, end
        switch lower(varargin{index})
         case 'title'
          set(hObject, 'Name', varargin{index+1});
         case 'string'
          set(handles.text1, 'String', varargin{index+1});
        end
    end
end

% Determine the position of the dialog - centered on the callback figure
% if available, else, centered on the screen
FigPos=get(0,'DefaultFigurePosition');
OldUnits = get(hObject, 'Units');
set(hObject, 'Units', 'pixels');
OldPos = get(hObject,'Position');
FigWidth = OldPos(3);
FigHeight = OldPos(4);
if isempty(gcbf)
    ScreenUnits=get(0,'Units');
    set(0,'Units','pixels');
    ScreenSize=get(0,'ScreenSize');
    set(0,'Units',ScreenUnits);

    FigPos(1)=1/2*(ScreenSize(3)-FigWidth);
    FigPos(2)=2/3*(ScreenSize(4)-FigHeight);
else
    GCBFOldUnits = get(gcbf,'Units');
    set(gcbf,'Units','pixels');
    GCBFPos = get(gcbf,'Position');
    set(gcbf,'Units',GCBFOldUnits);
    FigPos(1:2) = [(GCBFPos(1) + GCBFPos(3) / 2) - FigWidth / 2, ...
                   (GCBFPos(2) + GCBFPos(4) / 2) - FigHeight / 2];
end
FigPos(3:4)=[FigWidth FigHeight];
set(hObject, 'Position', FigPos);
set(hObject, 'Units', OldUnits);

% Show a question icon from dialogicons.mat - variables questIconData
% and questIconMap
load dialogicons.mat

IconData=questIconData;
questIconMap(256,:) = get(handles.figure1, 'Color');
IconCMap=questIconMap;

Img=image(IconData, 'Parent', handles.axes1);
set(handles.figure1, 'Colormap', IconCMap);

set(handles.axes1, ...
    'Visible', 'off', ...
    'YDir'   , 'reverse'       , ...
    'XLim'   , get(Img,'XData'), ...
    'YLim'   , get(Img,'YData')  ...
    );

% Make the GUI modal
set(handles.figure1,'WindowStyle','modal')

% UIWAIT makes ARMAX_close_dialog wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = ARMAX_close_dialog_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.figure1);

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.output = get(hObject,'String');

% Update handles structure
guidata(hObject, handles);

% Use UIRESUME instead of delete because the OutputFcn needs
% to get the updated handles structure.
uiresume(handles.figure1);

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.output = get(hObject,'String');

% Update handles structure
guidata(hObject, handles);

% Use UIRESUME instead of delete because the OutputFcn needs
% to get the updated handles structure.
uiresume(handles.figure1);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.output = 'No';
guidata(hObject, handles);
if isequal(get(handles.figure1, 'waitstatus'), 'waiting')
    % The GUI is still in UIWAIT, us UIRESUME
    uiresume(handles.figure1);
else
    % The GUI is no longer waiting, just close it
    delete(handles.figure1);
end

% --- Executes on key press over figure1 with no controls selected.
function figure1_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Check for "enter" or "escape"
if isequal(get(hObject,'CurrentKey'),'escape')
    % User said no by hitting escape
    handles.output = 'No'

    % Update handles structure
    guidata(hObject, handles);

    uiresume(handles.figure1);
end

if isequal(get(hObject,'CurrentKey'),'return')
    uiresume(handles.figure1);
end

author: bashtage

repo: mfe-toolbox


text:
function [W_G, IDX2, S, IDX_vote, S_vote] = RBF_affnity_knn_fast_outliers(data_noise,params,f_outliers)
%mydisp('Computing RBF affinity')
if nargin < 2
    neibours = 20;
else
    neibours = params.n_neighbors;
end
Sigma = params.Sigma;
N=size(data_noise,2);

%W_G=sparse(zeros(N,N));
W_G=spalloc(N,N,300*N);
%W_G=(zeros(N,N));
affinity_type = params.affinity_type;
retain = params.m*params.n_neighbors;

%S = pdist2imp(data_noise',data_noise','sqeuclidean');
switch  affinity_type
    case 'cosine'
    %S = pdist2imp(data_noise',data_noise','cosine');
    [IDX_vote,S]=knnsearch(data_noise',data_noise','k',retain+1,'distance','cosine');
    IDX_vote=IDX_vote(:,2:end);
    S=S(:,2:end);
    S_vote=S;
    IDX_vote=IDX_vote';
    S=S';
    S_vote=S;
    S = S(1:neibours,:);
    otherwise
%   S = single(pdist2imp(data_noise',data_noise','sqeuclidean'));
%   S = (pdist2imp(data_noise',data_noise','sqeuclidean'));
   %[IDX_vote,S]=knnsearch(data_noise',data_noise','k',neibours+1,'distance','euclidean');
    [IDX_vote,S]=knnsearch(data_noise',data_noise','k',retain+1,'distance','euclidean');

     IDX_vote=IDX_vote(:,2:end);
     S=S(:,2:end);
     IDX_vote=IDX_vote';
     S=S';

end
%% Clipping NN

IDX2 = IDX_vote(1:neibours,:);
S_vote=S;
S = S(1:neibours,:);
a1=1:N;
a2=repmat(a1,neibours);
a2=a2(1:neibours,1:N);
a3=a2(:);
I=a3;
J=IDX2(:);
% s = sqrt( sum( ( data_noise(:,I)-data_noise(:,J) ).^2, 1) );
% s=S(:);
% s = sqrt( sum( ( data_noise(:,I)-data_noise(:,J) ).^2, 1) );

Sigma=0.1;
% s = (1.1./ (1.1+ (f_outliers(I)-f_outliers(J)).^2)).^2;
% s1= ( sum( ( data_noise(:,I)-data_noise(:,J) ).^2, 1) )
% s2= (1.0./ (1+ s1)).^2;
% s2=s2';

s = (f_outliers(I)-f_outliers(J) ).^2;
w_temp=exp(-(s) ./ Sigma );
%s=s.*s2;
%w_temp=s;
waff = sparse(I,J,w_temp,N,N);
W_k=waff;
%W_k= (W_k + W_k')./2 ;
W_k= ((W_k + W_k')+ abs(W_k - W_k'))./2;
W_G=sparse(W_k);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function D = sampling_dep(f, depidx, domain, gridsize)
%SAMPLING_DEP Sampling a (multivariate) function
%	D = SAMPLING_DEP(f, depidx, domain, gridsize)
%
%	f        - function of a real^P vector (f depends on a subset of its input)
%	depidx   - index of the dependent elements
%	domain   - [min1 max1;... minP maxP] intervals for each element
%	gridsize - number of sampling grid points for each dependent element
%
%	D        - P-dimensional array of the sampled data
%
%	This function is used internally by the tptoool toolbox
%
%	eg.:   sampling_dep(@(x) x(2)+x(3), [2 3] [-1 1; 0 3], [7 5])
%
%	See also SAMPLING_LPV

% internal function (used by sampling_lpv)

%DIMENSIONS
P = length(gridsize);
p = zeros(max(depidx), 1);

%SAMPLING
if P > 0
	% allocation of sample array
	siz = prod(gridsize);
	D = zeros(siz,1);
	a = domain(:,1);
	b = domain(:,2);
	step = (b - a) ./ (gridsize - 1);

	% sampling
	z = ones(P,1);
	z(1) = 0;
	for k = 1:siz
		% next grid point index
		z = nexti(gridsize, z);
		% argument
		p(depidx) = a + step .* (z - 1);
		% sampling
		D(k) = f(p);
	end

	if P > 1
		% reshape to proper size
		D = reshape(D, gridsize');
	end
else
	D = f(p);
end

%SUBFUNCTION: Next grid point index
function n = nexti(gridsize, n)
i = 1;
n(i) = n(i) + 1;
while n(i) > gridsize(i)
	n(i) = 1;
	i = i + 1;
	n(i) = n(i) + 1;
end

author: Sable

repo: mcbench-benchmarks


text:
clc
close all
clear
% This file tests feasibility between MPC-DCBF and NMPC-DCBF.

%% System setup

timestep = 0.1;
system_param.A = [[1 ,timestep, 0]; [0, 1, timestep]; [0, 0, 1]];
system_param.B = [0; 0; timestep];
system_param.ul = -1;
system_param.uu = 1;
system_param.timestep = timestep;
x0 = [0;0;0]; % this value will be overrided
t0 = 0.0;

gammalist = [0.05, 0.1, 0.15, 0.2];
for gammaindex = 1:length(gammalist)
    %% Problem setup with different hyperparameters

    % MPC-DCBF simulator
    simulator_mpcdcbf = CBFDT(system_param, x0, t0);
    param_mpcdcbf = ParamMPCDCBF(8, gammalist(gammaindex), 10.0*eye(3), 10.0*eye(3), 1.0);
    simulator_mpcdcbf.setOpt('mpcdcbf', param_mpcdcbf);
    % NMPC-DCBF simulator
    simulator_nmpcdcbf = CBFDT(system_param, x0, t0);
    param_nmpcdcbf = ParamNMPCDCBF(8, 8, gammalist(gammaindex), 10.0*eye(3), 10.0*eye(3), 1.0, 10.0);
    simulator_nmpcdcbf.setOpt('nmpcdcbf', param_nmpcdcbf);
    % iterate over states
    x1list = linspace(-2,0,11);
    x2list = linspace(0,2,11);
    x3list = linspace(0,2,11);
    % data collection
    feas_points_mpcdcbf = [];
    feas_points_nmpcdcbf = [];
    infeas_points_all = [];

    %% Test feasibility with sampling states among controllers

    for k1 = 1:length(x1list)
        for k2 = 1:length(x2list)
            for k3 = 1:length(x3list)
                xfeas = [x1list(k1); x2list(k2); x3list(k3)];
                % solve the problem with MPC-DCBF
                simulator_mpcdcbf.xcurr = xfeas;
                [feas_mpcdcbf, ~, ~, ~] = simulator_mpcdcbf.solve;
                % solve the problem with NMPC-DCBF
                simulator_nmpcdcbf.xcurr = xfeas;
                [feas_nmpcdcbf, ~, ~, ~] = simulator_nmpcdcbf.solve;
                if feas_mpcdcbf == 1
                    feas_points_mpcdcbf = [feas_points_mpcdcbf, xfeas];
                end
                if feas_nmpcdcbf == 1
                    feas_points_nmpcdcbf = [feas_points_nmpcdcbf, xfeas];
                end
                if feas_mpcdcbf == 0 && feas_nmpcdcbf == 0
                    infeas_points_all = [infeas_points_all, xfeas];
                end
            end
        end
    end

    %% Plotting
    close all
    outer_color = [0.3010, 0.7450, 0.9330];
    inner_color = [1, 0, 0];
    figure('Renderer', 'painters', 'Position', [0 0 400 400]);
    set(gca,'LooseInset',get(gca,'TightInset'));
    hold on;
    scatter3(feas_points_mpcdcbf(1,:),feas_points_mpcdcbf(2,:),feas_points_mpcdcbf(3,:), 5, 'o','LineWidth',1.0, 'MarkerEdgeColor',inner_color,'MarkerFaceColor',inner_color);
    scatter3(feas_points_nmpcdcbf(1,:),feas_points_nmpcdcbf(2,:),feas_points_nmpcdcbf(3,:), 20, 'o','LineWidth',1.2, 'MarkerEdgeColor',outer_color);
    axis equal
    if gammaindex == 1
        h=get(gca,'Children');
        h_legend = legend(h([end, end-1]),...
            {'MPC-DCBF', 'NMPC-DCBF'}, 'Location', 'NorthEast');
        set(h_legend, 'Interpreter','latex');
    end
    set(gca,'LineWidth', 1.0, 'FontSize', 15);
    xlabel('$x (m)$','interpreter','latex','FontSize',20);
    ylabel('$v (m/s)$','interpreter','latex','FontSize',20);
    zlabel('$a (m/s^2)$','interpreter','latex','FontSize',20);
    xlim([-2, 0]);
    ylim([0, 2]);
    zlim([0, 2]);
    view(70, 4);
    grid on;
    figurename_eps = "feasibility-mpcdcbf" + gammaindex + ".eps";
    figurename_png = "feasibility-mpcdcbf" + gammaindex + ".png";
    dataname = "feasibility-mpcdcbf" + gammaindex + ".mat";
    % save data and generate figures
    print(gcf,strcat('figures/',figurename_eps), '-depsc');
    print(gcf,strcat('figures/',figurename_png), '-dpng', '-r800');
    save(strcat('data/',dataname));
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
A_bounds = [-eye(3);
            1,1,1];
b_bounds = [0;0;0;1];
lb = [0;0;0];
ub = [1;1;1];

[C, d] = iris.inner_ellipsoid.mosek_nofusion(A_bounds, b_bounds);
assert(det(C) > 0.01);

iris.drawing.draw_3d(A_bounds, b_bounds, C, d, [], lb, ub);

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function [ Z, H, dnorm ] = deep_seminmf ( X, layers, varargin )

% Process optional arguments
pnames = { ...
    'z0' 'h0' 'bUpdateH' 'bUpdateLastH' 'maxiter' 'TolFun', ...
    'verbose', 'bUpdateZ', 'cache', 'gnd' ...
};

X = bsxfun(@rdivide,X,sqrt(sum(X.^2,1)));

num_of_layers = numel(layers);

Z = cell(1, num_of_layers);
H = cell(1, num_of_layers);

dflts  = {0, 0, 1, 1, 500, 1e-5, 1, 1, 1, 0};

[z0, h0, bUpdateH, bUpdateLastH, maxiter, tolfun, verbose, bUpdateZ, cache, gnd] = ...
        internal.stats.parseArgs(pnames,dflts,varargin{:});

if  ~iscell(h0)
    for i_layer = 1:length(layers)
        if i_layer == 1
            % For the first layer we go linear from X to Z*H, so we use id
            V = X;
        else
            V = H{i_layer-1};
        end

        if verbose
            display(sprintf('Initialising Layer #%d with k=%d with size(V)=%s...', i_layer, layers(i_layer), mat2str(size(V))));
        end
        if ~iscell(z0)
            % For the later layers we use nonlinearities as we go from
            % g(H_{k-1}) to Z*H_k
            [Z{i_layer}, H{i_layer}, ~] = ...
                 seminmf(V, ...
                     layers(i_layer), ...
                     'maxiter', maxiter, ...
                     'bUpdateH', true, 'bUpdateZ', bUpdateZ, 'verbose', verbose, 'save', cache, 'fast', 1);
        else
            display('Using existing Z');
            [Z{i_layer}, H{i_layer}, ~] = ...
                 seminmf(V, ...
                     layers(i_layer), ...
                     'maxiter', 1, ...
                     'bUpdateH', true, 'bUpdateZ', 0, 'z0', z0{i_layer}, 'verbose', verbose, 'save', cache, 'fast', 1);
        end
    end

else
    Z=z0;
    H=h0;

    if verbose
        display('Skipping initialization, using provided init matrices...');
    end
end

dnorm0 = cost_function(X, Z, H);
dnorm = dnorm0 + 1;

if verbose
    display(sprintf('#%d error: %f', 0, dnorm0));
end

%% Error Propagation
if verbose
    display('Finetuning...');
end
H_err = cell(1, num_of_layers);

for iter = 1:maxiter
    H_err{numel(layers)} = H{numel(layers)};
    for i_layer = numel(layers)-1:-1:1
        H_err{i_layer} = Z{i_layer+1} * H_err{i_layer+1};
    end

    for i = 1:numel(layers)
        if bUpdateZ
            try
                if i == 1
                    Z{i} = X  * pinv(H_err{1});
                else
                    Z{i} = pinv(D') * X * pinv(H_err{i});
                end
            catch
                display(sprintf('Convergance error %f. min Z{i}: %f. max %f', norm(Z{i}, 'fro'), min(min(Z{i})), max(max(Z{i}))));
            end
        end

        if i == 1
            D = Z{1}';
        else
            D = Z{i}' * D;
        end

        if bUpdateH && (i < numel(layers) || (i == numel(layers) && bUpdateLastH))
            A = D * X;
            Ap = (abs(A)+A)./2;
            An = (abs(A)-A)./2;

            B = D * D';

            Bp = (abs(B)+B)./2;
            Bn = (abs(B)-B)./2;

            H{i} = H{i} .* sqrt((Ap + Bn * H{i}) ./ max(An + Bp * H{i}, 1e-10));
        end
    end

    assert(i == numel(layers));

    dnorm = cost_function(X, Z, H);

    if verbose
        display(sprintf('#%d error: %f', iter, dnorm));
    end

%     assert(dnorm <= dnorm0 + 0.01, ...
%         sprintf('Rec. error increasing! From %f to %f. (%d)', ...
%         dnorm0, dnorm, iter) ...
%     );

    if verbose && length(gnd) > 1
       if mod(iter, 50) == 0
          ac = evalResults(H{numel(H)}, gnd);
          fprintf(1, 'Clustering accuracy is %.2f\n', ac);
       end
    end
%
%     if dnorm0-dnorm <= tolfun*max(1,dnorm0)
%         if verbose
%             display( ...
%                 sprintf('Stopped at %d: dnorm: %f, dnorm0: %f', ...
%                     iter, dnorm, dnorm0 ...
%                 ) ...
%             );
%         end
%         break;
%     end

    dnorm0 = dnorm;
end
end

function error = cost_function(X, Z, H)
    error = norm(X - reconstruction(Z, H), 'fro');
end

function [ out ] = reconstruction( Z, H )

    out = H{numel(H)};

    for k = numel(H) : -1 : 1;
        out =  Z{k} * out;
    end

end

author: jwyang


text:
function [vertices, faces, facevertexcdata, renderer] = u3d_pre_surface(ax)
%U3D_PRE_SURFACE    Preprocess surface output to u3d.
%    U3D_PRE generates the input for the MESH_TO_LATEX function from
%    Alexandre Gramfort from your surface-graphs. The surface graphs 3d-model can be
%    displayed in u3d - format in pdf or if you don't delete the u3d-files
%    which are generated by MESH_TO_LATEX you can use deepview from
%    righthemisphere(http://www.righthemisphere.com/products/client-products/deep-view)
%    to embed the modell in Microsoft-Office products (Word, Excel, PowerPoint).
%
% usage
%   [vertices, faces, facevertexcdata] = U3D_PRE_SURFACE
%   [vertices, faces, facevertexcdata] = U3D_PRE_SURFACE(h)
%
% optional input
%   ax = axes object handle
%
% output
%   vertices = row vectors of point positions, as row cell array
%              for multiple surfaces
%            = {1 x #surfaces}
%            = {[#vertices x 3], ... }
%   faces = for each surface triangle face, indices of its 3 vertices,
%           these indices refer to the columns of matrix vertices,
%           as row cell array for multiple surfaces
%         = {1 x #surfaces}
%         = {[#faces x 3], ... }
%   facevertexcdata = RGB color information at each vertex,
%                     as row cell array for multiple surfaces
%                   = {1 x #surfaces}
%                   = {[#vertices x 3], ... }
%
% See also fig2idtf, u3d_pre_line, u3d_pre_patch, u3d_pre_quivergroup,
%          u3d_pre_contourgroup.
%
% File:      u3d_pre_surface.m
% Original Author: Sven Koerner, koerner(underline)sven(add)gmx.de
% Author:    Ioannis Filippidis, jfilippidis@gmail.com (added support for multiple surfaces)
% Date:      2012.06.10 -
% Language:  MATLAB R2012a
% Purpose:   preprocess surface children of axes for u3d export
% Copyright:
%
% License to use and modify this code is granted freely to all interested,
% as long as the original author is referenced and attributed as such.
% The original author maintains the right to be solely associated with this work.

%% input
if nargin < 1
    sh = findobj('type', 'surface');
else
    objs = get(ax, 'Children');
    sh = findobj(objs, 'type', 'surface');
end

if isempty(sh)
    disp('No surfaces found.');
    vertices            = [];
    faces               = [];
    facevertexcdata     = [];
    renderer = [];
    return
end

%% process each surface
N = size(sh, 1); % number of surfaces
vertices = cell(1, N);
faces = cell(1, N);
facevertexcdata = cell(1, N);
renderer = cell(1, N);
for i=1:N
    disp(['     Preprocessing surface No.', num2str(i) ] );
    h = sh(i, 1);

    [v, f, fvx, r] = single_surf_preprocessor(h);

    vertices{1, i} = v;
    faces{1, i} = f;
    facevertexcdata{1, i} = fvx;
    renderer{1, i} = r;
end

function [vertices, faces, facevertexcdata, renderer] = single_surf_preprocessor(h)
%% shading -> renderer in adobe reader
edgecolor = get(h, 'EdgeColor');
if strcmp(edgecolor, 'none')
    renderer = 'Solid';
else
    renderer = 'SolidWireframe';
end

%% get defined data-points
X = get(h, 'XData');
Y = get(h, 'YData');
Z = get(h, 'ZData');

realcolor = get_surf_color(h);

[faces, vertices, facevertexcdata] = surf2patch(X, Y, Z, realcolor, 'triangles');

%% remove nan vertices and faces
% remove faces using at least one vertex with some nan coordinate
nan_vertices = any(isnan(vertices), 2);
nan_faces = nan_vertices(faces);
nan_faces = any(nan_faces, 2);
nan_faces = ~nan_faces;
faces = faces(nan_faces, :);

% vertices with nan are not used anymore
% just make them contain numbers
% DO NOT REMOVE them! This would destroy face indexing
vertices(isnan(vertices) ) = 0;
vertices = vertices.';

%% surface concatenation (obsolete - although it reduces file size)
%{
[n, m] = size(X, 1);
Vi = n *m; % number of vertices
Fi = 2 *(n-1) *(m -1); % number of faces, due to wraping (closing of the surface)
disp(['     Object Vertex # = ', num2str(Vi) ] )
disp(['     Object Face # = ', num2str(Fi) ] )

tempfvc.faces = tempfvc.faces +V; % shift to account for previous vertices

% append to previous
fvc.vertices = [fvc.vertices; tempfvc.vertices];
fvc.faces = [fvc.faces; tempfvc.faces];
fvc.facevertexcdata = [fvc.facevertexcdata; tempfvc.facevertexcdata];

F = size(faces, 1);
[C1, C2] = size(facevertexcdata);

V = size(fvc.vertices, 1); %V = V +Vi;

disp(['     Vertex # = ', num2str(V) ] )
disp(['     Face # = ', num2str(F) ] )
disp(['     Vertex Colors = ', num2str(C1), ' x ', num2str(C2) ] )
%}

function [realcolor] = get_surf_color(h)
% texture mapping ?
facecolor = get(h, 'FaceColor');
if strcmp(facecolor, 'texturemap')
    realcolor = [];
    return
end

% texturemapping is when CData is an image and has different size than
% XData and the number of faces, so it is not treated here yet.

%% CData
cdata = get(h, 'CData');
k = size(cdata, 3);

% true color ?
if k == 3
    disp('CData is already True Color.')
    cdata(isnan(cdata) ) = 1;
    realcolor = cdata;
    return
end

% not indexed colors ?
if k ~= 1
    error('u3dsurf:cdata', 'size(CData, 3) \notin {1, 3}')
end

%% indexed color to RGB true color
ax = get(h, 'Parent');
realcolor = scaled_ind2rgb(cdata, ax);

function [realcolor] = scaled_ind2rgb(cdata, ax)
[n, m] = size(cdata);
cdata = double(cdata);

cmap = colormap(ax);
nColors = size(cmap, 1);
[cmin, cmax] = caxis(ax);

idx = (cdata -cmin) / (cmax -cmin) *nColors;
idx = ceil(idx);
idx(idx < 1) = 1;
idx(idx > nColors) = nColors;

%% handle nans in idx
nanmask = isnan(idx);
idx(nanmask) = 1; %temporarily replace w/ a valid colormap index

%% realcolor and output
realcolor = zeros(n, m);
for i = 1:3
    c = cmap(idx, i);
    c = reshape(c, n, m);
    realcolor(:, :, i) = c;
end


text:
function result=halfspacedepth(u,v,x,y)

%HALFSPACEDEPTH computes the halfspace depth of a (two-dimensional) point theta
% relative to a bivariate data set.
%
% The algoritm is described in:
%    Rousseeuw, P., Ruts, I. (1996),
%    "AS 307: Bivariate Location Depth",
%    Applied Statistics (JRSS-C), 45, 516-526.
%
% Required input arguments:
%            u : first coordinate of the point theta
%            v : second coordinate of the point theta
%            x : vector containing the first coordinates of all the data
%                points
%            y : vector containing the second coordinates of all the data
%                points
%
% I/O: [result] = halfspacedepth(u,v,x,y);
%
% This function is part of the Matlab Library for Robust Analysis,
% available at:
%           http://wis.kuleuven.be/stat/robust.html
%
% Written by Fabienne Verwerft
%
%
% Checking input
%
if not(length(x)==length(y))
    error('The vectors x and y must have the same length.')
end
if sum(isnan(x))>=1 || sum(isnan(y))>=1
    error('Missing values are not allowed')
end
%
% m is the number of data points that coincide with theta
%
eps=0.000001;
n=length(x);
norm=sqrt((x-u).^2 + (y-v).^2);
m=sum(norm<=eps);
ll=find(norm>eps);
norm=norm(ll);
xn=(x(ll)-u)./norm;
yn=(y(ll)-v)./norm;
%
% The vector containes the indices of the elements of x and y that satisfy
% the conditions.
%
k=find((abs(x(ll))>abs(y(ll))) & (xn>=0));
alfa(k)=asin(yn(k));
t=find(alfa(k)<0);
alfa(k(t))=2*pi+alfa(k(t));
k=find((abs(x(ll))>abs(y(ll))) & (xn<0));
alfa(k)=pi-asin(yn(k));
k=find((abs(x(ll))<=abs(y(ll))) & (yn>=0));
alfa(k)=acos(xn(k));
k=find((abs(x(ll))<=abs(y(ll))) & (yn<0));
alfa(k)=2*pi-acos(xn(k));
g=find(alfa>=(2*pi-eps));
alfa(g)=0;
%
nn=n-m;
if nn <= 1
    hdepth=m;
    result=hdepth;
    return
    % all data points coincide with theta
end
%
alfa=sort(alfa);
%
hoek=max(alfa(1)-alfa(nn)+2*pi,max(diff(alfa)));
if hoek > pi+eps
    hdepth=m;
    result=hdepth;
    return
    % hdepth=0 because theta lies outside the datacloud
end
%
% rotation around theta
% nu is the number of angles in the upper halfcircle
%
alfa=alfa-alfa(1);
nu=sum(alfa < (pi-eps));

if nu >= nn
    hdepth=m;
    result=hdepth;
    return
    % hdepth=0 every angle in the upper halfcircle so theta outside the
    % datacloud.
end
%
% construction of the array F
%
beta=alfa+pi;
alfatwee=alfa+2*pi;
A=[alfa,alfatwee,beta];
Aindex(1:2*nn)=1;
Aindex((2*nn+1):3*nn)=0;
[As,Asin]=sort(A);
Aindexs=Aindex(Asin);
pp=cumsum(Aindexs);
juisten=find(Aindexs==0);
F=pp(juisten);
%
% Adjust the array F for the angles that coincide with beta.
%
gelijkab=intersect(find(diff(As)<=eps)+1,juisten);
betagelijka=Asin(gelijkab);
if length(gelijkab)>0
    for i=1:length(gelijkab)
        aantal=sum((As(Aindexs==1)+eps)<As(gelijkab(i)));
        F(betagelijka(i)-2*nn)=aantal;
    end
end
%
gindex=find(diff(alfa)==0)+1;
G=0:(nn-1);
%
% Adjust the array G for angles alfa who coincide
%
if length(gindex)>0
    for i=1:length(gindex)
        aantal=sum(alfa(1:gindex(i))<alfa(gindex(i)));
        G(gindex(i))=aantal;
    end
end
%
k=F-G;
numh=min(min(k,(nn-k)));
result=(numh+m);

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function level = thresh_tool(im,cmap)
%THRESH_TOOL - GUI tool for selecting threshold level
%
%  Example 1 - nonblocking behavior:
%    x=imread('rice.png');
%    thresh_tool(x) %no return value, so MATLAB keeps running
%
%  Example 2 - blocking behavior:
%    x=imread('rice.png');
%    lev=thresh_tool(x) %MATLAB waits for GUI tool to finish first

% Copyright 2004-2010 The MathWorks, Inc.

max_colors=1000;    %practical limit

%calculate bins centers
color_range = double(limits(im));
if isinteger(im)
  %try direct indices first
  num_colors=diff(color_range)+1;
  if num_colors<max_colors %okay
    di=1;                               %inherent bins
  else %too many levels
    num_colors=max_colors;                  %practical limit
    di=diff(color_range)/(num_colors-1);
  end
else %noninteger
  %try infering discrete resolution first (intensities often quantized)
  di=min(diff(sort(unique(im(:)))));
  num_colors=diff(color_range)/di+1;
  if num_colors>max_colors %too many levels
    num_colors=max_colors;                  %practical limit
    di=diff(color_range)/(num_colors-1);
  end
end
bin_ctrs = [color_range(1):di:color_range(2)];
FmtSpec=['%.' num2str(ceil(-log10(di))) 'f'];

%new figure - interactive GUI tool for level segmenting
h_fig=figure;
if nargin>1 & isstr(cmap) & strmatch(lower(cmap),'gray')
  full_map=gray(num_colors);
else
  full_map=jet(num_colors);
end
setappdata(h_fig,'im',im)
setappdata(h_fig,'FmtSpec',FmtSpec)

%top left - input image
h_ax1=axes('unit','norm','pos',[0.05 0.35 0.4 0.60]);
subimage(im,full_map)
axis off, title('Input Image')

%top right - segmented (eventually)
h_ax2=axes('unit','norm','pos',[0.55 0.35 0.4 0.60]);
axis off
setappdata(h_fig,'h_ax2',h_ax2)

%next to bottom - intensity distribution
h_hist=axes('unit','norm','pos',[0.05 0.1 0.9 0.2]);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bin_ctrs=[0:double(max(im(:)))];
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=hist(double(im(:)),bin_ctrs);
bar(bin_ctrs,n)
axis([color_range limits(n(2:end))])
set(h_hist,'xtick',[],'ytick',[])
title('Intensity Distribution')

%very bottom - colorbar
h_cbar=axes('unit','norm','pos',[0.05 0.05 0.9 0.05],'tag','thresh_tool_cbar');
subimage(color_range,[0.5 1.5],1:num_colors,full_map)
set(h_cbar,'ytick',[],'xlim',color_range)
axis normal

%colorbar tick locations
set(h_cbar,'xtick',color_range)

%threshold level - initial guess (graythresh)
lo=double(color_range(1));
hi=double(color_range(2));
norm_im=(double(im)-lo)/(hi-lo);
norm_level=graythresh(norm_im);     %GRAYTHRESH assumes DOUBLE range [0,1]
my_level=norm_level*(hi-lo)+lo;

%display level as vertical line
axes(h_hist)
h_lev=vline(my_level,'-');
set(h_lev,'LineWidth',2,'color',0.5*[1 1 1],'UserData',my_level)
setappdata(h_fig,'h_lev',h_lev)

%attach draggable behavior for user to change level
move_vline(h_lev,@update_plot);

axes(h_cbar)
y_lim = get(h_cbar,'ylim');

% PLACE TEXT LOCATION ON COLORBAR (Laurens)
%h_text = text(my_level,mean(y_lim),num2str(round(my_level)));
h_text = text(my_level,mean(y_lim),'dummy','HorizontalAlignment','Center');
if nargin<2
  text_color=0.5*[1 1 1];
else
  text_color='m';
end
set(h_text,'FontWeight','Bold','color',text_color,'Tag','cbar_text')
movex_text(h_text,my_level)
%%%%%%%%%%%%%%%%%%%%%%%%

%segmented image
bw=(im>my_level);
axes(h_ax2)
hold on
subimage(bw), axis off, axis ij
hold off
title('Segmented')

update_plot

%add reset button (resort to initial guess)
h_reset = uicontrol('unit','norm','pos',[0.0 0.95 .1 .05]);
set(h_reset,'string','Reset','callback',@ResetOriginalLevel)

if nargout>0
  h_done=uicontrol('unit','norm','pos',[0.9 0.95 0.1 0.05]);
  set(h_done,'string','Done','callback',@DeleteDoneButton)
  set(h_fig,'WindowStyle','modal')
  waitfor(h_done)
  if ishandle(h_fig)
    h_lev=getappdata(gcf,'h_lev');
    level=mean(get(h_lev,'xdata'));
    delete(h_fig)
  else
    warning('User aborted; no return value.')
    level=[];
  end
end

function DeleteDoneButton(hObject,varargin)
delete(hObject)

function ResetOriginalLevel(hObject,varargin)
h_lev = getappdata(gcf,'h_lev');
init_level = get(h_lev,'UserData');
set(h_lev,'XData',init_level*[1 1])
text_obj = findobj('Type','Text','Tag','cbar_text');
movex_text(text_obj,init_level)
update_plot

function update_plot
im=getappdata(gcf,'im');
h_lev=getappdata(gcf,'h_lev');
my_level=mean(get(h_lev,'xdata'));
h_ax2=getappdata(gcf,'h_ax2');
h_im2=findobj(h_ax2,'type','image');
%segmented image
bw=(im>my_level);
rgb_version=repmat(double(bw),[1 1 3]);
set(h_im2,'cdata',rgb_version)


text:
% Add land and sea-ice grid boxes to the reconstructed fields
function Yrec_out = metoffice_add_land( Yrec, dataset )

[ mask, lat, lon ] = metoffice_get_mask( dataset );

% path from system
datapath = metoffice_get_path();
load( [ datapath, dataset ], 'time' )

% Add land and ice areas to the reconstruction fields
Yrec_out = NaN( length(lat)*length(lon), length(time) );
Yrec_out( logical(mask(:)), : ) = Yrec;

function Y = metoffice_add_land_old(X, maskfile)

S = load('/share/climate/data/UK_Met_Office/RecTest/mask');
Y = nan(numel(S.mask),size(X,2));
Y(S.mask(:)==1,:) = X;
% $$$ S = load('/share/climate/data/UK_Met_Office/RecTest/land_hadgem');
% $$$ land = (S.land(:) >= 0.5);
% $$$ Y = nan(length(land),size(X,2));
% $$$ Y(~land,:) = X;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function ABCtranslate = abc2num(S)
%ABC2NUM convert a spreadsheet label string to corresponding column number.
%   ABC2NUM(S) returns the number corresponding to the spreadsheet column
%   represented by the string, S.
%
%   The input string, S, may contain the ascii letters A-Z, in the format
%   commonly employed  by spreadsheet programs such as Calc and Microsoft
%   Excel. Other characters, including lower case letters, a-z, are
%   ignored.
%
%   Examples
%       S = 'CS';
%       abc2num(S)  returns  97
%
%     See also num2abc, num2str, str2num, char.

alphabetStrings = ['[A,B,C,D,E,F,G,'...
                   'H,I,J,K,L,M,N,'...
                   'O,P,Q,R,S,T,U,'...
                   'V,W,X,Y,Z]'];

n = regexp(S,alphabetStrings,'start');
ascii = zeros(size(n));
for i=1:length(n)
    ascii(i) = double(S(n(i)));
end
digits = ascii - 64;

ABCtranslate = polyval(flipdim(digits,1),26);

%==========================================================================
%abc2num.m
%Function that converts a spreadsheet column label to a number.
%--------------------------------------------------------------------------
%author:  David Sedarsky
%date:    February, 2012


text:
function linefit
zoom off
k = waitforbuttonpress;
point1 = get(gca,'CurrentPoint');    % button down detected
finalRect = rbbox;                   % return Figure units
point2 = get(gca,'CurrentPoint');    % button up detected
point1 = point1(1,1:2);              % extract x and y
point2 = point2(1,1:2);
p1 = min(point1,point2);             % calculate locations
offset = abs(point1-point2);         % and dimensions
x1 = p1(1);
x2 = p1(1)+offset(1);

zoom on
linehandle = findobj(gca, 'Type', 'line');
if ~isempty(linehandle)
    xdata = get(linehandle(1), 'XData');
    ydata = get(linehandle(1), 'YData');
    ind = find(xdata >=x1  & xdata <= x2);
    [p,s] = polyfit(xdata(ind), ydata(ind), 1);
    if p(2) >= 0
        msgbox(['y = ' num2str(p(1)) ' * x + ' num2str(p(2))], 'Line fit');
    else
        msgbox(['y = ' num2str(p(1)) ' * x - ' num2str(abs(p(2)))], 'Line fit');
    end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


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

    image([L R] * 255.0);
    w = numcols(L);
    colormap(gray(256))

    h = line(0, 0);

    while 1
        [x,y] = ginput(1);
        if isempty(x)
            break;
        end
        x
        y

        if x <= w
            disp('left image');
        else
            disp('right image')
        end

        p = [x y 1]';
        l = F * p;

        x1 = 1;
        x2 = w;

        y1 = (-l(1)*x1 - l(3)) / l(2);
        y2 = (-l(1)*x2 - l(3)) / l(2);

        set(h, 'Xdata', [x1+w x2+w], 'Ydata', [y1 y2]);
    end

author: petercorke

repo: machinevision-toolbox-matlab