text:
% Test code for irntv.m
%
% Legal:
%   irnTest.m is part of NUMIPAD (http://numipad.sf.net).
%
%   The NUMIPAD library is being developed under U.S. Government contract
%   W-7405-ENG-36 for Los Alamos National Laboratory.
%
% Authors
%   Paul Rodriguez    prodrig@pucp.edu.pe
%   Brendt Wohlberg   brendt@tmail.lanl.gov

%  example = 'l1deconv';
%  example = 'l2deconv';
%  example = 'l1denoise';
example = 'l2denoise';

%  Img = 'lena';
%  Img = 'peppers';
%  Img = 'goldhill';
Img = 'lena';
%  nmppath;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Input images        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch lower(Img)

  case{'lena'}
    Ig = double( imread('lena_gray_512.png') ) / 255;
    Ic = double( imread('lena_color_512.png') ) / 255;

  case{'peppers'}
    Ig = double( imread('peppers_gray.png') ) / 255;
    Ic = double( imread('peppers_color.png') ) / 255;

  case{'goldhill'}
    Ig = double( imread('goldhill_gray.png') ) / 255;
    Ic = double( imread('goldhill_color.png') ) / 255;

  case{'none'}
    disp(' ');
    disp('Select a test image (see code for an example)...');
    disp(' ');
    disp('Exiting irnTest code...');
    return;

  otherwise
    error('Not a valid image\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%
%     Normalize        %
%%%%%%%%%%%%%%%%%%%%%%%%

Normalize = @(x) (x - min(x(:)))/(max(x(:)) - min(x(:)));

%%%%%%%%%%%%%%%%%%%%%%%%
%       kernels        %
%%%%%%%%%%%%%%%%%%%%%%%%

kernel = fspecial('disk',3.2);

K = @(x) imfilter(x, kernel, 'symmetric','conv');
KT = @(x) K(x);
KC = {K, KT};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Blurred & noisy images        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch lower(example)

  case{'l1deconv'}
    IgBlur = K(Ig);
    Ig_01L1 = imnoise(IgBlur, 'salt & pepper', 0.1);
    Ig_03L1 = imnoise(IgBlur, 'salt & pepper', 0.3);

    IcBlur = K(Ic);
    Ic_01L1 = imnoise(IcBlur, 'salt & pepper', 0.1);
    Ic_03L1 = imnoise(IcBlur, 'salt & pepper', 0.3);

  case{'l1denoise'}

    Ig_01L1 = imnoise(Ig, 'salt & pepper', 0.1);
    Ig_03L1 = imnoise(Ig, 'salt & pepper', 0.3);

    Ic_01L1 = imnoise(Ic, 'salt & pepper', 0.1);
    Ic_03L1 = imnoise(Ic, 'salt & pepper', 0.3);

  case{'l2deconv'}

    IgBlur = K(Ig);
                                             %NOTE sigma^2 in imnoise
    Ig_01L2 = imnoise(IgBlur, 'gaussian', 0, 0.01*max(IgBlur(:)) );
    Ig_001L2 = imnoise(IgBlur,'gaussian', 0, 0.0001*max(IgBlur(:)) );

    IcBlur = K(Ic);
    Ic_01L2 = imnoise(IcBlur,'gaussian', 0, 0.01*max(IgBlur(:)) );
    Ic_001L2 = imnoise(IcBlur,'gaussian', 0, 0.0001*max(IgBlur(:)) );

  case{'l2denoise'}
    Ig_01L2 = imnoise(Ig,'gaussian', 0, 0.01*max(Ig(:)) );
    Ic_01L2 = imnoise(Ic,'gaussian', 0, 0.01*max(Ig(:)) );

    Ig_005L2 = imnoise(Ig,'gaussian', 0, 0.0025*max(Ig(:)) );
    Ic_005L2 = imnoise(Ic,'gaussian', 0, 0.0025*max(Ig(:)) );

end

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

if( strcmp(example,'l1deconv') || strcmp(example,'all') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L1 Deconvolved        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%noise level: 0.1

lambda = 0.45;

pars = irntvInputPars('l1tv');

%  pars.pcgtol_ini   = 1e-4;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 4;

%-----------------%
% -- Grayscale -- %

pars.U0           = Ig_01L1;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ig_01L1 = irntv(Ig_01L1, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ig_01L1) );
colormap gray; axis image; axis off;
title(sprintf('Deconvolved Image - Scalar IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_01L1), snr(Ig, Ig_01L1)));

%-------------%
% -- Color -- %

pars.U0           = Ic_01L1;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ic_01L1 = irntv(Ic_01L1, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ic_01L1) ); axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_01L1), snr(Ic, Ic_01L1)));

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

%noise level: 0.3

lambda = 0.95;

pars = irntvInputPars('l1tv');

%  pars.pcgtol_ini   = 1e-4;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 4;

%-----------------%
% -- Grayscale -- %

pars.U0           = Ig_03L1;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ig_03L1 = irntv(Ig_03L1, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ig_03L1) );
colormap gray; axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n', ...
               snr(Ig, IRN_Ig_03L1), snr(Ig, Ig_03L1)));

%-------------%
% -- Color -- %

pars.U0           = Ic_03L1;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ic_03L1 = irntv(Ic_03L1, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ic_03L1) ); axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB). \n', ...
               snr(Ic, IRN_Ic_03L1), snr(Ic, Ic_03L1)));

end % _END_ if( strcmp(example,'l1deconv') || strcmp(example,'all') )

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

if( strcmp(example,'l2deconv') || strcmp(example,'all') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L2 Deconvolved        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sigma (noise level) : 0.01

lambda  = 0.001;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini   = 1e-4;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 3;

%-----------------%
% -- Grayscale -- %

pars.U0           = Ig_001L2;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ig_001L2 = irntv(Ig_001L2, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ig_001L2) );
colormap gray; axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_001L2), snr(Ig, Ig_001L2)));

%-------------%
% -- Color -- %

pars.U0           = Ic_001L2;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ic_001L2 = irntv(Ic_001L2, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ic_001L2) ); axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_001L2), snr(Ic, Ic_001L2)));

% sigma (noise level) : 0.1

lambda  = 0.05;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini   = 1e-4;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.05;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.01;
pars.loops        = 3;

%-----------------%
% -- Grayscale -- %

pars.U0           = Ig_01L2;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ig_01L2 = irntv(Ig_01L2, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ig_01L2) );
colormap gray; axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_01L2), snr(Ig, Ig_01L2)));

%-------------%
% -- Color -- %

pars.U0           = Ic_01L2;  % initial solution

% >> Deconvolution via IRN <<
IRN_Ic_01L2 = irntv(Ic_01L2, KC, lambda, pars);

figure; imagesc( Normalize(IRN_Ic_01L2) );
colormap gray; axis image; axis off;
title(sprintf('Deconvolved Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_01L2), snr(Ic, Ic_01L2)));

end % _END_ if( strcmp(example,'l2deconv') || strcmp(example,'all') )

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

if( strcmp(example,'l1denoise') || strcmp(example,'all') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L1 Denoise        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%noise level: 10%

%-- IRN
lambda  = 1.1;

pars = irntvInputPars('l1tv');

pars.pcgtol_ini = 1e-4;
pars.epsf       = 1e-2;
pars.epsr       = 1e-4;
pars.loops      = 2;

%-----------------%
% -- Grayscale -- %

% >> Denoising via IRN <<
IRN_Ig_01L1 = irntv(Ig_01L1, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ig_01L1) );
colormap gray; axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_01L1), snr(Ig, Ig_01L1)));

%-------------%
% -- Color -- %

% >> Denoising via IRN <<
IRN_Ic_01L1 = irntv(Ic_01L1, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ic_01L1) );
colormap gray; axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_01L1), snr(Ic, Ic_01L1)));

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

%noise level: 30%

lambda  = 1.2;

pars = irntvInputPars('l1tv');

pars.pcgtol_ini   = 1e-4;
pars.adapt_epsR   = 1;
pars.epsR_cutoff  = 0.01;
pars.adapt_epsF   = 1;
pars.epsF_cutoff  = 0.05;
pars.loops        = 3;

%-----------------%
% -- Grayscale -- %

% >> Denoising via IRN <<
IRN_Ig_03L1 = irntv(Ig_03L1, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ig_03L1) );
colormap gray; axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_03L1), snr(Ig, Ig_03L1)));

%-------------%
% -- Color -- %

% >> Denoising via IRN <<
IRN_Ic_03L1 = irntv(Ic_03L1, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ic_03L1) ); axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_03L1), snr(Ic, Ic_03L1)));

end % _END_ if( strcmp(example,'l1denoise') || strcmp(example,'all') )

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

if( strcmp(example,'l2denoise') || strcmp(example,'all') )

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       L2 Denoise        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sigma (noise level) : 0.05

lambda  = 0.05;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini = 1e-4;
pars.epsf       = 1e-2;
pars.epsr       = 1e-5;
pars.loops      = 2;

%-----------------%
% -- Grayscale -- %

% >> Denoising via IRN <<
IRN_Ig_005L2 = irntv(Ig_005L2, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ig_005L2) );
colormap gray; axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_005L2), snr(Ig, Ig_005L2)));

%-------------%
% -- Color -- %

% >> Denoising via IRN <<
IRN_Ic_005L2 = irntv(Ic_005L2, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ic_005L2) ); axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_005L2), snr(Ic, Ic_005L2)));

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

% sigma (noise level) : 0.1

lambda  = 0.1;

pars = irntvInputPars('l2tv');

pars.pcgtol_ini = 1e-4;
pars.epsf       = 1e-2;
pars.epsr       = 1e-5;
pars.loops      = 2;

%-----------------%
% -- Grayscale -- %

% >> Denoising via IRN <<
IRN_Ig_01L2 = irntv(Ig_01L2, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ig_01L2) );
colormap gray; axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ig, IRN_Ig_01L2), snr(Ig, Ig_01L2)));

%-------------%
% -- Color -- %

% >> Denoising via IRN <<
IRN_Ic_01L2 = irntv(Ic_01L2, [], lambda, pars);

figure; imagesc( Normalize(IRN_Ic_01L2) ); axis image; axis off;
title(sprintf('Denoised Image - Vector IRN. \nSNR: %4.1fdB (noisy SNR: %4.1fdB).\n ', ...
               snr(Ic, IRN_Ic_01L2), snr(Ic, Ic_01L2)));

end % _END_ if( strcmp(example,'l2denoise') || strcmp(example,'all') )

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

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
% Based on https://se.mathworks.com/help/audio/ref/audiodevicereader-system-object.html
clear
clc
close all

% setup
samplingFreq = 8000; % Hz
recordingTime = 10; % s
segmentTime = 0.025; % s
segmentLength = round(segmentTime*samplingFreq);
deviceName = 'Default'; % use getAudioDevices(audioDeviceReader) to see available devices
micChannelNo = 1;

% initialise the recording object
% too see avai
AudioRecObj = audioDeviceReader(...
    'SampleRate', samplingFreq, ...
    'SamplesPerFrame', segmentLength, ...
    'NumChannels', length(micChannelNo), ...
    'Device', deviceName);

%
% AudioRecObj = audioDeviceReader(...
%     'SampleRate', samplingFreq, ...
%     'SamplesPerFrame', segmentLength, ...
%     'NumChannels', length(micChannelNo), ...
%     'Device', deviceName);
% Call setup to reduce the computational load of initialization in an
% audio stream loop.
AudioRecObj.setup();

% make a recording
nSegments = floor(recordingTime/segmentTime);
recordedData = nan(segmentLength, nSegments);
actFrameTimes = nan(1,nSegments);
for ii = 1:nSegments
%     tic
    [recordedSegment, numOverrun] = AudioRecObj();
%     toc
    % here we do some processing
    recordedData(:,ii) = recordedSegment(:);
    pause(0.01)
    % end of processing
    actFrameTimes(ii) = toc;
end
AudioRecObj.release();

