text:
function I=processvarargin(I,D)
%function I2=processvarargin(I,D)
%
%PURPOSE
%
%Subfunction of various functions in Icasso.
%
%INPUT
%
% I (cell array) input {'id1' value1 'id2' value2 ...}
% D (cell array) defaults {'id1 value1 'id2' value2 'id3' value3 ...}
%
%OUTPUT
%
% I2 (cell array) combination of inputs I and defaults D.
%
%DETAILS
%
%Checks that the format of input I is correct:
%Error if the input array I is not of form
% {string whatever string whatever ...}.
%That is, there are 2N elements in the array and elements with odd
%index are strings (identifiers).
%
%Adds defaults, i.e., adds those id, value pairs in
%cell array D whose id do _not_ appear in input I.

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

% ver 1.2 johan 100105

I=I(:);
D=D(:);
Ninput=length(I);

%%% Check number of pairs
if fix(Ninput/2)~=Ninput/2,
  error('Optional ''identifier'',value pairs mismatch!');
end

%%% Take default ids and values from the cell array

for i=1:2:Ninput,
  if ~ischar(I{i})
    error('Optional ''identifier'',value pairs mismatch!');
  end
  %% Check if some default is redefined, and drop the default away

  index2default=find(strcmp(lower(D(1:2:end)),lower(I{i})));
  if any(index2default);
    index2default=index2default(:)*2-1;
    index2default=[index2default;index2default+1];
    D(index2default)=[];
  end
end

I=[I;D]';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function val = objective(guess,IC)

%Runs a large number of simulations of the controlled simulation and then
%returns a combination of state and actuation error

duration = 3;
dt = 0.05;
n = round(duration/dt);

z = IC;

%Weighting between relative cost.
x_cost = 0.7^-2;
v_cost = 1.5^-2;
u_cost = 1.0^-2;

%Run simulation using RK4 integration in parallel
val = 0; %cost accumulator
for i=1:n
    [k1, u] = rhs(z,guess);
    k2 = rhs(z + 0.5*dt*k1,guess);
    k3 = rhs(z + 0.5*dt*k2,guess);
    k4 = rhs(z + k3,guess);

    z = z + (dt/6)*(k1 + 2*k2 + 2*k3 + k4);

    val = val + x_cost*sum(z(1,:).^2);
    val = val + v_cost*sum(z(2,:).^2);
    val = val + u_cost*sum(u.^2);
end

end

function [xdot, u] = rhs(z,k)

u = control(z,k);
xdot = dynamics(z,u);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% Minimize thermal noise power of an array with arbitrary 2-D geometry
% "Convex optimization examples" lecture notes (EE364) by S. Boyd
% "Antenna array pattern synthesis via convex optimization"
% by H. Lebret and S. Boyd
% (figures are generated)
%
% Designs an antenna array such that:
% - it has unit a sensitivity at some target direction
% - obeys constraint for minimum sidelobe level outside the beamwidth
% - minimizes thermal noise power in y (sigma*||w||_2^2)
%
% This is a convex problem described as:
%
%   minimize   norm(w)
%       s.t.   y(theta_tar) = 1
%              |y(theta)| <= min_sidelobe   for theta outside the beam
%
% where y is the antenna array gain pattern (complex function) and
% variables are w (antenna array weights or shading coefficients).
% Gain pattern is a linear function of w: y(theta) = w'*a(theta)
% for some a(theta) describing antenna array configuration and specs.
%
% Written for CVX by Almir Mutapcic 02/02/06

% select array geometry
ARRAY_GEOMETRY = '2D_RANDOM';
% ARRAY_GEOMETRY = '1D_UNIFORM_LINE';
% ARRAY_GEOMETRY = '2D_UNIFORM_LATTICE';

%********************************************************************
% problem specs
%********************************************************************
lambda = 1;           % wavelength
theta_tar = 60;       % target direction
half_beamwidth = 10;  % half beamwidth around the target direction
min_sidelobe = -20;   % maximum sidelobe level in dB

%********************************************************************
% random array of n antenna elements
%********************************************************************
if strcmp( ARRAY_GEOMETRY, '2D_RANDOM' )
  % set random seed to repeat experiments
  rand('state',0);

  % (uniformly distributed on [0,L]-by-[0,L] square)
  n = 36;
  L = 5;
  loc = L*rand(n,2);

%********************************************************************
% uniform 1D array with n elements with inter-element spacing d
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '1D_UNIFORM_LINE' )
  % (unifrom array on a line)
  n = 30;
  d = 0.45*lambda;
  loc = [d*[0:n-1]' zeros(n,1)];

%********************************************************************
% uniform 2D array with m-by-m element with d spacing
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '2D_UNIFORM_LATTICE' )
  m = 6; n = m^2;
  d = 0.45*lambda;

  loc = zeros(n,2);
  for x = 0:m-1
    for y = 0:m-1
      loc(m*y+x+1,:) = [x y];
    end
  end
  loc = loc*d;

else
  error('Undefined array geometry')
end

%********************************************************************
% construct optimization data
%********************************************************************
% build matrix A that relates w and y(theta), ie, y = A*w
theta = [1:360]';
A = kron(cos(pi*theta/180), loc(:,1)') + kron(sin(pi*theta/180), loc(:,2)');
A = exp(2*pi*i/lambda*A);

% target constraint matrix
[diff_closest, ind_closest] = min( abs(theta - theta_tar) );
Atar = A(ind_closest,:);

% stopband constraint matrix
ind = find(theta <= (theta_tar-half_beamwidth) | ...
           theta >= (theta_tar+half_beamwidth) );
As = A(ind,:);

%********************************************************************
% optimization problem
%********************************************************************
cvx_begin
  variable w(n) complex
  minimize( norm( w ) )
  subject to
    Atar*w == 1;
    abs(As*w) <= 10^(min_sidelobe/20);
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
end

fprintf(1,'The minimum norm of w is %3.2f.\n\n',norm(w));

%********************************************************************
% plots
%********************************************************************
figure(1), clf
plot(loc(:,1),loc(:,2),'o')
title('Antenna locations')

% plot array pattern
y = A*w;

figure(2), clf
ymin = -30; ymax = 0;
plot([1:360], 20*log10(abs(y)), ...
     [theta_tar theta_tar],[ymin ymax],'r--',...
     [theta_tar+half_beamwidth theta_tar+half_beamwidth],[ymin ymax],'g--',...
     [theta_tar-half_beamwidth theta_tar-half_beamwidth],[ymin ymax],'g--',...
     [0 theta_tar-half_beamwidth],[min_sidelobe min_sidelobe],'r--',...
     [theta_tar+half_beamwidth 360],[min_sidelobe min_sidelobe],'r--');
xlabel('look angle'), ylabel('mag y(theta) in dB');
axis([0 360 ymin ymax]);

% polar plot
figure(3), clf
zerodB = 50;
dBY = 20*log10(abs(y)) + zerodB;
plot(dBY.*cos(pi*theta/180), dBY.*sin(pi*theta/180), '-');
axis([-zerodB zerodB -zerodB zerodB]), axis('off'), axis('square')
hold on
plot(zerodB*cos(pi*theta/180),zerodB*sin(pi*theta/180),'k:') % 0 dB
plot( (min_sidelobe + zerodB)*cos(pi*theta/180), ...
      (min_sidelobe + zerodB)*sin(pi*theta/180),'k:')  % min level
text(-zerodB,0,'0 dB')
text(-(min_sidelobe + zerodB),0,sprintf('%0.1f dB',min_sidelobe));
theta_1 = theta_tar+half_beamwidth;
theta_2 = theta_tar-half_beamwidth;
plot([0 55*cos(theta_tar*pi/180)], [0 55*sin(theta_tar*pi/180)], 'k:')
plot([0 55*cos(theta_1*pi/180)], [0 55*sin(theta_1*pi/180)], 'k:')
plot([0 55*cos(theta_2*pi/180)], [0 55*sin(theta_2*pi/180)], 'k:')
hold off

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function value = lrline ( xu, yu, xv1, yv1, xv2, yv2, dv )

%*****************************************************************************80
%
%% LRLINE determines if a point is left of, right or, or on a directed line.
%
%  Discussion:
%
%    The directed line is parallel to, and at a signed distance DV from
%    a directed base line from (XV1,YV1) to (XV2,YV2).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe,
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, real XU, YU, the coordinates of the point whose
%    position relative to the directed line is to be determined.
%
%    Input, real XV1, YV1, XV2, YV2, the coordinates of two points
%    that determine the directed base line.
%
%    Input, real DV, the signed distance of the directed line
%    from the directed base line through the points (XV1,YV1) and (XV2,YV2).
%    DV is positive for a line to the left of the base line.
%
%    Output, integer VALUE, the result:
%    +1, the point is to the right of the directed line;
%     0, the point is on the directed line;
%    -1, the point is to the left of the directed line.
%
  tol = 100.0 * eps;

  dx = xv2 - xv1;
  dy = yv2 - yv1;
  dxu = xu - xv1;
  dyu = yu - yv1;

  tolabs = tol * max ( abs ( dx ), ...
                 max ( abs ( dy ), ...
                 max ( abs ( dxu ), ...
                 max ( abs ( dyu ), abs ( dv ) ) ) ) );

  t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy );

  if ( tolabs < t )
    value = 1;
  elseif ( -tolabs <= t )
    value = 0;
  else
    value = -1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% THRESHOLD_1D_WAVELET: Thresholded 1D wavelet transform.
