text:
function varargout = process_evt_detect_threshold( varargin )
% PROCESS_EVT_DETECT_THRESHOLD: Event detection based on a set threshold for a group of recordings file
%
% USAGE:  OutputFiles = process_evt_detect_threshold('Run', sProcess, sInputs)
%                 evt = process_evt_detect_threshold('Compute', F, TimeVector, OPTIONS)
%             OPTIONS = process_evt_detect_threshold('Compute')                         % Get the default options structure

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Elizabeth Bock, Francois Tadel, 2015-2016
%          Raymundo Cassani, 2021

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Detect events above threshold';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 45;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle?highlight=%28Detect+events+above+threshold%29#Artifact_correction_with_SSP';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data'};
    sProcess.OutputTypes = {'raw', 'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
	% Notice
    sProcess.options.notice.Comment = ['This process creates extended events to indicate all the time
' ...
                                       'samples where the signal is above a given amplitude threshold.
' ...
                                       'For creating simple events, use one of the following processes:
' ...
                                       '"Detect custom events" or "Detect analog triggers".

'];
    sProcess.options.notice.Type    = 'label';
    % Event name
    sProcess.options.eventname.Comment = 'Event name: ';
    sProcess.options.eventname.Type    = 'text';
    sProcess.options.eventname.Value   = 'artifact';
    % Separator
    sProcess.options.separator.Type = 'separator';
    sProcess.options.separator.Comment = ' ';
    % Channel name
    sProcess.options.channelname.Comment = 'Channel name: ';
    sProcess.options.channelname.Type    = 'channelname';
    sProcess.options.channelname.Value   = '';
    % Channel name comment
    sProcess.options.channelhelp.Comment = 'You can use the montage syntax here: "ch1, -ch2"';
    sProcess.options.channelhelp.Type    = 'label';
    % Time window
    sProcess.options.timewindow.Comment = 'Time window:';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    % Blanking period
    sProcess.options.blanking.Comment = 'Min duration between two events: ';
    sProcess.options.blanking.Type    = 'value';
    sProcess.options.blanking.Value   = {0.02, 'ms', []};
    % Threshold maximum
    sProcess.options.thresholdMAX.Comment = 'Maximum threshold: ';
    sProcess.options.thresholdMAX.Type    = 'value';
    sProcess.options.thresholdMAX.Value   = {0.00, '', 2};
    % units
    sProcess.options.label1.Comment = '
Threshold Units:';
    sProcess.options.label1.Type    = 'label';
    sProcess.options.units.Comment = {'None: 100', 'mV: 10-3', 'uV: 10-6', 'fT: 10-15', ''};
    sProcess.options.units.Type    = 'radio_line';
    sProcess.options.units.Value   = 1;
    % Filter
    sProcess.options.label2.Comment = '
Filter signal before detection:';
    sProcess.options.label2.Type    = 'label';
    sProcess.options.bandpass.Comment = 'Frequency band: ';
    sProcess.options.bandpass.Type    = 'range';
    sProcess.options.bandpass.Value   = {[], 'Hz', 2};
    % Use absolute value
    sProcess.options.isAbsolute.Comment = 'Use absolute value of signal';
    sProcess.options.isAbsolute.Type    = 'checkbox';
    sProcess.options.isAbsolute.Value   = 0;
    % Remove DC
    sProcess.options.isDCremove.Comment = 'Remove DC offset';
    sProcess.options.isDCremove.Type    = 'checkbox';
    sProcess.options.isDCremove.Value   = 0;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = ['Detect: ', sProcess.options.eventname.Value];
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % ===== GET OPTIONS =====
    % Event name
    evtName = strtrim(sProcess.options.eventname.Value);
    chanName = strtrim(sProcess.options.channelname.Value);
    if isempty(evtName) || isempty(chanName)
        bst_report('Error', sProcess, [], 'Event and channel names must be specified.');
        OutputFiles = {};
        return;
    end
    % Ignore bad segments
    % => Consider that if the event name contains "bad", we need to include the bad segments. If not, we ignore them.
    isIgnoreBad = panel_record('IsEventBad', evtName);
    % Prepare options structure for the detection function
    OPTIONS = Compute();
    OPTIONS.thresholdMAX = sProcess.options.thresholdMAX.Value{1};
    OPTIONS.isAbsolute   = sProcess.options.isAbsolute.Value;
    OPTIONS.isDCremove   = sProcess.options.isDCremove.Value;
    OPTIONS.bandpass     = sProcess.options.bandpass.Value{1};
    iUnits = sProcess.options.units.Value;
    switch iUnits
        case 1
            OPTIONS.thresholdMAX = OPTIONS.thresholdMAX * 1;
        case 2
            OPTIONS.thresholdMAX = OPTIONS.thresholdMAX * 1e-3;
        case 3
            OPTIONS.thresholdMAX = OPTIONS.thresholdMAX * 1e-6;
        case 4
            OPTIONS.thresholdMAX = OPTIONS.thresholdMAX * 1e-15;
    end

    % Time window to process
    if isfield(sProcess.options, 'timewindow') && isfield(sProcess.options.timewindow, 'Value') && iscell(sProcess.options.timewindow.Value) && ~isempty(sProcess.options.timewindow.Value)
        TimeWindow = sProcess.options.timewindow.Value{1};
    else
        TimeWindow = [];
    end

    % Blanking period between two events
    if isfield(sProcess.options,'blanking')
        OPTIONS.blanking = sProcess.options.blanking.Value{1};
    end

    % Option structure for function in_fread()
    ImportOptions = db_template('ImportOptions');
    ImportOptions.ImportMode      = 'Time';
    ImportOptions.UseCtfComp      = 1;
    ImportOptions.UseSsp          = 1;
    ImportOptions.EventsMode      = 'ignore';
    ImportOptions.DisplayMessages = 0;
    ImportOptions.RemoveBaseline  = 'no';

    % Get current progressbar position
    progressPos = bst_progress('get');
    nEvents = 0;
    nTotalOcc = 0;

    % For each file
    iOk = false(1,length(sInputs));
    for iFile = 1:length(sInputs)
        % ===== GET DATA =====
        % Progress bar
        bst_progress('text', 'Reading channel to process...');
        bst_progress('set', progressPos + round(iFile / length(sInputs) / 3 * 100));
        % Load the raw file descriptor
        isRaw = strcmpi(sInputs(iFile).FileType, 'raw');
        if isRaw
            DataMat = in_bst_data(sInputs(iFile).FileName, 'F', 'Time');
            sFile = DataMat.F;
        else
            DataMat = in_bst_data(sInputs(iFile).FileName, 'Time');
            sFile = in_fopen(sInputs(iFile).FileName, 'BST-DATA');
        end
        % Load channel file
        ChannelMat = in_bst_channel(sInputs(iFile).ChannelFile);
        % Process only continuous files
        if ~isempty(sFile.epochs)
            bst_report('Error', sProcess, sInputs(iFile), 'This function can only process continuous recordings (no epochs).');
            continue;
        end

        % Get channel to process: multiple channels
        if any(chanName == ',')
            [iChannels, iChanWeights] = process_evt_detect('ParseChannelMontage', chanName, {ChannelMat.Channel.Name});
            if isempty(iChannels)
                bst_report('Error', sProcess, sInputs(iFile), ['Montage "' chanName '" could not be interpreted. Please check channel names.']);
                continue;
            end
        % One channel
        else
            iChannels = find(strcmpi({ChannelMat.Channel.Name}, chanName));
            if isempty(iChannels)
                bst_report('Error', sProcess, sInputs(iFile), ['Channel "' chanName '" not found in the channel file.']);
                continue;
            elseif (length(iChannels) > 1)
                bst_report('Error', sProcess, sInputs(iFile), ['Found more than one channel with name "' chanName '" in the channel file.']);
                continue;
            end
            iChanWeights = 1;
        end

        % Read channel to process
        if ~isempty(TimeWindow)
            SamplesBounds = round(sFile.prop.times(1) .* sFile.prop.sfreq) + bst_closest(TimeWindow, DataMat.Time) - 1;
        else
            SamplesBounds = [];
        end
        [F, TimeVector] = in_fread(sFile, ChannelMat, 1, SamplesBounds, iChannels, ImportOptions);
        % Apply weights if reading multiple channels
        if (length(iChannels) > 1)
            F = iChanWeights * F;
        end
        % If nothing was read
        if isempty(F) || (length(TimeVector) < 2)
            bst_report('Error', sProcess, sInputs(iFile), 'Time window is not valid.');
            continue;
        end

        % ===== DETECT PEAKS =====
        % Progress bar
        bst_progress('text', 'Detecting peaks...');
        bst_progress('set', progressPos + round(2 * iFile / length(sInputs) / 3 * 100));
        % Perform detection
        detectedEvt = Compute(F, TimeVector, OPTIONS);

        % ===== CREATE EVENTS =====
        sEvent = [];
        % Basic events structure
        if ~isfield(sFile, 'events') || isempty(sFile.events)
            sFile.events = repmat(db_template('event'), 0);
        end
        % Process each event type separately
        for i = 1:length(detectedEvt)
            % Event name
            if (i > 1)
                newName = sprintf('%s%d', evtName, i);
            else
                newName = evtName;
            end
            % Get the event to create
            iEvt = find(strcmpi({sFile.events.label}, newName));
            % Existing event: reset it
            if ~isempty(iEvt)
                sEvent = sFile.events(iEvt);
                sEvent.epochs  = [];
                sEvent.times   = [];
                sEvent.reactTimes = [];
            % Else: create new event
            else
                % Initialize new event
                iEvt = length(sFile.events) + 1;
                sEvent = db_template('event');
                sEvent.label = newName;
                % Get the default color for this new event
                sEvent.color = panel_record('GetNewEventColor', iEvt, sFile.events);
            end
            % Times, samples, epochs
            sEvent.times    = detectedEvt{i};
            sEvent.epochs   = ones(1, size(sEvent.times,2));
            sEvent.channels = [];
            sEvent.notes    = [];
            % Add to events structure
            sFile.events(iEvt) = sEvent;
            nEvents = nEvents + 1;
            nTotalOcc = nTotalOcc + size(sEvent.times, 2);
        end

        % ===== SAVE RESULT =====
        % Progress bar
        bst_progress('text', 'Saving results...');
        bst_progress('set', progressPos + round(3 * iFile / length(sInputs) / 3 * 100));
        % Only save changes if something was detected
        if ~isempty(sEvent)
            % Report changes in .mat structure
            if isRaw
                DataMat.F = sFile;
            else
                DataMat.Events = sFile.events;
            end
            DataMat = rmfield(DataMat, 'Time');
            % Save file definition
            bst_save(file_fullpath(sInputs(iFile).FileName), DataMat, 'v6', 1);
            % Report number of detected events
            bst_report('Info', sProcess, sInputs(iFile), sprintf('%s: %d events detected in %d categories', chanName, nTotalOcc, nEvents));
        else
            bst_report('Warning', sProcess, sInputs(iFile), ['No event detected on channel "' chanName '". Please check the signal quality.']);
        end
        iOk(iFile) = true;
    end
    % Return all the input files
    OutputFiles = {sInputs(iOk).FileName};
end

%% ===== PERFORM DETECTION =====
% USAGE:      evt = Compute(F, TimeVector, OPTIONS=[])
%         OPTIONS = Compute()                            % Get the default options structure
function evt = Compute(F, TimeVector, OPTIONS)
    % Options structure
    defOptions = struct('thresholdMAX',  0, ...     % Maximum threshold value
                        'bandpass', [], ...         % bandpass filter data before detection
                        'isAbsolute', 0, ...        % If 1, the max point defines the event, else, first thresh crossing defines the event
                        'isDCremove', 0, ...        % If 1, remove the DC offset of the signal before detection
                        'blanking', 0.02);          % Min duration between two events. Default 0.02 s
    % Parse inputs
    if (nargin == 0)
        evt = defOptions;
        return;
    end
    % Copy the missing parameters
    OPTIONS = struct_copy_fields(OPTIONS, defOptions, 0);
    % Sampling frequency
    sFreq = 1 ./ (TimeVector(2) - TimeVector(1));
    % Initialize output
    evt = {};

    %% ===== DETECT EVENTS ABOVE THRESHOLD =====
     % Filter recordings
    if ~isempty(OPTIONS.bandpass)
        F = process_bandpass('Compute', F, sFreq, OPTIONS.bandpass(1), OPTIONS.bandpass(2), 'bst-hfilter-2019', 1);
    end

    % Remove DC offset with detrending
    if OPTIONS.isDCremove
        % Remove the linear trend
        nTime = length(F);
        iTime = 1:nTime;
        % Basis functions
        x = [ones(1,nTime); 0:nTime-1];
        % Estimate the contribution of the basis functions
        % beta = dat(:,iTime)/x(:,iTime); <-this leads to numerical issues, even in simple examples
        invxcov = inv(x(:,iTime) * x(:,iTime)');
        beta    = F(:,iTime) * x(:,iTime)' * invxcov;
        % Remove the estimated basis functions
        F = F - beta*x;
    end

    % Absolute valuse
    if OPTIONS.isAbsolute
        F = abs(F);
    end

    % Threshold mask
    threshMask = zeros(size(F));
    % Find all the indices that are above the threshold
    threshMask(F > OPTIONS.thresholdMAX) = 1;

    % Nothing detected: exit
    if sum(threshMask) < 1
        return;
    end

    % Group time points into extended events
    % group events that occur withing OPTIONS.blanking seconds of each other
    minNewEvent = round(sFreq * OPTIONS.blanking);
    % If minNewEvent samples is longer than the signal to process: exit
    if (minNewEvent >= length(F))
        bst_report('Warning', 'process_evt_detect_threshold', [], 'The minimun period between two events is longer than the signal. Cannot perform detection.');
        return;
    end
    diffThreshMask = diff([0 threshMask 0]);
    startEve = find(diffThreshMask == 1);
    endEve = find(diffThreshMask == -1) -1;
    iEve = 1;
    while iEve < length(startEve)-1
        if startEve(iEve + 1) - endEve(iEve) < minNewEvent
            % combine the two events
            endEve(iEve) = endEve(iEve + 1);
            startEve(iEve+1) = [];
            endEve(iEve+1) = [];
        else
            iEve = iEve + 1;
        end
    end

    % Find times of the events
    evt = {[TimeVector(startEve); TimeVector(endEve)]};

end

author: brainstorm-tools


text:
function [Vc,R]=circlefit(V)

% function [Vc,R]=circlefit(V)
% ------------------------------------------------------------------------
% This function fits a circle to 2 or more points. The input is an array
% containing the points (vertices) V in 2D or 3D. Although 3D data is
% supported the Z-coordinates should all be equal such that the fitting
% problem remains 2D.
% The output consists of the circle centre Vc and radius R.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
% 2019/11/21 Added to GIBBON
% ------------------------------------------------------------------------

%% Parse input

%Check if empty or not numeric
if isempty(V)
    error('Empty input')
elseif ~isnumeric(V)
    error('Input is not numeric');
end

%Force real
if ~isreal(V)
    V=real(V);
    warning('Complex data encountered. Complex parts ignored');
end

%Force double
if ~isa(V,'double')
    V=double(V); %Attempt conversion to double
end

%Check size
sizV=size(V);
if sizV(1)==1
    error('Only 1 point provided')
end

%Force 3D
if sizV(2)==2
    V(:,3)=0;
end

%% Fit cirlce
Vm=mean(V,1); %Get mean of point set
VD=V-Vm(ones(sizV(1),1),:); % shift points on mean
D=sqrt(sum(VD.^2,2)); %Compute distance to centre

if sizV(1)==2 %If only two points provided
    Vc=Vm; %Cirlc centre set to mean
    R=mean(D); %Radius set to mean distance
else %Fit circle to more than 2 points
    % Scale data
    scaleFactor=max(max(D,eps));
    VD=VD./scaleFactor; D=D./scaleFactor;

    % Solve sytem
    VF=[2*VD(:,1) 2*VD(:,2) ones(sizV(1),1)]\D.^2;

    % Unscale, collect output
    Q=[Vm(1) Vm(2) 0]+[VF(1) VF(2) sqrt(VF(3)+hypot(VF(1),VF(2)).^2)]*scaleFactor;

    Vc=[Q(1) Q(2) Vm(3)]; %Cirle centre
    R=Q(3); %Radius
end

if sizV(2)==2
    Vc=Vc(:,[1 2]);
end

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function val = get_param_val(param_list, group_name,param_name)

[col, row] = size(param_list);

for i = 1:row
   if strcmp(param_list(i).name, group_name)
      for j = 1:param_list(i).num_param
          if strcmp(param_list(i).param(j).name, param_name)
              val = param_list(i).param(j).val;
              return;
          end
      end
   end
end
% do not find parameter, return nan
val = NaN;

end

author: Firmament-Autopilot

repo: FMT-Model


text:
classdef HomogenizationTests < handle & matlab.unittest.TestCase

    properties (TestParameter)
            errorTests = {...
                'testNumericalConvergenceOfNumberOfLaminates';
                'testCommutingHomogPlaneStressWithZeroPoisson';
                'testCommutingVoigtHomog';
                'testAnisotropicPlaneStressbyEnergyEquivalence';
                'testStressInPlaneStress';
                'testStressRotationInVoigtNotationIn3D';
                'testStressRotationInVoigtNotationInPlaneStress';
                'testInverseSymmetricFourthOrderTensor' ;
                'testInverseNonSymmetricFourthOrderTensor';
                'testInverseOfInverseForStiffTensor';
                'testIsotropicFourthOrderTensor'
                'testSymmetrizeIsotropicFourthOrderTensor';
                'testSymmetryForIAniTensorInVoigt'
                'testMakeAnisotorpicTensorPlaneStressSymbolically';
                'testEnergyEquivalenceVoigtAndTensorNotationForIsoTensor';
                'testEnergyEquivalenceVoigtAndTensorNotationForIAniTensor';
                'testComplianceTensorThrougtVoigtComparingEnergy';
                'TestTwoRankSequentialLaminate';
                'testHorizontalTensorRotatedVsVPH';
                }
            passedTests = {...
                'testDiagonalLaminate';
                'testHorizontalLaminate';
                'TestGeneralTwoRankSequentialLaminate';
                'testNotCommutingHomogPlaneStress';
                'testSymmetrizeFourthOrderTensor';
                'testHorizontalTensorRotatedVsVHP';
                'testHorizontalTensorRotatedVsRank2';
                'testHorizontalTensorRotatedVsHVP';
                }
            elementdiffreact = {'testHorizontalTensorRotatedVsHVP'};
    end

%     methods (Test, TestTags = {'HomogenizationTests', 'Elementdiffreact'})
%
%         function testsElementdiffreact(testCase, elementdiffreact)
%             testCase.fixFolder();
%             test = eval(elementdiffreact);
%             passed = test.hasPassed();
%             verifyTrue(testCase, passed)
%         end
%
%     end

    methods (Test, TestTags = {'HomogenizationTests', 'ShowingError'})

        function testsError(testCase, errorTests)
            testCase.fixFolder();
            test = eval(errorTests);
            err = test.computeError();
            tol = test.tol;
            testCase.verifyLessThanOrEqual(err, tol)
        end

    end

    methods (Test, TestTags = {'HomogenizationTests', 'NotShowingError'})

        function testsPassed(testCase, passedTests)
            testCase.fixFolder();
            test = eval(passedTests);
            passed = test.hasPassed();
            verifyTrue(testCase, passed)
        end

    end

    methods (Access = private)

        function fixFolder(testCase)
            import matlab.unittest.fixtures.CurrentFolderFixture
            changeToFolder = '../../';
            testCase.applyFixture(CurrentFolderFixture(changeToFolder));
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%-------------------------------------------------------------------------
% This is a simulation demo to show multiple trackpaths for tracking the
% blood cell motion. In this demo, the diameter of the blood cell is 16
% pixels, so we use 9 trackpaths(in row 23) and the distance between each
% trackpath is 8 pixels(in row 25).
% Written by Yuan Chen,Nanjing university of aeronautic and astronautic(2010)
%-------------------------------------------------------------------------
clear all;
%% read the avi file, and calculate the mean image;
display('Reading and processing demo3.avi file...')
video=aviread('demo3.avi');             % read the avi file;
video = {video.cdata};                  % exracte the avi data(matrix);
z = 0;
for i=1:length(video);                  % change 'color' video to gray;
    video{i}=rgb2gray(video{i});
    z = z + double(video{i});
end
z = z/length(video);                    % the mean image;
[mz,nz] = size(z);
[mu1,mu2,v1x,v1y,v2x,v2y] = eigfunction(z,7,5);  % calculate the eigen values and eigen vectors of the mean image;

%% the trackpaths generation;
N = 9;                                  % THE number of trackpath can be set within [5,65];NOTE that,larger N will consume
                                        % much more processing times,we suggrest that N = 9, and no larger than 16;
t = 64/(N-1);                           % the distance between each trakcpath;
T = 64:t:128;                           % T are the radius of the trackpaths(the center(148,150)in the mean image z);
for k = 1:length(T)
    tmp = zeros(size(z));
    for i = 1:138
        for j = 1:300
            if sqrt((148-i)^2 + (150-j)^2 ) < T(k)+1 && sqrt((148-i)^2 + (150-j)^2 ) >= T(k);
               tmp(i,j) = 255;
            end
        end
    end
    tmp = im2bw(tmp);
    tmp = bwmorph(tmp,'thin',inf);
    trackpath{k} = tmp;
end
display([num2str(N),' trackpaths are selected to generate ',num2str(N),' ST images...Processing them may take several minutes'])
%% ST image generation, processing and traces extraction;
for i=1:size(trackpath,2)                                   % the number of the generated trackpath;
    map = [];
    [xx,yy] = startpoint(trackpath{i});                     % finding the starting point of a trackpath;
    ind_trackpath = ord_line_indx(trackpath{i},xx,yy);      % order the trackpath points;
    for j = 1:size(video,2)                                 % ST image generation;
        line = video{j}(ind_trackpath);
        map = [map;line];
    end
    map = map';
    Map{i} = map;                                           % store the ST image;
    display(['processing ST image',int2str(i)])
    [F,Theta,angle] = J4(map,2.5);                          % raw ST image enhanced by Jacob filter;
    R = nsf(F);                                             % applying noise supression function;
    K = angfilter(R,Theta,angle,20);                        % applying orientation filter function;
    th = graythresh(K);                                     % threshold for cell tracked trackpath;
    if th <= 0.01;                                          % threshold for none cell tracked trackpath;
        th = 0.01;
    end
    bw = im2bw(K,th);
    thin = bwmorph(bw,'thin',inf);                          % thinning extracted traces;
    trace{i} = bwareaopen(thin,10);                         % remove noise trace if its length smaller than 10 pixels;
end

%% calculate the trackpaths coordinate and grayscale;
h = fspecial('gaussian',7);
for i = 1:length(trackpath)                                 % denoise;
    Map_denoise{i} = conv2(Map{i},h,'same');
    Map_denoise{i} = conv2(Map_denoise{i},h,'same');
end
% -------find out the coordinate(X_trackpath,Y_trackpath) and grayscale of the trackpaths;
for i = 1:length(trackpath)
    [xx,yy] = startpoint(trackpath{i});                     % find the start coordinate of a trackpath;
    ind_trackpath = ord_line_indx(trackpath{i},xx,yy);      % search the trackpath's point from the start point orderly;
    for j = 1:length(ind_trackpath)
        [X_trackpath(j),Y_trackpath(j)] = ind2sub([mz,nz],ind_trackpath(j)); % change the index to (X,Y) coordinate;
    end
    Trackpath_value{i}(:,1) = X_trackpath';
    Trackpath_value{i}(:,2) = Y_trackpath';
    Trackpath_value{i}(:,3) = 0;
    [Y_trace,X_trace] = find(trace{i}==1);                  % find extracted trace coordinate;
    for k = 1:length(Y_trace)                               % store the grayscale of the trace;
        Trackpath_value{i}(Y_trace(k),3) = Map_denoise{i}(Y_trace(k),X_trace(k));
    end
    clear X_trackpath;clear Y_trackpath;
end
% -------backward mapping the extracted traces to the trackpaths;
for i = 1:length(Trackpath_value)
    for j = 1:size(Trackpath_value{i},1)
        z(Trackpath_value{i}(j,1),Trackpath_value{i}(j,2)) = Trackpath_value{i}(j,3);
    end
end
% figure,imshow(z,[]),title('The tracking state of trackpaths')

%% trackpaths alignment
k = round(length(T)/2);                                     % the centerline trackpath;
% -------turn all trackpaths and all ST images to the same length(or size) as the centreline trackpath and centreline ST image
% according to the eigen victor v1x,v1y;
for i = 1:length(T)
    p = findpoints(trackpath{k},trackpath{i},v1x,v1y);      % the corresponding relationship between the k and the i trackpaths.
    P{i} = round(smooth(p));
    if i~=k
        for j = 1:length(P{i})
            trace_aligned{i}(j,:) = trace{i}(P{i}(j),:);
            Points_aligned{i}(j,:) = Trackpath_value{i}(P{i}(j),:);
        end
    else
        trace_aligned{i} = trace{i};
        Points_aligned{i} = Trackpath_value{i};
    end
end

% -------base trace generation;
B_trace = 0;
for i = 1:length(trace_aligned)
    B_trace = B_trace + trace_aligned{i};
end
% figure,imshow(B_trace),title('the result of traces mapping')
B_trace = im2bw(B_trace);
B_trace = bwmorph(B_trace,'dilate');
B_trace = bwmorph(B_trace,'thin',inf);
% figure,imshow(B_trace),title('Base-trace')

%%  fusion of aligned points in trackpaths to calculate the trajectory;
l = 1;
for j = 1:size(Points_aligned{1},1)             % the number of the aligned points in the trackpath;
    Fusion_points = [];
    for i = 1:size(Points_aligned,2)            % the number of the trackpaths;
        if Points_aligned{i}(j,3) > 0           % >0 means the cell is tracked;
            X = Points_aligned{i}(j,1);         % X coordinate;
            Y = Points_aligned{i}(j,2);         % Y coordinate;
            V = Points_aligned{i}(j,3);         % grayvalue of (X,Y);
            Fusion_points = [Fusion_points,[X;Y;V]];
        end
    end
    if size(Fusion_points,2) > 2                % the cell is tracked by three or more than three trackpths;
        X_fused(l,1) = mean(Fusion_points(1,:));% the fused position is the center among the aligned points;
        Y_fused(l,1) = mean(Fusion_points(2,:));
        l = l+1;
    elseif size(Fusion_points,2) == 2           % the cell is tracked by two trackpaths(Eq.(15)in the paper);
        x1 = Fusion_points(1,1); x2 = Fusion_points(1,2);
        y1 = Fusion_points(2,1); y2 = Fusion_points(2,2);
        v1 = Fusion_points(3,1); v2 = Fusion_points(3,2);
        [X_fused(l,1),Y_fused(l,1)] = fusion2p(x1,y1,x2,y2,v1,v2,0.4); % the fused position calculation;
        l = l+1;
    elseif size(Fusion_points,2) ==1            % the cell is tracked by only one trackpath;
        X_fused(l,1) = Fusion_points(1,1);      % no fusion is needed;
        Y_fused(l,1) = Fusion_points(2,1);
        l = l+1;
    end
end

   X_trajectory = round(X_fused);               % calculated trajectory;
   Y_trajectory = round(Y_fused);
   X_smooth = round(smooth(X_fused,10));        % smooth process to prevent jagged trajectory;
   Y_smooth = round(smooth(Y_fused,10));
for i = 1:length(X_fused)
    z(X_trajectory(i),Y_trajectory(i)) = 255;
end
figure,imshow(z,[]);title(['Blood cell tracking by ',num2str(N), ' trackpaths(black line)'])
hold on;plot(Y_smooth,X_smooth,'r')

% -------the real trajectory of the cell;
x = 76:228;
y = 124*sin(pi/162*x+80.32);
y = round(-y+151);
hold on,plot(x,y);
h = legend('Tracked trajectory','Real trajectory',2);

%% calculating the error and variance;
for i = 1:length(X_trajectory)                  % by finding the closest points in the real trajectory to calculate the error;
    for j = 1:length(x)
        dis(j) = sqrt((X_trajectory(i) - y(j))^2 + (Y_trajectory(i) - x(j))^2);
        Error(i) = min(dis);
    end
end
ES = [mean(Error),std2(Error)];                 % mean error and standard deviation;
display(['The Error and Variance of the tracked tracjectory by using ',num2str(N),' trackpaths are'])
display([num2str(ES(1)) ' pixels and ' num2str(ES(2)),' pixels'])

author: Sable

repo: mcbench-benchmarks


text:
function [objfiles, timestamp] = cs_make (f, docomplex)
%CS_MAKE compiles CXSparse for use in MATLAB.
%   Usage:
%       cs_make
%       [objfiles, timestamp] = cs_make (f, docomplex)
%
%   With no input arguments, or with f=0, only those files needing to be
%   compiled are compiled (like the Unix/Linux/GNU "make" command, but not
%   requiring "make").  If f is a nonzero number, all files are compiled.
%   If f is a string, only that mexFunction is compiled.  For example,
%   cs_make ('cs_add') just compiles the cs_add mexFunction.  This option is
%   useful when developing a single new mexFunction.  This function can only be
%   used if the current directory is CXSparse/MATLAB/CSparse.  Returns a list of
%   the object files in CXSparse, and the latest modification time of any source
%   codes.
%
%   NOTE: if your compiler does not support the ANSI C99 complex type (most
%   notably Microsoft Windows), the CXSparse mexFunctions will not support
%   complex sparse matrices.  The complex case is not attempted if docomplex is
%   zero.
%
%   To add a new function and its MATLAB mexFunction to CXSparse:
%
%       (1) Create a source code file CXSparse/Source/cs_mynewfunc.c.
%       (2) Create a help file, CXSparse/MATLAB/CSparse/cs_mynewfunc.m.
%           This is very useful, but not strictly required.
%       (3) Add the prototype of cs_mynewfunc to CXSparse/Include/cs.h.
%       (4) Create its MATLAB mexFunction, CXSparse/MATLAB/cs_mynewfunc_mex.c.
%       (5) Edit cs_make.m, and add 'cs_mynewfunc' to the 'cs' and 'csm' lists.
%       (6) Type 'cs_make' in the CXSparse/MATLAB/CSparse directory.
%           If all goes well, your new function is ready for use in MATLAB.
%
%       (7) Optionally add 'cs_mynewfunc' to CXSparse/Source/Makefile
%           and CXSparse/MATLAB/CSparse/Makefile, if you want to use the
%           Unix/Linux/GNU make command instead of cs_make.m.  See where
%           'cs_add' and 'cs_add_mex' appear in those files, and add
%           'cs_mynewfunc' accordingly.
%       (8) Optionally add 'cs_mynewfunc' to Tcov/Makefile, and add additional
%           test code to cs_test.c, and add MATLAB test code to MATLAB/Test/*.
%
%   Example:
%       cs_make                  % compile everything
%       cs_make ('cs_chol') ;    % just compile cs_chol mexFunction
%
%   See also MEX.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

if (nargin < 1)
    f = 0 ;
end
if (nargin < 2)
    docomplex = 1 ;
end

try
    % ispc does not appear in MATLAB 5.3
    pc = ispc ;
catch
    % if ispc fails, assume we are on a Windows PC if it's not unix
    pc = ~isunix ;
end

if (pc)
    docomplex = 0 ;
end

if (docomplex == 0)
    % do not attempt to compile with complex matrices
    [objfiles, timestamp] = cs_make_helper (f, 0) ;
else
    try
        % try with complex support
        [objfiles, timestamp] = cs_make_helper (f, 1) ;
    catch
        % oops - that failed, try without complex support
        fprintf ('retrying without complex matrix support\n') ;
        [objfiles, timestamp] = cs_make_helper (f, 0) ;
    end
end

if (f > 0)
    fprintf ('CXSparse successfully installed.\n') ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Itheta,SigmaTheta,deltaMuTheta,suffStat] = VBA_Itheta(theta,y,posterior,suffStat,dim,u,options)
% Gauss-Newton update of the evolution parameters

if options.DisplayWin % Display progress
    set(options.display.hm(1),'string',...
        'VB Gauss-Newton on evolution parameters... ');
    set(options.display.hm(2),'string','0%');
    drawnow
end

%  Look-up which evolution parameter to update
indIn = options.params2update.theta;
indInx = options.params2update.x;

% Get precision parameters
alphaHat = posterior.a_alpha./posterior.b_alpha;

% Preallocate intermediate variables
iQx = options.priors.iQx;
Q = options.priors.SigmaTheta(indIn,indIn);
iQ = VBA_inv(Q);
muTheta0 = options.priors.muTheta;
Theta = muTheta0;
Theta(indIn) = theta;
dtheta0 = muTheta0-Theta;
dx = zeros(dim.n,dim.n_t);
div = 0;

%--- Initial condition ---%

% evaluate evolution function at current mode
[fx,dF_dX,dF_dTheta] = VBA_evalFun('f',posterior.muX0,Theta,u(:,1),options,dim,1);

% check infinite precision transition pdf
iQ2 = VBA_inv(iQx{1},indInx{1},'replace');

% posterior covariance matrix terms
d2fdx2 = dF_dTheta*iQ2*dF_dTheta';

% error terms
dx(:,1) = (posterior.muX(:,1) - fx);
dx2 = dx(:,1)'*iQ2*dx(:,1);
ddxdtheta = dF_dTheta*iQ2*dx(:,1);

%--- Loop over time series ---%
for t=1:dim.n_t-1

    % check infinite precision transition pdf
    iQ2 = VBA_inv(iQx{t+1},indInx{t+1},'replace');

    % evaluate evolution function at current mode
    [fx,dF_dX,dF_dTheta] = VBA_evalFun('f',posterior.muX(:,t),Theta,u(:,t+1),options,dim,t+1);

    % posterior covariance matrix terms
    d2fdx2 = d2fdx2 + dF_dTheta*iQ2*dF_dTheta';

    % error terms
    dx(:,t+1) = (posterior.muX(:,t+1) - fx);
    dx2 = dx2 + dx(:,t+1)'*iQ2*dx(:,t+1);
    ddxdtheta = ddxdtheta + dF_dTheta*iQ2*dx(:,t+1);

    % Display progress
    if options.DisplayWin && mod(t,dim.n_t./10) < 1
        set(options.display.hm(2),'string',[num2str(floor(100*t/dim.n_t)),'%']);
        drawnow
    end

    % Accelerate divergent update
    if VBA_isWeird ({dx2, dF_dX, dF_dTheta})
        div = 1;
        break
    end

end

if options.DisplayWin % Display progress
    set(options.display.hm(2),'string','OK');
    drawnow
end

% posterior covariance matrix
iSigmaTheta = iQ + alphaHat.*d2fdx2(indIn,indIn);
SigmaTheta = VBA_inv(iSigmaTheta);

% mode
tmp = iQ*dtheta0(indIn) + alphaHat.*ddxdtheta(indIn);
deltaMuTheta = SigmaTheta*tmp;

% variational energy
Itheta = -0.5.*dtheta0(indIn)'*iQ*dtheta0(indIn) -0.5*alphaHat.*dx2;
if VBA_isWeird ({Itheta, SigmaTheta}) || div
    Itheta = -Inf;
end

% update sufficient statistics
suffStat.Itheta = Itheta;
suffStat.dx = dx;
suffStat.dx2 = dx2;
suffStat.dtheta = dtheta0;
suffStat.div = div;

author: MBB-team

repo: VBA-toolbox


text:
function test_ft_sourcegrandaverage

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_sourcegrandaverage ft_selectdata ft_selectdata_new ft_datatype_source

%% format 1

source = [];
source.dim = [10 11 12];
source.transform = eye(4);
source.avg.pow = rand(10*11*12,1);
source.inside = 1:660;
source.outside = 661:1320;

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);

%% format 2

source = [];
source.transform = eye(4);
source.pos = rand(1320,3);
source.pow = rand(1320,1);
source.inside = 1:660;
source.outside = 661:1320;

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);

%% format 3

source = [];
source.dim = [10 11 12];
source.transform = eye(4);
source.pow = rand(10*11*12,1);
source.inside = 1:660;
source.outside = 661:1320;

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);

%% format 4

source = [];
source.pos = rand(1320,3);
source.time = 1:25;
source.avg.pow = rand(10*11*12,25);
source.inside = 1:660;
source.outside = 661:1320;

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'time'), 'time field is missing');

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'time'), 'time field is missing');

%% format 5

source = [];
source.pos = rand(1320,3);
source.freq = 1:6;
source.time = 1:5;
source.avg.pow = rand(10*11*12,6,5);
source.inside = 1:660;
source.outside = 661:1320;

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'freq'), 'freq field is missing');
assert(isfield(grandavg, 'time'), 'time field is missing');

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'freq'), 'freq field is missing');
assert(isfield(grandavg, 'time'), 'time field is missing');

%% format 6

source = [];
source.pos = rand(1320,3);
source.time = 1:5;
source.inside = 1:660;
source.outside = 661:1320;
source.avg.pow = cell(1320,1);
for i=source.inside
  source.avg.pow{i} = randn(3,5);
end
for i=source.outside
  source.avg.pow{i} = [];
end

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'time'), 'time field is missing');

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'time'), 'time field is missing');

%% format 7

source = [];
source.pos = rand(1320,3);
source.freq = 1:6;
source.time = 1:5;
source.inside = 1:660;
source.outside = 661:1320;
source.avg.pow = cell(1320,1);
for i=source.inside
  source.avg.pow{i} = randn(3,6,5);
end
for i=source.outside
  source.avg.pow{i} = [];
end

ft_checkdata(source, 'datatype', 'source');

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'freq'), 'freq field is missing');
assert(isfield(grandavg, 'time'), 'time field is missing');

cfg.keepindividual = 'yes';
grandavg = ft_sourcegrandaverage(cfg, source, source);
assert(isfield(grandavg, 'freq'), 'freq field is missing');
assert(isfield(grandavg, 'time'), 'time field is missing');


text:
function h = p25_h ( n, x )

%*****************************************************************************80
%
%% P25_H evaluates the Hessian for problem 25.
%
%  Discussion:
%
%    Note that, for P = 0, the Hessian matrix should be diagonal.
%    However, if it is estimated using finite differences, off diagonal
%    terms may appear.  This occurs when the argument increment dX is
%    too small to be significant in terms such as X**6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  for i = 1 : n
    h(i,i) = i * 12.0 * x(i)^2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function el = export_gsps_module_field(args)
%"export_gsps_module_field"
%   Given a single planC.structures struct and a tag in the ROI_contour
%   module, return a properly populated and formatted instance of
%   that tag.
%
%   For speed, tag must be a decimal representation of the 8 digit
%   hexidecimal DICOM tag desired, ie instead of '00100010', pass
%   hex2dec('00100010');
%
%   Arguments are passed in a structure, arg:
%       arg.tag         = decimal tag of field to fill
%       arg.data        = CERR structure(s) to fill from
%       arg.template    = an empty template of the module created by the
%                         function build_module_template.m
%
%   This function requires arg.data is a planC.structures.
%
%JRA 06/19/06
%NAV 07/19/16 updated to dcm4che3
%
%Usage:
%   dcmobj = export_roi_contour_module_field(args)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
gspsS       = args.data{1};
scanS       = args.data{2};
template    = args.template;

switch tag
    %Class 1 Tags -- Required, must have data.
    case  7340162  %0070,0082 Presentation Creation Date
        data = datestr(now,'yyyymmdd');
        el = data2dcmElement(data, tag);

    case  7340163   % 0070,0083 Presentation Creation Time
        data = datestr(now,'hhmmss');
        el = data2dcmElement(data, tag);

    case 542113824  % 2050,0020  Presentation LUT Shape
        data = 'IDENTITY'; % or 'INVERSE'
        el = data2dcmElement(data, tag);

    case 7340122  % 0070,005A  Displayed Area Selection Sequence
        templateEl  = template.getValue(tag);
        fHandle = @export_displayed_area_selection_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        % Only one item since all images use the displayed area selection
        % Might need to add additional items when generating annotations
        % belonging to different series.

        % for i=1:length(gspsS)
        i = 1; % only one item
        dcmobj = export_sequence(fHandle, templateEl, {gspsS,scanS});
        el.add(i-1, dcmobj);
        % end
        %get attribute to return
        el = el.getParent();

    case 528661   % 0008,1115  Referenced Series Sequence
        templateEl  = template.getValue(tag);
        fHandle = @export_referenced_series_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        % Only one item since all images use the displayed area selection
        % Might need to add additional items when generating annotations
        % belonging to different series.

        % for i=1:length(gspsS)
        i = 1; % only one item
        dcmobj = export_sequence(fHandle, templateEl, {gspsS,scanS});
        el.add(i-1, dcmobj);
        % end
        %get attribute to return
        el = el.getParent();

    case 7340033  %0070,0001 Graphic Annotation Sequence

        templateEl  = template.getValue(tag);
        fHandle = @export_graphic_annotation_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        scanSopInstanceUIDc = {scanS.scanInfo.sopInstanceUID};

        for i=1:length(gspsS)
            slcNum = strncmp(gspsS(i).referenced_SOP_instance_uid,...
                scanSopInstanceUIDc,length(gspsS(i).referenced_SOP_instance_uid));
            dcmobj = export_sequence(fHandle, templateEl, {gspsS(i),scanS.scanInfo(slcNum)});
            %dcmobj = export_sequence(fHandle, tag, {structS(i), i});
            el.add(i-1, dcmobj);
        end
        %get attribute to return
        el = el.getParent();

    case 7340128   %0070,0060 Graphic Layer sequence
        templateEl  = template.getValue(tag);
        fHandle = @export_graphic_layer_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        [uniqLabelC,indUniqV] = unique({gspsS.presentLabel});
        numLabels = length(uniqLabelC);

        for i = 1:numLabels
            presentColor = gspsS(indUniqV(i)).presentRGBColor;
            dcmobj = export_sequence(fHandle, templateEl, {uniqLabelC{i},i,presentColor});
            el.add(i-1, dcmobj);
        end

        %get attribute to return
        el = el.getParent();

    %Class 2 Tags -- Must be present, can be NULL.
    %Class 3 Tags -- presence is optional, currently undefined.
    %Class 1C Tags -- presence is required under special circumstances
    %Class 2C Tags -- presence is required under special circumstances

    otherwise
        warning(['No methods exist to populate DICOM GSPS module field ' dec2hex(tag,8) '.']);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [Wout,Phiout,Phioutbis]=dualwavelets(W,Phi,LP)
	%for jr=1:Jres
	%filters.psi=W;
	%filters.phi=Phi;
	%LP=littlewood_paley_ISCV(filters);
	J=size(W,2);
	L=size(W{end},2);
	eps=1e-12;
	tight=0;

	Wout = W;
	Phiout = Phi;

	for j=1:J
		if ~isempty(W{j})
			for l=1:L
				if tight
					% if frame is tight, LP = 1 for all omega, so we only need to conjugate
					Wout{j}{l}=conj(W{j}{l});
				else
					% if frame is not tight, we need to divide by LP so that
					% sum W_j*Wd_j = sum W_j*conj(W_j)/LP = LP/LP = 1
					% but this can be dangerous is LP becomes very small, so only
					% do this on the support of W_j
					supp=(abs(W{j}{l})>eps);
					Wout{j}{l} = supp.*(conj(W{j}{l})./(LP));
				end
			end
		end
	end
	eps = 2e-2;
	supp=(abs(Phi)>eps);
	Phiout=supp.*((Phi)./(LP));
	if tight
		Phiout=Phi;
		Phioutbis=supp.*min(8,1./Phi);
		Ic=find(supp==0);
		Phioutbis(Ic)=0;
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%PHASECORRELATE  Detect translational shifts that occur between two images
%
%     pshift = cv.phaseCorrelate(src1, src2)
%     [pshift,response] = cv.phaseCorrelate(src1, src2)
%     [...] = cv.phaseCorrelate(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __src1__ First source floating-point array (single-channel `single` or
%   `double`).
% * __src2__ Second source floating-point array (single-channel `single` or
%   `double`), of same size and type as `src1`.
%
% ## Output
% * __pshift__ detected phase shift (sub-pixel) between the two arrays `[x,y]`
% * __response__ Signal power within the 5x5 centroid around the peak, between
%   0 and 1 (optional).
%
% ## Options
% * __Window__ Floating-point array with windowing coefficients to reduce edge
%   effects (optional). Not set by default.
%
% The function is used to detect translational shifts that occur between two
% images.
%
% The operation takes advantage of the Fourier shift theorem for detecting the
% translational shift in the frequency domain. It can be used for fast image
% registration as well as motion estimation. For more information please see
% [Phase correlation](https://en.wikipedia.org/wiki/Phase_correlation).
%
% Calculates the cross-power spectrum of two supplied source arrays. The
% arrays are padded if needed with cv.getOptimalDFTSize.
%
% The function performs the following equations:
%
% * First it applies a
%   [Hanning window](https://en.wikipedia.org/wiki/Hann_function) to each image
%   to remove possible edge effects. This window is cached until the array
%   size changes to speed up processing time.
%
% * Next it computes the forward DFTs of each source array:
%
%       G_a = F{src1}, G_b = F{src2}
%
%   where `F` is the forward DFT.
%
% * It then computes the cross-power spectrum of each frequency domain array:
%
%       R = G_a * G_b^(*) / |G_a * G_b^(*)|
%
% * Next the cross-correlation is converted back into the time domain via the
%   inverse DFT:
%
%       r = F^(-1){R}
%
% * Finally, it computes the peak location and computes a 5x5 weighted
%   centroid around the peak to achieve sub-pixel accuracy.
%
%       (\Delta{x}, \Delta{y}) = weightedCentroid{argmax_(x,y){r}}
%
% * If non-zero, the response parameter is computed as the sum of the elements
%   of `r` within the 5x5 centroid around the peak location. It is normalized
%   to a maximum of 1 (meaning there is a single peak) and will be smaller
%   when there are multiple peaks.
%
% See also: cv.dft, cv.getOptimalDFTSize, cv.idft, cv.mulSpectrums,
%  cv.createHanningWindow, imregcorr
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function y=bcrecur(a, n)
%
% Computation of the fractional difference coefficients
% by the recurrence relation
%           bc(j)=(1-(a+1)/j)*bc(j-1),
%  a - order of the fractional difference
%  n - required number of coefficients
%
%  Computation of bcrecur(k,k) takes (3*k+1) flops.
%
%  Copyright (C) Igor Podlubny
%  15 Nov 1994
%
%  See also:
%  [1] Podlubny, I.: Fractional Differential Equations.
%      Academic Press, San Diego, 1999, 368 pages, ISBN 0125588402.

y=cumprod([1, 1 - ((a+1) ./ (1:n))]);


text:
function e = imEntropy(img, varargin)
%IMENTROPY Compute entropy of an image
%
%   H = imEntropy(IMG)
%   Computes the entropy of the image IMG.
%   The entropy is computed as:
%     H  = -sum(P .* log2(P));
%   where P is the density probability that image takes a given value.
%
%   Note that the computation of the entropy depends on the way the
%   histogram is computed. This should not be a concern for uint8 images,
%   but can give results inconsistent with imMutualInformation for double
%   images.
%
%   H = imEntropy(IMG, N)
%   H = imEntropy(IMG, BINS)
%   Computes the histogram using N bins, or the bins specified by BINS.
%
%
%   Example
%     % compute entropy on a sample image
%     img = imread('rice.png');
%     H = imEntropy(img)
%     H =
%         7.0115
%
%     % entropy is independent of pixel ordering
%     img2 = circshift(img, [30 40]);
%     H = imEntropy(img2)
%     H =
%         7.0115
%
%     % entropy computed on an histogram with fewer bins
%     imEntropy(img, 16)
%     ans =
%         3.1158
%
%   See also
%     imJointEntropy, imMutualInformation, imHistogram
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-08-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% first compute histogram
h = imHistogram(img, varargin{:});

% keep only positive values, and normalize by 1 to have density
h(h==0) = [];
h = h / sum(h);

% compute entropy
e = -sum(h .* log2(h));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function uvMat = sc_uvMat_from_uvMap(map, uvPixInd)

% Function: conver uvMap into matrix form

% Extract uv transform from the NNF
[imgH, imgW, nCh] = size(map);

offset = uint32((0:nCh-1)*imgH*imgW);
uvPixInd = bsxfun(@plus, uvPixInd, offset);
uvMat = map(uvPixInd);

% map(uvPixInd) = data;
%
%
%
% uvMat = zeros(nCh, numUvPixels, 'single');
%
% for i = 1: nCh
%     uvMat(i,:) = map(uvPix.ind + (i-1)*imgH*imgW);
% end

% Indices of valid uv pixels
% uvValid.ind = uvMapMask(uvPix.ind);
% uvValid.pos = find(uvValid.ind);

end

repo: StructCompletion


text:
function varargout = main(varargin)
% Main: callback functions for main window
%
%    ABOUT
%
%      -Created:     January 2004
%      -Last update: February 13th, 2004
%      -Revision:    0.0.1
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

% Edit the above text to modify the response to help_button register

% Last Modified by GUIDE v2.5 24-Mar-2004 17:39:06

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @register_OpeningFcn, ...
    'gui_OutputFcn',  @register_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin & isstr(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 register is made visible.
function register_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 register (see VARARGIN)

% Choose default command line output for register
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

if strcmp(get(hObject,'Visible'),'off')
    initialize_gui(hObject, handles);
end

% --- Outputs from this function are returned to the command line.
function varargout = register_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;

function initialize_gui(fig_handle, handles)

setappdata(0, 'test_data', '1');

function open_child_Callback(hObject, eventdata, handles)

child;

function show_params_Callback(hObject, eventdata, handles)

msgbox(getappdata(0, 'test_data'), 'Parameter Shown');


text:
function y = eq(X,Y)
%EQ (overloaded)

if isa(X,'blkvar')
    X = sdpvar(X);
end

if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

try
    y = constraint(X,'==',Y);
catch
    error(lasterr)
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function snakedemo2
% SNAKEDEMO - Example on how to build a Snake Toy, and twist
%             it into predefined shapes.
%

% Copyright (C) 2005-2011 The MathWorks Inc.

%%
  function b = oneblock(dark, parent)
  % Create one snake block.
  % DARK is boolean, true if this is a dark block.
  % PARENT is the object that is the parent of the new block.

    % The transform object is going to be our BLOCK representation.
    % The patches in the transform are the visible aspect of that block.
    b = hgtransform('parent',parent);

    % Create the patches.  Coordinate map.
    %
    %     SIDE         TOP/BACK
    %
    %  -1     1      -1     1
    % -1X-----X      X------X -1
    %   |    /       |      |
    %   |   /        |      |
    %   |  /         |      |
    %   | /          |      |
    %   |/           |      |
    % 1 X            X------X 1
    %
    %                  B (1 -1 -1)
    %               __+
    %(-1 -1 -1) __--   --__
    %  A    __--           --__
    %     +-                 __+  F  (1 -1 1)
    %     |--__          __-- /
    %     |    --__  __--    /
    %     |        +-       /
    %     |        |  E    /
    %     |        |      /  (-1 -1 1)
    %     |        |     /
    %     |        |    /
    %     |        |   /
    %     +        |  /
    %   C  --__    | /
    %          --__|/
    %(-1 1 -1)     +
    %                D  (-1 1 1)

    v = [ -1 -1 -1 ;
          1 -1 -1 ;
          -1 1 -1 ;
          -1 1 1 ;
          -1 -1 1 ;
          1 -1 1 ;
        ];

    A = 1;
    B = 2;
    C = 3;
    D = 4;
    E = 5;
    F = 6;

    f = [ A B C nan
          E F D nan
          A B F E
          A E D C
          C D F B
        ];

    %
    % Specify the face color for each block in the snake
    % as either a cream color, or brown.
    %
    if dark
      fc = [ .4 .4 .2 ];
    else
      fc = [ 1 1 .9 ];
    end

    patch('parent',b, ...
          'vertices', v, ...
          'faces', f, ...
          'facecolor',fc, ...
          'edgecolor','k');

    % Transformation function that moves each block
    % into the correct location within the parent.
    xform(b,0);

    % Remember the current rotation value for animation.
    setappdata(b,'rotation',0);
  end

%%
  function xform(block, rotation)
  % Apply a transform to BLOCK.
  % ROTATION is the amount of rotation off the default for that blcck.

    % Create the compound transform.
    tform = makehgtform('axisrotate',[0 1 0],rotation,...
                        'zrotate',pi/2,'yrotate',pi,...
                        'translate',[2 0 0]);

    % Apply the transformation to BLOCK.
    % Remember that BLOCK is an hgtransform object, and not
    % A patch object.  The patch just goes along for the ride.
    set(block,'matrix',tform);
  end

  %%
  % Create the Figure window and Axesin which we are going to create the
  % snake.

  figure('renderer','opengl');
  title('Snake Demo');
  % Force the axes to have sane aspect ratios.
  axis vis3d
  % Make the axes visible.  Allow objects to be drawn outside
  % the axes clip rectangle.
  set(gca,'vis','off','clip','off');
  set(gca,'pos',[0 0 1 1]);
  % Start with an off-angle view.
  view(3)

  % Create the first block.
  snake(1) = oneblock(false, gca);

  % We need 24 blocks to make a traditional snake toy.
  for i = 2:24
    snake(i) = oneblock(~mod(i,2),snake(i-1));
  end

  %%
  % We intend to animate snake transformations.
  % To do this so it is reliable on multiple platforms,
  % running at different speeds, we need to measure the time
  % it takes to draw, and then add a pause so the animation
  % is not too fast.
  function docam

    % Orbit the camera so there is lots of movement.
    % Good orbit values are usuall two numbers, where
    % the second value is slightly smaller than the first.
    camorbit(5,3);

    % Measure how long it takes to draw the snake.
    % We could use drawnow('expose'), but it is nice to
    % let someone use camera toolbar and menu items
    % while the snake is animating.
    tic
      drawnow;
    h=toc;

    % Calculate the delta between the desired pause of .04 seconds.
    delta = .04-h;

    % Add more delay
    if delta > 0
      pause(delta);
    end
  end

  %%
  % To animate the snake, we need a routine to permute the snake
  % to new angles.
  function config(in, conf)

    % Loop over all the blocks in order.
    for k=1:length(in)

      % Fetch the old rotation value from the input block.
      old = getappdata(in(k),'rotation');

      % Create an animation step value for 10 total steps.
      step = (conf(k)-old)/10;

      % Skip blocks already configured
      if step ~= 0

        % Perform the animation.
        for j = 1:10
          % Use the same transform function as we had before.
          xform(in(k),old+(step*j));
          docam;
        end
      end

      % Solve any round-off problems in the animation.
      xform(in(k),conf(k));
      % Store the new value into the block
      setappdata(in(k),'rotation',conf(k));
    end
    % Float the camera around a little bit so you can
    % Get the feel for what's going on.
    for k=1:50
      docam;
    end

  end

  %%
  % Perform some animations.
  %
  % By defining pre-defined rotation vectors, we can animate
  % the snake through different permutations.

  % Convenient name for 90 degrees.
  n=pi/2;

  %% Ball
  rv = [ 0 n -n -n n -n n n -n n -n -n ...
         n -n n n -n n -n -n n -n n n];
  config(snake, rv);
  %% Wave
  rv = [ 0 0 pi 0 pi 0 pi 0 pi 0 pi 0 ...
         pi 0 pi 0 pi 0 pi 0 pi 0 pi 0];
  config(snake, rv);
  %% Dna
  rv = [ n n n n n n n n n n n n ...
         n n n n n n n n n n n n];
  config(snake, rv);
  %% Spiral
  rv = [ 0 n 0 n 0 n 0 n 0 n 0 n ...
         0 n 0 n 0 n 0 n 0 n 0 n];
  config(snake, rv);
  %% Coil
  rv = [ 0 n pi n pi n pi n pi n pi n ...
         pi n pi n pi n pi n pi n pi n];
  config(snake, rv);
  %% box
  rv = [ 0 0 0 0 0 0 pi 0 0 0 0 ...
         pi 0 0 0 0 0 0 pi 0 0 0 0 pi];
  config(snake, rv);
  %% filledbox
  rv = [ 0 pi 0 0 0 pi pi 0 0 0 0 pi ...
         0 0 pi 0 0 0 0 pi pi 0 0 0];
  config(snake, rv);
  %% cross
  rv = [ 0 0 0 pi pi 0 0 0 pi 0 pi pi ...
         0 pi 0 pi pi 0 pi 0 pi pi 0 pi];
  config(snake, rv);
  %% heart
  rv = [ 0 0 pi 0 0 n 0 n 0 0 0 0 ...
         pi 0 0 0 0 -n 0 -n 0 0 pi 0];
  config(snake, rv);
  %% dog
  rv = [ 0 pi 0 0 0 pi pi 0 pi 0 0 ...
         pi 0 pi pi 0 0 0 pi 0 pi pi 0 pi ];
  config(snake, rv);
  %% chicken
  rv = [ 0 pi 0 0 0 0 0 pi pi 0 ...
         n 0 pi pi 0 pi pi 0 -n ...
         0 pi pi 0 0];
  config(snake, rv);
  %% swan
  rv = [ 0 pi pi 0 0 0 0 0 pi -n -n -n ...
         n -n n n -n n -n -n -n pi 0 pi];
  config(snake, rv);
  %% duck
  rv = [ 0 pi pi 0 0 0 n 0 n 0 0 -n ...
         pi n 0 0 n pi -n 0 0 n 0 n];
  config(snake, rv);
  %% robot
  rv = [ 0 0 n 0 -n n n n -n -n -n -n ...
         pi n n n n -n -n -n n 0 -n 0];
  config(snake, rv);
  %% symbol
  rv = [ 0 pi 0 pi n n -n n -n -n pi 0 ...
         pi pi 0 pi n n -n n -n -n pi 0];
  config(snake, rv);
  %% fluer
  rv = [ 0 0 n -n n 0 pi 0 -n n -n 0 ...
         pi 0 n -n n 0 pi 0 -n n -n 0];
  config(snake, rv);
  %% bowl
  rv = [ 0 0 pi pi 0 n -n pi n n -n -n ...
         pi n n -n -n pi n n -n -n pi n];
  config(snake, rv);
  %% spring
  rv = [ 0 -n n 0 n -n n 0 n -n n 0 ...
         n -n n 0 n -n n 0 n -n n 0];
  config(snake, rv);
  %% cobra
  rv = [ 0 pi -n 0 0 0 pi n n 0 pi 0 0 ...
         0 n 0 pi pi 0 pi pi 0 n 0];
  config(snake, rv);
  %% line
  rv = [ 0 0 0 0 0 0 0 0 0 0 0 0 ...
         0 0 0 0 0 0 0 0 0 0 0 0];
  config(snake, rv);

end


text:
function dicom_tarzip(varargin)
% tars all .dcm files in subdirectories you specify, and deletes the
% original .dcm files if successful.
%
% :Inputs:
%
%   none (default): all directories under the current one
%   directory wildcard, e.g., 'run*' to search only those subdirectories.
%   absolute path names should also be ok.
%
% :Examples:
% ::
%
%    dicom_tarzip
%
%    dicom_tarzip('18*')
%
% Or, let's say you have runs within subjects, and dicom files within run dirs.
% Subject dirs all start with NSF*.  run dirs all start with 18*
% Starting from the directory above subjects, you could do this to zip ALL subjects:
% ::
%
%    subjdirs = dir('NSF*');
%    for s = 1:length(subjdirs)
%        dicom_tarzip([subjdirs(s).name filesep '18*']);
%    end
%
% Needless to say, use this with extreme caution, as it deletes your original files.
% And test it thoroughly with a backup dataset first!!
%
% ..
%    Note: Ash's checksum program would be very useful to add here.
% ..

dirwildcard = '*';

if ~isempty(varargin), dirwildcard = varargin{1}; end

prevwd = pwd;

% list directories

d = dir(dirwildcard);

isd = cat(1, d.isdir);
d = d(isd);
d(strmatch('.', char(d.name))) = [];  % exclude . dirs

if isempty(d)
    fprintf('No matching directories for %s. Quitting.\n', dirwildcard);
    cd(prevwd);
    return
else
    disp('Working on directories:')
    disp(char(d.name))
    disp('  ')
end

basedir = fullfile(pwd, fileparts(dirwildcard));
% % dir returns relative paths, so go there if we have a path name for a
% % wildcard
% if ~isdir(basedir)
%     fprintf('%s is not a directory. Quitting.\n\n', basedir);
%     return
% end
if ~isempty(basedir), cd(basedir); end  % this would be the subject directory if a path name was entered as a wildcard
fprintf('Working in dir: %s\n', basedir)

for i = 1:length(d)

    dicomtarname = [d(i).name filesep d(i).name '_dicoms.tar'];

    fprintf('Working on: %s.\n', dicomtarname);

    if exist(dicomtarname, 'file')
        fprintf('%s already exists!!! Skipping.\n\n', dicomtarname);
        continue
    end

    str = ['!tar cvf ' dicomtarname ' ' d(i).name filesep '*dcm >& tarchecktmp.txt  > tarfilelisttmp.txt'];
    disp(str)
    eval(str)
    errchk = fileread('tarchecktmp.txt');
    isok = isempty(strfind(errchk, 'Error'));

    if ~isok
        fprintf('Error in creating: %s .\n', dicomtarname);
        disp(errchk)
        disp('Removing failed tar archive.')
        str = ['!rm ' dicomtarname];
        disp(str)
        eval(str)

    else
        % remove original files
        fprintf('%s created successfully. Checking number of files...\n', dicomtarname);

        % check length/files
        % TOO LONG str = ['!tar -t ' dicomtarname ' > tarchecktmp.txt'];
        %filelist = fileread('tarfilelisttmp.txt');
        fid = fopen('tarfilelisttmp.txt');
        filelist = textscan(fid, '%s');
        fclose(fid);
        filelist = filelist{1};
        nfiles = length(filelist);

        dd = dir([d(i).name filesep '*dcm']);
        nfiles2 = length(dd);

        if nfiles == nfiles2
            fprintf('tar archive contains %3.0f files.\n', nfiles);
            fprintf('Removing originals.\n');
            str = ['!rm ' d(i).name filesep '*dcm'];
            disp(str)
            eval(str)
            disp('   ')
        else
            fprintf('tar archive contains %3.0f files, but counted %3.0f that should be there.\n', nfiles, nfiles2);
            disp('Check this. Not removing originals.');
        end

    end

    !rm tarchecktmp.txt
    !rm tarfilelisttmp.txt

    disp('   ')

end

if ~isempty(basedir), cd(prevwd); end

end % main function

%dcmf = filenames([d(i).name filesep '*dcm'], 'char');

%%
% dcmf2 = [];
% if ~isempty(dcmf)
%     for ff = 1:size(dcmf, 1)
%        dcmf2 = strcat(dcmf2, [dcmf(ff, :) ' ']);
%     end
%
% end
% str = ['!tar cvf ' d(i).name '_dicoms.tar ' dcmf2];

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function att = coarsealign3(imudata,lat,smoothepoch)
% att = coarsealign3(imudata,lat,smoothepoch);
% initial coarse alignment based on the acceleration and gyro measurement
% given the initial latitude using Dr. Jekeli's book
% Author: Yudan Yi
% Aug. 2004
% Feb. 2005
%--------------------------------------------------------------------------
if (nargin<2) error('error in data'); end;
if (isempty(imudata)|isempty(lat)) error('error in data'); end;
[n,m] = size(imudata);
if (n<1 || m<7) error('error in data'); end;
if (isempty(lat)) error('error in data'); end;
if (nargin<3) smoothepoch = min(256,n); end;
if (isempty(smoothepoch) || smoothepoch<1) smoothepoch = min(134,n); end;
%--------------------------------------------------------------------------
we = 7.292115147e-05;  % Earth Rotate rate
att = [];
att_index = 0;
smoothepoch=fix(smoothepoch);
for index=1:smoothepoch:n
 	x=index:min(index+smoothepoch-1,n);
    curacc = sum(imudata(x,2:4))/size(imudata,1);
    curgyro = sum(imudata(x,5:7))/size(imudata,1);
    curaccgyro = cross(-curacc,curgyro);
%     gravmag = norm(curaccgyro);
    gravmag = norm(curacc);
    matrix1=[ 0 0 -gravmag; we*cos(lat) 0 -we*sin(lat); 0 gravmag*we*cos(lat) 0];
    matrix2=[curacc;curgyro;curaccgyro];
    c_nb = matrix1\matrix2;% see Jekeli, 2000, page 243, (8.12)
    curatt = Cbn2att(c_nb);

    att_index = att_index+1;
    att(att_index,:)= [imudata(index,1) curatt(:)'];
end
% [att1, pivot] = extract_pivot_angle(att(:,2:4)');
% meanatt = (mean(att1')'+pivot(:))*180/pi;
% stdatt = (std(att1'))*180/pi;
% fprintf('Mean:\t%15.12f\t%15.12f\t%15.12f\n',meanatt);
% fprintf('Std:\t%15.12f\t%15.12f\t%15.12f\n',stdatt);
%--------------------------------------------------------------------------
% figure
% plot(att(:,2:4)*180/pi,'-','Marker','.');
% legend(['\mu=' num2str(meanatt(1),'%5.2f') '\circ \sigma=' num2str(stdatt(1),'%5.2f') '\circ r'],...
%        ['\mu=' num2str(meanatt(2),'%5.2f') '\circ \sigma=' num2str(stdatt(2),'%5.2f') '\circ p'],...
%        ['\mu=' num2str(meanatt(3),'%5.2f') '\circ \sigma=' num2str(stdatt(3),'%5.2f') '\circ h']);
% grid
% xlabel('time');
% ylabel('\circ');
% title('initial coarse alignment for attitude (h,p,r): deg(\circ)');
% % roll, pitch, heading => heading, pitch and roll
% att(:,2) = meanatt(3)*pi/180; % heading
% att(:,3) = meanatt(2)*pi/180; % pitch
% att(:,4) = meanatt(1)*pi/180; % roll

repo: ekfmonoslam


text:
function [mps2] = ftps22mps2(ftps2)
% Convert acceleration from feet per square-second to meters per second
% squared.
% Chad A. Greene 2012
mps2 = ftps2*0.3048;
end


text:
function [ Y ] = classify2( Model, X )

    % Extract Features from the test data
    data = getImages(X);

    weights1 = Model.weights1;
    weights2 = Model.weights2;
    projection = Model.projection;

    % Project the test data on the new dimension
    xTest = getFeatures(data, projection);

    % Classify the images
    [Y, ~] = NNClassify(weights1, weights2, xTest);
end

function [Features] =  getImages(data)
    Features = [];
    for i = 1:size(data,1)
        image = reshape(data(i,:),[32,32,3]);
        image = imresize(image,4);
        feat = extract_feature(image);
        Features = horzcat(Features,feat);
    end
    Features = Features';
end

function [Feats] = getFeatures(Features, projection)
    featureSize = size(Features,2);

    % Perform Mean Normalization on the Feature Matrix.
    % Each pixel can have intensity values from 0 - 255
    for i = 1 : featureSize
        Features(:,i) = Features(:,i) - mean(Features(:,i));
    end

    % Project the Features on the new dimension
    Feats = projection * Features';
    Feats = Feats';
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [y]=fun(v,x)
y0=v(1);
A=v(2);
w=v(3);
x0=v(4);
y=y0+(2*A/pi).*(w./(4*(x-x0).^2+w.^2));


text:
function kernel = calibSPIRiT(kCalib, kSize, nCoils, CalibTyk)

% kernel = calibSPIRiT(kCalib, kSize, nCoils, CalibTyk)
%
% Function calibrates a SPIRiT kernel from a calibration area in k-space
%
%
% (c) Michael Lustig 2013
%

[AtA] = dat2AtA(kCalib,kSize);
for n=1:nCoils
	kernel(:,:,:,n) = calibrate(AtA,kSize,nCoils,n,CalibTyk);
end

author: thomaskuestner


text:
function pitchwatch(x,Ts)
% Plot the pitch keys.
% pitchwatch(x,[Ts])
%
% :: Syntax
%    The array x is the input signal and Ts is the (optional) sampling period.
%    Example on use: [x,Fs] = wavread('Hum.wav');
%                    pitchwatch(x,1/Fs);
%
% :: Information
%    Make your own wav-files with the Windows Sound Recorder. Choose the attributes
%    PCM 8000Hz, 16bit, Mono when saving the wav-file. For more information on pitch
%    keys, go to http://www.bookrags.com/wiki/Piano_key_frequencies.

% Set parameters.
if nargin < 2, Ts = 1/8000; end

% Set constants.
xlen = length(x);
wlen = 512;
wlag = 128;

% Zero-pad signal.
xpad = wlag-rem(xlen-wlen-1,wlag)-1;
x    = [x(:); zeros(xpad,1)];
L    = (xlen+xpad-wlen)/wlag+1;
L    = L-wlen/wlag/2;

% Calculate the AMDF and AMDF-fractional pitch periods.
for k1 = 1:L
   k2 = (k1-1)*wlag;
   for k3 = 1:wlen/2-1, c(k3) = sum(abs(x(k2+(1:wlen))-x(k2+k3+(1:wlen)))); end

   n = findpeaks(-c);
   if ~isempty(n), n(find(c(n) > mean(c(n))-sqrt(var(c(n))))) = []; end
   if ~isempty(n), t0 = n(1);
   else,           t0 = []; end

   if ~isempty(t0)
      u1    = x(k2+t0+(1:wlen))-x(k2+(1:wlen));
      u2    = x(k2+t0+(1:wlen))-x(k2+t0+1+(1:wlen));
      t1    = sum(u1.*u2)/sum(u2.*u2);
      y(k1) = 1/Ts/(t0+t1);
   else
      y(k1) = NaN;
   end
end

% Plot pitch keys.
keys = {'C','C#','D','D#','E','F','F#','G','G#','A','A#','B'};
key0 = 12*log2(y/440)+57;
u1   = max(0,min(floor(key0)));
u2   = min(96,max(ceil(key0)));
key1 = mod([u1:u2],12)+1;
key2 = floor([u1:u2]/12);
for k = [1:u2-u1+1], tick(k) = strcat(keys(key1(k)),num2str(key2(k))); end
figure, plot((0.5+[0:L-1])*wlag*Ts,key0,'b.');
set(gca,'FontSize',8,'XLim',[0,L*wlag]*Ts,'YLim',[u1,u2],'YTick',[u1:u2],'YTickLabel',tick);
xlabel('Time'); ylabel('Key');

% FUNCTIONS

function n = findpeaks(x)

n    = find(diff(diff(x) > 0) < 0);
u    = find(x(n+1) > x(n));
n(u) = n(u)+1;


text:
%MDS_STRESS - Sammon stress between dissimilarity matrices
%
% 	E = MDS_STRESS(Q,DS,D)
%
% INPUT
% 	Q					Indicator of the Sammon stress; Q = -2,-1,0,1,2
% 	DS				Original distance matrix
% 	D 				Approximated distance matrix
%
% OUTPUT
% 	E 				Sammon stress
%
% DESCRIPTION
% Computes the Sammon stress between the original distance matrix Ds
% and the approximated distance matrix D, expressed as follows:
%
%  E = 1/(sum_{i<j} DS_{ij}^(q+2)) sum_{i<j} (DS_{ij} - D_{ij})^2 * DS_{ij}^q
%

%
% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [e,alpha] = mds_stress (q,Ds,D,isratio)

	if nargin < 4
		isratio = 0;
	end

	[m,k]   = size(Ds);
	if any(size(D) ~= size(Ds)),
		error ('The sizes of matrices do not match.');
	end
	mk = m*k;

	D  = +D;
	Ds = +Ds;

	% I is  the index of non-zero (> eps) values to be included
	% for the computation of the stress

	I = 1:mk;
	nanindex = find(isnan(Ds(:)) | isnan(D(:)));
	if ~isempty(nanindex),
		I(nanindex) = [];
	end
	O = [];
	if m == k & (length(intersect(find(D(:) < eps), 1:m+1:(mk))) == m),
		O  = 1:m+1:mk;
		Ds(O) = 1;
		D (O) = 1;
    mm = m - 1;
	else
		mm = k;
	end

  if isratio,
    II = setdiff(I,O);
    alpha = sum((Ds(II).^q).*D(II).^2)/sum((Ds(II).^(q+1)).*D(II));
    Ds = alpha*Ds;
	else
		alpha = 1;
	end

	c = sum(Ds(I).^(q+2)) - length(O);
	if q ~= 0,
		e = sum(Ds(I).^q .* ((Ds(I)-D(I)).^2))/c;
	else
		e = sum(((Ds(I)-D(I)).^2))/c;
	end
return;

author: marianux

repo: ecg-kit

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


text:
% MAKEELEC - subroutine to make electrode file in EEGPLOT
%
% Usage:   >> makeelec(chans)
%          >> [channames] = makeelec(chans)
%
% Inputs:
%   chans - number of channels
%
% Author: Colin Humprhies, CNL / Salk Institute, 1996
%
% See also: EEGPLOT

% Copyright (C) Colin Humphries, CNL / Salk Institute, Aug, 1996
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

% 01-25-02 reformated help & license, added links -ad

function channames = makeelec(chans);

global is_done_pressed
is_done_pressed = 0;

lines = ceil(chans/2);

fh = 180 + chans*12;
fl = 400;

figure('Units','Pixels','Position',[400 200 fl fh])

fighandle = gcf;

axes('Units','Pixels','Position',[0 0 fl fh],'Visible','off','Xlim',[0 fl],'Ylim',[0 fh]);

text(fl/2,fh-20,'Make Electrode-File','HorizontalAlignment','Center','Fontsize',14)

titleb = uicontrol('Units','Pixels','Style','Edit','Position',...
                     [fl/2 fh-60 .4*fl 20],'UserData',0,'HorizontalAlignment','left');
text(fl/2-20,fh-50,'Filename:','HorizontalAlignment','Right')

text(fl/2,fh-75,'Electrodes','HorizontalAlignment','center')
%line([fl*.2 fl*.8],[fh-85 fh-85],'color','w')

for i = 1:chans

   k = 2*ceil(i/2);

   if (i==k)
      ede(i,:) = uicontrol('Units','Pixels','Style','Edit',...
'Position',[fl*.65 fh-90-k*12 fl*.15 20],'UserData',i,'HorizontalAlignment','left');
      text(fl*.65-15,fh-90-k*12+10,num2str(i),'HorizontalAlignment','right')
   else
      ede(i,:) = uicontrol('Units','Pixels','Style','Edit',...
'Position',[fl*.2 fh-90-k*12 fl*.15 20],'UserData',i,'HorizontalAlignment','left');
      text(fl*.2-15,fh-90-k*12+10,num2str(i),'HorizontalAlignment','right')
   end

end

TIMESTRING = ['chans1973 = get(gco,''UserData'');','OH1973 = findobj(''Style'',''edit'',''UserData'',0);','FILENAME1973 = get(OH1973,''string'');','if isempty(FILENAME1973);','fprintf(''Filename Missing'');','else;','channames1973 = zeros(chans1973,6);','for i1973 = 1:chans1973;','Aa1973 = findobj(''Style'',''edit'',''UserData'',i1973);','elabel1973 = get(Aa1973,''string'');','channames1973(i1973,6-length(elabel1973)+1:6) = elabel1973;','end;','for i1973 = 1:length(channames1973(:));','if channames1973(i1973) == 0;','channames1973(i1973) = ''.'';','end;','end;','FID1973 = fopen(FILENAME1973,''w'');','fprintf(FID1973,''%s'',channames1973'');','fclose(FID1973);','fprintf(''Saving file'');','end;','clear chans1973 OH1973 FILENAME1973  Aa1973 elabel1973 ii1973 FID1973 channames1973 i1973;'];
saveb = uicontrol('Units','Pixels','Style','PushButton','Position',[fl/10 20 fl*.15 25],'String','Save','UserData',chans,'Callback',TIMESTRING);

TIMESTRING = ['global is_done_pressed;','is_done_pressed = 1;','clear is_done_pressed'];
closeb = uicontrol('Units','Pixels','Style','PushButton','Position',...
               [fl*.45 20 fl*.15 25],'String','Done','Callback',TIMESTRING);

TIMESTRING = ['global is_done_pressed;','is_done_pressed = 2;','clear is_done_pressed'];
cancelb = uicontrol('Units','Pixels','Style','PushButton','Position',...
               [fl*.75 20 fl*.15 25],'String','Cancel','Callback',TIMESTRING);

while(is_done_pressed == 0)
   drawnow;
end

if is_done_pressed == 1
   channames = [];
   for i = 1:chans
      a = findobj('Style','edit','UserData',i);
      elabel = get(a,'string');
      channames = str2mat(channames,elabel);
   end
   channames = str2mat(channames,' ');
else

   channames = [];
end

delete(fighandle)
clear is_done_pressed

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function cvx_optval = sum_square_pos( x, varargin ) %#ok

%SUM_SQUARE_POS   Internal cvx version.

narginchk(1,2);
x2 = [];
cvx_begin
    variable x2( size( x ) );
    minimize( sum_square( x2, varargin{:} ) );
    x2 >= x; %#ok
cvx_end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% FUNCTION [idx,pnts] = kdtree_range(kdtree, range)
%
% AUTHOR:     Steven Michael
%             (smichael@ll.mit.edu)
%
% DATE:       2/17/05
%
% DESCRIPTION:
%
%  This function simply returns all the points of a kdtree that
%  are within an N-dimesional volume bounded by "range"
%
% INPUTS:
%
%   kdtree :    A KD Tree previously created with kdtree_create
%
%
%   range  :    An array (ndim X 2) of points that specify the bounds
%               of an "ndim" dimensional volume.  All the points of the
%               kdtree within this volume will be returned.
%               The range is boundary inclusive: [min,max]
%
%               Alternatively, range can be a (nsearches X ndim X 2)
%               array of points that describes multiple range boxes
%               -- the "nsearches" index references the range box.  The
%               output in this case will be put into "nsearches"
%               different cell arrays.
%
%
% OUTPUTS:
%
%   idx    :    A (1 X N) array of points, where N is the number of
%               points within the "ndim" dimensional volume. The value is
%               the index to the releavant point in the array from which
%               the tree was created.
%
%               If multiple range boxes are searches (a 3D array is input
%               for range), then "pnt" will be a cell array, where the
%               "nth" cell contains a (1XN) array of points ; N is the
%               number of points within the "ndim" dimensional volume
%               described by the "nth" index of the range input.
%
%
%   pnts:  :    A (N X ndim) array the actual points. If the actual point
%               instead of the index is desired, an optional second output
%               will hold the points.  Multiple range boxes are handled
%               as described above.
%
% Example:
%
%    % Create a list of 1000 random points in 3d space
%    r = rand(1000,3);
%
%    % Create a tree from this list
%    tree = kdtree(r);
%
%    % Find the point closest to the origin
%    [pntidx,pntval] = kdtree_closestpoint(tree,[0 0 0]);
%
%    % Create a list "r2" of 100 random points in 3d space and
%    % find the points in "r" that are closest to each point in "r2"
%    [pntidx,pntval] = kdtree_closestpoint(tree,r2);
%
%    % Find all the points within the cube defined by "rng"
%    rng = [ [.45 .55]; [.45 .55]; ; [.45 .55] ];
%    pntidx = kdtree_range(tree,rng);


text:
function [u,eqn,info] = fracLapP2P2(node,elem,pde,option)
%% FRACLAPP2P2 solves fractional Laplacian equation using P2-P2 element
%
% [u,eqn,info] = fracLapP1P1(node,elem,pde,option) solves the fractional
% Laplacian equation
%
%  (-\Delta^^s u = f in \Omega with u = 0 on \partial \Omega.

%% Options
if ~exist('option','var'), option = []; end
if ~isfield(option,'plotflag'), option.plotflag = 0; end

%% Parameters
global s;
alpha = 1-2*s;
if s == 0.5
    gamma = 1;
else
    gamma = 5/(2*s)+0.01; % a different grading factor
end

%% Mesh and data structure
[elem2dof,edge,bdDof] = dofP2(elem);
N = size(node,1);  NT = size(elem,1); NE = size(edge,1);
My = round(pde.L*sqrt(NT/2));
y = gradmap(0,pde.L,gamma,My);
Ny = length(y); % number of vertices in y-direction
NTy = Ny - 1;   % number of elements in y-direction
NTtotal = NT*NTy;
Nxdof = (N + NE); % dof in x  is number of vertices plus edges
Nydof = (Ny + NTy); % dof in y direction is number of vertices plus elements
Ndof = Nxdof*Nydof;

tic;
%% Stiffness matrix and mass matrix in the extended direction
% quantities in y direction
hy = diff(y);
a = zeros(NTy,5);
for i = 1:5
    a(:,i) = diff(y.^(alpha+i)/(alpha+i));
end
y1 = y(1:end-1);
y2 = y(2:end);
y12 = y1 + y2;
y1y2 = y1.*y2;
ym = y12/2;
% stiffness matrix in y direction
Ay = zeros(NTy,3,3);
Ay(:,1,1) = a(:,1)./hy.^2;
Ay(:,1,2) = -Ay(:,1,1);
Ay(:,2,1) = -Ay(:,1,1);
Ay(:,2,2) = Ay(:,1,1);
Ay(:,1,3) = 8*(a(:,2) - ym.*a(:,1))./hy.^3;
Ay(:,3,1) = Ay(:,1,3);
Ay(:,2,3) = -Ay(:,1,3);
Ay(:,3,2) = Ay(:,2,3);
Ay(:,3,3) = 64*(a(:,3) - 2*ym.*a(:,2) + ym.^2.*a(:,1))./hy.^4;
% mass matrix in y direction
My = zeros(NTy,3,3);
My(:,1,1) = (a(:,3) - 2*y1.*a(:,2) + y1.^2.*a(:,1))./hy.^2;
My(:,1,2) = (-a(:,3) + (y1+y2).*a(:,2) - y1.*y2.*a(:,1))./hy.^2;
My(:,2,1) = My(:,1,2);
My(:,2,2) = (a(:,3) - 2*y2.*a(:,2) + y2.^2.*a(:,1))./hy.^2;
My(:,1,3) = 4*(a(:,4) - (y12+y2).*a(:,3) + (y2.*y12 + y1y2).*a(:,2) ...
             - y2.*y1y2.*a(:,1))./hy.^3;
My(:,3,1) = My(:,1,3);
My(:,2,3) = 4*(-a(:,4) + (y12+y1).*a(:,3) - (y1.*y12 + y1y2).*a(:,2) ...
             + y1.*y1y2.*a(:,1))./hy.^3;
My(:,3,2) = My(:,2,3);
My(:,3,3) = 16*(a(:,5) - 2*y12.*a(:,4) + (y12.^2 + 2*y1y2).*a(:,3)...
               -2*y12.*y1.*y2.*a(:,2) + y1y2.^2.*a(:,1))./hy.^4;

%% Stiffness matrix and mass matrix in the original direction
[Dlambda,area] = gradbasis(node,elem);
% Compute a piecewise constant diffusion coefficient
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                                 % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)       % d is a function
    [lambda,weight] = quadpts(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
        K = K + weight(p)*pde.d(pxy);
   end
end
if ~isempty(pde.d) % build the coefficients into the scaled area
    areaK = K.*area;
else
    areaK = area;
end
At = zeros(NT,6,6);
Mt = zeros(NT,6,6);
% stiffness matrix
[lambda, w] = quadpts(2);
nQuad = size(lambda,1);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    for i = 1:6
        for j = 1:6
            At(:,i,j) = At(:,i,j) + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*areaK;
        end
    end
end
% mass matrix
[lambda, w] = quadpts(4);
nQuad = size(lambda,1);
phi(:,6) = 4*lambda(:,1).*lambda(:,2);
phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
phi(:,4) = 4*lambda(:,2).*lambda(:,3);
phi(:,5) = 4*lambda(:,3).*lambda(:,1);
for i = 1:6
    for j = 1:6
        for p = 1:nQuad
            Mt(:,i,j) = Mt(:,i,j) + w(p)*phi(p,i).*phi(p,j);
        end
        Mt(:,i,j) = Mt(:,i,j).*areaK;
    end
end
clear phi Dphip lambda Dlambda

%% Assemble stiffness matrix
dofMap = repmat(1:Nxdof,Nydof,1) + repmat((0:Nydof-1)'*Nxdof,1,Nxdof);
ii = zeros(18*18*NTtotal,1);
jj = zeros(18*18*NTtotal,1);
sA = zeros(18*18*NTtotal,1);
index = 0;
for m = 1:3
    for n = 1:3
        for i = 1:6
            for j = 1:6
                if m < 3
                    ii(index+1:index+NTtotal) = dofMap(m:Ny+m-2,elem2dof(:,i));
                else % m = 3: middle points of elements in y direction
                    ii(index+1:index+NTtotal) = dofMap(Ny+1:Nydof,elem2dof(:,i));
                end
                if n < 3
                    jj(index+1:index+NTtotal) = dofMap(n:Ny+n-2,elem2dof(:,j));
                else % n = 3: middle points of elements in y direction
                    jj(index+1:index+NTtotal) = dofMap(Ny+1:Nydof,elem2dof(:,j));
                end
                sA(index+1:index+NTtotal) = kron(At(:,i,j),My(:,m,n)) + ...
                                            kron(Mt(:,i,j),Ay(:,m,n));
                index = index + NTtotal;
            end
        end
    end
end
A = sparse(ii,jj,sA,Ndof,Ndof);
clear ii jj sA At Mt Ay My

%% Assemble the right hand side
b = zeros(Ndof,1);
u = zeros(Ndof,1);

%% Set up boundary conditions
% bdDof is found in dofP2
lateralbd = dofMap(:,bdDof);
topbd = dofMap(Ny,:);

% Modify the matrix to include the Dirichlet boundary condition
bdidx = zeros(Ndof,1);
bdidx(lateralbd) = 1;
bdidx(topbd) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*A*T + Tbd;

% Compute boundary integral over the original domain
if ~isfield(option,'fquadorder')
    option.fquadorder = 6;
end
[lambdaf,weightf] = quadpts(option.fquadorder);
phif(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
phif(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
phif(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
phif(:,4) = 4*lambda(:,2).*lambda(:,3);
phif(:,5) = 4*lambda(:,3).*lambda(:,1);
phif(:,6) = 4*lambda(:,1).*lambda(:,2);
nQuadf = size(lambdaf,1);
ft = zeros(NT,6);
for p = 1:nQuadf
    % quadrature points in the x-y coordinate
    pxy = lambdaf(p,1)*node(elem(:,1),:) ...
        + lambdaf(p,2)*node(elem(:,2),:) ...
        + lambdaf(p,3)*node(elem(:,3),:);
    fp = pde.f(pxy);
    for i = 1:6
        ft(:,i) = ft(:,i) + weightf(p)*phif(p,i)*fp;
    end
end
ft = ds*ft.*repmat(area,1,6);
b = b + accumarray(elem2dof(:),ft(:),[Ndof,1]);
% Neumann edfts are considered as open set. So the corner points should be
% set as Dirichlet boundary condition!
b(lateralbd) = 0;
clear ft

%% Record assembling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
freeNode = find(bdidx==0);
if isempty(freeNode), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for larft size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeNode) = A(freeNode,freeNode)\b(freeNode);
        residual = norm(b - A*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        option.x0 = u;
        option.solver = 'VCYCLE';
        option.freeDof = freeDof;
        [u,info] = mgfracLapP2P2(A,b,elem,option);
    case 'amg'
        option.solver = 'CG';
        [u(freeNode),info] = amg(A(freeNode,freeNode),b(freeNode),option);
end

%% Compute error using boundary integral
if isfield(pde,'exactu')
    err = zeros(NT,1);
    [lambda,weight] = quadpts(7);
    phi(:,6) = 4*lambda(:,2).*lambda(:,1);
    phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
    phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
    phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
    phi(:,4) = 4*lambda(:,2).*lambda(:,3);
    phi(:,5) = 4*lambda(:,1).*lambda(:,3);
    nQuad = size(lambda,1);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
             + lambda(p,2)*node(elem(:,2),:) ...
             + lambda(p,3)*node(elem(:,3),:);
        fp = pde.f(pxy);
        uhp = u(elem2dof(:,1)).*phi(p,1) + ...
              u(elem2dof(:,2)).*phi(p,2) + ...
              u(elem2dof(:,3)).*phi(p,3) + ...
              u(elem2dof(:,4)).*phi(p,4) + ...
              u(elem2dof(:,5)).*phi(p,5) + ...
              u(elem2dof(:,6)).*phi(p,6);
        up = pde.exactu(pxy);
        err = err + weight(p)*fp.*(up - uhp);
    end
else
    err = 0;
end
err = ds*sum(err.*area);
err = sqrt(err);
info.errH1 = err;

%% Output information
eqn = struct('A',A,'b',b,'freeNode',freeNode,'edge',edge);
info.assembleTime = assembleTime;

repo: ifem


text:
function [df d_arr] = richardson_based_derivative(eq, x0, n, h)
% July 1, 2013 by Ehsan Behnam.
% [df d_arr] = richardson_based_derivative(...)
% computes the derative of the equation "eq" % at "x0"
% using "n" iterations of Richardson extrapolation.
% "h" specifies the initial differential increment/decrement around x0
% (i.e. f(x+h) or f(x-h)).
% A couple of notes:
% 1) "eq" can be any string in the form of f(x). Note that it should be
% computable in Matlab so for example, use eq = 'exp(x)' instead of 'e^x'.
% 2) No matter of what is "h" at each step h <- h/2. Therefore if you want
% to compare the accuracy of this method to other methods (e.g. Tyler's
% series), you need to adjust the number of required steps properly.

% The first output "df" is the d(eq(x))/dx |@(x = x0). Regarding the
% educational purposes for producing the whole table, the function also
% outputs the whole array iteratively generated as the function progresses.

% Basic usage example:
% eq = 'x^3*exp(-x)';
% df = richardson_based_derivative(eq, 0.1, 10);
% Finding the difference between "df" and the Matlab output ...
% disp(subs(diff(sym(eq)), 0.1) - df);
% It will show "-4.6263e-07" on the workspace.

if (nargin < 2 || nargin > 4)
    help('richardson_based_derivative');
    df = '';
    d_arr = '';
    return;
end

if (nargin == 2)
    h = 1/2;
    n = 5;
end

if (nargin == 3)
    h = 1/2;
end

fx = sym(eq);

%Initialization:
d_arr = zeros(n);
for i = 1:n
    d_arr(i, 1) = subs(fx, 'x', x0 + h) - subs(fx, 'x', x0 - h);
    d_arr(i, 1) = d_arr(i, 1) / (2 * h);
    % Iteratively computing the coefficients:
    for j = 2:i
        d_arr(i, j) = d_arr(i, j - 1) + (d_arr(i, j - 1) - ...
                      d_arr(i - 1, j - 1)) / (4 ^j -1);
    end
    % next iteration "h":
    h = h/2;
end
df = d_arr(end, end);

repo: mcbench-benchmarks