text:
function y = quadratic(x, sigma, type)
%QUADRATIC   Quadratic "robust" function.
%   QUADRATIC(X, SIGMA, TYPE) evaluates the Quadratic "robust" function
%   with sigma SIGMA at point(s) X.
%   TYPE selects the evaluation type:
%    - 0: function value
%    - 1: first derivative
%    - 2: second derivative
%
%   This is a private member function of the class 'robust_function'.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date:  $
%   $Revision: $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  switch (type)
   case 0
    y = x.^2 / sigma^2;
   case 1
    y = 2 * x / sigma^2;
   case 2
    y = repmat(2 / sigma^2, size(x));
  end

author: kristinbranson

repo: JAABA


text:
  function [y, eo_fun] = embed_out(y, Mdim)
%|function [y, eo_fun] = embed_out(y, Mdim)
%|
%| Dual or adjoint of embed_in().  See embed_in() for details.
%| To be used around _back() routines that expect [(M) *L]
%| input and produce [(N) *L] output that then must be reshaped
%| to be either [(N) (L)] or [np *L]
%|
%| in
%|	y	[*M (L)] or [(M) (L)]	input data array(s), possibly as columns
%|	Mdim				(M)
%| out
%|	y	[(M) *L]		output data arrays
%|	eo_fun	strum object with methods:
%|	x = eo_fun.shape(x, mask, np)	reshape x from [(N) *L] to be either
%|					[np (L)] or [(N) (L)], depending on y
%|
%| Copyright 2006-12-9, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(y, 'test'), embed_out_test, clear y, return, end
if nargin < 2, Mdim = size(y,1); end % trick: try to guess...

if size(y,1) == prod(Mdim) % convert [*M (L)] to [(M) *L]
	state.column = true;
	dimy = size(y);
	state.diml = dimy(2:end); % [(L)]
	y = reshape(y, [Mdim prod(state.diml)]); % [(M) *L]

else % convert [(M) (L)] to [(M) *L]
	state.column = false;
	dimi = [size(y) 1];
	jf_equal(dimi(1:length(Mdim)), Mdim)
	state.diml = dimi(length(Mdim)+1:end); % (L) (possibly empty)
	y = reshape(y, [Mdim prod(state.diml) 1]); % [(N) *L]
end

eo_fun = strum(state, {'shape', @embed_out_shape});

%
% embed_out_shape()
%
function x = embed_out_shape(state, x, mask, np)

if nargin < 4, np = sum(mask(:)); end

diml = state.diml;

if state.column % column in yields column out, i.e., [(N) *L] to [np (L)]
	if any(diml > 1)
		x = reshape(x, numel(mask), []); % [*N *L]
		x = x(mask,:); % [np *L]
		x = reshape(x, [np diml]); % [np (L)]
	else
		x = x(mask); % [np,1]
	end

else % [(N) *L] to [(N) (L)]
	if any(diml > 1)
		x = reshape(x, [size(mask) diml]); % [(N) (L)]
	end
end

%
% embed_out_test
%
function embed_out_test
ig = image_geom('nx', 10, 'ny', 8, 'dx', 1);
ig.mask = ig.circ > 0;

Mdim = [6 7];
y0 = ones(Mdim); % [(M)]
dl = [2 3];

[y2 eo] = embed_out(y0(:), Mdim); % single column
x2 = eo.shape(ig.ones, ig.mask);
jf_equal(size(x2), [ig.np 1])

[y2 eo] = embed_out(y0, Mdim); % usual 2d
x2 = eo.shape(ig.ones, ig.mask);
jf_equal(size(x2), [ig.nx ig.ny])

yr = repmat(y0, [1 1 dl]); % [(M) (L)]
[y2 eo] = embed_out(yr, Mdim);
x2 = eo.shape(ones([ig.nx ig.ny prod(dl)]), ig.mask);
jf_equal(size(x2), [ig.nx ig.ny dl])

yr = repmat(y0(:), [1 dl]); % [*M (L)]
[y2 eo] = embed_out(yr, Mdim);
x2 = eo.shape(ones([ig.nx ig.ny prod(dl)]), ig.mask);
jf_equal(size(x2), [ig.np dl])

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function write_lines_with_prefix(obj, fn, decorator)
% write lines in mfile with prefix
%
% write_lines_with_prefix(obj, fn, decorator)
%
% Inputs:
%   obj                 MOcovMFile instance
%   fn                  output filename
%   decorator           function handle which takes a single argument with
%                       a line number, and returns a prefix string that
%                       forms an expression that can be evaluated. When
%                       called, each executable line is prefixed with its
%                       corresponding prefix expression and the results are
%                       written to the output file fn.
%
% Notes:
%   - the typical value for decorator is a function that returns
%       mocov_line_covered(fn,k)
%     with k the input line number. This provides functionality to keep
%     track of which lines are covered.

    orig_lines=get_lines(obj);
    n=numel(orig_lines);

    executable=get_lines_executable(obj);

    new_lines=cell(1,n);
    for k=1:n
        line=orig_lines{k};

        if executable(k)
            prefix=decorator(k);
        else
            prefix='';
        end

        new_lines{k}=sprintf('%s%s\n',prefix,line);
    end

    pth=fileparts(fn);
    mkdir_recursively(pth);

    fid=fopen(fn,'w');
    cleaner=onCleanup(@()fclose(fid));
    fprintf(fid,'%s',new_lines{:});

function mkdir_recursively(pth)
    if ~isempty(pth) && ~mocov_util_isfolder(pth)
        parent=fileparts(pth);
        mkdir_recursively(parent);
        mkdir(pth);
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% LOWPASSFILTER - Constructs a low-pass butterworth filter.
%
% usage: f = lowpassfilter(sze, cutoff, n)
%
% where: sze    is a two element vector specifying the size of filter
%               to construct [rows cols].
%        cutoff is the cutoff frequency of the filter 0 - 0.5
%        n      is the order of the filter, the higher n is the sharper
%               the transition is. (n must be an integer >= 1).
%               Note that n is doubled so that it is always an even integer.
%
%                      1
%      f =    --------------------
%                              2n
%              1.0 + (w/cutoff)
%
% The frequency origin of the returned filter is at the corners.
%
% See also: HIGHPASSFILTER, HIGHBOOSTFILTER, BANDPASSFILTER
%

% Copyright (c) 1999 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% October 1999
% August  2005 - Fixed up frequency ranges for odd and even sized filters
%                (previous code was a bit approximate)

function f = lowpassfilter(sze, cutoff, n)

    if cutoff < 0 | cutoff > 0.5
	error('cutoff frequency must be between 0 and 0.5');
    end

    if rem(n,1) ~= 0 | n < 1
	error('n must be an integer >= 1');
    end

    if length(sze) == 1
	rows = sze; cols = sze;
    else
	rows = sze(1); cols = sze(2);
    end

    % Set up X and Y matrices with ranges normalised to +/- 0.5
    % The following code adjusts things appropriately for odd and even values
    % of rows and columns.
    if mod(cols,2)
	xrange = [-(cols-1)/2:(cols-1)/2]/(cols-1);
    else
	xrange = [-cols/2:(cols/2-1)]/cols;
    end

    if mod(rows,2)
	yrange = [-(rows-1)/2:(rows-1)/2]/(rows-1);
    else
	yrange = [-rows/2:(rows/2-1)]/rows;
    end

    [x,y] = meshgrid(xrange, yrange);
    radius = sqrt(x.^2 + y.^2);        % A matrix with every pixel = radius relative to centre.
    f = ifftshift( 1.0 ./ (1.0 + (radius ./ cutoff).^(2*n)) );   % The filter

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [AtA,A] = corrMatrix4D(obj,i)
% calucate 4D correlation matrix
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

nCha = size(obj.kCalib{i},4);

% A = [];

