text:
function [p, f, t] = spm_mmtspec (x,Fs, freqs,timeres, timestep, NW)
% Moving multitaper based spectrogram
% FORMAT [p, f, t] = spm_mmtspec (x,Fs,freqs,timeres)
%
% x         input time series
% Fs        sampling frequency of input time series
% freqs     desired vector of frequencies for spectrogram eg. [6:1:30]
% timeres   desired time resolution for spectrogram, default T/16
%           where T is duration of x
%
% p         p(f, t) is estimate of power at freq f and time t
%
% Time series is split into a series of overlapping windows with 5% overlap.
% Desired frequency resolution is attained by zero padding
% as/if necessary. The taper approach is applied to each padded sample.
%
% Plot spectrogram using imagesc(t,f,p); axis xy
%___________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Partha Mitra, Ken Harris and Will Penny
% $Id: spm_mmtspec.m 4021 2010-07-28 12:43:16Z vladimir $

nChannels = size(x, 2);
nSamples = size(x,1);

if (nargin<4 || isempty(timeres))
    timeres = nSamples/(16*Fs);
end

if (nargin<5 || isempty(timestep))
    percent_overlap=0.05;
else
    percent_overlap= 1-timestep/timeres;
end

if (nargin<6 || isempty(NW))
   NW=3;
end

if length(unique(diff(freqs))) > 1
    error('Varying frequency resolution is not supported.');
end

df=freqs(2)-freqs(1);
nFFT=round(Fs/df);

WinLength=round(Fs*timeres);
nOverlap=ceil(percent_overlap*WinLength);

nTapers = 2*NW -1;

% Now do some computations that are common to all spectrogram functions

winstep = WinLength - nOverlap;

% check for column vector input
if nSamples == 1
    x = x';
    nSamples = size(x,1);
    nChannels = 1;
end;

if nSamples < WinLength
    disp('Error in spm_mmtspec: win length must be less than number of samples');
    return;
end

% calculate number of FFTChunks per channel
nFFTChunks = round(((nSamples-WinLength)/winstep));
% turn this into time, using the sample frequency
t = winstep*(0:(nFFTChunks-1))'/Fs;

% set up f and t arrays
if ~any(any(imag(x)))    % x purely real
    if rem(nFFT,2),    % nfft odd
        select = [1:(nFFT+1)/2];
    else
        select = [1:nFFT/2+1];
    end
    nFreqBins = length(select);
else
    select = 1:nFFT;
end
f = (select - 1)'*Fs/nFFT;

% allocate memory now to avoid nasty surprises later
y=complex(zeros(nFreqBins, nFFTChunks, nChannels, nChannels)); % output array
Periodogram = complex(zeros(nFreqBins, nTapers, nChannels, nFFTChunks)); % intermediate FFTs
Temp1 = complex(zeros(nFFT, nTapers, nFFTChunks));
Temp2 = complex(zeros(nFFT, nTapers, nFFTChunks));
Temp3 = complex(zeros(nFFT, nTapers, nFFTChunks));
eJ = complex(zeros(nFFT, nFFTChunks));

% calculate Slepian sequences.
%Tapers=spm_dpss(WinLength,NW);
Tapers=spm_dpss(max(nFFT,WinLength),NW);
Tapers=Tapers(:,1:nTapers);

% Vectorized alogrithm for computing tapered periodogram with FFT
TaperingArray = repmat(Tapers, [1 1 nChannels]);
for j=1:nFFTChunks
    Segment = x((j-1)*winstep+[1:WinLength], :);

    if WinLength<nFFT,
        % Zero pad sample to attain desired freq resolution
        Segment=[Segment;zeros(nFFT-WinLength, nChannels)];
    end

    SegmentsArray = permute(repmat(Segment, [1 1 nTapers]), [1 3 2]);
    TaperedSegments = TaperingArray .* SegmentsArray;

    fftOut = fft(TaperedSegments, nFFT);
    Periodogram(:,:,:,j) = fftOut(select,:,:); %fft(TaperedSegments,nFFT);
end

% Now make cross-products of them to fill cross-spectrum matrix
for Ch1 = 1:nChannels
        Temp1 = reshape(Periodogram(:,:,Ch1,:), [nFreqBins,nTapers,nFFTChunks]);
        Temp2 = Temp1;
        Temp2 = conj(Temp2);
        Temp3 = Temp1 .* Temp2;
        eJ=sum(Temp3, 2);
        p(: ,:, Ch1) = eJ/nTapers;
end

% Remove frequencies outside user-specified range (0.1 to correct for small
% numerical differences)
ind=find(f>=(freqs(1)-0.1) & f <=(freqs(end)+0.1));
f=f(ind);
p= p(ind,:, :);

author: spm


text:
function jed = epoch_to_jed_persian ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_PERSIAN returns the epoch of the Persian calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1952062.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function groupAssignment = prtUtilEquallySubDivideData(Y,nDivisions)
% prtUtilEquallySubDivideData  Equally sub-divide into several groups.
%   By equally sub-divide I mean make it so that each group has the same
%   number of data points of each class.
%
% Syntax: groupAssignment = prtUtilEquallySubDivideData(Y,nDivisions)
%
% Inputs:
%   Y - The class label vector for the dataset
%   nDivision - The number of division to make
%
% Outputs:
%   groupAssignment = Integer assignments specifying the group for each
%       datapoint. These are randomly drawn
%

if nDivisions > length(Y)
   warning('The number of requested divisions is larger than the amount of data. The number of divisions was changed to the length of the data.')
   nDivisions = length(Y);
end

% Leave One Out quick method...
if nDivisions == length(Y)
   groupAssignment = randperm(nDivisions)';
   return
end

% Y may contain NaNs. These are treated at missing data. All NaNs are
% distributed evenly across the folds. To do this we make a virtual class
% for nan and then use the same code below.
nanVals = isnan(Y);
if any(nanVals)
    Y(nanVals) = max(Y(~nanVals))+1;
end

sortedY = sort(Y);
nSamples = length(Y);
sortedGroupAssignment = repmat((1:nDivisions)',ceil(nSamples/nDivisions),1);
sortedGroupAssignment = sortedGroupAssignment(1:nSamples);

% Randomize within each class and revert back to the original order
groupAssignment = zeros(size(Y));
uY = unique(Y);
nClasses = length(uY);
for iClass = 1:nClasses
   cSortedGroupAssignment = sortedGroupAssignment(sortedY == uY(iClass));
   groupAssignment(Y==uY(iClass)) = cSortedGroupAssignment(randperm(sum(sortedY == uY(iClass))));
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


function vl_setupnn()
% VL_SETUPNN  Setup the MatConvNet toolbox
%    VL_SETUPNN() function adds the MatConvNet toolbox to MATLAB path.

% Copyright (C) 2014 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).

root = vl_rootnn() ;
addpath(fullfile(root, 'matlab')) ;
addpath(fullfile(root, 'matlab', 'mex')) ;
addpath(fullfile(root, 'matlab', 'xtest')) ;


text:
function YamlStruct = ReadYaml(yaml_file)
% This function reads Yaml file into struct
% Example
% >> yaml_file = 'EnaspolMain.yaml';
% >> YamlStruct = ReadYaml(yaml_file)
%
%  %======================================================================
%{
		Copyright (c) 2011
		This program is a result of a joined cooperation of Energocentrum
		PLUS, s.r.o. and Czech Technical University (CTU) in Prague.
        The program is maintained by Energocentrum PLUS, s.r.o. and
        licensed under the terms of MIT license. Full text of the license
        is included in the program release.

        Author(s):
		Jiri Cigler, Dept. of Control Engineering, CTU Prague
		Jan  Siroky, Energocentrum PLUS s.r.o.

        Implementation and Revisions:

        Auth  Date        Description of change
        ----  ---------   -------------------------------------------------
        jc    01-Mar-11   First implementation
        jc    02-Mar-11   .jar package initialization moved to external fun
        jc    18-Mar-11   Warning added when imported file not found
%}
%======================================================================

import prtExternal.yaml.*;

%InitYaml();
packagePath = fullfile(prtRoot,'+prtExternal','+yaml','YAMLMatlab','external-packages','snakeyaml');
jarPath = fullfile(packagePath,'snakeyaml-1.8.jar');

v = strfind(javaclasspath,packagePath);
if isempty(v) || isempty(v{1})
    javaaddpath(jarPath)
end
import('org.yaml.snakeyaml.Yaml');

yamlreader = Yaml();
yml = fileread(yaml_file);
if isempty(yml)
    YamlStruct = repmat(struct,0,0);
    return;
end
jymlobj = yamlreader.load(yml);

Data = Hash2Struct(jymlobj);
work_folder=fileparts(yaml_file);
if isfield(Data,'import')
    yml = '';
    for i=1:numel(Data.import)
        fToImport= fullfile(work_folder,Data.import{i});
        try
            fxx=fileread(fToImport);
            yml=sprintf('%s\n%s',yml,  fxx);
        catch %#ok<CTCH>
            warning('YAMLMatlab:FileNotFoundException','YAMLMatlab: File %s not found',fToImport);
        end
    end
    jymlobj = yamlreader.load(yml);
    YamlStruct = Hash2Struct(jymlobj);
