text:
%  Save back to the back to the original image a portion of slices that
%  were loaded by "load_untouch_nii". You can process those slices matrix
%  in any way, as long as their dimension is not changed.
%
%  Usage: save_untouch_slice(slice, filename, ...
%		slice_idx, [img_idx], [dim5_idx], [dim6_idx], [dim7_idx])
%
%  Example:
%	nii = load_nii('avg152T1_LR_nifti.nii');
%	save_nii(nii, 'test.nii');
%	view_nii(nii);
%	nii = load_untouch_nii('test.nii','','','','','',[40 51:53]);
%	nii.img = ones(91,109,4)*122;
%	save_untouch_slice(nii.img, 'test.nii', [40 51:52]);
%	nii = load_nii('test.nii');
%	view_nii(nii);
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function save_untouch_slice(slice, filename, slice_idx, img_idx, dim5_idx, dim6_idx, dim7_idx)

   if ~exist('filename','var') | ~exist('slice_idx','var')
      error('Usage: save_untouch_slice(slice, filename, slice_idx, [img_idx], [dim5_idx], [dim6_idx], [dim7_idx])');
   end

   if ~exist('img_idx','var') | isempty(img_idx)
      img_idx = [];
   end

   if ~exist('dim5_idx','var') | isempty(dim5_idx)
      dim5_idx = [];
   end

   if ~exist('dim6_idx','var') | isempty(dim6_idx)
      dim6_idx = [];
   end

   if ~exist('dim7_idx','var') | isempty(dim7_idx)
      dim7_idx = [];
   end

   %  Read the dataset header
   %
   [nii.hdr,nii.filetype,nii.fileprefix,nii.machine] = load_nii_hdr(filename);

   if nii.filetype == 0
      nii.hdr = load_untouch0_nii_hdr(nii.fileprefix,nii.machine);
   else
      nii.hdr = load_untouch_nii_hdr(nii.fileprefix,nii.machine,nii.filetype);
   end

   %  Save the dataset body
   %
   save_untouch_slice_img(slice,filename,nii.hdr,nii.filetype,nii.fileprefix, ...
	nii.machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx);

   return					% save_untouch_slice