%% plotting
figure(1)
plot((1:recordingTime*samplingFreq)'/samplingFreq, recordedData(:))
xlabel('time [s]')
ylabel('value [.]')
title('Recorded signal')

figure(2)
plot((1:nSegments)*segmentTime, cumsum(actFrameTimes),'.')
hold on
plot((1:nSegments)*segmentTime, (1:nSegments)*segmentTime, '--')
hold off
legend('Processed segment time', 'Real-time limit')
xlabel('sample recording time [s]')
ylabel('finished processing a segment [s]')
title('Processing time vs recorded time')

author: LimingShi


text:
function b = r8ge_mxv ( m, n, a, x )

%*****************************************************************************80
%
%% R8GE_MXV multiplies a R8GE matrix times a vector.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, real A(M,N), the R8GE matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(M), the product A * x.
%
  b(1:m) = a(1:m,1:n) * x(1:n)';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rad = complexPh2PositiveRad(complexPh);
%
% rad = complexPh2PositiveRad(complexPh);
%
% AUTHOR:  Wandell
% PURPOSE:
%   Convert the complex values values into radians that
% run from [0,2pi].  The Matlab angle function does the
% conversion into the range [-pi to pi].
%
rad = angle(complexPh);
l = find(rad < 0);
rad(l) = rad(l) + 2*pi;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% RMBASE - subtract basevector channel means from multi-epoch data matrix
%
% Usage:
%       >> [dataout] = rmbase(data); % remove whole-data channel means
%       >> [dataout datamean] = rmbase(data,frames,basevector);
%            % remove mean of basevector from each channel and epoch
% Inputs:
%   data       - data matrix (chans,frames*epochs) or (chans, frames, epochs);
%   frames     - data points per epoch {[]|0|default->data length}
%   basevector - vector of baseline frames per epoch
%                 Ex 1:128 {[]|0|default->whole epoch}
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 2-5-96

% Copyright (C) 2-5-96 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 07-30-97 converted to RMBASE -sm
% 09-30-97 fixed! -sm
% 05-10-01 caught empty input data -sm
% 06-03-01 added test for length-1 basevector, added [] defaults -sm
% 01-05-02 added check for negative basevector indices -Luca Finelli
% 01-25-02 reformated help & license -ad

function [dataout,datamean] = rmbase(data,frames,basevector)

	if nargin<3,
		basevector =0;
	end
    if isempty(basevector)
		basevector =0;
	end
    if length(basevector) == 1 && basevector(1) ~= 0
       fprintf('rmbase(): basevector should be a vector of frame indices.\n');
       return
    end

    if sum(basevector<0)~= 0
       fprintf('rmbase(): basevector should be 0 or a vector of positive frame indices.\n');
       return
    end

	if nargin < 2,
		frames = 0;
	end
    if isempty(frames)
		frames =0;
	end
	if nargin<1,
		help rmbase;
		fprintf('rmbase(): needs at least one argument.\n\n');
		return
	end
    if isempty(data)
		fprintf('rmbase(): input data is empty.\n\n');
		return
	end

    oridims = size(data);
	if ndims(data) == 3,
		data = reshape(data, size(data,1), size(data,2)*size(data,3));
	    reshape_flag=1;
	end

	[chans framestot]= size(data);
	if frames ==0,
		frames = framestot;
	end
    epochs = fix(framestot/frames);

	if length(basevector)>framestot,
		fprintf('rmbase(): length(basevector) > frames per epoch.\n\n');
		help rmbase;
		return
	end

    datamean = zeros(chans,epochs);
    % fprintf('removing epoch means for %d epochs\n',epochs);

    dataout = data;
    for e=1:epochs
        for c=1:chans
            if basevector(1)~=0,
                rmeans = nan_mean(double(data(c,(e-1)*frames+basevector)'));
            else
                rmeans = nan_mean(double(data(c,(e-1)*frames+1:e*frames)'));
                   %if e==1
                   %    fprintf('rmbase(): whole-data channel means removed. \n\n');
                   %end
            end
            datamean(c,e) = rmeans;
            dataout(c,(e-1)*frames+1:e*frames) = data(c,(e-1)*frames+1:e*frames) - rmeans;
        end
    end

    dataout = reshape(dataout, oridims);

% function out = nan_mean(in)
%
%     nans = find(isnan(in));
%     in(nans) = 0;
%     sums = sum(in);
%     nonnans = ones(size(in));
%     nonnans(nans) = 0;
%     nonnans = sum(nonnans);
%     nononnans = find(nonnans==0);
%     nonnans(nononnans) = 1;
%     out = sum(in)./nonnans;
%     out(nononnans) = NaN;
%

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function DICOM_fixMatlabBug
% When a DICOM I/O bug comes up, run this code once to fix.
%
% Details: Matlab's own DICOM toolbox is obviously written for Siemens
% scanner images. When running on GE images, it occassionally comes up with
% bugs. However, because the matlab source code is hidden in a private
% directory, we have to copy the file in VISTASOFT folder into that private
% directory, in order to overwrite the Matlab's problematic source code.
%
% Hence, this fixbug code may need to be expanded in the future if more
% bugs come up.
%
% Junjie Liu 2004/01/23
% 2004.04.08 Bob Dougherty: added check for pc- otherwise, this throws an
% error on linux due to file permissions.
% 2006.05.05 RFD: now politely backs up the original file. Is this patch
% still needed?
% 2006.05.22 RFD: disabled this fix- it's probably not needed and actually
% causes problems on newer matlabs.

matlabDir = fileparts(which('dicomread'));
if isempty(matlabDir);
    disp('Matlab version needs to be >6.5 to get DICOM toolbox');
    disp('In Linux please run matlab651 manually');
    error('dicomread function not found in Matlab. Your Matlab version is too old!');
end
% if(ispc)
%     matlabDir = fullfile(matlabDir,'private');
%
%     fixedDir = fileparts(which('DICOM_fixMatlabBug'));
%
%     disp('Overwritting your Matlab DICOM toolbox to fix possible bugs...');
%
%     fixedFile = fullfile(fixedDir,'dicom_read_attr.m');
%     matlabFile = fullfile(matlabDir,'dicom_read_attr.m');
%     [success message] = copyfile(matlabFile, fullfile(matlabDir,'dicom_read_attr_ORIGINAL.m'));
%     if(~success) warning(message); return; end;
%     [success message] = copyfile(fixedFile,matlabFile);
%     if ~success, warning(message); return; end;
%     disp('Overwritting Successful!');
% else
%     disp('Please ensure that your installation of matlab has been patched.');
% end
return;


text:
function [v,d,dref] = generate_doddbench(N)
% Generate DODDBENCH signal.
%
% Benchmark signal introduced in Dodd, T.J., Kadirkamanathan, V. and
% Harrison, R.F., "Function estimation in Hilbert space using sequential
% projections," Proc. of the IFAC Conf. on Intelligent Control Systems and
% Signal Processing, 113-118, 2003.
%
% Comment: copyright Cedric Richard, http://cedric-richard.fr/
%
% Input: N: number of data points
%
% Outputs: v: input sequence (2-dimensional sequence [v(:,1);v(:,2)])
%          d: noisy desired output (1-dimensional sequence)
%          dref: noise-free desired output
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

dref = zeros(1,N+2);
dref(1:2)=[0.1 0.1];

for t=3:N+2
    dref(t) = (0.8-0.5*exp(-dref(t-1)^2))*dref(t-1) - ...
        (0.3+0.9*exp(-dref(t-1)^2))*dref(t-2)+0.1*sin(pi*dref(t-1));
end
d = dref + 0.1*randn(1,N+2);
v = [d(1:N); d(2:N+1)]';

d(1:2)=[];
dref(1:2)=[];

repo: kafbox


text:
function [w] = simple_svm(X, y, lambda, opts)
% Simple function for solving an SVM with SDCA
% Used for local & global baselines

% Inputs
% X: input training data
% y: output training data
% lambda: regularization parameter

% Output
% w: the learned model

%% initialize
[n, d] = size(X);
w = zeros(d, 1);
alpha = zeros(n, 1);
primal_old = 0;

for iter=1:opts.max_sdca_iters
    % update coordinates cyclically
    for i=1:n
        % get current variables
        alpha_old = alpha(i);
        curr_x = X(i, :);
        curr_y = y(i);

        % calculate update
        grad = lambda * n * (1.0 - (curr_y * curr_x * w)) / (curr_x * curr_x') + (alpha_old * curr_y);

        % apply update
        alpha(i) = curr_y * max(0, min(1.0, grad));
        w = w + ((alpha(i) - alpha_old) * curr_x' * (1.0 / (lambda * n)));
    end

    % break if less than tol
    preds = y .* (X * w);
    primal_new = mean(max(0.0, 1.0 - preds)) + (lambda / 2.0) * (w' * w);
    if(abs(primal_old - primal_new) < opts.tol)
        break;
    end
    primal_old = primal_new;
end

end

author: gingsmith

repo: fmtl

save_path: github-repos/MATLAB/gingsmith-fmtl


text:
classdef RWMOP8 < PROBLEM
% <multi> <real> <constrained>
% Car side impact design problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 3;
            obj.D        = 7;
            obj.lower    = [0.5,0.45,0.5,0.5,0.875,0.4,0.4];
            obj.upper    = [1.5,1.35,1.5,1.5,2.625,1.2,1.2];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1);
            x2 = x(:,2);
            x3 = x(:,3);
            x4 = x(:,4);
            x5 = x(:,5);
            x6 = x(:,6);
            x7 = x(:,7);
            VMBP = 10.58-0.674.*x1.*x2-0.67275.*x2;
            VFD  = 16.45-0.489.*x3.*x7-0.843.*x5.*x6;
            % Objective function
            f(:,1) = 1.98+4.9.*x1.*6.67.*x2+6.98.*x3+4.01.*x4+1.78.*x5+1e-5.*x6+2.73.*x7;
            f(:,2) = 4.72-0.5.*x4-0.19.*x2.*x3;
            f(:,3) = 0.5.*(VMBP+VFD);
            % Constraints
            g(:,1) = -1+1.16-0.3717.*x2.*x4-0.0092928.*x3;
            g(:,2) = -0.32+0.261-0.0159.*x1.*x2-0.06486.*x1-0.019.*x2.*x7+0.0144.*x2.*x5+0.0154464.*x6;
            g(:,3) = -0.32+0.74-0.61.*x2-0.031296.*x3-0.031872.*x7+0.227.*x2.^2;
            g(:,4) = -0.32+0.214+0.00817.*x5-0.045195.*x1-0.0135168.*x1+0.03099.*x2.*x6-0.018.*x2.*x7+0.007176.*x3+0.023232.*x3-0.00364.*x5.*x6-0.018.*x2.^2;
            g(:,5) = -32+33.86+2.95.*x3-5.057.*x1.*x2-3.795.*x2-3.4431.*x7+1.45728;
            g(:,6) = -32+28.98+3.818.*x3-4.2.*x1.*x2+1.27296.*x6-2.68065.*x7;
            g(:,7) = -32+46.36-9.9.*x2-4.4505.*x1;
            g(:,8) = f(:,2)-4;
            g(:,9) = VMBP - 9.9;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
         %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [9.2596587e+01   4.0000000e+00   1.2699733e+01];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function q = affparaminv(p,q)
% function q = affparaminv(p[, q])
%
%    p(6,n) : [dx dy sc th sr phi]'
%    q(6,n) : [q(1) q(3) q(4); q(2) q(5) q(6)]

% Copyright (C) Jongwoo Lim and David Ross.  All rights reserved.

if (length(p) == 6)
  p = p(:);
end
if (nargin > 1)
  q = inv([p(3) p(4); p(5) p(6)]) * [q(1)-p(1) q(3:4); q(2)-p(2) q(5:6)];
else
  q = inv([p(3) p(4); p(5) p(6)]) * [-p(1) 1 0; -p(2) 0 1];
end
q = q([1,2,3,5,4,6]);

repo: Context-Aware-CF-Tracking


text:
function [A] = PoissonCstab1D();

% function [A] = PoissonCstab1D();
% Purpose: Set up symmetric Poisson matrix with estabilized central fluxes

Globals1D;
A = zeros(K*Np); g = zeros(K*Np,1);

% Build matrix -- one column at a time
for i=1:K*Np
    g(i) = 1.0;
    gmat = reshape(g,Np,K);
    Avec = PoissonCstabRHS1D(gmat);
    A(:,i) = reshape(Avec,K*Np,1);
    g(i)=0.0;
end
return

repo: Aero-matlab


text:
%% Plotting
%
% Below you find a list of plotting tools included in the MTEX toolbox. Of
% special importance is the command <saveFigure.html saveFigure> which
% allows saving plots in any kinds of image files.
%%
%

author: mtex-toolbox

repo: mtex

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


text:
function msg = er_displayParfiles(view);
%
% msg = er_displayParfiles(view);
%
% Pop up a message box with the scan group currently assigned
% to the present scan (if any), and the parfiles for each
% scan in the group. Returns the message text.
%
%
% ras, 02/06.
if notDefined('view'), view = getSelectedInplane;   end

mrGlobals;

dt = view.curDataType;
scan = view.curScan;
dtName = dataTYPES(dt).name;

if ~isfield(dataTYPES(dt).scanParams(scan), 'scanGroup') | ...
        ~isfield(dataTYPES(dt).scanParams(scan), 'parfile')
    msg = 'Parfiles and/or scan group not assigned for this scan.';
    myWarnDlg(msg);
    return
end

if ~isempty(dataTYPES(dt).scanParams(scan).scanGroup)
    [scans dt] = er_getScanGroup(view);
    scanGroup = sprintf('%s scans %s', dataTYPES(dt).name, num2str(scans));
else
    % just use cur scan / dt
    scans = scan;
    scanGroup = 'none';
end

msg = [sprintf('%s %s scan %i \n', view.name, dtName, scan) ...
       sprintf('Scan Group: %s \n\n', scanGroup) ...
       sprintf('Parfiles: \n')];
for s = scans
    msg = [msg sprintf('%i: %s \n', s, dataTYPES(dt).scanParams(s).parfile)];
end

msgbox(msg, 'Event-Related Info');

return


text:
function quiverSection(ori,varargin)
% plot orientations to ODF sections
%
% Input
%  ori - @orientation
%
% Options
%  sections - number of sections
%  points   - number of orientations to be plotted
%  all      - plot all orientations
%
% Flags
%  phi2      - phi2 sections (default)
%  phi1      - phi1 sections
%  sigma     - sigma = phi1 ~ phi2 sections
%  axisAngle - rotational angle sections
%
% See also
% vector3d/scatter saveFigure

if check_option(varargin,{'contour','contourf','smooth'})
  warning('Using the options contour, contourf or smooth in orientation sections plots is not recommented. Computing an ODF and plotting it is usually much better');
end

if ori.antipodal, varargin = [varargin,'antipodal']; end
oS = newODFSectionPlot(ori.CS,ori.SS,varargin{:});

oS.quiver(ori,varargin{:});


text:
function calpak_test ( )

%*****************************************************************************80
%
%% CALPAK_TEST tests the CALPAK library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the CALPAK library.\n' );

  calpak_test0005 ( );
  calpak_test0007 ( );

  calpak_test001 ( );
  calpak_test002 ( );
  calpak_test003 ( );
  calpak_test004 ( );
  calpak_test005 ( );
  calpak_test006 ( );
  calpak_test0065 ( );
  calpak_test007 ( );
  calpak_test0075 ( );
  calpak_test00755 ( );
  calpak_test00756 ( );
  calpak_test0076 ( );
  calpak_test0006 ( );
  calpak_test008 ( );
  calpak_test009 ( );

  calpak_test010 ( );
  calpak_test011 ( );
  calpak_test012 ( );
  calpak_test0125 ( );
  calpak_test013 ( );
  calpak_test014 ( );
  calpak_test015 ( );
  calpak_test016 ( );
  calpak_test017 ( );
  calpak_test0175 ( );
  calpak_test018 ( );
  calpak_test019 ( );

  calpak_test020 ( );

  calpak_test165 ( );
  calpak_test17 ( );
  calpak_test175 ( );
  calpak_test18 ( );
  calpak_test185 ( );
  calpak_test19 ( );
  calpak_test195 ( );

  calpak_test20 ( );
  calpak_test21 ( );
  calpak_test215 ( );
  calpak_test22 ( );
  calpak_test23 ( );
  calpak_test24 ( );
  calpak_test25 ( );
  calpak_test255 ( );
  calpak_test26 ( );
  calpak_test265 ( );
  calpak_test27 ( );
  calpak_test275 ( );
  calpak_test28 ( );
  calpak_test29 ( );

  calpak_test30 ( );
  calpak_test31 ( );
  calpak_test315 ( );
  calpak_test32 ( );
  calpak_test325 ( );
  calpak_test326 ( );
  calpak_test327 ( );
  calpak_test328 ( );
  calpak_test33 ( );
  calpak_test335 ( );
  calpak_test336 ( );
  calpak_test337 ( );
  calpak_test34 ( );
  calpak_test344 ( );
  calpak_test345 ( );
  calpak_test35 ( );
  calpak_test36 ( );
  calpak_test365 ( );
  calpak_test37 ( );
  calpak_test373 ( );
  calpak_test375 ( );
  calpak_test376 ( );
  calpak_test38 ( );
  calpak_test389 ( );
  calpak_test39 ( );
  calpak_test394 ( );
  calpak_test395 ( );

  calpak_test40 ( );
  calpak_test41 ( );
  calpak_test415 ( );
  calpak_test42 ( );
  calpak_test43 ( );
  calpak_test435 ( );
  calpak_test44 ( );
  calpak_test445 ( );
  calpak_test45 ( );
  calpak_test46 ( );
  calpak_test47 ( );
  calpak_test48 ( );
  calpak_test49 ( );
  calpak_test492 ( );
  calpak_test493 ( );
  calpak_test495 ( );

  calpak_test50 ( );
  calpak_test501 ( );
  calpak_test502 ( );
  calpak_test503 ( );
  calpak_test51 ( );
  calpak_test515 ( );
  calpak_test5153 ( )
  calpak_test51535 ( );
  calpak_test5154 ( );
  calpak_test5155 ( );
  calpak_test5156 ( );
  calpak_test52 ( );
  calpak_test525 ( );
  calpak_test53 ( );
  calpak_test535 ( );
  calpak_test54 ( );
  calpak_test555 ( );
  calpak_test56 ( );
  calpak_test565 ( );
  calpak_test566 ( );
  calpak_test57 ( );
  calpak_test58 ( );
  calpak_test585 ( );
  calpak_test59 ( );

  calpak_test60 ( );
  calpak_test605 ( );
  calpak_test61 ( );
  calpak_test615 ( );
  calpak_test616 ( );
  calpak_test62 ( );
  calpak_test621 ( );
  calpak_test622 ( );
  calpak_test623 ( );
  calpak_test624 ( );
  calpak_test63 ( );
  calpak_test635 ( );
  calpak_test636 ( );
  calpak_test64 ( );
  calpak_test65 ( );
  calpak_test66 ( );
  calpak_test67 ( );
  calpak_test675 ( );
  calpak_test68 ( );
  calpak_test685 ( );
  calpak_test686 ( );
  calpak_test687 ( );
  calpak_test688 ( );
  calpak_test689 ( );
  calpak_test69 ( );
  calpak_test695 ( );

  calpak_test70 ( );
  calpak_test71 ( );
  calpak_test72 ( );
  calpak_test73 ( );
  calpak_test74 ( );
  calpak_test75 ( );
  calpak_test76 ( );
  calpak_test77 ( );
  calpak_test775 ( );
  calpak_test78 ( );
  calpak_test79 ( );
  calpak_test795 ( );

  calpak_test80 ( );
  calpak_test805 ( );
  calpak_test81 ( );
  calpak_test82 ( );
  calpak_test83 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function h=semiaudplot(x,y,varargin)
%SEMIAUDPLOT  2D plot on auditory scale
%   Usage: h=semiaudplot(x,y);
%
%   `semiaudplot(x,y)` plots the data $(x,y)$ on an auditory scale. By
%   default the values of the x-axis will be shown on the Erb-scale.
%
%   `semiaudplot` takes the following parameters at the end of the line of input
%   arguments:
%
%     'x'       Make the x-axis use the auditory scale. This is the default.
%
%     'y'       Make the y-axis use the auditory scale.
%
%     'opts',c  Pass options stored in a cell array onto the plot
%               function.
%
%   In addition to these parameters, the auditory scale can be
%   specified. All scales supported by |freqtoaud| are supported. The default
%   is to use the erb-scale.
%
%   See also: freqtoaud

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import       = {'ltfattranslate','freqtoaud'};
definput.flags.plotdir= {'x','y'};
definput.keyvals.tick = [0,100,250,500,1000,2000,4000,8000,16000];
definput.keyvals.res  = 500;
definput.keyvals.opts = {};

[flags,kv]=ltfatarghelper({},definput,varargin);

n=500;
tickpos=freqtoaud(kv.tick,flags.audscale);

if flags.do_x
  xmin=min(x);
  xmax=max(x);
  audminmax=freqtoaud([xmin,xmax],flags.audscale);

  plotval=spline(x,y,audspace(xmin,xmax,n,flags.audscale));
  plot(linspace(audminmax(1),audminmax(2),n),plotval,kv.opts{:});
  set(gca,'XTick',tickpos);
  set(gca,'XTickLabel',num2str(kv.tick(:)));
  xlabel(sprintf('%s (Hz)',kv.frequency));

end;

if flags.do_y
  ymin=min(y);
  ymax=max(y);
  audminmax=freqtoaud([ymin,ymax],flags.audscale);

  plot(x,freqtoerb(y),kv.opts{:});
  set(gca,'YTick',tickpos);
  set(gca,'YTickLabel',num2str(tick(:)));

  ylabel(sprintf('%s (Hz)',kv.frequency));
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
clear all; close all; clc

addpath('apm')

s = 'http://byu.apmonitor.com';
b = 'mhe';

% Connect to Arduino
tclab;

% Run time in minutes
run_time = 5.0;

% Number of cycles (1 cycle per 2 seconds)
loops = round(30*run_time);

% Temperature (degC)
T1 = ones(1,loops) * T1C(); % measured T
T1mhe = ones(1,loops) * T1C(); % measured T
Kp = ones(1,loops) * 0.3598;
tau = ones(1,loops) * 47.73;
TC_ss = ones(1,loops) * 23;
time = zeros(1,loops);

% milli-volts input
Q1 = zeros(1,loops);
Q2 = zeros(1,loops);
Q1(6:end) = 100.0;
Q1(50:end) = 20.0;
Q1(100:end) = 80.0;
Q1(150:end) = 10.0;
Q1(200:end) = 95.0;
Q1(250:end) = 0.0;

% time
tm = zeros(1,loops);

% moving horizon estimation
mhe_init();

start_time = clock;
prev_time = start_time;

% dynamic plot (note: subplots needs to be declared here first)
figure(1)
subplot(3,1,1)
hold on, grid on
anexp1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anpred1 = animatedline('LineStyle','--','Color', 'r','LineWidth', 2);
ylabel('Temperature \circC')
legend('T_1 Measured', 'T_1 Predicted', ...
    'Location', 'northwest')
title('Temperature Estimation')
subplot(3,1,2)
hold on, grid on
anQ1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anQ2 = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ylabel('Power Level Q (%)')
legend('Q_1', 'Q_2', 'Location', 'northwest')
subplot(3,1,3)
hold on, grid on
anK = animatedline('LineStyle','-', 'Color', 'r', 'LineWidth', 2);
antau = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ylabel('Parameters')
legend('K x 100', 'tau', 'Location', 'northwest')
xlabel('Time (sec)')

for ii = 1:loops
    % adjust power level
    h1(Q1(ii));
    h2(Q2(ii));

    % Pause Sleep time
    pause_max = 2.0;
    pause_time = pause_max - etime(clock,prev_time);
    if pause_time >= 0.0
        pause(pause_time - 0.01)
    else
        pause(0.01)
    end

    % Record time and change in time
    t = clock;
    dt = etime(t,prev_time);
    if ii>=2
        time(ii) = time(ii-1) + dt;
    end
    prev_time = t;

    % read and record from temperature controller
    T1(ii) = T1C();
    T2(ii) = T2C();

    % non-linear energy balance
    jj = ii+1;
    params = mhe(T1(ii),Q1(ii));
    Kp(jj) = params(1);
    tau(jj) = params(2);
    T1mhe(jj) = params(4);

    % plot
    addpoints(anexp1,time(ii),T1(ii))
    addpoints(anpred1,time(ii),T1mhe(ii))
    addpoints(anQ1,time(ii),Q1(ii))
    addpoints(anQ2,time(ii),Q2(ii))
    addpoints(anK,time(ii),100*Kp(ii))
    addpoints(antau,time(ii),tau(ii))
    drawnow

    if ii==10
        apm_web(s,b);
    end
end

h1(0);
h2(0);
disp('Heaters off')
% turn off heater but keep LED on if T > 50
if (T1C() || T2C()) > 50
    led(1)
    disp(['Warning, heater temperature 1 =', num2str(T1C())])
    disp(['Warning, heater temperature 2 =', num2str(T2C())])
else
    led(0)
end

% save txt file with data
data = [time',Q1',Q2',T1',T2'];
csvwrite('data.txt',data);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returs the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

for k=1:K
    % use logical arrays for indexing
    % see http://www.mathworks.com/help/matlab/math/matrix-indexing.html#bq7egb6-1
    indexes = idx == k;
    centroids(k, :) = mean(X(indexes, :));
end;

% =============================================================

end

author: zsiciarz

repo: ml-coursera


text:
% Copyright 2014-2015 The MathWorks, Inc.
%% Load sample frames
load sampleFrames.mat
subplot(1,3,1)
imshow(vidFrame1)

%% Threshold image
% Convert RGB image to chosen color space
I = rgb2hsv(vidFrame1);

% Define thresholds for channel 1 based on histogram settings
channel1Min = 0.379;
channel1Max = 0.496;

% Define thresholds for channel 2 based on histogram settings
channel2Min = 0.436;
channel2Max = 1.000;

% Define thresholds for channel 3 based on histogram settings
channel3Min = 0.000;
channel3Max = 1.000;

% Create mask based on chosen histogram thresholds
BW = (I(:,:,1) >= channel1Min ) & (I(:,:,1) <= channel1Max) & ...
    (I(:,:,2) >= channel2Min ) & (I(:,:,2) <= channel2Max) & ...
    (I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max);

subplot(1,3,2)
imshow(BW)

%% Remove disturbances
diskElem = strel('disk',3);
Ibwopen = imopen(BW,diskElem);
subplot(1,3,3)
imshow(Ibwopen)

%% Blob Analysis
hBlobAnalysis = vision.BlobAnalysis('MinimumBlobArea',200,...
    'MaximumBlobArea',5000);
[objArea,objCentroid,bboxOut] = step(hBlobAnalysis,Ibwopen);

%% Annotate image
Ishape = insertShape(vidFrame1,'rectangle',bboxOut,'Linewidth',4);
figure
subplot(1,2,1)
imshow(Ishape)

numObj = numel(objArea);
hTextIns = vision.TextInserter('%d','Location',[20 20],'Color',...
    [255 255 0],'FontSize',30);
Itext = step(hTextIns,Ishape,int32(numObj));
subplot(1,2,2)
imshow(Itext)

%% Clean up
release(hBlobAnalysis)
release(hTextIns)

author: TheAlgorithms

repo: MATLAB-Octave


text:
function PPI = spm_peb_ppi(varargin)
% Bold deconvolution to create physio- or psycho-physiologic interactions
% FORMAT PPI = spm_peb_ppi(SPMname,ppiflag,VOI,Uu,ppiname,showGraphics)
%
% SPM          - Structure containing generic details about the analysis or
%                the fully qualified filename of such a structure.
% ppiflag      - Type of analysis. Must be one of:
%                  'simple deconvolution'          or 'sd'
%                  'psychophysiologic interaction' or 'ppi'
%                  'physiophysiologic interaction' or 'phipi'
% VOI          - Structure containing details about a VOI (as produced by
%                spm_regions) or the fully qualified filename of such a
%                structure. If a structure, then VOI should be of size 1x1
%                in the case of simple deconvolution, and psychophysiologic
%                interactions) or 1x2, in the case of physiophysiologic
%                interactions. If a file name it should be 1xN or 2xN.
% Uu           - Matrix of input variables and contrast weights. This is an
%                [n x 3] matrix. The first column indexes SPM.Sess.U(i). The
%                second column indexes the name of the input or cause, see
%                SPM.Sess.U(i).name{j}. The third column is the contrast
%                weight. Unless there are parametric effects the second
%                column will generally be a 1.
% ppiname      - Basename of the PPI file to save. The saved file will be:
%                <PATH_TO_SPM.MAT>/PPI_<ppiname>.mat
% showGraphics - empty or 1 = yes, 0 = no.
%
%
% PPI.ppi      - (PSY*xn  or xn1*xn2) convolved with the HRF
% PPI.Y        - Original BOLD eigenvariate. Use as covariate of no interest
% PPI.P        - PSY convolved with HRF for psychophysiologic interactions,
%                or in the case of physiophysologic interactions contains
%                the eigenvariate of the second region.
% PPI.name     - Name of PPI
% PPI.xY       - Original VOI information
% PPI.xn       - Deconvolved neural signal(s)
% PPI.psy.u    - Psychological variable or input function (PPIs only)
% PPI.psy.w    - Contrast weights for psychological variable (PPIs only)
% PPI.psy.name - Names of psychological conditions (PPIs only)
%__________________________________________________________________________
%
% This routine is effectively a hemodynamic deconvolution using full priors
% and EM to deconvolve the HRF from a hemodynamic time series to give a
% neuronal time series [that can be found in PPI.xn].  This deconvolution
% conforms to Wiener filtering. The neuronal process is then used to form
% PPIs. See help text within function for more details.
%__________________________________________________________________________
% Copyright (C) 2002-2014 Wellcome Trust Centre for Neuroimaging

% Darren Gitelman
% $Id: spm_peb_ppi.m 6556 2015-09-15 15:42:04Z guillaume $

% SETTING UP A PPI THAT ACCOUNTS FOR THE HRF
% =========================================================================
% PPI's were initially conceived as a means of identifying regions whose
% reponses can be explained in terms of an interaction between activity in
% a specified source (the physiological factor) and some experimental
% effect (the psychological factor). However, a problem in setting up PPI's
% is that in order to derive a proper estimate of the interaction between
% a psychological variable (P) and measured hemodynamic signal (x), one
% cannot simply convolve the psychological variable with the hrf (HRF) and
% multiply by the signal. Thus:
%
%                  conv(P,HRF).* x ~= conv((P.*xn),HRF)
%
% P   = psychological variable
% HRF = hemodynamic response function
% xn  = underlying neural signal which in fMRI is convolved with the hrf to
%       give the signal one measures -- x.
% x   = measured fmri signal
%
% It is actually the right hand side of the equation one wants.
% Thus one has to work backwards, in a sense, and deconvolve the hrf
% from x to get xn. This can then be multiplied by P and the resulting
% vector (or matrix) reconvolved with the hrf.
%
% This algorithm uses a least squares strategy to solve for xn.
%
% The source's hemodynamics are x = HRF*xn;
%
% Using the constraint that xn should have a uniform spectral density
% we can expand x in terms of a discrete cosine set (xb)
%
%      xn  = xb*B
%       B  = parameter estimate
%
% The estimator of x is then
%
%       x  = HRF(k,:)*xn
%       x  = HRF(k,:) * xb * B
%
% This accounts for different time resolutions between our hemodynamic
% signal and the discrete representation of the psychological variable. In
% this case k is a vector representing the time resolution of the scans.
%
% Conditional estimates of B allow for priors that ensure uniform variance
% over frequencies.
%
% PPI STATISTICAL MODEL
% =========================================================================
% Once the PPI.ppi interaction term has been calculated a new GLM must be
% setup to search for the interaction effects across the brain. This is
% done using a standard, first level, fMRI model, which must include 3
% covariates, PPI.ppi (interaction), PPI.Y (main effect: source region bold
% signal) and PPI.P (main effect: "psychological" condition), plus any
% nuisance regressors according to the particular design.
%
% NB: Designs that include only the interaction term without the main
% effects are not proper as inferences on the interaction will include a
% mixture of both main and interaction effects.
%
% Once the model has been setup and run, a contrast of [1 0 0] over the
% PPI.ppi, PPI.Y and PPI.P columns respectively, will show regions with a
% positive relationship to the interaction term, discounting any main
% effects. Negative regressions can be examined with [-1 0 0]. A PPI random
% effects analysis would involve taking the con* files from the [1 0 0]
% t-contrast for each subject and forwarding them to a second level
% analysis.

% Shortcut for PPI display
%--------------------------------------------------------------------------
if nargin && ischar(varargin{1}) && strcmpi(varargin{1},'display')
    if nargin > 1
        PPI = varargin{2};
    else
        [PPI, sts] = spm_select(1,'^PPI.*\.mat$');
        if ~sts, return; end
    end
    if ischar(PPI), load(PPI); end
    display_ppi(PPI);
    return;
end

% Set up the graphical interface
%--------------------------------------------------------------------------
Finter = spm_figure('FindWin','Interactive');
if ~nargin, Finter = spm_figure('GetWin','Interactive'); end
header = get(Finter,'Name');
spm_clf(Finter); set(Finter,'name','PPI Setup');

% Check inputs
%--------------------------------------------------------------------------
if nargin && isstruct(varargin{1})
    SPM = varargin{1};
    try, swd = SPM.pwd; catch, swd = pwd; end
else
    try
        P            = varargin{1};
    catch
        [P, sts]     = spm_select(1,'^SPM\.mat$','Select SPM.mat');
        if ~sts, PPI = []; return; end
    end
    swd = spm_file(P,'fpath');
    load(fullfile(swd,'SPM.mat'));
end
SPM.swd = spm_file(swd,'cpath');
cwd     = pwd;
cd(SPM.swd)

% Ask whether to perform physiophysiologic or psychophysiologic interactions
%--------------------------------------------------------------------------
try
    ppiflag = varargin{2};
catch
    ppiflag = {'simple deconvolution',...
               'psychophysiologic interaction',...
               'physiophysiologic interaction'};
    i       = spm_input('Analysis type?',1,'m',ppiflag);
    ppiflag = ppiflag{i};
end

switch lower(ppiflag)

    case  {'simple deconvolution','sd'}
    %======================================================================
    if nargin>2 && isstruct(varargin{3})
        p.xY = varargin{3};
    else
        try
            VOI = varargin{3};
            p   = load(deblank(VOI(1,:)),'xY');
        catch
            spm_input('physiological variable:...  ',2,'d');
            [VOI,sts] = spm_select(1,'^VOI.*\.mat$',{'select VOI'});
            if ~sts, PPI = []; return; end
            p   = load(VOI,'xY');
        end
    end
    xY(1) = p.xY;
    Sess  = SPM.Sess(xY(1).Sess);

    case  {'physiophysiologic interaction','phipi'}
    %======================================================================
    if nargin>2 && isstruct(varargin{3})
        xY = varargin{3};
        xY = xY(:)';
        if numel(xY) ~= 2
            error('Must include 2 VOI structures for physiophysiologic interactions')
        end
    else
        try
            VOI = varargin{3};
            if size(VOI,1) ~= 2
                error('Must include 2 VOI filenames for physiophysiologic interactions')
            end
            for i = 1:2
                p     = load(deblank(VOI(i,:)),'xY');
                xY(i) = p.xY;
            end
        catch
            spm_input('physiological variables:...  ',2,'d');
            [VOI,sts] = spm_select(2,'^VOI.*\.mat$',{'select VOIs'});
            if ~sts, PPI = []; return; end
            for i = 1:2
                p     = load(deblank(VOI(i,:)),'xY');
                xY(i) = p.xY;
            end
        end
    end
    Sess = SPM.Sess(xY(1).Sess);

    case  {'psychophysiologic interaction','ppi'}
    %======================================================================
    if nargin>2 && isstruct(varargin{3})
        p.xY = varargin{3};
    else
        try
            VOI = varargin{3};
            p   = load(deblank(VOI(1,:)),'xY');
        catch
            spm_input('physiological variable:...  ',2,'d');
            [VOI,sts] = spm_select(1,'^VOI.*\.mat$',{'select VOI'});
            if ~sts, PPI = []; return; end
            p   = load(VOI,'xY');
        end
    end
    xY(1) = p.xY;
    Sess  = SPM.Sess(xY(1).Sess);

    % get 'causes' or inputs U
    %----------------------------------------------------------------------
    U.name = {};
    U.u    = [];
    U.w    = [];
    try
        Uu = varargin{4};
        for i = 1:size(Uu,1)
            U.u           = [U.u Sess.U(Uu(i,1)).u(33:end,Uu(i,2))];
            U.name{end+1} = Sess.U(Uu(i,1)).name{Uu(i,2)};
            U.w           = [U.w Uu(i,3)];
        end
    catch
        spm_input('Psychological variable:...  ',2,'d');
        u      = length(Sess.U);
        for  i = 1:u
            for  j = 1:length(Sess.U(i).name)
                str   = ['include ' Sess.U(i).name{j} '?'];
                if spm_input(str,3,'y/n',[1 0])
                    str             = 'Contrast weight';
                    tmpw            = spm_input(str,4,'e',[],1);
                    % if tmpw==0 then don't include the column in the
                    % design. This takes care of the possibility that the
                    % user would select to include the column but then give
                    % it a 0 weight.
                    %------------------------------------------------------
                    if tmpw ~= 0
                        U.w             = [U.w tmpw];
                        U.u             = [U.u Sess.U(i).u(33:end,j)];
                        U.name{end + 1} = Sess.U(i).name{j};
                    end
                end
            end
        end
    end

    otherwise
    %======================================================================
    error('Unknown type of analysis');

end % (switch setup)

% Name of PPI file to be saved
%--------------------------------------------------------------------------
try
    PPI.name = varargin{5};
catch
    PPI.name = spm_input('Name of PPI',3,'s','PPI');
end

% Check if Graphical output should be shown
%--------------------------------------------------------------------------
try
    showGraphics = varargin{6};
catch
    showGraphics = 1;
end

% Setup variables
%--------------------------------------------------------------------------
RT      = SPM.xY.RT;
dt      = SPM.xBF.dt;
NT      = round(RT/dt);
fMRI_T0 = SPM.xBF.T0;
N       = length(xY(1).u);
k       = 1:NT:N*NT;                       % microtime to scan time indices

% Setup other output variables
%--------------------------------------------------------------------------
PPI.xY = xY;
PPI.RT = RT;
PPI.dt = dt;

% Create basis functions and hrf in scan time and microtime
%--------------------------------------------------------------------------
hrf = spm_hrf(dt);

% Create convolved explanatory {Hxb} variables in scan time
%--------------------------------------------------------------------------
xb  = spm_dctmtx(N*NT + 128,N);
Hxb = zeros(N,N);
for i = 1:N
    Hx       = conv(xb(:,i),hrf);
    Hxb(:,i) = Hx(k + 128);
end
xb = xb(129:end,:);

% Get confounds (in scan time) and constant term
%--------------------------------------------------------------------------
X0 = xY(1).X0;
M  = size(X0,2);

% Get response variable
%--------------------------------------------------------------------------
for i = 1:size(xY,2)
    Y(:,i) = xY(i).u;
end

% Remove confounds and save Y in ouput structure
%--------------------------------------------------------------------------
Yc    = Y - X0*inv(X0'*X0)*X0'*Y;
PPI.Y = Yc(:,1);
if size(Y,2) == 2
    PPI.P = Yc(:,2);
end

% Specify covariance components; assume neuronal response is white
% treating confounds as fixed effects
%--------------------------------------------------------------------------
Q = speye(N,N)*N/trace(Hxb'*Hxb);
Q = blkdiag(Q, speye(M,M)*1e6  );

% Get whitening matrix (NB: confounds have already been whitened)
%--------------------------------------------------------------------------
W = SPM.xX.W(Sess.row,Sess.row);

% Create structure for spm_PEB
%--------------------------------------------------------------------------
clear P
P{1}.X = [W*Hxb X0];        % Design matrix for lowest level
P{1}.C = speye(N,N)/4;      % i.i.d assumptions
P{2}.X = sparse(N + M,1);   % Design matrix for parameters (0's)
P{2}.C = Q;

switch ppiflag

    case  {'simple deconvolution','sd'}
    %======================================================================
    C  = spm_PEB(Y,P);
    xn = xb*C{2}.E(1:N);
    xn = spm_detrend(xn);

    % Save variables (NOTE: xn is in microtime and does not account for
    % slice timing shifts). To convert to BOLD signal convolve with a hrf.
    % Use a microtime to scan time index to convert to scan time: e.g.,
    % k = 1:NT:N*NT; where NT = number of bins per TR = TR/dt or SPM.xBF.T
    % and N = number of scans in the session. Finally account for slice
    % timing effects by shifting the index accordingly.
    %----------------------------------------------------------------------
    PPI.xn = xn;

    case  {'physiophysiologic interaction','phipi'}
    %======================================================================
    C    = spm_PEB(Y(:,1),P);
    xn1  = xb*C{2}.E(1:N);
    C    = spm_PEB(Y(:,2),P);
    xn2  = xb*C{2}.E(1:N);
    xn1  = spm_detrend(xn1);
    xn2  = spm_detrend(xn2);
    xnxn = xn1.*xn2;

    % Convolve, convert to scan time, and account for slice timing shift
    %----------------------------------------------------------------------
    ppi = conv(xnxn,hrf);
    ppi = ppi((k-1) + fMRI_T0);

    % Save variables
    %----------------------------------------------------------------------
    PPI.xn  = [xn1 xn2];
    PPI.ppi = spm_detrend(ppi);

    case  {'psychophysiologic interaction','ppi'}
    %======================================================================

    % COMPUTE PSYCHOPHYSIOLOGIC INTERACTIONS
    % use basis set in microtime
    %----------------------------------------------------------------------
    % get parameter estimates and neural signal; beta (C) is in scan time
    % This clever trick allows us to compute the betas in scan time which
    % is much quicker than with the large microtime vectors. Then the betas
    % are applied to a microtime basis set generating the correct neural
    % activity to convolve with the psychological variable in microtime
    %----------------------------------------------------------------------
    C  = spm_PEB(Y,P);
    xn = xb*C{2}.E(1:N);
    xn = spm_detrend(xn);

    % Setup psychological variable from inputs and contrast weights
    %----------------------------------------------------------------------
    PSY = zeros(N*NT,1);
    for i = 1:size(U.u,2)
        PSY = PSY + full(U.u(:,i) * U.w(i));
    end
    PSY = spm_detrend(PSY);

    % Multiply psychological variable by neural signal
    %----------------------------------------------------------------------
    PSYxn = PSY.*xn;

    % Convolve, convert to scan time, and account for slice timing shift
    %----------------------------------------------------------------------
    ppi = conv(PSYxn,hrf);
    ppi = ppi((k-1) + fMRI_T0);

    % Convolve psych effect, convert to scan time, and account for slice
    % timing shift
    %----------------------------------------------------------------------
    PSYHRF = conv(PSY,hrf);
    PSYHRF = PSYHRF((k-1) + fMRI_T0);

    % Save psychological variables
    %----------------------------------------------------------------------
    PPI.psy = U;
    PPI.P   = PSYHRF;
    PPI.xn  = xn;
    PPI.ppi = spm_detrend(ppi);

end % (switch)

% Display
%--------------------------------------------------------------------------
if showGraphics
    display_ppi(PPI);
end

% Save
%--------------------------------------------------------------------------
str    = ['PPI_' PPI.name '.mat'];

save(fullfile(SPM.swd,str),'PPI', spm_get_defaults('mat.format'))

cmd = 'spm_peb_ppi(''display'',''%s'')';
fprintf('   PPI saved as %s\n',spm_file(fullfile(SPM.swd,str),'link',cmd));

% Clean up
%--------------------------------------------------------------------------
set(Finter,'name',header);
cd(cwd);

%==========================================================================
% function display_ppi(PPI)
%==========================================================================
function display_ppi(PPI)

Fgraph = spm_figure('GetWin','PPI');
FS     = spm('FontSizes');
spm_clf(Fgraph);

try
    RT = PPI.RT;
catch
    RT = spm_get_defaults('stats.fmri.t')*PPI.dt; % backward compatibility
end
N  = length(PPI.xY(1).u);
NT = round(RT/PPI.dt);
t  = PPI.RT*[1:N];
T  = PPI.dt*[1:(N*NT)];

%-Simple deconvolution
%--------------------------------------------------------------------------
if numel(PPI.xY) == 1 && ~isfield(PPI,'psy')

    ax = subplot(1,1,1);
    plot(t,PPI.Y,T,PPI.xn)
    title('hemodynamic and neuronal responses')
    xlabel(['time (secs)'])
    axis tight square
    grid on
    legend('BOLD','neuronal')

    str = sprintf('Simple Deconvolution: %s\n',PPI.name);
    str = [str sprintf('VOI file: %s',PPI.xY.name)];

%-Physiophysiologic interaction
%--------------------------------------------------------------------------
elseif numel(PPI.xY) == 2

    ax = subplot(2,1,1);
    plot(t,PPI.ppi)
    title('PPI')
    xlabel(['time (secs)'])
    axis tight square
    grid on

    subplot(2,2,3)
    plot(t,PPI.Y(:,1),T,PPI.xn(:,1))
    title('hemodynamic and neuronal responses (1st)')
    xlabel(['time (secs)'])
    axis tight square
    grid on
    legend('BOLD','neuronal')

    subplot(2,2,4)
    plot(t,PPI.P,T,PPI.xn(:,2))
    title('hemodynamic and neuronal responses (2nd)')
    xlabel(['time (secs)'])
    axis tight square
    grid on
    legend('BOLD','neuronal')

    str = sprintf('Physiophysiologic Interaction: %s\n',PPI.name);
    str = [str, sprintf('VOI File 1: %s\n',PPI.xY(1).name)];
    str = [str, sprintf('VOI File 2: %s',PPI.xY(2).name)];

%-Psychophysiologic interaction
%--------------------------------------------------------------------------
elseif isfield(PPI,'psy')

    PSY = zeros(N*NT,1);
    for i = 1:size(PPI.psy.u,2)
        PSY = PSY + full(PPI.psy.u(:,i)*PPI.psy.w(:,i));
    end

    ax = subplot(2,1,1);
    plot(t,PPI.Y,T,PPI.xn(:,1))
    title('hemodynamic and neuronal responses')
    xlabel(['time (secs)'])
    axis tight square
    grid on
    legend('BOLD','neuronal')

    subplot(2,2,3)
    plot(T,PSY,'LineStyle','--','Color',[0 .65 0]);
    hold on
    plot(t,PPI.P,'LineStyle','-','LineWidth',1,'Color','b');
    hold off
    title('[convolved] psych. variable')
    xlabel(['time (secs)'])
    axis tight square
    grid on

    subplot(2,2,4)
    plot(t,PPI.ppi)
    title('PPI')
    xlabel(['time (secs)'])
    axis tight square
    grid on

    str = sprintf('Psychophysiologic Interaction: %s\n',PPI.name);
    str = [str, sprintf('VOI File: %s\n',PPI.xY(1).name)];
    str = [str, sprintf('Factors: ')];
    for i = 1:numel(PPI.psy.name)
        str = [str, sprintf('%s [%0.0f]',PPI.psy.name{i},PPI.psy.w(i))];
        if i < numel(PPI.psy.name)
            str = [str, sprintf('; ')];
        end
    end

else
    error('Unknown PPI type.');
end

hAx = axes('Position',[0 0.90 1 0.1],...
    'DefaultTextFontSize',FS(8),...
    'DefaultTextInterpreter','None',...
    'DefaultTextVerticalAlignment','Baseline',...
    'Parent',Fgraph,...
    'Units','points',...
    'Visible','off');
AxPos = get(hAx,'Position');
set(hAx,'XLim',[0,AxPos(3)]); set(hAx,'YLim',AxPos(2)+[0,AxPos(4)])
dy    = FS(9);
h     = text(dy,floor(AxPos(2)+AxPos(4))-2*dy,str);
set(h,'Parent',hAx);

author: spm


text:
%% Image Pyramids
%
% In this tutorial:
%
% * We will learn about Image Pyramids
% * We will use Image pyramids to create a new fruit, "Orapple"
% * We will see these functions: |cv.pyrUp|, |cv.pyrDown|, |cv.buildPyramid|
%
% Sources:
%
% * <https://docs.opencv.org/3.1.0/dc/dff/tutorial_py_pyramids.html>
%

%% Theory
% Normally, we used to work with an image of constant size. But on some
% occasions, we need to work with (the same) images in different resolutions.
% For example, while searching for something in an image, like face, we are
% not sure at what size the object will be present in said image. In that
% case, we will need to create a set of the same image with different
% resolutions and search for object in all of them. These set of images with
% different resolutions are called *Image Pyramids* (because when they are
% kept in a stack with the highest resolution image at the bottom and the
% lowest resolution image at top, it looks like a pyramid).
%
% There are two kinds of Image Pyramids:
%
% # Gaussian Pyramid
% # Laplacian Pyramids
%
% Higher level (Low resolution) in a Gaussian Pyramid is formed by removing
% consecutive rows and columns in Lower level (higher resolution) image.
% Then each pixel in higher level is formed by the contribution from 5 pixels
% in underlying level with gaussian weights. By doing so, a |MxN| image
% becomes |M/2xN/2| image. So area reduces to one-fourth of original area. It
% is called an Octave. The same pattern continues as we go upper in pyramid
% (ie, resolution decreases). Similarly while expanding, area becomes 4 times
% in each level. We can find Gaussian pyramids using |cv.pyrDown| and
% |cv.pyrUp| functions.
%
% Below is the 4 levels in an image pyramid:
%

% compute Gaussian pyramid
img = cv.imread(fullfile(mexopencv.root(),'test','fruits.jpg'), 'ReduceScale',2);
p = cv.buildPyramid(img, 'MaxLevel',4);

% combine all levels in one image
[rows,cols,~] = cellfun(@size, p);
dst = zeros([rows(1) cols(1)+cols(2) size(img,3)], class(img));
r_start = [0 0 cumsum(rows(2:end-1))] + 1;
r_end = [rows(1) cumsum(rows(2:end))];
c_start = [0 repmat(cols(1),1,numel(p)-1)] + 1;
c_end = [cols(1) cols(1)+cols(2:end)];
for i=1:numel(p)
    dst(r_start(i):r_end(i), c_start(i):c_end(i), :) = p{i};
end
figure, imshow(dst)

%%
% Now you can go down the image pyramid with |cv.pyrUp| function.
%
% Remember, |higher_reso2| is not equal to |higher_reso|, because once you
% decrease the resolution, you loose the information.
%

higher_reso = img;
lower_reso = cv.pyrDown(higher_reso);
higher_reso2 = cv.pyrUp(lower_reso);
figure, imshow(cat(2, higher_reso, higher_reso2))

%%
% Laplacian Pyramids are formed from the Gaussian Pyramids. There is no
% exclusive function for that. Laplacian pyramid images are like edge images
% only. Most of its elements are zeros. They are used in image compression.
% A level in Laplacian Pyramid is formed by the difference between that level
% in Gaussian Pyramid and expanded version of its upper level in Gaussian
% Pyramid. The three levels of a Laplacian level will look like below
% (contrast is adjusted to enhance the contents):
%

% compute Gaussian pyramid
img = imread(fullfile(mexopencv.root(),'test','butterfly.jpg'));
p = cv.Blender.createLaplacePyr(img, 4);

% combine all levels in one image
[rows,cols,~] = cellfun(@size, p);
dst = zeros([rows(1)+rows(2) cols(1) size(img,3)], class(img));
r_start = [0 repmat(rows(1),1,numel(p)-1)] + 1;
r_end = [rows(1) rows(1)+rows(2:end)];
c_start = [0 0 cumsum(cols(2:end-1))] + 1;
c_end = [cols(1) cumsum(cols(2:end))];
for i=1:numel(p)-1
    dst(r_start(i):r_end(i),  c_start(i):c_end(i), :) = uint8(p{i});
end
figure, imshow(dst)
%imshow(cv.CLAHE(rgb2gray(dst)))

%% Image Blending using Pyramids
% One application of Pyramids is Image Blending. For example, in image
% stitching, you will need to stack two images together, but it may not look
% good due to discontinuities between images. In that case, image blending
% with Pyramids gives you seamless blending without leaving much data in the
% images. One classical example of this is the blending of two fruits, Orange
% and Apple.
%
% Please check first reference in additional resources, it has full
% diagramatic details on image blending, Laplacian Pyramids etc. Simply it is
% done as follows:
%
% * Load the two images of apple and orange
% * Find the Gaussian Pyramids for apple and orange (in this particular
%   example, number of levels is 6)
% * From Gaussian Pyramids, find their Laplacian Pyramids
% * Now join the left half of apple and right half of orange in each levels
%   of Laplacian Pyramids
% * Finally from this joint image pyramids, reconstruct the original image.
%
% Below is the full code. (For sake of simplicity, each step is done
% separately which may take more memory. You can optimize it if you want so).
%

% a pair of images of same size
A = imread(fullfile(mexopencv.root(),'test','apple.jpg'));
B = imread(fullfile(mexopencv.root(),'test','orange.jpg'));

% generate Laplacian Pyramids
pA = cv.Blender.createLaplacePyr(A, 5);
pB = cv.Blender.createLaplacePyr(B, 5);

% add left and right halves of images in each level
p = pB;
for i=1:numel(p)
    idx = round(size(pA{i},2)/2);
    p{i}(:,1:idx,:) = pA{i}(:,1:idx,:);
end
%{
% add top and bottom halves of images in each level
p = pB;
for i=1:numel(p)
    idx = round(size(pA{i},1)/2);
    p{i}(1:idx,:,:) = pA{i}(1:idx,:,:);
end
%}

% reconstruct
C = cv.Blender.restoreImageFromLaplacePyr(p);
C = uint8(C);

% image with direct connecting each half
idx = round(size(A,2)/2);
D = B;
D(:,1:idx,:) = A(:,1:idx,:);

% show images
figure
subplot(221), imshow(A), title('Apple')
subplot(222), imshow(B), title('Orange')
subplot(223), imshow(C), title('Pyramid Blending')
subplot(224), imshow(D), title('Direct Connection')

%% Additional Resources
%
% <http://pages.cs.wisc.edu/~csverma/CS766_09/ImageMosaic/imagemosaic.html Image Blending>
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% MORLET_2D_PYRAMID computes the 2-D elliptic Morlet filter given a set of
%    parameters in spatial domain
%
% Usage
%    gab = MORLET_2D_PYRAMID(N, M, sigma, slant, xi, theta, offset)
%
% Input
%    N (numeric): width of the filter
%    M (numeric): height of the filter
%    sigma (numeric): standard deviation of the envelope
%    slant (numeric): excentricity of the elliptic envelope
%            (the smaller slant, the larger angular resolution)
%    xi (numeric):  the frequency peak
%    theta (numeric): orientation in radians of the filter
%    offset (numeric): 2-by-1 Vvector index of the row and column of the
%       center of the filter (if offset is [0,0] the filter is centered in
%       [1,1])
%
% Output
%    gab (numeric): N-by-M matrix representing the gabor filter in spatial
%       domain.
%
% Description
%    Compute a Morlet wavelet in spatial domain.
%
%    Morlet wavelets have a 0 DC component.
%
% See also
%    GABOR_2D, MORLET_2D_NODC

function gab = morlet_2d_pyramid(N, M, sigma, slant, xi, theta, offset)

	if ~exist('offset', 'var')
		offset = [floor(N/2), floor(M/2)];
	end

	[x , y] = meshgrid(1:M, 1:N);

	x = x - offset(2) - 1;
	y = y - offset(1) - 1;

	Rth = rotation_matrix_2d(theta);
	A = Rth\  [1/sigma^2, 0 ; 0 slant^2/sigma^2] * Rth ;
	s = x.* ( A(1,1)*x + A(1,2)*y) + y.*(A(2,1)*x + A(2,2)*y ) ;

	%normalize sucht that the maximum of fourier modulus is 1
	gaussian_envelope = exp( - s/2);
	oscilating_part = gaussian_envelope .* exp(1i*(x*xi*cos(theta) + y*xi*sin(theta)));
	K = sum(oscilating_part(:)) ./ sum(gaussian_envelope(:));
	gabc = oscilating_part - K.*gaussian_envelope;

	gab=1/(2*pi*sigma^2/slant^2)*(gabc);

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function component = partition(model, solver, algorithm)
% swiftcc++ and fastcc++ augment swiftcc and fastcc by this preprocess
%
% USAGE:
%
%    component = partition(model, solver, algorithm)
%
% INPUTS:
%    model:        the metabolic network reconstruction
%                    * .S - the associated sparse stoichiometric matrix
%                    * .lb - feasible flux distribution lower bound
%                    * .ub - feasible flux distribution uppper bound
%                    * .rxns - cell array of reaction abbreviations
%                    * .rev - the 0-1 indicator vector of the reversible reactions
%    solver:       the LP solver to be used; the currently available options
%                  are 'gurobi', 'linprog', and 'cplex' with the default value
%                  of 'linprog'. It fallbacks to the COBRA LP solver interface
%                  if another supported solver is called.
%    algorithm:    the backend algorithm to be utilized between 'swift' and 'fast'
%
% OUTPUT:
%    component:    the index set of the reactions constituting the maximum
%                  flux consistent metabolic subnetwork
%
% NOTE:
%
%    requires bioinformatics toolbox
%
% .. Authors:
%       - Mojtaba Tefagh, Stephen P. Boyd, 2019, Stanford University

    assert(license('test', 'bioinformatics_toolbox') == 1, ...
        'The required Bioinformatics toolbox is not available!');
    S = model.S;
    rev = model.rev;
    lb = model.lb;
    ub = model.ub;
    c = model.c;
    rxns = model.rxns;
    [m, n] = size(S);

    %% constructing the directed graph
    DG = zeros(m+1);
    for i = 1:n
        head = S(:, i) > 0;
        if ~any(head)
            head = m+1;
        end
        tail = S(:, i) < 0;
        if ~any(tail)
            tail = m+1;
        end
        DG(head, tail) = 1;
        if rev(i)
            DG(tail, head) = 1;
        end
    end

    %% finding strongly or weakly connected components in the graph
    DG = sparse(DG);
    [~, C] = graphconncomp(DG, 'Directed', true);
    C = C(1:end-1);

    %% finding weakly connected components in the graph
    if range(C) == 0
        [~, C] = graphconncomp(max(DG, DG.'), 'Directed', false);
        C = C(1:end-1);
    end

    %% partitioning the metabolic network
    component = zeros(n, 1);
    for i = 1:n
        v = C(S(:, i) ~= 0);
        component(i) = all(v == v(1))*v(1);
    end
    if range(component) == 0
        component = zeros(n, 1);
        if strcmp(algorithm, 'fast')
            component(fastcc(model, getCobraSolverParams('LP', 'feasTol')*100)) = 1;
        elseif strcmp(algorithm, 'swift')
            component(swiftcc(S, rev, solver)) = 1;
        end
    else
        newcomponent = zeros(n, 1);
        for i = unique(component).'
            core = component == i;
            if sum(core) > 1
                model.S = S(:, core);
                model.rev = rev(core);
                model.lb = lb(core);
                model.ub = ub(core);
                model.c = c(core);
                model.rxns = rxns(core);
                newcomponent(core) = partition(model, solver, algorithm);
            end
        end
        component = newcomponent;
    end
end

author: opencobra

repo: cobratoolbox


text:
classdef Anodal2gausComputer < handle

    properties (Access = public)
        A_nodal_2_gauss
    end

    properties (Access = private)
        nnode
        nelem
        npnod
        ngaus
        connec
        shape
    end

    methods (Access = public)
        function obj = Anodal2gausComputer(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.computeA();
        end

        function intX = integrateP1FunctionWithShapeFunction(obj,cParams)
            ndof = size(obj.A_nodal_2_gauss{1},2);
            intX = zeros(ndof,1);
            for igaus = 1:obj.ngaus
                dVG  = cParams.dV(:,igaus);
                xG   = cParams.x(:,igaus);
                A    = obj.A_nodal_2_gauss{igaus};
                intX = intX + A'*(xG.*dVG);
            end
        end
    end

    methods (Access = private)
        function init(obj,cParams)
            obj.nnode  = cParams.nnode;
            obj.nelem  = cParams.nelem;
            obj.npnod  = cParams.npnod;
            obj.ngaus  = cParams.ngaus;
            obj.connec = cParams.connec;
            obj.shape  = cParams.shape;
        end

        function computeA(obj)
            A0    = sparse(obj.nelem,obj.npnod);
            A2g   = cell(obj.ngaus,1);
            nodes = obj.connec;
            for igaus = 1:obj.ngaus
                A2g{igaus} = A0;
                for inode = 1:obj.nnode
                    node   = nodes(:,inode);
                    shapeN = obj.shape(inode,igaus);
                    Ni = ones(obj.nelem,1)*shapeN;
                    A  = sparse(1:obj.nelem,node,Ni,obj.nelem,obj.npnod);
                    A2g{igaus} = A2g{igaus} + A;
                end
            end
            obj.A_nodal_2_gauss = A2g;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function updateBins = checkThresholds(data,params,binVals)
% function updateBins = checkThresholds(data,params,binVals)
% Checks if the binVals are still accurate.

updateBins = true;
if isempty(binVals)|| size(binVals,1)~=params.numBins,
  return;
end

nex = size(data,1);
numDim = size(data,2);
numBins = params.numBins;
temp = linspace(0,100,numBins+2);
prcValues = temp(2:end-1);
mpoint = round(numel(prcValues)/2);
% selpts = [1 mpoint numel(prcValues)];
selpts = [1 numel(prcValues)];
selValues = prcValues(selpts);
binVals = permute(binVals,[3 2 1]);
curbvals = binVals(:,:,selpts);

kk1 = bsxfun(@lt,data,curbvals);
kk2 = bsxfun(@le,data,curbvals);

tt1 = sum(kk1,1)/nex*100;
tt2 = sum(kk2,1)/nex*100;
tt1 = squeeze(tt1);
tt2 = squeeze(tt2);
mm = abs(tt1-tt2);
dd = abs(bsxfun(@minus,tt1,selValues));

dd(mm>0.2) = 0;

if max(dd(:))<(100/numBins),
  updateBins = false;
end

author: kristinbranson

repo: JAABA


text:
function mrk= mrk_evenlyInBlocks(blk, msec, varargin)

% MRK_EVENLYINBLOCKS - inserts additional markers between the existing
% markers, starting msec after the existing markers.
%
% Synopsis:
%   MRK= mrk_evenlyInBlocks(BLK, MSEC, <OPT>)
%
% Arguments:
%   BLK:  [STRUCT] structure defining blocks. It must have a field 'ival',
%         with each row defining a time interval ('block') in msec.
%
%   MSEC: [DOUBLE] distance of markers within each block in milliseconds
%
% OPT - struct or property/value list of optional fields/properties:
%      .OffsetStart -  specify offset in milliseconds after which the first
%                      after an existing marker block is to be set (default 0)
%      .OffsetEnd -    minimum length between block and next marker
%                      (default msec)
%
% Returns:
%   MRK: marker structure

% Author: Benjamin B
% 7-2010: Documented, extended, cleaned up (Matthias T)
% 5-2015 adapted to the new toolbox (Laura A, Benjamin B, Markus W)

props= {'OffsetEnd'     msec   '!DOUBLE[1]'
        'OffsetStart'   0      '!DOUBLE[1]'
       };

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props, 1);
opt_checkProplist(opt, props);

misc_checkType(blk, 'STRUCT(ival)');
misc_checkType(blk.ival, 'DOUBLE[- 2]', 'blk.ival');
misc_checkType(msec, '!DOUBLE[1]');

mrk= struct('time',[], 'event',struct);
mrk.event= struct('blkno',[]);

if isfield(blk, 'y'),
  [nClasses, nBlocks]= size(blk.y);
  mrk.y= zeros(nClasses,0);
  mrk.className= blk.className;
end

% fields from blk.event will be adapted and added to mrk.event. Let's prepare:
if isfield(blk, 'event'),
  eventFields= fieldnames(blk.event);
else
  eventFields= {};
end
for Fld= eventFields,
  fld= Fld{1};
  mrk.event.(fld)= [];
end

nBlocks= size(blk.ival,1);
for bb= 1:nBlocks,
  new_time= blk.ival(bb,1)+opt.OffsetStart:msec:blk.ival(bb,2)-opt.OffsetEnd;
  nMrk= length(new_time);
  mrk.time= cat(2, mrk.time, new_time);
  mrk.event.blkno= cat(1, mrk.event.blkno, bb*ones(nMrk,1)); % blkno has to be column vector for later functions
  % adapt fields from blk.event and add to mrk.event
  for Fld= eventFields,
    fld= Fld{1};
    val= blk.event.(fld)(bb);
    mrk.event.(fld)= cat(1, mrk.event.(fld), repmat(val, [nMrk 1]));
  end
  if isfield(blk, 'y'),
    new_y= zeros(nClasses, nMrk);
    iClass= find(blk.y(:,bb));
    new_y(iClass,:)= 1;
    mrk.y= cat(2, mrk.y, new_y);
  end
end

author: bbci


text:
function v = vec(a)

v = a(:);

author: andrewssobral

repo: lrslibrary


text:
function [yu,dyudtheta] = u_GaussianBumps(Theta,u,inF)

% input function for free-form deterministic DCM

spread = exp(Theta(1));
scale = exp(Theta(inF.indscale));
centres = exp(Theta(inF.indcentres));
n = length(centres);
yu = 0;
dyudtheta = zeros(size(Theta));
for i=1:n
   dt = centres(i)-u;
   bump = scale(i)*exp(-0.5*dt.^2./spread);
   dyudtheta(1) = ...
       dyudtheta(1) + scale(i)*dt^2./(2*spread*exp(dt^2/(2*spread)));
   dyudtheta(inF.indscale(i)) = bump;
   dyudtheta(inF.indcentres(i)) = ...
       -scale(i)*centres(i)*dt./(spread*exp(dt^2/(2*spread)));
   yu = yu + bump;
end

author: MBB-team

repo: VBA-toolbox


text:
function [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr)
% Multi-taper time-frequency coherence - two point processes given as times
%
% Usage:
%
% [C,phi,S12,S1,S2,t,f,zerosp,confC,phistd,Cerr]=cohgrampt(data1,data2,movingwin,params,fscorr)
% Input:
% Note units have to be consistent. Thus, if movingwin is in seconds, Fs
% has to be in Hz. see chronux.m for more information.
%
%       data1  (structure array of spike times with dimension trials; also accepts 1d array of spike times) -- required
%       data2  (structure array of spike times with dimension trials; also accepts 1d array of spike times) -- required
%       movingwin (in the form [window winstep] -- required
%       params: structure with fields tapers, pad, Fs, fpass, err, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%                    Note that T has to be equal to movingwin(1).
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%           trialave (average over trials when 1, don't average when 0) - optional. Default 0
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
%                 1 (use finite size corrections) - optional
%                (available only for spikes). Defaults 0.
% Output:
%       C (magnitude of coherency time x frequencies x trials for trialave=0;
%              time x frequency for trialave=1)
%       phi (phase of coherency time x frequencies x trials for no trial averaging;
%              time x frequency for trialave=1)
%       S12 (cross spectrum - time x frequencies x trials for no trial averaging;
%              time x frequency for trialave=1)
%       S1 (spectrum 1 - time x frequencies x trials for no trial averaging;
%              time x frequency for trialave=1)
%       S2 (spectrum 2 - time x frequencies x trials for no trial averaging;
%              time x frequency for trialave=1)
%       t (time)
%       f (frequencies)
%       zerosp (1 for windows and trials where spikes were absent (in either channel),zero otherwise)
%       confC (confidence level for C at 1-p %) - only for err(1)>=1
%       phistd - theoretical/jackknife (depending on err(1)=1/err(1)=2) standard deviation for phi
%                Note that phi + 2 phistd and phi - 2 phistd will give 95% confidence
%                bands for phi - only for err(1)>=1
%       Cerr  (Jackknife error bars for C - use only for Jackknife - err(1)=2)

if nargin < 3; error('Need data1 and data2 and window parameters'); end;
if nargin < 4; params=[]; end;

if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
end

[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
if nargin < 5 || isempty(fscorr); fscorr=0; end;

if nargout > 10 && err(1)~=2;
    error('Cerr computed only for Jackknife. Correct inputs and run again');
end;
if nargout > 8 && err(1)==0;
    error('Errors computed only if err(1) is not equal to zero');
end;

[N,Ch]=check_consistency(data1,data2);
[mintime1,maxtime1]=minmaxsptimes(data1);
[mintime2,maxtime2]=minmaxsptimes(data2);
mintime=min(mintime1,mintime2);
maxtime=max(maxtime1,maxtime2);

tn=mintime+movingwin(1)/2:movingwin(2):maxtime-movingwin(1)/2;
Nwin=round(Fs*movingwin(1)); % number of samples in window
% Nstep=round(movingwin(2)*Fs); % number of samples to step through
nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
f=getfgrid(Fs,nfft,fpass); Nf=length(f);
params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers
nw=length(tn);
if trialave;
   C=zeros(nw,Nf);
   S12=zeros(nw,Nf);
   S1=zeros(nw,Nf);
   S2=zeros(nw,Nf);
   phi=zeros(nw,Nf);
   Cerr=zeros(2,nw,Nf);
%    phierr=zeros(2,nw,Nf);
   phistd=zeros(nw,Nf);
else
   C=zeros(nw,Nf,Ch);
   S12=zeros(nw,Nf,Ch);
   S1=zeros(nw,Nf,Ch);
   S2=zeros(nw,Nf,Ch);
   phi=zeros(nw,Nf,Ch);
   Cerr=zeros(2,nw,Nf,Ch);
%    phierr=zeros(2,nw,Nf,Ch);
   phistd=zeros(nw,Nf,Ch);
end;
zerosp=zeros(nw,Ch);

for n=1:nw;
   t=linspace(tn(n)-movingwin(1)/2,tn(n)+movingwin(1)/2,Nwin);
   datawin1=extractdatapt(data1,[t(1) t(end)]);datawin2=extractdatapt(data2,[t(1) t(end)]);
   if nargout==11;
     [c,ph,s12,s1,s2,f,zsp,confc,phie,cerr]=coherencypt(datawin1,datawin2,params,fscorr,t);
%      phierr(1,n,:,:)=squeeze(phie(1,:,:));
%      phierr(2,n,:,:)=squeeze(phie(2,:,:));
     phistd(n,:,:)=phie;
     Cerr(1,n,:,:)=squeeze(cerr(1,:,:));
     Cerr(2,n,:,:)=squeeze(cerr(2,:,:));
   elseif nargout==10;
     [c,ph,s12,s1,s2,f,zsp,confc,phie]=coherencypt(datawin1,datawin2,params,fscorr,t);
%      phierr(1,n,:,:)=squeeze(phie(1,:,:));
%      phierr(2,n,:,:)=squeeze(phie(2,:,:));
     phistd(n,:,:)=phie;
   else
     [c,ph,s12,s1,s2,f,zsp]=coherencypt(datawin1,datawin2,params,fscorr,t);
   end;
   C(n,:,:)=c;
   phi(n,:,:)=ph;
   S12(n,:,:)=s12;
   S1(n,:,:)=s1;
   S2(n,:,:)=s2;
   zerosp(n,:)=zsp;
end;
t=tn;
C=squeeze(C); phi=squeeze(phi);S12=squeeze(S12); S1=squeeze(S1); S2=squeeze(S2);zerosp=squeeze(zerosp);
if nargout > 9; confC=confc; end;
if nargout==11;Cerr=squeeze(Cerr);end;
% if nargout==10; phierr=squeeze(phierr);end
if nargout==10; phistd=squeeze(phistd);end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ total_bboxes ] = apply_correction( total_bboxes, corrections, add1 )
%APPLY_CORRECTION Summary of this function goes here
%   Detailed explanation goes here

    % Perform correction based on regression values
    bbw = total_bboxes(:,3) - total_bboxes(:,1);
    bbh = total_bboxes(:,4) - total_bboxes(:,2);

    % TODO is this needed?
    if(add1)
        bbw = bbw + 1;
        bbh = bbh + 1;
    end

    new_min_x = total_bboxes(:,1) + corrections(:,1) .* bbw;
    new_min_y = total_bboxes(:,2) + corrections(:,2) .* bbh;
    new_max_x = total_bboxes(:,3) + corrections(:,3) .* bbw;
    new_max_y = total_bboxes(:,4) + corrections(:,4) .* bbh;
    score = total_bboxes(:,5);
    total_bboxes = [new_min_x, new_min_y, new_max_x, new_max_y, score];

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function savetSeries(tSeries,vw,scan,slice,nii)
% Save time series to a nifti file
%
%   savetSeries(tSeries,vw,scan,[slice],[nii])
%
% This function should be called everytime you save a tSeries. Makes the
% tSeries directory and scan subdirectory if they don't already exist.
%
% Nifti functionality has been added. In the directory
% '[sessionDir]/[dataTYPE]/TSeries/tSeriesScan[Scan#]Slice[Slice#].nii.gz'
%
% tSeries: matrix of N time samples by M pixels
% vw, scan, & slice: Used only to determine the full path for the tSeries file.
%
% djh, 2/17/2001
mrGlobals;

if notDefined('scan'),  scan  = viewGet(vw, 'curScan'); end
if notDefined('slice'), slice = viewGet(vw, 'curSlice'); end

viewType = viewGet(vw,'View Type');

if strcmp(viewType,'Inplane')

    tseriesdir = viewGet(vw,'tSeriesDir', 1);

    if ~exist(fullfile(tseriesdir),'dir')
        mkdir(tseriesdir);
    end

    pathStr = fullfile(tseriesdir,['tSeriesScan',num2str(scan),'.nii.gz']);

    %Let's also make sure that the dataTYPES are updated properly
    %First, check if this datatype exists
    %Then, check if this scan exists
    %Then, change the path for this scan to the new path

    %This is the same code as in duplicateDataType
    curDt = viewGet(vw,'Cur Dt');
    if curDt > numel(dataTYPES)
        newName = [];
        dlg = 'Please give a new name for the new datatype';
        while isempty(newName)
            newName = inputdlg(dlg,'Name new datatype');
            if iscell(newName) && ~isempty(newName)
                newName = deblank(newName{1}); % remove possible blanks
                if existDataType(newName)
                    dlg = 'What you just input is already in dataTYPES. Try again !!!';
                    newName = []; disp(['Warning: ',dlg]);
                end
            else
                dlg = 'You MUST give a NEW name for the new datatype !!!';
                disp(['Warning: ',dlg]);
            end %if
        end %while
        addDataType(newName);
    end %if

    numScans = dtGet(dataTYPES(curDt),'N Scans');

    if numScans < scan
        for newScan = (numScans+1):(scan-numScans)
             %Copy over the scan information from the previous scan
             vw = initScan(vw,curDt,newScan);
        end %for
    end %if

    %We now have a path to save the nifti to, let's make a nifti!
    if ~exist('nii','var')
        % We have no Nifti, we need to create our own
        % We can take the data from loadtSeries
        [~,nii] = loadtSeries(vw,scan,slice);

        % Now we need to change the data, the dimensions and the filepath
        % to reflect the new data that was passed in the tSeries

        % Let's first reshape the data into the proper size

        sizeDim = viewGet(vw,'Functional Slice Dim');
        totalSize = [size(tSeries,1),sizeDim,size(tSeries,3)];
        keepFrames = zeros(numScans, 2);
        keepFrames(1:numScans,2) = -1;
        newtSeries = reshape(tSeries,totalSize);
        newtSeries = permute(newtSeries, [2 3 4 1]); %Put the time at the end
        %Now let's update the nifti
        dim = size(newtSeries);
        nii = niftiSet(nii,'Dim',dim);
        nii = niftiSet(nii,'Data',newtSeries);

    else %if
        dim = niftiGet(nii, 'dim');
    end

    %Now we need to add another dimension to the pixdim if it only has 3
    if numel(niftiGet(nii,'Pix Dim')) < numel(dim)
        %Add another pixdim
        nii = niftiSet(nii, 'Pix Dim',[niftiGet(nii,'Pix Dim') 1]);
    end

	nii = niftiSet(nii,'File Path',pathStr);

    niftiWrite(nii,pathStr);

    dataTYPES(curDt) = dtSet(dataTYPES(curDt),'Inplane Path',pathStr,scan);
    dataTYPES(curDt) = dtSet(dataTYPES(curDt),'Size',dim(1:2),scan);

    %Add in keepFrames information
    dataTYPES(curDt) = dtSet(dataTYPES(curDt),'Keep Frames',keepFrames,scan);

    verbose = prefsVerboseCheck;
    if verbose > 1		% starting to use graded levels of feedback
        fprintf('Saved time series %s. (%s)\n', pathStr, datestr(now));
    end %if

    %Save the new dataTypes variable back to the overall session.
    saveSession;

else

    tseriesdir = tSeriesDir(vw, 1);
    scandir = ['Scan',num2str(scan)];
    if ~exist(fullfile(tseriesdir,scandir),'dir')
        mkdir(tseriesdir,scandir);
    end
    pathStr = fullfile(tseriesdir,scandir,['tSeries',num2str(slice)]);

    if strcmp(viewType,'Flat')
        % In the flat view, we will be getting multiple slices, so go
        % through each slice and save down the tSeries
        % We may sometimes encounter the case where we have 3 or 4
        % dimensions, and we want to save them all. We need to account for
        % that:

        numEle = numel(size(tSeries));
        numSlices = size(tSeries,numEle);

        if numEle == 4

            for i = 1:numSlices
                pathStr = fullfile(tseriesdir,scandir,['tSeries',num2str(i)]);
                tSeriesTmp = tSeries(:,:,i,:); %The slices will be on the 3rd dimension

                tSeriesTmp = single(tSeriesTmp);

                %disp(['Saving: ',pathStr]);
                save(pathStr,'tSeriesTmp');
            end %for

        elseif numEle ==3

            for i = 1:numSlices
                pathStr = fullfile(tseriesdir,scandir,['tSeries',num2str(i)]);
                tSeriesTmp = tSeries(:,:,i);

                tSeriesTmp = single(tSeriesTmp);

                %disp(['Saving: ',pathStr]);
                save(pathStr,'tSeriesTmp');
            end %for

        else
            %We don't know what to do
            error('Incorrect number of elements in tSeries.');
        end %if

    else
        %This is what happens in the Gray view, since only 1 slice
        tSeries = single(tSeries); %#ok<NASGU>

        %disp(['Saving: ',pathStr]);
        save(pathStr,'tSeries');
    end %if

    verbose = prefsVerboseCheck;
    if verbose > 1		% starting to use graded levels of feedback
        fprintf('Saved time series %s. (%s)\n', pathStr, datestr(now));
    end %if

end %if

return

author: vistalab

repo: vistasoft


text:
function [ num_int, pint ] = halfspace_normal_triangle_int_3d ( p, normal, t )

%*****************************************************************************80
%
%% HALFSPACE_NORMAL_TRIANGLE_INT_3D: intersection ( normal halfspace, triangle ) in 3D.
%
%  Discussion:
%
%    The normal form of a halfspace in 3D may be described as the set
%    of points (X,Y,Z) on or "above" a plane described in normal form:
%
%      P is a point on the plane,
%      NORMAL is the unit normal vector, pointing "out" of the
%      halfspace.
%
%    The triangle is specified by listing its three vertices.
%
%    The intersection may be described by the number of vertices of the
%    triangle that are included in the halfspace, and by the location of
%    points between vertices that separate a side of the triangle into
%    an included part and an unincluded part.
%
%    0 vertices, 0 separators    (no intersection)
%    1 vertex, 0 separators      (point intersection)
%    2 vertices, 0 separators    (line intersection)
%    3 vertices, 0 separators    (triangle intersection)
%
%    1 vertex, 2 separators,     (intersection is a triangle)
%    2 vertices, 2 separators,   (intersection is a quadrilateral).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(3,3), the vertices of the triangle.
%
%    Input, real P(3,1), a point on the bounding plane
%    that defines the halfspace.
%
%    Input, real NORMAL(3,1), the components of the normal vector
%    to the bounding plane that defines the halfspace.  By convention, the
%    normal vector points "outwards" from the halfspace.
%
%    Output, integer NUM_INT, the number of intersection points returned,
%    which will always be between 0 and 4.
%
%    Output, real PINT(3,4), the coordinates of the NUM_INT
%    intersection points.  The points will lie in sequence on the triangle.
%    Some points will be vertices, and some may be separators.
%
  dim_num = 3;
%
%  Compute the signed distances between the vertices and the plane.
%
  d = - normal(1:dim_num,1)' * p(1:dim_num,1);
%
%  Compute the signed distances between the vertices and the plane.
%
  dist1 = d + normal(1:dim_num,1)' * t(1:dim_num,1);
  dist2 = d + normal(1:dim_num,1)' * t(1:dim_num,2);
  dist3 = d + normal(1:dim_num,1)' * t(1:dim_num,3);
%
%  Now we can find the intersections.
%
  [ num_int, pint ] = halfspace_triangle_int_3d ( t, dist1, dist2, dist3 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Housekeeping
clc
clear
close all

%% Using the Grobner solver: constant-parameter case
m0=rise('rbcs0');
m1=solve(m0);
print_solution(m1)

% Finding all solutions
%-----------------------
m2=solve(m0,'solver','dsge_groebner','solve_check_stability',false);
print_solution(m2)
print_solution(m1)

%% Using the Grobner solver: regime switching example
clc

m00=rise('rbcs','solve_perturbation_type',{'frwz',{'mu'}});

% Finding one solution
%----------------------
m01=solve(m00);

% Finding all solutions
%-----------------------
m02=solve(m00,'solver','dsge_groebner','solve_check_stability',false);
print_solution(m01)
print_solution(m02)

%% Using the dsge_udc solver: regime switching example
clc

m000=set(m00,'solve_perturbation_type','mw');

% Finding one solution
%----------------------
m001=solve(m000);

% Finding all solutions Groebner
%-------------------------------
m002=solve(m000,'solver','dsge_groebner','solve_check_stability',false);

% Finding all solutions UDC
%--------------------------
refine=false;checkStab=false;allSols=true;msvOnly=true;
xplosRoots=false;debug=false;
solver={'dsge_udc',refine,checkStab,allSols,msvOnly,xplosRoots,debug};
m003=solve(m000,'solver',solver,'solve_check_stability',false);

% Finding all solutions UDC with refinement
%------------------------------------------
refine=true;checkStab=false;allSols=true;msvOnly=true;
xplosRoots=false;debug=false;
solver={'dsge_udc',refine,checkStab,allSols,msvOnly,xplosRoots,debug};

m004=solve(m000,'solver',solver,'solve_check_stability',false);

author: jmaih


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

%*****************************************************************************80
%
%% CN_LEG_03_XIU implements the Xiu precision 3 rule for region CN_LEG.
%
%  Discussion:
%
%    The rule has order
%
%      O = 2 * N.
%
%    The rule has precision P = 3.
%
%    CN_LEG is the cube [-1,+1]^N with the Legendre weight function
%
%      w(x) = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongbin Xiu,
%    Numerical integration formulas of degree two,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 1515-1520.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  o = 2 * n;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  expon = 0;
  volume = c1_leg_monomial_integral ( expon );
  volume = volume ^ n;

  for j = 1 : o

    i = 0;
    for r = 1 : floor ( n / 2 )
      arg = ( 2 * r - 1 ) * j * pi / n;
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * cos ( arg ) / sqrt ( 3.0 );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * sin ( arg ) / sqrt ( 3.0 );
    end

    if ( i < n )
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * r8_mop ( j ) / sqrt ( 3.0 );
      if ( n == 1 )
        x(i,j) = x(i,j) / sqrt ( 2.0 );
      end
    end

  end

  w(1:o) = volume / o;

  return
end