else
    YamlStruct = Data;
end


text:
function pass = test_integral(pref)
% Test path integral for 3D functions.

if ( nargin == 0)
    pref = chebfunpref;
end
tol = 1e3*pref.cheb3Prefs.chebfun3eps;

f = chebfun3(@(x,y,z) 2*x.*y+3*z);
curve = chebfun(@(t) [2*t 5*t 4*t], [0, 1]);
exact = 38*sqrt(5);
pass(1) = abs(integral(f, curve) - exact)/exact < tol;

f = chebfun3(@(x,y,z) x.^2.*z);
curve = chebfun(@(t) [-t 6+2*t 2+5*t], [0 1]);
exact = 23*sqrt(30)/12;
pass(2) = abs(integral(f, curve) - exact )/exact < tol;

% http://tutorial.math.lamar.edu/Classes/CalcIII/LineIntegralsPtI.aspx
dom = [-1 1 -1 1 0 12*pi];
f = chebfun3(@(x,y,z) x.*y.*z, dom);
curve = chebfun(@(t) [cos(t) sin(t) 3*t], [0 4*pi]);
exact = -3*sqrt(10)*pi;
pass(3) = abs(integral(f, curve) - exact )/exact < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function cache_dir = resolve_cache_dir(model, model_dir, opts)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

cache_dir_base = fullfile(model_dir, 'cache_dir');
mkdir_if_missing(cache_dir_base);
if ~isempty(opts.cache_dir)
    assert(isempty(opts.cache_dir_name));
    cache_dir = opts.cache_dir;
elseif ~isempty(opts.cache_dir_name)
    assert(isempty(opts.cache_dir));
    cache_dir = fullfile(cache_dir_base, opts.cache_dir_name);
else
    if opts.use_detection_svms
        cache_dir = fullfile(cache_dir_base, 'detection_svms');
    else
        cache_dir = fullfile(cache_dir_base, 'softmax');
    end
end
mkdir_if_missing(cache_dir);
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function [ i_min, i_max, j_min, j_max, m, n, nst ] = st_header_read ( ...
  input_filename )

%*****************************************************************************80
%
%% ST_HEADER_READ reads the header of an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the ST file.
%
%    Output, integer I_MIN, I_MAX, the minimum and maximum rows.
%
%    Output, integer J_MIN, J_MAX, the minimum and maximum columns.
%
%    Output, integer M, the number of rows.
%
%    Output, integer N, the number of columns.
%
%    Output, integer NST, the number of nonzeros.
%
  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ST_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'ST_HEADER_READ - Fatal error!' );
    return
  end

  nst = 0;
  i_min = + Inf;
  i_max = - Inf;
  j_min = + Inf;
  j_max = - Inf;

  string = '  %d  %d  %f';

  while ( 1 )

    line = fgets ( input_unit );

    if ( line == -1 )
      break
    end

    [ x, count ] = sscanf ( line, string );

    i = x(1);
    j = x(2);
    aij = x(3);

    nst = nst + 1;
    i_min = min ( i_min, i );
    i_max = max ( i_max, i );
    j_min = min ( j_min, j );
    j_max = max ( j_max, j );

  end

  fclose ( input_unit );

  m = i_max - i_min + 1;
  n = j_max - j_min + 1;

  return
end


text:
function out = isdelta(f)
%ISDELTA   Always returns false, since CLASSICFUNs don't have delta functions.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

    out = false;
end


text:
function van_der_corput_seed_set ( seed )

%*****************************************************************************80
%
%% VAN_DER_CORPUT_SEED_SET sets the seed of the van der Corput sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Johannes van der Corput,
%    Verteilungsfunktionen I & II,
%    Nederl. Akad. Wetensch. Proc.,
%    Volume 38, 1935, pages 813-820, pages 1058-1066.
%
%  Parameters:
%
%    Input, integer SEED, the new value for the seed.  SEED should be
%    nonnegative.  Only the integer part of SEED is used.  SEED = 0
%    is allowed, and returns R = 0.
%
%    Global, integer VAN_DER_CORPUT_SEED, the current seed.
%
  global van_der_corput_BASE
  global van_der_corput_SEED
%
  if ( seed < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VAN_DER_CORPUT_SEED_SET - Fatal error!\n' );
    fprintf ( 1, '  Input value of SEED < 0!\n' );
    fprintf ( 1, '  SEED = %d\n', seed );
    return
  end

  van_der_corput_SEED = floor ( seed );

  return
end


text:
function grid_weight = sparse_grid_weights_cfn ( dim_num, level_max, rule, ...
  point_num, grid_index )

%*****************************************************************************80
%
%% SPARSE_GRID_WEIGHTS_CFN computes sparse grid weights based on a CFN 1D rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Input, integer POINT_NUM, the total number of points in
%    the grids.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of
%    point indices, representing a subset of the product grid of level
%    LEVEL_MAX, representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%
%    Output, real GRID_WEIGHT(POINT_NUM), the weights
%    associated with the sparse grid points.
%
  if ( level_max == 0 )
    grid_weight(1:point_num) = 2.0^dim_num;
    return
  end

  grid_weight(1:point_num) = 0.0;

  level_min = max ( 0, level_max + 1 - dim_num );

  for level = level_min : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 );

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_closed ( dim_num, level_1d );
%
%  The product of the 1D orders gives us the number of points in this grid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  Generate the indices of the points corresponding to the grid.
%
      grid_index2 = multigrid_index_cfn ( dim_num, order_1d, order_nd );
%
%  Compute the weights for this grid.
%
      grid_weight2 = product_weights ( dim_num, order_1d, order_nd, rule );
%
%  Adjust the grid indices to reflect LEVEL_MAX.
%
      grid_index2 = multigrid_scale_closed ( dim_num, order_nd, ...
        level_max, level_1d, grid_index2 );
%
%  Now determine the coefficient.
%
      coeff = r8_mop ( level_max - level ) ...
        * r8_choose ( dim_num - 1, level_max - level );

      for point2 = 1 : order_nd

        for point = 1 : point_num

          if ( all ( ...
            grid_index2(1:dim_num,point2) == grid_index(1:dim_num,point) ...
          ) )
            grid_weight(point) = grid_weight(point) ...
              + coeff * grid_weight2(point2);
            break
          end

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function mfx = spm_cfg_mfx
% SPM Configuration file for MFX
%__________________________________________________________________________
% Copyright (C) 2010-2014 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_cfg_mfx.m 6239 2014-10-13 14:53:48Z guillaume $

%--------------------------------------------------------------------------
% dir Directory
%--------------------------------------------------------------------------
dir         = cfg_files;
dir.tag     = 'dir';
dir.name    = 'Directory';
dir.help    = {'Select a directory where the SPM.mat file containing the specified design matrix will be written.'};
dir.filter = 'dir';
dir.ufilter = '.*';
dir.num     = [1 1];

%--------------------------------------------------------------------------
% spmmat Select SPM.mat
%--------------------------------------------------------------------------
spmmat         = cfg_files;
spmmat.tag     = 'spmmat';
spmmat.name    = 'Select SPM.mat files';
spmmat.help    = {...
    'Select the SPM.mat files that contains first-level designs.'
    'They must have the same number of parameters for each session.'
    ['These are assumed to represent session-specific realisations of ' ...
    '2nd-level effects.']};
spmmat.filter  = 'mat';
spmmat.ufilter = '^SPM\.mat$';
spmmat.num     = [1 Inf];

%--------------------------------------------------------------------------
% ffx Create first-level design
%--------------------------------------------------------------------------
ffx       = cfg_exbranch;
ffx.tag   = 'ffx';
ffx.name  = 'FFX Specification';
ffx.val   = {dir spmmat};
ffx.help  = {'Create FFX multi-session first-level design'};
ffx.prog  = @spm_local_ffx;
ffx.vout  = @vout_ffx;

%--------------------------------------------------------------------------
% spmmat Select SPM.mat
%--------------------------------------------------------------------------
spmmat         = cfg_files;
spmmat.tag     = 'spmmat';
spmmat.name    = 'Select SPM.mat';
spmmat.help    = {...
    'Design and estimation structure after a 1st-level analysis.'};
spmmat.filter  = 'mat';
spmmat.ufilter = '^SPM\.mat$';
spmmat.num     = [1 1];

%--------------------------------------------------------------------------
% con Contrast
%--------------------------------------------------------------------------
con         = cfg_entry;
con.tag     = 'contrast';
con.name    = 'Contrast';
con.help    = {
    'Contrast used to define 2nd level design matrix.'
    'E.g. ones(n,1) contrast where n is the number of sessions/subjects.'
    ['The specification of a contrast that is not ones(n,1) allows, ' ...
    'for example, specified sessions/subjects to be ignored.']
    ['If left empty (default), a contrast ones(n,1) will be specified ' ...
    'automatically at run time.']
    };