% Usage
%    [x_phi, x_psi, meta_phi, meta_psi] = threshold_wavelet_1d(x, ...
%       filters, options)
% Input
%    x: The signal to be transformed.
%    filters: The filters of the wavelet transform.
%    options: Various options for the transform. options.oversampling controls
%       the oversampling factor when subsampling. options.threshold specifies
%       the threshold used.
% Output
%    x_phi: x filtered by lowpass filter phi
%    x_psi: cell array of x filtered by wavelets psi
%    meta_phi, meta_psi: meta information on x_phi and x_psi, respectively
% See also  WAVELET_1D

function [x_phi, x_psi, meta_phi, meta_psi] = threshold_wavelet_1d(x, ...
    filters, options)

if nargin < 3
    options = struct();
end

options = fill_struct(options, 'threshold', 1e-4);

[x_phi, x_psi, meta_phi, meta_psi] = wavelet_1d(x, filters, options);

for p1 = find(~cellfun(@isempty, x_psi))
    res = meta_psi.resolution(p1);

    x_psi{p1 ...
        } = T_theta(x_psi{p1}, options.threshold*2^(res/2));
end
end

function y = T_theta(x, theta)
y=sign(x).*max(0,abs(x)-theta);

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function high_card_simulation_test02 ( )

%*****************************************************************************80
%
%% HIGH_CARD_SIMULATION_TEST02 plots the results for a deck of 100 cards.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2012
%
%  Author:
%
%    John Burkardt
%
  deck_size = 100;
  trial_num = 1000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HIGH_CARD_SIMULATION_TEST02\n' );
  fprintf ( 1, '  Using %d cards and %d trials, compute the chances of\n', deck_size, trial_num );
  fprintf ( 1, '  picking the high card with a skip of 0 through 99.\n' );

  p = zeros ( deck_size, 1 );

  for skip_num = 0 : deck_size - 1

    p(skip_num+1) = high_card_simulation ( deck_size, skip_num, trial_num );

  end

  plot ( 0:deck_size-1, p, 'b', 'Linewidth', 2 )
  grid on
  title ( 'Estiamted chance of winning per given skip number' )
  xlabel ( 'Number of cards to skip before choice' )
  ylabel ( 'Chance of correct choice.' )

  print ( '-dpng', 'high_card_simulation_test02.png' )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data plotted in "high_card_simulation_test02.png".\n' );

  return
end


text:
function dtiVisualizeClusterOfFibers(clusterID, clusterlabels, fg)
%Display a group of fibers from fg.
%dtiVisualizeClusterOfFibers(clusterID, clusterlabels, fg)
%e.g., dtiVisualizeClusterOfFibers(1, fg.subgroup, fg)
%clusterlabels is a vector denoting, for each fiber, which cluster the
%fiber belongs to.
%clusterID is which cluster we want to visualize.
%ER 11/2007
if length(clusterlabels)~=1
    clusterlabels=clusterlabels';
end

figure;
for fbindex = find(clusterlabels==clusterID)'
curve=fg.fibers{fbindex};

tubeplot(curve(1, :), curve(2, :), curve(3, :), 1);
hold on;
end
axis([ -80    80  -120    90   -60    90]); xlabel('L<->R'); ylabel('P<->A'); zlabel('I<->S'); grid on;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function x = blend_102 ( r, s, x00, x01, x10, x11 )

%*****************************************************************************80
%
%% BLEND_102 extends scalar point data into a square.
%
%  Diagram:
%
%    01------------11
%     |      .      |
%     |      .      |
%     |.....rs......|
%     |      .      |
%     |      .      |
%    00------------10
%
%  Formula:
%
%    Written in terms of R and S, the map has the form:
%
%      X(R,S) =
%               1     * ( + x00 )
%             + r     * ( - x00 + x10 )
%             + s     * ( - x00       + x01 )
%             + r * s * ( + x00 - x10 - x01 + x11 )
%
%    Written in terms of the coefficients, the map has the form:
%
%      X(R,S) = x00 * ( 1 - r - s + r * s )
%             + x01 * (         s - r * s )
%             + x10 * (     r     - r * s )
%             + x11 * (             r * s )
%
%             = x00 * ( 1 - r ) * ( 1 - s )
%             + x01 * ( 1 - r ) *       s
%             + x10 *       r   * ( 1 - s )
%             + x11 *       r           s
%
%    The nonlinear term ( r * s ) has an important role:
%
%      If ( x01 + x10 - x00 - x11 ) is zero, then the input data lies in
%      a plane, and the mapping is affine.  All the interpolated data
%      will lie on the plane defined by the four corner values.  In
%      particular, on any line through the square, data values at
%      intermediate points will lie between the values at the endpoints.
%
%      If ( x01 + x10 - x00 - x11 ) is not zero, then the input data does
%      not lie in a plane, and the interpolation map is nonlinear.  On
%      any line through the square, data values at intermediate points
%      may lie above or below the data values at the endpoints.  The
%      size of the coefficient of r * s will determine how severe this
%      effect is.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Gordon,
%    Blending-Function Methods of Bivariate and Multivariate Interpolation
%    and Approximation,
%    SIAM Journal on Numerical Analysis,
%    Volume 8, Number 1, March 1971, pages 158-177.
%
%    William Gordon and Charles Hall,
%    Transfinite Element Methods: Blending-Function Interpolation over
%    Arbitrary Curved Element Domains,
%    Numerische Mathematik,
%    Volume 21, Number 1, 1973, pages 109-129.
%
%    William Gordon and Charles Hall,
%    Construction of Curvilinear Coordinate Systems and Application to
%    Mesh Generation,
%    International Journal of Numerical Methods in Engineering,
%    Volume 7, 1973, pages 461-477.
%
%    Joe Thompson, Bharat Soni, Nigel Weatherill,
%    Handbook of Grid Generation,
%    CRC Press, 1999.
%
%  Parameters:
%
%    Input, real R, S, the coordinates where an
%    interpolated value is desired.
%
%    Input, real X00, X01, X10, X11, the data values
%    at the corners.
%
%    Output, real X, the interpolated data value at (R,S).
%
  x =             + x00 ...
      + r *     ( - x00 + x10 ) ...
      + s *     ( - x00       + x01 ) ...
      + r * s * ( + x00 - x10 - x01 + x11 );

  return
end


text:
function[labelImage] = imageInsertBlobLabels(labelImage, labelMap, labelNames, varargin)
% [labelImage] = imageInsertBlobLabels(labelImage, labelMap, labelNames, varargin)
%
% Uses insertText (from Vision Toolbox) to insert a cell of strings into
% the image. Each label is positioned at the center of mass of its blob.
%
% Test case (should output a red 'c' on a white image):
% labelImage = ones(256, 256, 3);
% labelMap = zeros(256, 256);
% labelMap(100, 100) = 3;
% labelMap(101, 100) = 3;
% labelNames = {'a', 'b', 'c'};
% outImage = imageInsertBlobLabels(labelImage, labelMap, labelNames, 'fontColor', [255, 0, 0], 'minComponentSize', 1);
% imshow(outImage);
%
% Copyright by Holger Caesar, 2016

% Parse input
p = inputParser;
addParameter(p, 'fontSize', 15);
addParameter(p, 'fontWidth', 8);
addParameter(p, 'fontColor', [26, 232, 222]);
addParameter(p, 'minComponentSize', 100); % at least 2px
addParameter(p, 'skipLabelInds', []); % labels that are ignored
parse(p, varargin{:});

fontSize = p.Results.fontSize;
fontWidth = p.Results.fontWidth;
fontColor = p.Results.fontColor;
minComponentSize = p.Results.minComponentSize;
skipLabelInds = p.Results.skipLabelInds;

% Check inputs
assert(~isa(labelMap, 'gpuArray'));

% Get unique list of labels
labelList = double(unique(labelMap(:)));
labelList(labelList == 0) = [];
labelListCount = numel(labelList);
if labelListCount == 0
    % Don't do anything if there are no labels
    return;
end

% Init
pixelIdxLists = cell(labelListCount, 1);
pixelIdxLabels = cell(labelListCount, 1);
usedMap = false(size(labelMap));

% Get the indices for all blobs
for labelMapUnIdx = 1 : labelListCount
    labelIdx = labelList(labelMapUnIdx);

    % Get connected components of that label
    components = bwconncomp(labelMap == labelIdx);
    pixelIdxList = components.PixelIdxList(:);
    sel = cellfun(@(x) numel(x), pixelIdxList) >= minComponentSize;
    pixelIdxList = pixelIdxList(sel);
    pixelIdxLabels{labelMapUnIdx} = ones(numel(pixelIdxList), 1) .* labelIdx;
    pixelIdxLists{labelMapUnIdx} = pixelIdxList;
end

% Convert cell to matrix
pixelIdxLists = flattenCellArray(pixelIdxLists);
pixelIdxLabels = cell2mat(pixelIdxLabels);
assert(numel(pixelIdxLists) == numel(pixelIdxLabels));
compCount = numel(pixelIdxLists);

% Sort by size
[~, sortOrder] = sort(cellfun(@numel, pixelIdxLists), 'descend');
pixelIdxLists = pixelIdxLists(sortOrder);
pixelIdxLabels = pixelIdxLabels(sortOrder);

