text:
function [Eft, Varft, lpyt, Eyt, Varyt] = gpmc_loopred(gp, x, y, varargin)
%GPMC_LOOPRED  Leave-one-out predictions with Gaussian Process MCMC approximation.
%
%  Description
%    [EFT, VARFT, LPYT, EYT, VARYT] = GPMC_LOOPRED(RECGP, X, Y, OPTIONS)
%    takes a Gaussian processes record structure RECGP (returned by
%    gp_mc) together with a matrix X of training inputs and vector
%    Y of training targets, and evaluates the leave-one-out
%    predictive distribution at inputs X and returns means EFT and
%    variances VARFT of latent variables, the logarithm of the
%    predictive densities PYT, and the predictive means EYT and
%    variances VARYT of observations at input locations X.
%
%    OPTIONS is optional parameter-value pair
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      is     - defines if importance sampling weighting is used 'on'
%               (default). If set to 'off', MCMC samples from the
%               full data posterior are used.
%
%    Given Gaussian likelihood or non-Gaussian likelihood and
%    latent method Laplace or EP, LOO-posterior given
%    hyperparameters is computed analytically or with analytic
%    approximation and LOO-posterior of the hyperparameters is
%    approximated using importance sampling. Optionally full data
%    posterior for hyperparameters can be used (by setting option
%    'is' to 'off'). Given non-Gaussian likelihood and latent
%    method MCMC, LOO-posterior of the hyperparameters and latent
%    values is approximated using importance sampling.
%
%  References:
%    Aki Vehtari and Jouko Lampinen (2002). Bayesian model
%    assessment and comparison using cross-validation predictive
%    densities. Neural Computation, 14(10):2439-2468.
%
%  See also
%   GP_LOOPRED, GP_MC, GP_PRED
%

% Copyright (c) 2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.FunctionName = 'GPMC_LOOPRED';
ip.addRequired('gp',@isstruct);
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('is', 'on', @(x) ismember(x,{'on' 'off'}))
ip.parse(gp, x, y, varargin{:});
z=ip.Results.z;
is=ip.Results.is;

if isfield(gp,'meanf') && ~isempty(gp.meanf)
  error('GPMC_LOOPRED: Mean functions not yet supported');
end

nmc=size(gp.jitterSigma2,1);
[n,nin]=size(x);

if strcmp(gp.type, 'PIC_BLOCK') || strcmp(gp.type, 'PIC')
  ind = gp.tr_index;           % block indeces for training points
  gp = rmfield(gp,'tr_index');
end

for i1=1:nmc
  % compute leave-one-out predictions for each hyperparameter sample
  % if latent method is MCMC, then these samples are for latent values, too
  Gp = take_nth(gp,i1);
  switch gp.type
    case 'FULL'

    case {'FIC' 'CS+FIC'}
      % Reformat the inducing inputs
      u = reshape(Gp.X_u,length(Gp.X_u)/nin,nin);
      Gp.X_u = u;

    case {'PIC' 'PIC_BLOCK'}
      % Reformat the inducing inputs
      u = reshape(Gp.X_u,length(Gp.X_u)/nin,nin);
      Gp.X_u = u;
      Gp.tr_index = ind;
  end
  if nargout <= 3
    [Efts(:,i1), Varfts(:,i1), lpyts(:,i1)] = gp_loopred(Gp, x, y, 'z', z);
  else
    [Efts(:,i1), Varfts(:,i1), lpyts(:,i1), Eyts(:,i1), Varyts(:,i1)] = gp_loopred(Gp, x, y, 'z', z);
  end
end

if isequal(is,'off')
  w=1/nmc;
  lw=-log(nmc);
else
  % log importance sampling weights
  lw=-lpyts;
  % normalize weights
  for i2=1:n
    % this works even when lw have large magnitudes
    lw(i2,:)=lw(i2,:)-sumlogs(lw(i2,:));
  end
  % importance sampling weights
  w=exp(lw);
  % check the effective sample size
  m_eff=1./sum(w.^2,2);
  if min(m_eff)<n/10
    mn=sum(m_eff<(n/10));
    if mn==1
      warning(sprintf('For %d data point the effective sample size in IS is less than n/10',mn))
    else
      warning(sprintf('For %d data points the effective sample size in IS is less than n/10',mn))
    end
  end
end

% IS-LOO predictions
Eft = sum(Efts.*w, 2);
Varft = sum(Varfts.*w,2) + sum(bsxfun(@minus,Efts,Eft).^2.*w,2);
if nargout > 2
  lpyt = log(sum(exp(lpyts+lw),2)); % same as lpy=log(sum(exp(lpys).*w,2));
  if nargout > 3
    Eyt = sum(Eyts.*w,2);
    Varyt = sum(Varyts.*w,2) + sum(bsxfun(@minus,Eyts,Eyt).^2.*w,2);
  end
end


text:
%demoL1iPotts_DecImp
% Reconstruction of a blurred jump-sparse signal from incomplete measurements under
% impulsive noise using the inverse L1-Potts functional

% load signal
groundTruth = loadPcwConst('sampleDec');
n = numel(groundTruth);

% create Gaussian kernel
K = convkernel('gaussian', 51, 6);

% create measurement matrix
Afull = spconvmatrix(K, numel(groundTruth));
fraction = 0.5;
idx = sort(randidx(n, fraction)) ;
A = Afull(idx, :);

% create blurred signal
fBlurry = A * groundTruth(:);

% impulsive noise (noiseFraction = number of pixels destroyed)
noiseFraction = 0.3;
ridx = randidx(numel(fBlurry), noiseFraction);
f = fBlurry;
f(ridx) =  (rand(size(ridx)));

% Solve inverse L1-Potts problem
gamma = 0.4;
u = minL1iPotts(f, gamma, A);

% show result
showPotts(f, u, groundTruth, 'L^1-iPotts')

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
%get_file_info2.m
%
%Replacement for get_file_info.m:
%   Numerous efficiency fixes, smaller filesize, easier readability,
%   and structural/organizational changes to keep returned info together