con.strtype = 'r';
con.num     = [Inf Inf];
con.val     = {[]};

%--------------------------------------------------------------------------
% spec MFX Specification
%--------------------------------------------------------------------------
spec       = cfg_exbranch;
spec.tag   = 'spec';
spec.name  = 'MFX Specification';
spec.val   = {spmmat con};
spec.help  = {'MFX Specification'};
spec.prog  = @spm_local_mfx;
spec.vout  = @vout_mfx;

%--------------------------------------------------------------------------
% mfx Mixed-effects (MFX) analysis
%--------------------------------------------------------------------------
mfx         = cfg_choice;
mfx.tag     = 'mfx';
mfx.name    = 'Mixed-effects (MFX) analysis';
mfx.help    = {'Mixed-effects (MFX) analysis'};
mfx.values  = {ffx spec};

%==========================================================================
% function out = spm_local_ffx(job)
%==========================================================================
function out = spm_local_ffx(job)
spmmat = job.spmmat;
SPMS = cell(size(spmmat));
for i=1:numel(spmmat)
    load(spmmat{i},'SPM');
    SPMS{i} = SPM;
end

matlabbatch{1}.spm.stats.fmri_spec.dir            = cellstr(job.dir);
matlabbatch{1}.spm.stats.fmri_spec.timing.units   = SPMS{1}.xBF.UNITS;
matlabbatch{1}.spm.stats.fmri_spec.timing.RT      = SPMS{1}.xY.RT;
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t  = SPMS{1}.xBF.T;
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t0 = SPMS{1}.xBF.T0;
switch SPMS{1}.xBF.name
    case 'hrf'
        matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [0 0];
    case 'hrf (with time derivative)'
        matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [1 0];
    case 'hrf (with time and dispersion derivatives)'
        matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [1 1];
    case 'Fourier set'
        matlabbatch{1}.spm.stats.fmri_spec.bases.fourier.length = SPMS{1}.xBF.length;
        matlabbatch{1}.spm.stats.fmri_spec.bases.fourier.order  = SPMS{1}.xBF.order;
    case 'Fourier set (Hanning)'
        matlabbatch{1}.spm.stats.fmri_spec.bases.fourier_han.length = SPMS{1}.xBF.length;
        matlabbatch{1}.spm.stats.fmri_spec.bases.fourier_han.order  = SPMS{1}.xBF.order;
    case 'Gamma functions'
        matlabbatch{1}.spm.stats.fmri_spec.bases.gamma.length = SPMS{1}.xBF.length;
        matlabbatch{1}.spm.stats.fmri_spec.bases.gamma.order  = SPMS{1}.xBF.order;
    case 'Finite Impulse Response'
        matlabbatch{1}.spm.stats.fmri_spec.bases.fir.length = SPMS{1}.xBF.length;
        matlabbatch{1}.spm.stats.fmri_spec.bases.fir.order  = SPMS{1}.xBF.order;
end
matlabbatch{1}.spm.stats.fmri_spec.volt   = SPMS{1}.xBF.Volterra;
matlabbatch{1}.spm.stats.fmri_spec.global = SPMS{1}.xGX.iGXcalc;
if isfield(SPMS{1}.xM,'gMT')
    matlabbatch{1}.spm.stats.fmri_spec.mthresh = SPMS{1}.xM.gMT;
end
if ~isempty(SPMS{1}.xM.VM)
    matlabbatch{1}.spm.stats.fmri_spec.mask = cellstr(SPMS{1}.xM.VM.fname); % can be intersection
end
if strncmp('AR',SPMS{1}.xVi.form,2)
    matlabbatch{1}.spm.stats.fmri_spec.cvi  = 'AR(1)';
elseif strcmp(SPMS{1}.xVi.form,'FAST')
    matlabbatch{1}.spm.stats.fmri_spec.cvi  = 'FAST';
else
    %matlabbatch{1}.spm.stats.fmri_spec.cvi  = 'none';
    disp('Forcing serial correlation modelling to be AR(1).');
    matlabbatch{1}.spm.stats.fmri_spec.cvi  = 'AR(1)';
end

k  = 1;
np = zeros(1,numel(SPMS));
disp(' ');
for i=1:numel(SPMS)
    fprintf('Subject %d:\n',i);
    np(i) = 0;
    n     = cumsum([1 SPMS{i}.nscan]);
    nSess = numel(SPMS{i}.Sess);
    for j=1:nSess
        fprintf(' Session %d:\n',j);
        matlabbatch{1}.spm.stats.fmri_spec.sess(k).scans = cellstr(SPMS{i}.xY.P(n(j):n(j+1)-1,:));
        nCond=numel(SPMS{i}.Sess(j).U);
        for l=1:nCond
            matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).name = SPMS{i}.Sess(j).U(l).name{1};
            matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).onset = SPMS{i}.Sess(j).U(l).ons;
            matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).duration = SPMS{i}.Sess(j).U(l).dur;
            o    = 1;
            nReg = numel(SPMS{i}.Sess(j).U(l).P);
            for m=1:nReg
                fprintf('  Condition %d: Columns %d\n',l,nReg);
                np(i) = np(i)+nReg;
                switch SPMS{i}.Sess(j).U(l).P(m).name
                    case 'time'
                        matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).tmod = SPMS{i}.Sess(j).U(l).P(m).h;
                    case 'none'
                    otherwise
                        matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).pmod(o).name  = SPMS{i}.Sess(j).U(l).P(m).name;
                        matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).pmod(o).param = SPMS{i}.Sess(j).U(l).P(m).P;
                        matlabbatch{1}.spm.stats.fmri_spec.sess(k).cond(l).pmod(o).poly  = SPMS{i}.Sess(j).U(l).P(m).h;
                        o = o + 1;
                end
            end
        end
        for l=1:numel(SPMS{i}.Sess(j).C.name)
            matlabbatch{1}.spm.stats.fmri_spec.sess(k).regress(l).name = SPMS{i}.Sess(j).C.name{l};
            matlabbatch{1}.spm.stats.fmri_spec.sess(k).regress(l).val  = SPMS{i}.Sess(j).C.C(:,l);
        end
        matlabbatch{1}.spm.stats.fmri_spec.sess(k).hpf = SPMS{i}.xX.K(j).HParam;
        k = k + 1;
    end
    disp(' ');
end

% Output number of sessions/conditions/columns for each subject

if any(np-np(1))
    disp('Error in FFX specification: all subject should have same number of parameters');
    return
end

spm_jobman('run',matlabbatch);

out.spmmat{1} = fullfile(job.dir{1},'SPM.mat');

%==========================================================================
% function dep = vout_ffx(job)
%==========================================================================
function dep = vout_ffx(job)
dep = cfg_dep;
dep.sname      = 'SPM.mat File';
dep.src_output = substruct('.','spmmat');
dep.tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});

%==========================================================================
% function out = spm_local_mfx(job)
%==========================================================================
function out = spm_local_mfx(job)
load(job.spmmat{1},'SPM');
c = job.contrast;
if isempty(c)
    n = length(SPM.Sess);
    c = ones(n,1);
end
spm_mfx(SPM,c);
out.spmmat{1} = fullfile(fileparts(job.spmmat{1}),'mfx','SPM.mat');

%==========================================================================
% function dep = vout_mfx(job)
%==========================================================================
function dep = vout_mfx(job)
dep = cfg_dep;
dep.sname      = 'SPM.mat File';
dep.src_output = substruct('.','spmmat');
dep.tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});


text:
function prepare_voc12val_det_imgs(output_img_dir, bbox_mat_filelist, resize_dim)
% PREPARE_VOC12VAL_DET
%   Prepare testing images from detector outputs on VOC12 val set.
% input:
%   output_img_dir: output image folder name, output will be like <output_img_dir>/<category>/<imgs>
%   bbox_mat_filenames: txt files containing a list of bbox mat filenames
%   resize_dim (D): if >0 resize images to DxD.
%

addpath(fullfile(mfilename('fullpath'), '../../'));
global_variables;
cls_names = g_cls_names;

% paths
annotation_path = fullfile(g_pascal3d_root_folder, 'Annotations');
image_path = fullfile(g_pascal3d_root_folder, 'PASCAL/VOCdevkit/VOC2012/JPEGImages');
addpath(fullfile(g_pascal3d_root_folder, 'VDPM'));
addpath(fullfile(g_pascal3d_root_folder, 'PASCAL/VOCdevkit/VOCcode'));

% read ids of train/val set images
VOCinit;
img_set = 'val';
ids = textread(sprintf(VOCopts.imgsetpath, img_set), '%s');
M = numel(ids);

bbox_mat_filenames = importdata(bbox_mat_filelist);
for k = 1:numel(bbox_mat_filenames)
    bbox_mat_filenames{k} = fullfile(g_detection_results_folder, bbox_mat_filenames{k});
end

% avoid duplication
assert(exist(output_img_dir,'dir')==0);
mkdir(output_img_dir);
for cls_idx = 1:numel(cls_names)
    mkdir([output_img_dir '/' cls_names{cls_idx}]);