for compIdx = 1 : compCount
    labelIdx = pixelIdxLabels(compIdx);
    if ismember(labelIdx, skipLabelInds)
        continue;
    end
    labelName = labelNames{labelIdx};

    % Get list of relevant pixels
    pixInds = pixelIdxLists{compIdx};
    pixInds = setdiff(pixInds, find(usedMap(:)));
    if isempty(pixInds)
        continue;
    end

    % Compute center of mass
    [y, x] = ind2sub(size(labelMap), pixInds);
    yCenter = median(y);
    xCenter = median(x);

    height = fontSize;
    width = fontWidth * numel(labelName);
    yStart = max(1, round(yCenter - height / 2));
    xStart = max(1, round(xCenter - width / 2));
    yEnd = min(yStart + height - 1, size(usedMap, 1));
    xEnd = min(xStart + width  - 1, size(usedMap, 2));

    if any(any(usedMap(yStart:yEnd, xStart:xEnd)))
        continue;
    end

    % Place label here
    labelImage = insertText(labelImage, [xStart, yStart], labelName, 'Font',  'LucidaSansDemiBold', 'FontSize', fontSize, 'TextColor', fontColor, 'BoxOpacity', 0.0, 'AnchorPoint', 'LeftTop');

    % Mark those pixels as used
    usedMap(yStart:yEnd, xStart:xEnd) = true;
end;

author: nightrome


text:
function [ meta ] = read_ceos_img_meta( filename )
%READ_CEOS_IMG_META Read CEOS SAR image file
%
% Written by: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Open file
fid=fopen(filename,'r','b');

% Read volume descriptor records
meta.rec_num=fread(fid,1,'uint32'); % Record sequence number
meta.rec_subtype1=fread(fid,1,'uint8'); % Record subtype code
meta.rec_type=fread(fid,1,'uint8'); % Record type
meta.rec_subtype2=fread(fid,1,'uint8'); % Record subtype code
meta.rec_subtype3=fread(fid,1,'uint8'); % Record subtype code
meta.rec_length=fread(fid,1,'uint32'); % Record length
meta.ascii_ebcdic=fread(fid,2,'uint8=>char').'; % Record length
fseek(fid,2,'cof'); % Blanks
meta.doc_id=fread(fid,12,'uint8=>char').'; % Superstructure format control document ID
meta.doc_rev=fread(fid,2,'uint8=>char').'; % Superstructure format control document revision level
meta.rec_rev=fread(fid,2,'uint8=>char').'; % Superstructure record format revision level
meta.soft_rel_rev=fread(fid,12,'uint8=>char').'; % Software release and revision level
meta.file_num=fread(fid,4,'uint8=>char').'; % File number
meta.file_id=fread(fid,16,'uint8=>char').'; % File ID
meta.rec_seq_loc_type_flag=fread(fid,4,'uint8=>char').'; % Record sequency and location type flag
meta.seq_num_loc=fread(fid,8,'uint8=>char').'; % Sequence number of location
meta.fld_len_seq=fread(fid,4,'uint8=>char').'; % Field length of sequence number
meta.rec_code_loc_type_flag=fread(fid,4,'uint8=>char').'; % Records code and location type flag
meta.loc_rec_code=fread(fid,8,'uint8=>char').'; % Location of record code
meta.fld_len_code=fread(fid,4,'uint8=>char').'; % Field length of record code
meta.rec_len_loc_type_flag=fread(fid,4,'uint8=>char').'; % Record length and location type flag
meta.loc_rec_len=fread(fid,8,'uint8=>char').'; % Location of record length
meta.len_rec_len=fread(fid,4,'uint8=>char').'; % Location of record length
fseek(fid,68,'cof'); % Blanks
meta.num_data_rec=str2double(fread(fid,6,'uint8=>char')); % Number of data set records
meta.data_len=str2double(fread(fid,6,'uint8=>char')); % Data set summary record length
fseek(fid,24,'cof'); % Blanks
% Sample group data
meta.sample_len=str2double(fread(fid,4,'uint8=>char')); % Bit length per sample
meta.num_samples=str2double(fread(fid,4,'uint8=>char')); % Number of samples per data group
meta.num_bytes=str2double(fread(fid,4,'uint8=>char')); % Number of bytes per data group
meta.just_order=fread(fid,4,'uint8=>char').'; % Justification and order of samples
% SAR related data
meta.num_chan=str2double(fread(fid,4,'uint8=>char')); % Number of SAR channels
meta.num_lines=str2double(fread(fid,8,'uint8=>char')); % Number of lines per data set
meta.num_left=str2double(fread(fid,4,'uint8=>char')); % Number of left border pixels per line
meta.num_pixels=str2double(fread(fid,8,'uint8=>char')); % Number of data groups (or pixels) per line
meta.num_right=str2double(fread(fid,4,'uint8=>char')); % Number of right border pixels per line
meta.num_top=str2double(fread(fid,4,'uint8=>char')); % Number of top border lines
meta.num_bottom=str2double(fread(fid,4,'uint8=>char')); % Number of bottom border lines
meta.interleave=fread(fid,4,'uint8=>char').'; % Interleaving ID
% Record data
meta.phys_rec_line=str2double(fread(fid,2,'uint8=>char')); % Number of physical records per line
meta.phys_rec_multi_chan=str2double(fread(fid,2,'uint8=>char')); % Number of physical records per multi-channel line
meta.prefix_bytes=str2double(fread(fid,4,'uint8=>char')); % Number of bytes of prefix per record
meta.sar_data_bytes=str2double(fread(fid,8,'uint8=>char')); % Number of bytes of SAR data per record
meta.suffix_bytes=str2double(fread(fid,4,'uint8=>char')); % Number of bytes of suffix data per record
meta.pre_suf_rpt_flg=fread(fid,4,'uint8=>char').'; % Prefix/suffix repeat flag
% Prefix/suffix data locations
meta.loc_sar_data=fread(fid,8,'uint8=>char').'; % Sample data line number locator
meta.loc_sar_chan_num=fread(fid,8,'uint8=>char').'; % SAR channel number locator
meta.loc_time=fread(fid,8,'uint8=>char').'; % Time of SAR data line locator
meta.loc_leftfill=fread(fid,8,'uint8=>char').'; % Left-fill count locator
meta.loc_rightfill=fread(fid,8,'uint8=>char').'; % Right-fill count locator
meta.pad_pixels=fread(fid,4,'uint8=>char').'; % Pad pixels present indictor
fseek(fid,28,'cof'); % Blanks
meta.loc_data_qual=fread(fid,8,'uint8=>char').'; % SAR data line quality code locator
meta.loc_cal_info=fread(fid,8,'uint8=>char').'; % Calibration information field locator
meta.loc_gain=fread(fid,8,'uint8=>char').'; % Gain values field locator
meta.loc_bias=fread(fid,8,'uint8=>char').'; % Bias values field lcoator
meta.sar_datatype=fread(fid,28,'uint8=>char').'; % SAR data format type indicator
meta.sar_datatype_code=fread(fid,4,'uint8=>char').'; % SAR data format type code
meta.num_leftfill=fread(fid,4,'uint8=>char').'; % Number of left fill bits within pixel
meta.num_rightfill=fread(fid,4,'uint8=>char').'; % Number of right fill bits within pixel
meta.max_data_range=fread(fid,8,'uint8=>char').'; % Maximum data range of pixel
meta.scansar_num_bursts=fread(fid,4,'uint8=>char').'; % ScanSAR, number of burst data in this file
meta.scansar_num_lines=fread(fid,4,'uint8=>char').'; % ScanSAR, number of lines per one burst
meta.scansar_num_overlap=fread(fid,4,'uint8=>char').'; % ScanSAR, number of overlap lines with adjacent bursts
fseek(fid,260,'cof'); % Blanks
sig_start = ftell(fid);