if(isreal(obj.kCalib))
    A = zeros(prod(obj.calibSize{i} - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision);
else
    A = complex(zeros(prod(obj.calibSize{i} - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision),zeros(prod(obj.calibSize{i} - obj.kernelSize + 1), prod(obj.kernelSize)*nCha,obj.measPara.precision));
end
if(isempty(A))
    error('corrMatrix4D(): Unable to create GRAPPA kernel. Check kernel and calibration size dimensionality');
end
counter = 1;
for n=1:nCha
    if(isreal(obj.kCalib{i}(:,:,:,n)))
        if(strcmp(obj.measPara.precision,'single'))
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colRSingle(obj.kCalib{i}(:,:,:,n),obj.kernelSize).';
        else
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colR(obj.kCalib{i}(:,:,:,n),obj.kernelSize).'; % before: tmp =
        end
    else
        if(strcmp(obj.measPara.precision,'single'))
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colCSingle(obj.kCalib{i}(:,:,:,n),obj.kernelSize).';
        else
            A(:,counter:counter+prod(obj.kernelSize)-1) = im3colC(obj.kCalib{i}(:,:,:,n),obj.kernelSize).';
        end
    end
    counter = counter + prod(obj.kernelSize);
end

AtA = A'*A;

end

author: thomaskuestner


text:
% Maximize stopband attenuation of a lowpass IIR filter
% "Linear Programming Design of IIR Digital Filters with Arbitrary
%  Magnitude Functions" by L.R. Rabiner, N.Y. Graham, and H.D. Helms
% (figures are generated)
%
% Designs a lowpass IIR filter using spectral factorization method where we:
% - minimize maximum stopband attenuation
% - have a constraint on the maximum passband ripple
%
%   minimize   max |H(w)|                      for w in the stopband
%       s.t.   1/delta <= |H(w)| <= delta      for w in the passband
%
% where we now have a rational frequency response function:
%
%   H(w) = sum_{m=0}^{M-1} b_m exp{-jwm} / sum_{n=1}^{N-1} a_n exp{-jwn}
%
% We change variables via spectral factorization method and get:
%
%   minimize   max R(w)                        for w in the stopband
%       s.t.   (1/delta)^2 <= R(w) <= delta^2  for w in the passband
%              R(w) >= 0                       for all w
%
% where R(w) is the squared magnited of the frequency response
% (and the Fourier transform of the autocorrelation coefficients r).
% We represent R(w) = N_hat(w)/D_hat(w), where now R(w) is a rational
% function since we deal with IIR filter (see the reference paper).
%
% Variables are coeffients of the numerator, denoted as c, and
% denominator, denoted as d. delta is the allowed passband ripple.
% This is a quasiconvex problem and can be solved using bisection.
%
% Written for CVX by Almir Mutapcic 02/02/06

%*********************************************************************
% user's filter specs (for a low-pass filter example)
%*********************************************************************
% number of coefficients for the IIR filter (including the zeroth one)
% (also without loss of generality we can assume that d_0 = 1, which
%  is the zeroth coefficient of the autocorrelation denominator)
M = 4; % nominator
N = 4; % denominator

wpass = 0.12*pi;   % end of the passband
wstop = 0.24*pi;   % start of the stopband
delta = 1;         % maximum passband ripple in dB (+/- around 0 dB)

%*********************************************************************
% create optimization parameters
%*********************************************************************
% rule-of-thumb discretization (from Cheney's Approx. Theory book)
sample_order = 30;
m = 15*(sample_order);
w = linspace(0,pi,m)'; % omega

% A's are matrices used to compute the power spectrum
Anum = [ones(m,1) 2*cos(kron(w,[1:M-1]))];
Aden = [ones(m,1) 2*cos(kron(w,[1:N-1]))];

% passband 0 <= w <= w_pass
ind = find((0 <= w) & (w <= wpass));    % passband
Ap_num  = Anum(ind,:);
Ap_den  = Aden(ind,:);

% transition band is not constrained (w_pass <= w <= w_stop)

% stopband (w_stop <= w)
ind = find((wstop <= w) & (w <= pi));   % stopband
As_num  = Anum(ind,:);
As_den  = Aden(ind,:);

%********************************************************************
% optimization
%********************************************************************

% use bisection (on the log of vars) to solve for the min stopband atten
Us_top = 1e-0; % 0 dB
Us_bot = 1e-6; % -60 dB (in original variables)

while( 20*log10(Us_top/Us_bot) > 1)
  % try to find a feasible design for given specs
  Us_cur = sqrt(Us_top*Us_bot);

  % formulate and solve the magnitude design problem
  cvx_begin quiet
    variable c(M,1)
    variable d(N-1,1)

    % feasibility problem
    % passband constraints
    (Ap_num*c) <= (10^(+delta/20))^2*(Ap_den*[1;d]); % upper constr
    (Ap_num*c) >= (10^(-delta/20))^2*(Ap_den*[1;d]); % lower constr
    % stopband constraint
    (As_num*c) <= (Us_cur)*(As_den*[1;d]); % upper constr
    % nonnegative-real constraint
    Anum*c >= 0;
    Aden*[1;d] >= 0;
  cvx_end

  % bisection
  if ~any(isnan(c)) % feasible
    fprintf(1,'Problem is feasible for stopband atten = %3.2f dB\n', ...
               10*log10(Us_cur));
    Us_top = Us_cur;
    b = spectral_fact(c);
    a = spectral_fact([1;d]);
  else % not feasible
    fprintf(1,'Problem not feasible for stopband atten = %3.2f dB\n', ...
               10*log10(Us_cur));
    Us_bot = Us_cur;
  end
end

% display the max attenuation in the stopband (convert to original vars)
fprintf(1,'\nOptimum min stopband atten is between %3.2f and %3.2f dB.\n',...
        10*log10(Us_bot),10*log10(Us_top));
disp('Optimal IIR filter coefficients are: ')
disp('Numerator: '), b
disp('Denominator: '), a

%*********************************************************************
% plotting routines
%*********************************************************************
% frequency response of the designed filter, where j = sqrt(-1)
H = ([exp(-j*kron(w,[0:M-1]))]*b)./([exp(-j*kron(w,[0:N-1]))]*a);

% magnitude plot
figure(1)
subplot(2,1,1)
plot(w,20*log10(abs(H)), ...
     [0 wpass],[delta delta],'r--', ...
     [0 wpass],[-delta -delta],'r--', ...
     [wstop pi],[10*log10(Us_top) 10*log10(Us_top)],'r--')
xlabel('w')
ylabel('mag H(w) in dB')
axis([0 pi -55 10]);

% phase plot
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w'), ylabel('phase H(w)')

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function varargout = im2metaimage( str, res, scale, crop, ext, file )
% IM2METAIMAGE  Read a batch of image files, collate them and save as a
% single MetaImage "file" (actually, one .mha and one .raw file)
%
% IM2METAIMAGE(STR)
%
%   STR is a string with the path and output mha header file name. For
%   example,
%
%     STR='/home/john/data/study01/study01.mha'; % linux
%     STR='C:\data\study01\study01.mha';         % windows
%
%   It is assumed that the MetaImage file (.mha) will be created
%   from the image files found in the target directory.
%
% IM2METAIMAGE(STR, RES, SCALE, CROP, EXT, FILE)
%
%   RES is a 3-vector with the pixel size in the x-, y- and z-coordinates.
%   By default RES=[1 1 1];
%
%   SCALE is a scalar factor to reduce the size of each frame (the number
%   of frames doesn't change). By default, SCALE=1.
%
%   CROP is a vector with values [xmin, xmax, ymin, ymax, zmin, zmax].
%   These values allow to crop the data volume. It follows the same
%   convention as Seg3D [1]. For example, if xmin=0, xmax=3, then voxels 1,
%   2 and 3 are selected.
%   Note that the output MetaImage file will have an offset header so that
%   the data retains its true coordinates.
%
%   EXT is a string with the file extension of the images. By default,
%   EXT='tif'.
%
%   FILE is the list of image files. By default, all the files in the
%   target directory are read, but FILE allows to change the system listing
%   order, select a subset of files, etc. This is useful, for example, when
%   the system listing order does not correspond to the order in which you
%   want to collate the slices. FILE is expected to have the same struct
%   format as the output of function DIR().
%
%   If you have filenames like
%
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.10.2009.10.23.16.02.38.562500.9699412.BMP
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.11.2009.10.23.16.02.38.562500.9699451.BMP
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.1.2009.10.23.16.02.38.562500.9699056.BMP
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.12.2009.10.23.16.02.38.562500.9699490.BMP
%
%   where the slice order is given by the field between the "6" and the
%   "2009", then you can use Gerardus function sortdirbynum() to order the
%   slices correctly. Example:
%
%   >> str1 = 'MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.*.BMP';
%   >> str2 = 'MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.mha';
%   >> file = sortdirbynum( str1, 5, '.' );
%   >> im2metaimage( str2, [], [], [], 'BMP', file );
%
% IM = IM2METAIMAGE(...);
%
%   IM is the image volume.
%
% [1] http://www.sci.utah.edu/cibc/software/seg3d.html

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2009 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error( nargchk( 1, 6, nargin, 'struct' ) );
error( nargoutchk( 0, 1, nargout, 'struct' ) );

% defaults
if ( nargin < 2 || isempty( res ) )
    res = [ 1, 1, 1 ];
end
if ( nargin < 3 || isempty( scale ) )
    scale = 1;
end
if ( nargin < 4 || isempty( crop ) )
    crop = [];
end
if ( nargin < 5 || isempty( ext ) )
    ext = 'tif';
end
if ( nargin < 6 || isempty( file ) )
    file = [];
end

% check arguments
if ( scale <= 0 || scale > 1 )
    error( 'SCALE must be in (0, 1]' );
end

% get path to the data files
[ dirdata, name ] = fileparts( str );

% get list of image files
if isempty( file )
    file = dir( [ dirdata filesep '*.' ext ] );
end

% adjust the image resolution to the scaling factor
res( 1:2 ) = res( 1:2 ) / scale;

% load first slice
frame = imread( [ dirdata filesep file( 1 ).name ] );

% remove alpha channel or other layers if present
frame = frame( :, :, 1 );

% init volume to load data
im = zeros( [ size( frame )*scale length( file ) ], class( frame ) );

% load slices of whole volume
im( :, :, 1 ) = imresize( frame, scale, 'bilinear' );
for I = 2:length( file )
    frame = imread( [ dirdata filesep file( I ).name ] );
    frame = frame( :, :, 1 );
    % resize frames
    im( :, :, I ) = imresize( frame, scale, 'bilinear' );
end

% crop volume (this is less efficient than initialising the volume as
% cropped, and then cropping each frame, but it makes the code simpler)
%
% The "+1" is necessary to follow the Seg3D convention
if ~isempty( crop )
    im = im( ...
        crop( 1 )+1:crop( 2 ), ...
        crop( 3 )+1:crop( 4 ), ...
        crop( 5 )+1:crop( 6 ) ...
        );
end

% compute offset in the coordinates of the cropped volume
if isempty( crop )
    offset = [ 0.0 0.0 0.0 ];
else
    offset = crop( [ 1 3 5 ] ) .* res;
end

% write mha file
writemetaimagefile([dirdata filesep name '.mha'], im, res, offset);

% avoid outputing the image volume unless the user has requested it
% explicitly
if ( nargout == 0 )
    varargout = [];
else
    varargout = { im };
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% MAHAL2CONF - Translates a Mahalanobis distance into a confidence
%              interval.  Consider a multivariate Gaussian
%              distribution of the form
%
%   p(x) = 1/sqrt((2 * pi)^d * det(C)) * exp((-1/2) * MD(x, m, inv(C)))
%
%              where MD(x, m, P) is the Mahalanobis distance from x
%              to m under P:
%
%                 MD(x, m, P) = (x - m) * P * (x - m)'
%
%              A particular Mahalanobis distance k identifies an
%              ellipsoid centered at the mean of the distribution.
%              The confidence interval associated with this ellipsoid
%              is the probability mass enclosed by it.
%
%              If X is an d dimensional Gaussian-distributed vector,
%              then the Mahalanobis distance of X is distributed
%              according to the Chi-squared distribution with d
%              degrees of freedom.  Thus, the confidence interval is
%              determined by integrating the chi squared distribution
%              up to the Mahalanobis distance of the measurement.
%
% Usage:
%
%   c = mahal2conf(m, d);
%
% Inputs:
%
%   m    - the Mahalanobis radius of the ellipsoid
%   d    - the number of dimensions of the Gaussian distribution
%
% Outputs:
%
%   c    - the confidence interval, i.e., the fraction of
%          probability mass enclosed by the ellipsoid with the
%          supplied Mahalanobis distance
%
% See also: CONF2MAHAL

% Copyright (C) 2002 Mark A. Paskin
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
% USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = mahal2conf(m, d)

c = chi2cdf(m, d);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = linearCompressibility(C,varargin)
% computes the linear compressibility of an elasticity tensor
%
% Description
%
% $$\beta(x) = S_{ijkk} x_i x_j$$
%
% Input
%  C - elastic @stiffnessTensor
%  x - list of @vector3d
%
% Output
%  beta - linear compressibility in directions v
%

% compute linear compressibility from complience
[varargout{1:nargout}] = linearCompressibility(inv(C),varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function oki = straightCIv1(actionstr)

%	STRAIGHT command interpreter

%       STRAIGHT --- Speech Transformation and Representation using
%       Adaptive Interpolation of weiGHTed spectrum
%	Interactive interface for trial use
%
%       (c) Copyright ATR Human Information Processing Research
%       Laboratories, 1996,1997
%       Author and Inventor: Hideki Kawahara
%       02/July/1996
%       04/July/1996
%       07/July/1996
%	07/Sep./1996
%	01/Nov./1996
%	25/Nov./1996 TEMPO enhancement
%	30/Nov./1996 This version is with AG window trick
%	25/Dec./1996 This version is with AG window trick
%	27/Dec./1996 using fundamentalness as U/V measure
%	02/Feb./1997 without V/UV decision
%	05/Feb./1997 Bug Fix
%	10/Feb./1997 Maximum liklihood F0 tracking
%	15/Feb./1997 frequency range division
%	15/Feb./1997 Big bug fix in spectral calculation
%	22/Feb./1997 No need for temporal smoothing (thanks for the new window)
%	11/Mar./1997 New wavelet for F0 extraction
%	03/May /1997 Quick bug fix for Matlab v5
%	30/May /1997 Quick bug fix for Matlab v5
%	21/June/1997 Quick bug fix for new fundamentalness
%	12/Aug./1997 Revised optimum smoother (minor revision)
%	13/Aug./1997 Revised SPIKES scaling property and TD enhancement
%	19/Dec./1997 Revised for mixed mode excitation
%	09/Jan./1998 GUI command line interpreter
%	29/Jan./1998 Integrated with mixed mode and GUI
%	02/Feb./1998 Minor inconvenience fix
%	03/Fev./1998 Minor bug fix
%	13/April/1999 New pitch extractor
%	19/April/1999 Minor modificatio to display information
%	22/April/1999 Minor modificatio for unvoiced spectrum
%	05/May  /1999 New periodic/aperiodic mixing ratio
%	21/July/1999  reduced version with crisp V/UV decision
%	17/Feb./2001 Bug fix for multi channel sound files
%	30/May/2001 New aperiodicity measure and control
%	8/April/2002 modified to remove magic numbers
%	26/June/2002 minor modification: two globals are added.
%	08/Dec./2002 Default parameters were made user definable.
%   04/Feb./2003 Aperiodicity correction and bug fix in STRAIGHTbody proc.
%   31/Aug./2004 Bug fix for f0 abnormal values
%	30/April/2005 modification for Matlab v7.0 compatibility

global n2sgram nsgram n3sgram n2sgrambk nwsgram xold x f0floor f0ceil fs framem shiftm f0shiftm ...
  fftl eta pc framel fftl2 acth pwth pcnv fconv sconv delsp gdbw cornf fname ofname delfracind ...
  tpath cpath ... %%paraminitialized
  mag delfrac hr f0raw f0l f0var f0varL sy pcorr pecorr ... %%gobjlist ...
  upsampleon hhb pwt pwh amp defaultendian indefaultendian outdefaultendian f0varbak

global maphandles

global bv % 02/Sept./1999
global apv dpv % 21/Sept./1999
global defaultch % 16/Feb./2001
global apve dpve % 22/Oct./2001

global f0v vrv % 25/June/2002 for debug
global ecrt % 03/Feb./2003 for C/N based correction of dpv

oki = true;

switch(actionstr)
  %-----------------------------------------------
  %	Initialization part
  %-----------------------------------------------
  case 'GUIinitialize'
    clear global
    straightCIv1 initializeparams;
    straightPanel98bak;
    syncgui;
    straightCIv1 syncbuttons;
  case 'resetparamsbtn'
    straightCIv1 initializeparams;
    syncgui;
  case 'initialize'
    clear global n2sgram nsgram n3sgram n2sgrambk nwsgram xold x f0raw f0l sy hhb
    straightCIv1 initializeparams;
    syncgui;
    straightCIv1 syncbuttons;
  case 'initializeparams'
    if exist('defaultparams.m','file')==0 %	08/Dec./2002
      f0floor=40;
      f0ceil=800;
      fs=22050;	% sampling frequency (Hz)
      framem=40;	% default frame length for pitch extraction (ms)
      shiftm=1;       % default frame shift (ms) for spectrogram
      f0shiftm=1;     % default frame shift (ms) for F0 information
      fftl=1024;	% default FFT length
      eta=1.4;        % time window stretch factor
      pc=0.6;         % exponent for nonlinearity
      mag=0.2;      % This parameter should be revised.
      framel=framem*fs/1000;

      if fftl < framel
        fftl=2^ceil(log(framel)/log(2));
      end;
      fftl2=fftl/2;
      defaultch=1; % 17/Feb./2001

      %-------------- Decision parameter for source information

      acth=0.5;	% Threshold for normalized correlation (dimension less)
      pwth=32;	% Threshold for instantaneous power below maximum (dB)

      %-----------------------------------------------------
      %       Synthesis parameters
      %-----------------------------------------------------

      pcnv=1.0; 	% pitch stretch
      fconv=1.0; 	% frequency stretch
      sconv=1.0; 	% time stretch

      %	  delsp=2; 	%  standard deviation of random group delay in ms
      delsp=0.5; 	%  standard deviation of random group delay in ms 26/June/2002
      gdbw=70; 	% smoothing window length of random group delay (in Hz)
      %	  cornf=3000;  	% corner frequency for random phase (Hz)
      cornf=4000;  	% corner frequency for random phase (Hz) 26/June 2002
      delfrac=0.2;  % This parameter should be revised.
      delfracind=0;

      %-----------------------------------------------------
      %	file parameters
      %-----------------------------------------------------

      fname='none';		% input data file name

      hr='on';
      tpath=pwd;
      if strcmp(computer,'MAC2')==0
        tpath=[tpath '/'];
      end;
      upsampleon=0;
    else %	08/Dec./2002
      defaultparams;
    end; % of if exist('defaultparams.m','file')==0 %	08/Dec./2002
    defaultendian=chkdefaultendian;
    indefaultendian=defaultendian;
    outdefaultendian=defaultendian;
  case 'resetvalues'
    straightCIv1 initializeparams
    syncgui;

    %-----------------------------------------------------
    %	file I/O part
    %-----------------------------------------------------
  case 'bininputformat'
    hh=findobj('Tag','bininputformat');
    indefaultendian=get(hh,'Value');
  case 'binoutputformat'
    hh=findobj('Tag','binoutputformat');
    outdefaultendian=get(hh,'Value');
  case 'readfile'
    [fname,cpath]=uigetfile(...
      {'*.wav';'*.aif';'*.WAV';'*.aiff';'*.dat';'*.dat'},...
      'sound file input');
    if fname(1)~=0
      tcpath=[char(39) cpath char(39)];
      eval(['cd ' tcpath]);
      if ~isempty(strfind(lower(fname),'.wav')) % 16/Feb./2001
        [x,fs]=audioread(fname);
        x=x*32768;
      elseif ~isempty(strfind(lower(fname),'.aif')) % 16/Feb./2001
        [x,fs]=aiffread(fname);
      else
        if indefaultendian==1
          fid=fopen(fname,'r','ieee-le');
        else
          fid=fopen(fname,'r','ieee-be');
        end;
        x=fread(fid,'short')';
        fclose(fid);
      end;
      [tnn,tmm]=size(x); % 16/Feb./2001
      if min(tnn,tmm)>1 % 16/Feb./2001
        switch tnn>tmm
          case 1,
            x=x(:,defaultch);
          case 0,
            x=x(defaultch,:);
        end;
      end;
      x=x(:)'; % make sure that the vector is row vector
      xold=x;
      x=xold+std(x)/1000*randn(size(x));  % 03/Feb./2001
    else
      disp('file input is cancelled. ');
      disp(' ');
    end;
    syncgui;
    straightCIv1 syncbuttons;

  case 'savefile'
    tcpath=[char(39) tpath char(39)];
    eval(['cd ' tcpath]);
    tsy=sy; tfs=fs;
    if upsampleon
      switch fs
        case {8000, 10000, 11025, 12000}
          tfs=fs*4; tsy=interp(sy,4);
        case {16000, 20000, 22050, 24000}
          tfs=fs*2; tsy=interp(sy,2);
      end;
    end;
    [ofname,tpath]=uiputfile('*','sound file output');
    if ofname(1)~=0
      if ~isempty(strfind(ofname,'.wav'))
        audiowrite([tpath ofname], tsy/32768,tfs);
      elseif ~isempty(strfind(ofname,'.aif'))
        ok=aiffwrite(tsy,tfs,16,ofname);
        if isempty(ok)
          disp(['File output is failed. ' ofname ' was not written.']);
        end;
      else
        if outdefaultendian==1
          fid2=fopen([tpath ofname],'w','ieee-le');
        else
          fid2=fopen([tpath ofname],'w','ieee-be');
        end;
        fwrite(fid2,tsy,'short');
        fclose(fid2);
      end;
      disp(['% Saved successfully as: ' ofname]);
      disp(' ');
    else
      disp('file output is cancelled. ');
      disp(' ');
    end;

    %-----------------------------------------------------
    %	Display spectrogram group
    %-----------------------------------------------------
  case 'dispnsgram'
    mxil=max(max(20*log10(nsgram+0.001)));
    [nsy,nsx]=size(nsgram);
    figure;
    imagesc((0:nsx-1)*shiftm,(0:nsy-1)/nsy*fs/2, ...
      max(20*log10(nsgram+0.001),mxil-50));
    axis('xy'); colormap(1-gray);
    title(['Equal resolution spectrum of ' fname ' ' date ' ' mktstr]);
    xlabel('time (ms)');
    ylabel('frequency (Hz)');
  case 'dispnwsgram'
    mxil=max(max(20*log10(nwsgram+0.001)));
    [nsy,nsx]=size(nwsgram);
    figure;
    imagesc((0:nsx-1)*shiftm,(0:nsy-1)/nsy*fs/2, ...
      max(20*log10(nwsgram+0.001),mxil-50));
    axis('xy'); colormap(1-gray);
    title(['Wide band spectrum of ' fname ' ' date ' ' mktstr]);
    xlabel('time (ms)');
    ylabel('frequency (Hz)');
  case 'dispn2sgram'
    mxil=max(max(20*log10(n2sgram+0.001)));
    [nsy,nsx]=size(n2sgram);
    figure;
    imagesc((0:nsx-1)*shiftm,(0:nsy-1)/nsy*fs/2,...
      max(real(20*log10(n2sgram+0.001)),mxil-50));
    axis('xy'); colormap(1-gray);
    title(['enhanced spectrum of ' fname ' ' date ' ' mktstr]);
    xlabel('time (ms)');
    ylabel('frequency (Hz)');
  case 'dispn2sgrambk'
    mxil=max(max(20*log10(n2sgrambk+0.001)));
    [nsy,nsx]=size(n2sgrambk);
    figure;
    imagesc((0:nsx-1)*shiftm,(0:nsy-1)/nsy*fs/2,...
      max(real(20*log10(n2sgrambk+0.001)),mxil-50));
    axis('xy'); colormap(1-gray);
    title(['Interpolated spectrum of ' fname ' ' date ' ' mktstr]);
    xlabel('time (ms)');
    ylabel('frequency (Hz)');
  case 'dispn3sgram'
    mxil=max(max(20*log10(n3sgram+0.001)));
    [nsy,nsx]=size(n3sgram);
    figure;
    imagesc((0:nsx-1)*shiftm,(0:nsy-1)/nsy*fs/2,...
      max(real(20*log10(n3sgram+0.001)),mxil-50));
    axis('xy'); colormap(1-gray);
    title(['enhanced spectrum (without 2nd strctr) of ' fname ' ' date ' ' mktstr]);
    xlabel('time (ms)');
    ylabel('frequency (Hz)');
  case 'disphhbspectrograme'
    bb=1:length(n2sgram(1,:));
    mmx=hhb(:,bb).*n3sgram+(1-hhb(:,bb)).*nwsgram;
    mxil=max(max(20*log10(mmx+0.001)));
    [nsy,nsx]=size(mmx);
    figure;
    imagesc((0:nsx-1)*shiftm,(0:nsy-1)/nsy*fs/2,...
      max(real(20*log10(mmx+0.001)),mxil-50));
    axis('xy'); colormap(1-gray);
    title(['final composite spectrum of ' fname ' ' date ' ' mktstr]);
    xlabel('time (ms)');
    ylabel('frequency (Hz)');
  case 'showF0'
    figure(gcf+1); subplot(111);plot((1:length(f0l))*f0shiftm,f0l); grid on;
    title(['F0 of ' fname ' ' date ' ' mktstr]);
    ylabel('frequency (Hz)'); xlabel('time (ms)');

    %-----------------------------------------------
    %	audio display part
    %-----------------------------------------------
  case 'playoriginal'
    straightsound(xold,fs);
  case 'playsynth'
    straightsound(sy,fs);

    %-----------------------------------------------
    %	parameter modification part
    %-----------------------------------------------
  case 'peekvars'  % This is the most poweful interaction
    keyboard;
    syncgui;
    straightCIv1 syncbuttons;
  case 'getfsmenu'
    fs=getfsfrommenu(gco);
    syncgui;
  case 'editf0ceil'
    f0ceil=getvalufromedit(gco,800);
    syncgui;
  case 'editf0floor'
    f0floor=getvalufromedit(gco,40);
    syncgui;
  case 'editshiftm'
    f0shiftm=getvalufromedit(gco,1);
    shiftm=f0shiftm;
    syncgui;
  case 'fftledit'
    syncgui;
  case 'wndwstrtchedit'
    eta=getvalufromedit(gco,1.4);
    syncgui;
  case 'pwrcnstntedit'
    pc=getvalufromedit(gco,0.6);
    syncgui;
  case 'magfactoredit'
    mag=getvalufromedit(gco,0.2);
    syncgui;
  case 'delfracedit'
    delfrac=getvalufromedit(gco,0.2);
    syncgui;
  case 'delspedit'
    delsp=getvalufromedit(gco,2);
    syncgui;
  case 'cornfedit'
    cornf=getvalufromedit(gco,2400);
    syncgui;
  case 'gdbwedit'
    gdbw=getvalufromedit(gco,70);
    syncgui;
  case 'pcnvedit'
    pcnv=getvalufromedit(gco,1);
    syncgui;
  case 'fconvedit'
    fconv=getvalufromedit(gco,1);
    syncgui;
  case 'sconvedit'
    sconv=getvalufromedit(gco,1);
    syncgui;
  case 'tpathedit'
    tpath=get(gco,'Value');
    syncgui;
  case 'pcnvslider'
    pcnv=10.0.^get(gco,'Value');
    syncgui;
  case 'fconvslider'
    fconv=3.0.^get(gco,'Value');
    syncgui;
  case 'sconvslider'
    sconv=10.0.^get(gco,'Value');
    syncgui;
  case 'delfracradio'
    delfracind=~delfracind;
    syncgui;
  case 'delspradio'
    delfracind=~delfracind;
    syncgui;
  case 'upsamplebtn'
    upsampleon=~upsampleon;
    syncgui;

    %--------------------------------------------------
    %	non-linear manipulations
    %--------------------------------------------------
  case 'FqNLbtn'
    hh=findobj('Tag','FqNLbtn');
    if get(hh,'UserData') ==0 ||  isempty(get(hh,'UserData'))
      set(hh,'UserData',1);
      set(hh,'BackgroundColor',[0.9 0.33333 0.33333]);
      bendline initialize;
    else
      set(hh,'UserData',0);
      set(hh,'BackgroundColor',[0.733333 0.733333 0.733333]);
      bendline close;
    end;

    %	This part is obsolate. This part will be revised completely.
  case 'interactSGRAM'
    figure;
    mxil=max(max(20*log10(n2sgram)));
    imagesc(max(20*log10(n2sgram),mxil-45));
    axis('xy'); colormap(1-gray);
    title(['Interpolated spectrum of ' fname ' ' date ' ' mktstr]);
    disp('%	Now you have to define trajectory using mouse');
    disp('%	Please type "return", if you are ready.');
    disp('%	It is recommended for you to select important portion ');
    disp('%	using "zoom on" command.');
    disp('%	Please do not forget to issue "zoom off" before continue.');
    disp('%	In graphical input interaction, click defines point and return');
    disp('%	notifies it is the last point.');
    keyboard;
    disp('%	Interaction started. Put the cursor inside the graphics.');
    zoom off;
    getTrace;
    disp('%	You can modify spectrum using the following command.');
    disp('%       n2sgram=n2sgrambk.*(1+nsgm).^X;');
    disp('%		X*6 dB amplification is made.');
    disp('%		Default 6dB amplification was already done.');
    disp('%	If you are OK, type "return". ');
    disp('%	Otherwise, please change.');
    n2sgram=n2sgrambk.*(1+nsgm);

    %-------------------------------------------
    %	mapping control part
    %-------------------------------------------
  case 'frequencymapmod'
    [nii,~]=size(n2sgram);
    vx=(0:nii-1)/(nii-1);
    idcv=vx*maphandles(20)+sin(vx*pi)*maphandles(21)+sin(2*pi*vx)*maphandles(22);
    fconv=max(1,min(nii,idcv*(nii-1)+1));

    %-------------------------------------------
    %	synthesis part
    %-------------------------------------------
  case 'synthesizechar' % This part is useless now
    disp('%-------- Current Synthesis parameters ------');
    disp(['%	delsp=' num2str(delsp) ...
      ';     % standard deviation of random group delay in ms']);
    disp(['%	gdbw=' num2str(gdbw) ...
      ';    % smoothing window length of random group delay (in Hz)']);
    disp(['%	cornf=' num2str(cornf) ...
      ';  % corner frequency for random phase (in Hz)']);
    disp(['%	pcnv=' num2str(pcnv) ...
      ';    % pitch stretch']);
    disp(['%	fconv=' num2str(fconv) ...
      ';   % frequency stretch']);
    disp(['%	sconv=' num2str(sconv) ...
      ';   % time stretch']);
    disp('%	');
    disp('%	If you are happy with these parameters please type "return".');
    disp('%	You can change these setting using Matlab command(s)');
    disp('%	If you want to restore default parameters please type');
    disp('%	"default22kparams"  There are similar prog. for 12k,16k files.');
    keyboard;
    disp('%	Now, synthesis is in progress. Please wait a moment.');
    syncgui;
    straightCIv1 synthesize

  case 'synthesizegradedqqq'  % OBSOLATE!!!
    hh=findobj('Tag','FqNLbtn');
    if ~isempty(hh)
      if ~isempty(get(hh,'UserData'))
        if get(hh,'UserData')==1
          straightCIv1 frequencymapmod
        end;
      end;
    end;
    sy=straightSynthTC01(n3sgram,nwsgram,f0raw,hhb,shiftm,fs, ...
      pcnv,fconv,sconv,gdbw,delfrac,delsp,cornf,delfracind);
    dBsy=powerchk(sy,fs,15); % 23/Sept./1999
    cf=(20*log10(32768)-22)-dBsy;
    sy=sy*(10.0.^(cf/20));
    disp('%	Done!');
    straightCIv1 syncbuttons;

  case 'synthesizegraded'
    hh=findobj('Tag','FqNLbtn');
    if ~isempty(hh)
      if ~isempty(get(hh,'UserData') )
        if get(hh,'UserData')==1
          straightCIv1 frequencymapmod
        end;
      end;
    end;
    sy=straightSynthTB07ca(n3sgram,f0raw,shiftm,fs, ...
      pcnv,fconv,sconv,gdbw,delfrac,delsp,cornf,delfracind, ...
      aperiodiccomp(apv,dpv,5,f0raw,f0shiftm),1); % 8/April/2002
    dBsy=powerchk(sy,fs,15); % 23/Sept./1999
    cf=(20*log10(32768)-22)-dBsy;
    sy=sy*(10.0.^(cf/20));
    disp('%	Done!');
    straightCIv1 syncbuttons;

  case 'synthesize'
    hh=findobj('Tag','FqNLbtn');
    if ~isempty(hh)
      if ~isempty(get(hh,'UserData') )
        if get(hh,'UserData')==1
          straightCIv1 frequencymapmod
        end;
      end;
    end;
    sy=straightSynthTB06(n3sgram,f0raw,f0var,f0varL,shiftm,fs, ...
      pcnv,fconv,sconv,gdbw,delfrac,delsp,cornf,delfracind);
    dBsy=powerchk(sy,fs,15); % 23/Sept./1999
    cf=(20*log10(32768)-22)-dBsy;
    sy=sy*(10.0.^(cf/20));
    disp('%	Done!');
    straightCIv1 syncbuttons;

    %------------------------------------------------------
    %	analysis part
    %	This part is modified to introduce a new F0 and
    %	source information extraction method. (19/April/1999)
    %------------------------------------------------------
  case 'source'
    nvo=24;
    nvc=ceil(log(f0ceil/f0floor)/log(2)*nvo);
    [f0v,vrv,dfv,~,aav]=fixpF0VexMltpBG4(xold,fs,f0floor,nvc,nvo,1.2,1,shiftm,1,5,0.5,1);
    title([fname '  ' datestr(now,0)]);
    %drawnow;
    [~,~]=size(f0v);
    subplot(614);
    [pwt,pwh]=plotcpower(xold,fs,shiftm);drawnow;

    [f0raw,irms,~,amp]=f0track5(f0v,vrv,dfv,pwt,pwh,aav,shiftm);
    f0t=f0raw;avf0=mean(f0raw(f0raw>0));
    f0t(f0t==0)=f0t(f0t==0)*NaN;tt=1:length(f0t);
    %	  keyboard;
    subplot(615);plot(tt*shiftm,f0t,'g');grid on;
    if ~isnan(avf0)
      axis([1 max(tt)*shiftm ...
        min(avf0/sqrt(2),0.95*min(f0raw(f0raw>0)))  ...
        max(avf0*sqrt(2),1.05*max(f0raw(f0raw>0)))]);
    end;
    ylabel('F0 (Hz)');
    %----------- 31/July/1999
    hold on;
    dn=floor(fs/(f0ceil*3*2)); % fix by H.K. at 28/Jan./2003
    [f0raw,ecr]=refineF06(decimate(xold,dn),fs/dn,f0raw,1024,1.1,3,f0shiftm,1,length(f0raw)); % 31/Aug./2004
    f0t=f0raw;%%avf0=mean(f0raw(f0raw>0));
    f0t(f0t==0)=f0t(f0t==0)*NaN;tt=1:length(f0t);
    subplot(615);plot(tt*shiftm,f0t,'k');hold off;
    drawnow
    %----------- 31/July/1999

    tirms=irms;
    tirms(f0raw==0)=tirms(f0raw==0)*NaN;
    tirms(f0raw>0)=-20*log10(tirms(f0raw>0));
    ecrt=ecr;
    ecrt(f0raw==0)=ecrt(f0raw==0)*NaN;
    subplot(616);hrms=plot(tt*shiftm,tirms,'g',tt*shiftm,20*log10(ecrt),'r'); %31/July/1999
    set(hrms,'LineWidth',2);hold on
    plot(tt*shiftm,-10*log10(vrv),'k.');
    grid on;hold off
    axis([1 max(tt)*shiftm -10 60]);
    xlabel('time (ms)');ylabel('C/N (dB)');
    drawnow;
    irmsz=irms*0;

    %---------- This part is for maintaining compatibility with old synthesis routine ----
    f0var=max(0.00001,irms-irmsz).^2;
    f0var(f0var>0.99)=f0var(f0var>0.99)*0+100;
    f0var(f0raw==0)=f0var(f0raw==0)*0+100;
    f0varbak = f0var;  % backup for f0var (18/July/1999)
    f0var=f0var/2;  %  2 is a magic number. If everything is OK, it should be 1.
    f0var=(f0var>0.9);  % This modification is to make V/UV decision crisp  (18/July/1999)
    f0varL=f0var;
    %-------------------------------------------------------------------------------------
    f0raw(f0raw<=0)=f0raw(f0raw<=0)*0; % safeguard 31/August/2004
    f0raw(f0raw>f0ceil)=f0raw(f0raw>f0ceil)*0+f0ceil; % safeguard 31/August/2004

    straightCIv1 syncbuttons;

    %--------------------------------------------------------------
    %	classic STRAIGHT with a single V/UV measure
    %--------------------------------------------------------------
  case 'straightcore'
    disp('% Now, adaptive window analysis has started. Please wait a moment.');
    [n2sgrambk,nsgram]=straightBodyC03ma(xold,fs,shiftm,fftl,f0raw,f0var,f0varL,eta,pc); %%
    if mag>0
      n2sgram=specreshape(fs,n2sgrambk,eta,pc,mag,f0raw);
    else
      n2sgram=n2sgrambk;
    end;
    straightCIv1 syncbuttons;

    %--------------------------------------------------------------
    %	revised STRAIGHT with a multi band graded V/UV decision (OBSOLATE!!)
    %--------------------------------------------------------------
  case 'bandcorrbtnqqq'
    [n2sgrambk,nsgram,nwsgram]= ...
      straightBodyB04m(xold,fs,shiftm,fftl,f0raw,eta,pc);
    straightCIv1 syncbuttons;
    if mag>0
      n2sgram=specreshape(fs,n2sgrambk,eta,pc,mag,f0raw);
    else
      n2sgram=n2sgrambk;
    end;

    [pcorr,pecorr]=BcorrMap(xold,fs,f0raw,shiftm);

    wvm3=wfromMap4(pcorr,pecorr,n2sgram,fs);
    emap=max(pcorr,pecorr);
    hh=wvm3'*emap;

    a=0.32;b=15;c=0.15;   % blending parameter; this is very tentative
    hhb=max(0,(1.0./(1+exp(-(hh-a)*b))-1.0/(1+exp(-(c-a)*b))) ...
      /(1.0/(1+exp(-(1-a)*b))-1.0/(1+exp(-(c-a)*b))));
    straightCIv1 syncbuttons;

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

    %   MBE type analysis  2/Sept./1999

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

  case 'bandcorrbtn'
        [n2sgrambk,nsgram]=straightBodyC03ma(xold,fs,shiftm,fftl,f0raw,f0var,f0varL,eta,pc);
    if mag>0
      n2sgram=specreshape(fs,n2sgrambk,eta,pc,mag,f0raw);
    else
      n2sgram=n2sgrambk;
    end;

    [apvq,dpvq,apve,dpve]=aperiodicpartERB2(xold,fs,f0raw,f0shiftm,5,fftl/2+1); % 10/April/2002
    apv=10*log10(apvq); % for compatibility
    dpv=10*log10(dpvq); % for compatibility
    %- ---------
    %   Notes on aperiodicity estimation: The previous implementation of
    %   aperiodicity estimation was sensitive to low frequency noise. It is a
    %   bad news, because environmental noise usually has its power in the low
    %   frequency region. The following corrction uses the C/N information
    %   which is the byproduct of fixed point based F0 estimation.
    %   by H.K. 04/Feb./2003
    %- ---------
    dpv=correctdpv(apv,dpv,5,f0raw,ecrt,f0shiftm,fs); % Aperiodicity correction 04/Feb./2003 by H.K.

    bv=boundmes2(apv,dpv,fs,f0shiftm,5,fftl/2+1);

    figure;

    semilogy((0:length(bv)-1)*f0shiftm,0.5./10.0.^(bv));grid on;
    straightCIv1 syncbuttons;

  case 'remove2ndstructue'
    n3sgram=rmv2nd(n2sgram,f0raw,fs);
    straightCIv1 syncbuttons;

  case 'bypassbtn'
    n3sgram=n2sgram;
    straightCIv1 syncbuttons;

    %-----------------------------------------------------------------
    %	suppress buttons which are nor appropriate
    %-----------------------------------------------------------------
  case 'syncbuttons'
    hh=findobj('Tag','analyzesrcbtn');
    if length(xold) >fftl
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','bandcorrbtn');
    if length(f0raw) >5
      set(hh,'Enable','on');  % Enabled again 02/Sept./1999

    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','analyzespcbtn');
    if length(f0raw) >5
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','bypassbtn');
    if length(n2sgram) >2
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','remove2ndbtn');
    if length(n2sgram) >2
      %		  set(hh,'Enable','on');  % commented out at 21/July/1999
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','synthgradbtn');
    if (length(bv) >2) && (length(n3sgram)>2)
      set(hh,'Enable','on');  % enabled again at 02/Sept./1999
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','synthesizebtn');
    if length(n3sgram) >2
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','savetobtn');
    if length(sy) >fftl
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','adaptivespecbtn');
    if length(nsgram) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','widespecbtn');
    if length(nwsgram) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','smoothedspecbtn');
    if length(n2sgrambk) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','dispn2sgrambtn');
    if length(n2sgram) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','removedspecbtn');
    if length(n3sgram) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','cmpstspecgrambtn');
    if length(hhb) >1
      %		  set(hh,'Enable','on');  %  commented out on 21/July/1999
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','playorgbtn');
    if length(xold) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;
    %----------------------------------------------------
    hh=findobj('Tag','playsynthbtn');
    if length(sy) >1
      set(hh,'Enable','on');
    else
      set(hh,'Enable','off');
    end;

end;
end

function defaultendian=chkdefaultendian
%	defaultendian	: 1-littel endian, 2-big endian

gg=computer;
switch gg(1:3)
  case 'PCW'
    defaultendian=1;
  case 'MAC'
    defaultendian=2;
  case 'SUN'
    defaultendian=2;
  case 'SOL'
    defaultendian=2;
  case 'HP7'
    defaultendian=1;
  case 'SGI'
    defaultendian=2;
  case 'ALP'
    defaultendian=1;
  case 'AXP'
    defaultendian=1;
  case 'LNX'
    defaultendian=1;
  otherwise
    defaultendian=2;
end;
end

author: HidekiKawahara


text:
classdef  (SharedTestFixtures={matlab.unittest.fixtures.PathFixture(...
        fullfile(fileparts(mfilename('fullpath')), '..', '..', '..', 'processingTools'), ...
        'IncludeSubfolders', true)})  ...
        testop_CSIApodize < matlab.unittest.TestCase
    properties
        MRSIStruct
        kMaxX
        kMaxY
        kx
        ky

    end
    methods(TestMethodSetup)
        function createPhantom(testCase)
            MRSI = struct();
            fovX = 200;
            fovY = 200;
            fovZ = 10;
            voxelSizeX = 10;
            voxelSizeY = 10;
            voxelSizeZ = 10;

            MRSI.fov.x = fovX;
            MRSI.fov.y = fovY;
            MRSI.fov.z = fovZ;
            MRSI.voxelSize.x = voxelSizeX;
            MRSI.voxelSize.y = voxelSizeY;
            MRSI.voxelSize.z = voxelSizeZ;
            x = -fovX/2 + voxelSizeX/2:voxelSizeX:fovX/2 - voxelSizeX/2;
            y = -fovY/2 + voxelSizeY/2:voxelSizeY:fovY/2 - voxelSizeY/2;
            z = -fovZ/2 + voxelSizeZ/2:voxelSizeZ:fovZ/2 - voxelSizeZ/2;
            MRSI.coordinates.x = x;
            MRSI.coordinates.y = y;
            MRSI.coordinates.z = z;

            kxFov = 1/(x(2) - x(1));
            kyFov = 1/(y(2) - y(1));
            deltaKx = 1/fovX;
            deltaKy = 1/fovY;
            testCase.kMaxX = kxFov/2;
            testCase.kMaxY = kyFov/2;
            testCase.kx = -kxFov/2 + deltaKx/2:deltaKx:kxFov/2 - deltaKx/2;
            testCase.ky = -kyFov/2 + deltaKy/2:deltaKy:kyFov/2 - deltaKy/2;
            MRSI.data = ones(1024, 20, 20);
            MRSI.dims.t = 1;
            MRSI.flags.apodized = 0;
            testCase.MRSIStruct = MRSI;

        end
    end

    methods(Test)

        function testApodizeCosine(testCase)
            MRSILocal = testCase.MRSIStruct;
            MRSILocal.flags.spatialFT = 0;
            MRSILocal.dims.kx = 2;
            MRSILocal.dims.ky = 3;
            apodized = op_CSIApodize(MRSILocal, "functionType", 'cosine');

            W = getWeights('cosine', testCase.kx, testCase.ky, testCase.kMaxX, testCase.kMaxY);
            realData = applyApodization(MRSILocal, W, 0);
            testCase.verifyEqual(apodized.data, realData);
        end

        function testApodizeHamming(testCase)
            MRSILocal = testCase.MRSIStruct;
            MRSILocal.flags.spatialFT = 0;
            MRSILocal.dims.kx = 2;
            MRSILocal.dims.ky = 3;
            apodized = op_CSIApodize(MRSILocal, "functionType", 'hamming');

            W = getWeights('hamming', testCase.kx, testCase.ky, testCase.kMaxX, testCase.kMaxY);
            realData = applyApodization(MRSILocal, W, 0);
            testCase.verifyEqual(apodized.data, realData);

        end

        function testApodizeGaussian(testCase)
            MRSILocal = testCase.MRSIStruct;
            MRSILocal.flags.spatialFT = 0;
            MRSILocal.dims.kx = 2;
            MRSILocal.dims.ky = 3;
            apodized = op_CSIApodize(MRSILocal, "functionType", 'gaussian');

            W = getWeights('gaussian', testCase.kx, testCase.ky, testCase.kMaxX, testCase.kMaxY);
            realData = applyApodization(MRSILocal, W, 0);

            testCase.verifyEqual(apodized.data, realData);
        end

        function testApodizeCosineSpatial(testCase)
            MRSILocal = testCase.MRSIStruct;
            MRSILocal.flags.spatialFT = 1;
            MRSILocal.dims.x = 2;
            MRSILocal.dims.y = 3;
            apodized = op_CSIApodize(MRSILocal, "functionType", 'cosine');

            W = getWeights('cosine', testCase.kx, testCase.ky, testCase.kMaxX, testCase.kMaxY);
            W = fourierTransformWeights(W);
            realData = applyApodization(MRSILocal, W, 1);

            testCase.verifyEqual(apodized.data, realData, 'RelTol', 1e-10);
        end

        function testApodizeHammingSpatial(testCase)
            MRSILocal = testCase.MRSIStruct;
            MRSILocal.flags.spatialFT = 1;
            MRSILocal.dims.x = 2;
            MRSILocal.dims.y = 3;
            apodized = op_CSIApodize(MRSILocal, "functionType", 'hamming');

            W = getWeights('hamming', testCase.kx, testCase.ky, testCase.kMaxX, testCase.kMaxY);
            W = fourierTransformWeights(W);
            realData = applyApodization(MRSILocal, W, 1);

            testCase.verifyEqual(apodized.data, realData, 'RelTol', 1e-10);
        end

        function testApodizeGaussianSpatial(testCase)
            MRSILocal = testCase.MRSIStruct;
            MRSILocal.flags.spatialFT = 1;
            MRSILocal.dims.x = 2;
            MRSILocal.dims.y = 3;
            apodized = op_CSIApodize(MRSILocal, "functionType", 'gaussian');

            W = getWeights('gaussian', testCase.kx, testCase.ky, testCase.kMaxX, testCase.kMaxY);
            W = fourierTransformWeights(W);
            realData = applyApodization(MRSILocal, W, 1);

            testCase.verifyEqual(apodized.data, realData, 'RelTol', 1e-10);
        end

        function testImageAndKSpace(testCase)
            MRSILocalK = testCase.MRSIStruct;
            MRSILocalK.flags.spatialFT = 0;
            MRSILocalK.dims.kx = 2;
            MRSILocalK.dims.ky = 3;
            apodizedK = op_CSIApodize(MRSILocalK, 'functionType', 'hamming');
            apodizedK.data = fourierTransformData(apodizedK.data, apodizedK.dims.kx, apodizedK.dims.ky);

            MRSILocalImage = testCase.MRSIStruct;
            MRSILocalImage.flags.spatialFT = 1;

            MRSILocalImage.dims.x = 2;
            MRSILocalImage.dims.y = 3;
            MRSILocalImage.data = fourierTransformData(MRSILocalImage.data, MRSILocalImage.dims.x, MRSILocalImage.dims.y);

            apodizedImage = op_CSIApodize(MRSILocalImage, 'functionType', 'hamming');

            testCase.verifyEqual(apodizedImage.data, apodizedK.data, 'AbsTol', 1e-10);
        end
    end
end

function realData = applyApodization(MRSILocal, W, isSpatial)
    realData = MRSILocal.data;
    realData = permute(realData, [2, 3, 1]);
    if(isSpatial)
        for iIndex = 1:size(realData, 3)
            realData(:, :, iIndex) = conv2(realData(:, :, iIndex), W, 'same');
        end
    else
        realData = realData .* W;
    end
    realData = permute(realData, [3, 1, 2]);
end

function wFT = fourierTransformWeights(W)
    if(mod(size(W, 1), 2) == 1)
        W = circshift(W, 1, 1);
    end
    if(mod(size(W, 2), 2) == 1)
        W = circshift(W, 1, 2);
    end

    wFT = fftshift(fftshift(fft2(fftshift(fftshift(W, 1), 2)), 1), 2);
    wFT = wFT/numel(wFT);
end

function W = getWeights(apodizationType, kX, kY, kMaxX, kMaxY)
    switch(apodizationType)
        case("cosine")
            W1 = cos(pi*(kX)/(2*kMaxX));
            W2 = cos(pi*(kY)/(2*kMaxY));
        case("gaussian")
            W1 = exp(-4*(kX/kMaxX).^2);
            W2 = exp(-4*(kY/kMaxY).^2);
        case('hamming')
            W1 = 0.54 + 0.46*cos(pi*kX/kMaxY);
            W2 = 0.54 + 0.46*cos(pi*kX/kMaxY);
        otherwise
            error('No function found!');
    end
    W = W1' * W2;
end

function data = fourierTransformData(data, kxDim, kyDim)
    if(mod(size(data, kxDim), 2) == 1)
        data = circshift(data, 1, kxDim);
    end
    if(mod(size(data, kyDim), 2) == 1)
        data = circshift(data, 1, kyDim);
    end

    data = fftshift(fft(fftshift(data, kxDim), [], kxDim), kxDim);
    data = fftshift(fft(fftshift(data, kyDim), [], kyDim), kyDim);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [varargout] = openDbConnection(view);
% [View,OK] = openDbConnection(view);
% OK = openDbConnection(view);
% openDbConnection(view);
%
% First overload is used for login, others - for  simple connection.
%
% In login mode: Check for session  existence. If  not exist - ask operator for his
% username, password etc. Then connect. If name doesn't exist - ask for  it  again and create
% a  new record  in the database.
% In simple mode: Just try to connect. If it fails - bring an error message.

defautValues = {'snarp.stanford.edu','vista','',''};

if(nargout<2) % SIMPLE MODE
  if(isfield(view,'mysqlsession')) % We are logged in already
    mysql('open',view.mysqlsession.hostname,view.mysqlsession.username,...
      view.mysqlsession.password);
    mysql('use','mrDataDB');
    varargout{1} = 1; % OK button
    return;
  else % Error Message
    myWarnDlg('You''re not logged in yet. Please, go to File/Login to database and log in.');
    varargout{1} = 0; % Cancel button
    return;
  end
else % ADVANCED MODE
  varargout{2} = 1; % OK button
  if(isfield(view, 'mysqlsession')) %  Connection was  already open
    Answer = questdlg('You have selected your database,username and password already. Do you wish to continue?',...
      'Alert','Cancel','OK','Cancel');
    if(strcomp(Answer,'Cancel'))
       varargout{1} = view;
       varargout{2} = 0;
       return;
    end
  end

  prompt={'Hostname:','Login:','Password:','Username (how are you listed in a database'};
  answer=inputdlg(prompt,'Connecting to database...',1,defautValues);

  if ~isempty(answer)
    mysqlsession = struct('hostname',answer{1},'username',answer{2},...
      'password',answer{3},'userId',0);
    view = setfield(view,'mysqlsession',mysqlsession);
    Flag = mysql('open',mysqlsession.hostname,mysqlsession.username,...
      mysqlsession.password);
    if(~exist('Flag')) % An error occured when opening a connection
      varargout{1} = view;
      varargout{2} = 0;
      return;
    end
    mysql('use','mrDataDB');
    view.mysqlsession.userId = checkForUser(answer{4}); % Transform nick name to an Id (see below)
    view.mysqlsession.userId
    varargout{1} = view;
    return;
  else % Cancel button
    varargout{1} = view;
    varargout{2} = 0;
    return;
  end
end

%%-------------------------------
function userId = checkForUser(userNick);
userId = mysql(['SELECT id FROM users WHERE username="',userNick,'"']);
if(max(size(userId))>1) % More than one user found
  userId = userId(1);
  myWarnDlg('More than one user with this username found. Using  the first one.');
  return;
elseif(min(size(userId))==0) % No such user in a database
  Answer = questdlg('There''re no users with this username in the database. Do you want to create a new record?',...
      'New record','Cancel','OK','Cancel');
  if(strcomp(Answer,'Cancel'))
    userId = 0;
    myWarnDlg('You''re listed as a GUEST user');
  else
    ResultFlag = 0;
    while(ResultFlag==0)
      answer=inputdlg({'First name:','Last name:','Organization:','E-mail:','Username:'}...
        ,'Create a new user',1,{'','','','',userNick});
      if(isempty(answer)) % Cancel button
        userId = 0;
        myWarnDlg('You''re listed as a GUEST user');
        return;
      end
      listedEmail = mysql(['SELECT id FROM users WHERE email="',answer{4},'"']);
      listedNick = mysql(['SELECT id FROM users WHERE username="',answer{5},'"']);
      if(~isempty(listedEmail))
        myWarnDlg('This e-mail is already in the database. Try again.');
      elseif(~isempty(listedNick))
        myWarnDlg('This username is already in the database. Try again.');
      else
        ResultFlag=1;
      end;
    end
    mysql(['INSERT INTO users (firstName,lastName,organization,email,username,notes) VALUES("',...
      answer{1},'","',answer{2},'","',answer{3},'","',answer{4},'","',answer{5},...
      '","mrLoadRet creation")']);
    disp(['New user ',answer{5},' might be created.']);
  end
  return;
else % Normal situation - one user found
  return;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = runtests_xunit(varargin)
%runtests Run unit tests
%   runtests runs all the test cases that can be found in the current directory
%   and summarizes the results in the Command Window.
%
%   Test cases can be found in the following places in the current directory:
%
%       * An M-file function whose name starts or ends with "test" or
%         "Test" and that returns no output arguments.
%
%       * An M-file function whose name starts or ends with "test" or
%         "Test" and that contains subfunction tests and uses the
%         initTestSuite script to return a TestSuite object.
%
%       * An M-file defining a subclass of TestCase.
%
%   runtests(dirname) runs all the test cases found in the specified directory.
%
%   runtests(packagename) runs all the test cases found in the specified
%   package. (This option requires R2009a or later).
%
%   runtests(mfilename) runs test cases found in the specified function or class
%   name. The function or class needs to be in the current directory or on the
%   MATLAB path.
%
%   runtests('mfilename:testname') runs the specific test case named 'testname'
%   found in the function or class 'name'.
%
%   Multiple directories or file names can be specified by passing multiple
%   names to runtests, as in runtests(name1, name2, ...) or
%   runtests({name1, name2, ...}, ...)
%
%   runtests(..., '-verbose') displays the name and result, result, and time
%   taken for each test case to the Command Window.
%
%   runtests(..., '-logfile', filename) directs the output of runtests to
%   the specified log file instead of to the Command Window.
%
%   out = runtests(...) returns a logical value that is true if all the
%   tests passed.
%
%   Examples
%   --------
%   Find and run all the test cases in the current directory.
%
%       runtests
%
%   Find and run all the test cases in the current directory. Display more
%   detailed information to the Command Window as the test cases are run.
%
%       runtests -verbose
%
%   Save verbose runtests output to a log file.
%
%       runtests -verbose -logfile my_test_log.txt
%
%   Find and run all the test cases contained in the M-file myfunc.
%
%       runtests myfunc
%
%   Find and run all the test cases contained in the TestCase subclass
%   MyTestCase.
%
%       runtests MyTestCase
%
%   Run the test case named 'testFeature' contained in the M-file myfunc.
%
%       runtests myfunc:testFeature
%
%   Run all the tests in a specific directory.
%
%       runtests c:\Work\MyProject\tests
%
%   Run all the tests in two directories.
%
%       runtests c:\Work\MyProject\tests c:\Work\Book\tests

%   Steven L. Eddins
%   Copyright 2009-2010 The MathWorks, Inc.

verbose = false;
logfile = '';
if nargin < 1
    suite = TestSuite.fromPwd();
else
    [name_list, verbose, logfile] = getInputNames(varargin{:});
    if numel(name_list) == 0
        suite = TestSuite.fromPwd();
    elseif numel(name_list) == 1
        suite = TestSuite.fromName(name_list{1});
    else
        suite = TestSuite();
        for k = 1:numel(name_list)
            suite.add(TestSuite.fromName(name_list{k}));
        end
    end
end

if isempty(suite.TestComponents)
    error('xunit:runtests:noTestCasesFound', 'No test cases found.');
end

if isempty(logfile)
    logfile_handle = 1; % File handle corresponding to Command Window
else
    logfile_handle = fopen(logfile, 'w');
    if logfile_handle < 0
        error('xunit:runtests:FileOpenFailed', ...
            'Could not open "%s" for writing.', logfile);
    else
        cleanup = onCleanup(@() fclose(logfile_handle));
    end
end

fprintf(logfile_handle, 'Test suite: %s\n', suite.Name);
if ~strcmp(suite.Name, suite.Location)
    fprintf(logfile_handle, 'Test suite location: %s\n', suite.Location);
end
fprintf(logfile_handle, '%s\n\n', datestr(now));

if verbose
    monitor = VerboseTestRunDisplay(logfile_handle);
else
    monitor = TestRunDisplay(logfile_handle);
end
did_pass = suite.run(monitor);

if nargout > 0
    out = did_pass;
end

function [name_list, verbose, logfile] = getInputNames(varargin)
name_list = {};
verbose = false;
logfile = '';
k = 1;
while k <= numel(varargin)
    arg = varargin{k};
    if iscell(arg)
        name_list = [name_list; arg];
    elseif ~isempty(arg) && (arg(1) == '-')
        if strcmp(arg, '-verbose')
            verbose = true;
        elseif strcmp(arg, '-logfile')
            if k == numel(varargin)
                error('xunit:runtests:MissingLogfile', ...
                    'The option -logfile must be followed by a filename.');
            else
                logfile = varargin{k+1};
                k = k + 1;
            end
        else
            warning('runtests:unrecognizedOption', 'Unrecognized option: %s', arg);
        end
    else
        name_list{end+1} = arg;
    end
    k = k + 1;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [hfig,finalim,final_pforefly] = PlotSampleFramesColor(trx,readframe,bkgdim0,predictions,mainfly,otherflies,ts,varargin)

colorpos = [.7,0,0];
colorneg = [0,0,.7];
border = 20;

[colorpos,colorneg,border,hfig,figpos,...
  cm,fg_thresh,bg_thresh,sigma_bkgd,wmah,frac_a_back,dist_epsilon,...
  ncolors_reference,solidbkgdcolor,...
  flipim,intensityoff] = ...
  myparse(varargin,'colorpos',colorpos,...
  'colorneg',colorneg,...
  'border',border,...
  'hfig',[],...
  'figpos',[],...
  'colormap',{@jet,@jet},...
  'fg_thresh',115/255,...
  'bg_thresh',10/255,...
  'sigma_bkgd',115/255-10/255,...
  'wmah',.5,...
  'frac_a_back',1,...
  'dist_epsilon',.1,...
  'ncolors_reference',[],...
  'solidbkgdcolor',[.85,.85,.85],...
  'flipim',false,...
  'intensityoff',0);

if isempty(otherflies),
  flygroups = {mainfly};
else
  flygroups = {mainfly,otherflies};
end
nflygroups = numel(flygroups);
nframes = numel(ts);
allts = ts(1):ts(end);
nframes_total = numel(allts);
allflies = [flygroups{:}];

[nr0,nc0,ncolors] = size(bkgdim0);

%% grab out relevant time points
% not rotating currently, maybe rotate in the future?

x = nan(trx.nflies,nframes_total);
y = nan(trx.nflies,nframes_total);
a = nan(trx.nflies,nframes_total);
b = nan(trx.nflies,nframes_total);
theta = nan(trx.nflies,nframes_total);
imcenter = [(nc0+1)/2,(nr0+1)/2]; %#ok<NASGU>
for fly = 1:trx.nflies,
  idx = allts+trx(fly).off;
  i0 = find(idx>=1,1);
  i1 = find(idx<=trx(fly).nframes,1,'last');
  if isempty(i0) || isempty(i1),
    continue;
  end
  x(fly,i0:i1) = trx(fly).x(idx(i0:i1));
  y(fly,i0:i1) = trx(fly).y(idx(i0:i1));
  a(fly,i0:i1) = trx(fly).a(idx(i0:i1));
  b(fly,i0:i1) = trx(fly).b(idx(i0:i1));
  theta(fly,i0:i1) = trx(fly).theta(idx(i0:i1));
end

bkgdim = bkgdim0;

%% figure out the axes limits
xlim = [inf,-inf];
ylim = [inf,-inf];
for fly = allflies,
  for j = 1:nframes_total,
    [x1,x2,y1,y2] = ellipse_to_bounding_box(x(fly,j),y(fly,j),a(fly,j),b(fly,j),theta(fly,j));
    xlim(1) = min([xlim(1),x1,x2]);
    xlim(2) = max([xlim(2),x1,x2]);
    ylim(1) = min([ylim(1),y1,y2]);
    ylim(2) = max([ylim(2),y1,y2]);
  end
end
xlim = [floor(xlim(1))-border,ceil(xlim(2))+border];
ylim = [floor(ylim(1))-border,ceil(ylim(2))+border];
bkgdim = bkgdim(ylim(1):ylim(2),xlim(1):xlim(2),:);
x = x - xlim(1)+1;
y = y - ylim(1)+1;

%% colors for each frame
colors_all = nan(ts(end)-ts(1)+1,3,nflygroups);
for i = 1:nflygroups,
  cmcurr = cm{min(numel(cm),i)};
  if isempty(cmcurr),
    colors_all(:,:,i) = 1;
  elseif isnumeric(cmcurr),
    colors_all(:,:,i) = cmcurr;
  elseif ~isempty(ncolors_reference),
    colors_reference = cmcurr(ncolors_reference);
    idxcurr = round(linspace(1,ncolors_reference,ts(end)-ts(1)+1));
    colors_all(:,:,i) = colors_reference(idxcurr,:);
  else
    colors_all(:,:,i) = cmcurr(ts(end)-ts(1)+1);
  end
end
colors = colors_all(ts-ts(1)+1,:,:);

%% read in the frames, convert to grayscale
for i = 1:nframes,
  t = ts(i);
  imcurr = readframe(t);
  imcurr = im2double(imcurr(ylim(1):ylim(2),xlim(1):xlim(2),:));
  if ncolors > 1,
    grayim_curr = rgb2gray(imcurr);
  else
    grayim_curr = imcurr;
    imcurr = repmat(imcurr,[1,1,3]);
  end
  if i == 1,
    ims = repmat(imcurr,[1,1,1,nframes]);
    gray_ims = repmat(grayim_curr,[1,1,nframes]);
  else
    ims(:,:,:,i) = imcurr;
    gray_ims(:,:,i) = grayim_curr;
  end
end

%% probability that each pixel is foreground
dbkgd = permute(sum(abs(bsxfun(@minus,ims,bkgdim)),3),[1,2,4,3])/3;
pback = exp(-dbkgd/sigma_bkgd^2/2);
pback = pback / max(pback(:));
minv = exp(-fg_thresh/sigma_bkgd^2/2);
maxv = exp(-bg_thresh/sigma_bkgd^2/2);
pback = min(1,max(0,(pback-minv) / (maxv-minv)));
pfore = 1-pback;

%% make a connectivity graph for computing distances within foreground
% pixels
nr = ylim(2)-ylim(1)+1;
nc = xlim(2)-xlim(1)+1;
np = nr*nc;

[cgrid,rgrid] = meshgrid(1:nc,1:nr);
rconn1 = repmat(rgrid,[1,1,4]);
cconn1 = repmat(cgrid,[1,1,4]);
rconn2 = nan(size(rconn1));
cconn2 = nan(size(cconn1));

rconn2(:,:,1) = rgrid+1; % above
cconn2(:,:,1) = cgrid;
rconn2(:,:,2) = rgrid-1; % below
cconn2(:,:,2) = cgrid;
rconn2(:,:,3) = rgrid; % right
cconn2(:,:,3) = cgrid+1;
rconn2(:,:,4) = rgrid; % left
cconn2(:,:,4) = cgrid-1;
badidx = rconn2 < 1 | rconn2 > nr | cconn2 < 1 | cconn2 > nc;
rconn1(badidx) = [];
cconn1(badidx) = [];
rconn2(badidx) = [];
cconn2(badidx) = [];
conn1 = sub2ind([nr,nc],rconn1,cconn1);
conn2 = sub2ind([nr,nc],rconn2,cconn2);

%isdy = conn2 == conn1+1 | conn2 == conn1-1;
%w = min(maxdist,-log(pfore)+dist_epsilon)/maxdist;
w = reshape((pback+dist_epsilon)/(1+dist_epsilon),[np,nframes]);

%% assign pixels to flies

% probability that each pixel belongs to each fly
pfly = zeros(nr*nc,trx.nflies,nframes);
%allflies = [flygroups{:}];

for ii = 1:nframes,
  t = ts(ii);
  i = t-ts(1)+1;
  for fly = 1:trx.nflies,
    if t > trx(fly).endframe || t < trx(fly).firstframe,
      continue;
    end
    mu = [x(fly,i)-frac_a_back*a(fly,i)*cos(theta(fly,i)),...
      y(fly,i)-frac_a_back*a(fly,i)*sin(theta(fly,i))];
%     mu = [trx(fly).x(j)-frac_a_back*trx(fly).a(j)*cos(trx(fly).theta(j)),...
%       trx(fly).y(j)-frac_a_back*trx(fly).a(j)*sin(trx(fly).theta(j))];

    s = sub2ind([nr,nc],min(nr,max(1,round(mu(2)))),...
      min(nc,max(1,round(mu(1)))));
    wcurr = w(conn2,ii);
    S = axes2cov(a(fly,i),b(fly,i),theta(fly,i));
    diffs = bsxfun(@minus,[cgrid(:),rgrid(:)],mu);
    c = chol(S);
    dmah = sum((diffs/c).^2',1); %#ok<UDIM>
%     ratio = S(1,1)/S(2,2);
%     rationorm = 2 / (ratio+1);
%     wcurr(isdy) = wcurr(isdy)*ratio*rationorm;
%     wcurr(~isdy) = wcurr(~isdy)*rationorm;
    G = sparse(conn1,conn2,wcurr,np,np);
    d = graphshortestpath(G,s,'Directed',true);
    %pfly(:,fly,i) = exp(-d.^2/(trx(fly).a(j)*dist_epsilon)^2);
    pfly(:,fly,ii) = exp(-(1-wmah)*d.^2/(a(fly,i)*dist_epsilon)^2 - wmah*dmah);
  end
end
Z = sum(pfly,2);
Z(Z==0) = 1;
pfly = bsxfun(@rdivide,pfly,Z);

gray_ims = reshape(gray_ims,[nr*nc,nframes]);
pfore = reshape(pfore,[nr*nc,nframes]);
psomefly = reshape(1-prod(1-pfly(:,[flygroups{:}],:),2),[nr*nc,nframes]);
pforefly = pfore.*psomefly;
%pfore_any = 1 - prod(1-pfore.*psomefly,2);

%% color the images

im = zeros(nr*nc,3,nframes);
if flipim,
  if intensityoff > 0,
    gray_ims = intensityoff+(1 - gray_ims)*(1-intensityoff);
  end
end
for i = 1:nframes,
  imcurr = zeros(nr*nc,3);
  for flygroupi = 1:nflygroups,
    flygroup = flygroups{flygroupi};
    pcurr = 1 - prod(1-pfly(:,flygroup,i),2);
    colorcurr = colors(i,:,flygroupi);
    imcurr = imcurr + bsxfun(@times,colorcurr,pcurr);
  end
  if ~isempty(solidbkgdcolor),
    imcurr = bsxfun(@times,imcurr,gray_ims(:,i).*pforefly(:,i)) + ...
      bsxfun(@times,1-pforefly(:,i),repmat(reshape(solidbkgdcolor,[1,3]),[nr*nc,1]));
  else
    imcurr = bsxfun(@plus,bsxfun(@times,imcurr,gray_ims(:,i).*pforefly(:,i)),...
      (1-pforefly(:,i)).*reshape(bkgdim,[nr*nc,1]));
  end
  im(:,:,i) = imcurr;
end
finalim = reshape(im,[nr,nc,3,nframes]);

final_pforefly = reshape(pforefly,[nr,nc,nframes]);

%%

if isempty(hfig),
  hfig = figure;
end

figure(hfig);
clf;
hax = createsubplots(1,numel(ts),.01);

idxpos = predictions{mainfly}(ts(1):ts(end))>0;
idxneg = predictions{mainfly}(ts(1):ts(end))==0;
idxunlabeled = predictions{mainfly}(ts(1):ts(end))<0;

imid = ceil(numel(ts)/2);
try
  timestamps = trx(mainfly).timestamps(ts+trx(mainfly).off)-...
    trx(mainfly).timestamps(ts(imid)+trx(mainfly).off);
catch
  timestamps = [0,cumsum(trx(mainfly).dt)];
  timestamps = timestamps(ts+trx(mainfly).off)-...
    timestamps(ts(imid)+trx(mainfly).off);
end

for i = 1:numel(ts),
  t = ts(i);
  image(finalim(:,:,:,i),'Parent',hax(i));
  axis(hax(i),'image','off');
  hold(hax(i),'on');
  plot(hax(i),x(mainfly,:),y(mainfly,:),'k.-');
  plot(hax(i),x(mainfly,idxunlabeled),y(mainfly,idxunlabeled),'.','Color','w');
  plot(hax(i),x(mainfly,idxneg),y(mainfly,idxneg),'.','Color',colorneg);
  plot(hax(i),x(mainfly,idxpos),y(mainfly,idxpos),'.','Color',colorpos);
  if predictions{mainfly}(t)>0
    colorcurr = colorpos;
  elseif predictions{mainfly}(t)==0,
    colorcurr = colorneg;
  else
    colorcurr = [1,1,1];
  end
  plot(hax(i),x(mainfly,t-ts(1)+1),y(mainfly,t-ts(1)+1),'o','color',colorcurr,'markerfacecolor',colorcurr);
  text(1,1,sprintf('t = %.2fs',timestamps(i)),'HorizontalAlignment','left','VerticalAlignment','top','Parent',hax(i));
end

axes(hax(end));
colormap(colors_all(:,:,1)*.75);
tlim = timestamps(end)-timestamps(1);
set(hax,'CLim',[0,tlim]);
hcb = colorbar('East');
tticks = [0,tlim];
realylim = get(hcb,'YLim');
yticks = realylim(1)+(tticks/tlim)*(realylim(2)-realylim(1));
set(hcb,'YTick',yticks,'YTickLabel',tticks');

if isempty(figpos),
  truesize;
else
  set(hfig,'Units','pixels','Position',figpos);
end


text:
function jac = p07_jac ( option, nvar, x )

%*****************************************************************************80
%
%% P07_JAC evaluates the jacobian for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the jacobian.
%
%    Output, real JAC(NVAR-1,NVAR), the jacobian matrix evaluated
%    at X.  The NVAR-th row is not set by this routine.
%
  jac = zeros ( nvar, nvar );

  for i = 1 : nvar - 1
    jac(i,i) = 100.0 * ( 1.0 - x(i) * x(i) ) ...
      / ( 1.0 + x(i) + x(i) * x(i) )^2;
  end

  jac(1,1) = jac(1,1) + 2.0;
  jac(1,2) = jac(1,2) - 1.0;
  jac(1,nvar) = jac(1,nvar) - 1.0;

  for i = 2 : nvar-2
    jac(i,i-1) = jac(i,i-1) - 1.0;
    jac(i,i) = jac(i,i) + 3.0;
    jac(i,i+1) = jac(i,i+1) - 1.0;
    jac(i,nvar) = jac(i,nvar) - 1.0;
  end

  jac(nvar-1,nvar-2) = jac(nvar-1,nvar-2) - 1.0;
  jac(nvar-1,nvar-1) = jac(nvar-1,nvar-1) + 2.0;
  jac(nvar-1,nvar) = jac(nvar-1,nvar) - 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function drgui
%DRGUI Runs the graphical user interface of the DR toolbox
%
%   drgui
%
% Runs the graphical user interface of the dimensionality reduction toolbox.

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    drtool;

repo: cluster


text:
function sF = log(sF, varargin)
% log of a function
%
% Syntax
%   sF = log(sF)
%   sF = log(sF, 'bandwidth', bandwidth)
%
% Input
%  sF - @S2FunHarmonic
%
% Output
%  sF - @S2FunHarmonic
%
% Options
%  bandwidth - minimal degree of the spherical harmonic
%

%sF = sF.quadrature(@(v) log(sF.eval(v)),varargin{:});

sF = S2FunHandle(@(v) log(sF.eval(v)));

end


text:
function varargout = pm_pad(varargin)
% Pads a (partially) unwrapped phasemap such that the phase
% at a non-unwrapped location is a weighted average of unwrapped
% neighbouring phase-values.
% FORMAT [pm,wmap] = pm_pad(pm,wmap,kernel)
%
% Input:
% pm     : 2 or 3D phasemap where some voxels have been unwrapped
%          and some not.
% wmap   : Wrap-map, where a non-zero value indicates corresponding
%          phase-value in pm has been unwrapped.
% kernel : kernel used to generate a weighted average of surrounding
%          voxels.
%
% Output:
% pm     : Same as pm in, but where some previously unwrapped
%          phase-values have now been replaced.
% wmap   : Same as wmap in, but where values that was replaced
%          by weighted average in pm have now been set.
%__________________________________________________________________________
% Copyright (C) 2008-2015 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: pm_pad.m 6501 2015-07-17 14:32:09Z spm $

error('mex-function pm_pad.c not compiled');

author: spm


text:
function result = regularizeGray(input,nodes,edges,lambda,tol,maxiters)
%
% result = regularizeMap(view,input,lambda)
%
% input: data from one scan, e.g.,
%    co(:,1) or  z = amp(:,1).*exp(i*ph(:,1))
% can have NaNs indicating no data (e.g., for voxels below cothresh)
%
% nodes:  8xN array of (x,y,z,num_edges,edge_offset,layer,dist,pqindex).
%
% edges:  1xM array of node indices.  The edge_offset of
%    each node points into the starting location of its set
%    of edges.
% where N, M are the number of nodes, edges in the graph.
%
% lambda: smoothness coefficient
%
% result^{i+1}[x] = c2[x] s^i[x]               data missing
% 		  = c1[x] (input[x] + lambda s^i[x])     otherwise
%
% c_1[x] = 1/(1 + numNeighbors * lambda)
% c_2[x] = 1/numNeighbors
% s[x] = sumNeighbors

% Defaults for tol and maxiters
%
if ~exist('tol','var')
  tol = 1e-2;
end
if ~exist('maxiters','var')
  maxiters = 100;
end

% Get indices for missing data
%
NaNs=find(isnan(input));
notNaNs=find(~isnan(input));

% Initialize result
%

numNodes = length(input);
result = zeros(size(input));
result(notNaNs) = input(notNaNs);
newresult = zeros(size(input));

% Get numNeighbors and compute c1 and c2
%
numNeighbors = double(nodes(4,:));
edges        = double(edges);

edgeOffsets = double(nodes(5,:));
c1 = 1 ./ (1 + lambda*numNeighbors);
c2 = 1 ./ numNeighbors;
c2(find(numNeighbors==0)) = 0;
c = zeros(size(input));
c(notNaNs) = c1(notNaNs);
c(NaNs) = c2(NaNs);

% Initialize iterations
%
inputSD=std(input(notNaNs));
snr=Inf;
iter=0;

waitHandle = mrvWaitbar(0,'Smoothing data.  Please wait...');
while ((snr>tol) & (iter<maxiters))
   mrvWaitbar(iter/maxiters)
   iter = iter+1;

   % Compute sumNeighbors
   sumNeighbors = sumOfNeighbors(real(result),edges,edgeOffsets,numNeighbors) + ...
      j*sumOfNeighbors(imag(result),edges,edgeOffsets,numNeighbors);

   newresult(NaNs) = c(NaNs) .* sumNeighbors(NaNs);
   newresult(notNaNs) = c(notNaNs) .* (input(notNaNs) + lambda*sumNeighbors(notNaNs));

   % Compute snr (stopping criterion)

   snr=sqrt(mean(abs(newresult-result).^2))/inputSD;
   result=newresult;
end
close(waitHandle);

if (iter >= maxiters)
   disp(['Warning: maximum number of iterations exceeded: ',num2str(maxiters)]);
end

return;

%%% Test code

mrLoadRet

% open volume window
% switch to gray mode
% load anatomy
% view phase

% coronal slice 60

curScan = getCurScan(VOLUME{1});
amp = VOLUME{1}.amp(:,curScan);
co = VOLUME{1}.co(:,curScan);
ph = VOLUME{1}.ph(:,curScan);
z = co.*exp(i*ph);

zthresh = z;
belowThreshIndices = find(co<.2);
zthresh(belowThreshIndices) = NaN;

result = regularizeGray(zthresh,VOLUME{1}.nodes,VOLUME{1}.edges,.1,1e-2);
newPh = angle(result);
newPh(newPh<0) = newPh(newPh<0)+pi*2;
VOLUME{1}.ph(:,curScan) = newPh;
VOLUME{1}.ph(:,curScan) = ph;

%%% Compile MEX file for sumOfNeighbors
% mcc -ir sumOfNeighbors

%%% Test sumOfNeighbors

nodes = VOLUME{1}.nodes;
edges = VOLUME{1}.edges;
numNeighbors = nodes(4,:);
sumNeighbors = sumOfNeighbors(real(z),nodes,edges,numNeighbors);
figure(3)
plot(sumNeighbors,'.')


text:
function value = year_length_months_coptic ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_COPTIC returns the number of months in a Coptic year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  value = 12;

  return
end


text:
function [res] = mne_transpose_named_matrix(mat)
%
% [res] = mne_transpose_named_matrix(mat)
%
% Transpose a named matrix
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.2  2006/04/23 15:29:41  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/18 23:21:22  msh
%   Added mne_transform_source_space_to and mne_transpose_named_matrix
%
%

res.nrow = mat.ncol;
res.ncol = mat.nrow;
res.row_names = mat.col_names;
res.col_names = mat.row_names;
res.data      = mat.data';

return;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x,fval,exitflag,info,Opt] = opti_mintprog(f,A,b,Aeq,beq,lb,ub,int,sos,opts)
%OPTI_MINTPROG Solve a MILP using an OPTI MILP Solver (Matlab-Like Overload)
%
%   [x,fval,exitflag,info] = opti_mintprog(f,A,b,Aeq,beq,lb,ub,int) solves
%   the linear program min f'x where A,b are the inequality constraints,
%   Aeq,beq are the equality constraints, lb,ub are the decision variable
%   bounds and int is a string of integer variables ('C', 'I', 'B').
%
%   [x,fval,exitflag,info] = opti_minprog(f,...,int,sos) allows Special
%   Ordered Sets (SOS) to be specified via the sos structure. Fields must
%   include 'type' ('1' or '2'), 'index' (vector) and 'weight' (vector).
%   Multiple SOS can be supplied as cell arrays.
%
%   [x,fval,exitflag,info] = opti_mintprog(f,...,sos,opts) allows the user
%   to specify optiset options. This includes specifying a solver via the
%   'solver' field of optiset.
%
%   [x,...,info,Opt] = opti_mintprog(f,...) returns the internally built
%   OPTI object.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

% Handle missing arguments
if nargin < 10, opts = optiset; end
if nargin < 9, sos = []; end
if nargin < 8, int = []; end
if nargin < 7, ub = []; end
if nargin < 6, lb = []; end
if nargin < 5, beq = []; end
if nargin < 4, Aeq = []; end
if nargin < 3, error('You must supply at least 3 arguments to opti_mintprog'); end

%Check SOS
if(~isempty(sos))
    if(~isfield(sos,'type')), error('The SOS structure must have the field type!'); end
    if(~isfield(sos,'index')), error('The SOS structure must have the field index!'); end
    if(~isfield(sos,'weight')), error('The SOS structure must have the field weight!'); end
end

%Build OPTI Object
Opt = opti('f',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'int',int,'sos',sos,'options',opts);

%Solve
[x,fval,exitflag,info] = solve(Opt);


text:
function featAugm = get_augm_spatial_features_diff_neighbour_locref(feat)

% relative coord of 2 detections
delta = feat(:, 1:2);
a = compute_angle(delta(:,1), delta(:,2));

delta_forward = feat(:, 5:6);
a_forward = compute_angle(delta_forward(:,1), delta_forward(:,2));

delta1 = delta - delta_forward;
dist1 = sqrt(delta1(:,1).^2 + delta1(:,2).^2);
a1 = wrap_angle(a - a_forward);
abs_a1 = abs(a1);
%abs_a1(:) = 0;

% now the same for backward
delta = feat(:, 3:4);
a = compute_angle(-delta(:,1), -delta(:,2));

delta_backward = feat(:, 7:8);
a_backward = compute_angle(delta_backward(:,1), delta_backward(:,2));

delta2 = delta - delta_backward;
dist2 = sqrt(delta2(:,1).^2 + delta2(:,2).^2);
a2 = wrap_angle(a - a_backward);
abs_a2 = abs(a2);

%featAugm = cat(2, dist, exp(-dist), dx, dy, abs(dx), abs(dy), dxdy, dx.^2, dy.^2, abs_a, a.^2, exp(-abs_a));

featAugm = cat(2, dist1, abs_a1, dist2, abs_a2);
featAugm = cat(2, featAugm, exp(-featAugm));

end

function angle = compute_angle(deltaX, deltaY)
angle = atan2(deltaY,deltaX);
angle = wrapMinusPiPifast(angle);
assert(all(angle <= pi) && all(angle >= -pi));
end

function a = wrap_angle(a)
larger = a > pi;
smaller = a < -pi;
a(larger)  = a(larger) - 2*pi;
a(smaller) = a(smaller)+ 2*pi;
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [ shape2D, global_params, local_params, final_lhood, landmark_lhoods, view_used ] = Fitting_from_bb( Image, DepthImage, bounding_box, PDM, patchExperts, clmParams, varargin)
%FITTING Summary of this function goes here
%   Detailed explanation goes here

    % the bounding box format is [minX, minY, maxX, maxY];

    % the mean model shape
    M = PDM.M;

    num_points = numel(M) / 3;

    if(any(strcmp(varargin,'orientation')))
        orientation = varargin{find(strcmp(varargin, 'orientation'))+1};
        rot = Euler2Rot(orientation);
    else
        rot = eye(3);
        orientation = [0;0;0];
    end

    rot_m = rot * reshape(M, num_points, 3)';
    width_model = max(rot_m(1,:)) - min(rot_m(1,:));
    height_model = max(rot_m(2,:)) - min(rot_m(2,:));

    a = (((bounding_box(3) - bounding_box(1)) / width_model) + ((bounding_box(4) - bounding_box(2))/ height_model)) / 2;

    tx = (bounding_box(3) + bounding_box(1))/2;
    ty = (bounding_box(4) + bounding_box(2))/2;

    % correct it so that the bounding box is just around the minimum
    % and maximum point in the initialised face
    tx = tx - a*(min(rot_m(1,:)) + max(rot_m(1,:)))/2;
    ty = ty - a*(min(rot_m(2,:)) + max(rot_m(2,:)))/2;

    % visualisation of the initial state
    %hold off;imshow(Image);hold on;plot(a*rot_m(1,:)+tx, a*rot_m(2,:)+ty,'.r');hold on;rectangle('Position', [bounding_box(1), bounding_box(2), bounding_box(3)-bounding_box(1), bounding_box(4)-bounding_box(2)]);
    global_params = [a, 0, 0, 0, tx, ty]';
    global_params(2:4) = orientation;

    local_params = zeros(numel(PDM.E), 1);

    if(any(strcmp(varargin,'gparam')))
        global_params = varargin{find(strcmp(varargin, 'gparam'))+1};
    end

    if(any(strcmp(varargin,'lparam')))
        local_params = varargin{find(strcmp(varargin, 'lparam'))+1};
    end

    scale = clmParams.startScale;

    if(size(Image, 3) == 1)
        GrayImage = Image;
    else
        GrayImage = rgb2gray(Image);
    end

    [heightImg, widthImg] = size(GrayImage);

    % Some predefinitions for faster patch extraction
    [xi, yi] = meshgrid(0:widthImg-1,0:heightImg-1);
    xi = double(xi);
    yi = double(yi);

    GrayImageDb = double(GrayImage);

    clmParams_old = clmParams;

    % In case there are no iterations initialize thes to empty
    view = GetView(patchExperts(scale).centers, global_params(2:4));
    final_lhood = 0;
    [shape2D] = GetShapeOrtho(M, PDM.V, local_params, global_params);
    landmark_lhoods = zeros(size(shape2D,1),1);

    % multi iteration refinement using NU-RLMS in each one
    for i=1:clmParams.numPatchIters

        current_patch_scaling = patchExperts(scale).trainingScale;
        visibilities = patchExperts(scale).visibilities;

        view = GetView(patchExperts(scale).centers, global_params(2:4));

        % The shape fitting is performed in the reference frame of the
        % patch training scale
        refGlobal = [current_patch_scaling, 0, 0, 0, 0, 0]';

        % the reference shape
        refShape = GetShapeOrtho(M, PDM.V, local_params, refGlobal);

        % shape around which the patch experts will be evaluated in the original image
        [shape2D] = GetShapeOrtho(M, PDM.V, local_params, global_params);
        shape2D_img = shape2D(:,1:2);

        % Create transform using a slightly modified version of Kabsch that
        % takes scaling into account as well, in essence we get a
        % similarity transform from current estimate to reference shape
        [A_img2ref, T_img2ref, ~, ~] = AlignShapesWithScale(shape2D_img(:,1:2),refShape(:,1:2));

        % Create a transform, from shape in image to reference shape
        T = maketform('affine', [A_img2ref;T_img2ref]);

        shape_2D_ref = tformfwd(T, shape2D_img);

        % transform the current shape to the reference one, so we can
        % interpolate
        shape2D_in_ref = (A_img2ref * shape2D_img')';

        sideSizeX = (clmParams.window_size(i,1) - 1)/2;
        sideSizeY = (clmParams.window_size(i,2) - 1)/2;

        patches = zeros(size(shape2D_in_ref,1), clmParams.window_size(i,1) * clmParams.window_size(i,2));

        Ainv = inv(A_img2ref);

        % extract patches on which patch experts will be evaluted
        for l=1:size(shape2D_in_ref,1)
            if(visibilities(view,l))

                xs = (shape2D_in_ref(l,1)-sideSizeX):(shape2D_in_ref(l,1)+sideSizeX);
                ys = (shape2D_in_ref(l,2)-sideSizeY):(shape2D_in_ref(l,2)+sideSizeY);

                [xs, ys] = meshgrid(xs, ys);

                pairs = [xs(:), ys(:)];

                actualLocs = (Ainv * pairs')';

                actualLocs(actualLocs(:,1) < 0,1) = 0;
                actualLocs(actualLocs(:,2) < 0,2) = 0;
                actualLocs(actualLocs(:,1) > widthImg - 1,1) = widthImg - 1;
                actualLocs(actualLocs(:,2) > heightImg - 1,2) = heightImg - 1;

                [t_patch] = interp2_mine(xi, yi, GrayImageDb, actualLocs(:,1), actualLocs(:,2), 'bilinear');
                t_patch = reshape(t_patch, size(xs));

                patches(l,:) = t_patch(:);

            end
        end

        % Calculate patch responses, either SVR or CCNF
        if(strcmp(patchExperts(scale).type, 'SVR'))
            responses = PatchResponseSVM_multi_modal( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), patchExperts(scale).normalisationOptionsCol, clmParams, clmParams.window_size(i,:));
        elseif(strcmp(patchExperts(scale).type, 'CCNF'))
            responses = PatchResponseCCNF( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), patchExperts(scale), clmParams.window_size(i,:));
        elseif(strcmp(patchExperts(scale).type, 'CEN'))

            % responses_o = PatchResponseCEN( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), patchExperts(scale), clmParams.window_size(i,:));

            % A faster (and very slightly less accurate) version of patch
            % response computation
            responses = PatchResponseCEN_sparse( patches, patchExperts(scale).patch_experts(view,:), visibilities(view,:), clmParams.window_size(i,:), patchExperts(scale).normalisationOptionsCol.patchSize);
        end

        % If a depth image is provided compute patch experts around it as
        % well (unless it's the final iteration)
        if(~isempty(DepthImage) && (i ~= clmParams.numPatchIters))

            % Extracting the depth patches here
            patches_depth = zeros(size(shape2D_in_ref,1), clmParams.window_size(i,1) * clmParams.window_size(i,2));

            % extract patches on which patch experts will be evaluted
            for l=1:size(shape2D_in_ref,1)
                if(visibilities(view,l))

                    xs = (shape2D_in_ref(l,1)-sideSizeX):(shape2D_in_ref(l,1)+sideSizeX);
                    ys = (shape2D_in_ref(l,2)-sideSizeY):(shape2D_in_ref(l,2)+sideSizeY);

                    [xs, ys] = meshgrid(xs, ys);

                    pairs = [xs(:), ys(:)];

                    actualLocs = (Ainv * pairs')';

                    actualLocs(actualLocs(:,1) < 1,1) = 1;
                    actualLocs(actualLocs(:,2) < 1,2) = 1;
                    actualLocs(actualLocs(:,1) > widthImg,1) = widthImg;
                    actualLocs(actualLocs(:,2) > heightImg,2) = heightImg;

                    % use nearest neighbour interpolation as bilinear would
                    % produce artefacts in depth image (when missing data
                    % is there)
                    [t_patch] = interp2_mine(xi, yi, DepthImage, actualLocs(:,1), actualLocs(:,2), 'nearest');
                    t_patch = reshape(t_patch, size(xs));

                    patches_depth(l,:) = t_patch(:);

                end
            end

            old_mm = clmParams.use_multi_modal;
            clmParams.use_multi_modal = 0;
            responses_depth = PatchResponseSVM_multi_modal( patches_depth, patchExperts(scale).patch_experts_depth(view,:), visibilities(view,:), patchExperts(scale).normalisationOptionsDepth, clmParams, clmParams.window_size(i,:));
            clmParams.use_multi_modal = old_mm;

            % Combining the patch responses from different channels here
            for l=1:size(shape2D_in_ref,1)
                responses{l} = responses{l} + responses_depth{l};
            end
        end

        % the better the correlation in training the more reliable the feature
        % the reliabilities are independent for every modality in SVR, so
        % combine them (also correlation is inverse to variance)
        if(strcmp(patchExperts(scale).type, 'SVR'))
            if(clmParams.use_multi_modal)
                reliabilities = patchExperts(scale).patch_experts{1,1}(end).correlations{1};
            else
                reliabilities = patchExperts(scale).patch_experts{1,1}(1).correlations{1};
            end
        else
            % for CCNF the modalities work together
            reliabilities = patchExperts(scale).correlations;
        end
        reliabilities = reliabilities(view,:);

        % deal with the fact that params might be different for different
        % scales
        if(numel(clmParams_old.regFactor) > 1)
            clmParams.regFactor = clmParams_old.regFactor(scale);
        end
        if(numel(clmParams_old.sigmaMeanShift) > 1)
            clmParams.sigmaMeanShift = clmParams_old.sigmaMeanShift(scale);
        end
        if(numel(clmParams_old.tikhonov_factor) > 1)
            clmParams.tikhonov_factor = clmParams_old.tikhonov_factor(scale);
        end

        % The actual NU-RLMS step

        % first the rigid transform
        [global_params, local_params] = NU_RLMS(global_params, local_params, PDM, responses, visibilities, view, reliabilities, shape2D_img, T, true, clmParams, []);

        % second the combined transform
        [global_params, local_params, final_lhood, landmark_lhoods] = ...
            NU_RLMS(global_params, local_params, PDM, responses, visibilities, view, reliabilities, shape2D_img, T, false, clmParams, []);

        % Clamp orientation and make sure it doesn't get out of hand
        orientation = global_params(2:4);
        orientation(orientation < -pi/2) = -pi/2;
        orientation(orientation > pi/2) = pi/2;
        global_params(2:4) = orientation;

        % move up a scale if possible
        if(clmParams.useMultiScale && scale ~= numel(patchExperts))

            % only go up a scale if we don't need to upsample
            if(0.9 * patchExperts(scale+1).trainingScale < global_params(1))
                scale = scale + 1;
            else
                break;
            end
        end
    end

    % the view in last iteration
    view_used = view;

    % See how good the tracking was in the end
    [shape2D] = GetShapeOrtho(M, PDM.V, local_params, global_params);

    % Moving to matlab format
    shape2D = shape2D(:,1:2) + 1;

end

function [id] = GetView(centers, rotation)

    [~,id] = min(sum((centers * pi/180 - repmat(rotation', size(centers,1), 1)).^2,2));

end

author: TadasBaltrusaitis

repo: OpenFace