end

for cls_idx = 1:numel(cls_names)
    cls = cls_names{cls_idx};

    imgfile = fopen(fullfile(output_img_dir, sprintf('%s.txt', cls)),'w');

    % load detection bboxes of all images for class
    detection_filename = bbox_mat_filenames{cls_idx};
    assert(~isempty(strfind(detection_filename, cls)));
    mat = load(detection_filename);
    det_bboxes_all = mat.boxes;

    for i = 1:M
        if mod(i,100) == 0
            fprintf('%s: %d/%d\n', cls, i, M);
        end

        % bboxes for image i, class cls
        det_bboxes = det_bboxes_all{i};

        % load i-th image in VOC val
        img_filename = fullfile(image_path, sprintf('%s.jpg', ids{i}));
        im = imread(img_filename);

        for k = 1:size(det_bboxes,1)
            box = det_bboxes(k,:);
            box = box(1:4);
            w = box(3)-box(1)+1;
            h = box(4)-box(2)+1;
            rect = [box(1),box(2), w, h];

            % write cropped img
            cropped_im = imcrop(im, rect);
            if resize_dim > 0
                copped_im = imresize(cropped_im, [resize_dim, resize_dim]);
            end
            cropped_im_filename = sprintf('%s/%s/%s_%s_%s.jpg', output_img_dir, cls, cls, ids{i}, num2str(k));
            imwrite(cropped_im, cropped_im_filename);
            fprintf(imgfile, '%s\n', cropped_im_filename);
        end
    end
    fclose(imgfile);
end

author: ShapeNet

repo: RenderForCNN


text:
function geometry_test0183 ( )

%*****************************************************************************80
%
%% TEST0183 tests CIRCLE_LLR2IMP_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  p_hi =  10.0;
  p_lo = -10.0;
  test_num = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0183\n' );
  fprintf ( 1, '  CIRCLE_LLR2IMP_2D is given:\n' );
  fprintf ( 1, '  a line through P1 and P2,\n' );
  fprintf ( 1, '  a line through Q1 and Q2,\n' );
  fprintf ( 1, '  and a radius R,\n' );
  fprintf ( 1, '  and determines the centers C of 4 circles\n' );
  fprintf ( 1, '  of the given radius, tangent to both lines.\n' );

  seed = 123456789;

  for test = 1 : test_num

    [ p1(1:2,1), seed ] = r8vec_uniform ( 2, p_lo, p_hi, seed );
    [ p2(1:2,1), seed ] = r8vec_uniform ( 2, p_lo, p_hi, seed );
    [ q1(1:2,1), seed ] = r8vec_uniform ( 2, p_lo, p_hi, seed );
    [ q2(1:2,1), seed ] = r8vec_uniform ( 2, p_lo, p_hi, seed );

    r_lo = 1.0;
    r_hi = 5.0;
    [ r, seed ] = r8_uniform ( r_lo, r_hi, seed );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Radius R = %f\n', r );

    fprintf ( 1, '  Point #P1: %f  %f\n', p1(1:2,1) );
    fprintf ( 1, '  Point #P2: %f  %f\n', p2(1:2,1) );
    fprintf ( 1, '  Point #Q1: %f  %f\n', q1(1:2,1) );
    fprintf ( 1, '  Point #Q2: %f  %f\n', q2(1:2,1) );

    pc = circle_llr2imp_2d ( p1, p2, q1, q2, r );

    fprintf ( 1, '  Center #1: %f  %f\n', pc(1:2,1) );
    fprintf ( 1, '  Center #2: %f  %f\n', pc(1:2,2) );
    fprintf ( 1, '  Center #1: %f  %f\n', pc(1:2,3) );
    fprintf ( 1, '  Center #2: %f  %f\n', pc(1:2,4) );

    d1 = line_exp_point_dist_2d ( p1, p2, pc(1:2,1) );
    d2 = line_exp_point_dist_2d ( p1, p2, pc(1:2,2) );
    d3 = line_exp_point_dist_2d ( p1, p2, pc(1:2,3) );
    d4 = line_exp_point_dist_2d ( p1, p2, pc(1:2,4) );

    fprintf ( 1, '  %f  %f  %f  %f\n', d1, d2, d3, d4 );

    d1 = line_exp_point_dist_2d ( q1, q2, pc(1:2,1) );
    d2 = line_exp_point_dist_2d ( q1, q2, pc(1:2,2) );
    d3 = line_exp_point_dist_2d ( q1, q2, pc(1:2,3) );
    d4 = line_exp_point_dist_2d ( q1, q2, pc(1:2,4) );

    fprintf ( 1, '  %f  %f  %f  %f\n', d1, d2, d3, d4 );

  end

  return
end


text:
function [ y, m, d ] = day_borrow_roman ( y, m, d )

%*****************************************************************************80
%
%% DAY_BORROW_ROMAN borrows days from months in a Roman date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, a year, month, and day
%    representing a date.  On input, D might be negative.  On output,
%    M should have decreased by one month, and D gone up by the
%    number of days in the month we "cashed in".  Y may be affected
%    if the input value of M was 1.
%
  while ( d <= 0 )

    m = m - 1;

    [ y, m ] = month_borrow_roman ( y, m );

    days = month_length_roman ( y, m );

    d = d + days;

  end

  return
end


text:
function dat=nj_tseries(ncRef,var,loni,lati,varargin)
% NJ_TSERIES Interpolate time series at lon,lat points from CF model output
% Usage:  dat=nj_tseries(uri,var,lon,lat,varargin)
% Inputs:
%    ncRef = netcdf file name, dods url, or ncgeodataset object
%    var = variable name (string) (e.g. 'salt');
%    lon = vector of longitudes where time series are desired
%    lat = vector of latitudes  where time series are desired
% Optional additional parameter pairs
%    'method' =   'nearest' or 'linear'
%    'layer' =  number of layer to extract
%    'itime' =   indices of time to extract
%    'ele' =    connectivity array (required for non structured grids)
%
% Unstructured Grid Example:  surface layer salinity for 1st 100 time steps
% url='http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3'
% dat=nj_tseries(url,'salinity',-70.5,42.5,'method','nearest',...
%     'layer',1,'itime',1:100,'ele',1);
%
% Structured Grid Example:  surface layer salinity for 1st 100 time steps
% url='http://geoport.whoi.edu/thredds/dodsC/coawst_2_2/fmrc/coawst_2_2_best.ncd'
% dat=nj_tseries(url,'salt',-70.5,42.5,'method','nearest',...
%     'layer',16,'itime',1:100);

% Outputs:
%      dat.vals = data values
%      dat.time = time values
%      dat.ii= i value
%      dat.jj= j value

% Rich Signell (rsignell@usgs.gov)

if (isa(ncRef, 'ncgeodataset')) %check for ncgeodataset Object
  nc = ncRef;
else
  % open CF-compliant NetCDF File as a Common Data Model (CDM) "Grid Dataset"
  nc = ncgeodataset(ncRef);
end
vart=nc.geovariable(var);
lon=nc.data(vart.getlonname);
lat=nc.data(vart.getlatname);

while ~isempty(varargin),
  switch lower(varargin{1}),
    case 'method',
      method=varargin{2};
    case 'layer',
      layer=varargin{2};
    case 'itime',
      itime=varargin{2};
    case 'ele',
      ele=varargin{2};
    otherwise,
      error(['Can''t understand property:' varargin{1}]);
  end
  varargin([1 2])=[];
end
if ~exist('method','var'),method='nearest';end
jdmat=nj_time(ncRef,var); % returns [] if no time variable exists
if ~exist('itime','var'),
  if ~isempty(jdmat)
    itime=1:length(jdmat);
  end
end

ntime=length(itime);
npoint=length(loni);
% initialize
u=zeros(ntime,npoint);
lono=zeros(npoint);
lato=zeros(npoint);
inear=zeros(npoint);
jnear=zeros(npoint);