% Pixel data
switch meta.file_id(8)
    case 'B'
        % Signal data records
        % recs = 1:meta.num_data_rec;
        % Only need the first (and maybe last) record, since these
        % per-record fields are mostly all identical except for a few
        % fields we don't need.
        recs = [1, meta.num_data_rec];
        for i = 1:numel(recs)
            fseek(fid, sig_start + (meta.prefix_bytes + ...
                (meta.num_pixels*meta.num_bytes)) * (recs(i)-1), 'bof');
            meta.signal(i).rec_num=fread(fid,1,'uint32'); % Record sequence number
            meta.signal(i).rec_subtype1=fread(fid,1,'uint8'); % Record subtype code
            meta.signal(i).rec_type=fread(fid,1,'uint8'); % Record type
            meta.signal(i).rec_subtype2=fread(fid,1,'uint8'); % Record subtype code
            meta.signal(i).rec_subtype3=fread(fid,1,'uint8'); % Record subtype code
            meta.signal(i).rec_length=fread(fid,1,'uint32'); % Record length
            % Prefix data-general information
            meta.signal(i).line_num=fread(fid,1,'uint32'); % SAR image data line number
            meta.signal(i).sar_rec_ind=fread(fid,1,'uint32'); % SAR image data record index
            meta.signal(i).left_fill=fread(fid,1,'uint32'); % Actual count of left-fill pixels
            meta.signal(i).num_pixels=fread(fid,1,'uint32'); % Actual count of data pixels
            meta.signal(i).right_fill=fread(fid,1,'uint32'); % Actual count of right-fill pixels
            % Prefix data-sensor parameters
            meta.signal(i).update_flg=fread(fid,1,'uint32'); % Sensor parameters update flag
            meta.signal(i).year=fread(fid,1,'uint32'); % Sensor acquisition year
            meta.signal(i).day=fread(fid,1,'uint32'); % Sensor acquisition day of year
            meta.signal(i).msec=fread(fid,1,'uint32'); % Sensor acquisition milli-seconds of day
            meta.signal(i).chan_id=fread(fid,1,'uint16'); % SAR channel ID
            meta.signal(i).chan_code=fread(fid,1,'uint16'); % SAR channel code
            meta.signal(i).tx_pol=fread(fid,1,'uint16'); % Transmitted polarization
            meta.signal(i).rcv_pol=fread(fid,1,'uint16'); % Received polarization
            meta.signal(i).prf=fread(fid,1,'uint32'); % PRF [mHz]
            meta.signal(i).scan_id=fread(fid,1,'uint32'); % Scan ID
            meta.signal(i).rng_comp_flg=fread(fid,1,'uint16'); % Onboard range compressed flag
            meta.signal(i).chirp_type=fread(fid,1,'uint16'); % Chirp type designator
            meta.signal(i).chirp_length=fread(fid,1,'uint32'); % Chirp length (pulse width) [nsec]
            meta.signal(i).chirp_const=fread(fid,1,'int32'); % Chirp constant coefficient [Hz]
            meta.signal(i).chirp_lin=fread(fid,1,'int32'); % Chirp linear coefficient [Hz/usec]
            meta.signal(i).chirp_quad=fread(fid,1,'int32'); % Chirp quadratic coefficient [Hz/usec^2]
            meta.signal(i).usec=fread(fid,1,'uint64'); % Sensor acquisition micro-second of day (rounded or floored)
            meta.signal(i).gain=fread(fid,1,'uint32'); % Receiver gain [dB]
            meta.signal(i).invalid_flg=fread(fid,1,'uint32'); % Invalid line flag
            meta.signal(i).elec_ele=fread(fid,1,'int32'); % Electronic elevation angle at nadir of antenna [deg]
            meta.signal(i).mech_ele=fread(fid,1,'int32'); % Mechanical elevation angle at nadir of antenna [deg]
            meta.signal(i).elec_squint=fread(fid,1,'int32'); % Electronic antenna squint angle [deg]
            meta.signal(i).mech_squint=fread(fid,1,'int32'); % Mechanical antenna squint angle [deg]
            meta.signal(i).slant_rng=fread(fid,1,'uint32'); % Slant range to 1st data sample [m]
            meta.signal(i).wind_pos=fread(fid,1,'uint32'); % Data record window position (SAMPLE DELAY [nsec])
            fseek(fid,4,'cof'); % Blanks
            % Prefix data-platform reference information
            meta.signal(i).pos_update_flg=fread(fid,1,'uint32'); % Platform position parameters update flag
            meta.signal(i).plat_lat=fread(fid,1,'int32'); % Platform latitude [1/1,000,000 deg]
            meta.signal(i).plat_lon=fread(fid,1,'int32'); % Platform longitude [1/1,000,000 deg]
            meta.signal(i).plat_alt=fread(fid,1,'int32'); % Platform altitude [m]
            meta.signal(i).grnd_spd=fread(fid,1,'int32'); % Platform ground speed [cm/sec]
            meta.signal(i).vel_x=fread(fid,1,'int32'); % Platform velocity X [cm/sec]
            meta.signal(i).vel_y=fread(fid,1,'int32'); % Platform velocity Y [cm/sec]
            meta.signal(i).vel_z=fread(fid,1,'int32'); % Platform velocity Z [cm/sec]
            meta.signal(i).acc_x=fread(fid,1,'int32'); % Platform acceleration X [cm/sec^2]
            meta.signal(i).acc_y=fread(fid,1,'int32'); % Platform acceleration X [cm/sec^2]
            meta.signal(i).acc_z=fread(fid,1,'int32'); % Platform acceleration X [cm/sec^2]
            meta.signal(i).track=fread(fid,1,'int32'); % Platform track angle [1/1,000,000 deg]
            meta.signal(i).true_track=fread(fid,1,'int32'); % Platform rue track angle [1/1,000,000 deg]
            meta.signal(i).pitch=fread(fid,1,'int32'); % Platform pitch angle [1/1,000,000 deg]
            meta.signal(i).roll=fread(fid,1,'int32'); % Platform roll angle [1/1,000,000 deg]
            meta.signal(i).yaw=fread(fid,1,'int32'); % Platform yaw angle [1/1,000,000 deg]
            % Prefix data-sensor/facility specific auxiliary data
            meta.signal(i).lat_first=fread(fid,1,'int32'); % Latitude of 1st pixel [1/1,000,000 deg]
            meta.signal(i).lat_center=fread(fid,1,'int32'); % Latitude of center pixel [1/1,000,000 deg]
            meta.signal(i).lat_last=fread(fid,1,'int32'); % Latitude of last pixel [1/1,000,000 deg]
            meta.signal(i).lon_first=fread(fid,1,'int32'); % Longitude of 1st pixel [1/1,000,000 deg]
            meta.signal(i).lon_center=fread(fid,1,'int32'); % Longitude of center pixel [1/1,000,000 deg]
            meta.signal(i).lon_last=fread(fid,1,'int32'); % Longitude of last pixel [1/1,000,000 deg]
            % ScanSAR burst data parameters
            meta.signal(i).burst_num=fread(fid,1,'uint32'); % ScanSAR, burst number
            meta.signal(i).line_num=fread(fid,1,'uint32'); % ScanSAR, line number in this burst
            fseek(fid,60,'cof'); % Blanks
            meta.signal(i).frame_num=fread(fid,1,'uint32'); % ALOS2 frame number
            % fseek(fid,256,'cof'); % PALSAR auxiliary data
            % fseek(fid,meta.num_pixels*meta.num_bytes,'cof'); % SAR data
        end
    case {'C','D'}
        % Processed data records
        disp('Processed data');
end

% Close file
fclose(fid);

end

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

author: ngageoint


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

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

  switch (type)
   case 0
    y = log(1 + x.^2 / (2 * sigma^2));
   case 1
    y = 2 * x ./ (2 * sigma^2 + x.^2);
   case 2
    y = 2 ./ (2 * sigma^2 + x.^2);      % not second order, but first order/x  (Deqing Sun 24 Nov 2007)
  end

author: kristinbranson

repo: JAABA


text:
function owen_values_test ( )

%*****************************************************************************80
%
%% OWEN_VALUES_TEST demonstrates the use of OWEN_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'OWEN_VALUES_TEST:\n' );
  fprintf ( 1, '  OWEN_VALUES stores values of\n' );
  fprintf ( 1, '  the Owen T function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '           H             A             T\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, h, a, t ] = owen_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %24.16f\n',  h, a, t );

  end

  return
end


text:
function boolVal=isPermutation(possiblePerm,startIdx)
%%ISPERMUTATION Given an array, determine whether it is a valid
%        permutation of the values startIdx:(startIdx+permLength-1).
%
%INPUTS: possiblePerm A length permLength vector to test for being a
%                     permutation.
%            startIdx The starting index of the permutation values to
%                     consider. If this input is omitted or an empty matrix
%                     is passed, then the default of 1 is used.
%
%OUTPUTS: boolVal This is true if possiblePerm is a permutation vector and
%                 false otherwise. Empty matrices are considered
%                 permutation vectors.
%
%After initial checks, the function just allocates a vector having the
%length of possiblePerm and adds 1 to each index given in possiblePerm. If
%any index is repeated. Afterwards, it checks whether all of the indices
%contain a "1".
%
%EXAMPLE:
% boolVal1=isPermutation([1;4;3;5])
% boolVal2=isPermutation([1;4;3;2;5])
% boolVal3=isPermutation([13;16;15;17;14],13)
%One will find boolVal1=false, boolVal2=true, boolVal3=true.
%
%September 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(startIdx))
    startIdx=1;
end

possiblePerm=possiblePerm(:);

if(~isnumeric(possiblePerm)||~isreal(possiblePerm)||any(~isfinite(possiblePerm))||any(possiblePerm~=fix(possiblePerm)))
    boolVal=false;
    return;
end

numEls=length(possiblePerm);

possiblePerm=possiblePerm-startIdx+1;
if(any(possiblePerm>numEls)||any(possiblePerm<1))
    boolVal=false;
    return;
end

numValsPresent=zeros(numEls,1);
for k=1:numEls
    idx=possiblePerm(k);
    numValsPresent(idx)=numValsPresent(idx)+1;
end

boolVal=all(numValsPresent==1);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function headers = checkPLXData(data, headers)
% checkPLXData - Check that PLX data and update/create PLX file headers.
%
% headers = checkPLXData(data, headers)
%
% checkPLXData is designed as a companion function to writePLXFile and
% readPLXHeaders. It takes spiking data that is destined for a PLX file,
% and (optionally) a header produced by readPLXHeaders, and checks that the
% data and header are ready for writing. This function is called by
% writePLXFile, so there is no need to call it first (although this won't
% harm anything). If headers is not provided, common reasonable defaults
% are used instead.
%
% WARNING: Continuous (slow channel) and event data is not yet supported by
% this function. All related data is thrown away and all related fields are
% zeroed out.
%
% Input:
%   data - A structure containing the spiking data
%     data.ts    - [n x 1] array of timestamps (seconds)
%     data.chans - [n x 1] array of channel assignments (1 based)
%     data.units - [n x 1] array of unit assignments (0 = unsorted)
%     data.waves - [n x npw] matrix of waveforms
%   headers - PLX file header information (optional). See readPLXHeaders.
%
% Output:
%   headers - PLX file header information, incorporating information from
%             data structure about spike and unit counts.
%
% See also writePLXFile, readPLXHeaders.
%
% AUTHOR: Benjamin Kraus (bkraus@bu.edu, ben@benkraus.com)
% Copyright (c) 2011, Benjamin Kraus
% $Id$

idstr = 'BJK MATLAB R';
matlabrel = version('-release');

if(length(matlabrel)<=(18-length(idstr)))
    idstr = [idstr matlabrel];
elseif(length(matlabrel)>=15)
    idstr = ['BJK' matlabrel(1:15)];
else
    idstr = [idstr(1:18-length(matlabrel)) matlabrel];
end

% First verify that all fields within the data structure have the same
% number of entries.
if(size(data.ts,1) ~= size(data.chans,1) ||...
        size(data.ts,1) ~= size(data.units,1) ||...
        size(data.ts,1) ~= size(data.waves,1))
    error('All fields in PLX data must have same number of rows.');