%--------------------------------------------------------------------------
function save_untouch_slice_img(slice,filename,hdr,filetype,fileprefix,machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx)

   if ~exist('hdr','var') | ~exist('filetype','var') | ~exist('fileprefix','var') | ~exist('machine','var')
      error('Usage: save_untouch_slice_img(slice,filename,hdr,filetype,fileprefix,machine,slice_idx,[img_idx],[dim5_idx],[dim6_idx],[dim7_idx]);');
   end

   if ~exist('slice_idx','var') | isempty(slice_idx) | hdr.dime.dim(4)<1
      slice_idx = [];
   end

   if ~exist('img_idx','var') | isempty(img_idx) | hdr.dime.dim(5)<1
      img_idx = [];
   end

   if ~exist('dim5_idx','var') | isempty(dim5_idx) | hdr.dime.dim(6)<1
      dim5_idx = [];
   end

   if ~exist('dim6_idx','var') | isempty(dim6_idx) | hdr.dime.dim(7)<1
      dim6_idx = [];
   end

   if ~exist('dim7_idx','var') | isempty(dim7_idx) | hdr.dime.dim(8)<1
      dim7_idx = [];
   end

   %  check img_idx
   %
   if ~isempty(img_idx) & ~isnumeric(img_idx)
      error('"img_idx" should be a numerical array.');
   end

   if length(unique(img_idx)) ~= length(img_idx)
      error('Duplicate image index in "img_idx"');
   end

   if ~isempty(img_idx) & (min(img_idx) < 1 | max(img_idx) > hdr.dime.dim(5))
      max_range = hdr.dime.dim(5);

      if max_range == 1
         error(['"img_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"img_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim5_idx
   %
   if ~isempty(dim5_idx) & ~isnumeric(dim5_idx)
      error('"dim5_idx" should be a numerical array.');
   end

   if length(unique(dim5_idx)) ~= length(dim5_idx)
      error('Duplicate index in "dim5_idx"');
   end

   if ~isempty(dim5_idx) & (min(dim5_idx) < 1 | max(dim5_idx) > hdr.dime.dim(6))
      max_range = hdr.dime.dim(6);

      if max_range == 1
         error(['"dim5_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim5_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim6_idx
   %
   if ~isempty(dim6_idx) & ~isnumeric(dim6_idx)
      error('"dim6_idx" should be a numerical array.');
   end

   if length(unique(dim6_idx)) ~= length(dim6_idx)
      error('Duplicate index in "dim6_idx"');
   end

   if ~isempty(dim6_idx) & (min(dim6_idx) < 1 | max(dim6_idx) > hdr.dime.dim(7))
      max_range = hdr.dime.dim(7);

      if max_range == 1
         error(['"dim6_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim6_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim7_idx
   %
   if ~isempty(dim7_idx) & ~isnumeric(dim7_idx)
      error('"dim7_idx" should be a numerical array.');
   end

   if length(unique(dim7_idx)) ~= length(dim7_idx)
      error('Duplicate index in "dim7_idx"');
   end

   if ~isempty(dim7_idx) & (min(dim7_idx) < 1 | max(dim7_idx) > hdr.dime.dim(8))
      max_range = hdr.dime.dim(8);

      if max_range == 1
         error(['"dim7_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim7_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check slice_idx
   %
   if ~isempty(slice_idx) & ~isnumeric(slice_idx)
      error('"slice_idx" should be a numerical array.');
   end

   if length(unique(slice_idx)) ~= length(slice_idx)
      error('Duplicate index in "slice_idx"');
   end

   if ~isempty(slice_idx) & (min(slice_idx) < 1 | max(slice_idx) > hdr.dime.dim(4))
      max_range = hdr.dime.dim(4);

      if max_range == 1
         error(['"slice_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"slice_idx" should be an integer within the range of [' range '].']);
      end
   end

   write_image(slice,filename,hdr,filetype,fileprefix,machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx);

   return					% save_untouch_slice_img

%---------------------------------------------------------------------
function write_image(slice,fn,hdr,filetype,fileprefix,machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx)

   fid = fopen(fn,'r+',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   end

   %  Set bitpix according to datatype
   %
   %  /*Acceptable values for datatype are*/
   %
   %     0 None                     (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN
   %     1 Binary                         (ubit1, bitpix=1) % DT_BINARY
   %     2 Unsigned char         (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8
   %     4 Signed short                  (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16
   %     8 Signed integer                (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32
   %    16 Floating point    (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32
   %    32 Complex, 2 float32      (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
   %    64 Double precision  (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64
   %   128 uint8 RGB                 (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24
   %   256 Signed char            (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8
   %   511 Single RGB              (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
   %   512 Unsigned short               (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16
   %   768 Unsigned integer             (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32
   %  1024 Signed long long              (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
   %  1280 Unsigned long long           (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64
   %  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
   %  1792 Complex128, 2 float64  (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
   %  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
   %
   switch hdr.dime.datatype
   case   2,
      hdr.dime.bitpix = 8;  precision = 'uint8';
   case   4,
      hdr.dime.bitpix = 16; precision = 'int16';
   case   8,
      hdr.dime.bitpix = 32; precision = 'int32';
   case  16,
      hdr.dime.bitpix = 32; precision = 'float32';
   case  64,
      hdr.dime.bitpix = 64; precision = 'float64';
   case 128,
      hdr.dime.bitpix = 24; precision = 'uint8';
   case 256
      hdr.dime.bitpix = 8;  precision = 'int8';
   case 511
      hdr.dime.bitpix = 96; precision = 'float32';
   case 512
      hdr.dime.bitpix = 16; precision = 'uint16';
   case 768
      hdr.dime.bitpix = 32; precision = 'uint32';
   case 1024
      hdr.dime.bitpix = 64; precision = 'int64';
   case 1280
      hdr.dime.bitpix = 64; precision = 'uint64';
   otherwise
      error('This datatype is not supported');
   end

   hdr.dime.dim(find(hdr.dime.dim < 1)) = 1;

   %  move pointer to the start of image block
   %
   switch filetype
   case {0, 1}
      fseek(fid, 0, 'bof');
   case 2
      fseek(fid, hdr.dime.vox_offset, 'bof');
   end

   if hdr.dime.datatype == 1 | isequal(hdr.dime.dim(4:8),ones(1,5)) | ...
	(isempty(img_idx) & isempty(dim5_idx) & isempty(dim6_idx) & isempty(dim7_idx) & isempty(slice_idx))

      msg = [char(10) char(10) '   "save_untouch_slice" is used to save back to the original image a' char(10)];
      msg = [msg '   portion of slices that were loaded by "load_untouch_nii". You can' char(10)];
      msg = [msg '   process those slices matrix in any way, as long as their dimension' char(10)];
      msg = [msg '   is not changed.'];
      error(msg);
   else

      d1 = hdr.dime.dim(2);
      d2 = hdr.dime.dim(3);
      d3 = hdr.dime.dim(4);
      d4 = hdr.dime.dim(5);
      d5 = hdr.dime.dim(6);
      d6 = hdr.dime.dim(7);
      d7 = hdr.dime.dim(8);

      if isempty(slice_idx)
         slice_idx = 1:d3;
      end

      if isempty(img_idx)
         img_idx = 1:d4;
      end

      if isempty(dim5_idx)
         dim5_idx = 1:d5;
      end

      if isempty(dim6_idx)
         dim6_idx = 1:d6;
      end

      if isempty(dim7_idx)
         dim7_idx = 1:d7;
      end

      %ROMAN: begin
      roman = 1;
      if(roman)

         %  compute size of one slice
         %
         img_siz = prod(hdr.dime.dim(2:3));

         %  For complex float32 or complex float64, voxel values
         %  include [real, imag]
         %
         if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
            img_siz = img_siz * 2;
         end

         %MPH: For RGB24, voxel values include 3 separate color planes
         %
         if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
            img_siz = img_siz * 3;
         end

      end; %if(roman)
      % ROMAN: end

      for i7=1:length(dim7_idx)
         for i6=1:length(dim6_idx)
            for i5=1:length(dim5_idx)
               for t=1:length(img_idx)
               for s=1:length(slice_idx)

                  %  Position is seeked in bytes. To convert dimension size
                  %  to byte storage size, hdr.dime.bitpix/8 will be
                  %  applied.
                  %
                  pos = sub2ind([d1 d2 d3 d4 d5 d6 d7], 1, 1, slice_idx(s), ...
			                    img_idx(t), dim5_idx(i5),dim6_idx(i6),dim7_idx(i7)) -1;
                  pos = pos * hdr.dime.bitpix/8;

                  % ROMAN: begin
                  if(roman)
                      % do nothing
                  else
                     img_siz = prod(hdr.dime.dim(2:3));

                     %  For complex float32 or complex float64, voxel values
                     %  include [real, imag]
                     %
                     if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
                        img_siz = img_siz * 2;
                     end

                     %MPH: For RGB24, voxel values include 3 separate color planes
                     %
                     if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
                        img_siz = img_siz * 3;
                     end
                  end; % if (roman)
                  % ROMAN: end

                  if filetype == 2
                     fseek(fid, pos + hdr.dime.vox_offset, 'bof');
                  else
                     fseek(fid, pos, 'bof');
                  end

                  %  For each frame, fwrite will write precision of value
                  %  in img_siz times
                  %
                  fwrite(fid, slice(:,:,s,t,i5,i6,i7), sprintf('*%s',precision));

               end
               end
            end
         end
      end
   end

   fclose(fid);

   return						% write_image


text:
function B = rot3d90(A, X, K)
% ROT3D90 - rotate 3D array 90 degrees around a specific axis
%
%   B = rot3D90(A, X) is the 90 degree counterclockwise rotation of the 3D
%   array A along a specific axis. X denotes the plane of rotation,
%   according to the table below. For instance, when X is 1 the plane of
%   rotation is formed by the 2nd and 3rd dimension.
%
%   rot3d90(A, X, K) is the K*90 degree rotation of A, K = +-1,+-2,...
%
%   If A has a dimension of [N, M, L], the dimensions of the rotation
%   depend on X. For value of X, the subindex in the X-th dimension will be
%   unaffected. Therefore, specific vectors will be unaffected.
%   Note that rot3d90(A,3) and rot90(A) produce the same result.
%
%      X    Rotation plane    Dimension of B   Unaffected vectors
%     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%      1    dims 2 & 3        [N, L, M]        N-by-1 (row)
%      2    dims 1 & 3        [L, M, N]        1-by-N (column)
%      3    dims 1 & 2        [M, N, L]        1-by-1-by-N (plane)
%
%   Example:
%      A = cat(3,[1 2 ; 3 4],[5 6 ; 7 8])
%      B = rot3d90(A, 2) % rotation in plane formed by 1st and 3rd dimension
%      % B(:,:,1) = 5  6
%      %            1  2
%      % B(:,:,2) = 7  8
%      %            3  4
%
%   See also ROT90, FLIP

% version 3.1, jan 2018
% (c) Jos van der Geest
% email: samelinoa@gmail.com
% FEX: http://www.mathworks.nl/matlabcentral/fileexchange/authors/10584

% 1.0 (jan 2018) created, using rot90 and permute
% 2.0 (jan 2018) replaced calls to rot90 by direct permuting and flipping
% 3.0 (jan 2018) added K, specifiying the number of rotations
% 3.1 (jan 2018) fixed an error for 2D arrays inputs

narginchk(2, 3)

if nargin == 2
    K = 1 ;
else
    if ~isscalar(K) || fix(K) ~= K
        error('rot3d90:KInvalid', 'K should be a scalar integer.');
    end
    K = mod(K, 4) ; % after 4 rotations we are back to A
end

B = A ;
for j = 1:K
    switch X
        case 1
            % rotation in the plane formed by the 2nd and 3rd dimension
            % N-by-1-by-1 (column) vectors stay the same
            % for all elements, the subindex in the 1st dim does not change
            B = permute(B, [1 3 2 4:ndims(B)]);
            B = flip(B,2);
        case 2
            % rotation in the plane formed by 1st and 3rd dimension
            % 1-by-N-by-1 (row) vectors stay the same
            % for all elements, the subindex in the 2nd dim does not change
            B = permute(B, [2 3 1 4:ndims(B)]) ;
            B = flip(B, 2);
            B = permute(B, [2 1 3:ndims(B)]);
        case 3
            % rotate in the plane formed by 1st and 2nd dimension
            % 1-by-1-by-N vectors stay the same
            % for all elements, the subindex in the 3rd dim does not change
            B = flip(B, 2);
            B = permute(B, [2 1 3:ndims(B)]);
        otherwise
            error('rot3d90:XInvalid', 'Invalid rotation axis. X should be 1, 2, or 3.') ;
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function C = gsp_jtv_mat2vec(C)
%GSP_JTV_MAT2VEC vector to matrix transform
%   Usage:  C  = gsp_jtv_mat2vec(C);
%
%   Input parameters:
%       C       : Data
%
%   Ouput parameter
%       C       : Data
%
%   Reshape the data from the matrix form to the vector form

[N,T,Nf,Ns] = size(C);

C = reshape(C,N*T*Nf,Ns);

end

repo: gspbox


text:
function [ know, x ] = p08_sol ( m, know )

%*****************************************************************************80
%
%% P08_SOL returns the solution for problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution, if known.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = 0.0;
  else
    know = 0;
    x(1:m,1) = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function status = test_hrtf_extrapolation(modus,hrtf_set)
%TEST_HRTF_EXTRAPOLATION tests the HRTF extrapolation functions
%
%   Usage: status = test_hrtf_extrapolation(modus,hrtf_set)
%
%   Input parameters:
%       modus    - 0: numerical
%                  1: visual
%       hrtf_set - HRTFs for testing extrapolation. Can be
%                    'QU_KEMAR'
%                    'FABIAN_3D'
%
%   Output parameters:
%       status - true or false

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

% TODO: add mode to save data as reference data
status = false;

%% ===== Checking of input  parameters ===================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);

%% ===== Configuration ===================================================
conf = SFS_config;
conf.showprogress = true;
if strcmp('QU_KEMAR',hrtf_set)
    % QU KEMAR horizontal HRTFs
    disp('Extrapolate QU KEMAR anechoic 3m ...');
    % check if HRTF data set is available, download otherwise
    basepath = get_sfs_path();
    hrtf_file = [basepath '/data/HRTFs/QU_KEMAR_anechoic_3m.sofa'];
    if ~exist(hrtf_file,'file')
        url = ['https://raw.githubusercontent.com/sfstoolbox/data/master/', ...
               'HRTFs/QU_KEMAR_anechoic_3m.sofa'];
        download_file(url,hrtf_file);
    end
    % load HRTF data set
    sofa = SOFAload(hrtf_file);
    x0 = SOFAcalculateAPV(sofa);
    % do the extrapolation
    sofa_pw = extrapolate_farfield_hrtfset(sofa,conf);
    % plot the original HRTF data set
    if modus
        figure;
        imagesc(x0(:,1),1:size(sofa.Data.IR,3),squeeze(sofa.Data.IR(:,1,:))');
        title('QU KEMAR anechoic 3m');
        xlabel('phi / deg');
        % plot the interplated HRTF data set
        figure;
        imagesc(x0(:,1),1:size(sofa_pw.Data.IR,3),squeeze(sofa_pw.Data.IR(:,1,:))');
        title('QU KEMAR anechoic extrapolated');
        xlabel('phi / deg');
    end
    % ILD of both HRTF sets
    ild1 = db(rms(squeeze(sofa.Data.IR(:,1,:))')) - ...
           db(rms(squeeze(sofa.Data.IR(:,2,:))'));
    ild2 = db(rms(squeeze(sofa_pw.Data.IR(:,1,:))')) - ...
           db(rms(squeeze(sofa_pw.Data.IR(:,2,:))'));
    if modus
        figure;
        plot(x0(:,1),ild1,'-b', ...
             x0(:,1),ild2,'-r');
        legend('original','extrapolated');
        title('Interaural Level Differences');
        xlabel('phi / deg');
        ylabel('Amplitude difference / dB');
    end
elseif strcmp('FABIAN_3D',hrtf_set)
    error(['%s: the FABIAN 3D data set is not publicly available at the ' ...
    'moment.'],upper(mfilename));
    % SEACEN FABIAN 3D HRTFs
    disp('Extrapolate SEACEN FABIAN anechoic ...');
    % extrapolation settings
    conf.dimension = '3D';
    conf.usetapwin = false;
    conf.secondary_sources.geometry = 'sphere';
    hrtf_file = 'FABIAN_3d_anechoic.sofa';
    % do the extrapolation
    irs_pw = extrapolate_farfield_hrtfset(hrtf_file,conf);
    %save FABIAN_3D_extrapolated.sofa irs_pw
    % plot the original HRTF data set
    if modus
        figure;
        imagesc(x0(:,1),1:size(sofa.Data.IR,3),squeeze(sofa.Data.IR(:,1,:))');
        title('SEACEN FABIAN anechoic 1.7m');
        xlabel('phi / deg');
        % plot the interplated HRTF data set
        figure;
        imagesc(x0(:,1),1:size(sofa_pw.Data.IR,3),squeeze(sofa.Data_pw.IR(:,1,:))');
        title('SEACEN FABIAN anechoic extrapolated');
        xlabel('phi / deg');
    end
    % ILD of both HRTF sets
    % ILD of both HRTF sets
    ild1 = db(rms(squeeze(sofa.Data.IR(:,1,:))')) - ...
           db(rms(squeeze(sofa.Data.IR(:,2,:))'));
    ild2 = db(rms(squeeze(sofa_pw.Data.IR(:,1,:))')) - ...
           db(rms(squeeze(sofa_pw.Data.IR(:,2,:))'));
    if modus
        figure;
        plot(x0(:,1),ild1,'-b', ...
             x0(:,1),ild2,'-r');
        legend('original','extrapolated');
        title('Interaural Level Differences');
        xlabel('phi / deg');
        ylabel('Amplitude difference / dB');
    end
end

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
function pout = polydern(p,diffvar)
% polydern: partial derivatives of a n-dimensional polynomial model created by polyfitn.
% usage: pout = polydern(p,diffvar)
%
% arguments: (input)
%  p - a structure containing a polynomial regression model.  See polyfitn
%      for details.
%
%  diffvar - integer, which variable to compute the derivative with respect
%      to. If diffvar is a character string, then it refers to an explicit
%      variable in the model that is found in the VarNames field. If so,
%      the variable name must be an exact match for one of the variables
%      listed in p.VarNames.
%
% arguments: (output)
% pout - a structure containing a polynomial regression model which
%           describes the partial derivative with respect to variable #i.
%
% Note:
% Polydern modifies the "ModelTerms" and "Coefficients" fields in creating
% the output polynomial. ParameterVar, parameterStd are also updated to
% have the proper shape, as well as scaling them to reflect the derivative.
%
% The other fields (R2, RMSE, VarNames) are copied without change, even
% though the error measures are no longer meaningful for the derivative
% polynomial.
%
% This function uses polynomial structure variables created by polyfitn,
% which is part of the PolyfitnTools package by John D'Errico:
%     http://www.mathworks.com/matlabcentral/fileexchange/10065
%
% See also: polyder, polyfitn, polyvaln
%
% Polydern author: Jason Goodman
% Error checks, ParameterStd/Var updates and allowance for the variable
% name itself instead of an index added by John D'Errico

% check the number of input args
if nargin ~= 2
  error('POLYDERN:argumentcount','Exactly two arguments are required')
end

% is p a struct, created by polyfitn?
if ~isstruct(p) || ~isfield(p,'ModelTerms') || ~isfield(p,'Coefficients')
  error('POLYDERN:invalidp','p must be a struct as created by polyfitn')
end

% check diffvar too. scalar, numeric, integer, real, in the proper range,
% or it must be character, an exact match for one of the variables listed
% in p.VarNames.
if ischar(diffvar)
  ind = ismember(p.VarNames,diffvar);
  if sum(ind) == 1
    % we have a hit
    diffvar = find(ind);
  else
    error('POLYDERN:invaliddiffvar','diffvar must be a scalar numeric real integer value, or a valid variable name')
  end
elseif ~isnumeric(diffvar) || (numel(diffvar) ~= 1) || (round(diffvar) ~= diffvar) || ~isreal(diffvar)
  % diffvar must be scalar, numeric, etc. Also verify that it is real to
  % be complete
  error('POLYDERN:invaliddiffvar','diffvar must be a scalar numeric real integer value, or a valid variable name')
elseif (diffvar < 1) || (diffvar > size(p.ModelTerms,2))
  % test that diffvar is at least 1, and is not larger than the
  % number of variables in the nmodel
  error('POLYDERN:invaliddiffvar', ...
    ['diffvar (=',num2str(diffvar), ...
    ') cannot represent the variable index to be differentiated, as it is out of range for this model'])
end

pout = p;
pout.ModelTerms = [];
pout.Coefficients = [];

remainingterms = true(1,size(p.ModelTerms,1));
jout = 1;
for jin=1:length(p.Coefficients)
    if (p.ModelTerms(jin,diffvar) ~= 0)
        pout.Coefficients(jout) = p.Coefficients(jin).*p.ModelTerms(jin,diffvar);
        pout.ModelTerms(jout,:) = p.ModelTerms(jin,:);
        pout.ModelTerms(jout,diffvar) = pout.ModelTerms(jout,diffvar)-1;
        jout = jout + 1;
    else
        % this term got dropped from the model
        remainingterms(jin) = false;
    end
end

if (jout == 1)   % Polynomial has no terms in this variable: deriv is zero
    [ncoeff, nvars] = size(p.ModelTerms);
    pout.ModelTerms(1,:) = zeros(1,nvars);
    pout.Coefficients = 0;

    pout.ParameterVar = 0;
    pout.ParameterStd = 0;
else
    % at least some terms remain, so update the parameter variances and
    % standard deviations. First, drop those terms that died.
    pout.ParameterVar = reshape(p.ParameterVar(remainingterms),1,[]);
    pout.ParameterStd = reshape(p.ParameterStd(remainingterms),1,[]);

    % scale them appropriately. Thus if std(A) is sigma, then std(k*A)
    % is k*sigma, but multiply the corresponding variances by k.^2.
    scale = p.ModelTerms(remainingterms,diffvar).';
    pout.ParameterStd = pout.ParameterStd.*scale;
    pout.ParameterVar = pout.ParameterVar.*scale.^2;
end

author: Sable

repo: mcbench-benchmarks


text:
function lobatto_polynomial_test ( )

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

  lobatto_polynomial_value_test ( );
  lobatto_polynomial_derivative_test ( );
  lobatto_polynomial_plot_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LOBATTO_POLYNOMIAL_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%TEST_bezierFunction.m
%
% This script tests the use of a bezier curve to represent a vector
% function, rather than an arbitrary space curve
%

order = 4;

p = rand(1,order+1);

tSpan = [1,4];
t = linspace(tSpan(1),tSpan(2),100);
pGrid = linspace(tSpan(1),tSpan(2),order+1);

tic
x = bezierCurve(p,t,tSpan);
toc

figure(1); clf;
plot(t,x); hold on;
plot(pGrid,p,'x');

author: MatthewPeterKelly

repo: dscTutorials


text:
bn = '/biac3/wandell4/data/reading_longitude/dti_y4/zs070717/raw/dti_g13_b800';
bn = '/biac3/wandell4/data/reading_longitude/dti_y4/zs070717/raw/dti_g87_b800';

dwRaw = niftiRead([bn '.nii.gz']);
sz = size(dwRaw.data);
d = single(dwRaw.data);
clear dwRaw;
bvecs = dlmread([bn '.bvecs']);
bvals = dlmread([bn '.bvals']);

% Apply eddy/motion correction
load([bn '_ecXform.mat']);

[X,Y,Z] = ndgrid([1:sz(1)],[1:sz(2)],[1:sz(3)]);
x = [X(:)'; Y(:)'; Z(:)']; clear X Y Z;
% Run the first one separate to cache the voxel list (x).
d(:,:,:,1) = reshape(mrAnatFastInterp3(d(:,:,:,1), x, [xform(1).ecParams xform(1).phaseDir]),sz(1:3));
for(ii=2:sz(4))
    if(mod(ii,30)==0), fprintf('Processing vol %d of %d...\n',ii,sz(4)); end
    % Save a few cycles by not resending the voxel list (x) on each iteration.
    d(:,:,:,1) = reshape(mrAnatFastInterp3(d(:,:,:,ii), [], [xform(ii).ecParams xform(ii).phaseDir]),sz(1:3));
end
d = double(d);

bv = dtiRawReorientBvecs(bvecs,xform,false);
tau = 40;
q = [bv.*sqrt(repmat(bvals./tau,3,1))]';
X = [ones(size(q,1),1) -tau.*q(:,1).^2 -tau.*q(:,2).^2 -tau.*q(:,3).^2 -2*tau.*q(:,1).*q(:,2) -2*tau.*q(:,1).*q(:,3) -2*tau.*q(:,2).*q(:,3)];

[x,y,z]=meshgrid([5:10],[5:10],[30:40]);
noiseRegion = sub2ind(sz(1:3),x,y,z);
for(ii=1:sz(4))
   tmp = d(:,:,:,ii);
   n(ii,:) = double(tmp(noiseRegion(:)));
end
sd = std(n(:));
signalSigma = 1.5267*sd;

tic;[dt,pdd] = dtiFitTensor(d,X); toc
makeMontage3(abs(pdd));

[fa,md] = dtiComputeFA(dt(:,:,:,2:7));
md(md>5)=5; md(md<0)=0;
showMontage(md);
fa(fa<0)=0; fa(fa>1)=1;
showMontage(fa);

sl = squeeze(md(6,:,:));
figure;imagesc(sl);axis image; colormap gray; colorbar
roi = roipoly;
t = 25;
fprintf('Self-diffusion @ %0.1f C = %0.2f, measured diffusivity = %0.2f.\n',t,dtiSelfDiffusionOfWater(t),mean(sl(roi(:))));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
addpath /n/fs/vision/www/pvt/

FrameHeight = 630;
TemplateFileName = 'template.html';
TemplateString = file2string(TemplateFileName);
amazonHead = '<?xml version="1.0"?><HTMLQuestion xmlns="http://mechanicalturk.amazonaws.com/AWSMechanicalTurkDataSchemas/2011-11-11/HTMLQuestion.xsd"><HTMLContent><![CDATA[<!DOCTYPE html><title>HIT</title><meta http-equiv="Content-Type" content="text/html; charset=UTF-8"/><form name="mturk_form" method="post" id="mturk_form" action="https://www.mturk.com/mturk/externalSubmit"><input type="hidden" value="" name="assignmentId" id="assignmentId" />';
amazonTail = ['</form>]]></HTMLContent><FrameHeight>' num2str(FrameHeight) '</FrameHeight></HTMLQuestion>'];
QuestionTemplate = [amazonHead TemplateString amazonTail];

% put your stuff here
for questionID=1:1000
    Question = QuestionTemplate;
    Question = strrep(Question,'${questionID}',num2str(questionID));
    Question = strrep(Question,'${data}','yinda fillin');
    Question = strrep(Question,'${examples}','yinda fillin');
    Question = strrep(Question,'${definition}','yinda fillin');
    Question = strrep(Question,'${question}','yinda fillin');
    Question = strrep(Question,'${listname}','yinda fillin');
end

author: jianxiongxiao

repo: ProfXkit


text:
clc
clear all
seed = 601;
randn('state',seed); rand('state',seed);
tol = 5e-6; % optimality tolerance for stopping_type 1
DB=20; % SNR of the video
frame_array=[35,100,125]; % frames that will be shown at the end

global D X S
seed = 602;
randn('state',seed); rand('state',seed);
load Hall_airport_1000_1496_497_144_176_gray.mat;
D = images(:,1:201);
n1=144*176; n2=201;
stdev = norm(D,'fro')/(sqrt(144*176*201)*10^(DB/20))
D = D+stdev*randn(144*176,201);
[X,S]=nsa_v1(D,stdev,tol);
figure
plot_data(frame_array,D,X,S,144,176)

author: andrewssobral

repo: lrslibrary


text:
function arglist = hlp_getModelingApproaches(varargin)
% return a cell array (arglist) defining the available modeling approaches
% each entry of the cell array is a cell array of the form
% {'Approach Name' @function_handle}
% where @function_handle is the handle to the function implementing the
% approach
%
% hlp_getModelingApproaches('defNameOnly') returns only the name of the
%  default modeling approach (as a string)
%
% hlp_getModelingApproaches('mfileNameOnly',approachName) where approachName
%   is a string with the name of a valid approach, returns the m-file name
%   of the function implementing the specified approach
%
% Author: Tim Mullen, 2012 SCCN/INC/UCSD

defNameOnly = false;
approachName = '';

if nargin==1 && strcmp(varargin{1},'defNameOnly')
    defNameOnly = true;
end

if nargin==2 && ismember_bc('mfileNameOnly',varargin);
    approachName = varargin{2};
    if ~ischar(approachName)
        error('hlp_getModelingApproaches:badInput','Bad argument pair for ''mfileNameOnly''');
    end
end

% construct the arglist
arglist = { ...
            {'Segmentation VAR' @est_fitMVAR}, ...   % {'Kalman Filtering' {arg('arg_dummy',[],[],'replace this with @est_fitMVARKalman')}}, ...
            };

% determine the output
if defNameOnly
    arglist = arglist{1}{1};
elseif ~isempty(approachName)
    % get the function handle matching the desired approach name
    tmp = cellfun(@(x) fastif(isequal(x{1},approachName),x{2},''),arglist,'UniformOutput',false);
    tmp(cellfun(@isempty,tmp))=[];
    if isempty(tmp)
        error('hlp_getModelingApproachMfileName:badApproachName','Unknown modeling approach ''%s''',approachName);
    else
        % return function name as a string
        arglist = char(tmp{1});
    end
end

repo: REST


text:
function [u,x]=SFrontTrack(u0,x0,T,u_flux,f_flux);
% Doing front tracking on the initial fronts defined by (u0,x0) in the
% time interval [0 T]. The piecewise linear flux function is in (u_flux,
% f_flux). u0 must take values in the set {u_flux}.
% This version plots the fronts as the Riemann problmes are solved.

[u s t_coll,x]=Sinitialize(u0,x0,T,u_flux,f_flux);
% Initializes the collision times.
t_start=zeros(size(x));
asm=mean(abs(s));
xlower=min([x(1)-asm*T,x(1)]);
xhigher=max([x(length(x))+asm*T,x(length(x))]);
if xlower==xhigher
  xlower=xhigher-0.5;
  xhigher=xlower+1;
end;
axis([xlower xhigher 0 T]);
set(gca,'Box','on');
%axes('XLim',[xlower xhigher],'YLim',[0 T],'Box','on');
hold on;
[t i]=min(t_coll);
while t<T&~isempty(s),  % The main loop until no more collisions ...
  [ncr nu]=size(u);
  nx=length(x);
  if nx~=nu-1,
    [nx nu];
    error('Wrong length of x vs u in FrontTrack');
  end;
  xc=x(i)+s(i)*(t-t_start(i));
  line([x(i) xc],[t_start(i) t]);       % Plotting the lines to the
  line([x(i-1) xc],[t_start(i-1) t]);   % collision about to be solved
  [ur sr]=SRiemannsol(u(:,i-1),u(:,i+1),u_flux,f_flux);  % Solving the RP
  [nc nr]=size(ur);
  ns=length(sr);
  if ns>0,
    if i>2,
      t_coll(i-1)=collision_time([s(i-2) sr(1)],[t_start(i-2),t],...
                                 [x(i-2),xc],T+1);
      if t_coll(i-1)<=t,     % The collision is not after present time.
        [-1 t_coll(i-1) t]   % Fatal error.
        [s(i-2) sr(1) t_start(i-2) t]
        error('not increasing t_coll : a');
      end;
    end;
    if i<nu-1,
      t_coll(i+1)=collision_time([sr(ns),s(i+1)],[t,t_start(i+1)],...
                                 [xc,x(i+1)],T+1);
      if t_coll(i+1)<=t,  % The collision is not after present time.
        [1 t_coll(i+1) t] % Fatal error.
        [sr(ns) s(i+1)]
        error('not increasing t_coll: b ');
      end;
    end;
  else
    if i>2&i<nu-1,
      t_coll(i-1)=collision_time([s(i-2) s(i+1)],...
                                 [t_start(i-2),t_start(i+1)],...
                                 [x(i-2),x(i+1)],T+1);
    end;
  end;
  u=[u(:,1:i-1) ur u(:,i+1:nu)];    % Managing the lists ...
  s=[s(1:i-2) sr s(i+1:nu-1)];
  hone=ones([1,ns]);
  x=[x(1:i-2) xc*hone x(i+1:nu-1)];
  t_start=[t_start(1:i-2) t*hone t_start(i+1:nu-1)];
  t_coll=[t_coll(1:i-1) (T+1)*ones([1,nr]) t_coll(i+1:nu)];
  [t i]=min(t_coll);
end;
n=length(x);
for i=1:n  % Drawing the fronts up to final time T
  line([x(i),x(i)+s(i)*(T-t_start(i))],[t_start(i) T]);
end;
x=x+s.*(T-t_start);
hold off;

repo: Aero-matlab


text:
% NMFLAB for Signal Processing written by A. Cichocki and R. Zdunek
% in cooperation with other members of Laboratory for Advanced Brain Signal
% Processing, BSI, RIKEN, Saitama, JAPAN

function [A,X,Distance_output]=nmf_conv(Y,r,Index_norm, A_true, S, mc_on, NoAlts,restart_mc_on,AL,Y_true,type_alg_A, type_alg_X, max_restart, no_conv)
%
% Non-negative Matrix Factorization (NMF) with convolutive algorithms
%
% [A,X]=nmf_IP(Y, r, Index_norm, A_true, S, mc_on, NoAlts, restart_mc_on, AL, Y_true, type_alg_A, type_alg_X, ...
%       max_restart, no_iter, alpha0, tau)
%       produces mixing matrix A of dimension [m by r],
%       and source matrix X of dimension [r by T], for the linear mixing model: AX = Y,
%       where Y is an observation matrix [m by T].
% Note: > m: number of sensors,
%       > r: number of sources,
%       > T: number of samples,
%
% INPUTS:
%       > Index_norm:    vector of 13 binary entries indicating which number
%                        divergence measures are turned on in the View Options,
%
%       > A_true:        true mixing matrix (only for synthetic data),
%       > S:             true source matrix (only for synthetic data),
%       > mc_on:         1 - Monte Carlo analysis enabled, 0 - otherwise,
%       > NoAlts:        number of alternating steps (only for Monte Carlo
%                        analysis and the option "Fixed Alternatings" is selected)
%       > restart_mc_on: 1 - restarts in Monte Carlo analysis are enabled,
%                        0 - otherwise,
%       > AL:            mixing matrix estimaed from the preceeding layer,
%       > Y_true:        the first layer mixed signals (mixtures),
%       > type_alg_A:    indicates the selected algorithm for computation
%                        of the mixing matrix,
%       > type_alg_X:    indicates the selected algorithm for computation of the sources,
%       > max_restart:   number of restarts,
%       > no_iter:       number of inner iterations,
%       > alpha0:        initial magnitude in the exponential model for regularization parameter,
%       > tau:           damping factor in the exponential model for regularization parameter,
%
% OUTPUTS:
%       > A:               estimated mixing matrix,
%       > X:               estimated source matrix,
%       > Distance_output: structures of different divergences measured between "Y" and estimated "AX" versus iterations,
%
% #########################################################################
A = [];
X = [];
if (nargin < 14) | isempty(no_conv) | (no_conv < 0) | max(size(no_conv) > 1)
   disp('Incorrect number of convolutive components');
   return
end
if (nargin < 13) | isempty(max_restart)  | (max_restart < 0) | max(size(max_restart) > 1)
   disp('Number of restarts must be given correctly');
   return
end
if (nargin < 12) | isempty(type_alg_X) | (type_alg_X < 1) | max(size(type_alg_X) > 1)
   disp('Incorrect algorithm for X');
   return
end
if (nargin < 11) | isempty(type_alg_A) | (type_alg_A < 1) | max(size(type_alg_A) > 1)
   disp('Incorrect algorithm for A');
   return
end
if (nargin < 10) | isempty(Y_true)
   disp('The first layer mixed signals are unknown');
   Y_true = zeros(size(Y_true));
end
if (nargin < 9) | isempty(AL)
   disp('Mixing matrix from the preceeding layer unknown');
   AL = eye(size(Y,1));
end
if (nargin < 8) | isempty(restart_mc_on) | max(size(restart_mc_on) > 1)
   disp('Index od restarts in MC analysis unknown');
   restart_mc_on = 0;
end
if (nargin < 7) | isempty(NoAlts) | max(size(NoAlts) > 1)
   disp('Adjustable number of alternatings');
   NoAlts = [];
end
if (nargin < 6) | isempty(mc_on) | max(size(mc_on) > 1)
   disp('No Monte Carlo Analysis');
   mc_on = 0;
end
if (nargin < 5) | isempty(S)
   disp('X_true not given');
end
if (nargin < 4) | isempty(A_true)
   disp('A_true not given');
   index_fixed_A = 1;
else
   index_fixed_A = 0;
end
if (nargin < 3) | isempty(Index_norm)
   '"Index_norm" must be specified'
   return
end
if (nargin < 2) | isempty(r)
   'Rank of factorization must be given'
   return
end

if isempty(Y) | isnan(Y)
   error('No data');
   return
end
% test for negative values in Y
if min(min(Y)) < 0
    disp('Some matrix entries are changed from negative to small positive');
    Y(Y< 0) = eps;
end
if min(sum(Y,2)) == 0
    disp('Not all entries in a row can be zero');
    return
end

[m,T]=size(Y);
niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
niter_sample = 30; % maximum number of iterations for each random sample
epsil_normA = 1E-14; % tolerance for alternating

if mc_on & ~restart_mc_on
    max_restart = 0;
end
if ~isempty(NoAlts)
    niter_selected = NoAlts;
end

if (type_alg_A == 3) & (size(A_true,1) ~= size(Y,1))
   disp('Multilayer technique cannot be used with A fixed');
   Distance_output = [];
   return
end

% Declaration for A and X
A=zeros(m,r,no_conv+1);
Ap = A;
X=zeros(r,T);
Ainit = A;
Xinit = X;
Z = zeros(m,T);
KL_outer_temp = 0;
Z_outer = 0;
nr = 0; restart_on = 0; norm_A = 10; nr_best = -1;
m_sx = 1:m; r_sx = 1:r; T_sx = 1:T;
s_dist = 0;
delta = .1;

while (nr <= max_restart)

   % Initialize random A and X
       if ~nr & (~mc_on | restart_mc_on)
          Ainit_1(m_sx',r_sx) = abs(repmat(.1*sin(2*pi*.1*m_sx'),1,r) + repmat(.1*cos(2*pi*.1*r_sx),m,1) + repmat(cos(2*pi*.471*m_sx'),1,r) + repmat(sin(2*pi*.471*r_sx),m,1));
          Ainit_1 = Ainit_1/max(max(Ainit_1));
          Ainit = repmat(Ainit_1,[1 1 no_conv+1]);

          Xinit(r_sx',T_sx) = abs(repmat(.1*sin(2*pi*.1*r_sx'),1,T) + repmat(.1*cos(2*pi*.1*T_sx),r,1) + repmat(cos(2*pi*.471*r_sx'),1,T) + repmat(sin(2*pi*.471*T_sx),r,1));
          Xinit = Xinit/max(max(Xinit));
       else
          Ainit=rand(m,r,no_conv+1);
          Xinit=rand(r,T);
       end

        % Normalization of initial guess
       Ainit = Ainit.*repmat((1./sum(Ainit,1)),[m,1,1]);

        if (nr == max_restart)&(max_restart > 0)
           A = A_best;
           X = X_best;
        else
           A = Ainit;
           X = Xinit;
        end % initial guess assignment

    Yx = zeros(m,T);
    n = 0; k = 0;

while ((k <= niter_sample)&(nr < max_restart)) | ((k <= niter_selected)&(nr == max_restart)&(norm_A > epsil_normA)& isempty(NoAlts)) | ((k <= niter_selected)&(nr == max_restart)& (NoAlts > 0))

k = k + 1;

l = mod(k-1,no_conv+1)+1;
%Xs = (circshift(X',l-1))';
Xs = X(:, (no_conv-l+2):(T-l+1));
Ys = Y(:, (no_conv-l+2):(T-l+1));

        switch type_alg_A

           case 1 % ALS

              Ap = A;
              Xpinv = Xs'*pinv(Xs*Xs');
              A(:,:,l) = max(1E6*eps, Ys*Xpinv);
              A = A.*repmat((1./sum(A,1)),[m,1,1]);

           case 2 % EMML

              Ap = A;
              A(:,:,l) = A(:,:,l).*((Y./(Yz + eps))*Xs')./repmat(sum(Xs,2)',m,1);
              A = A.*repmat((1./sum(A,1)),[m,1,1]);

           case 3 % Fixed A

              Ap = A;
              niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
              A = A_true + eps;
              A = A.*repmat((1./sum(A,1)),[m,1,1]);

        end % type_alg_A

        switch type_alg_X

            case 1 % ALS

               Xs = max(1E6*eps,pinv(squeeze(A(:,:,l))'*squeeze(A(:,:,l)))*squeeze(A(:,:,l))'*Ys);

            case 2 % EMML

              Zy = (circshift((Y./(Yz + eps))',-l+1))';
               if l > 1
                  Zy(:,(T-l+2):T) = zeros(m,l-1);
               end
              X = X.*(A(:,:,l)'*Zy)./repmat(sum(A(:,:,l),1)',1,T);

            case 3 % Fixed X

              X = S + eps;
              niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)

        end % type_alg_X

      %  X = (circshift(Xs',-l+1))';

                if (nr == max_restart)&(mod(k,50)==0)& (~mc_on | restart_mc_on)
                   norm_A = sqrt(sum(sum(sum((A - Ap).^2))));
                   fprintf(1, 'Restart %d,  %d-th alternating step\n',nr_best+1,k);
                end

                if sum(Index_norm)
                   if (nr == max_restart) & (((k < 50) & (mod(k,5)==0)) | ((k>49) & ((mod(k,50)==0))))

                       s_dist = s_dist + 1;
                       k_select(s_dist) = k;
                       Z = A*X + eps;
                       Z = diag(1./sqrt(var(Z')))*Z;
                       Y = Y + eps;

                        dist_Fro(s_dist) = norm(Y - Z,'fro');
                        dist_KL(s_dist) = sum(sum(Y.*log(Y./Z + eps) - Y + Z));
                        dist_KL2(s_dist) = sum(sum(Z.*log(Z./Y + eps) + Y - Z));
                        dist_Pearson(s_dist) = sum(sum( ((Y - Z).^2)./Z ));
                        dist_Hellinger(s_dist) = sum(sum( (sqrt(Z) - sqrt(Y)).^2 ));
                        dist_JS_rel(s_dist) = sum(sum(2*Y.*log(2*Y./(Y + Z) + eps) + Z - Y));
                        dist_JS_rel2(s_dist) = sum(sum(2*Z.*log(2*Z./(Y + Z) + eps) - Z + Y));
                        Zy = Y + Z;
                        dist_JS(s_dist) = sum(sum(Y.*log(2*Y./Zy + eps) + Z.*log(2*Z./Zy + eps) ));
                        dist_AG_rel(s_dist) = sum(sum(Zy.*log(.5*Zy./Y + eps) + Y - Z));
                        dist_AG(s_dist) = sum(sum(.5*Zy.*log(.5*Zy./sqrt(Y.*Z) + eps)));
                        dist_J(s_dist) = sum(sum( .5*(Y - Z).*log(Y./Z + eps) ));
                        dist_Chi(s_dist) = sum(sum( ((Y + Z).*(Y - Z).^2)./(Y.*Z) ));
                        dist_Tria(s_dist) = sum(sum( ((Y - Z).^2)./(Y + Z) ));

                    end % if multiple
                end % if sum

end % while (k)

% Outer KL divergence
if ndims(A) > 2
   Yz = zeros(m,T);
   for t = 1:size(A,3)
       Ytmp = squeeze(A(:,:,t))*X;
       Yz = Yz + Ytmp;
   end
   Z = AL*Yz;
else
    Z = AL*A*X;
end
Z_outer = norm(Z,'fro') + eps;
KL_outer = sum(sum(Y_true.*log((Y_true + eps)./(Z + eps)) - Y_true + Z))/Z_outer;

           if (nr == 0) | (KL_outer < KL_outer_temp)
              A_best = A; X_best = X; KL_outer_temp = KL_outer; nr_best = nr;
           end % multi-conditions

   nr = nr + 1;

   if nr <=max_restart
      fprintf(1, 'Restart %d, Kullback-Leibler divergence = %e\n',	nr, KL_outer);
   end

end % while (restarts)

% One-Variance scaling
X(X <= 0) = eps;

if ndims(A) > 2
   A =  sum(A,3);
end

Distance_output = cell(length(s_dist),1);
Distance_output(1) = {[]};
Distance_output(2) = {[]};
Distance_output(3) = {[]};
Distance_output(4) = {[]};
Distance_output(5) = {[]};
Distance_output(6) = {[]};
Distance_output(7) = {[]};
Distance_output(8) = {[]};
Distance_output(9) = {[]};
Distance_output(10) = {[]};
Distance_output(11) = {[]};
Distance_output(12) = {[]};
Distance_output(13) = {[]};
Distance_output(14) = {[]};

if sum(Index_norm)
   Distance_output(1) = {k_select};
   Distance_output(2) = {dist_Fro};
   Distance_output(3) = {dist_KL};
   Distance_output(4) = {dist_KL2};
   Distance_output(5) = {dist_Pearson};
   Distance_output(6) = {dist_Hellinger};
   Distance_output(7) = {dist_JS_rel};
   Distance_output(8) = {dist_JS_rel2};
   Distance_output(9) = {dist_JS};
   Distance_output(10) = {dist_AG_rel};
   Distance_output(11) = {dist_AG};
   Distance_output(12) = {dist_J};
   Distance_output(13) = {dist_Chi};
   Distance_output(14) = {dist_Tria};
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [segm_params, filepath_params, other_params] = ...
        internal_params(segm_params, filepath_params, other_params, ...
                       extra_params)
% the internal parameters which are not usually of severe importance to a
% practictioner. Nevertheless, feel free to tinker if you know what you are
% doing :)
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    fp = filepath_params;                                                   % Parameters related to directory locations
    sp = segm_params;                                                       % Parameters for FastSeg's technical internals
    op = other_params;                                                      % Misc. parameters
    ep = extra_params;                                                      % Params passed as an argument for overriding parameters here

    % set the filepath parameters
    fp = get_param_val(fp, ep, 'extern_src_dir', ...
                           fullfile(fp.code_root_dir, 'extern_src'));       % External dependencies directory
    fp = get_param_val(fp, ep, 'fuxin_lib_dir', ...
                           fullfile(fp.code_root_dir, 'extern_src/fuxin_lib_src'));    % External dependencies to Fuxin's library
    fp = get_param_val(fp, ep, 'extern_codes_dir', ...
                           fullfile(fp.code_root_dir, 'extern_src'));       % External dependencies packaged with this code
    fp = get_param_val(fp, ep, 'utils_src_dir', ...
                           fullfile(fp.code_root_dir, 'utils'));            % Directory to our utility functions
    fp = get_param_val(fp, ep, 'seg_save_dirpath', ...
                           fullfile(fp.data_save_dirpath, 'MySegmentsMat'));% Directory where segments computed (results) are saved
   fp = get_param_val(fp, ep, 'trained_models_dirpath', ...
                           fullfile(fp.code_root_dir, 'trained_models'));   % All the trained models (for instance, for boundary energies) reside in this directory
    fp = get_param_val(fp, ep, 'scores_dirpath', ...
                           fullfile(fp.seg_save_dirpath, 'scores'));        % Used by compute_print_results() to save quantitative results for the segments computed
    fp = get_param_val(fp, ep, 'boundaries_parent_dirpath', ...
                           fp.data_save_dirpath);                           % The main directory where boundary detection results are saved
    fp = get_param_val(fp, ep, 'debug_parent_dirpath', ...
                           fullfile(fp.data_save_dirpath, 'Debug'));        % Where any debug information is output (usually when other_params.debug = true)

    % set the Segmenter parameters
    sp = get_param_val(sp, ep, 'bg_frame_width',        1);                 % Decides the no. of pixels to use at the boundary, which is then used to generate the bg seed. A large number would force more superpixels
                                                                            % to be used background seeds. Invoked at Segmenter.precompute_unary_data().
    sp = get_param_val(sp, ep, 'filter_max_energy',     5000);              % Threshold used by Segmenter.filter_segments to filter segments with too high of a cut ratio. The energy is computed in compute_energies().
    sp = get_param_val(sp, ep, 'graph_pairwise_multiplier', ...
            {         1000           ,          1000           });          % Used in AbstractGraph.set_pairwise_graph() as a scaling factor for all pairwise capacities. A larger value will increase pairwise
                                                                            % capacities and hence increase the likelihood for larger segments. Each cell gives a value used for each graph method.
    sp = get_param_val(sp, ep, 'graph_sub_methods_cut_param', ...
            {[     1    ,     1     ], [    1     ,     1     ]});          % For graph-cut settings. Invoked in GraphProcessor.multiseed_param_min_st_cut() as <graphcut_params>. A 0 value here would schedule
                                                                            % parametric lambdas in reverse i.e. if you have lambda values of 0,1,2,5 for paramertic min-cut, a 0 parameter here would compute min-cuts
                                                                            % in the order 5,2,1,0. A parameter 1 schedules it in the normal ascending direction. Both methods give the same results but might effect the
                                                                            % computation time. This can possibly be extended to support other parameter settings for graph-cuts.
    sp = get_param_val(sp, ep, 'graph_filter_segs', ...
            {         true           ,           true          });          % Indicates whether you want to filter the segments produced by a particular graph method. Invoked in Segmenter.filter_segments().
    sp = get_param_val(sp, ep, 'graph_sol_upper_bp', ...
            [         20             ,           300           ]);          % Decides the upper limit for the lambdas used for the parameteric min-cut. This is the <u> value in GraphProcessor.generate_param_lambdas()
                                                                            % See the function to find the exact algorithm used to generate the list of lambdas for parametric min-cut.

    % set miscellaneous parameters
    hash_params.Format = 'base64';
    hash_params.Method = 'SHA-1';
    op = get_param_val(op, ep, 'hash_params', hash_params);                 % Hashing parameters. This is used to create a hash of the settings used - this is useful for figuring out if the same set of parameters were
                                                                            % used for computing a set of segments, without comparing the parameters directly.

    % compute the hash for the parameters (acts as a unique ID for params)
%     addpath(fullfile(fp.extern_codes_dir, 'DataHash'));
%     op.filepath_params_hash = DataHash(fp, hash_params);
%     op.segm_params_hash = DataHash(sp, hash_params);
    op.filepath_params_hash = num2str(sum(getByteStreamFromArray(fp)));
    op.segm_params_hash = num2str(sum(getByteStreamFromArray(sp)));

    filepath_params = fp;
    segm_params = sp;
    other_params = op;
end

author: Cloud-CV

repo: object-proposals


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function pc_all = loadTestingMatFiles(path, mat_files)
    % load pc
    for i = 1:size(mat_files,2)
        testing_pc = loadPointCloud(path, mat_files(i));

        % randomly find a scan to project
        num_scan = size(testing_pc, 1);
        scan = randi([1 num_scan],1);

        points = testing_pc(scan, :, 1:3);
        points = (reshape(points, size(points, 2),[]))';
        for j = 1:size(points, 2)
            if ~any(all(points(:, j), 2))
                    pc_all(i).mat_pc = [points(:, 1:j-1); ...
                        ones(1, size(points(:, 1:j-1), 2))];
                break;
            end
        end
    end
end

author: UMich-BipedLab


text:
function r8sd_print ( n, ndiag, offset, a, title )

%*****************************************************************************80
%
%% R8SD_PRINT prints a R8SD matrix.
%
%  Discussion:
%
%    The R8SD storage format is for symmetric matrices whose only nonzero entries
%    occur along a few diagonals, but for which these diagonals are not all
%    close enough to the main diagonal for band storage to be efficient.
%
%    In that case, we assign the main diagonal the offset value 0, and
%    each successive superdiagonal gets an offset value 1 higher, until
%    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
%
%    Assuming there are NDIAG nonzero diagonals (ignoring subdiagonals!),
%    we then create an array B that has N rows and NDIAG columns, and simply
%    "collapse" the matrix A to the left:
%
%  Example:
%
%    The "offset" value is printed above each column.
%
%    Original matrix               New Matrix
%
%       0   1   2   3   4   5       0   1   3   5
%
%      11  12   0  14   0  16      11  12  14  16
%      21  22  23   0  25   0      22  23  25  --
%       0  32  33  34   0  36      33  34  36  --
%      41   0  43  44  45   0      44  45  --  --
%       0  52   0  54  55  56      55  56  --  --
%      61   0  63   0  65  66      66  --  --  --
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer NDIAG, the number of diagonals of the matrix
%    that are stored in the array.
%    NDIAG must be at least 1, and no more than N.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal storage.
%
%    Input, real A(N,NDIAG), the R8SD matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8sd_print_some ( n, ndiag, offset, a, 1, 1, n, n, title );

  return
end


text:
%kkclrgetargs 'Clear the main_get_args_call field in the main of a kroutine or xvroutine'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kclrgetargs.pane file
%
% Parameters:
% String: tb 'Toolbox Name', default: ' ': 'Toolbox to search'
% String: oname 'Object Name ', default: ' ': 'Software object name'
% Toggle: rgw 'Run Ghostwriter?', default: 0: 'Run ghostwriter after stripping the get args section'
%
% Example: kkclrgetargs( {'tb',' ';'oname',' ';'rgw',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kclrgetargs - Clear the main_get_args_call field in the main of a kroutine or xvroutine
%
%  DESCRIPTION
% This program can be used to clear the -main_get_args_call field in the
% main routine of kroutines and xvroutines.  The optional flag -rgw tells
% kclrgetargs to run the code generators after clearing this field, thereby
% putting the latest version of the call to the clui and pane parsing API.
% It was developed primarily to help migrate kroutines that were developed
% before the Khoros code generators supported the generation of continuous
% run processes.  All new kroutines and xvroutines use the new API clui
% initialize API that supports continuous run, and with this program it
% is easy to update existing kroutines and xvroutines.
%
%
%
%  EXAMPLES
%
% % kclrgetargs -tb testtoolbox -oname myobject -rgw
% The -main_get_args_call section of myobject has been cleared
% ghostwriter: done generating myobject man page
% ghostwriter: done generating myobject.h
% ghostwriter: done generating myobject.c
% ghostwriter: done generating usage.c
% ghostwriter: done generating myobject help page
% ghostwriter: updating object database [TESTTOOLBOX::myobject]
%
%
%  "SEE ALSO"
% ghostwriter(1), kcodegen(3)
%
%  RESTRICTIONS
% Since the -main_get_args_call section of a kroutine or xvroutine is
% one of the fields that users are allowed to edit, running this
% program can potentially remove special code the user put in that
% section so that their program would do some special handling of
% the command line parameters.
%
%  REFERENCES
% None
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkclrgetargs(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kkclrgetargs(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'tb', ' ';'oname', ' ';'rgw', 0};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'String','String','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'kclrgetargs"  '],Inputs,narglist);


text:
%% compile
%mex -O ba_interp2.cpp

%% The mandrills eye
clear
IMG=load('mandrill');
IMG = ind2rgb(IMG.X, IMG.map);
[Dx Dy] = meshgrid(-25:0.1:25, -25:0.1:25);

R = [cos(pi/18) sin(pi/18); -sin(pi/18) cos(pi/18)];
RD = R * [Dx(:)'; Dy(:)'];
RDx = reshape(RD(1,:), size(Dx))+320;
RDy = reshape(RD(2,:), size(Dy))+42;

methods = {'nearest', 'linear', 'cubic'};
la=nan(1,3);
for i=1:3
  tic;
  IMG_R = ba_interp2(IMG, RDx, RDy, methods{i});
  elapsed=toc;

  IMG_R_M = nan(size(IMG_R));
  tic;
  for c=1:size(IMG,3)
    IMG_R_M(:,:,c) = interp2(IMG(:,:,c), RDx, RDy, ['*' methods{i}]);
  end
  elapsed_matlab=toc;

  la(i) = subplot(2,2,i);
  imshow(IMG_R);
  difference = max(max(max(abs(IMG_R-IMG_R_M))));
  title({...
    sprintf('Rotation and zoom using %s interpolation took %gs.', methods{i}, elapsed), ...
    sprintf('Speedup against inbuilt method by a factor of %g', elapsed_matlab/elapsed), ...
    sprintf('The maximum difference with the built in method was %g', difference)});
  end
linkaxes(la(:));

%% The mandrills eye single
clear
IMG=load('mandrill');
IMG = single(ind2rgb(IMG.X, IMG.map));
[Dx Dy] = meshgrid(single(-25:0.1:25), single(-25:0.1:25));

R = single([cos(pi/18) sin(pi/18); -sin(pi/18) cos(pi/18)]);
RD = R * [Dx(:)'; Dy(:)'];
RDx = reshape(RD(1,:), size(Dx))+320;
RDy = reshape(RD(2,:), size(Dy))+42;

methods = {'nearest', 'linear', 'cubic'};
la=nan(1,3);
for i=1:3
  tic;
  IMG_R = ba_interp2(IMG, RDx, RDy, methods{i});
  elapsed=toc;

  IMG_R_M = nan(size(IMG_R));
  tic;
  for c=1:size(IMG,3)
    IMG_R_M(:,:,c) = interp2(IMG(:,:,c), RDx, RDy, ['*' methods{i}]);
  end
  elapsed_matlab=toc;

  la(i) = subplot(2,2,i);
  imshow(IMG_R);
  difference = max(max(max(abs(IMG_R-IMG_R_M))));
  title({...
    sprintf('Rotation and zoom using %s interpolation took %gs.', methods{i}, elapsed), ...
    sprintf('Speedup against inbuilt method by a factor of %g', elapsed_matlab/elapsed), ...
    sprintf('The maximum difference with the built in method was %g', difference)});
  end
linkaxes(la(:));

%% Cumulative effect of resampling.
clear
IMG=load('mandrill');
IMG = ind2rgb(IMG.X, IMG.map);
[Dx Dy] = meshgrid(-349.5:1:349.5, -349.5:1:349.5);

R = [cos(pi/18) sin(pi/18); -sin(pi/18) cos(pi/18)];
RD = R * [Dx(:)'; Dy(:)'];
RDx = reshape(RD(1,:), size(Dx))+350;
RDy = reshape(RD(2,:), size(Dy))+350;

methods = {'nearest', 'linear', 'cubic'};
la=nan(1,numel(methods));
for i=1:numel(methods)
  IMG_R{i} = IMG;
  IMG_R{i}(end+1:700,:,:) = 0;
  IMG_R{i}(:,end+1:700,:) = 0;
  IMG_R{i} = circshift(IMG_R{i}, ceil([[700-size(IMG,1) 700-size(IMG,2)]/2, 0]));
  imshow(IMG_R{i});
end
[X Y] = meshgrid(1:700, 1:700);
mask = ((X-350).^2+(Y-350).^2)<350^2;
for j=1:36
  for i=1:numel(methods)
    tic;
    IMG_R{i} = ba_interp2(IMG_R{i}, RDx, RDy, methods{i});
    IMG_R{i}(repmat(~mask, [1 1 3])) = 0;
    elapsed=toc;

    la(i) = subplot(1,3,i);
    imshow(IMG_R{i});
    title(sprintf('Rotating a mandrill %g times by 10degree using %s interpolation.', j, methods{i}));
    axis([100 600 100 600]);
  end
  drawnow
end
linkaxes(la(:));

for i=1:numel(methods)
  imwrite(IMG_R{i}, sprintf('mandril_%s.png', methods{i}), 'Alpha', double(mask));
end

author: shuochsu

repo: DeepVideoDeblurring


text:
%QUESTION NO:2

%Using the Perceptron Learning Law design a classifier for the following
%problem:
% Class C1 : [-2 2]', [-2 1.5]', [-2 0]', [1 0]' and [3 0]'
% Class C2 : [ 1 3]', [3 3]', [1 2]', [3 2]', and [10 0]'

inp=[-2 -2 -2 1 3 1 3 1 3 10;2 1.5 0 0 0 3 3 2 2 0];
out=[1 1 1 1 1 0 0 0 0 0];
network=newp([-2 10;0 3],1);
network.iw{1}=[0.5 0.5];
network.b{1}=0.5;
y=sim(network,inp);
figure,plot(inp,out,inp,y,'o'),title('Before Training');
axis([-10 20 -2.0 2.0]);
network.trainParam.epochs = 20;
network=train(network,inp,out);
y=sim(network,inp);
figure,plot(inp,out,inp,y,'o'),title('After Training');
axis([-10 20 -2.0 2.0]);
display('Final weight vector and bias values : \n');
Weights=network.iw{1};
Bias=network.b{1};
Weights
Bias
Actual_Desired=[y' out'];
Actual_Desired


text:
% LFUtilCalLensletCam - calibrate a lenslet-based light field camera
%
% Usage:
%
%     LFUtilCalLensletCam
%     LFUtilCalLensletCam( Inputpath )
%     LFUtilCalLensletCam( InputPath, CalOptions )
%     LFUtilProcessCalibrations( [], CalOptions )
%
% All parameters are optional and take on default values as set in the "Defaults" section at the top
% of the implementation. As such, this can be called as a function or run directly by editing the
% code. When calling as a function, pass an empty array "[]" to omit a parameter.
%
% This function recursively crawls through a folder identifying decoded light field images and
% performing a calibration based on those images. It follows the calibration procedure described in:
%
% D. G. Dansereau, O. Pizarro, and S. B. Williams, "Decoding, calibration and rectification for
% lenslet-based plenoptic cameras," in Computer Vision and Pattern Recognition (CVPR), IEEE
% Conference on. IEEE, Jun 2013.
%
% Minor differences from the paper: camera parameters are automatically initialized, so no prior
% knowledge of the camera's parameters is required; the free intrinsic parameters have been reduced
% by two: H(3:4,5) were previously redundant with the camera's extrinsics, and are now automatically
% centered; and the light field indices [i,j,k,l] are 1-based in this implementation, and not
% 0-based as described in the paper.
%
% The calibration produced by this process is saved in a JSON file, by default named 'CalInfo.json',
% in the top-level folder of the calibration images. The calibration includes a 5x5 homogeneous
% intrinsic matrix EstCamIntrinsicsH, a 5-element vector EstCamDistortionV describing its distortion
% parameters, and the lenslet grid model formed from the white image used to decode the
% checkerboard images. Taken together, these relate light field indices [i,j,k,l] to spatial rays
% [s,t,u,v], and can be utilized to rectify light fields, as demonstrated in LFUtilDecodeLytroFolder
% / LFCalRectifyLF.
%
% The input light fields are ideally of a small checkerboard (order mm per square), taken at close
% range, and over a diverse set of poses to maximize the quality of the calibration. The
% checkerboard light field images should be decoded without rectification, and colour correction is
% not required.
%
% Calibration is described in more detail in LFToolbox.pdf, including links to example datasets and
% a walkthrough of a calibration process.
%
% Inputs -- all are optional, see code below for default values :
%
%     InputPath : Path to folder containing decoded checkerboard images -- note the function
%                 operates recursively, i.e. it will search sub-folders.
%
%     CalOptions : struct controlling calibration parameters
%          .ExpectedCheckerSize : Number of checkerboard corners, as recognized by the automatic
%                                 corner detector; edge corners are not recognized, so a standard
%                                 8x8-square chess board yields 7x7 corners
%     .ExpectedCheckerSpacing_m : Physical extents of the checkerboard, in meters
%           .LensletBorderSize : Number of pixels to skip around the edges of lenslets; a low
%                                 value of 1 or 0 is generally appropriate, as invalid pixels are
%                                 automatically skipped
%                   .SaveResult : Set to false to perform a "dry run"
%                  .ShowDisplay : Enables various displays throughout the calibration process
%       .ForceRedoCornerFinding : Forces the corner finding procedure to run, overwriting existing
%                                 results
%                .ForceRedoInit : Forces the parameter initialization step to run, overwriting
%                                 existing results
%                   .OptTolX : Determines when the optimization process terminates. When the
%                              estimted parameter values change by less than this amount, the
%                              optimization terminates. See the Matlab documentation on lsqnonlin,
%                              option `TolX' for more information. The default value of 5e-5 is set
%                              within the LFCalRefine function; a value of 0 means the optimization
%                              never terminates based on this criterion.
%                 .OptTolFun : Similar to OptTolX, except this tolerance deals with the error value.
%                              This corresponds to Matlab's lsqnonlin option `TolFun'. The default
%                              value of 0 is set within the LFCalRefine function, and means the
%                              optimization never terminates based on this criterion.
%
% Output takes the form of saved checkerboard info and calibration files.
%
% Example :
%
%   LFUtilCalLensletCam('.', ...
%     struct('ExpectedCheckerSize', [8,6], 'ExpectedCheckerSpacing_m', 1e-3*[35.1, 35.0]))
%
%   Run from within the top-level path of a set of decoded calibration images of an 8x6 checkerboard
%   of dimensions 35.1x35.0 mm, will carry out all the stages of a calibration. See the toolbox
%   documentation for a more complete example.
%
% User guide: LFToolbox.pdf
% See also:  LFCalFindCheckerCorners, LFCalInit, LFCalRefine, LFUtilDecodeLytroFolder, LFSelectFromDatabase

% Copyright (c) 2013-2020 Donald G. Dansereau

function LFUtilCalLensletCam( InputPath, CalOptions )

%---Tweakables---
InputPath = LFDefaultVal('InputPath', '.');

CalOptions = LFDefaultField( 'CalOptions', 'ExpectedCheckerSize', [19, 19] );
CalOptions = LFDefaultField( 'CalOptions', 'ExpectedCheckerSpacing_m', [3.61, 3.61] * 1e-3 );
CalOptions = LFDefaultField( 'CalOptions', 'LensletBorderSize', 1 );
CalOptions = LFDefaultField( 'CalOptions', 'SaveResult', true );
CalOptions = LFDefaultField( 'CalOptions', 'ForceRedoCornerFinding', false );
CalOptions = LFDefaultField( 'CalOptions', 'ForceRedoInit', false );
CalOptions = LFDefaultField( 'CalOptions', 'ShowDisplay', true );
CalOptions = LFDefaultField( 'CalOptions', 'CalInfoFname', 'CalInfo.json' );

%---Check for previously started calibration---
CalInfoFname = fullfile(InputPath, CalOptions.CalInfoFname);
if( ~CalOptions.ForceRedoInit && exist(CalInfoFname, 'file') )
    fprintf('---File %s already exists\n   Loading calibration state and options\n', CalInfoFname);
    CalOptions = LFStruct2Var( LFReadMetadata(CalInfoFname), 'CalOptions' );
else
    CalOptions.Phase = 'Start';
end

RefineComplete = false; % always at least refine once

%---Step through the calibration phases---
while( ~strcmp(CalOptions.Phase, 'Refine') || ~RefineComplete )
    switch( CalOptions.Phase )

        case 'Start'
            %---Find checkerboard corners---
            CalOptions.Phase = 'Corners';
            CalOptions = LFCalFindCheckerCorners( InputPath, CalOptions );

        case 'Corners'
            %---Initialize calibration process---
            CalOptions.Phase = 'Init';
            CalOptions = LFCalInit( InputPath, CalOptions );

            if( CalOptions.ShowDisplay )
                LFFigure(2);
                clf
                LFCalDispEstPoses( InputPath, CalOptions, [], [0.7,0.7,0.7] );
            end

        case 'Init'
            %---First step of optimization process will exclude distortion---
            CalOptions.Phase = 'NoDistort';
            CalOptions = LFCalRefine( InputPath, CalOptions );

            if( CalOptions.ShowDisplay )
                LFFigure(2);
                LFCalDispEstPoses( InputPath, CalOptions, [], [0,0.7,0] );
            end

        case 'NoDistort'
            %---Next step of optimization process adds distortion---
            CalOptions.Phase = 'WithDistort';
            CalOptions = LFCalRefine( InputPath, CalOptions );
            if( CalOptions.ShowDisplay )
                LFFigure(2);
                LFCalDispEstPoses( InputPath, CalOptions, [], [0,0,1] );
            end

        otherwise
            %---Subsequent calls refine the estimate---
            CalOptions.Phase = 'Refine';
            CalOptions = LFCalRefine( InputPath, CalOptions );
            RefineComplete = true;
            if( CalOptions.ShowDisplay )
                LFFigure(2);
                LFCalDispEstPoses( InputPath, CalOptions, [], [1,0,0] );
            end
    end
end

repo: LFToolbox


text:
%  INTERNAL FUNCTION: makes to vectors have the same minimum and maximum
%
%  ::
%
%    new=apply_range(type,refv,newv)
%
%  Args:
%
%     - **refv** [vector]: reference vector
%     - **newv** [vector]: vector to modify
%
%  Returns:
%     :
%
%     - **newv** [vector]: modified vector
%
%

author: jmaih


text:
function x = jacobi2 ( n, a, b, x )

%*****************************************************************************80
%
%% JACOBI2 carries out one step of the Jacobi iteration.
%
%  Discussion:
%
%    MATLAB's vector notation is used to express the operation in one line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of equations.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real B(N), the right hand side.
%
%    Input, real X(N), the estimated solution.
%
%    Output, real X_NEW(N), the improved estimate of the solution.
%
  x = ( b - a * x + ( diag ( a ) .* x ) ) ./ diag ( a );

  return
end