text:
classdef StructuredEdgeDetection < handle
    %STRUCTUREDEDGEDETECTION  Class implementing edge detection algorithm
    %
    % As described in [Dollar2013].
    %
    % ### Structured forests for fast edge detection
    %
    % This module contains implementations of modern structured edge detection
    % algorithms, i.e. algorithms which somehow takes into account pixel
    % affinities in natural images.
    %
    % ![image01](https://docs.opencv.org/3.3.1/01.jpg)
    % ![image02](https://docs.opencv.org/3.3.1/02.jpg)
    % ![image03](https://docs.opencv.org/3.3.1/03.jpg)
    % ![image04](https://docs.opencv.org/3.3.1/04.jpg)
    % ![image05](https://docs.opencv.org/3.3.1/05.jpg)
    % ![image06](https://docs.opencv.org/3.3.1/06.jpg)
    % ![image07](https://docs.opencv.org/3.3.1/07.jpg)
    % ![image08](https://docs.opencv.org/3.3.1/08.jpg)
    % ![image09](https://docs.opencv.org/3.3.1/09.jpg)
    % ![image10](https://docs.opencv.org/3.3.1/10.jpg)
    % ![image11](https://docs.opencv.org/3.3.1/11.jpg)
    % ![image12](https://docs.opencv.org/3.3.1/12.jpg)
    %
    % ## References
    % [Dollar2013]:
    % > Piotr Dollar and C Lawrence Zitnick. "Structured forests for fast edge
    % > detection". In IEEE International Conference on Computer Vision (ICCV)
    % > 2013, pages 1841-1848.
    %
    % [Lim2013]:
    % > Joseph J Lim, C Lawrence Zitnick, and Piotr DollAr. "Sketch tokens: A
    % > learned mid-level representation for contour and object detection".
    % > In IEEE Conference on Computer Vision and Pattern Recognition (CVPR),
    % > 2013, pages 3158-3165.
    %
    % See also: cv.StructuredEdgeDetection.StructuredEdgeDetection
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = StructuredEdgeDetection(model, varargin)
            %STRUCTUREDEDGEDETECTION  The only constructor
            %
            %     obj = cv.StructuredEdgeDetection(model)
            %     obj = cv.StructuredEdgeDetection(model, howToGetFeatures)
            %
            % ## Input
            % * __model__ name of the file where the model is stored.
            % * __howToGetFeatures__ optional name of MATLAB M-function that
            %   implements custom feature extractor. A helper function for
            %   training part of: "P. Dollar and C. L. Zitnick. Structured
            %   Forests for Fast Edge Detection, 2013". You need it only if
            %   you would like to train your own forest, otherwise leave it
            %   unspecified for the default implementation. See example below.
            %
            % ## Example
            % The following is an example of a custom feature extractor
            % MATLAB function:
            %
            %     % This function extracts feature channels from src. The
            %     % StructureEdgeDetection uses this feature space to detect
            %     % edges.
            %     function features = myRFFeatureGetter(src, opts)
            %         % src: source image to extract features
            %         % features: output n-channel floating-point feature matrix
            %         % opts: struct of options
            %         gnrmRad = opts.normRad;    % gradientNormalizationRadius
            %         gsmthRad = opts.grdSmooth; % gradientSmoothingRadius
            %         shrink = opts.shrink;      % shrinkNumber
            %         outNum = opts.nChns;       % numberOfOutputChannels
            %         gradNum = opts.nOrients;   % numberOfGradientOrientations
            %
            %         nsize = [size(src,1) size(src,2)] ./ shrink;
            %         features = zeros([nsize outNum], 'single');
            %         % ... here your feature extraction code
            %     end
            %
            % TODO: Custom extractor is not internally used in the current
            % cv.StructuredEdgeDetection implementation. See this
            % [tutorial](https://docs.opencv.org/3.3.1/d2/d59/tutorial_ximgproc_training.html)
            % for more information about training your own structured forest
            % (it uses an external MATLAB toolbox for the training part).
            %
            % See also: cv.StructuredEdgeDetection.detectEdges
            %
            this.id = StructuredEdgeDetection_(0, 'new', model, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.StructuredEdgeDetection
            %
            if isempty(this.id), return; end
            StructuredEdgeDetection_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.StructuredEdgeDetection.empty,
            %  cv.StructuredEdgeDetection.load
            %
            StructuredEdgeDetection_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.StructuredEdgeDetection.clear,
            %  cv.StructuredEdgeDetection.load
            %
            b = StructuredEdgeDetection_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.StructuredEdgeDetection.load
            %
            StructuredEdgeDetection_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.StructuredEdgeDetection.save
            %
            StructuredEdgeDetection_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.StructuredEdgeDetection.save,
            %  cv.StructuredEdgeDetection.load
            %
            name = StructuredEdgeDetection_(this.id, 'getDefaultName');
        end
    end

    %% StructuredEdgeDetection
    methods
        function dst = detectEdges(this, src)
            %DETECTEDGES  The function detects edges in src and draw them to dst
            %
            %     dst = obj.detectEdges(src)
            %
            % ## Input
            % * __src__ source image (RGB, float, in [0;1]) to detect edges.
            %
            % ## Output
            % * __dst__ destination image (grayscale, float, in [0;1]) where
            %   edges are drawn.
            %
            % The algorithm underlies this function is much more robust to
            % texture presence, than common approaches, e.g. cv.Sobel.
            %
            % See also: cv.Sobel, cv.Canny
            %
            dst = StructuredEdgeDetection_(this.id, 'detectEdges', src);
        end

        function orientation_image = computeOrientation(this, edge_image)
            %COMPUTEORIENTATION  Computes orientation map from edge image
            %
            %     orientation_image = obj.computeOrientation(edge_image)
            %
            % ## Input
            % * **edge_image** edge image from `detectEdges` function.
            %
            % ## Output
            % * **orientation_image** orientation image.
            %
            % See also: cv.StructuredEdgeDetection.detectEdges
            %
            orientation_image = StructuredEdgeDetection_(this.id, 'computeOrientation', edge_image);
        end

        function dst = edgesNms(this, edge_image, orientation_image, varargin)
            %EDGESNMS  Suppress edges (nonmaximum suppression)
            %
            %     dst = obj.edgesNms(edge_image, orientation_image)
            %     dst = obj.edgesNms(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * **edge_image** edge image from `detectEdges` function.
            % * **orientation_image** orientation image from `computeOrientation` function.
            %
            % ## Output
            % * __dst__ suppressed image (grayscale, float, in [0;1]).
            %
            % ## Options
            % * __R__ radius for NMS suppression. default 2
            % * __S__ radius for boundary suppression. default 0
            % * __M__ multiplier for conservative suppression. default 1.0
            % * __IsParallel__ enables/disables parallel computing.
            %   default true
            %
            % The function suppresses edges where edge is stronger in
            % orthogonal direction.
            %
            % See also: cv.StructuredEdgeDetection.detectEdges,
            %  cv.StructuredEdgeDetection.computeOrientation
            %
            dst = StructuredEdgeDetection_(this.id, 'edgesNms', edge_image, orientation_image, varargin{:});
        end
    end

    %% Static functions
    methods (Static)
        function features = getFeatures(src, opts)
            %GETFEATURES  Extracts features from image
            %
            %     features = cv.StructuredEdgeDetection.getFeatures(src, opts)
            %
            % ## Input
            % * __src__ source image to extract features (RGB float in [0;1]).
            % * __opts__ a scalar struct of random forest options
            %   (feature params), with the following fields:
            %   * __normRad__ `gradientNormalizationRadius` gradient
            %     normalization radius.
            %   * __grdSmooth__ `gradientSmoothingRadius` radius for smoothing
            %     of gradients (using convolution with triangle filter).
            %   * __shrink__ `shrinkNumber` amount to shrink channels.
            %   * __nChns__ `numberOfOutputChannels` number of edge
            %     orientation bins for output.
            %   * __nOrients__ `numberOfGradientOrientations` number of
            %     orientations per gradient scale.
            %
            % ## Output
            % * __features__ extracted features.
            %
            % Extracted features are appropriate for StructuredEdgeDetection
            % training.
            %
            % See also: cv.StructuredEdgeDetection.StructuredEdgeDetection
            %
            features = StructuredEdgeDetection_(0, 'getFeatures', src, opts);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function hB = plugin_board(board)
% Use Plugin API to create board plugin object

hB = hdlcoder.Board;

% Target Board Information
hB.BoardName    = sprintf('AnalogDevices ADRV9371 %s', upper(board));

% FPGA Device
hB.FPGAVendor   = 'Xilinx';

% Determine the device based on the board
switch(upper(board))
    case 'ZC706'
        hB.FPGAFamily   = 'Zynq';
        hB.FPGADevice   = sprintf('xc7%s', 'z045');
        hB.FPGAPackage  = 'ffg900';
        hB.FPGASpeed    = '-2';
    case 'ZCU102'
        hB.FPGAFamily   = 'Zynq UltraScale+';
        hB.FPGADevice   = sprintf('xc%s', 'zu9eg-ffvb1156-2-e');
        hB.FPGAPackage  = '';
        hB.FPGASpeed    = '';
    otherwise
        hB.FPGAFamily   = 'Zynq UltraScale+';
        hB.FPGADevice   = sprintf('xc%s', 'zu9eg-ffvb1156-2-e');
        hB.FPGAPackage  = '';
        hB.FPGASpeed    = '';
end

% Tool Info
hB.SupportedTool = {'Xilinx Vivado'};

% FPGA JTAG chain position
hB.JTAGChainPosition = 2;

%% Add interfaces
% Standard "External Port" interface

author: analogdevicesinc


text:
function [dat] = read_neuralynx_sdma(dataset, begsample, endsample, chanindx);

% READ_NEURALYNX_SDMA read specified channels and samples from a Neuralynx splitted DMA dataset
%
% Use as
%    [hdr] = read_neuralynx_sdma(dataset)
%    [dat] = read_neuralynx_sdma(dataset, begsample, endsample, chanindx)
%
% The splitted DMA dataset is not a formal Neuralynx format, but at
% the FCDC we use it in conjunction with SPIKEDOWNSAMPLE. The dataset
% directory contains files, one for each channel, each containing a
% 8-byte header followed by the binary values for all samples. Commonly
% the binary values are represented as int32, but it is possible to use
% int16 or other numeric representations. The 8-byte header specifies the
% numeric representation and the bitshift that should be applied (in case
% of integer representations).
%
% This function returns the electrophysiological data in AD units
% and not in uV. You should look up the details of the headstage and
% the Neuralynx amplifier and scale the values accordingly.

% Copyright (C) 2006-2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_neuralynx_sdma.m 945 2010-04-21 17:41:20Z roboos $

needhdr = (nargin==1);
needdat = (nargin>=2);

if ~isdir(dataset)
  error('dataset should be a directory');
end

% determine the content of the dataset directory
dirlist = dir(dataset);

% determine the correspondence between files and channels
label         = cell(length(dirlist),1);
filesel       = zeros(length(dirlist),1);
headerfile    = [];
[pd, nd, xd]  = fileparts(dataset);
for i=1:length(dirlist)
  [pf, nf, x1] = fileparts(dirlist(i).name);
  [pf, nf, x2] = fileparts(nf);
  % compare the name of the channel with the name of the dataset and look at the file extention
  % the file should be called "dataset/dataset.chanlabel" or "dataset/dataset.chanlabel.bin"
  if strcmp(nf, nd) && ~isempty(x1)
    if ~isempty(x2)
      label{i}   = x2(2:end);
      filesel(i) = 1;
    elseif strcmp(x1, '.txt')
      % this is not a proper channel, but an ascii file with header info
      label{i}   = [];
      filesel(i) = 0;
      headerfile = fullfile(dataset, dirlist(i).name);
    else
      label{i}   = x1(2:end);
      filesel(i) = 1;
    end
  else
    label{i}   = [];
    filesel(i) = 0;
  end
end
label = label(filesel==1);
bytes = cell2mat({dirlist(filesel==1).bytes}); % number of bytes in each file
filelist = {dirlist(filesel==1).name};         % filename excluding the path
clear dirlist filesel pd nd xd

for i=1:length(filelist)
  % include the full path in the filename
  filelist{i} = fullfile(dataset, filelist{i});
end

% there is a preferred ordering, which corresponds with the channel order in the non-splitted DMA file
preferred = {
  'stx'
  'pid'
  'siz'
  'tsh'
  'tsl'
  'cpu'
  'ttl'
  'x01'
  'x02'
  'x03'
  'x04'
  'x05'
  'x06'
  'x07'
  'x08'
  'x09'
  'x10'
  'csc001'
  'csc002'
  'csc003'
  'csc004'
  'csc005'
  'csc006'
  'csc007'
  'csc008'
  'csc009'
  'csc010'
  'csc011'
  'csc012'
  'csc013'
  'csc014'
  'csc015'
  'csc016'
  'csc017'
  'csc018'
  'csc019'
  'csc020'
  'csc021'
  'csc022'
  'csc023'
  'csc024'
  'csc025'
  'csc026'
  'csc027'
  'csc028'
  'csc029'
  'csc030'
  'csc031'
  'csc032'
  'csc033'
  'csc034'
  'csc035'
  'csc036'
  'csc037'
  'csc038'
  'csc039'
  'csc040'
  'csc041'
  'csc042'
  'csc043'
  'csc044'
  'csc045'
  'csc046'
  'csc047'
  'csc048'
  'csc049'
  'csc050'
  'csc051'
  'csc052'
  'csc053'
  'csc054'
  'csc055'
  'csc056'
  'csc057'
  'csc058'
  'csc059'
  'csc060'
  'csc061'
  'csc062'
  'csc063'
  'csc064'
  'csc065'
  'csc066'
  'csc067'
  'csc068'
  'csc069'
  'csc070'
  'csc071'
  'csc072'
  'csc073'
  'csc074'
  'csc075'
  'csc076'
  'csc077'
  'csc078'
  'csc079'
  'csc080'
  'csc081'
  'csc082'
  'csc083'
  'csc084'
  'csc085'
  'csc086'
  'csc087'
  'csc088'
  'csc089'
  'csc090'
  'csc091'
  'csc092'
  'csc093'
  'csc094'
  'csc095'
  'csc096'
  'csc097'
  'csc098'
  'csc099'
  'csc100'
  'csc101'
  'csc102'
  'csc103'
  'csc104'
  'csc105'
  'csc106'
  'csc107'
  'csc108'
  'csc109'
  'csc110'
  'csc111'
  'csc112'
  'csc113'
  'csc114'
  'csc115'
  'csc116'
  'csc117'
  'csc118'
  'csc119'
  'csc120'
  'csc121'
  'csc122'
  'csc123'
  'csc124'
  'csc125'
  'csc126'
  'csc127'
  'csc128'
  'csc129'
  'csc130'
  'csc131'
  'csc132'
  'csc133'
  'csc134'
  'csc135'
  'csc136'
  'csc137'
  'csc138'
  'csc139'
  'csc140'
  'csc141'
  'csc142'
  'csc143'
  'csc144'
  'csc145'
  'csc146'
  'csc147'
  'csc148'
  'csc149'
  'csc150'
  'csc151'
  'csc152'
  'csc153'
  'csc154'
  'csc155'
  'csc156'
  'csc157'
  'csc158'
  'csc159'
  'csc160'
  'csc161'
  'csc162'
  'csc163'
  'csc164'
  'csc165'
  'csc166'
  'csc167'
  'csc168'
  'csc169'
  'csc170'
  'csc171'
  'csc172'
  'csc173'
  'csc174'
  'csc175'
  'csc176'
  'csc177'
  'csc178'
  'csc179'
  'csc180'
  'csc181'
  'csc182'
  'csc183'
  'csc184'
  'csc185'
  'csc186'
  'csc187'
  'csc188'
  'csc189'
  'csc190'
  'csc191'
  'csc192'
  'csc193'
  'csc194'
  'csc195'
  'csc196'
  'csc197'
  'csc198'
  'csc199'
  'csc200'
  'csc201'
  'csc202'
  'csc203'
  'csc204'
  'csc205'
  'csc206'
  'csc207'
  'csc208'
  'csc209'
  'csc210'
  'csc211'
  'csc212'
  'csc213'
  'csc214'
  'csc215'
  'csc216'
  'csc217'
  'csc218'
  'csc219'
  'csc220'
  'csc221'
  'csc222'
  'csc223'
  'csc224'
  'csc225'
  'csc226'
  'csc227'
  'csc228'
  'csc229'
  'csc230'
  'csc231'
  'csc232'
  'csc233'
  'csc234'
  'csc235'
  'csc236'
  'csc237'
  'csc238'
  'csc239'
  'csc240'
  'csc241'
  'csc242'
  'csc243'
  'csc244'
  'csc245'
  'csc246'
  'csc247'
  'csc248'
  'csc249'
  'csc250'
  'csc251'
  'csc252'
  'csc253'
  'csc254'
  'csc255'
  'csc256'
  'crc'
  };

[sel1, sel2] = match_str(preferred, label);
if length(sel2)==length(label)
  % all reorder the labels/files
  label = label(sel2);
  bytes = bytes(sel2);
  filelist = filelist(sel2);
else
  % not all files in this SDMA dataset could be accounted for in the
  % preference list, hence do not attempt to apply the preferred ordering
end

if needhdr
  % construct a general header for the complete dataset
  if ~isempty(headerfile)
    orig             = neuralynx_getheader(headerfile);
    hdr.Fs           = orig.SamplingFrequency;
    hdr.nSamples     = [];        % see below
    hdr.nSamplesPre  = 0;         % number of pre-trigger samples in each trial
    hdr.nTrials      = 1;         % number of trials
    hdr.label        = label;
    hdr.nChans       = length(label);
    hdr.orig.dataset = orig;       % keep the header details
  else
    % some parts of the header have to be hardcoded, since the splitted dataset does not contain all header information
    hdr.Fs           = 32556;     % sampling frequency
    hdr.nSamples     = [];        % see below
    hdr.nSamplesPre  = 0;         % number of pre-trigger samples in each trial
    hdr.nTrials      = 1;         % number of trials
    hdr.label        = label;
    hdr.nChans       = length(label);
  end

  % read the header of each individual file, i.e. for each variable
  clear orig
  for i=1:length(filelist)
    orig(i) = read_neuralynx_bin(filelist{i});
  end

  % determine the number of samples for each channel
  nsamples = cell2mat({orig.nSamples});
  if any(nsamples~=nsamples(1))
    error('different number of samples over channels are not supported');
  else
    hdr.nSamples = nsamples(1);
  end

  % determine the sampling frequency for each channel
  fsample = cell2mat({orig.Fs});
  if any(diff(fsample))
    error('different sampling rates over channels are not supported');
  elseif any(fsample~=hdr.Fs)
    error('inconsistent sampling rates');
  end

  % determine the first and last timestamp, by reading them from the timestamp channels
  tslfile = filelist{find(strcmp('tsl', label))};
  tshfile = filelist{find(strcmp('tsh', label))};
  beg_tsl = read_neuralynx_bin(tslfile, 1, 1);
  beg_tsh = read_neuralynx_bin(tshfile, 1, 1);
  end_tsl = read_neuralynx_bin(tslfile, hdr.nSamples, hdr.nSamples);
  end_tsh = read_neuralynx_bin(tshfile, hdr.nSamples, hdr.nSamples);
  hdr.FirstTimeStamp = timestamp_neuralynx(beg_tsl, beg_tsh);
  hdr.LastTimeStamp  = timestamp_neuralynx(end_tsl, end_tsh);
  hdr.TimeStampPerSample = double(hdr.LastTimeStamp-hdr.FirstTimeStamp)./(hdr.nSamples-1);  % this should be double, since it can be fractional

  % also remember the original header details
  hdr.orig.chan = orig;

  % only return the header information
  dat = hdr;

else
  % allocate memory for all data
  dat = zeros(length(chanindx), endsample-begsample+1);

  % read all channels, one small chunk at at time, and write it to seperate files
  for i=1:length(chanindx)
    j = chanindx(i);
    dat(i,:) = double(read_neuralynx_bin(filelist{j}, begsample, endsample));
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function[]=makefigs_dlines
%MAKEFIGS_DLINES  Makes a sample figure for DLINES.

figure
subplot(2,2,1)
axis([-1 0 0 1]),dlines(-1,'b'),dlines(1,'r')
subplot(2,2,2),
axis([0 1 0 1]),dlines(-1,'b'),dlines(1,'r')
subplot(2,2,3)
axis([-1 0 -1 0]),dlines(-1,'b'),dlines(1,'r')
subplot(2,2,4)
axis([0 1 -1 0]),dlines(-1,'b'),dlines(1,'r')

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% ACWMF algorithm implementation
%function [y, noise_matrix] = ACWMF(I);
function output = ACWMF(I);

% Threshold parameters
%delta = [40 25 10 5];   %% for uniform impulse noise
delta = [55,40,25,15];  %% for salt and pepper noise
%delta = (delta1+delta2)/2;   %% for mixed impulse noise

x = double(I);
image_size = size(I);
B = im2col(padarray(x,[1 1],'symmetric','both'),[3 3],'sliding');

% Compute filter output
m = medfilt2(x,[3 3],'symmetric');

%Compute differences
d = abs(x-m);
d0 = d(:)';
clear d;
B(10,:) = x(:)';
B(11,:) = x(:)';
m1 = median(B);
d1 = abs(x(:)'-m1);
B(12,:) = x(:)';
B(13,:) = x(:)';
m1=median(B);
d2 = abs(x(:)'-m1);
B(14,:) = x(:)';
B(15,:) = x(:)';
m1=median(B);
d3 = abs(x(:)'-m1);
clear B;

% Compute MAD values
B_x = im2col(padarray(x,[1 1],'symmetric','both'),[3 3],'sliding');

for i = 1:9
    B(i,:) = abs(B_x(i,:) - m(:)');
end
MAD = median(B);

clear B;

% Compute threshold values
s = 0.1;
T1 = (MAD * s) + delta(1);
T2 = (MAD * s) + delta(2);
T3 = (MAD * s) + delta(3);
T4 = (MAD * s) + delta(4);

x2 = x(:);

% Detect noisy pixels
F = find((d0>T1)|(d1>T2)|(d2>T3)|(d3>T4));

%noise_matrix = zeros(image_size);
%noise_matrix(F) = 1;

% Replace noisy pixels
x2(F) = m(F);
clear F m d0 d1 d2 d3 T x;
output = uint8(col2im(x2,[1 1],image_size,'sliding'));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% NOMAD Install for OPTI Toolbox

% This file will help you compile NOMAD for use with MATLAB.

% My build platform:
% - Windows 8 x64
% - Visual Studio 2012

% To recompile you will need to get / do the following:

% 1) Get NOMAD
% NOMAD is available from http://www.gerad.ca/NOMAD/PHP_Forms/Download.php.
% Complete the download form then download the latest version. Once you
% have installed NOMAD, locate the /src/ directory.

% 2) Compile NOMAD
% The easiest way to compile NOMAD is to use the Visual Studio Project
% Builder included with OPTI. Use the following commands, substituting the
% required path on your computer:

% %% Visual Studio Builder Commands
% path = 'full path to NOMAD here'; %e.g. 'C:\Solvers\NOMAD'
% sdir = [path '\src'];
% name = 'libnomad';
% opts.exPP = {'_CRT_SECURE_NO_WARNINGS'};
% opts.exclude = 'nomad.cpp';
% VS_WriteProj(sdir,name,[],opts)
% %%

% Once complete, you will have a directory called NOMAD\libnomad. Open the
% Visual Studio 2012 project file, then complete the following steps:
%   a) In defines.hpp change the following:
%       - line 194 change MAX_DIMENSION to 5000.
%   b) Build a Win32 or x64 Release to compile the code.
%   c) Copy the generated .lib file to the following folder:
%
%   OPTI/Solvers/nomad/Source/lib/win32 or win64
%
%   And rename it to libnomad.lib.
%
%   You will also need to copy all header files from NOMAD/src to
%   the following folder:
%
%   OPTI/Solvers/nomad/Source/Include

% 3) NOMAD MEX Interface
% % The NOMAD MEX Interface is a simple MEX interface I wrote to use NOMAD.
% Note the MEX interface contains two versions, a GERAD version (for the
% authors of NOMAD) and an OPTI version (for easier integration into this
% toolbox). By default the GERAD version is built, however to build the
% OPTI version simply add the preprocessor "OPTI_VERSION" to the compiler.

% 4) Compile the MEX File
% The code below will automatically include all required libraries and
% directories to build the NOMAD MEX file. Once you have completed all the
% above steps, simply run this file to compile NOMAD! You MUST BE in the
% base directory of OPTI!

clear nomad

% Get Arch Dependent Library Path
libdir = opti_GetLibPath();

fprintf('\n------------------------------------------------\n');
fprintf('NOMAD MEX FILE INSTALL\n\n');

%Get NOMAD Libraries
post = [' -IInclude -L' libdir ' -llibnomad -llibut -output nomad -DOPTI_VERSION'];

%CD to Source Directory
cdir = cd;
cd 'Solvers/nomad/Source';

%Compile & Move
pre = 'mex -v -largeArrayDims nomadmex.cpp';
try
    eval([pre post])
    movefile(['nomad.' mexext],'../','f')
    fprintf('Done!\n');
catch ME
    cd(cdir);
    error('opti:nomad','Error Compiling NOMAD!\n%s',ME.message);
end
cd(cdir);
fprintf('------------------------------------------------\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function value = p26_autonomous ( )

%*****************************************************************************80
%
%% P26_AUTONOMOUS reports whether problem P26 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = conex1_inverse ( alpha )

%*****************************************************************************80
%
%% CONEX1_INVERSE returns the inverse of the CONEX1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar defining A.
%    A common value is 100.0.
%
%    Output, real A(4,4), the matrix.
%
  a(1,1) =  1.0;
  a(1,2) =  1.0 - alpha;
  a(1,3) =        alpha;
  a(1,4) =  2.0;

  a(2,1) =  0.0;
  a(2,2) =  1.0 + alpha;
  a(2,3) =      - alpha;
  a(2,4) =  0.0;

  a(3,1) =  0.0;
  a(3,2) = -1.0;
  a(3,3) =  1.0;
  a(3,4) =  1.0 / alpha;

  a(4,1) = 0.0;
  a(4,2) = 0.0;
  a(4,3) = 0.0;
  a(4,4) = 1.0 / alpha;

  return
end


text:
function flag = checkinterfacemesh3(node, elem, interfaceData)
%% CHECKINTERFACEMESH3 check the validity of the polyhedra mesh generated
%   by interfacemesh3.
%
%   Here we check  validity of the polyhedra mesh by the Euler formla:
%            F - E + V = 2
%   where F, E and V are the number of faces, edges and vertices of a simply
%   closed polyhedron.
%
%   Author: Huayi Wei <weihuayi@xtu.edu.cn>, based on discussion with Long Chen.
%
%   See also: interfacemesh3
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

face = interfaceData.face;
face2elem = interfaceData.face2elem;
N = size(node,1);
NF = size(face,1);

NP = max(face2elem);
isCutPoly = false(NP, 1);
isCutPoly(face2elem) = true;
cutPolyIdx = zeros(NP, 1);
NP = sum(isCutPoly);
cutPolyIdx(isCutPoly) = 1:NP;
face2elem = cutPolyIdx(face2elem);

isTriFace = face(:,4) == 0;
edge = 4*ones(NF,1);
edge(isTriFace) = 3;
E = accumarray(face2elem, edge)/2;
F = accumarray(face2elem, ones(NF,1));

poly2node = sparse(face2elem(isTriFace)*ones(1,3), face(isTriFace, 1:3), 1, NP,N)...
+ sparse(face2elem(~isTriFace)*ones(1,4), face(~isTriFace, 1:4), 1, NP,N);
poly2node = poly2node > 0;
NV = poly2node*ones(N,1);

isBdPoly = ((F - E + NV) ~= 2);

if ~isempty(find(isBdPoly,1))
    figure
    showmesh(node, face(isTriFace & isBdPoly(face2elem),1:3));
    sface = face(~isTriFace & isBdPoly(face2elem),1:4);
    hold on
    tsface = [sface(:, 1:3);sface(:,[3, 4,1])];
    showmesh(node, tsface, 'Facecolor','y');

%     bdPolyIdx = find(isBdPoly);
%     figure
%     for i = 1:size(bdPolyIdx)
%         sface = face(~isTriFace & face2elem == bdPolyIdx(i),1:4);
%         tsface = [sface(:, 1:3);sface(:,[3, 4,1])];
%         showmesh(node, [face(isTriFace & face2elem == bdPolyIdx(i),1:3);tsface]);
%     end
    flag = false;
else
    flag = true;
end

repo: ifem


text:
% Isoparametric Formulation Implementation
% clear memory
clear all; %close all; clc;

%% Physical parameters
% E: modulus of elasticity
% A: area of cross section
% L: length of bar
E=2e11; A=12.5e-4; L=[0.75,0.75];

%% Build Elements
% numberElements: number of elements
numberElements=2;
% numberNodes: number of nodes
numberNodes=3;
% generation of coordinates and connectivities
elementNodes=[1 2;2 3];
nodeCoordinates=[0 0.75 1.5];
% for structure:
   % displacements: displacement vector
   % force : force vector
   % stiffness: stiffness matrix
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
% computation of the system stiffness matrix and force vector
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:) ;
  detJacobian=L(e)/2;
  invJacobian=1/detJacobian;
  ngp = 2;
  [w,xi]=gauss1d(ngp);
  xc=0.5*(nodeCoordinates(elementDof(1))+nodeCoordinates(elementDof(2)));
  for ip=1:ngp;
      [shape,naturalDerivatives]=shapeFunctionL2(xi(ip));
      B=naturalDerivatives*invJacobian;
      stiffness(elementDof,elementDof)=...
      stiffness(elementDof,elementDof)+ B'*B*w(ip)*detJacobian*E*A;
      force(elementDof) = force(elementDof)+...
          60000*shape'*detJacobian*w(ip);
  end
end

%% BCs and solution
% prescribed dofs
prescribedDof = [1];
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness, ...
    numberNodes,prescribedDof,force)

%% Stress
stress = zeros(numberElements,1);
for e=1:numberElements;
    stress(elementNodes(e,:)) = E*B*displacements(elementNodes(e,:));
end

%% Exact solutions
x = 0:0.1:1.5;
displacements_exact=(90000*x-30000*x.^2)/(E*A);
stress_exact=-((60000*x)-90000)/A;

%% plot figures
figure(1)

% Displacements
subplot(1,2,1); hold on;
plot(nodeCoordinates,displacements,'-.sb')
plot(x,displacements_exact,'-r'); hold off;
xlabel('x'); ylabel('displacement, u');
legend('FEM','Exact Solution',2);

% Stress
subplot(1,2,2); hold on;
stairs(nodeCoordinates,stress,'-.sb')
plot(x,stress_exact,'-r'); hold off;
xlabel('x'); ylabel('stress, \sigma');
legend('FEM','Exact Solution',1);

repo: Aero-matlab


text:
classdef(Abstract) AbstractSensorState < matlab.mixin.SetGet & matlab.mixin.Heterogeneous & matlab.mixin.Copyable
    %AbstractSensorState Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        sensor = getSensor(obj);

        activeTf = getSensorActiveState(obj);

        setSensorActiveState(obj, activeTf);

        steeringModel = getSensorSteeringMode(obj);

        setSteeringModel(obj, steeringModel);

        range = getSensorMaxRange(obj);

        setSensorMaxRange(obj, range);

        angle = getMaxAngle(obj);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [] = test11()
% function test11()
%
% Try out centroid function in Nelder-Mead solver.
%

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    % Pick the manifold
    n = 10;
    k = 2;
%     M = euclideanfactory(n);
    M = obliquefactory(n, k);

    % Generate a bunch of points
    m = 25;
    x = cell(m, 1);
    for i = 1 : m
        x{i} = M.rand();
    end

    y = centroid(M, x);

    % This is for Euclidean only
    if strfind(M.name(), 'Euclidean')
        xx = zeros(n, 1);
        for i = 1 : m
            xx = xx + x{i}/m;
        end
        fprintf('Distance between mean point and centroid found: %e.\n',...
                 norm(xx-y));
    end

end

author: NicolasBoumal

repo: manopt


text:
function [CP, ppf] = getControlPoints( input_A, input_B, maxppf, ransac)
    disp('Detecting and Matching SIFT features...');
    fileListA = dir(input_A);
    fileListA = fileListA(3:length(fileListA));
    fileListB = dir(input_B);
    fileListB = fileListB(3:length(fileListB));
    nFrames = min(length(fileListA), length(fileListB));
    CP = zeros(nFrames, maxppf, 4);
    ppf = zeros(nFrames, 1);
    trackerA = vision.PointTracker('MaxBidirectionalError', 1);
    trackerB = vision.PointTracker('MaxBidirectionalError', 1);
    for frameIndex = 1:nFrames
        fprintf('%5d', frameIndex);
        if mod(frameIndex, 20) == 0
            fprintf('\n') ;
        end

        fileNameA = fileListA(frameIndex).name;
        fileNameB = fileListB(frameIndex).name;
        IA = imread([input_A fileNameA]);
        IB = imread([input_B fileNameB]);

        [H, W, ~] = size(IA);

        if frameIndex > 1
            setPoints(trackerA, trackA);
            [trackAcont, validityA] = step(trackerA, IA);
            setPoints(trackerB, trackB);
            [trackBcont, validityB] = step(trackerB, IB);
            trackAcont = trackAcont(validityA & validityB, :);
            trackBcont = trackBcont(validityA & validityB, :);
        end
        if ransac
            [trackAsurf, trackBsurf] = SURF(IA, IB);
        else
            [trackAsurf, trackBsurf] = SURF2(IA, IB);
        end

        trackA = trackAsurf;
        trackB = trackBsurf;

        valid = trackA(:, 1) > 0 & trackA(:, 1) < W & trackA(:, 2) > 0 & trackA(:, 2) < H ...
            & trackB(:, 1) > 0 & trackB(:, 1) < W & trackB(:, 2) > 0 & trackB(:, 2) < H;

        trackA = trackA(valid, :);
        trackB = trackB(valid, :);

        valid = filtermask(IA, trackA);
        trackA = trackA(valid, :);
        trackB = trackB(valid, :);
        valid = filtermask(IB, trackB);
        trackA = trackA(valid, :);
        trackB = trackB(valid, :);

        if frameIndex == 1
            initialize(trackerA, trackA, IA);
            initialize(trackerB, trackB, IB);
        else
            if length(trackA) + length(trackAcont) > maxppf
                ordering = randperm(length(trackAcont));
                trackAcont = trackAcont(ordering(1:maxppf - length(trackA)), :);
                trackBcont = trackBcont(ordering(1:maxppf - length(trackA)), :);
            end
            trackA = [trackA ; trackAcont];
            trackB = [trackB ; trackBcont];
        end

        IA = insertMarker(IA, trackA, 'o', 'color', 'red');
        IB = insertMarker(IB, trackB, 's', 'color', 'yellow');
        figure(1);
        imshow(IA);
        figure(2);
        imshow(IB);

        if length(trackA) > maxppf
            ppf(frameIndex) = maxppf;
            CP(frameIndex, :, 1:2) = trackA(1:maxppf, :);
            CP(frameIndex, :, 3:4) = trackB(1:maxppf, :);
        else
            ppf(frameIndex) = length(trackA);
            CP(frameIndex, 1:ppf(frameIndex), 1:2) = trackA;
            CP(frameIndex, 1:ppf(frameIndex), 3:4) = trackB;
        end
%         CP(frameIndex, :, :) = [featuresA featuresB; zeros(maxppf - ppf(frameIndex), 4)];

    end
end

function valid = filtermask(frame, points_)
    [H, W, ~] = size(frame);
    valid = points_(:, 1) > 0 & points_(:, 1) < W & points_(:, 2) > 0 & points_(:, 2) < H;
    points_ = points_(valid, :);
    mask = frame(:, :, 1) < 20 & frame(:, :, 2) < 20 & frame(:, :, 3) < 20;
    mask = imgaussfilt(double(mask), 50);
    videoH = size(frame, 1);
    mask(mask > 0.2) = 1;
    mask(mask ~= 1) = 0;
    valid = mask(round(points_(:, 1) - 1) * videoH + round(points_(:, 2))) == 0;
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function [dcos] = triang2(pp,tt)
%TRIANG2 calc. enclosed angles for a 2-simplex triangulation
%embedded in the two-dimensional plane.
%   [ADEG] = TRIANG2(VERT,TRIA) returns the enclosed angles
%   associated with each triangle, where ADEG is a T-by-3
%   array of the angles subtended at each vertex, VERT is a
%   V-by-2 array of XY coordinates, and TRIA is a T-by-3 ar-
%   ray of vertex indexing, where each row defines a triang-
%   le, such that VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and
%   VERT(TRIA(II,3),:) are the coordinates of the II-TH tri-
%   angle. Angles are returned in degrees.
%
%   See also TRISCR2, TRIAREA, TRIBAL2

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 08/07/2018

%---------------------------------------------- basic checks
    if (~isnumeric(pp) || ~isnumeric(tt) )
        error('triang2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(pp) ~= +2 || ndims(tt) ~= +2 )
        error('triang2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(pp,2)~= +2 || size(tt,2) < +3 )
        error('triang2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nnod = size(pp,1) ;

%---------------------------------------------- basic checks
    if (min(min(tt(:,1:3))) < +1 || ...
            max(max(tt(:,1:3))) > nnod )
        error('triang2:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%----------------------------------- compute enclosed angles
    dcos = zeros(size(tt,1),3) ;

    ev12 = pp(tt(:,2),:)-pp(tt(:,1),:) ;
    ev23 = pp(tt(:,3),:)-pp(tt(:,2),:) ;
    ev31 = pp(tt(:,1),:)-pp(tt(:,3),:) ;

    lv11 = sqrt(sum(ev12.^2,2));
    lv22 = sqrt(sum(ev23.^2,2));
    lv33 = sqrt(sum(ev31.^2,2));

    ev12 = ev12 ./ ...
        lv11(:,ones(1,size(pp,2)));
    ev23 = ev23 ./ ...
        lv22(:,ones(1,size(pp,2)));
    ev31 = ev31 ./ ...
        lv33(:,ones(1,size(pp,2)));

    dcos(:,1) = sum(-ev12.*ev23,2);
    dcos(:,2) = sum(-ev23.*ev31,2);
    dcos(:,3) = sum(-ev31.*ev12,2);

    dcos(:,1) = max(-1.,dcos(:,1));
    dcos(:,1) = min(+1.,dcos(:,1));
    dcos(:,2) = max(-1.,dcos(:,2));
    dcos(:,2) = min(+1.,dcos(:,2));
    dcos(:,3) = max(-1.,dcos(:,3));
    dcos(:,3) = min(+1.,dcos(:,3));

    dcos = acos(dcos) * 180. / pi ;

end

author: dengwirda


text:
%% Off-policy learning for a turbocharged diesel engine
% This is the Numerical Example 2.2 in Chapter 2 of book
% Robust Adaptive Dynamic Programming
% by Yu Jiang and Zhong-Ping Jiang

%% Parameters configuration
xn = 6;
un = 2;

% Set the weighting matrices for the cost function
Q = diag([100 0 0 0 0 100]);
R = diag([1 1]);

% Initialize the feedback gain matrix
Kinit  = zeros(un,xn);  %Only if A is Hurwitz, K can be set as zero.
K = Kinit;
N  = 100;           %Length of the window, should be at least xn^2+2*xn*un
MaxIteration = 10;  %Max iteration times
T  = 0.01;          %Length of each integration interval

x0 = [20;5;10;2;-1;-2]; %Initial condition

expl_noise_freq = (rand(un,100)-.5)*100; % Exploration noise frequencies

% Matrices to collect online data and perform learning
Dxx=[];
Ixx=[];
Ixu=[];

% Initial condition of the augmented system
X=[x0;kron(x0',x0')';kron(x0,zeros(un,1))]';

%% Run the simulation and obtain the data matrices \delta_{xx},
%I_{xx}, and I_{xu}

ifLearned = 0;

x_save=[];
t_save=[];

for i=1:N
    % Simulation the system and at the same time collect online info.
    [t,X] = ode45(@(t,x)aug_sys(t,x,K,ifLearned,expl_noise_freq), ...
        [(i-1)*T,i*T],X(end,:));

    %Append new data to the data matrices
    Dxx=[Dxx;kron(X(end,1:xn),X(end,1:xn))-kron(X(1,1:xn),X(1,1:xn))];
    Ixx=[Ixx;X(end,xn+1:xn+xn^2)-X(1,xn+1:xn+xn^2)];
    Ixu=[Ixu;X(end,xn+xn^2+1:end)-X(1,xn+xn^2+1:end)];

    % Keep track of the system trajectories
    x_save=[x_save;X];
    t_save=[t_save;t];
end

ifLearned = 1; % Mark learning finished

P_old = zeros(xn);
P = eye(xn)*10;    % Initialize the previous cost matrix
it = 0;            % Counter for iterations
p_save = [];       % Track the cost matrices in all the iterations
k_save = [];       % Track the feedback gain matrix in each iterations

% The system matrices are copied here only for the purpose of analyzing the
% results
A = [-0.4125    -0.0248        0.0741     0.0089   0          0;
    101.5873    -7.2651        2.7608     2.8068   0          0;
    0.0704       0.0085       -0.0741    -0.0089   0          0.0200;
    0.0878       0.2672        0         -0.3674   0.0044     0.3962;
    -1.8414      0.0990        0          0       -0.0343    -0.0330;
    0            0             0       -359      187.5364   -87.0316];

B = [-0.0042  0.0064
     -1.0360  1.5849
      0.0042  0;
      0.1261  0;
      0      -0.0168;
      0       0];

[Kopt,Popt] = lqr(A,B,Q,R); % Calculate the ideal solution for comparion purpose
k_save  = norm(K-Kopt);  % keep track of the differences between the actual K
% and the idea valu

%% Off-policy learning using the collected online data
while norm(P-P_old)>1e-8 & it<MaxIteration
    it = it+1;                        % Update and display the # of iters
    P_old = P;                        % Update the previous cost matrix
    QK = Q+K'*R*K;                    % Update the Qk matrix

    Theta = [Dxx,-Ixx*kron(eye(xn),K')-Ixu];  % Left-hand side of
    % the key equation
    Xi = -Ixx*QK(:);                 % Right-hand side of the key equation
    pp = pinv(Theta)*Xi;             % Solve the equations in the LS sense
    P = reshape(pp(1:xn*xn), [xn, xn]);  % Reconstruct the symmetric matrix
    P = (P + P')/2;

    BPv = pp(end-(xn*un-1):end);
    K = inv(R)*reshape(BPv,un,xn)/2;% Get the improved gain matrix
    p_save = [p_save,norm(P-Popt)];   % Keep track of the cost matrix
    k_save = [k_save,norm(K-Kopt)];    % Keep track of the control gains

    disp(['K_', num2str(it), '=']);
    disp(K);
end

%% Post-learning simulations
[tt,xx]=ode45(@(t,x)aug_sys(t,x,K,ifLearned,expl_noise_freq), ...
    [t(end) 10],X(end,:)');

% Keep track of the post-learning trajectories
t_final = [t_save;tt];
x_final = [x_save;xx];

% For comparson, also get the unlearned simulation
[ttt,xxx]=ode45(@(t,x)aug_sys(t,x,Kinit,ifLearned,expl_noise_freq), ...
    [t(end) 10],X(end,:)');
t_unlearned = [t_save;ttt];
x_unlearned = [x_save;xxx];

%% Plotting results
figure(1)
subplot(211)
plot(0:length(p_save)-1,p_save,'o',0:length(p_save)-1,p_save,'Linewidth',2);
ylim([-max(p_save)*0.2 max(p_save)*1.2]);
legend('||P_k-P^*||')
xlabel('Number of iterations')

subplot(212)
plot(0:length(k_save)-1,k_save,'^',0:length(k_save)-1,k_save,'Linewidth',2)
ylim([-max(k_save)*0.2 max(k_save)*1.2]);
legend('||K_k-K^*||')
xlabel('Number of iterations')

% Uncomment to save figure
% print('Ch2_ex2_fig1_pk','-depsc')

figure(2)
plot(t_final,x_final(:,1:6),'Linewidth',2)
legend('x_1','x_2','x_3','x_4','x_5','x_6')
xlabel('Time (sec)')

figure(3)
plot(t_final, 3.6*x_final(:,6),'k-', ...
     ttt, 3.6*xxx(:,6),'r--', ...
    'Linewidth',2)
ylim([-400 150])

legend('MAF (Under ADP)', 'MAF (Unlearned)')
xlabel('Time (sec)','FontSize',12)

% Create textarrow
annotation(figure(3),'textarrow',[0.2630173564753 0.218958611481976],...
	[0.166023166023166 0.198841698841699],'String',{'Controller Updated'},...
	'FontSize',14);

% Uncomment to save figure
% print('Ch2_ex2_fig2_y','-depsc')

%% Display results
disp('Approximate Cost Matrix')
P
disp('Optimal Cost Matrix')
Popt
disp('Approximate Gain Matrix')
K
disp('Optimal Gain Matrix')
Kopt

author: yu-jiang

repo: radpbook

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


text:
function tor_ihb_TalSpace()
% :Usage:
% ::
%
%    tor_ihb_TalSpace
%
% Main function to view SPM99 cluster's (contours & 3D) in Talariach space
%
%==================================================================================
%   Graphic objects and their 'Tags':
%==================================================================================
%   type        handle          Tag                     Description
%----------------------------------------------------------------------------------
%   figure      hFigMain        ihb_TalSpaceMain_fig    Main figure
%   figure      hFig3D          ihb_TalSpace3D_fig      3D view figure
%   figure      hFigLeg         ihb_TalSpaceLeg_fig     Legend figure
%
%   axes        hAxis3D         ihb_axesTalSpace3D      Axes for 3D view
%   light       hLight          ihb_lightTalSpace3D     Light object for 3D view
%
%   patch                       ihb_cntAxial            contours patch handles
%   patch                       ihb_cntFrontal              in 3D view window
%   patch                       ihb_cntSaggitalL
%   patch                       ihb_cntSaggitalR
%
%   patch                       ihb_surfaceCluster      isosurface patch
%
%   uicontrol                   ihb_CursorCoord         coordinates under cursor
%   uicontrol                   ihb_ClusterLevelProb    cluster probability
%   uicontrol                   ihb_VoxelLevelProb      extr. probability
%   uicontrol                   ihb_VolumeInVox         cluster volume
%   uicontrol                   ihb_ClusterPopUp        select cluster combo box
%
%   uicontrol                   ihb_clr_Orig            pushbuttons for
%   uicontrol                   ihb_clr_Symm            color change in
%   uicontrol                   ihb_clr_Surface         legend window
%   uicontrol                   ihb_clr_Axial
%   uicontrol                   ihb_clr_Frontal
%   uicontrol                   ihb_clr_Saggital
%
%   line                        ihb_LineOriginal        lines for original clusters
%   line                        ihb_LineSymmetrical     lines for symmetrical clusters
%
%   uimenu                      ihb_VoxSizeMenu111
%   uimenu                      ihb_VoxSizeMenu222
%   uimenu                      ihb_VoxSizeMenu444
%
%   uimenu                      ihb_RendererMenuZbuffer
%   uimenu                      ihb_RendererMenuOpenGL
%
%   uimenu                      ihb_TransparencyMenu
%   uimenu                      ihb_TransparencyMenuOpaque
%   uimenu                      ihb_TransparencyMenuLow
%   uimenu                      ihb_TransparencyMenuMedium
%   uimenu                      ihb_TransparencyMenuHigh
%
%   uimenu                      ihb_Axes3DMenuClusterOnly
%   uimenu                      ihb_Axes3DMenuWholeBrain
%
%   uimenu                      ihb_LineWidthMenuThin
%   uimenu                      ihb_LineWidthMenuNormal
%   uimenu                      ihb_LineWidthMenuThick
%
%==================================================================================
%  Application data (name used with setappdata the same as handle or variable name
%==================================================================================
% for hFigMain:
%
%   axes            hAxAxial            axial
%   axes            hAxFront            frontal
%   axes            hAxSagL             left saggital
%   axes            hAxSagR             right saggital
%   strucure array  clusters            array of selected clusters
%   integer         indClusterToView    index of currently viewed cluster in clusters array
%   bool            drawSymLS           == 1 if draw symmetrical on left saggital
%   bool            drawSymRS           == 1 if draw symmetrical on right saggital
%
%----------------------------------------------------------------------------------
%
% for hAxis (where hAxis is one of hAxAxial, hAxFront, hAxSagL, hAxSagR):
%   (see ihb_LoadSlice for detail)
%
%   strucure array  sliceInfo           array of Talarich slice information (ihbdfl_ax_info etc.)
%   integer array   rngArray            array of local indexes in sliceInfo
%   integer         sliceIndex          index of slice in the local range array rngArray
%   double          sliceDist           slice distance from origin
%   string          axType              axis type: 'ax' 'fr' 'sl' 'sr'
%   1x2 vector      xLim                limit values for X direction
%   1x2 vector      yLim                limit values for Y direction
%
%
%==================================================================================
% Default information in ihb_TalSpaceDfl.mat file (see ihb_ResetDefaults)
%==================================================================================
% To get default variable with name var use:
%
% load('ihb_TalSpaceDfl.mat', 'var');
%
% To save:
%
% save('ihb_TalSpaceDfl.mat', 'var', '-append');
%----------------------------------------------------------------------------------
% Figures tags for windows created by SPM during call spm_getSPM
%----------------------------------------------------------------------------------
% ihbdfl_spm_fig_Interactive
% ihbdfl_spm_fig_SelFileWin
% ihbdfl_spm_fig_ConMan
%----------------------------------------------------------------------------------
% Minimal size of clusters to draw
%----------------------------------------------------------------------------------
% ihbdfl_min_size_to_draw
%----------------------------------------------------------------------------------
% Sizes of voxels used in Talariach space
%----------------------------------------------------------------------------------
% ihbdfl_tal_x_vox
% ihbdfl_tal_y_vox
% ihbdfl_tal_z_vox
%----------------------------------------------------------------------------------
% Sizes of main window and axis for ihb_TalSpace
%----------------------------------------------------------------------------------
% ihbdfl_main_width
% ihbdfl_main_height
% ihbdfl_main_gap
% ihbdfl_main_x
% ihbdfl_main_y
% ihbdfl_main_z
% ihbdfl_main_info_h
%----------------------------------------------------------------------------------
% Talariach slices information
%----------------------------------------------------------------------------------
% ihbdfl_ax_info
% ihbdfl_fr_info
% ihbdfl_sl_info
% ihbdfl_sr_info
%----------------------------------------------------------------------------------
% Axis limits to draw slice
%----------------------------------------------------------------------------------
% ihbdfl_xLimD
% ihbdfl_yLimD
% ihbdfl_zLimD
% ihbdfl_xLimI
% ihbdfl_yLimI
% ihbdfl_zLimI
%----------------------------------------------------------------------------------
% 3D view Renderer ('zbuffer' or 'OpenGL')
% 3D view transparency (available only for OpenGL
%----------------------------------------------------------------------------------
% ihbdfl_renderer
% ihbdfl_transparency
%----------------------------------------------------------------------------------
% Axes for 3D view type. If ihbdfl_bAxesIsClusterOnly == 1 - axes have cluster range
%                                                        0 - axes have brain range
%----------------------------------------------------------------------------------
% ihbdfl_bAxesIsClusterOnly
%----------------------------------------------------------------------------------
% Colors defaults
%----------------------------------------------------------------------------------
% ihbdfl_color_clOrig
% ihbdfl_color_clSymm
% ihbdfl_color_surface
% ihbdfl_color_cntAx
% ihbdfl_color_cntFr
% ihbdfl_color_cntSag
%----------------------------------------------------------------------------------
% Line width
%----------------------------------------------------------------------------------
% ihbdfl_line_width
%==================================================================================
%
% ..
%   09.04.01    Sergey Pakhomov
%   01.08.01    last modified
% ..

% ..
%    Create Default information file if necessary
% ..
dflFileName = ihb_FileFolderName('dfl');
fid = fopen(dflFileName);
if fid == -1
    ihb_ResetDefaults;
else
    fclose(fid);
end
%==================================================================================
% Delete previous figures
%==================================================================================
ihb_ResetFigure('ihb_TalSpaceMain_fig');
ihb_ResetFigure('ihb_TalSpace3D_fig');
ihb_ResetFigure('ihb_TalSpaceLeg_fig');
%==================================================================================
% Create main figure
%==================================================================================
screenUnits = get(0, 'units');      % Save current screen units
set(0, 'units', 'pixels');          % Set pixels units for screen
scrSize = get(0,'ScreenSize');      % Get screen size in pixels
set(0, 'units', screenUnits);       % Restore current screen units

load(dflFileName, 'ihbdfl_main_width');
load(dflFileName, 'ihbdfl_main_height');
pos(3) = ihbdfl_main_width;
pos(4) = ihbdfl_main_height;
pos(1) = (scrSize(3) - pos(3))/2;
pos(2) = (scrSize(4) - pos(4))/2;

hFigMain = figure(...
    'Visible', 'off',...
    'NumberTitle', 'off',...
    'MenuBar', 'none',...
    'Name', 'Talairach slices',...
    'Tag', 'ihb_TalSpaceMain_fig',...
    'CloseRequestFcn', 'ihb_CloseTalSpaceCB',...
    'ResizeFcn', 'ihb_ResizeTalSpaceCB',...
    'units', 'pixels',...
    'Position', pos);
%==================================================================================
% Set menu
%==================================================================================
%----------------------------------------------------------------------------------
% File menu
%----------------------------------------------------------------------------------
hFileMenu = uimenu('Label', 'File');
uimenu(hFileMenu, 'Label', 'New clusters', 'Callback', 'ihb_GetClusterSet');
uimenu(hFileMenu, 'Label', 'Exit', 'Callback', 'close', 'Separator', 'on');
%----------------------------------------------------------------------------------
% 3D View menu
%----------------------------------------------------------------------------------
h3DViewMenu = uimenu('Label', '3D View', 'Callback', 'ihb_TalSpaceView3D');
%----------------------------------------------------------------------------------
% Color Legend menu
%----------------------------------------------------------------------------------
hLegendMenu = uimenu('Label', 'Colors', 'Callback', 'ihb_ColorsTalSpace');
%----------------------------------------------------------------------------------
% Option menu
%----------------------------------------------------------------------------------
hOptionMenu = uimenu('Label', 'Options');
	%------------------------------------------------------------------------------
	%     Voxel size submenu
	%------------------------------------------------------------------------------
    hVoxSizeMenu = uimenu(hOptionMenu, 'Label', 'Voxel size');
        uimenu(hVoxSizeMenu, 'Tag', 'ihb_VoxSizeMenu111', 'Label', '1x1x1', 'Callback', 'ihb_TalSpaceVoxSizeMenuCB(1)');
        uimenu(hVoxSizeMenu, 'Tag', 'ihb_VoxSizeMenu222', 'Label', '2x2x2', 'Callback', 'ihb_TalSpaceVoxSizeMenuCB(2)');
        uimenu(hVoxSizeMenu, 'Tag', 'ihb_VoxSizeMenu444', 'Label', '4x4x4', 'Callback', 'ihb_TalSpaceVoxSizeMenuCB(4)');
	%------------------------------------------------------------------------------
	%     Renderer submenu
	%------------------------------------------------------------------------------
    hRendererMenu  = uimenu(hOptionMenu, 'Label', '3D renderer');
        uimenu(hRendererMenu, 'Tag', 'ihb_RendererMenuZbuffer', 'Label', 'Zbuffer', 'Callback', 'ihb_TalSpaceRendererMenuCB(1)');
        uimenu(hRendererMenu, 'Tag', 'ihb_RendererMenuOpenGL', 'Label', 'OpenGL', 'Callback', 'ihb_TalSpaceRendererMenuCB(2)');
	%------------------------------------------------------------------------------
	%     Transparency submenu
	%------------------------------------------------------------------------------
    hTransparencyMenu  = uimenu(hOptionMenu, 'Tag', 'ihb_TransparencyMenu', 'Label', '3D Transparency');
        uimenu(hTransparencyMenu, 'Tag', 'ihb_TransparencyMenuOpaque', 'Label', 'Opaque', 'Callback', 'ihb_TalSpaceTransparencyMenuCB(1.0)');
        uimenu(hTransparencyMenu, 'Tag', 'ihb_TransparencyMenuLow', 'Label', 'Low', 'Callback', 'ihb_TalSpaceTransparencyMenuCB(0.8)');
        uimenu(hTransparencyMenu, 'Tag', 'ihb_TransparencyMenuMedium', 'Label', 'Medium', 'Callback', 'ihb_TalSpaceTransparencyMenuCB(0.5)');
        uimenu(hTransparencyMenu, 'Tag', 'ihb_TransparencyMenuHigh', 'Label', 'High', 'Callback', 'ihb_TalSpaceTransparencyMenuCB(0.2)');
	%------------------------------------------------------------------------------
	%     3D axes submenu
	%------------------------------------------------------------------------------
    hAxes3DMenu  = uimenu(hOptionMenu, 'Label', '3D view axes');
        uimenu(hAxes3DMenu, 'Tag', 'ihb_Axes3DMenuClusterOnly', 'Label', 'Cluster Only', 'Callback', 'ihb_TalSpaceAxes3DMenuCB(1)');
        uimenu(hAxes3DMenu, 'Tag', 'ihb_Axes3DMenuWholeBrain', 'Label', 'Whole Brain', 'Callback', 'ihb_TalSpaceAxes3DMenuCB(0)');
	%------------------------------------------------------------------------------
	%     Smoothing submenu
	%------------------------------------------------------------------------------
    % hSmoothMenu = uimenu(hOptionMenu, 'Label', 'Smoothing');
	%------------------------------------------------------------------------------
	%     Line width submenu
	%------------------------------------------------------------------------------
    hLineWidthMenu = uimenu(hOptionMenu, 'Label', 'Line width');
        uimenu(hLineWidthMenu, 'Tag', 'ihb_LineWidthMenuThin', 'Label', 'Thin', 'Callback', 'ihb_TalSpaceLineWidthMenuCB(1)');
        uimenu(hLineWidthMenu, 'Tag', 'ihb_LineWidthMenuNormal', 'Label', 'Normal', 'Callback', 'ihb_TalSpaceLineWidthMenuCB(2)');
        uimenu(hLineWidthMenu, 'Tag', 'ihb_LineWidthMenuThick', 'Label', 'Thick', 'Callback', 'ihb_TalSpaceLineWidthMenuCB(4)');
    hResetMenu = uimenu(hOptionMenu, 'Label', 'Reset Defaults', 'Callback', 'ihb_TalSpaceResetMenuCB');
%----------------------------------------------------------------------------------
% Help menu
%----------------------------------------------------------------------------------
hHelpMenu  = uimenu('Label', 'Help');
    uimenu(hHelpMenu, 'Label', 'Web site', 'Callback', 'ihb_TalSpaceHelp');
    uimenu(hHelpMenu, 'Label', 'About', 'Callback', 'ihb_AboutTalSpace');
%----------------------------------------------------------------------------------
ihb_TalSpaceSetMenuCheck;
%----------------------------------------------------------------------------------
% Set axis and information UI
%----------------------------------------------------------------------------------
ihb_SetAxis(hFigMain);
ihb_SetInfoUI(hFigMain);
%----------------------------------------------------------------------------------
% Set mouse move callback
%----------------------------------------------------------------------------------
set(hFigMain, 'WindowButtonMotionFcn', 'ihb_MotionTalSpaceCB');
%----------------------------------------------------------------------------------
% Select cluster set and cluster to view    while
%----------------------------------------------------------------------------------
bIsEmpty = 0;
while bIsEmpty == 0
    bIsEmpty = tor_ihb_GetClusterSet;
    if bIsEmpty == 0
        h = msgbox('Cluster to view should be selected', 'Warning!', 'warn', 'modal');
        waitfor(h);
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function asa144_test ( )

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

  asa144_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA144_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function centroid = triangle_centroid_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_CENTROID_2D computes the centroid of a triangle in 2D.
%
%  Discussion:
%
%    The centroid of a triangle can also be considered the
%    center of gravity, or center of mass, assuming that the triangle
%    is made of a thin uniform sheet of massy material.
%
%    The centroid of a triangle is the intersection of the medians.
%
%    A median of a triangle is a line connecting a vertex to the
%    midpoint of the opposite side.
%
%    In barycentric coordinates, in which the vertices of the triangle
%    have the coordinates (1,0,0), (0,1,0) and (0,0,1), the centroid
%    has coordinates (1/3,1/3,1/3).
%
%    In geometry, the centroid of a triangle is often symbolized by "G".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real CENTROID(2,1), the coordinates of the centroid.
%
  dim_num = 2;

  for i = 1 : dim_num
    centroid(i,1) = sum ( t(i,1:3) ) / 3.0;
  end

  return
end


text:
function test_failed = test_wpfbt(verbose)
%TEST_WFBTPR
%
% Checks perfect reconstruction of the general wavelet transform of different
% filters
%
disp('========= TEST WPFBT ============');
global LTFAT_TEST_TYPE;
tolerance = 1e-8;
if strcmpi(LTFAT_TEST_TYPE,'single')
   tolerance = 2e-6;
end

test_failed = 0;
if(nargin>0)
   verbose = 1;
else
   verbose = 0;
end

type = {'dec'};
ext = {'per','zero','odd','even'};
format = {'pack','cell'};

J = 3;
%! Mild tree
wt1 = wfbtinit({'db10',6,'full'});
wt1 = wfbtremove(2,1,wt1,'force');
wt1 = wfbtremove(2,3,wt1,'force');

%! Hardcore tree
wt2 = wfbtinit({'db3',1});
wt2 = wfbtput(1,1,'mband1',wt2);
wt2 = wfbtput(2,2,'mband1',wt2);
wt2 = wfbtput(3,3,'mband1',wt2);
wt2 = wfbtput(3,1,'db10',wt2);
wt2 = wfbtput(4,1,'dgrid2',wt2);
wt2 = wfbtput(5,1,'db3',wt2);

%! Another tree
wt3 = wfbtinit();
wt3 = wfbtput(0,0,'cmband4',wt3);
wt3 = wfbtput(1,0,'cmband6',wt3);
wt3 = wfbtput(1,1,'cmband6',wt3);
wt3 = wfbtput(2,0,'cmband4',wt3);
wt3 = wfbtput(2,1,'cmband4',wt3);
wt3 = wfbtput(3,1,'cmband4',wt3);
wt3 = wfbtput(3,2,'cmband4',wt3);
wt3 = wfbtput(3,3,'cmband4',wt3);
wt3 = wfbtput(3,4,'cmband4',wt3);

% wt2 = wfbtinit();
% wt2 = wfbtput(0,0,{'db',4},wt2);
% wt2 = wfbtput(1,0,{'algmband',1},wt2);
% wt2 = wfbtput(1,1,{'hden',3},wt2);
% wt2 = wfbtput(2,0,{'dgrid',2},wt2);
% wt2 = wfbtput(2,1,{'dgrid',2},wt2);

test_filters = {
               {'algmband2',J} % 4 filters, uniform, crit. sub.
               {'db4',J}
               {'algmband1',J} % 3 filters, uniform, crit. sub.
               %{{'hden',3},J} % 3 filters, non-uniform, no crit. sub. no correct
               {'dgrid1',J} % 4 filters. sub. fac. 2
               wt1
               wt2
               wt3
               };

scaling = {'intscale','intsqrt','intnoscale'};
scalingInv = scaling(end:-1:1);

%testLen = 4*2^7-1;%(2^J-1);
testLen = 53;
f = tester_rand(testLen,1);
for scIdx = 1:numel(scaling)
for extIdx=1:length(ext)
   extCur = ext{extIdx};

   for typeIdx=1:length(type)
     for tt=1:length(test_filters)
        actFilt = test_filters{tt};
         if verbose, if(~isstruct(actFilt))fprintf('J=%d, filt=%s, ext=%s, inLen=%d \n',actFilt{2},actFilt{1},extCur,size(f,1)); else disp('Custom'); end; end;

        [c,info] = wpfbt(f,actFilt,extCur,scaling{scIdx});
        fhat = iwpfbt(c,actFilt,size(f,1),extCur,scalingInv{scIdx});

        %MSE
            err = norm(f-fhat,'fro');
            [test_failed,fail]=ltfatdiditfail(err,test_failed,tolerance);
            if(~verbose)
              if(~isstruct(actFilt))fprintf('J=%d, %5.5s, ext=%4.4s, %s, L=%d, err=%.4e %s \n',actFilt{2},actFilt{1},extCur,scaling{scIdx},size(f,1),err,fail); else fprintf('Custom, %s, err=%.4e %s \n',scaling{scIdx},err,fail); end;
            end
            if strcmpi(fail,'FAILED')
               if verbose
                if(~isstruct(actFilt)) fprintf('err=%d, filt=%s, ext=%s, inLen=%d \n',err,actFilt{1},extCur,testLen); else disp('Custom'); end;
                 figure(1);clf;stem([f,fhat]);
                 figure(2);clf;stem([f-fhat]);
                 break;
               end
            end

           fhat2 = iwpfbt(c,info);
           err = norm(f-fhat2,'fro');
           [test_failed,fail]=ltfatdiditfail(err,test_failed,tolerance);
           if(~isstruct(actFilt))
                fprintf('INFO J=%d, %5.5s, ext=%4.4s, %s, L=%d, err=%.4e %s \n',actFilt{2},actFilt{1},extCur,scaling{scIdx},size(f,1),err,fail);
           else
               fprintf('INFO Custom, %s, err=%.4e %s \n',scaling{scIdx},err,fail);
           end;

            if test_failed && verbose, break; end;

     end
     if test_failed && verbose, break; end;
   end
   if test_failed && verbose, break; end;
end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function a = wishart_sample_inverse ( m, df, sigma )

%*****************************************************************************80
%
%% WISHART_SAMPLE_INVERSE returns the inverse of a sample Wishart matrix.
%
%  Discussion:
%
%    This function requires functions from the PDFLIB and RNGLIB libraries.
%
%    The "initialize()" function from RNGLIB must be called before using
%    this function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Odell, Alan Feiveson,
%    A numerical procedure to generate a sample covariance matrix,
%    Journal of the American Statistical Association,
%    Volume 61, Number 313, March 1966, pages 199-203.
%
%    Stanley Sawyer,
%    Wishart Distributions and Inverse-Wishart Sampling,
%    Washington University,
%    30 April 2007, 12 pages.
%
%  Parameters:
%
%    Input, integer M, the order of the matrix.
%
%    Input, integer DF, the number of degrees of freedom.
%    M <= DF.
%
%    Input, real SIGMA(M,M), the covariance matrix, which should be
%    a symmetric positive definite matrix.
%
%    Output, real A(M,M), the inverse of a sample matrix from the Wishart
%    distribution.
%
  if ( df < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WISHART_SAMPLE_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  DF = %d < M = %d.\n', df, m );
    error ( 'WISHART_SAMPLE_INVERSE - Fatal error!\n' );
  end
%
%  Get R, the upper triangular Cholesky factor of SIGMA.
%
  r = chol ( sigma );
%
%  Get S, the inverse of R.
%
  s = r8ut_inverse ( m, r );
%
%  Get UA, the inverse of a sample from the unit Wishart distribution.
%
  ua = wishart_unit_sample_inverse ( m, df );
%
%  Construct the matrix A = S * UA * S'.
%
  a = s * ua * s';

  return
end


text:
function [ value, isterminal, direction ] = event_function ( t, y )

  value = y(1) - 0.99;
  isterminal = 1;
  direction = 0;

  return
end


text:
% fmri_data: Data class for storing and analyzing neuroimaging datasets and associated information
%
% 'fmri_data' is a data class containing information about generic fmri
% datasets stored in a structure-like object.  Using this has the
% advantages that the properties and methods are standardized and controlled.
% It also keeps track of the history of what was done to the dataset.
%
% -------------------------------------------------------------------------
% Features and philosophy:
% -------------------------------------------------------------------------
%
% - Store image data in a flat (2-d), space-efficient voxels x images matrix
% - Save space: Tools to remove out-of-image voxels and empty (zero/NaN) voxels, store data in single-precision format
% - Analysis-friendly: 2-d matrices can be read and analyzed in multiple packages/algorithms
% - Meta-data included to convert back to 3-d image volume space,
%   with easy tools (methods) to reconstruct and visualize images
%   Built-in resampling makes it easy to compare/combine datasets with different voxel sizes and image bounding boxes
%   Reduces overhead for statisticians/data scientists unfamiliar with neuroimaging to apply their algorithms.
% - Multiple images can be stored in a single object
% - Methods have short, intuitive names, and perform high-level functions specialized for neuroimaging:
% - Visualization (plot, orthviews, surface, montage, histogram, isosurface methods)
% - Image manipulation (apply_mask, get_wh_image, resample_space, compare_space, flip, threshold methods)
% - Data extraction (apply_atlas, apply_parcellation, extract_gray_white_csf, extract_roi_averages)
% - Analysis (ica, mahal, image_math, and many more in the fmri_data subclass)
% - Provenance: Ability to track and update history of changes to objects
%
% fmri_data is a subclass of the image_vector object and inherits its
% properties and methods.
%
% -------------------------------------------------------------------------
% To construct/create a new instance of an object:
% -------------------------------------------------------------------------
%
% Basic usage:
% obj = fmri_data(image_names, [maskinput], [other optional inputs])
%
% maskinput     :   [optional] name of mask image to use.  Default: 'brainmask_canlab.nii', a
%                   brain mask that is distributed with SPM software and in
%                   the CANlab Core Tools
%                   Alternative in CANlab tools: which('gray_matter_mask.img')
% 'noverbose'   :   Suppress verbose output during image loading
% 'sample2mask' :   Sample images to mask space. Default: Sample mask to
%                   image space, use native image space
%
% Creating class instances
% -----------------------------------------------------------------------
% You can create an empty object by using:
% fmri_dat = fmri_data
% - fmri_dat is the object.
% - It will be created with a standard brain mask, brainmask_canlab.nii
% - This image should be placed on your Matlab path
% - The space information is stored in fmri_dat.volInfo
% - Data is stored in fmri_dat.dat, in a [voxels x images] matrix
% - You can replace or append data to the fmri_dat.dat field.
%
% You can create an object by assembling an image_vector object from parts
% (entering fields) and converting using fmri_obj = fmri_data(image_vec_obj)
% This can also be a structure with field names that match properties of
% the new intended fmri_data object. e.g.,
%   new_obj.volInfo = MC_Setup.volInfo;
%   new_obj.dat = MC_Setup.unweighted_study_data(:, sym_cons);
%
% You can create an fmri_data object with extacted image data.
% - Let "imgs" be a string array or cell array of image names
% - This command creates an object with your (4-D) image data:
% - fmri_dat = fmri_data(imgs);
% - Images can be zipped (.gz) as well. fmri_data() will unpack them.
% - Only values in the standard brain mask, brainmask_canlab.nii, will be included.
% - This saves memory by reducing the number of voxels saved dramatically.
%
% You can specify any mask you'd like to extract data from.
% - Let "maskimagename" be a string array with a mask image name.
% - this command creates the object with data saved in the mask:
% - fmri_dat = fmri_data(imgs, maskimagename);
% - The mask information is saved in fmri_dat.mask
%
% e.g., this extracts data from images within the standard brain mask:
% dat = fmri_data(imgs, which('brainmask_canlab.nii'));
%
% Defining the space of the extracted data
% -----------------------------------------------------------------------
% Note: There are two options for defining the space (i.e., coordinates/voxels)
% that the data is mapped to.
% By default, the mask is resliced to the same space as the first image in the
% input image name set (not coregistered; just resliced to the same voxel sizes.
% The images are assumed to be in register.)
%
% Reampling to mask space: YOU CAN ALSO map the image data to the space of the mask, by entering
% 'sample2mask' as in input argument.
% For loading images in different spaces together in one object, use the 'sample2mask' option.
%
% Attaching additional data
% -----------------------------------------------------------------------
% The fmri_data object has a number of fields for appending specific types of data.
%
% - You can replace or append data to the fmri_dat.dat field.
% - The fmri_data object will also store predictor data (.X) also outcome data (.Y)
% - There are many fields for descriptions, notes, etc., like "dat_descrip" and "source_notes"
% - Attach custom descriptions in these fields to document your object.
% - The "history" field stores a cell array of strings with the processing
% history of the object. Some methods add to this history automatically.
%
% -----------------------------------------------------------------------
% Properties and methods
% -----------------------------------------------------------------------
% Properties are data fields associated with an object.
% Type the name of an object (class instance) you create to see its
% properties, and a link to its methods (things you can run specifically
% with this object type). For example: After creating an fmri_data object
% called fmri_dat, as above, type fmri_dat to see its properties.
%
% There are many other methods that you can apply to fmri_data objects to
% do different things.
% - Try typing methods(fmri_data) for a list.
% - You always pass in an fmri_data object as the first argument.
% - Methods include utilities for many functions - e.g.,:
% - resample_space(fmri_dat) resamples the voxels
% - write(fmri_dat) writes an image file to disk (careful not to overwrite by accident!)
% - regress(fmri_dat) runs multiple regression
% - predict(fmri_dat) runs cross-validated machine learning/prediction algorithms
%
% Key properties and methods (a partial list; type doc fmri_data for more):
% -------------------------------------------------------------------------
% fmri_data Properties (a partial list; type doc fmri_data for more):
%   dat                     - Image data, a [voxels x images] matrix, single-format
%   fullpath                - List of image names loaded into object with full paths
%   history                 - History of object processing, for provenance
%   image_names             - List of image names loaded into object, no paths
%   removed_images          - Vector of images that have been removed (saves space; see remove_empty.m, replace_empty.m)
%   removed_voxels          - Vector of empty in-mask voxels that have been removed (saves space; see remove_empty.m, replace_empty.m)
%   volInfo                 - Structure with info on brain mask (saves space) and mapping voxels to brain space
%
% fmri_data Methods (a partial list; type doc fmri_data for more):
%   General:
% . 	descriptives        -  Get descriptives for an fmri_data or other image_vector object
%       enforce_variable_types	-  Re-casts variables in objects into standard data types, which can save
%       flip                - Flips images stored in an object left-to-right
%     	history             - Display history for image_vector object
%       write               - Write an image_vector object to hard drive as an Analyze image (uses .fullpath field for image names)
%   Data extraction:
%       apply_atlas         - Computes the mean value or pattern expression for each reference region specified in an atlas object
%       apply_mask          - Apply a mask image (image filename or fmri_mask_image object) to an image_vector object
%       apply_parcellation  - Computes the mean value or pattern expression for each parcel specified in a data object
%       extract_gray_white_csf	- Extracts mean values (values) and top 5 component scores (components)
%       extract_roi_averages	- This image_vector method a extracts and averages data stored in an fmri_data object
%
%   Handling brain space and image selection:
%       compare_space           - Compare spaces of two image_vector objects
%       get_wh_image            - For an image_vector with multiple images (cases, contrasts, etc.), select a subset.
%       reconstruct_image       - Reconstruct a 3-D or 4-D image from image_vector object obj
%       remove_empty            - remove vox: logical vector of custom voxels to remove, VOX x 1
%       reparse_contiguous      - Re-construct list of contiguous voxels in an image based on in-image
%       replace_empty           - Replace empty/missing values in an image data object
%       resample_space          - Resample the images in an fmri_data object (obj) to the space of another
%
%   Display and visualization:
%   	display_slices      - Creates 3 separate montage views - ax, cor, sagg in a special figure window
%       histogram           - Create a histogram of image values or a series of histograms for each
%       image_similarity_plot - Associations between images in object and set of 'spatial basis function' images (e.g., 'signatures' or pre-defined maps)
%       isosurface          - Create and visualize an isosurface created from the boundaries in an image object.
%       montage             - Create a montage of an image_vector (or statistic_image or fmri_data)
%       orthviews               - display SPM orthviews for CANlab image_vector (or fmri_data, statistic_image) object
%       pattern_surf_plot_mip	- axial maximum intensity projection pattern surface plot
%   	sagg_slice_movie	- Movie of successive differences (sagittal slice)
%       slices              - Create a montage of single-slice results for every image in an image_vector object
%       surface             - Render image data on brain surfaces; options for cutaways and canonical surfaces
%       wedge_plot_by_atlas	- Plot a data object or 'signature' pattern divided into local regions
%
%   Data processing and analysis:
%   	ica                 - Spatial ICA of an fmri_data object
%       image_math          - Perform simple mathematical and boolean operations on image objects (see also plus, minus, power)
%       mahal               - Mahalanobis distance for each image in a set compared to others in the set
%       mean                - Mean across a set of images. Returns a new image_vector object.
%       preprocess          - Preprocesses data in an image_vector (e.g., fmri_data) object; many options for filtering and outlier id
%       qc_metrics_second_level	- Quality metrics for a 2nd-level analysis (set of images from different subjects)
%       searchlight         - Run searchlight multivariate prediction/classification on an image_vector
%       threshold           - Threshold image_vector (or fmri_data or fmri_obj_image) object based on raw threshold values
%       union               - ...and intersection masks for two image_vector objects
%
% -------------------------------------------------------------------------
% Examples and help:
% -------------------------------------------------------------------------
%
% To list properties and methods for this object, type:
% doc fmri_data, methods(fmri_data)
%
% Example 1: Load images (and run a simple analysis)
% % --------------------------------------------------------------------
%
% Load a sample dataset into an fmri_data object (subclass of image_vector)
% This loads one of a set of named image collections used in demos/help:
% data_obj = load_image_set('emotionreg');
%
% You can load the same images manually, by locating the files, listing
% their names in a character array (or 1 x n cell array of strings), and
% then passing those into fmri_data:
%
% data_obj = fmri_data(which('Wager_2008_emo_reg_vs_look_neg_contrast_images.nii.gz'));
%
% filedir = what(fullfile('CanlabCore', 'Sample_datasets', 'Wager_et_al_2008_Neuron_EmotionReg'));
% image_names = filenames(fullfile(filedir.path, '*img'));
% data_obj = fmri_data(image_names);
%
% Now you can interact with the object.  Try, e.g.,:
% methods(data_obj)                               % List methods for object type
% descriptives(data_obj);                         % Print summary of descriptive statistics for the dataset
% plot(data_obj)                                  % Custom fmri_data specific plots
% t = ttest(data_obj);                            % Perform a voxel-wise one-sample t-test across images
% t = threshold(t, .005, 'unc', 'k', 10);         % Re-threshold with extent threshold of 10 contiguous voxels
% r = region(t);                                  % Turn t-map into a region object with one element per contig region
%
% Example 2: Extract data averaged over regions of interest:
% % --------------------------------------------------------------------
%
% First run Example 1.  Now you have a thresholded t-statistic map.
% Extract averages (across voxels) for each subject in each contiguous
% region by typing:
%
% r = extract_roi_averages(data_obj, t);
%
% This returns r, which is another object--a "region"-class object.
% Region objects contain vectors, one element per pre-defined region in the
% image (in this case, significant blobs from our analysis).
% Each element contains info that describes the region, including the voxels included,
% and average and voxel-by-voxel data if extracted from images and attached.
% Type "doc region" for more info.
% r has properties that hold multiple types of data:
% - .dat holds generic extracted data
% - .all_data holds voxel-by-voxel extracted data
%
% % --------------------------------------------------------------------
% % Group one-sample t-test:
% % A simple, complete example of a group analysis
% % For more, see walkthroughs on canlab.github.io
% % --------------------------------------------------------------------
% % Load sample images, creating and fmri_data class object with 30 images
%     imgs = load_image_set('emotionreg');
%
% % Display a slice from each image in a montage:
%     slices(imgs);
%
% % Display some useful summary plots of the dataset:
%     plot(imgs);
%
% % Perform a t-test on each voxel, returning a statistic_image object
% % containing t-stats and p-values:
%     t = ttest(imgs);
%
% % Display the unthresholded results in a quick-render montage of image
% % values only:
%     display_slices(t, 'axial'); colormap summer; colorbar;
%
% % Display the unthresholded results over an anatomical underlay,
% % on a combination of slices and surfaces,
% % returning an fmridisplay class object with registered handles
%     o2 = canlab_results_fmridisplay(t, 'full');
%
% % Remove colors from slices and surfaces registered in the o2 object:
%     o2 = removeblobs(o2);
%
% % Threshold the t-statistic_image object at p < 0.005
%     t = threshold(t, .005, 'unc');
%
% % Re-display the thresholded images on slices/surfaces registered in o2:
%     o2 = addblobs(o2, region(t), 'nolegend');
%
% % Display the thresholded t-map with orthviews:
%     orthviews(t);
%
% % display on a slice montage:
%    create_figure('slices'); axis off;
%    montage(t);
%
% % Re-threshold at q < 0.05 FDR, and re-display on orthviews:
%     t = threshold(t, .05, 'fdr');
%     orthviews(t);
%
% % Print a table of results, and return a region-class object r with labels
% % from a default atlas (an atlas-class object):
%     r = table(t);
%
% % Display a slice montage showing each activation blob, with labels:
%   montage(r, 'regioncenters', 'colormap');
%
% --------------------------------------------------------------------
% For more examples and walkthroughs, see
% walkthroughs on canlab.github.io
% These are also found in the CANlab_help_examples
% repository at https://github.com/canlab/CANlab_help_examples
%
% Some example tutorials:
% canlab_help_1_installing_tools
% canlab_help_2_load_a_sample_dataset
% canlab_help_3_voxelwise_t_test_walkthrough
% canlab_help_4_write_data_to_image_file_format
% canlab_help_5_regression_walkthrough
% ... and more
% % --------------------------------------------------------------------

% Programmers' notes:
% Tor Wager, 1/14/17 : Previously, if you stack names of images in different spaces and load them,
% fmri_data did not return an error, but it returned distorted/incorrectly
% loaded images.  Now, it returns an error. For loading images in different
% spaces together, use the 'sample2mask' option.
%
% Stephan Geuter, 5/16/17: Processing of masks in line 270 had been
% changed, masks were ignored. modified check in line 270 to include
% masking again

classdef fmri_data < image_vector

    properties
        % also inherits the properties of image_vector.

        source_notes = 'Source notes...';

        X % legacy; temporary, so we can load old objects

        mask = fmri_mask_image;
        mask_descrip = 'mask is an fmri_mask_image object that defines the mask.';

        images_per_session

        Y = [];
        Y_names;
        Y_descrip = 'Behavioral or outcome data matrix.';

        covariates;
        covariate_names = {''};
        covariates_descrip = 'Nuisance covariates associated with data';

        history_descrip = 'Cell array: names of methods applied to this data, in order';

        additional_info = struct('');

        metadata_table = table(); % a table for storing image-level metadata. Numbers of rows should be the number of images in .dat

    end % properties

    methods

        % Class constructor
        function obj = fmri_data(image_names, maskinput, varargin)
            %
            % [obj, cl_with_averages] = fmri_data(image_names, mask_image, varargin)
            %
            % Reads a set of image files and a mask image, and returns
            % an fmri_data object with data for all in-mask voxels.

            % ---------------------------------
            % Create empty fmri_data object, and return if no additional
            % arguments
            % ---------------------------------

            obj.source_notes = 'Info about image source here';
            obj.mask = fmri_mask_image;
            obj.mask_descrip = 'Volume and in-area mask info from iimg_read_img';

            obj.X = []; % legacy; temporary, so we can load old objects
            obj.Y = [];
            obj.Y_names;
            obj.Y_descrip = 'Behavioral or outcome data matrix.';
            obj.covariates;
            obj.covariate_names = {''};
            obj.covariates_descrip = 'Nuisance covariates associated with data';

            obj.images_per_session = [];

            obj.history = {''};
            obj.history_descrip = 'Cell array of names of methods applied to this data, in order';
            obj.additional_info = struct('');

            % DEFAULT INPUTS
            % -----------------------------------

            verbose = 1;
            verbosestr = 'verbose';
            sample2mask = 0;

            % SET UP OPTIONAL INPUTS
            % -----------------------------------

            for i = 1:length(varargin)
                if ischar(varargin{i})
                    switch varargin{i}

                        case 'verbose', varargin{i} = []; % nothing else needed
                        case 'noverbose', verbose = 0; verbosestr = 'noverbose'; varargin{i} = [];
                        case 'sample2mask', sample2mask = 1; varargin{i} = [];

                        case 'native_image_space' % do nothing, for convenience in calling scripts

                        otherwise, warning(['Unknown input string option:' varargin{i}]);
                    end
                end
            end

            if nargin == 0
                % SPECIAL DEFAULT METHOD: Standard empty mask
                % -----------------------------------

                % Empty: Define with standard default mask
                [image_names, maskinput] = deal(which('brainmask_canlab.nii'));

                if isempty(maskinput)
                    disp('Warning: Cannot find brainmask_canlab.nii, creating without mask info.');
                    return
                end

                verbose = 0;
                verbosestr = 'noverbose';

            elseif isempty(image_names)
               % We must have image names to run
               disp('fmridata: image_names is empty. No images to load.');
               return

            end

            % ---------------------------------
            % Special: if existing image_vector
            % map into space of fmri_data and return
            % ---------------------------------
            if iscell(image_names), image_names = char(image_names{:}); end

            if isstruct(image_names) || isa(image_names, 'image_vector')
                % Map fields of input object into fmri_data structure

                warning off
                obj2 = struct(image_names); %image_names;  % tor: struct not needed i think. nope, is needed.
                warning on

                N = fieldnames(obj);
                for i = 1:length(N)
                    if isfield(obj2, (N{i}))
                        obj.(N{i}) = obj2.(N{i});
                    end
                end

                obj.mask.volInfo = obj2.volInfo;

                obj = run_checks_and_fixes(obj, verbosestr);

                return

            else
                % It's a file. Unzip if needed and check that image_names exist

                % Handle .gz by unzipping if needed
                [image_names, was_gzipped] = gunzip_image_names_if_gz(image_names);

            end

            % ---------------------------------
            % define mask object
            % ---------------------------------

            % if all keywords, like 'noverbose', then varargin will be empty
            % SG changed check, maskinput can be defined and varargin empty.
            % maskinput is not part of varargin. 5/16/17

            % 5/16/17 Tor: if all varargin cells are empty, varargin may not
            % be... also, deal with special case where 2nd input is not
            % mask, but keyword.

            % Special case: 2nd argument is keyword, not mask.  This is
            % improper usage, but allow it for legacy reasons and
            % usability.
            if exist('maskinput', 'var') && ~isempty(maskinput) && ischar(maskinput)

                switch maskinput

                    case 'verbose', maskinput = []; % nothing else needed
                    case 'noverbose', verbose = 0; verbosestr = 'noverbose'; maskinput = [];
                    case 'sample2mask', sample2mask = 1; maskinput = [];
                end
            end

            % Empty mask: use default
            if (nargin < 2 || isempty(maskinput)) % && isempty(varargin)

                maskinput = which('brainmask_canlab.nii');
                if verbose, fprintf('Using default mask: %s\n', maskinput); end
                if isempty(maskinput), error('Cannot find mask image!'); end

            end

            switch class(maskinput)
                case 'char' % string file name
                    maskobj = fmri_mask_image(maskinput);

                case {'fmri_mask_image'} % , 'fmri_data', 'statistic_image'} % others will not work with resample_to_space_defining... below

                    maskobj = maskinput;

                otherwise
                    error('region class constructor: unknown mask input type.')
            end
            clear maskinput

            % Either extract the image data in the space of the mask, or
            % vice versa
            %sample2mask = strmatch('sample2mask', varargin);

            if sample2mask
                % Read data in mask space; map images to mask
                % ------------------------------------------------

                if verbose, fprintf('Expanding image filenames if necessary\n'); end

                for i = 1:size(image_names, 1)

                    iinames{i} = expand_4d_filenames(image_names(i, :));
                end
                iinames = char(iinames{:});
                imgdat = zeros(length(maskobj.volInfo.wh_inmask), size(iinames, 1), 'single');

                if isempty(iinames)
                    disp('Images do not exist!'); disp(image_names); error('Exiting');
                end

                % Now extract the actual data from the mask
                if verbose
                    fprintf('Sampling %3.0f images to mask space: %04d', size(iinames, 1), 0)
                end

                for i = 1:size(iinames, 1)

                    if verbose, fprintf('\b\b\b\b%04d', i); end
                    idat = scn_map_image(iinames(i, :), maskobj);
                    imgdat(:, i) = idat(maskobj.volInfo.wh_inmask);

                end

                if verbose, fprintf('\n'); end

                [dd, ff, ee] = fileparts(maskobj.volInfo.fname);
                maskobj.space_defining_image_name = [ff ee];

                obj.volInfo = maskobj.volInfo;

            else % Not sample2mask, so resample mask to data image space

                % Read data in image space; map mask to images
                % ------------------------------------------------

                % resample to image space if necessary
                space_defining_image = deblank(image_names(1, :));
                maskobj = resample_to_image_space(maskobj, space_defining_image);

                % Now extract the actual data from the mask
                switch spm('Ver')

                    case {'SPM12','SPM8', 'SPM5'}
                        imgdat = iimg_get_data(maskobj.volInfo, image_names, 'single', verbosestr, 'noexpand');

                    case {'SPM2', 'SPM99'}
                        % legacy, for old SPM
                        imgdat = iimg_get_data(maskobj.volInfo, image_names, 'single', verbosestr);

                    otherwise
                        error('Unknown version of SPM! Update code, check path, etc.');
                end

                imgdat = imgdat';

            end % read data, depending on mask sampling

            % re-zip images if they were originally zipped
            % add .gz back to file names.
            if any(was_gzipped)

                image_names = re_zip_images(image_names, was_gzipped);

            end

            % add mask object to fmri_data object
            obj = create(obj, 'mask', maskobj, verbosestr);

            % append data
            obj = create(obj, 'dat', imgdat, 'image_names', image_names, verbosestr);
            clear imgdat

            % append description info
            [~, ff, ee] = fileparts(maskobj.volInfo.fname);
            mask_image_name = [ff ee];

            %obj = create(obj, 'dat_descrip', sprintf('Data from %s: %s', mask_image_name, obj.dat_descrip));
            obj = create(obj, 'mask_descrip', mask_image_name, verbosestr);

            obj.volInfo = maskobj.volInfo;

            obj.history(end+1) = {sprintf('Sampled to space of %s', maskobj.space_defining_image_name)};
            obj.history(end+1) = {['Masked with ' mask_image_name]};

            % Checks and fixes
            % -------------------------------------------------------------

            obj = run_checks_and_fixes(obj, verbosestr);

        end % constructor function

    end % methods

end

% -------------------------------------------------------------
% -------------------------------------------------------------

% Sub-functions

% -------------------------------------------------------------
% -------------------------------------------------------------

function obj = run_checks_and_fixes(obj, verbosestr)

obj = check_image_filenames(obj, verbosestr);

if isempty(obj.volInfo)
    error('obj.volInfo cannot be empty');
end

if isempty(obj.mask) && ~isempty(obj)  % isempty is an object method here
    % fix/create mask

    obj.mask.dat = any(obj.dat, 2);
    obj.mask.removed_voxels = obj.removed_voxels;
    obj.mask.volInfo_descrip = 'Generic mask built from .dat. Any voxel with a value in .dat is in-mask';

end

if issparse(obj.dat)
    obj.dat = full(obj.dat);
end

if ~isfield(obj.volInfo, 'cluster')
    obj = reparse_contiguous(obj);
end

if strcmp(verbosestr, 'verbose')

    % Check data bit rate
    % ---------------------------------------------------------------------

    databitrate = length(unique(obj.dat(:)));

    fprintf('Number of unique values in dataset: %d  Bit rate: %3.2f bits\n', databitrate, log2(databitrate));

    if databitrate < 2^10
        fprintf('Warning: Number of unique values in dataset is low, indicating possible restriction of bit rate. For comparison, Int16 has 65,536 unique values\n');
    end

end

end % function

% -------------------------------------------------------------

function image_names_out = re_zip_images(image_names, was_gzipped)
% in : char, out: char

image_names_out = cellstr(image_names);

for i = 1:size(image_names, 1)

    if was_gzipped(i)

        % Use system to remove unzipped version after zipping.
        % Will wait for input, and not overwrite, if images exist
        %         [status, result] = system(['gzip ' image_names(i, :)]);

        try
            gzip(deblank(image_names(i, :)));
            image_names_out{i} = [deblank(image_names_out{i}) '.gz'];

        catch
            warning('Error writing .gz images. Check permissions (or maybe using Git Annex?');

        end
    end

end

image_names_out = char(image_names_out{:});

end