end

headers.magic = 'PLEX'; % Required at the head of every PLX file.
headers.version = 106;  % This code is designed to write version 106.
headers.processingsoftware =  idstr;

headers.numPointsWave = size(data.waves,2);

% If headers are not provided, make them up.
if(nargin < 2)
    headers.ADFrequency = 40000;
    headers.lasttimestamp = round(max(data.ts)*headers.ADFrequency);
    headers.numDSPChannels = max(data.chans);
    headers.numPointsPreThr = 0;
    headers.date = now();
    headers.fastread = 0;
    headers.waveformfreq = 40000;
    headers.trodalness = 4;             % Default for version < 103
    headers.datatrodalness = 1;         % Default for version < 103
    headers.bitsperspikesample = 12;    % Default for version < 103
    headers.bitsperslowsample = 12;     % Default for version < 103
    headers.spikeMaxMagnitudeMV = 3000; % Default for version < 103
    headers.slowMaxMagnitudeMV = 5000;  % Default for version < 103
    headers.spikePreAmpGain = 1000;     % Default for version < 105
    headers.acquiringsoftware = '';     % Default for version < 106
    headers.comment='';
    % Channel header information common all channels.
    chanheader.name = '';
    chanheader.signame = '';
    chanheader.channel = 0;
    chanheader.wfrate = 0;
    chanheader.sig = 0;
    chanheader.ref = 0;
    chanheader.gain = 1;
    chanheader.filter = 0;
    chanheader.threshold = 0;
    chanheader.method = 1;
    chanheader.nunits = 0;
    chanheader.template = zeros(5,64);
    chanheader.fit = zeros(5,1);
    chanheader.sortwidth = headers.numPointsWave;
    chanheader.boxes = zeros(5,2,4);
    chanheader.sortbeg = 0;
    chanheader.comment = ''; % Default for version < 105
    chanheader.srcid = '';   % Default for version < 106
    chanheader.chanid = 0;   % Default for version < 106

    % Create channel specific header values.
    for ch = 1:headers.numDSPChannels;
        DSPname = sprintf('sig%03u',ch);
        headers.chans(ch) = chanheader;
        headers.chans(ch).name = DSPname;
        headers.chans(ch).signame = DSPname;
        headers.chans(ch).channel = ch;
        headers.chans(ch).sig = ch;
    end
end

if isfield(headers,'numEventChannels') && isfield(headers,'numSlowChannels')&& ((headers.numEventChannels ~= 0 || headers.numSlowChannels ~= 0))
    wstr = ['Continuous (slow channel) and event data is not yet ',...
        'supported by this function. All related data is thrown away ',...
        'and all related fields are zeroed out.'];
    warning('checkPLXData:unsupported',wstr);
end
headers.numEventChannels = 0; % Event channels not supported
headers.numSlowChannels = 0;  % Slow channels not supported
headers.evcounts = zeros(1, 512);

% Make sure that numDSPChannels matches the number of channel headers.
headers.numDSPChannels = length(headers.chans);

% Shorthand for variables that will be used later.
nDSP = headers.numDSPChannels;
trodalness = headers.trodalness;

switch trodalness
    case 1; elstr = 'electrode';
    case 2; elstr = 'stereotrode';
    case 4; elstr = 'tetrode';
    otherwise; error('Unrecognized trodalness.');
end

% Initialize the timestamp, event, and waveform counts.
% tscounts and wvcounts must be at least [5 x 130]
tscounts = zeros(max(4,max(data.units))+1, max(nDSP,128)+2);

if(max(data.chans)>nDSP)
    error('More spike data channels than header spike channels.');
end

if(mod(nDSP,trodalness)~=0)
    error('Number of spike channel headers must be multiple of trodalness');
end

% Count the number of units on each channel and spikes in each unit.
nspks = zeros(nDSP,1);
tunits = [];
tref = 0;
for ch = 1:nDSP;
    tt = ceil(ch/trodalness);
    units = data.units(data.chans == ch);
    if(isempty(units))
        nspks(ch) = 0;
        headers.chans(ch).nunits = 0;
        tscounts(:,ch+1) = 0;
    else
        % Check if we have a unit assignment for this trode yet. If not,
        % save this unit assignment. If so, check that this is the same
        % unit assignment.
        if(tref ~= tt); tunits = units; tref = tt;
        elseif(length(tunits)~=length(units))
            error('All channels within a %s must have the same number of spikes.',elstr);
        elseif(~all(tunits == units))
            error('All channels within a %s must have the same unit assignments.',elstr);
        end
        nspks(ch) = length(units);
        headers.chans(ch).nunits = max(units);
        for jj = 0:max(data.units);
            tscounts(jj+1,ch+1) = sum(units == jj);
        end
    end
end
headers.tscounts = tscounts;
headers.wfcounts = tscounts;

% Check that each channel on a stereotrode or tetrode has the same number
% of spikes.
if(trodalness >= 2)
    nspks = reshape(nspks,trodalness,[])';
    maxspks = repmat(max(nspks,[],2),1,trodalness);
    nspks(nspks==0) = maxspks(nspks==0);
    if(any(min(nspks,[],2)~=max(nspks,[],2)))
        error('All channels within a %s must have the same number of spikes.',elstr);
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out = drawFromIG(theta, chi)
    [m, n] = size(theta);
    chisq1 = randn(m, n).^2;
    out = theta + 0.5*theta./chi .* (theta.*chisq1 - sqrt(4*theta.*chi.*chisq1 + theta.^2.*chisq1.^2) );
    l = (rand(m, n) >= theta./(theta+out));
    out(l) = theta(l).^2 ./ out(l);
end

author: andrewssobral

repo: lrslibrary


text:
function [hdr] = read_ns_hdr(filename)

% READ_NS_HDR read the header from a NeuroScan 3.x or 4.x AVG/EEG/CNT File
%
% [hdr] = read_ns_hdr(filename)
%
% The output data structure hdr has the fields:
%   hdr.label       - electrode labels
%   hdr.nchan       - number of channels
%   hdr.npnt        - number of samplepoints in ERP waveform
%   hdr.rate        - sample rate (Hz)
%   hdr.xmin        - prestimulus epoch start (e.g., -100 msec)
%   hdr.xmax        - poststimulus epoch end (e.g., 900 msec)
%   hdr.nsweeps     - number of accepted trials/sweeps
%   hdr.domain      - time (0) or frequency (1) domain

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

fid = fopen(filename,'r','ieee-le');

if fid<0,
    error(['cannot open ', filename]);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% these structures are from Neuroscan sethead.h