switch method
  case 'nearest'  % find closest point
    for k=1:npoint
      if exist('ele','var');  % unstructured grid
        ind=nearxy(double(lon(:)),double(lat(:)),loni(k),lati(k));
        lono(k)=lon(ind);
        lato(k)=lat(ind);
        fprintf('extracting point %d:lon=%f,lat=%f\n',k,lono(k),lato(k));

        if ~exist('itime','var');
          u(:,k)=nc{var}(ind);  % no time variable, e.g. bathymetry
        elseif (exist('itime','var') && ~exist('layer','var'))
          u(:,k)=squeeze(nc{var}(itime,ind));
        elseif (exist('itime','var') &&  exist('layer','var'))
          u(:,k)=squeeze(nc{var}(itime,layer,ind));
          dat.layer=layer;
        else
          disp('invalid arguments');return
        end
        dat.ii(k)=ind;
      else   % structured grid
        if isvector(lon), [lon,lat]=meshgrid(lon,lat);end;
        ind=nearxy(double(lon(:)),double(lat(:)),loni(k),lati(k));
        [jj,ii]=ind2ij(double(lon),ind);
        lono(k)=lon(jj,ii);
        lato(k)=lat(jj,ii);
        fprintf('extracting point %d:lon=%f,lat=%f\n',k,lono(k),lato(k));
        if ~exist('itime','var');
          u(:,k)=nc{var}(jj,ii);  % no time variable, e.g. bathymetry
        elseif (exist('itime','var') && ~exist('layer','var'))
          u(:,k)=squeeze(nc{var}(itime,jj,ii));
        elseif (exist('itime','var') &&  exist('layer','var'))
          u(:,k)=squeeze(nc{var}(itime,layer,jj,ii));
          dat.layer=layer;
        else
          disp('invalid arguments');return
        end
        dat.ii(k)=ii;
        dat.jj(k)=jj;
      end
    end
  case 'linear'  % interpolate from surrounding points
    if exist('ele','var');  % unstructured grid
      error(' ''linear'' interpolation not working yet for unstructured grid. Use ''nearest'' ')
    else % structured grid
      lono=loni;
      lato=lati;
      if isvector(lon), [lon,lat]=meshgrid(lon,lat);end;
      [n,m]=size(lon);
      [ii,jj]=meshgrid(1:m,1:n);
      F1=TriScatteredInterp(double(lon(:)),double(lat(:)),ii(:));
      ivar=F1(loni(:),lati(:));
      F2=TriScatteredInterp(double(lon(:)),double(lat(:)),jj(:));
      jvar=F2(loni(:),lati(:));
      for k=1:length(loni);
        fprintf('interpolating at point %d:lon=%f,lat=%f\n',k,loni(k),lati(k));
        i0=floor(ivar(k));
        ifrac=ivar(k)-i0;
        j0=floor(jvar(k));
        jfrac=jvar(k)-j0;
        if ~exist('itime','var'); % 2D (no time dimension)
          uall=nc{var}(j0:j0+1,i0:i0+1);
          u1=uall(1,1);
          u2=uall(1,2);
          u3=uall(2,1);
          u4=uall(2,2);
        elseif (exist('itime','var') && ~exist('layer','var'))
          uall=squeeze(nc{var}(itime,j0:j0+1,i0:i0+1));
          u1=uall(:,1,1);
          u2=uall(:,1,2);
          u3=uall(:,2,1);
          u4=uall(:,2,2);
        elseif (exist('itime','var') &&  exist('layer','var'))
          uall=squeeze(nc{var}(itime,layer,j0:j0+1,i0:i0+1));
          u1=uall(:,1,1);
          u2=uall(:,1,2);
          u3=uall(:,2,1);
          u4=uall(:,2,2);
          dat.layer=layer;
        else               % 3D (2D with time dimension)
          disp('invalid arguments');return
        end
        dat.ii=ivar;
        dat.jj=jvar;
        ua=u1+(u2-u1)*ifrac;
        ub=u3+(u4-u3)*ifrac;
        u(:,k)=ua+jfrac*(ub-ua);
      end
    end
end
dat.vals=u;
dat.lon=lono;
dat.lat=lato;
if ~isempty(jdmat)
  dat.time=jdmat(itime);
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function hRD = plugin_rd_tx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9364z7020.common.plugin_rd('ccbob_cmos','Tx');
AnalogDevices.adrv9364z7020.common.add_io(hRD, 'Tx');

author: analogdevicesinc


text:
function X = arrange(X,foo)
%ARRANGE Arranges the rank-1 components of a ktensor.
%
%   ARRANGE(X) normalizes the columns of the factor matrices and then sorts
%   the ktensor components by magnitude, greatest to least.
%
%   ARRANGE(X,N) absorbs the weights into the Nth factor matrix instead of
%   lambda.
%
%   ARRANGE(X,P) rearranges the components of X according to the
%   permutation P. P should be a permutation of 1 to NCOMPOMENTS(X).
%
%   See also KTENSOR, NCOMPONENTS.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Just rearrange and return if second argument is a permutation
if exist('foo','var') && (length(foo) > 1)
    X.lambda = X.lambda(foo);
    for i = 1 : ndims(X)
        X.u{i} = X.u{i}(:,foo);
    end
    return;
end

%% Ensure that matrices are normalized
X = normalize(X);

%% Sort
[X.lambda, idx] = sort(X.lambda, 1, 'descend');
for i = 1 : ndims(X)
    X.u{i} = X.u{i}(:,idx);
end

%% Absorb the weight into one factor, if requested
if exist('foo','var')
    r = length(X.lambda);
    X.u{end} = X.u{end} * spdiags(X.lambda,0,r,r);
    X.lambda = ones(size(X.lambda));
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function aux_depthcontour(params, hParentFigure)
% function aux_FMD(params, hParentFigure);
%-------------------------------------------
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% Thomas van Stiphout, thomas@sed.ethz.ch
% last update: 7.9.2005

report_this_filefun(mfilename('fullpath'));

% prepare Value that will be illustrated on the 3D surface defined above
 params.mPlotValues=ones(length(params.vUsedNodes),1)*nan;
 params.mPlotValues(params.vUsedNodes)=params.mValueGrid(:,1);
 vSel=~isnan(params.mValueGrid(:,1));
 vTmp=params.mValueGrid(:,1);
 vTmp(~vSel)=nan;
 params.mPlotValues(params.vUsedNodes)=vTmp;
 %       params.mPlotValues(params.mPlotValues == 0)=nan;
 params.mPlotValues=reshape(params.mPlotValues,length(params.vY),length(params.vX));
 mZ_=params.mZ;
 mZ_(isnan(params.mPlotValues))=nan;
 V_=-10:-20:-190;
 hold on;contour3(params.mX,params.mY,-mZ_,V_,'k--');
 V_=-0:-20:-200;
 hold on;contour3(params.mX,params.mY,-mZ_,V_,'k-');
 hold on;contour(params.mX,params.mY,-mZ_,'ShowText','on');
%  figure;surf(params.mX,params.mY,-mZ_,params.mPlotValues);
      view(3)

author: CelsoReyes


text:
function [ D ] = gsp_grad_mat( G )
%GSP_GRAD_MAT Gradient sparse matrix of the graph G
%   Usage:  D = gsp_gradient_mat(G);
%
%   Input parameters:
%       G   : Graph structure
%
%   Output parameters:
%       D   : Gradient sparse matrix
%
%   This function return the gradient matrix. To be more effiecient, call
%   the function::
%
%           G = gsp_adj2vec(G)
%
%   before this function.
%
%   Example:::
%
%       N = 40;
%       G = gsp_sensor(N);
%       G = gsp_adj2vec(G);
%       D = gsp_grad_mat(G);
%

% Author: Nathanael Perraudin
% Date  : 14 Mai 2014
% Testing: test_operators

if ~isfield(G,'v_in')
    G = gsp_adj2vec(G);
    warning(['GSP_GRADIENT_MAT: To be more efficient you should run: ',...
        'G = gsp_adj2vec(G); before using this proximal operator.']);
end

% if isfield(G,'Diff');
%     D = G.Diff;
%     return;
% end

% In case the graph has logical in the weights matrix
G.weights = double(G.weights);

if strcmp(G.lap_type,'combinatorial')
    n = G.Ne;
    Dr = [1:n 1:n];
    Dc(1:n) = G.v_in;
    Dc(n+1:2*n) = G.v_out;
    Dv(1:n) = sqrt(G.weights);
    Dv(n+1:2*n) = -sqrt(G.weights);
elseif strcmp(G.lap_type,'normalized')
    n = G.Ne;
    Dr = [1:n 1:n];
    Dc(1:n) = G.v_in;
    Dc(n+1:2*n) = G.v_out;
    Dv(1:n) = sqrt(G.weights./G.d(G.v_in));
    Dv(n+1:2*n) = -sqrt(G.weights./G.d(G.v_out));
%     [v_i, v_j, weights] = find(G.W);
%     n = length(v_i);
%     Dr = [1:n 1:n];
%     Dc(1:n) = v_i;
%     Dc(n+1:2*n) = v_j;
%     Dv(1:n) = sqrt(weights./G.d(v_i));
%     Dv(n+1:2*n) = -sqrt(weights./G.d(G.v_j));
else
    error('Not implemented yet!')
end

D = sparse(Dr,Dc,Dv,n,G.N);

end

repo: gspbox


text:
function [gradS,gradV,gradAlpha] = refineShape(S,V,alpha,stateFiles,lambda,shapePrior,deformationPrior,shapePriorInstance,iter)

% V is (numpoints*3 X K)
% alpha is K X length(fnames)

%% Initializations
params = get_params();

%numcores = max(matlabpool('size'),1);

poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    numcores = 0;
else
    numcores = poolobj.NumWorkers;
end

gradICPS = zeros(size(S));
gradICPV = zeros(size(V));
gradICPAlpha = zeros(size(alpha));
numpoints = size(S,1);
K = size(V,2);
numimages = length(stateFiles);

gOS = zeros([size(S) numcores]);
gOV = zeros([size(V) numcores]);
gOA = zeros([size(alpha) numcores]);
cO = zeros([numpoints,1,numcores]);