function [file_info] = get_file_info2(file, handles, override)
    [file_info.path, file_info.name, file_info.ext] = fileparts(file);

    try
        image0    = imread(file,2);
        file_info.type = 'multipage';
        int_prec=0;
    catch
        image0    = imread(file);
        file_info.type = 'singlepage';
    end

    class_type    = class(image0);
    file_info.bitdepth = log2(double(intmax(class_type))-double(intmin(class_type)));
    file_info.width    = size(image0,2);
    file_info.height   = size(image0,1);

    if ~isnan(str2double(get(handles.prec_edit,'String')))
        int_prec = str2double(get(handles.prec_edit,'String'));
        file_info.prec = ['%',num2str(int_prec),'.',num2str(int_prec),'d'];
    else
        length_name = length(file_info.name);
        for x = length_name:-1:1
            if isnan(str2double(file_info.name(x)))
                int_prec  = length_name-x;
                file_info.prec = ['%',num2str(int_prec),'.',num2str(int_prec),'d'];
                break
            end
        end
    end

    file_info.part_name = file_info.name(1:length(file_info.name)-int_prec);

    if get(handles.custom_frame,'Value')
        file_info.start = str2double(get(handles.from_edit,'String'));
        file_info.stop  = str2double(get(handles.to_edit,'String'));
    elseif ~get(handles.custom_frame,'Value')
        if strcmp(file_info.type,'multipage')
            file_info.start = 1;

            if ~exist('override','var') %Allow function inputs to override autodetect of end frame (saves on time)
                file_imf = imfinfo(file);
                file_info.stop  = length(file_imf);
            end
        elseif strcmp(file_info.type,'singlepage')
            file_info.start = 1;

            if ~exist('override','var') %Allow function inputs to override autodetect of end frame (saves on time)
                index = 1;
                while exist([file_info.path,'\',file_info.part_name,num2str(index,file_info.prec),file_info.ext],'file')
                    index = index + 1;
                end
                file_info.stop = index - 1;
            end
        end
    end

    if strcmp(file_info.type,'multipage') %Override prec and partial name for multipage images
        file_info.part_name = file_info.name;
        file_info.prec  = '%0.0d';
    end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef SettingsHomogenizedVarComputerFactory < handle

    methods (Access = public, Static)

        function obj = create(cParams)
            switch cParams.type
                case 'ByInterpolation'
                    obj = SettingsHomogenizedVarComputerFromInterpolation(cParams);
                case 'ByVademecum'
                    obj = SettingsHomogenizedVarComputerFromVademecum(cParams);
                otherwise
                    error('Invalid HomogenizerVarComputer type.');
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function vl_compilenn(varargin)
%VL_COMPILENN Compile the MatConvNet toolbox.
%   The `vl_compilenn()` function compiles the MEX files in the
%   MatConvNet toolbox. See below for the requirements for compiling
%   CPU and GPU code, respectively.
%
%   `vl_compilenn('OPTION', ARG, ...)` accepts the following options:
%
%   `EnableGpu`:: `false`
%      Set to true in order to enable GPU support.
%
%   `Verbose`:: 0
%      Set the verbosity level (0, 1 or 2).
%
%   `Debug`:: `false`
%      Set to true to compile the binaries with debugging
%      information.
%
%   `CudaMethod`:: Linux & Mac OS X: `mex`; Windows: `nvcc`
%      Choose the method used to compile the CUDA code. There are two
%      methods:
%
%      * The **`mex`** method uses the MATLAB MEX command with the
%        configuration file
%        `<MatConvNet>/matlab/src/config/mex_CUDA_<arch>.[sh/xml]`
%        This configuration file is in XML format since MATLAB 8.3
%        (R2014a) and is a Shell script for earlier versions. This
%        is, principle, the preferred method as it uses the
%        MATLAB-sanctioned compiler options.
%
%      * The **`nvcc`** method calls the NVIDIA CUDA compiler `nvcc`
%        directly to compile CUDA source code into object files.
%
%        This method allows to use a CUDA toolkit version that is not
%        the one that officially supported by a particular MATALB
%        version (see below). It is also the default method for
%        compilation under Windows and with CuDNN.
%
%   `CudaRoot`:: guessed automatically
%      This option specifies the path to the CUDA toolkit to use for
%      compilation.
%
%   `EnableImreadJpeg`:: `true`
%      Set this option to `true` to compile `vl_imreadjpeg`.
%
%   `ImageLibrary`:: `libjpeg` (Linux), `gdiplus` (Windows), `quartz` (Mac)
%      The image library to use for `vl_impreadjpeg`.
%
%   `ImageLibraryCompileFlags`:: platform dependent
%      A cell-array of additional flags to use when compiling
%      `vl_imreadjpeg`.
%
%   `ImageLibraryLinkFlags`:: platform dependent
%      A cell-array of additional flags to use when linking
%      `vl_imreadjpeg`.
%
%   `EnableCudnn`:: `false`
%      Set to `true` to compile CuDNN support.
%
%   `CudnnRoot`:: `'local/'`
%      Directory containing the unpacked binaries and header files of
%      the CuDNN library.
%
%   ## Compiling the CPU code
%
%   By default, the `EnableGpu` option is switched to off, such that
%   the GPU code support is not compiled in.
%
%   Generally, you only need a C/C++ compiler (usually Xcode, GCC or
%   Visual Studio for Mac, Linux, and Windows respectively). The
%   compiler can be setup in MATLAB using the
%
%      mex -setup
%
%   command.
%
%   ## Compiling the GPU code
%
%   In order to compile the GPU code, set the `EnableGpu` option to
%   `true`. For this to work you will need:
%
%   * To satisfy all the requirement to compile the CPU code (see
%     above).
%
%   * A NVIDIA GPU with at least *compute capability 2.0*.
%
%   * The *MATALB Parallel Computing Toolbox*. This can be purchased
%     from Mathworks (type `ver` in MATLAB to see if this toolbox is
%     already comprised in your MATLAB installation; it often is).
%
%   * A copy of the *CUDA Devkit*, which can be downloaded for free
%     from NVIDIA. Note that each MATLAB version requires a
%     particular CUDA Devkit version:
%
%     | MATLAB version | Release | CUDA Devkit |
%     |----------------|---------|-------------|
%     | 2013b          | 2013b   | 5.5         |
%     | 2014a          | 2014a   | 5.5         |
%     | 2014b          | 2014b   | 6.0         |
%
%     A different versions of CUDA may work using the hack described
%     above (i.e. setting the `CudaMethod` to `nvcc`).
%
%   The following configurations have been tested successfully:
%
%   * Windows 7 x64, MATLAB R2014a, Visual C++ 2010 and CUDA Toolkit
%     6.5 (unable to compile with Visual C++ 2013).
%   * Windows 8 x64, MATLAB R2014a, Visual C++ 2013 and CUDA
%     Toolkit 6.5.
%   * Mac OS X 10.9 and 10.10, MATLAB R2013a and R2013b, Xcode, CUDA
%     Toolkit 5.5.
%   * GNU/Linux, MATALB R2014a, gcc, CUDA Toolkit 5.5.
%
%   Furthermore your GPU card must have ComputeCapability >= 2.0 (see
%   output of `gpuDevice()`) in order to be able to run the GPU code.
%   To change the compute capabilities, for `mex` `CudaMethod` edit
%   the particular config file.  For the 'nvcc' method, compute
%   capability is guessed based on the GPUDEVICE output. You can
%   override it by setting the 'CudaArch' parameter (e.g. in case of
%   multiple GPUs with various architectures).
%
%   See also: [Compliling MatConvNet](../install.md#compiling),
%   [Compiling MEX files containing CUDA
%   code](http://mathworks.com/help/distcomp/run-mex-functions-containing-cuda-code.html),
%   `vl_setup()`, `vl_imreadjpeg()`.

% Copyright (C) 2014-15 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% Get MatConvNet root directory
root = fileparts(fileparts(mfilename('fullpath'))) ;
addpath(fullfile(root, 'matlab')) ;

% --------------------------------------------------------------------
%                                                        Parse options
% --------------------------------------------------------------------

opts.enableGpu        = false;
opts.enableImreadJpeg = true;
opts.enableCudnn      = false;
opts.imageLibrary = [] ;
opts.imageLibraryCompileFlags = {} ;
opts.imageLibraryLinkFlags = [] ;
opts.verbose          = 0;
opts.debug            = false;
opts.cudaMethod       = [] ;
opts.cudaRoot         = [] ;
opts.cudaArch         = [] ;
opts.defCudaArch      = [...
  '-gencode=arch=compute_20,code=\"sm_20,compute_20\" '...
  '-gencode=arch=compute_30,code=\"sm_30,compute_30\"'];
opts.cudnnRoot        = 'local/cudnn' ;
opts = vl_argparse(opts, varargin);

% --------------------------------------------------------------------
%                                                     Files to compile
% --------------------------------------------------------------------

arch = computer('arch') ;
if isempty(opts.imageLibrary)
  switch arch
    case 'glnxa64', opts.imageLibrary = 'libjpeg' ;
    case 'maci64', opts.imageLibrary = 'quartz' ;
    case 'win64', opts.imageLibrary = 'gdiplus' ;
  end
end
if isempty(opts.imageLibraryLinkFlags)
  switch opts.imageLibrary
    case 'libjpeg', opts.imageLibraryLinkFlags = {'-ljpeg'} ;
    case 'quartz', opts.imageLibraryLinkFlags = {'LDFLAGS=$LDFLAGS -framework Cocoa -framework ImageIO'} ;
    case 'gdiplus', opts.imageLibraryLinkFlags = {'-lgdiplus'} ;
  end
end

lib_src = {} ;
mex_src = {} ;

% Files that are compiled as CPP or CU depending on whether GPU support
% is enabled.
if opts.enableGpu, ext = 'cu' ; else, ext='cpp' ; end
lib_src{end+1} = fullfile(root,'matlab','src','bits',['data.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['datamex.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnconv.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnfullyconnected.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnsubsample.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnpooling.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnnormalize.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnbnorm.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnbias.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnconv.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnconvt.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnpool.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnnormalize.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnbnorm.' ext]) ;

% CPU-specific files
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','im2row_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','subsample_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','copy_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','pooling_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','normalize_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','bnorm_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','tinythread.cpp') ;

% GPU-specific files
if opts.enableGpu
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','im2row_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','subsample_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','copy_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','pooling_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','normalize_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','bnorm_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','datacu.cu') ;
end

% cuDNN-specific files
if opts.enableCudnn
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnconv_cudnn.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnbias_cudnn.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnpooling_cudnn.cu') ;
end

% Other files
if opts.enableImreadJpeg
  mex_src{end+1} = fullfile(root,'matlab','src', ['vl_imreadjpeg.' ext]) ;
  lib_src{end+1} = fullfile(root,'matlab','src', 'bits', 'impl', ['imread_' opts.imageLibrary '.cpp']) ;
end

% --------------------------------------------------------------------
%                                                   Setup CUDA toolkit
% --------------------------------------------------------------------

if opts.enableGpu
  opts.verbose && fprintf('%s: * CUDA configuration *\n', mfilename) ;

  % Find the CUDA Devkit
  if isempty(opts.cudaRoot), opts.cudaRoot = search_cuda_devkit(opts) ; end
  opts.verbose && fprintf('%s:\tCUDA: using CUDA Devkit ''%s''.\n', ...
                          mfilename, opts.cudaRoot) ;

  opts.nvccPath = fullfile(opts.cudaRoot, 'bin', 'nvcc') ;
  switch arch
    case 'win64', opts.cudaLibDir = fullfile(opts.cudaRoot, 'lib', 'x64') ;
    case 'maci64', opts.cudaLibDir = fullfile(opts.cudaRoot, 'lib') ;
    case 'glnxa64', opts.cudaLibDir = fullfile(opts.cudaRoot, 'lib64') ;
    otherwise, error('Unsupported architecture ''%s''.', arch) ;
  end

  % Activate the CUDA Devkit
  cuver = activate_nvcc(opts.nvccPath) ;
  opts.verbose && fprintf('%s:\tCUDA: using NVCC ''%s'' (%d).\n', ...
                          mfilename, opts.nvccPath, cuver) ;

  % Set the CUDA arch string (select GPU architecture)
  if isempty(opts.cudaArch), opts.cudaArch = get_cuda_arch(opts) ; end
  opts.verbose && fprintf('%s:\tCUDA: NVCC architecture string: ''%s''.\n', ...
                          mfilename, opts.cudaArch) ;
end

if opts.enableCudnn
  opts.cudnnIncludeDir = fullfile(opts.cudnnRoot, 'include') ;
  switch arch
    case 'win64', opts.cudnnLibDir = fullfile(opts.cudnnRoot, 'lib', 'x64') ;
    case 'maci64', opts.cudnnLibDir = fullfile(opts.cudnnRoot, 'lib') ;
    case 'glnxa64', opts.cudnnLibDir = fullfile(opts.cudnnRoot, 'lib64') ;
    otherwise, error('Unsupported architecture ''%s''.', arch) ;
  end
end

% --------------------------------------------------------------------
%                                                     Compiler options
% --------------------------------------------------------------------

% Build directories
mex_dir = fullfile(root, 'matlab', 'mex') ;
bld_dir = fullfile(mex_dir, '.build');
if ~exist(fullfile(bld_dir,'bits','impl'), 'dir')
  mkdir(fullfile(bld_dir,'bits','impl')) ;
end

% Compiler flags
flags.cc = {} ;
flags.link = {} ;
if opts.verbose > 1
  flags.cc{end+1} = '-v' ;
  flags.link{end+1} = '-v' ;
end
if opts.debug
  flags.cc{end+1} = '-g' ;
else
  flags.cc{end+1} = '-DNDEBUG' ;
end
if opts.enableGpu, flags.cc{end+1} = '-DENABLE_GPU' ; end
if opts.enableCudnn,
  flags.cc{end+1} = '-DENABLE_CUDNN' ;
  flags.cc{end+1} = ['-I' opts.cudnnIncludeDir] ;
end
flags.link{end+1} = '-lmwblas' ;
switch arch
  case {'maci64', 'glnxa64'}
  case {'win64'}
    % VisualC does not pass this even if available in the CPU architecture
    flags.cc{end+1} = '-D__SSSE3__' ;
end

if opts.enableImreadJpeg
  flags.cc = horzcat(flags.cc, opts.imageLibraryCompileFlags) ;
  flags.link = horzcat(flags.link, opts.imageLibraryLinkFlags) ;
end

if opts.enableGpu
  flags.link{end+1} = ['-L' opts.cudaLibDir] ;
  flags.link{end+1} = '-lcudart' ;
  flags.link{end+1} = '-lcublas' ;
  switch arch
    case {'maci64', 'glnxa64'}
      flags.link{end+1} = '-lmwgpu' ;
    case 'win64'
      flags.link{end+1} = '-lgpu' ;
  end
  if opts.enableCudnn
    flags.link{end+1} = ['-L' opts.cudnnLibDir] ;
    flags.link{end+1} = '-lcudnn' ;
  end
end

% For the MEX command
flags.link{end+1} = '-largeArrayDims' ;
flags.mexcc = flags.cc ;
flags.mexcc{end+1} = '-largeArrayDims' ;
flags.mexcc{end+1} = '-cxx' ;
if strcmp(arch, 'maci64') && opts.enableGpu && cuver < 70000
  % CUDA prior to 7.0 on Mac require GCC libstdc++ instead of the native
  % Clang libc++. This should go away in the future.
  flags.mexcc{end+1} = 'CXXFLAGS=$CXXFLAGS -stdlib=libstdc++' ;
  flags.link{end+1} = 'LDFLAGS=$LDFLAGS -stdlib=libstdc++' ;
  if  ~verLessThan('matlab', '8.5.0')
    % Complicating matters, MATLAB 8.5.0 links to Clang's libc++ by
    % default when linking MEX files overriding the option above. We
    % force it to use GCC libstdc++
    flags.link{end+1} = 'LINKLIBS=$LINKLIBS -L"$MATLABROOT/bin/maci64" -lmx -lmex -lmat -lstdc++' ;
  end
end
if strcmp(arch, 'maci64') && opts.enableGpu
  % Mac OS X 10.11 disables LD_LIBRARY_PATH for security reason. To
  % address this issue, we need to rpath the required CUDA and
  % cuDNN libraries.
  flags.link{end+1} = sprintf('LDFLAGS=$LDFLAGS -Wl,-rpath -Wl,"%s"', opts.cudaLibDir) ;
  if opts.enableCudnn
    flags.link{end+1} = sprintf('LDFLAGS=$LDFLAGS -Wl,-rpath -Wl,"%s"', opts.cudnnLibDir) ;
  end
end
if opts.enableGpu
  flags.mexcu = flags.cc ;
  flags.mexcu{end+1} = '-largeArrayDims' ;
  flags.mexcu{end+1} = '-cxx' ;
  flags.mexcu(end+1:end+2) = {'-f' mex_cuda_config(root)} ;
  flags.mexcu{end+1} = ['NVCCFLAGS=' opts.cudaArch '$NVCC_FLAGS'] ;
end

% For the cudaMethod='nvcc'
if opts.enableGpu && strcmp(opts.cudaMethod,'nvcc')
  flags.nvcc = flags.cc ;
  flags.nvcc{end+1} = ['-I"' fullfile(matlabroot, 'extern', 'include') '"'] ;
  flags.nvcc{end+1} = ['-I"' fullfile(matlabroot, 'toolbox','distcomp','gpu','extern','include') '"'] ;
  if opts.debug
    flags.nvcc{end+1} = '-O0' ;
  end
  flags.nvcc{end+1} = '-Xcompiler' ;
  switch arch
    case {'maci64', 'glnxa64'}
      flags.nvcc{end+1} = '-fPIC' ;
    case 'win64'
      flags.nvcc{end+1} = '/MD' ;
      check_clpath(); % check whether cl.exe in path
  end
  flags.nvcc{end+1} = opts.cudaArch;
end

if opts.verbose
  fprintf('%s: * Compiler and linker configurations *\n', mfilename) ;
  fprintf('%s: \tintermediate build products directory: %s\n', mfilename, bld_dir) ;
  fprintf('%s: \tMEX files: %s/\n', mfilename, mex_dir) ;
  fprintf('%s: \tMEX compiler options: %s\n', mfilename, strjoin(flags.mexcc)) ;
  fprintf('%s: \tMEX linker options: %s\n', mfilename, strjoin(flags.link)) ;
end
if opts.verbose & opts.enableGpu
  fprintf('%s: \tMEX compiler options (CUDA): %s\n', mfilename, strjoin(flags.mexcu)) ;
end
if opts.verbose & opts.enableGpu & strcmp(opts.cudaMethod,'nvcc')
  fprintf('%s: \tNVCC compiler options: %s\n', mfilename, strjoin(flags.nvcc)) ;
end
if opts.verbose & opts.enableImreadJpeg
  fprintf('%s: * Reading images *\n', mfilename) ;
  fprintf('%s: \tvl_imreadjpeg enabled\n', mfilename) ;
  fprintf('%s: \timage library: %s\n', mfilename, opts.imageLibrary) ;
  fprintf('%s: \timage library compile flags: %s\n', mfilename, strjoin(opts.imageLibraryCompileFlags)) ;
  fprintf('%s: \timage library link flags: %s\n', mfilename, strjoin(opts.imageLibraryLinkFlags)) ;
end

% --------------------------------------------------------------------
%                                                              Compile
% --------------------------------------------------------------------

% Intermediate object files
srcs = horzcat(lib_src,mex_src) ;
for i = 1:numel(horzcat(lib_src, mex_src))
  [~,~,ext] = fileparts(srcs{i}) ; ext(1) = [] ;
  if strcmp(ext,'cu')
    if strcmp(opts.cudaMethod,'nvcc')
      nvcc_compile(opts, srcs{i}, toobj(bld_dir,srcs{i}), flags.nvcc) ;
    else
      mex_compile(opts, srcs{i}, toobj(bld_dir,srcs{i}), flags.mexcu) ;
    end
  else
    mex_compile(opts, srcs{i}, toobj(bld_dir,srcs{i}), flags.mexcc) ;
  end
end

% Link into MEX files
for i = 1:numel(mex_src)
  [~,base,~] = fileparts(mex_src{i}) ;
  objs = toobj(bld_dir, {mex_src{i}, lib_src{:}}) ;
  mex_link(opts, objs, mex_dir, flags.link) ;
end

% Reset path adding the mex subdirectory just created
vl_setupnn() ;

% --------------------------------------------------------------------
%                                                    Utility functions
% --------------------------------------------------------------------

% --------------------------------------------------------------------
function objs = toobj(bld_dir,srcs)
% --------------------------------------------------------------------
str = fullfile('matlab','src') ;
multiple = iscell(srcs) ;
if ~multiple, srcs = {srcs} ; end
for t = 1:numel(srcs)
  i = strfind(srcs{t},str);
  objs{t} = fullfile(bld_dir, srcs{t}(i+numel(str):end)) ;
end
if ~multiple, objs = objs{1} ; end
objs = strrep(objs,'.cpp',['.' objext]) ;
objs = strrep(objs,'.cu',['.' objext]) ;
objs = strrep(objs,'.c',['.' objext]) ;

% --------------------------------------------------------------------
function objs = mex_compile(opts, src, tgt, mex_opts)
% --------------------------------------------------------------------
mopts = {'-outdir', fileparts(tgt), src, '-c', mex_opts{:}} ;
opts.verbose && fprintf('%s: MEX CC: %s\n', mfilename, strjoin(mopts)) ;
mex(mopts{:}) ;

% --------------------------------------------------------------------
function obj = nvcc_compile(opts, src, tgt, nvcc_opts)
% --------------------------------------------------------------------
nvcc_path = fullfile(opts.cudaRoot, 'bin', 'nvcc');
nvcc_cmd = sprintf('"%s" -c "%s" %s -o "%s"', ...
                   nvcc_path, src, ...
                   strjoin(nvcc_opts), tgt);
opts.verbose && fprintf('%s: NVCC CC: %s\n', mfilename, nvcc_cmd) ;
status = system(nvcc_cmd);
if status, error('Command %s failed.', nvcc_cmd); end;

% --------------------------------------------------------------------
function mex_link(opts, objs, mex_dir, mex_flags)
% --------------------------------------------------------------------
mopts = {'-outdir', mex_dir, mex_flags{:}, objs{:}} ;
opts.verbose && fprintf('%s: MEX LINK: %s\n', mfilename, strjoin(mopts)) ;
mex(mopts{:}) ;

% --------------------------------------------------------------------
function ext = objext()
% --------------------------------------------------------------------
% Get the extension for an 'object' file for the current computer
% architecture
switch computer('arch')
  case 'win64', ext = 'obj';
  case {'maci64', 'glnxa64'}, ext = 'o' ;
  otherwise, error('Unsupported architecture %s.', computer) ;
end

% --------------------------------------------------------------------
function conf_file = mex_cuda_config(root)
% --------------------------------------------------------------------
% Get mex CUDA config file
mver = [1e4 1e2 1] * sscanf(version, '%d.%d.%d') ;
if mver <= 80200, ext = 'sh' ; else ext = 'xml' ; end
arch = computer('arch') ;
switch arch
  case {'win64'}
    config_dir = fullfile(matlabroot, 'toolbox', ...
                          'distcomp', 'gpu', 'extern', ...
                          'src', 'mex', arch) ;
  case {'maci64', 'glnxa64'}
    config_dir = fullfile(root, 'matlab', 'src', 'config') ;
end
conf_file = fullfile(config_dir, ['mex_CUDA_' arch '.' ext]);
fprintf('%s:\tCUDA: MEX config file: ''%s''\n', mfilename, conf_file);

% --------------------------------------------------------------------
function check_clpath()
% --------------------------------------------------------------------
% Checks whether the cl.exe is in the path (needed for the nvcc). If
% not, tries to guess the location out of mex configuration.
status = system('cl.exe -help');
if status == 1
  warning('CL.EXE not found in PATH. Trying to guess out of mex setup.');
  cc = mex.getCompilerConfigurations('c++');
  if isempty(cc)
    error('Mex is not configured. Run "mex -setup".');
  end
  prev_path = getenv('PATH');
  cl_path = fullfile(cc.Location, 'VC','bin','x86_amd64');
  setenv('PATH', [prev_path ';' cl_path]);
  status = system('cl.exe');
  if status == 1
    setenv('PATH', prev_path);
    error('Unable to find cl.exe');
  else
    fprintf('Location of cl.exe (%s) successfully added to your PATH.\n', ...
      cl_path);
  end
end

% -------------------------------------------------------------------------
function paths = which_nvcc(opts)
% -------------------------------------------------------------------------
switch computer('arch')
  case 'win64'
    [~, paths] = system('where nvcc.exe');
    paths = strtrim(paths);
    paths = paths(strfind(paths, '.exe'));
  case {'maci64', 'glnxa64'}
    [~, paths] = system('which nvcc');
    paths = strtrim(paths) ;
end

% -------------------------------------------------------------------------
function cuda_root = search_cuda_devkit(opts)
% -------------------------------------------------------------------------
% This function tries to to locate a working copy of the CUDA Devkit.

opts.verbose && fprintf(['%s:\tCUDA: seraching for the CUDA Devkit' ...
                    ' (use the option ''CudaRoot'' to override):\n'], mfilename);

% Propose a number of candidate paths for NVCC
paths = {getenv('MW_NVCC_PATH')} ;
paths = [paths, which_nvcc(opts)] ;
for v = {'5.5', '6.0', '6.5', '7.0'}
  switch computer('arch')
    case 'glnxa64'
      paths{end+1} = sprintf('/usr/local/cuda-%s/bin/nvcc', char(v)) ;
    case 'maci64'
      paths{end+1} = sprintf('/Developer/NVIDIA/CUDA-%s/bin/nvcc', char(v)) ;
    case 'win64'
      paths{end+1} = sprintf('C:\\Program Files\\NVIDIA GPU Computing Toolkit\\CUDA\\v%s\\bin\\nvcc.exe', char(v)) ;
  end
end
paths{end+1} = sprintf('/usr/local/cuda/bin/nvcc') ;

% Validate each candidate NVCC path
for i=1:numel(paths)
  nvcc(i).path = paths{i} ;
  [nvcc(i).isvalid, nvcc(i).version] = validate_nvcc(paths{i}) ;
end
if opts.verbose
  fprintf('\t| %5s | %5s | %-70s |\n', 'valid', 'ver', 'NVCC path') ;
  for i=1:numel(paths)
    fprintf('\t| %5d | %5d | %-70s |\n', ...
            nvcc(i).isvalid, nvcc(i).version, nvcc(i).path) ;
  end
end

% Pick an entry
index = find([nvcc.isvalid]) ;
if isempty(index)
  error('Could not find a valid NVCC executable\n') ;
end
nvcc = nvcc(index(1)) ;
cuda_root = fileparts(fileparts(nvcc.path)) ;

if opts.verbose
  fprintf('%s:\tCUDA: choosing NVCC compiler ''%s'' (version %d)\n', ...
          mfilename, nvcc.path, nvcc.version) ;
end

% -------------------------------------------------------------------------
function [valid, cuver]  = validate_nvcc(nvccPath)
% -------------------------------------------------------------------------
[status, output] = system(sprintf('"%s" --version', nvccPath)) ;
valid = (status == 0) ;
if ~valid
  cuver = 0 ;
  return ;
end
match = regexp(output, 'V(\d+\.\d+\.\d+)', 'match') ;
if isempty(match), valid = false ; return ; end
cuver = [1e4 1e2 1] * sscanf(match{1}, 'V%d.%d.%d') ;

% --------------------------------------------------------------------
function cuver = activate_nvcc(nvccPath)
% --------------------------------------------------------------------

% Validate the NVCC compiler installation
[valid, cuver] = validate_nvcc(nvccPath) ;
if ~valid
  error('The NVCC compiler ''%s'' does not appear to be valid.', nvccPath) ;
end

% Make sure that NVCC is visible by MEX by setting the MW_NVCC_PATH
% environment variable to the NVCC compiler path
if ~strcmp(getenv('MW_NVCC_PATH'), nvccPath)
  warning('Setting the ''MW_NVCC_PATH'' environment variable to ''%s''', nvccPath) ;
  setenv('MW_NVCC_PATH', nvccPath) ;
end

% In some operating systems and MATLAB versions, NVCC must also be
% available in the command line search path. Make sure that this is%
% the case.
[valid_, cuver_] = validate_nvcc('nvcc') ;
if ~valid || cuver_ ~= cuver
  warning('NVCC not found in the command line path or the one found does not matches ''%s''.', nvccPath);
  nvccDir = fileparts(nvccPath) ;
  prevPath = getenv('PATH') ;
  switch computer
    case 'PCWIN64', separator = ';' ;
    case {'GLNXA64', 'MACI64'}, separator = ':' ;
  end
  setenv('PATH', [nvccDir separator prevPath]) ;
  [valid_, cuver_] = validate_nvcc('nvcc') ;
  if ~valid_ || cuver_ ~= cuver
    setenv('PATH', prevPath) ;
    error('Unable to set the command line path to point to ''%s'' correctly.', nvccPath) ;
  else
    fprintf('Location of NVCC (%s) added to your command search PATH.\n', nvccDir) ;
  end
end

% --------------------------------------------------------------------
function cudaArch = get_cuda_arch(opts)
% --------------------------------------------------------------------
opts.verbose && fprintf('%s:\tCUDA: determining GPU compute capability (use the ''CudaArch'' option to override)\n', mfilename);
try
  gpu_device = gpuDevice();
  arch_code = strrep(gpu_device.ComputeCapability, '.', '');
  cudaArch = ...
      sprintf('-gencode=arch=compute_%s,code=\\\"sm_%s,compute_%s\\\" ', ...
              arch_code, arch_code, arch_code) ;
catch
  opts.verbose && fprintf(['%s:\tCUDA: cannot determine the capabilities of the installed GPU;' ...
                      'falling back to default\n'], mfilename);
  cudaArch = opts.defCudaArch;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function Offspring = Operator(Problem,Particle,Pbest,Gbest)
% Particle swarm optimization in MMOPSO

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Parameter setting
    ParticleDec = Particle.decs;
    PbestDec    = Pbest.decs;
    GbestDec    = Gbest.decs;
    [N,D]       = size(ParticleDec);
    ParticleVel = Particle.adds(zeros(N,D));

    %% Particle swarm optimization
    W  = repmat(unifrnd(0.1,0.5,N,1),1,D);
    r1 = repmat(rand(N,1),1,D);
    r2 = repmat(rand(N,1),1,D);
    C1 = repmat(unifrnd(1.5,2,N,1),1,D);
    C2 = repmat(unifrnd(1.5,2,N,1),1,D);
    OffVel        = W.*ParticleVel;
    temp          = repmat(rand(N,1)<0.7,1,D);
    OffVel(temp)  = OffVel(temp) + C1(temp).*r1(temp).*(PbestDec(temp)-ParticleDec(temp));
    OffVel(~temp) = OffVel(~temp) + C2(~temp).*r2(~temp).*(GbestDec(~temp)-ParticleDec(~temp));
    OffDec        = ParticleDec + OffVel;
    Offspring     = Problem.Evaluation(OffDec,OffVel);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function view = cmapSetStandard(view, mode, phInMiddle)
%
%   view = cmapSetStandard(view,mode,phInMiddle)
%
% Author: JL, BW
% Purpose:
%     Put in place a standard color map for the different modes.

if ieNotDefined('mode'), mode = 'ph'; end

switch mode
    case 'ph'
    case 'co'
    case 'amp'
    otherwise
end

% ??? view = refreshView(view);

return;

%---------------------------
function cmap = cmapCenter(mp,centerPhase)

[mp, numGrays,numColors] = getColorMap(FLAT{2},'ph',1);

step = numColors/(2*pi);
horPhase = 2;

% We want the part of the default map from pi to 2pi to be centered around
% the phase that represents the horizontal midline.  That means we want to
% shift the map so that 3pi/2 in the default color map shifts down to the
% horizontal position.

sz = round(( (3*pi/2) - horPhase) *step);
cmap = circshift(mp,sz)

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [m, C, W] = gaussian_CPD_params_given_dps(CPD, domain, evidence)
% GAUSSIAN_CPD_PARAMS_GIVEN_EV_ON_DPS Extract parameters given evidence on all discrete parents
% function [m, C, W] = gaussian_CPD_params_given_ev_on_dps(CPD, domain, evidence)

ps = domain(1:end-1);
dps = ps(CPD.dps);
if isempty(dps)
  m = CPD.mean;
  C = CPD.cov;
  W = CPD.weights;
else
  odom = domain(~isemptycell(evidence(domain)));
  dops = myintersect(dps, odom);
  dpvals = cat(1, evidence{dops});
  if length(dops) == length(dps)
    dpsizes = CPD.sizes(CPD.dps);
    dpval = subv2ind(dpsizes, dpvals(:)');
    m = CPD.mean(:, dpval);
    C = CPD.cov(:, :, dpval);
    W = CPD.weights(:, :, dpval);
  else
    map = find_equiv_posns(dops, dps);
    index = mk_multi_index(length(dps), map, dpvals);
    m = CPD.mean(:, index{:});
    C = CPD.cov(:, :, index{:});
    W = CPD.weights(:, :, index{:});
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Configurations for the algorithsm
% Input:
%       configName:  Name of the config, corresponding to the desired
%       dataset
function config = readConfig(configName)

    if (nargin==0); configName = 'synthetic'; end

    if (strcmp(configName, 'synthetic'))
        config = readConfig_synthetic();

     % To be added
%     elseif (strcmp(configName, 'redwood'))
%         config = readConfigRedWood();
    end

end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function fixed = p11_fixed_points ( m, fixed_num )

%*****************************************************************************80
%
%% P11_FIXED_POINTS returns the fixed points in problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer FIXED_NUM, the number of fixed points.
%
%    Output, real FIXED(M,FIXED_NUM), the fixed points.
%
  fixed = [    ...
    0.0,  0.0; ...
    1.0,  0.0; ...
    1.0,  0.5; ...
    0.5,  0.5; ...
    0.5,  1.0; ...
    0.0,  1.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [xDistSort,yDistSort]=sortPointMinDist(x,y)
%% function [xDistSort,yDistSort]=sortPointMinDist(x,y)
% Sorts arbitrary points into a Clock Wise contour.
%
%% Syntax
% [xDistSort,yDistSort]=sortPointMinDist(x,y)
%
%% Description
% This functions goal is to order points efficintly to allow plotting a contoure with
% them. While this problem can be hard to solve, and even sometimes considered unsolvable,
% it is solved here basing on the assumption that "nearest point" is the nect point. This
% results in a solution of O(N^2) complexity (calculating all distamnces between points),
% but im most cases it will be of significantly lower complexity.
%
%% Input arguments:
% x- x coordinates of the points
%
% y- x coordinates of the points
%
%% Output arguments
% xDistSort- sorted points x coordinates.
%
% yDistSort- sorted points y coordinates.
%
%% Issues & Comments (None)
% The resulting shape will have points connected,accordung to a predfined metric, so it
% will not always result in a shape user wished for. Add additonal points to fix this
% issue, when needed.
% The funciton is pretty damandig computationally (~20 slower than sortPoint2ContourCW),
% so should be used not too frequently.
%
%% Example
% N=11;
% x=10*rand(1,N);
% y=10*rand(1,N);
% [xCW,yCW]=sortPointMinDist(x,y);
% figure;
% plot(x,y,'.-b');
% hold on;
% plot(xCW,yCW,'.-r','LineWidth',2);
% hold off;
% axis equal;
% title('Sorting arbitrary Points to form a Contour','FontSize',14);
% legend('Unsorted points contour', 'Sorted points contour');
%
%% See also
% sortPoint2ContourCW;  % Custom function
% mask2poly;            % Custom function
% poly2mask;            % Matlab function
%
%% Revision history
% First version: Nikolay S. 2011-07-24.
% Last update:   Nikolay S. 2011-07-25.
%
% *List of Changes:*
%

%% convert to column vectors and sort
nPoints=length(x);
iDistSort=zeros(1,nPoints);
distMat=Inf(nPoints,nPoints); % col index- distance from, row index- distance to

for iDistRow=1:nPoints
   ind2=(iDistRow+1):nPoints;
   distMat(iDistRow, ind2)=sqrt((x(ind2)-x(iDistRow)).^2+(y(ind2)-y(iDistRow)).^2);
   % distance between a->b eqauls distance between b->a
   distMat(ind2, iDistRow)=transpose(distMat(iDistRow, ind2));
end

% find pair of closest points- first and second points
[~,pointLinInd] = min(distMat(:));
[iDistSort(1),nextPoint] = ind2sub(size(distMat), pointLinInd);
for iPoint=1:nPoints-1
   currPoint=nextPoint;
   [minDist,nextPoint] = min(distMat(currPoint, :)); % find next point- closest to current
   distMat(currPoint,:)=Inf; % delete distances from current point
   distMat(:,currPoint)=Inf; % delete distances to   current point
   iDistSort(iPoint+1)=nextPoint; % store sorted points indexes

   if isinf(minDist) % this will b true in case of an error
      break;
   end
end
xDistSort=x(iDistSort);
yDistSort=y(iDistSort);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef ViewGridTypeEnum < matlab.mixin.SetGet
    %ViewGridTypeEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Major('Major Grid Lines','on')
        Minor('Major and Minor Grid Lines','minor')
        Off('No Grid Lines','off')
    end

    properties
        name(1,:) char = '';
        gridStr(1,:) char = '';
    end

    methods
        function obj = ViewGridTypeEnum(name, gridStr)
            obj.name = name;
            obj.gridStr = gridStr;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('ViewGridTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('ViewGridTypeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('ViewGridTypeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Abar = computeAbar(A, ny)
    Abar_tmp = cell(ny,1);
    tmp = A;

    for i=1:ny
        Abar_tmp{i,1} = tmp;
        tmp = tmp*A;
    end

    Abar = cell2mat(Abar_tmp);

end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
clear;
addpath('../../CVPR16_VDSR/');
addpath('./matlabPyrTools-master/');
%% test
% load('DRRN_22C128_dataset_1_x3');
load('VDSR_20C64_dataset_1_x3');
path = ['../../SR Test dataset/Set5/'];
d = dir([path '*.bmp']);
up_scale = 3;
sum = 0 ;
for i=1:length(d)
    im = imread([path d(i).name]);
    %% work on illuminance only
    if size(im,3)>1
        im = rgb2ycbcr(im);
    end
    im = im2double(im(:,:,1));
    im_gnd = modcrop(im, up_scale);
    [hei,wid] = size(im_gnd);
    im_b = imresize(imresize(im_gnd,1/up_scale,'bicubic'),up_scale,'bicubic');
    %% remove border
    im_b1 = shave(uint8(single(im_b) * 255), [up_scale, up_scale]);
    im_gnd1 = shave(uint8(single(im_gnd) * 255), [up_scale, up_scale]);
    im_h1 = im_h_set{i};
    sum = sum+ifcvec(double(im_gnd1),double(im_h1));
end
ifc = sum/length(d)

author: tyshiwo


text:
function test_anls(varargin)
%
% demonstration file for NMFLibrary.
%
% This file illustrates how to use this library.
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on June 23, 2022

    if nargin < 1
        clc;
        clear;
        close all;
        rng('default')

        m = 500;
        n = 100;
        V = rand(m,n);
        rank = 20;
        options = [];
        options.verbose = 2;
        options.max_epoch = 100;
        health_check_mode = false;
    else
        V = varargin{1};
        rank = varargin{2};
        options = varargin{3};
        health_check_mode = true;
    end

    % initialize factor matrices
    [x_init, ~] = generate_init_factors(V, rank, []);
    inner_max_epoch = 10;
    rescale = true;

    %% Scale initialization so that argmin_a ||a * WH - X||_F = 1
    W = x_init.W;
    H = x_init.H;
    if rescale
        [W, H] = normalize_WH(V, W, H, rank, 'type1');
    end

    x_init.W = W;
    x_init.H = H;

    %% ANLS
    %algorithm = 'anls_bpp';
    %algorithm = 'anls_asgroup';
    algorithm = 'anls_asgivens';
    % standard
    options.momentum_h = 0;
    options.momentum_w = 0;
    options.inner_max_epoch = inner_max_epoch;
    options.alg = algorithm;

    [w_nmf_anls, infos_nmf_anls] = anls_nmf(V, rank, options);

    % acc only H
    options.momentum_h = 3;
    options.momentum_w = 0;
    options.inner_max_epoch = inner_max_epoch;
    options.alg = algorithm;
    [w_nmf_anls_momemtum_h, infos_nmf_anls_momemtum_h] = anls_nmf(V, rank, options);

    options.momentum_h = 3;
    options.momentum_w = 2;
    options.inner_max_epoch = inner_max_epoch;
    options.alg = algorithm;
    [w_nmf_anls_momemtum_hw, infos_nmf_anls_momemtum_hw] = anls_nmf(V, rank, options);

    options.momentum_h = 3;
    options.momentum_w = 2;
    options.inner_max_epoch = inner_max_epoch;
    options.scaling = false;
    options.alg = algorithm;
    [w_nmf_anls_momemtum_hw_inner, infos_nmf_anls_momemtum_hw_inner] = anls_nmf(V, rank, options);

    if ~health_check_mode

        %% plot
        display_graph('iter','cost', {'ANLS-Standard', 'ANLS-MOMENTUM (H)', 'ANLS-MOMENTUM (W,H)', 'ANLS-MOMENTUM-nonscale (W,H)'}, ...
                                {w_nmf_anls, w_nmf_anls_momemtum_h, w_nmf_anls_momemtum_hw, w_nmf_anls_momemtum_hw_inner}, ...
                                {infos_nmf_anls, infos_nmf_anls_momemtum_h, infos_nmf_anls_momemtum_hw, infos_nmf_anls_momemtum_hw_inner});
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function obj=addFromSuite(obj, s)
% Add MOxUnitTestSuite tests to the suite
%
% obj=addFromSuite(obj, s)
%
% Inputs:
%   obj             MoxUnitTestSuite instance
%   s               MOxUnitTestSuite instance with tests to be added to
%                   the suite.
%
% Output:
%   obj             MoxUnitTestSuite instance with all tests in s
%                   added to the suite
%
% See also: initTestSuite
%
% NNO 2015

    n_test_nodes=countTestNodes(s);

    for k=1:n_test_nodes
        test_node=getTestNode(s, k);
        obj=addTest(obj, test_node);
    end


text:
k=0;
while k<4
	k=k+1;
	B(k)=k.^2
end

author: Sable

repo: mcbench-benchmarks


text:
function match_files = getRegexFiles(path, regex)
% files = dir([path '/*.mat']);
files = dir(path);
names = {files.name};
match_indices = ~cellfun(@isempty, regexp(names, regex, 'match', 'once'));
match_files = names(match_indices);
end

author: ramintoosi

repo: ROSS

save_path: github-repos/MATLAB/ramintoosi-ROSS


text:
function [P, Transf] = cs_convert(sMri, src, dest, P, isNanAllowed)
% CS_CONVERT: Convert 3D points between coordinates systems.
%
% USAGE:       P = cs_convert(sMri, src, dest, P, isNanAllowed=0)
%         Transf = cs_convert(sMri, src, dest)
%
% INPUT:
%     - sMri         : Brainstorm MRI structure
%     - src          : Current coordinate system {'voxel','mri','scs','mni','acpc','world','captrak'}
%     - dest         : Target coordinate system {'voxel','mri','scs','mni','acpc','world','captrak'}
%     - P            : a Nx3 matrix of point coordinates to convert
%     - isNanAllowed : If 0, throw an error and if the conversion of some points P leads to NaN,
%                      and if there is no linear MNI registration to be used instead
%                      (case of non-linear MNI registration with points outside of the defined area)
%
% DESCRIPTION:   https://neuroimage.usc.edu/brainstorm/CoordinateSystems
%     - voxel : X=left>right,  Y=posterior>anterior,   Z=bottom>top
%               Coordinate of the center of the first voxel at the bottom-left-posterior of the MRI volume: (1,1,1)
%     - mri   : Same as 'voxel' but in millimeters instead of voxels:  mriXYZ = voxelXYZ * Voxsize
%     - scs   : Based on: Nasion, left pre-auricular point (LPA), and right pre-auricular point (RPA).
%               Origin: Midway on the line joining LPA and RPA
%               Axis X: From the origin towards the Nasion (exactly through)
%               Axis Y: From the origin towards LPA in the plane defined by (NAS,LPA,RPA), and orthogonal to X axis
%               Axiz Z: From the origin towards the top of the head
%     - mni   : MNI coordinates based on SPM affine or non-linear registration
%     - acpc  : Based on: Anterior commissure (AC), Posterior commissure (PC) and an interhemisperic point (IH)
%               Origin: AC
%               Axis X: From the origin towards the right
%               Axis Y: Negative y-axis is passing from AC through PC
%               Axis Z: Passing through a mid-hemispheric point in the superior direction
%     - world : Transformation available in the initial file loaded as the default MRI (vox2ras/qform/world transformation)
%     - captrak: RAS orientation and the origin approximately between LPA and RPA
%               Axis X: From LPA through RPA exactly
%               Axis Y: Orthogonal to the X-axis through the nasion (NAS)
%               Axis Z: Orthogonal to the XY-plane through the vertex of the head

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

% Throw errors by default if NaN
if (nargin < 5) || isempty(isNanAllowed)
    isNanAllowed = 0;
end
% Check matrices orientation
if (nargin < 4) || isempty(P)
    P = [];
elseif (size(P,2) ~= 3) && (size(P,1) == 3)
    P = P';
elseif (size(P,2) ~= 3)
    error('P must have 3 columns (X,Y,Z).');
end
% If the coordinate system didn't change
if strcmpi(src, dest)
    return;
end
% Keep track of points that cannot be transformed
iMissing = [];
isApplied = 0;
Porig = P;
% Transform to homogeneous coordinates
if ~isempty(P)
    P = [P'; ones(1,size(P,1))];
end

% ===== GET MRI=>WORLD TRANSFORMATION =====
if strcmpi(src, 'world') || strcmpi(dest, 'world') || (strcmpi(src, 'mni') && isfield(sMri,'NCS') && isfield(sMri.NCS,'y') && ~isempty(sMri.NCS.y))
    % Get the vox2ras transformation
    if isempty(sMri) || isempty(sMri.InitTransf)
        P = [];
        return;
    end
    iTransf = find(strcmpi(sMri.InitTransf(:,1), 'vox2ras'));
    if isempty(iTransf)
        P = [];
        return;
    end
    vox2ras = sMri.InitTransf{iTransf,2};
    % 2nd operation: Change reference from (0,0,0) to (1,1,1)
    vox2ras = vox2ras * [1 0 0 -1; 0 1 0 -1; 0 0 1 -1; 0 0 0 1];
    % 1st operation: Convert from MRI(mm) to voxels
    vox2ras = vox2ras * diag(1 ./ [sMri.Voxsize, 1]);
    % Convert to meters: transformation MRI=>WORLD (in meters)
    mri2world = vox2ras .* [ones(3,3), 1e-3.*ones(3,1); ones(1,4)];
    % Compute inverse transformation WORLD=>MRI (in meters)
    world2mri = inv(mri2world);
end

% ===== COMPUTE ACPC TRANSFORMATION =====
if strcmpi(src, 'acpc') || strcmpi(dest, 'acpc')
    tACPC = cs_compute(sMri, 'acpc');
    if isempty(tACPC) || isempty(tACPC.R)
        P = [];
        return;
    end
    scs2acpc = [tACPC.R, tACPC.T; 0 0 0 1];
end

% ===== COMPUTE CAPTRAK TRANSFORMATION =====
if strcmpi(src, 'captrak') || strcmpi(dest, 'captrak')
    tCapTrak = cs_compute(sMri, 'captrak');
    if isempty(tCapTrak) || isempty(tCapTrak.R)
        P = [];
        return;
    end
    scs2captrak = [tCapTrak.R, tCapTrak.T; 0 0 0 1];
end

% ===== CONVERT SRC => MRI =====
% Evaluate the transformation to apply
switch lower(src)
    case 'voxel'
        RT1 = diag([sMri.Voxsize(:) ./ 1000; 1]);
    case 'mri'
        RT1 = eye(4);
    case 'scs'
        if ~isfield(sMri,'SCS') || ~isfield(sMri.SCS,'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS,'T') || isempty(sMri.SCS.T)
            P = [];
            return;
        end
        RT1 = inv([sMri.SCS.R, sMri.SCS.T./1000; 0 0 0 1]);
    case 'mni'
        % Transformation of each point by indirection in the deformation field y
        if ~isempty(P) && isfield(sMri,'NCS') && isfield(sMri.NCS,'y') && ~isempty(sMri.NCS.y)
            % Convert MNI => voxel space of the registration matrix
            P_reg = inv(sMri.NCS.y_vox2ras) * (P .* [1000;1000;1000;1]);
            % Convert from 0-based to 1-based??
            % => This solution was obtained empirically by minimizing:
            %    sqrt(sum((cs_convert(sMri, 'mri', 'mni', cs_convert(sMri, 'mni', 'mri', P)) - P).^2)).*1000 => around 0.003 with this adjustment
            P_reg = P_reg + [1;1;1;0];
            % Convert Voxel => World
            P_world = [...
                interp3(sMri.NCS.y(:,:,:,1), P_reg(2,:), P_reg(1,:), P_reg(3,:), 'linear', NaN); ...
                interp3(sMri.NCS.y(:,:,:,2), P_reg(2,:), P_reg(1,:), P_reg(3,:), 'linear', NaN); ...
                interp3(sMri.NCS.y(:,:,:,3), P_reg(2,:), P_reg(1,:), P_reg(3,:), 'linear', NaN)] ./ 1000;
            % Check if some points could not be converted
            if ~isNanAllowed
                iMissing = find(any(isnan(P_world),1));
            end
            % Convert World => MRI
            P = world2mri * [double(P_world); ones(1, size(P_world,2))];
            RT1 = eye(4);
        elseif isfield(sMri,'NCS') && isfield(sMri.NCS,'R') && ~isempty(sMri.NCS.R) && isfield(sMri.NCS,'T') && ~isempty(sMri.NCS.T)
            RT1 = inv([sMri.NCS.R, sMri.NCS.T./1000; 0 0 0 1]);
        else
            P = [];
            return;
        end
    case 'acpc'
        % ACPC => SCS => MRI
        RT1 = inv(scs2acpc);
    case 'captrak'
        % CapTrak => SCS => MRI
        RT1 = inv(scs2captrak);
    case 'world'
        RT1 = world2mri;
    otherwise
        error(['Invalid coordinate system: ' src]);
end

% ===== CONVERT MRI => DEST =====
% Evaluate the transformation to apply
switch lower(dest)
    case 'voxel'
        RT2 = diag([1000 ./ sMri.Voxsize(:); 1]);
    case 'mri'
        RT2 = eye(4);
    case 'scs'
        if ~isfield(sMri,'SCS') || ~isfield(sMri.SCS,'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS,'T') || isempty(sMri.SCS.T)
            P = [];
            return;
        end
        RT2 = [sMri.SCS.R, sMri.SCS.T./1000; 0 0 0 1];
    case 'mni'
        % Using non-linear MNI normalization: Transformation of each point by indirection in the deformation field iy
        if ~isempty(P) && isfield(sMri,'NCS') && isfield(sMri.NCS,'iy') && ~isempty(sMri.NCS.iy)
            % Convert: src => MRI => voxel
            P_vox = diag([1000 ./ sMri.Voxsize(:); 1]) * RT1 * P;
            % Get values from the iy volumes
            Pmni = [interp3(sMri.NCS.iy(:,:,:,1), P_vox(2,:), P_vox(1,:), P_vox(3,:), 'linear', NaN); ...
                 interp3(sMri.NCS.iy(:,:,:,2), P_vox(2,:), P_vox(1,:), P_vox(3,:), 'linear', NaN); ...
                 interp3(sMri.NCS.iy(:,:,:,3), P_vox(2,:), P_vox(1,:), P_vox(3,:), 'linear', NaN)] ./ 1000;
            % Check if some points could not be converted
            if ~isNanAllowed
                iMissing = find(any(isnan(Pmni),1));
            end
            % Transpose the matrix back
            P = double(Pmni');
            Transf = [];
            isApplied = 1;
        elseif isfield(sMri,'NCS') && isfield(sMri.NCS,'R') && ~isempty(sMri.NCS.R) && isfield(sMri.NCS,'T') && ~isempty(sMri.NCS.T)
            RT2 = [sMri.NCS.R, sMri.NCS.T./1000; 0 0 0 1];
        else
            P = [];
            return;
        end
    case 'acpc'
        % MRI => SCS => ACPC
        RT2 = scs2acpc;
    case 'captrak'
        % MRI => SCS => CapTrak
        RT2 = scs2captrak;
    case 'world'
        RT2 = mri2world;
    otherwise
        error(['Invalid coordinate system: ' dest]);
end

% If the final transformation is not already applied
if ~isApplied
    % Compute the final transformation matrix
    Transf = RT2 * RT1;
    % Apply the transformation matrix to the points
    if ~isempty(P)
        % Apply rotation-translation
        P = Transf * P;
        % Remove the last coordinate and transpose the matrix back
        P = P(1:3,:)';
    else
        P = Transf;
    end
end

% If there are NaN points from non-linear MNI: Try to fix them with linear MNI
if ~isempty(iMissing)
    % If there is a linear transformation available as well, use it for missing points (same code as below)
    if isfield(sMri,'NCS') && isfield(sMri.NCS,'R') && ~isempty(sMri.NCS.R) && isfield(sMri.NCS,'T') && ~isempty(sMri.NCS.T)
        % Remove the non-linear transformation, so that it forces using the linear one
        sMri.NCS.iy = [];
        sMri.NCS.y = [];
        P(iMissing,:) = cs_convert(sMri, src, dest, Porig(iMissing,:), isNanAllowed);
    % If NaN not allowed and no linear transformation: check for missing values
    elseif ~isNanAllowed
        error(['Some points are outside the definition of the non-linear MNI registration.' 10 'Please compute a linear MNI normalization to convert these points.']);
    end
end

author: brainstorm-tools


text:
%% kabschRotationMatrix
% Below is a demonstration of the features of the |kabschRotationMatrix| function

%%
clear; close all; clc;

%% Syntax
% |[Q]=kabschRotationMatrix(V1,V2);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function colormask = wbmask(m,n,wbmults,align)
% COLORMASK = wbmask(M,N,WBMULTS,ALIGN)
%
% Makes a white-balance multiplicative mask for an image of size m-by-n
% with RGB while balance multipliers WBMULTS = [R_scale G_scale B_scale].
% ALIGN is string indicating Bayer arrangement: 'rggb','gbrg','grbg','bggr'
colormask = wbmults(2)*ones(m,n); %Initialize to all green values
switch align
case 'rggb'
colormask(1:2:end,1:2:end) = wbmults(1); %r
colormask(2:2:end,2:2:end) = wbmults(3); %b
case 'bggr'
colormask(2:2:end,2:2:end) = wbmults(1); %r
colormask(1:2:end,1:2:end) = wbmults(3); %b
case 'grbg'
colormask(1:2:end,2:2:end) = wbmults(1); %r
colormask(2:2:end,1:2:end) = wbmults(3); %b
case 'gbrg'
colormask(2:2:end,1:2:end) = wbmults(1); %r
colormask(1:2:end,2:2:end) = wbmults(3); %b
end
end

author: AbdoKamel

repo: sidd-ground-truth-image-estimation


text:
function ori = discreteSample(odf,npoints,varargin)
% draw a random sample
%
%

% preallocate orientations
q = quaternion.id(npoints,1);

% which component
if numel(odf.weights) == 1
  icmp = ones(size(q));
else
  icmp = discretesample(odf.weights,npoints);
end

% compute discrete sample for each component seperately
for ic = 1:length(odf.components)
  q(icmp == ic) = discreteSample(odf.components{ic},sum(icmp==ic),varargin{:});
end

% take random symmetrically equivalent samples
qcs = reshape(quaternion(odf.CS.properGroup),[],1);
qss = reshape(quaternion(odf.SS.properGroup),[],1);
ics = discretesample(length(qcs),npoints,1);
iss = discretesample(length(qss),npoints,1);

ori = orientation(qss(iss(:)) .* q .* qcs(ics(:)),odf.CS,odf.SS);

% the antipodal case
if  odf.antipodal
  ori.antipodal =true;
  doInv = 1==discretesample(2,npoints,1);
  ori(doInv)=inv(ori(doInv));
end

author: mtex-toolbox

repo: mtex

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


text:
function varargout = cstructdecode(buf, varargin);

% CSTRUCTDECODE decodes a structure from a uint8 buffer
%
% See READ_NEURALYNX_NEV for an example

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

if ~isa(buf, 'uint8')
  error('incorrect type of input data, should be uint8');
end

nbytes = numel(buf);
nfield = length(varargin);

wordsize = zeros(1,nfield);
for i=1:nfield
  switch varargin{i}
  case 'uint8'
    wordsize(i) = 1;
  case 'int8'
    wordsize(i) = 1;
  case 'uint16'
    wordsize(i) = 2;
  case 'int16'
    wordsize(i) = 2;
  case 'uint32'
    wordsize(i) = 4;
  case 'int32'
    wordsize(i) = 4;
  case 'uint64'
    wordsize(i) = 8;
  case 'int64'
    wordsize(i) = 8;
  case {'float32' 'single'}
    varargin{i} = 'single';
    wordsize(i) = 4;
  case {'float64' 'double'}
    varargin{i} = 'double';
    wordsize(i) = 8;
  otherwise
    if strncmp(varargin{i}, 'char', 4)
      if length(varargin{i})>4
        % assume a string like 'char128' which means 128 characters
        wordsize(i) = str2num(varargin{i}(5:end));
        varargin{i} = 'char';
      else
        wordsize(i) = 1;
      end
    else
      error('incorrect type specification');
    end
  end
end

pklen = sum(wordsize);
pknum = nbytes/sum(wordsize);

buf = reshape(buf, pklen, pknum);

for i=1:nfield
  rowbeg = sum(wordsize(1:(i-1)))+1;
  rowend = sum(wordsize(1:(i-0)))+0;
  sel = buf(rowbeg:rowend,:);
  if strcmp(varargin{i}, 'char')
    varargout{i} = char(sel)';
  else
    varargout{i} = typecast(sel(:), varargin{i});
  end
end


text:
% dipfit_erpeeg - fit multiple component dipoles using DIPFIT
%
% Usage:
%         >> [ dipole model EEG] = dipfit_erpeeg(data, chanlocs, 'key', 'val', ...);
%
% Inputs:
%  data      - input data [channel x point]. One dipole per point is
%              returned.
%  chanlocs  - channel location structure (returned by readlocs()).
%
% Optional inputs:
%  'settings'  - [cell array] dipfit settings (arguments to the
%                pop_dipfit_settings() function). Default is none.
%  'dipoles'   - [1|2] use either 1 dipole or 2 dipoles contrain in
%                symetry. Default is 1.
%  'dipplot'   - ['on'|'off'] plot dipoles. Default is 'off'.
%  'plotopt'   - [cell array] dipplot() 'key', 'val' options. Default is
%                'normlen', 'on', 'image', 'fullmri'
%
% Outputs:
%  dipole      - dipole structure ('posxyz' field is the position; 'momxyz'
%                field is the moment and 'rv' the residual variance)
%  model       - structure containing model information ('vol.r' field is
%                radius, 'vol.c' conductances, 'vol.o' the 3-D origin and
%                'chansel', the selected channels).
%  EEG         - faked EEG structure containing erp activation at the place
%                of ICA components but allowing to plot ERP dipoles.
%
% Note: residual variance is set to NaN if Dipfit does not converge
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, Nov. 2003

% Copyright (C) 10/2003 Arnaud Delorme, SCCN/INC/UCSD, arno@sccn.ucsd.edu
%
% 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 [dipoles, model, EEG] = dipfit_erpeeg(DATA, chanlocs, varargin);

    if nargin < 1
        help dipfit_erpeeg;
        return;
    end;

    ncomps = size(DATA,2);
    if size(DATA,1) ~= length(chanlocs)
        error('# of row in ''DATA'' must equal # of channels in ''chanlocs''');
    end;

    % faking an EEG dataset
    % ---------------------
    EEG          = eeg_emptyset;
    EEG.data     = rand(size(DATA,1), 1000);
    EEG.nbchan   = size(DATA,1);
    EEG.pnts     = 1000;
    EEG.trials   = 1;
    EEG.chanlocs = chanlocs;
    EEG.icawinv    = [ DATA DATA ];
    EEG.icaweights = zeros(size([ DATA DATA ]))';
    EEG.icasphere  = zeros(size(DATA,1), size(DATA,1));
    %EEG            = eeg_checkset(EEG);
    EEG.icaact     = EEG.icaweights*EEG.icasphere*EEG.data(:,:);
    EEG.icaact     = reshape( EEG.icaact, size(EEG.icaact,1), size(EEG.data,2), size(EEG.data,3));

    % uses mutlifit to fit dipoles
    % ----------------------------
    EEG            = pop_multifit(EEG, [1:ncomps], varargin{:});

    % process outputs
    % ---------------
    dipoles = EEG.dipfit.model;
    if isfield(dipoles, 'active')
        dipoles = rmfield(dipoles, 'active');
    end;
    if isfield(dipoles, 'select')
        dipoles = rmfield(dipoles, 'select');
    end;
    model   = EEG.dipfit;
    if isfield(model, 'model')
        model   = rmfield(model, 'model');
    end;
    return;

author: PatternRecognition

repo: OpenBMI


text:
function prob_test079 ( )

%*****************************************************************************80
%
%% TEST079 tests GENLOGISTIC_CDF, GENLOGISTIC_CDF_INV, GENLOGISTIC_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST079\n' );
  fprintf ( 1, '  For the Generalized Logistic PDF:\n' );
  fprintf ( 1, '  GENLOGISTIC_PDF evaluates the PDF.\n' );
  fprintf ( 1, '  GENLOGISTIC_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  GENLOGISTIC_CDF_INV inverts the CDF.\n' );

  a = 1.0;
  b = 2.0;
  c = 3.0;

  check = genlogistic_check ( a, b, c );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST079 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF parameter C =             %14f\n', c );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = genlogistic_sample ( a, b, c, seed );

    pdf = genlogistic_pdf ( x, a, b, c );

    cdf = genlogistic_cdf ( x, a, b, c );

    x2 = genlogistic_cdf_inv ( cdf, a, b, c );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function test_libphaseret_overlaynthframe
f = greasy;
a = 128;
M = 1024;
M2 = floor(M/2) + 1;
gl = 1024;
L = dgtlength(numel(f),a,M);
g = firwin('hann',gl);
gamma = 0.25645*gl^2;
gd = long2fir(gabdual(g,a,M),gl);
gshift = fftshift(g);
N = 10;
idx = 9;
g2 = repmat(g.*gd,1,N);
g1 = repmat(fftshift(g.*gd),1,N);
cout = zeros(gl,1);
coutPtr = libpointer('doublePtr',cout);

calllib('libphaseret','phaseret_overlaynthframe_d',g1,gl,N,a,idx,coutPtr);

[~,out2] = comp_overlayframes(g2,a,gl,idx);
[~,out3orig] = overlayframes(g2,a,gl,idx);
out1orig = M*coutPtr.Value;
out3orig = out3orig*M;

out1 = out1orig;
out3 = out3orig;
out1(out1orig==0) = 1;
out1(out1<1e-6) = 1e-6;
out3(out3orig==0) = 1;
out3(out3<1e-6) = 1e-6;

figure(1); plot([gshift./out1,gshift./out3])

function [partrec,frame] = overlayframes(cframes,a,M,n)

N = size(cframes,2);
bufLen = N*a - (a-1) + M-1;
partrec = zeros(bufLen,1);

startidx = ceil(M/2)-1;
idxrange = startidx + [0:floor(M/2),-ceil(M/2)+1:-1];
for ii=0:N-1
    idx = ii*a + idxrange + 1;
    partrec(idx) = partrec(idx) + cframes(:,ii+1);
end

frame = partrec(a*n+1:a*n + M);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function this = t_filter(this, cutoffSeconds)
% high-pass filters temporally (4th dimension of the image) as SPM
%
%   MrImage = t_filter(MrImage)
%
% This is a method of class MrImage.
%
% IN
%   cutoffSeconds   slower drifts than this will be filtered out
%
% OUT
%
% EXAMPLE
%   t_filter
%
%   See also MrImage spm_filter

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-02
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% convert to 2D
nVoxel = this.geometry.nVoxels;
Y = reshape(this.data, [], nVoxel(4))'; % Y = [nVolumes, nVoxel]
nVoxel3D = prod(nVoxel(1:3));

% create K for spm_filter and do it
K.RT = this.geometry.TR_s;
K.HParam = cutoffSeconds;
K.row = 1:nVoxel(4);

% spm_filter assumes Y = [nVolumes, nVoxel] dimensions
% K.row is specified to enable different filtering for different time
% frames e.g. sessions, to not filter drifts between session time gaps
Y = spm_filter(K, Y);

% back-conversion to 4D image
this.data = reshape(Y', nVoxel);

author: translationalneuromodeling

repo: tapas


text:
function [pvec, pstruct] = tapas_softmax_mu3_wld_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2018 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1)       = ptrans(1);         % la_wd
pstruct.la_wd = pvec(1);
pvec(2)       = ptrans(2);         % la_ld
pstruct.la_ld = pvec(2);

return;


text:
function res = glp(im,thresh)

% inputs
% im is the fourier transform of the image
% thresh is the cutoff circle radius

%outputs
% res is the filtered image

[r,c]=size(im);
d0=thresh;

d=zeros(r,c);
h=zeros(r,c);

for i=1:r
    for j=1:c
     d(i,j)=  sqrt( (i-(r/2))^2 + (j-(c/2))^2);
    end
end

for i=1:r
    for j=1:c
      h(i,j)= exp ( -( (d(i,j)^2)/(2*(d0^2)) ) );
    end
end

for i=1:r
    for j=1:c
    res(i,j)=(h(i,j))*im(i,j);

    end
end


text:
function [meanCoherehce,meanAmps,stdCos,semCos] = meanCos(view)
%
% [meanCos,meanAmps] = vectorMeans(view)
%
% Calculates mean correlations and amplitudes (ignoring phase)
% for all scans, for pixels that are in the currently selected
% ROI.
%
% djh 7/98

nscans = numScans(view);

% Get selpts from current ROI
ROIcoords = getCurROIcoords(view);

% Compute vector mean for each scan
meanCoherehce = zeros(1,nscans);
meanAmps = zeros(1,nscans);
stdCos=zeros(1,nscans);

for scanNum = 1:nscans
  [meanCoVM,meanAmp,stdCo,semCo] = meanCo(view,scanNum,ROIcoords);
  meanCoherehce(scanNum) = meanCoVM;
  meanAmps(scanNum) = meanAmp;
  stdCos(scanNum) = stdCo;
  semCos(scanNum) = semCo;

end

author: vistalab


text:
function [adfreq, n, ts, fn, d] = nex_cont(filename, varname)
% nex_cont(filename, varname): Read continuous variable from a .nex file
%
% [adfreq, n, ts, fn, d] = nex_cont(filename, varname)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   varname - variable name
%
%           continuous (a/d) data come in fragments. Each fragment has a timestamp
%           and a number of a/d data points. The timestamp corresponds to
%           the time of recording of the first a/d value in this fragment.
%           All the data values stored in the vector d.
% OUTPUT:
%   n - total number of data points
%   ts - array of fragment timestamps (one timestamp for fragment, in seconds)
%   fn - number of data points in each fragment
%   d - array of a/d values (in millivolts)

% original from Plexon, download from http://www.plexoninc.com (8/4/02)
% modifications by Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

n = 0;
adfreq = 0;
ts = 0;
fn = 0;
d = 0;

if(nargin ~= 2)
   disp('2 input arguments are required')
   return
end

if(ischar(filename) == 0)
   disp('input arguments should be character arrays')
   return
end

if(ischar(varname) == 0)
   disp('input arguments should be character arrays')
   return
end

if(length(filename) == 0)
   [fname, pathname] = uigetfile('*.nex', 'Select a Nex file');
    filename = strcat(pathname, fname);
end

fid = fopen(filename, 'r', 'ieee-le');
if(fid == 0)
   return
end

disp(strcat('file = ', filename));
magic = fread(fid, 1, 'int32');
version = fread(fid, 1, 'int32');
comment = fread(fid, 256, 'char');
freq = fread(fid, 1, 'double');
tbeg = fread(fid, 1, 'int32');
tend = fread(fid, 1, 'int32');
nvar = fread(fid, 1, 'int32');
fseek(fid, 260, 'cof');
name = zeros(1, 64);
found = 0;
for i=1:nvar
    type = fread(fid, 1, 'int32');
    var_version = fread(fid, 1, 'int32');
    name = fread(fid, [1 64], 'char');
    offset = fread(fid, 1, 'int32');
    nf = fread(fid, 1, 'int32');
    dummy = fread(fid, 32, 'char');
    adfreq = fread(fid, 1, 'double');
    adtomv = fread(fid, 1, 'double');
    n = fread(fid, 1, 'int32');
    name = char(name);
    name = deblank(name);
    k = strcmp(name, deblank(varname));
    if(k == 1)
        if type ~= 5
            disp(sprintf('%s is not a continuous variable', deblank(varname)));
            return;
        end
        found = 1;
        fseek(fid, offset, 'bof');
        ts = fread(fid, [1 nf], 'int32');
        fn = fread(fid, [1 nf], 'int32');
        d = fread(fid, [1 n], 'int16');
        break
    end
    dummy = fread(fid, 76, 'char');
end

fclose(fid);

if found == 0
    disp('did not find variable in the file');
else
    ts = ts/freq;
    d = d*adtomv;
    fn(nf+1) = n;
    fn = diff(fn);
    disp(strcat('number of data points = ', num2str(n)));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test_suite=test_assert_exception_thrown()
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

% Test cases where exceptions are thrown and that is OK
function test_assert_exception_thrown_passes
    assertExceptionThrown(@()error('Throw w/o ID'));
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        'moxunit:error');
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        'moxunit:error','message');
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        'moxunit:error');
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        '*','message');                                    % Same as above
    assertExceptionThrown(@()error('Throw w/o ID'),...
        '*','message');                                    % Same as above
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        '*');                                              % Any error OK
    assertExceptionThrown(@()error('Throw w/o ID'),...
        '*');                                              % Same as above

    % Explicitly assert that an error is thrown without an ID
    assertExceptionThrown(@()error('Throw w/o ID'),...
        '');

    % Allow id to be a cellstr
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        {'moxunit:error'},'message');
    assertExceptionThrown(@()error('Throw w/o ID'),...
        '*','message');
    assertExceptionThrown(@()error('moxunit:error','msg'),...
        {'moxunit:foo','moxunit:error','moxunit:foo'});
    assertExceptionThrown(@()error('Throw w/o ID'),...
        {'','moxunit:baz'});

function test_assert_exception_thrown_illegal_arguments
    args_cell={ ...
                {@()error('foo'),'message'},...
                {@()error('foo'),'not:id:_entifier','message'},...
                {@()error('foo'),struct},...
                {@()error('foo'),9},...
                {@()error('foo'),'error:id',9},...
                {@()error('foo'),'error:id',struct},...
                {@()error('foo'),'error_id','message'},...
                {'foo'},...
              };

    for k=1:numel(args_cell)
        args=args_cell{k};
        try
            assertExceptionThrown(args{:})

            error_exception_not_thrown('moxunit:illegalParameter');
        catch
            [unused,error_id]=lasterr();
            error_if_wrong_id_thrown('moxunit:illegalParameter',error_id);
        end
    end

% Test cases where func throws exceptions and we need to throw as well
function test_assert_exception_thrown_wrong_exception
    % Verify that when func throws but the wrong exception comes out, we respond
    % with the correct exception (assertExceptionThrown:wrongException)
    exception_id_cell = {'moxunit:failed',...
                         '',...
                         {'','moxunit:failed'},...
                        };

    for k=1:numel(exception_id_cell)
        exception_id=exception_id_cell{k};

        try
            assertExceptionThrown(@()error('moxunit:error','msg'),...
                exception_id,'msg');
            error_exception_not_thrown('moxunit:wrongExceptionRaised');
        catch
            [unused,error_id]=lasterr();
            error_if_wrong_id_thrown('moxunit:wrongExceptionRaised',...
                                                            error_id);
        end
    end

% Test cases where func does not throw but was expected to do so
function test_assert_exception_thrown_exceptions_not_thrown

    % For all combination of optional arguments we expect the same
    % behavior. We will loop, testing all combinations here
    args_cell = {...
                    {@do_nothing},...                   % No arguments
                    {@do_nothing,'moxunit:failed'},...  % Identifier only
                    {@do_nothing,'*','message'},...     % Wildcard
                    {@do_nothing,'a:b','msg'},...       % All arguments
                    {@do_nothing,{'a:b','c:d'},'msg'}   % Cell str
                 };

   for k=1:numel(args_cell)
        args=args_cell{k};

        % Run the test
        try
            assertExceptionThrown(args{:});
            error_exception_not_thrown('moxunit:exceptionNotRaised');
        catch
            [unused,error_id]=lasterr();
            error_if_wrong_id_thrown('moxunit:exceptionNotRaised',...
                                        error_id);
        end
   end

function error_exception_not_thrown(error_id)
    error('moxunit:exceptionNotRaised', 'Exception ''%s'' not thrown', error_id);

function error_if_wrong_id_thrown(expected_error_id, thrown_error_id)
    if ~strcmp(thrown_error_id, expected_error_id)
        error('moxunit:wrongExceptionRaised',...
              'Exception raised with id ''%s'' expected id ''%s''',...
              thrown_error_id,expected_error_id);
    end

function do_nothing
    % do nothing


text:
function test_suite = test_modelassesment2

%   Run specific demo and save values for comparison.
%
%   See also
%     TEST_ALL, DEMO_MODELASSESMENT2

% Copyright (c) 2011-2012 Ville Tolvanen

initTestSuite;

  function testDemo
    % Set random number stream so that test failing isn't because randomness.
    % Run demo & save test values.
    prevstream=setrandstream(0);

    disp('Running: demo_modelassesment2')
    demo_modelassesment2
    path = which('test_modelassesment2.m');
    path = strrep(path,'test_modelassesment2.m', 'testValues');
    if ~(exist(path, 'dir') == 7)
      mkdir(path)
    end
    path = strcat(path, '/testModelAssesment2');
    save(path, 'DIC', 'DIC2', 'DIC_latent', 'p_eff', ...
      'p_eff2', 'p_eff_latent', 'p_eff_latent2', 'mlpd_cv', 'mlpd_loo');

    % Set back initial random stream
    setrandstream(prevstream);
    drawnow;clear;close all

  % Compare test values to real values.

  function testDICParameters
    values.real = load('realValuesModelAssesment2.mat', 'DIC');
    values.real.DIC(isnan(values.real.DIC)) = 0;
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'),'DIC');
    values.test.DIC(isnan(values.test.DIC)) = 0;
    assertVectorsAlmostEqual(values.real.DIC, values.test.DIC, 'relative', 0.05);

  function testDICAll
    values.real = load('realValuesModelAssesment2.mat', 'DIC2');
    values.real.DIC2(isnan(values.real.DIC2)) = 0;
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'),'DIC2');
    values.test.DIC2(isnan(values.test.DIC2)) = 0;
    assertVectorsAlmostEqual(values.real.DIC2, values.test.DIC2, 'relative', 0.05);

  function testDICLatent
    values.real = load('realValuesModelAssesment2.mat', 'DIC_latent');
    values.real.DIC_latent(isnan(values.real.DIC_latent)) = 0;
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'),'DIC_latent');
    values.test.DIC_latent(isnan(values.test.DIC_latent)) = 0;
    assertVectorsAlmostEqual(values.real.DIC_latent, values.test.DIC_latent, 'relative', 0.05);

  function testPeffLatentsMarginalized
    values.real = load('realValuesModelAssesment2.mat', 'p_eff');
    values.real.p_eff(isnan(values.real.p_eff)) = 0;
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'),'p_eff');
    values.test.p_eff(isnan(values.test.p_eff)) = 0;
    assertVectorsAlmostEqual(values.real.p_eff, values.test.p_eff, 'relative', 0.4);

	function testPeffLatent
    values.real = load('realValuesModelAssesment2.mat', 'p_eff_latent');
    values.real.p_eff_latent(isnan(values.real.p_eff_latent)) = 0;
     values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'),'p_eff_latent');
     values.test.p_eff_latent(isnan(values.test.p_eff_latent)) = 0;
    assertVectorsAlmostEqual(values.real.p_eff_latent, values.test.p_eff_latent, 'relative', 0.05);

	function testPeffLatent2
    values.real = load('realValuesModelAssesment2.mat', 'p_eff_latent2');
    values.real.p_eff_latent2(isnan(values.real.p_eff_latent2)) = 0;
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'),'p_eff_latent2');
    values.test.p_eff_latent2(isnan(values.test.p_eff_latent2)) = 0;
    assertVectorsAlmostEqual(values.real.p_eff_latent2, values.test.p_eff_latent2, 'relative', 0.25);

	function testLogPredDensity10foldCV
    values.real = load('realValuesModelAssesment2.mat', 'mlpd_cv');
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'), 'mlpd_cv');
    assertVectorsAlmostEqual(values.real.mlpd_cv, values.test.mlpd_cv, 'relative', 0.05);

  function testLogPredDensityLOOPRED
    values.real = load('realValuesModelAssesment2.mat', 'mlpd_loo');
    values.test = load(strrep(which('test_modelassesment2.m'), 'test_modelassesment2.m', 'testValues/testModelAssesment2.mat'), 'mlpd_loo');
    assertVectorsAlmostEqual(values.real.mlpd_loo, values.test.mlpd_loo, 'relative', 0.05);


text:
%     Steffen Urban email: steffen.urban@kit.edu
%     Copyright (C) 2014  Steffen Urban
%
%     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.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

% 04.03.2014 by Steffen Urban
% this is a modified file from
% Davide Scaramuzzas Toolbox OcamCalib
% original filename: findcenter.m

function findcenterUrban(calib_data)

if isempty(calib_data.ima_proc) | isempty(calib_data.Xp_abs)
    fprintf(1,'\nNo corner data available. Extract grid corners before calibrating.\n\n');
    return;
end

fprintf(1,'\nComputing center coordinates.\n\n');

if isempty(calib_data.taylor_order),
    calib_data.taylor_order = calib_data.taylor_order_default;
end

%% ================
%  added code
options = optimset('Display','off','MaxIter',10000,'LargeScale','off');
x0   = [calib_data.ocam_model.xc, calib_data.ocam_model.yc];

[x0,~,~,~] = fminsearch(@errCenterUrban, x0, options, calib_data);

calib_data.ocam_model.xc = x0(1);
calib_data.ocam_model.yc = x0(2);
calib_data.ocam_model.c=1;
calib_data.ocam_model.d=0;
calib_data.ocam_model.e=0;

%% do calibration again
[calib_data.RRfin,calib_data.ocam_model.ss] = calibrate(calib_data.Xt, calib_data.Yt, calib_data.Xp_abs, calib_data.Yp_abs, calib_data.ocam_model.xc, calib_data.ocam_model.yc, calib_data.taylor_order, calib_data.ima_proc);
calib_data.calibrated = 1; %This flag i s1 when the camera has been calibrated

% reproject
M=[calib_data.Xt,calib_data.Yt,zeros(size(calib_data.Xt))];
[allerr,rms] = reprojectpoints_advUrban(calib_data.ocam_model, calib_data.RRfin, calib_data.ima_proc, calib_data.Xp_abs, calib_data.Yp_abs, M);
calib_data.rmsAfterCenter = rms;

ss = calib_data.ocam_model.ss;
ss

figure(3);
set(3,'Name','Calibration results','NumberTitle','off');
subplot(2,1,1);
plot(0:floor(calib_data.ocam_model.width/2),polyval([ss(end:-1:1)],[0:floor(calib_data.ocam_model.width/2)])); grid on; axis equal;
xlabel('Distance ''rho'' from the image center in pixels');
ylabel('f(rho)');
title('Forward projection function');
%
subplot(2,1,2);
plot(0:floor(calib_data.ocam_model.width/2),180/pi*atan2(0:floor(calib_data.ocam_model.width/2),-polyval([ss(end:-1:1)],[0:floor(calib_data.ocam_model.width/2)]))-90); grid on;
xlabel('Distance ''rho'' from the image center in pixels');
ylabel('Degrees');
title('Angle of optical ray as a function of distance from circle center (pixels)');

calib_data.calibrated = 1; %This flag is 1 when the camera has been calibrated

author: urbste

repo: ImprovedOcamCalib