% the first two columns are the size in bytes and the offset
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% GENERAL SETUP STRUCTURE, total 900 bytes
%    0     20  |  char   rev[20];         /* Revision string                         */
%   20      1  |  char   type;            /* File type AVG=1, EEG=0                  */
%   21     20  |  char   id[20];          /* Patient ID                              */
%   41     20  |  char   oper[20];        /* Operator ID                             */
%   61     20  |  char   doctor[20];      /* Doctor ID                               */
%   81     20  |  char   referral[20];    /* Referral ID                             */
%  101     20  |  char   hospital[20];    /* Hospital ID                             */
%  121     20  |  char   patient[20];     /* Patient name                            */
%  141      2  |  short  int age;         /* Patient Age                             */
%  143      1  |  char   sex;             /* Patient Sex Male='M', Female='F'        */
%  144      1  |  char   hand;            /* Handedness Mixed='M',Rt='R', lft='L'    */
%  145     20  |  char   med[20];         /* Medications                             */
%  165     20  |  char   class[20];       /* Classification                          */
%  185     20  |  char   state[20];       /* Patient wakefulness                     */
%  205     20  |  char   label[20];       /* Session label                           */
%  225     10  |  char   date[10];        /* Session date string                     */
%  235     12  |  char   time[12];        /* Session time strin                      */
%  247      4  |  float  mean_age;        /* Mean age (Group files only)             */
%  251      4  |  float  stdev;           /* Std dev of age (Group files only)       */
%  255      2  |  short int n;            /* Number in group file                    */
%  257     38  |  char   compfile[38];    /* Path and name of comparison file        */
%  295      4  |  float  SpectWinComp;    // Spectral window compensation factor
%  299      4  |  float  MeanAccuracy;    // Average respose accuracy
%  303      4  |  float  MeanLatency;     // Average response latency
%  307     46  |  char   sortfile[46];    /* Path and name of sort file              */
%  353      4  |  int    NumEvents;       // Number of events in eventable
%  357      1  |  char   compoper;        /* Operation used in comparison            */
%  358      1  |  char   avgmode;         /* Set during online averaging             */
%  359      1  |  char   review;          /* Set during review of EEG data           */
%  360      2  |  short unsigned nsweeps;      /* Number of expected sweeps               */
%  362      2  |  short unsigned compsweeps;   /* Number of actual sweeps                 */
%  364      2  |  short unsigned acceptcnt;    /* Number of accepted sweeps               */
%  366      2  |  short unsigned rejectcnt;    /* Number of rejected sweeps               */
%  368      2  |  short unsigned pnts;         /* Number of points per waveform           */
%  370      2  |  short unsigned nchannels;    /* Number of active channels               */
%  372      2  |  short unsigned avgupdate;    /* Frequency of average update             */
%  374      1  |  char  domain;           /* Acquisition domain TIME=0, FREQ=1       */
%  375      1  |  char  variance;         /* Variance data included flag             */
%  376      2  |  unsigned short rate;    /* D-to-A rate                             */
%  378      8  |  double scale;           /* scale factor for calibration            */
%  386      1  |  char  veogcorrect;      /* VEOG corrected flag                     */
%  387      1  |  char  heogcorrect;      /* HEOG corrected flag                     */
%  388      1  |  char  aux1correct;      /* AUX1 corrected flag                     */
%  389      1  |  char  aux2correct;      /* AUX2 corrected flag                     */
%  390      4  |  float veogtrig;         /* VEOG trigger percentage                 */
%  394      4  |  float heogtrig;         /* HEOG trigger percentage                 */
%  398      4  |  float aux1trig;         /* AUX1 trigger percentage                 */
%  402      4  |  float aux2trig;         /* AUX2 trigger percentage                 */
%  406      2  |  short int heogchnl;     /* HEOG channel number                     */
%  408      2  |  short int veogchnl;     /* VEOG channel number                     */
%  410      2  |  short int aux1chnl;     /* AUX1 channel number                     */
%  412      2  |  short int aux2chnl;     /* AUX2 channel number                     */
%  414      1  |  char  veogdir;          /* VEOG trigger direction flag             */
%  415      1  |  char  heogdir;          /* HEOG trigger direction flag             */
%  416      1  |  char  aux1dir;          /* AUX1 trigger direction flag             */
%  417      1  |  char  aux2dir;          /* AUX2 trigger direction flag             */
%  418      2  |  short int veog_n;       /* Number of points per VEOG waveform      */
%  420      2  |  short int heog_n;       /* Number of points per HEOG waveform      */
%  422      2  |  short int aux1_n;       /* Number of points per AUX1 waveform      */
%  424      2  |  short int aux2_n;       /* Number of points per AUX2 waveform      */
%  426      2  |  short int veogmaxcnt;   /* Number of observations per point - VEOG */
%  428      2  |  short int heogmaxcnt;   /* Number of observations per point - HEOG */
%  430      2  |  short int aux1maxcnt;   /* Number of observations per point - AUX1 */
%  432      2  |  short int aux2maxcnt;   /* Number of observations per point - AUX2 */
%  434      1  |  char   veogmethod;      /* Method used to correct VEOG             */
%  435      1  |  char   heogmethod;      /* Method used to correct HEOG             */
%  436      1  |  char   aux1method;      /* Method used to correct AUX1             */
%  437      1  |  char   aux2method;      /* Method used to correct AUX2             */
%  438      4  |  float  AmpSensitivity;  /* External Amplifier gain                 */
%  442      1  |  char   LowPass;         /* Toggle for Amp Low pass filter          */
%  443      1  |  char   HighPass;        /* Toggle for Amp High pass filter         */
%  444      1  |  char   Notch;           /* Toggle for Amp Notch state              */
%  445      1  |  char   AutoClipAdd;     /* AutoAdd on clip                         */
%  446      1  |  char   baseline;        /* Baseline correct flag                   */
%  447      4  |  float  offstart;        /* Start point for baseline correction     */
%  451      4  |  float  offstop;         /* Stop point for baseline correction      */
%  455      1  |  char   reject;          /* Auto reject flag                        */
%  456      4  |  float  rejstart;        /* Auto reject start point                 */
%  460      4  |  float  rejstop;         /* Auto reject stop point                  */
%  464      4  |  float  rejmin;          /* Auto reject minimum value               */
%  468      4  |  float  rejmax;          /* Auto reject maximum value               */
%  472      1  |  char   trigtype;        /* Trigger type                            */
%  473      4  |  float  trigval;         /* Trigger value                           */
%  477      1  |  char   trigchnl;        /* Trigger channel                         */
%  478      2  |  short int trigmask;     /* Wait value for LPT port                 */
%  480      4  |  float trigisi;          /* Interstimulus interval (INT trigger)    */
%  484      4  |  float trigmin;          /* Min trigger out voltage (start of pulse)*/
%  488      4  |  float trigmax;          /* Max trigger out voltage (during pulse)  */
%  492      1  |  char  trigdir;          /* Duration of trigger out pulse           */
%  493      1  |  char  Autoscale;        /* Autoscale on average                    */
%  494      2  |  short int n2;           /* Number in group 2 (MANOVA)              */
%  496      1  |  char  dir;              /* Negative display up or down             */
%  497      4  |  float dispmin;          /* Display minimum (Yaxis)                 */
%  501      4  |  float dispmax;          /* Display maximum (Yaxis)                 */
%  505      4  |  float xmin;             /* X axis minimum (epoch start in sec)     */
%  509      4  |  float xmax;             /* X axis maximum (epoch stop in sec)      */
%  513      4  |  float AutoMin;          /* Autoscale minimum                       */
%  517      4  |  float AutoMax;          /* Autoscale maximum                       */
%  521      4  |  float zmin;             /* Z axis minimum - Not currently used     */
%  525      4  |  float zmax;             /* Z axis maximum - Not currently used     */
%  529      4  |  float lowcut;           /* Archival value - low cut on external amp*/
%  533      4  |  float highcut;          /* Archival value - Hi cut on external amp */
%  537      1  |  char  common;           /* Common mode rejection flag              */
%  538      1  |  char  savemode;         /* Save mode EEG AVG or BOTH               */
%  539      1  |  char  manmode;          /* Manual rejection of incomming data      */
%  540     10  |  char  ref[10];          /* Label for reference electode            */
%  550      1  |  char  Rectify;          /* Rectification on external channel       */
%  551      4  |  float DisplayXmin;      /* Minimun for X-axis display              */
%  555      4  |  float DisplayXmax;      /* Maximum for X-axis display              */
%  559      1  |  char  phase;            /* flag for phase computation              */
%  560     16  |  char  screen[16];       /* Screen overlay path name                */
%  576      2  |  short int CalMode;      /* Calibration mode                        */
%  578      2  |  short int CalMethod;    /* Calibration method                      */
%  580      2  |  short int CalUpdate;    /* Calibration update rate                 */
%  582      2  |  short int CalBaseline;  /* Baseline correction during cal          */
%  584      2  |  short int CalSweeps;    /* Number of calibration sweeps            */
%  586      4  |  float CalAttenuator;    /* Attenuator value for calibration        */
%  590      4  |  float CalPulseVolt;     /* Voltage for calibration pulse           */
%  594      4  |  float CalPulseStart;    /* Start time for pulse                    */
%  598      4  |  float CalPulseStop;     /* Stop time for pulse                     */
%  602      4  |  float CalFreq;          /* Sweep frequency                         */
%  606     34  |  char  taskfile[34];     /* Task file name                          */
%  640     34  |  char  seqfile[34];      /* Sequence file path name                 */
%  674      1  |  char  SpectMethod;      // Spectral method
%  675      1  |  char  SpectScaling;     // Scaling employed
%  676      1  |  char  SpectWindow;      // Window employed
%  677      4  |  float SpectWinLength;   // Length of window %
%  681      1  |  char  SpectOrder;       // Order of Filter for Max Entropy method
%  682      1  |  char  NotchFilter;      // Notch Filter in or out
%  683     11  |  char  unused[11];       // Free space
%  694      2  |  short  FspStopMethod;   /* FSP - Stoping mode                      */
%  696      2  |  short  FspStopMode;     /* FSP - Stoping mode                      */
%  698      4  |  float FspFValue;        /* FSP - F value to stop terminate         */
%  702      2  |  short int FspPoint;     /* FSP - Single point location             */
%  704      2  |  short int FspBlockSize; /* FSP - block size for averaging          */
%  706      2  |  unsigned short FspP1;   /* FSP - Start of window                   */
%  708      2  |  unsigned short FspP2;   /* FSP - Stop  of window                   */
%  710      4  |  float FspAlpha;         /* FSP - Alpha value                       */
%  714      4  |  float FspNoise;         /* FSP - Signal to ratio value             */
%  718      2  |  short int FspV1;        /* FSP - degrees of freedom                */
%  720     40  |  char  montage[40];      /* Montage file path name                  */
%  760     40  |  char  EventFile[40];    /* Event file path name                    */
%  800      4  |  float fratio;           /* Correction factor for spectral array    */
%  804      1  |  char  minor_rev;        /* Current minor revision                  */
%  805      2  |  short int eegupdate;    /* How often incomming eeg is refreshed    */
%  807      1  |  char   compressed;      /* Data compression flag                   */
%  808      4  |  float  xscale;          /* X position for scale box - Not used     */
%  812      4  |  float  yscale;          /* Y position for scale box - Not used     */
%  816      4  |  float  xsize;           /* Waveform size X direction               */
%  820      4  |  float  ysize;           /* Waveform size Y direction               */
%  824      1  |  char   ACmode;          /* Set SYNAP into AC mode                  */
%  825      1  |  unsigned char   CommonChnl;      /* Channel for common waveform             */
%  826      1  |  char   Xtics;           /* Scale tool- 'tic' flag in X direction   */
%  827      1  |  char   Xrange;          /* Scale tool- range (ms,sec,Hz) flag X dir*/
%  828      1  |  char   Ytics;           /* Scale tool- 'tic' flag in Y direction   */
%  829      1  |  char   Yrange;          /* Scale tool- range (uV, V) flag Y dir    */
%  830      4  |  float  XScaleValue;     /* Scale tool- value for X dir             */
%  834      4  |  float  XScaleInterval;  /* Scale tool- interval between tics X dir */
%  838      4  |  float  YScaleValue;     /* Scale tool- value for Y dir             */
%  842      4  |  float  YScaleInterval;  /* Scale tool- interval between tics Y dir */
%  846      4  |  float  ScaleToolX1;     /* Scale tool- upper left hand screen pos  */
%  850      4  |  float  ScaleToolY1;     /* Scale tool- upper left hand screen pos  */
%  854      4  |  float  ScaleToolX2;     /* Scale tool- lower right hand screen pos */
%  858      4  |  float  ScaleToolY2;     /* Scale tool- lower right hand screen pos */
%  862      2  |  short int port;         /* Port address for external triggering    */
%  864      4  |  long  NumSamples;       /* Number of samples in continous file     */
%  868      1  |  char  FilterFlag;       /* Indicates that file has been filtered   */
%  869      4  |  float LowCutoff;        /* Low frequency cutoff                    */
%  873      2  |  short int LowPoles;     /* Number of poles                         */
%  875      4  |  float HighCutoff;       /* High frequency cutoff                   */
%  879      2  |  short int HighPoles;    /* High cutoff number of poles             */
%  881      1  |  char  FilterType;       /* Bandpass=0 Notch=1 Highpass=2 Lowpass=3 */
%  882      1  |  char  FilterDomain;     /* Frequency=0 Time=1                      */
%  883      1  |  char  SnrFlag;          /* SNR computation flag                    */
%  884      1  |  char  CoherenceFlag;    /* Coherence has been  computed            */
%  885      1  |  char  ContinousType;    /* Method used to capture events in *.cnt  */
%  886      4  |  long  EventTablePos;    /* Position of event table                 */
%  890      4  |  float ContinousSeconds; // Number of seconds to displayed per page
%  894      4  |  long  ChannelOffset;    // Block size of one channel in SYNAMPS
%  898      1  |  char  AutoCorrectFlag;  // Autocorrect of DC values
%  899      1  |  unsigned char DCThreshold; // Auto correct of DC level