gSS = zeros([size(S) numcores]);
gSV = zeros([size(V) numcores]);
gSA = zeros([size(alpha) numcores]);

gPS = zeros([size(S) numcores]);
gPV = zeros([size(V) numcores]);
gPA = zeros([size(alpha) numcores]);

gKS = zeros([size(S) numcores]);
gKV = zeros([size(V) numcores]);
gKA = zeros([size(alpha) numcores]);

normS = sqrt(sum(S.^2,2));
directions = S./repmat(normS,1,3);

%% Parallely computing gradients

%h = tic;
localParams = params;
parfor c = 1:numcores
    [gKS(:,:,c),gKV(:,:,c),gKA(:,:,c),gOS(:,:,c),gOV(:,:,c),gOA(:,:,c),cO(:,:,c),...
        gSS(:,:,c),gSV(:,:,c),gSA(:,:,c),gPS(:,:,c),gPV(:,:,c),gPA(:,:,c)]...
        =computeGrads(S,V,alpha,stateFiles,lambda,shapePriorInstance,c,numcores,localParams);
end

gradOcclusionS = sum(gOS,3);
gradOcclusionV = sum(gOV,3);
gradOcclusionAlpha = sum(gOA,3);
occlusionCount = sum(cO,3);

gradSilS = sum(gSS,3);
gradSilV = sum(gSV,3);
gradSilAlpha = sum(gSA,3);

gradShapePriorS = sum(gPS,3);
gradShapePriorV = sum(gPV,3);
gradShapePriorAlpha = sum(gPA,3);

gradKeypointS = sum(gKS,3);
gradKeypointV = sum(gKV,3);
gradKeypointAlpha = sum(gKA,3);

%% Computing soft occlusion grad for the shape
gradOcclusionS = repmat(double((occlusionCount/numimages-0.02)>0),1,3).*gradOcclusionS; %hacking here

%% Computing gradients from priors
%h = tic;
gradPriorS = lambda(3)*shapePrior(S);
[gradPriorV,deformationPriorGradAlpha] = deformationPrior(V,alpha);
gradPriorV = lambda(4)*gradPriorV;
deformationPriorGradAlpha = lambda(4)*deformationPriorGradAlpha;
%fprintf('Prior computation: %.3fs\n',toc(h));

%% Display all gradients
fprintf('%d\t%3.3f\t%3.3f\t%3.3f\t%3.3f\t%3.3f\t%3.3f\t',...
iter, mean(sqrt(sum(mean(gKS,3).^2,2))),mean(sqrt(sum(mean(gOS,3).^2,2)))...
    ,mean(sqrt(sum(mean(gSS,3).^2,2))),mean(sqrt(sum(mean(gPS,3).^2,2))),...
    mean(sqrt(sum(gradPriorS.^2,2))),mean(sqrt(sum(gradPriorV.^2,2))));

%% Computing overall change in shape
gradS = gradICPS + gradOcclusionS + gradPriorS + gradSilS + gradShapePriorS + gradKeypointS;
gradV = gradICPV + gradOcclusionV + gradPriorV + gradSilV + gradShapePriorV + gradKeypointV;
gradAlpha = gradICPAlpha + gradOcclusionAlpha + deformationPriorGradAlpha + gradSilAlpha...
    + gradShapePriorAlpha + gradKeypointAlpha;

%% Computing grad along radial direction to preserve mesh
if(params.opt.spherical)
    projections = sum(gradS.*directions,2);
    gradS = directions.*repmat(projections,1,3);
    for k=1:K
        deformation = reshape(gradV(:,k),numpoints,3);
        projections = sum(deformation.*directions,2);
        deformation = directions.*repmat(projections,1,3);
        gradV(:,k)=deformation(:);
    end
end

end

repo: CategoryShapes


text:
function A = makehatch_plus(hatch,n,m)
%MAKEHATCH_PLUS Predefined hatch patterns
%
% Modification of MAKEHATCH to allow for selection of matrix size. Useful whe using
%   APPLYHATCH_PLUS with higher resolution output.
%
% input (optional) N    size of hatch matrix (default = 6)
% input (optional) M    width of lines and dots in hatching (default = 1)
%
%  MAKEHATCH_PLUS(HATCH,N,M) returns a matrix with the hatch pattern for HATCH
%   according to the following table:
%      HATCH        pattern
%     -------      ---------
%        /          right-slanted lines
%        \          left-slanted lines
%        |          vertical lines
%        -          horizontal lines
%        +          crossing vertical and horizontal lines
%        x          criss-crossing lines
%        .          square dots
%        c          circular dots
%        w          Just a blank white pattern
%        k          Just a totally black pattern
%
%  See also: APPLYHATCH, APPLYHATCH_PLUS, APPLYHATCH_PLUSCOLOR, MAKEHATCH

%  By Ben Hinkle, bhinkle@mathworks.com
%  This code is in the public domain.

% Modified Brian FG Katz    8-aout-03
% Modified David M Kaplan    19-fevrier-08

if ~exist('n','var'), n = 6; end
if ~exist('m','var'), m = 1; end
n=round(n);