% ELECTRODE STRUCTURE, total 75 bytes
%   10     0   |  char  lab[10];          /* Electrode label - last bye contains NULL */
%    1    10   |  char  reference;        /* Reference electrode number               */
%    1    11   |  char  skip;             /* Skip electrode flag ON=1 OFF=0           */
%    1    12   |  char  reject;           /* Artifact reject flag                     */
%    1    13   |  char  display;          /* Display flag for 'STACK' display         */
%    1    14   |  char  bad;              /* Bad electrode flag                       */
%    2    15   |  unsigned short int n;   /* Number of observations                   */
%    1    17   |  char  avg_reference;    /* Average reference status                 */
%    1    18   |  char  ClipAdd;          /* Automatically add to clipboard           */
%    4    19   |  float x_coord;          /* X screen coord. for 'TOP' display        */
%    4    23   |  float y_coord;          /* Y screen coord. for 'TOP' display        */
%    4    27   |  float veog_wt;          /* VEOG correction weight                   */
%    4    31   |  float veog_std;         /* VEOG std dev. for weight                 */
%    4    35   |  float snr;              /* signal-to-noise statistic                */
%    4    39   |  float heog_wt;          /* HEOG Correction weight                   */
%    4    43   |  float heog_std;         /* HEOG Std dev. for weight                 */
%    2    47   |  short int baseline;     /* Baseline correction value in raw ad units*/
%    1    49   |  char  Filtered;         /* Toggel indicating file has be filtered   */
%    1    50   |  char  Fsp;              /* Extra data                               */
%    4    51   |  float aux1_wt;          /* AUX1 Correction weight                   */
%    4    54   |  float aux1_std;         /* AUX1 Std dev. for weight                 */
%    4    59   |  float sensitivity;      /* electrode sensitivity                    */
%    1    63   |  char  Gain;             /* Amplifier gain                           */
%    1    64   |  char  HiPass;           /* Hi Pass value                            */
%    1    65   |  char  LoPass;           /* Lo Pass value                            */
%    1    66   |  unsigned char Page;     /* Display page                             */
%    1    67   |  unsigned char Size;     /* Electrode window display size            */
%    1    68   |  unsigned char Impedance;/* Impedance test                           */
%    1    69   |  unsigned char PhysicalChnl; /* Physical channel used                    */
%    1    70   |  char  Rectify;          /* Free space                               */
%    4    71   |  float calib;            /* Calibration factor                       */

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read general parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
offset_type       =  20;
offset_nsweeps    = 364;
offset_pnts       = 368;
offset_nchans     = 370;
offset_domain     = 374;
offset_variance   = 375;
offset_rate       = 376;
offset_xmin       = 505;
offset_xmax       = 509;
fseek(fid, offset_type,     'bof');        type      = fread(fid, 1, 'ushort');
fseek(fid, offset_nsweeps,  'bof');    hdr.nsweeps   = fread(fid, 1, 'ushort');
fseek(fid, offset_pnts,     'bof');    hdr.npnt      = fread(fid, 1, 'ushort');
fseek(fid, offset_nchans,   'bof');    hdr.nchan     = fread(fid, 1, 'ushort');
fseek(fid, offset_domain,   'bof');    hdr.domain    = fread(fid, 1, 'uchar');
fseek(fid, offset_variance, 'bof');    hdr.variance  = fread(fid, 1, 'uchar');
fseek(fid, offset_rate,     'bof');    hdr.rate      = fread(fid, 1, 'ushort');
fseek(fid, offset_xmin,     'bof');    hdr.xmin      = fread(fid, 1, 'float32') * 1000;
fseek(fid, offset_xmax,     'bof');    hdr.xmax      = fread(fid, 1, 'float32') * 1000;
fseek(fid, offset_xmax,     'bof');    hdr.xmax      = fread(fid, 1, 'float32') * 1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read electrode configuration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hdr.label = {};
for elec = 1:hdr.nchan,
  offset_label    = 900+(elec-1)*75 + 0;
  offset_baseline = 900+(elec-1)*75 + 47;
  offset_sens     = 900+(elec-1)*75 + 59;
  offset_calib    = 900+(elec-1)*75 + 71;
  fseek(fid, offset_label,    'bof'); lab                   = fread(fid, 10, 'uchar'); lab(find(lab==0)) = ' ';
  fseek(fid, offset_baseline, 'bof'); hdr.baseline(elec)    = fread(fid, 1, 'ushort');
  fseek(fid, offset_sens,     'bof'); hdr.sensitivity(elec) = fread(fid, 1, 'float32');
  fseek(fid, offset_calib,    'bof'); hdr.calib(elec)       = fread(fid, 1, 'float32');
  hdr.label{elec}  = fliplr(deblank(fliplr(deblank(char(lab')))));
  hdr.factor(elec) = hdr.calib(elec) * hdr.sensitivity(elec) / 204.8;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to infer whether time or frequency data is represented
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if hdr.domain==0
  % data in time domain
elseif hdr.domain==1
  % data in frequency domain
else
  % probably old datafile, assume the data to be in time domain
  warning('assuming the data to be in time domain (domain was %d)', hdr.domain);
  hdr.domain=0;
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function ops = buildRegOps(ops)

% number of channels in db.expts that aren't db.expred
ops.nchannels      = getOr(ops, 'nchannels', 1);
ops.nplanes        = getOr(ops, 'nplanes', 1);

% --- registration options --- %
ops.kriging         = getOr(ops, 'kriging', 1);  % kriging subpixel alignment
ops.PhaseCorrelation = getOr(ops, 'PhaseCorrelation', 1); % set to 0 for non-whitened cross-correlation
if ~ops.kriging
    ops.registrationUpsample = getOr(ops, {'registrationUpsample'}, 1);  % upsampling factor during registration
    % 1 for no upsampling is much faster, 2 may give better subpixel accuracy
end
ops.NimgFirstRegistration  = getOr(ops, 'NimgFirstRegistration', 500); % number of images to include in the first registration pass
ops.NiterPrealign   = getOr(ops, {'NiterPrealign'}, 10); % number of iterations for initial mean image
ops.showTargetRegistration = getOr(ops, 'showTargetRegistration', 1); % shows the image targets for all planes to be registered

% --- non-rigid registration splits FOV into multiple blocks, registers
% each block separately, and then smooths shifts over blocks --- %
% numBlocks = [y-blocks x-blocks]
if getOr(ops, 'nonrigid', 0)
    % default is 8 blocks in the Y direction (1/6 pixels each)
    % ops.numBlocks(1) = # of blocks in Y, ops.numBlocks(2) = # of blocks in X
    ops.numBlocks      = getOr(ops, {'numBlocks'}, [8 1]);
    disp('non-rigid registration chosen');
elseif isfield(ops, 'numBlocks') && ~isempty(ops.numBlocks)
    % if numBlocks is set to a value greater than 1, then run non-rigid
    if numel(ops.numBlocks) == 1
        ops.numBlocks = [ops.numBlocks 1];
    end
    if sum(ops.numBlocks) > 2
        ops.nonrigid               = 1;
        disp('non-rigid registration chosen');
    end
else
    ops.numBlocks = [1 1];
    ops.nonrigid = 0;
end

% --- keep frames at beginning and end of expts (to check for movement)
ops.nimgbegend     = getOr(ops, {'nimgbegend'}, 0);

% --- write tiffs to RegFileTiffLocation if field is not empty --- %
ops.RegFileTiffLocation = getOr(ops, {'RegFileTiffLocation'}, []);
% --- copy binary file to regfilebinlocation if field is not empty --- %
ops.RegFileBinLocation = getOr(ops, {'RegFileBinLocation'}, []);

% --- for aligning across planes --- %
ops.alignTargetImages  = getOr(ops, {'alignTargetImages'}, false); % if true, align target images to each other
ops.interpolateAcrossPlanes = getOr(ops, {'interpolateAcrossPlanes'}, false); %if true, similar looking planes will be averaged together to generate final movie
ops.planesToInterpolate = getOr(ops, {'planesToInterpolate'}, 1:ops.nplanes); % these planes will be considered for interpolation
ops.alignAcrossPlanes  = getOr(ops, {'alignAcrossPlanes'}, false); % at each time point, frame will be aligned to best matching target image (from different planes)

% --- for within plane splitting (if tiffs are too large in Y and X) --- %
ops.splitFOV           = getOr(ops, {'splitFOV'}, [1 1]); % split FOV into chunks if memory issue
% ops.splitFOV(1) = # of subsets in Y, ops.splitFOV(2) = # of subsets in X
ops.smooth_time_space  = getOr(ops, 'smooth_time_space', []);

% --- bidirectional phase offset computation --- %
% assumes same bidirectional offset for each plane
ops.dobidi             = getOr(ops, {'dobidi'}, 1); % compute bidiphase?
% if set to a value by user, do not recompute
if isfield(ops, 'BiDiPhase')
    ops.dobidi         = 0;
end
ops.BiDiPhase          = getOr(ops, {'BiDiPhase'}, 0); % set to default 0

author: cortex-lab


text:
function [s,success] = dmc_xml2struct(file)

%Convert xml file into a MATLAB structure
%  [S,SUCCESS]= DMC_XML2STRUCT(SOURCE) Reads xml-formatted SOURCE and returns a
%  structure S. S contains the attributes in SOURCE, nested to reflect the
%  xml children as appropriate. If SOURCE cannot be parsed for any reason
%  the s is returned as an empty matrix [] and SUCCESS is set to 0.
%  Otherwise SUCCESS is 1.
%
%  EXAMPLE
%  A file containing:
%    <XMLname attrib1="Some value">
%      <Element>Some text</Element>
%      <DifferentElement attrib2="2">Some more text</Element>
%      <DifferentElement attrib3="2" attrib4="1">Even more text</DifferentElement>
%    </XMLname>
%
%  Will produce:
%    s.XMLname.Attributes.attrib1 = "Some value";
%    s.XMLname.Element.Text = "Some text";
%    s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2";
%    s.XMLname.DifferentElement{1}.Text = "Some more text";
%    s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2";
%    s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1";
%    s.XMLname.DifferentElement{2}.Text = "Even more text";
%
%  Note the characters : - and . are not supported in structure fieldnames and
%  are replaced by _

%  Modified from function xml2struct (Matlab exchange File ID: #28518)
%  Copyright by W. Falkena. distributed under BSD license.
%  Last accessed Feb. 3, 2011
%  Modifications here allow xml loading from URLs and return null
%  result when file or URL is invalid.
%
%  MODIFIED BY: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
%  $Date$
%  $Revision$

%check for URL for file name
if (nargin~=1)
    error('xml2struct:IncorrectNumberOfArguments','Incorrect number of arguments');
end
% VALID = 0;
% if strfind(file,'http://')
%     VALID = 1;
% elseif exist(file,'file')
%     VALID = 1;
% elseif exist([file '.xml'],'file')
%     VALID = 1;
% else
%     warning('xml2struct:argumentIsNotFileOrURL',['URL or filename not recognized: ' file]);
% end

try
    xDoc = xmlread(file);
    s = parseChildNodes(xDoc);
    success = 1;
catch
    %warning('xml2struct:InvalidFileOrURL',['No xml code could be parsed from URL (or filename): ' file]);
    s = [];
    success = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Subfunction parseChildNodes
function [children,ptext] = parseChildNodes(theNode)
% Recurse over node children.
children = struct;
ptext = [];
if theNode.hasChildNodes
    childNodes = theNode.getChildNodes;
    numChildNodes = childNodes.getLength;

    for count = 1:numChildNodes
        theChild = childNodes.item(count-1);
        [text,name,attr,childs] = getNodeData(theChild);

        if (~strcmp(name,'#text') && ~strcmp(name,'#comment'))
            %XML allows the same elements to be defined multiple times,
            %put each in a different cell
            if (isfield(children,name))
                if (~iscell(children.(name)))
                    %put existsing element into cell format
                    children.(name) = {children.(name)};
                end
                index = length(children.(name))+1;
                %add new element
                children.(name){index} = childs;
                if(~isempty(text))
                    children.(name){index}.('Text') = text;
                end
                if(~isempty(attr))
                    children.(name){index}.('Attributes') = attr;
                end
            else
                %add previously unkown new element to the structure
                children.(name) = childs;
                if(~isempty(text))
                    children.(name).('Text') = text;
                end
                if(~isempty(attr))
                    children.(name).('Attributes') = attr;
                end
            end
        elseif (strcmp(name,'#text'))
            %this is the text in an element (i.e. the parentNode)
            if (~isempty(regexprep(text,'[\s]*','')))
                if (isempty(ptext))
                    ptext = text;
                else
                    %what to do when element data is as follows:
                    %<element>Text <!--Comment--> More text</element>

                    %put the text in different cells:
                    % if (~iscell(ptext)) ptext = {ptext}; end
                    % ptext{length(ptext)+1} = text;

                    %just append the text
                    ptext = [ptext text];
                end
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GET NODE DATA
function [text,name,attr,childs] = getNodeData(theNode)
% Create structure of node info.

%make sure name is allowed as structure name
name = regexprep(char(theNode.getNodeName),'[-:.]','_');

attr = parseAttributes(theNode);
if (isempty(fieldnames(attr)))
    attr = [];
end

%parse child nodes
[childs,text] = parseChildNodes(theNode);

if (isempty(fieldnames(childs)))
    %get the data of any childless nodes
    try
        %faster then if any(strcmp(methods(theNode), 'getData'))
        text = char(theNode.getData);
    catch
        %no data
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARSE THE ATTRIBURES
function attributes = parseAttributes(theNode)
% Create attributes structure.

attributes = struct;
if theNode.hasAttributes
    theAttributes = theNode.getAttributes;
    numAttributes = theAttributes.getLength;

    for count = 1:numAttributes
        attrib = theAttributes.item(count-1);
        attr_name = regexprep(char(attrib.getName),'[-:.]','_');
        attributes.(attr_name) = char(attrib.getValue);
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function z=fisherz(r);
% function z=fisherz(r);
%
% transforms r scores (correlations) into fisher z scores by the
% transform r = 0.5*log(1+r)/(1-r);
%
% rmk, 1/14/99

z=0.5*log((1+r)./(1-r));


text:
function linpack_c_test14 ( )

%*****************************************************************************80
%
%% TEST14 tests CHPCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST14\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  Hermitian matrix using packed storage (HP),\n' );
  fprintf ( 1, '  CHPCO factors the matrix and estimates\n' );
  fprintf ( 1, '  the reciprocal condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  k = 0;
  seed = 123456789;

  for j = 1 : n

    for i = 1 : j-1
      k = k + 1;
      [ a(k), seed ] = c4_uniform_01 ( seed );
      a_save(i,j) = a(k);
      a_save(j,i) = conj ( a(k) );
    end

    k = k + 1;
    [ a(k), seed ] = r4_uniform_01 ( seed );
    a_save(j,j) = a(k);

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a_save(i,j) ), imag ( a_save(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, rcond ] = chpco ( a, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition RCOND = %f\n', rcond );

  return
end


text:
function [ vX, paramLambda ] = SolveLsNormSquaredConst( mA, vB, normSquaredConst, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsNormConst( mA, vB, normConst )
% Solves norm constrained Least Squares problem by finding the optimal Dual
% Variable (paramLambda) of the KKT Conditions by succesively solving
% Tikhonov Regularized Least Squares problems.
% The objective Funciton is given by:
% \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2}
% subject to {\left\| x \right\|}_{2}^{2} \leq normSquaredConst
% Input:
%   - mA                -   Input Matrix.
%                           The given matrix of the problem 0.5 || A x - b|| ^ 2.
%                           Structure: Vector (Column Vector).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The given vector of the problem 0.5 || A x - b|| ^ 2.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - normConst         -   Norm Constraint.
%                           The solution must obey || x || <= normConst.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
% Output:
%   - vX                -   Solution Vector.
%                           The optimal solution of the problem 0.5 || A x
%                           - b || ^ 2 subject to || x || <= normConst.
%                           Structure: Vector (Column Vector).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  See - https://stats.stackexchange.com/questions/401212.
% Remarks:
%   1.  T
% TODO:
%   1.  U
% Release Notes:
%   -   1.0.000     14/04/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

paramLambda = 0; %<! Initialization

vX = pinv(mA) * vB;
mI = eye(size(vX, 1));

mAA = mA.' * mA;
mAb = mA.' * vB;

if((norm(vX, 2) ^ 2) <= normSquaredConst)
    return;
end

hObjFun     = @(paramLambda) (norm((mAA + (2 * paramLambda * mI)) \ mAb, 2) ^ 2) - normSquaredConst;
hObjFunGrad = @(paramLambda) -4 * vB.' * mA * inv(mAA + (2 * paramLambda * mI)) * inv(mAA + (2 * paramLambda * mI)) * inv(mAA + (2 * paramLambda * mI)) * mA.' * vB;

for ii = 1:numIterations
    paramLambda = paramLambda - (hObjFun(paramLambda) / hObjFunGrad(paramLambda));
end

% paramLambda = fzero(hObjFun, 1e-6); %<! Doesn't work for some reason
% paramLambda = FindZero(hObjFun, paramLambda, 1e6);

% hObjFun     = @(paramLambda) (((norm((mAA + (2 * paramLambda * mI)) \ mAb, 2) ^ 2) - normSquaredConst) ^ 2);
% % paramLambda = fminbnd(hObjFun, 0, 2000);
% paramLambda = fminsearch(hObjFun, paramLambda);

vX = (mAA + (2 * paramLambda * mI)) \ mAb;

end

author: RoyiAvital

repo: StackExchangeCodes