switch (hatch)
  case '\'
    [B,C] = meshgrid( 0:n-1 );
    B = B-C;
    clear C
    A = abs(B) <= m/2;
    A = A | abs(B-n) <= m/2;
    A = A | abs(B+n) <= m/2;
  case '/'
    A = fliplr(makehatch_plus('\',n,m));
  case '|'
    A=zeros(n);
    A(:,1:m) = 1;
  case '-'
    A = makehatch_plus('|',n,m);
    A = A';
  case '+'
    A = makehatch_plus('|',n,m);
    A = A | A';
  case 'x'
    A = makehatch_plus('\',n,m);
    A = A | fliplr(A);
  case '.'
    A=zeros(n);
    A(1:2*m,1:2*m)=1;
  case 'c'
    [B,C] = meshgrid( 0:n-1 );
    A = sqrt(B.^2+C.^2) <= m;
    A = A | fliplr(A) | flipud(A) | flipud(fliplr(A));
  case 'w'
    A = zeros(n);
  case 'k'
    A = ones(n);
  otherwise
    error(['Undefined hatch pattern "' hatch '".']);
end

author: Sable

repo: mcbench-benchmarks


text:
%% Convex Hull demo
% We learn how to get hull contours and draw them.
%
% In this sample you will learn how to use the OpenCV function
% <matlab:doc('cv.convexHull') cv.convexHull>.
%
% Sources:
%
% * <https://docs.opencv.org/3.4.0/d7/d1d/tutorial_hull.html>
% * <https://github.com/opencv/opencv/blob/3.4.0/samples/cpp/tutorial_code/ShapeDescriptors/hull_demo.cpp>
%

function varargout = hull_demo_gui(im)
    % load source image
    if nargin < 1
        src = cv.imread(fullfile(mexopencv.root(),'test','stuff.jpg'));
    elseif ischar(im)
        src = imread(im);
    else
        src = im;
    end

    % Convert image to gray and blur it
    if size(src,3) == 3
        src = cv.cvtColor(src, 'RGB2GRAY');
    end
    src = cv.blur(src, 'KSize',[5 5]);

    % create the UI
    h = buildGUI(src);
    if nargout > 0, varargout{1} = h; end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    thresh = round(get(h.slid, 'Value'));
    set(h.txt, 'String',sprintf('Threshold: %3d',thresh));

    % Detect edges using Threshold
    threshold_output = cv.threshold(h.src, thresh, ...
        'MaxValue',255, 'Type','Binary');

    % Find contours
    contours = cv.findContours(threshold_output, ...
        'Mode','Tree', 'Method','Simple');

    % Find the convex hull object for each contour
    hull = cell(size(contours));
    for i=1:numel(contours)
        hull{i} = cv.convexHull(contours{i}, 'Clockwise',false);
    end

    % Draw contours + hull results
    drawing = zeros([size(threshold_output) 3], 'uint8');
    for i=1:numel(contours)
        clr = randi([0 255], [1 3], 'uint8');
        drawing = cv.drawContours(drawing, contours, ...
            'ContourIdx',i-1, 'MaxLevel',0, ...
            'Color',clr, 'Thickness',1, 'LineType',8);
        drawing = cv.drawContours(drawing, hull, ...
            'ContourIdx',i-1, 'MaxLevel',0, ...
            'Color',clr, 'Thickness',2, 'LineType',8);
    end

    % show result
    set(h.img, 'CData',drawing);
    drawnow;
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    thresh = 100;
    max_thresh = 255;
    sz = size(img);
    sz(2) = max(sz(2), 250);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Convex Hull Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+29]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 30 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    h.txt = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[5 5 130 20], 'String',sprintf('Threshold: %3d',thresh));
    h.slid = uicontrol('Parent',h.fig, 'Style','slider', 'Value',thresh, ...
        'Min',0, 'Max',max_thresh, 'SliderStep',[1 10]./(max_thresh-0), ...
        'Position',[135 5 sz(2)-135-5 20]);

    % hook event handlers, and trigger default start
    set(h.slid, 'Callback',{@onChange,h}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [res,K] = dotk(x,y,n)
%DOTK Dot product in K-fold precision.
%   dotk(x,y) computes the scalar product x'*y of the vectors x and y as in
%   K-fold precision, where K is chosen adaptively so that the result is
%   accurate up to machine precision. If x and y are N-D arrays, dotk(x,y)
%   operates along their first non-singleton dimension.
%
%   dotk(x,y,n) returns the dot product of x and y in the dimension n.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] T. Ogita, S. M. Rump, S. Oishi, "Accurate sum and dot product",
%       SIAM J. Sci. Comp., Vol. 26, No. 6, 2005, pp. 1955-1988.

% Permute x and y so that the n-th mode is the first one.
sx = size(x);
sy = size(y);
if any(sx ~= sy), error('dotk:xy','x and y should have equal size.'); end
if nargin < 3, n = find(sx > 1,1); end
if isempty(n) || size(x,n) == 1, res = conj(x).*y; K = 0; return; end

% Call recursively if either x or y is complex.
if ~isreal(x) && isreal(y)
    [resr,Kr] = dotk(real(x),y,n);
    [resi,Ki] = dotk(-imag(x),y,n);
    res = resr+resi*1i;
    K = Kr+Ki*1i;
    return;
elseif isreal(x) && ~isreal(y)
    [resr,Kr] = dotk(x,real(y),n);
    [resi,Ki] = dotk(x,imag(y),n);
    res = resr+resi*1i;
    K = Kr+Ki*1i;
    return;
elseif ~isreal(x) && ~isreal(y)
    x = cat(n,real(x),imag(x));
    [resr,Kr] = dotk(x,cat(n,real(y),imag(y)),n);
    [resi,Ki] = dotk(x,cat(n,imag(y),-real(y)),n);
    res = resr+resi*1i;
    K = Kr+Ki*1i;
    return;
end
if n > 1
    x = permute(x,[n 1:n-1 n+1:ndims(x)]);
    y = permute(y,[n 1:n-1 n+1:ndims(y)]);
end

% Initialize the algorithm.
res = zeros([2*sx(n) sx(1:n-1) sx(n+1:length(sx))]);
factor = 2^ceil(log2(2/eps(class(x)))/2)+1;

% Adaptive K-fold precision dot product.
[res(end,:),res(1,:)] = twoproduct(x(1,:),y(1,:),factor);
for i = 2:size(x,1)
    [tmp,res(i,:)] = twoproduct(x(i,:),y(i,:),factor);
    [res(end,:),res(i+size(x,1)-1,:)] = twosum(res(end,:),tmp);
end
[res,K] = sumk(res);
if n > 1, res = permute(res,[2:n 1 n+1:length(sx)]); end
if isreal(K), K = K+1; else K = K+1+1*1i; end

end

function [x,y] = twoproduct(a,b,factor)

tmp = factor*a;
a1 = tmp-(tmp-a);
a2 = a-a1;
tmp = factor*b;
b1 = tmp-(tmp-b);
b2 = b-b1;

x = a.*b;
y = a2.*b2-(((x-a1.*b1)-a2.*b1)-a1.*b2);

end

function [x,y] = twosum(a,b)

x = a+b;
tmp = x-a;
y = (a-(x-tmp))+(b-tmp);

end

repo: lrslibrary


text:
function  CMETAA  = readCMETAA(fid)
%
% <<<<<<<<<< CLASSIFICATION: UNCLASSIFIED >>>>>>>>>>
%
% READCMETAA Create structure for metadata fields.
%   CMETAA = READCMETAA(FID) returns a structure of fields and their
%   associated values.

% Written by Matt Donath, NGA, matthew.b.donath@nga.ic.gov
% References:
%  > MIL-STD-2500A, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.0
%  > MIL-STD-2500C, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.1
%  > STDI-0001,     NATIONAL SUPPORT DATA EXTENSIONS (SDE)(VERSION 1.3/CN2)
%                   FOR THE NATIONAL IMAGERY TRANSMISSION FORMAT (NITF)
%  > STDI-0002,     THE COMPENDIUM OF CONTROLLED EXTENSIONS (CE) FOR THE
%                   NATIONAL IMAGERY TRANSMISSION FORMAT (NITF) VERSION 2.1
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Create the metadata structure for CMETAA.

CMETAA.CETAG = fread(fid,6,'uint8=>char')';
CMETAA.CEL = str2double(fread(fid,5,'uint8=>char')');

CMETAA.RELATED_TRES = strtrim(fread(fid,2,'uint8=>char')');
CMETAA.RELATED_TRES = strtrim(fread(fid,120,'uint8=>char')');
CMETAA.RD_PRC_NO = strtrim(fread(fid,12,'uint8=>char')');
CMETAA.IF_PROCESS = strtrim(fread(fid,4,'uint8=>char')');
CMETAA.RD_CEN_FREQ = strtrim(fread(fid,4,'uint8=>char')');
CMETAA.RD_MODE = strtrim(fread(fid,5,'uint8=>char')');
CMETAA.RD_PATCH_NO = str2double(fread(fid,4,'uint8=>char')');

CMETAA.CMPLX_DOMAIN = strtrim(fread(fid,5,'uint8=>char')');
CMETAA.CMPLX_MAG_REMAP_TYPE = strtrim(fread(fid,4,'uint8=>char')');
CMETAA.CMPLX_LIN_SCALE = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CMPLX_AVG_POWER = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CMPLX_LINLOG_TP = str2double(fread(fid,5,'uint8=>char')');
CMETAA.CMPLX_PHASE_QUANT_FLAG = strtrim(fread(fid,3,'uint8=>char')');
CMETAA.CMPLX_PHASE_QUANT_BIT_DEPTH = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_SIZE_1 = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_IC_1 = strtrim(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_SIZE_2 = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_IC_2 = strtrim(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_IC_BPP = str2double(fread(fid,5,'uint8=>char')');
CMETAA.CMPLX_WEIGHT = strtrim(fread(fid,3,'uint8=>char')');
CMETAA.CMPLX_AZ_SLL = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_RNG_SLL = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_AZ_TAY_NBAR = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_RNG_TAY_NBAR = str2double(fread(fid,2,'uint8=>char')');
CMETAA.CMPLX_WEIGHT_NORM = strtrim(fread(fid,3,'uint8=>char')');
CMETAA.CMPLX_SIGNAL_PLANE = strtrim(fread(fid,1,'uint8=>char')');

CMETAA.IF_DC_SF_ROW = str2double(fread(fid,6,'uint8=>char')');
CMETAA.IF_DC_SF_COL = str2double(fread(fid,6,'uint8=>char')');

for lp = 1:4
    CMETAA.IF_PATCH_ROW(lp) = str2double(fread(fid,6,'uint8=>char')');
    CMETAA.IF_PATCH_COL(lp) = str2double(fread(fid,6,'uint8=>char')');
end

CMETAA.IF_DC_IS_ROW = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_DC_IS_COL = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_IMG_ROW_DC = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_IMG_COL_DC = str2double(fread(fid,8,'uint8=>char')');

for lp = 1:4
    CMETAA.IF_TILE_ROW(lp) = str2double(fread(fid,6,'uint8=>char')');
    CMETAA.IF_TILE_COL(lp) = str2double(fread(fid,6,'uint8=>char')');
end

CMETAA.IF_RD = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_RGWLK = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_KEYSTN = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_LINSFT = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_SUBPATCH = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_GEODIST = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_RGFO = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_BEAM_COMP = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_RGRES = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_AZRES = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_RSS = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_AZSS = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_RSR = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_AZSR = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_RFFT_SAMP = str2double(fread(fid,7,'uint8=>char')');
CMETAA.IF_AZFFT_SAMP = str2double(fread(fid,7,'uint8=>char')');
CMETAA.IF_RFFT_TOT = str2double(fread(fid,7,'uint8=>char')');
CMETAA.IF_AZFFT_TOT = str2double(fread(fid,7,'uint8=>char')');
CMETAA.IF_SUBP_ROW = str2double(fread(fid,6,'uint8=>char')');
CMETAA.IF_SUBP_COL = str2double(fread(fid,6,'uint8=>char')');
CMETAA.IF_SUB_RG = str2double(fread(fid,4,'uint8=>char')');
CMETAA.IF_SUB_AZ = str2double(fread(fid,4,'uint8=>char')');
CMETAA.IF_RFFTS = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_AFFTS = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_RANGE_DATA = strtrim(fread(fid,7,'uint8=>char')');
CMETAA.IF_INCPH = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.IF_SR_NAME1 = strtrim(fread(fid,8,'uint8=>char')');
CMETAA.IF_SR_AMOUNT1 = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_SR_NAME2 = strtrim(fread(fid,8,'uint8=>char')');
CMETAA.IF_SR_AMOUNT2 = str2double(fread(fid,8,'uint8=>char')');
CMETAA.IF_SR_NAME3 = strtrim(fread(fid,8,'uint8=>char')');
CMETAA.IF_SR_AMOUNT3 = str2double(fread(fid,8,'uint8=>char')');

for lp = 1:3
    CMETAA.AF_TYPE{lp} = strtrim(fread(fid,5,'uint8=>char')');
end

CMETAA.POL_TR = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_RE = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_REFERENCE = strtrim(fread(fid,40,'uint8=>char')');
CMETAA.POL = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_REG = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_ISO_1 = str2double(fread(fid,5,'uint8=>char')');
CMETAA.POL_BAL = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_BAL_MAG = str2double(fread(fid,8,'uint8=>char')');
CMETAA.POL_BAL_PHS = str2double(fread(fid,8,'uint8=>char')');
CMETAA.POL_HCOMP = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_HCOMP_BASIS = strtrim(fread(fid,10,'uint8=>char')');

for lp = 1:3
    CMETAA.POL_HCOMP_COEF(lp) = str2double(fread(fid,9,'uint8=>char')');
end

CMETAA.POL_AFCOMP = strtrim(fread(fid,1,'uint8=>char')');
CMETAA.POL_SPARE_A = strtrim(fread(fid,15,'uint8=>char')');
CMETAA.POL_SPARE_N = strtrim(fread(fid,9,'uint8=>char')');

CMETAA.T_UTC_YYYYMMMDD = strtrim(fread(fid,9,'uint8=>char')');
CMETAA.T_HHMMSSUTC = strtrim(fread(fid,6,'uint8=>char')');
CMETAA.T_HHMMSSLOCAL = strtrim(fread(fid,6,'uint8=>char')');

CMETAA.CG_SRAC = str2double(fread(fid,11,'uint8=>char')');
CMETAA.CG_SLANT_CONFIDENCE = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CG_CROSS = str2double(fread(fid,11,'uint8=>char')');
CMETAA.CG_CROSS_CONFIDENCE = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CG_CAAC = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_CONE_CONFIDENCE = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_GPSAC = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_GPSAC_CONFIDENCE = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_SQUINT = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_GAAC = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CG_GAAC_CONFIDENCE = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_INCIDENT = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CG_SLOPE = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CG_TILT = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_LD = fread(fid,1,'uint8=>char')';
CMETAA.CG_NORTH = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_NORTH_CONFIDENCE = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_EAST = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_RLOS = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_LOS_CONFIDENCE = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_LAYOVER = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_SHADOW = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_OPM = str2double(fread(fid,7,'uint8=>char')');
CMETAA.CG_MODEL = strtrim(fread(fid,5,'uint8=>char')');

CMETAA.CG_AMPT_X = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_AMPT_Y = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_AMPT_Z = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_AP_CONF_XY = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_AP_CONF_Z = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_APCEN_X = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_APCEN_Y = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_APCEN_Z = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_APER_CONF_XY = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_APER_CONF_Z = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_FPNUV_X = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_FPNUV_Y = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_FPNUV_Z = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_IDPNUVX = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_IDPNUVY = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_IDPNUVZ = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_SCECN_X = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_SCECN_Y = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_SCECN_Z = str2double(fread(fid,13,'uint8=>char')');
CMETAA.CG_SC_CONF_XY = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_SC_CONF_Z = str2double(fread(fid,6,'uint8=>char')');
CMETAA.CG_SWWD = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_SNVEL_X = str2double(fread(fid,10,'uint8=>char')');
CMETAA.CG_SNVEL_Y = str2double(fread(fid,10,'uint8=>char')');
CMETAA.CG_SNVEL_Z = str2double(fread(fid,10,'uint8=>char')');
CMETAA.CG_SNACC_X = str2double(fread(fid,10,'uint8=>char')');
CMETAA.CG_SNACC_Y = str2double(fread(fid,10,'uint8=>char')');
CMETAA.CG_SNACC_Z = str2double(fread(fid,10,'uint8=>char')');
CMETAA.CG_SNATT_ROLL = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_SNATT_PITCH = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_SNATT_YAW = str2double(fread(fid,8,'uint8=>char')');
CMETAA.CG_GTP_X = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_GTP_Y = str2double(fread(fid,9,'uint8=>char')');
CMETAA.CG_GTP_Z = str2double(fread(fid,9,'uint8=>char')');

CMETAA.CG_MAP_TYPE = strtrim(fread(fid,4,'uint8=>char')');

if strcmp(CMETAA.CG_MAP_TYPE,'GEOD')
    CMETAA.CG_PATCH_LATCEN = str2double(fread(fid,11,'uint8=>char')');
    CMETAA.CG_PATCH_LNGCEN = str2double(fread(fid,12,'uint8=>char')');
    CMETAA.CG_PATCH_LTCORUL = str2double(fread(fid,11,'uint8=>char')');
    CMETAA.CG_PATCH_LGCORUL = str2double(fread(fid,12,'uint8=>char')');
    CMETAA.CG_PATCH_LTCORUR = str2double(fread(fid,11,'uint8=>char')');
    CMETAA.CG_PATCH_LGCORUR = str2double(fread(fid,12,'uint8=>char')');
    CMETAA.CG_PATCH_LTCORLR = str2double(fread(fid,11,'uint8=>char')');
    CMETAA.CG_PATCH_LGCORLR = str2double(fread(fid,12,'uint8=>char')');
    CMETAA.CG_PATCH_LTCORLL = str2double(fread(fid,11,'uint8=>char')');
    CMETAA.CG_PATCH_LNGCOLL = str2double(fread(fid,12,'uint8=>char')');
    CMETAA.CG_PATCH_LAT_CONFIDENCE = str2double(fread(fid,9,'uint8=>char')');
    CMETAA.CG_PATCH_LONG_CONFIDENCE = str2double(fread(fid,9,'uint8=>char')');
elseif strcmp(CMETAA.CG_MAP_TYPE,'MGRS')
    CMETAA.CG_MGRS_CENT = strtrim(fread(fid,23,'uint8=>char')');
    CMETAA.CG_MGRSCORUL = strtrim(fread(fid,23,'uint8=>char')');
    CMETAA.CG_MGRSCORUR = strtrim(fread(fid,23,'uint8=>char')');
    CMETAA.CG_MGRSCORLR = strtrim(fread(fid,23,'uint8=>char')');
    CMETAA.CG_MGRCORLL = strtrim(fread(fid,23,'uint8=>char')');
    CMETAA.CG_MGRS_CONFIDENCE = str2double(fread(fid,7,'uint8=>char')');
    CMETAA.CG_MGRS_PAD = strtrim(fread(fid,11,'uint8=>char')');
elseif strcmp(CMETAA.CG_MAP_TYPE,'NA')
    CMETAA.CG_MAP_TYPE_BLANK = strtrim(fread(fid,133,'uint8=>char')');
end
CMETAA.CG_SPARE_A = strtrim(fread(fid,144,'uint8=>char')');
CMETAA.CA_CALPA = str2double(fread(fid,7,'uint8=>char')');
CMETAA.WF_SRTFR = str2double(fread(fid,14,'uint8=>char')');
CMETAA.WF_ENDFR = str2double(fread(fid,14,'uint8=>char')');
CMETAA.WF_CHRPRT = str2double(fread(fid,10,'uint8=>char')');
CMETAA.WF_WIDTH = str2double(fread(fid,9,'uint8=>char')');
CMETAA.WF_CENFRQ = str2double(fread(fid,13,'uint8=>char')');
CMETAA.WF_BW = str2double(fread(fid,13,'uint8=>char')');
CMETAA.WF_PRF = str2double(fread(fid,7,'uint8=>char')');
CMETAA.WF_PRI = str2double(fread(fid,9,'uint8=>char')');
CMETAA.WF_CDP = str2double(fread(fid,7,'uint8=>char')');
CMETAA.WF_NUMBER_OF_PULSES = str2double(fread(fid,9,'uint8=>char')');

CMETAA.VPH_COND = strtrim(fread(fid,1,'uint8=>char')');

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function f = goldstein_price ( x )

%*****************************************************************************80
%
%% GOLDSTEIN_PRICE evaluates the Goldstein-Price polynomial.
%
%  Discussion:
%
%    The minimizer is
%
%      X* = [ 0.0, -1.0 ]
%      F(X*) = 3.0
%
%    Suggested starting point:
%
%      X init = [ -0.5, 0.25 ] (easy convergence)
%      X init = [ -4.0, 5.00 ] (harder convergence)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, real X(2), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  if ( length ( x ) ~= 2 )
    error ( 'Error: function expects a two dimensional input\n' );
  end

  a = x(1) + x(2) + 1.0;

  b = 19.0 - 14.0 * x(1) + 3.0 * x(1) * x(1) - 14.0 * x(2) ...
    + 6.0 * x(1) * x(2) + 3.0 * x(2) * x(2);

  c = 2.0 * x(1) - 3.0 * x(2);

  d = 18.0 - 32.0 * x(1) + 12.0 * x(1) * x(1) + 48.0 * x(2) ...
    - 36.0 * x(1) * x(2) + 27.0 * x(2) * x(2);

  f = ( 1.0 + a * a * b ) * ( 30.0 + c * c * d );

  return
end