text:
function [c,r,k]=pdfmoments(t,m,b,a)
%PDFMOMENTS convert between central moments, raw moments and cumulants [C,R,K]=(T,M,B,A)
%  Inputs: t  text string containing:
%               'm','r','k'  Imput is central moments, raw moments, cumulants [default 'm']
%               'M','R','K'  Ouptut c is central moments, raw moments, cumulants [default 'M']
%          m    vector of input moments; m(r) is moment r, m(1) is always the mean
%          a,b  If input moments are for x, output moments are for a*x+b [defaulta=1, b=0]
%
% Outputs: c    central moments (or as determined by 'R' or 'K' options)
%          r    raw moments
%          k    cumulants
%
% (a) For all formats, the first element is the mean (i.e. not the first central moment or cumulant which equal zero).
% (b) pdfmoments('k',[m,s^2,0,0,0])=pdfmoments('k',[0,1,0,0,0],m,s) gives a Gaussian with mean m and std dev s

%   Copyright (c) 1998 Mike Brookes,  mike.brookes@ic.ac.uk
%      Version: $Id: pdfmoments.m 9678 2017-04-04 14:28:37Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent n0 bc mk fa
if isempty(n0)
    n0=3;  % order 3
    bc=[1 1 0 0; 1 2 1 0; 1 3 3 1]; % binomial coefficients
    mk={[1 1 1];[0 1 1 1]};  %  cumulant coefficients: one row per term, powers of k(2:i+1), coef k->m, coef m->k
    mn=[1 0; 1 1]; % [i,j] = number of terms in m(i+1) whose lowest moment is >= j+1
    fa=1; % factorial list
end
% check arguments
if nargin<4 || isempty(a)
    a=1;
end
if nargin<3 || isempty(b)
    b=0;
end
if isempty(t)
    t='';
end
n=length(m);                                    % number of moments required
if n>n0                                         % check if need to update coefficient arrays
    if fix(n/2-1)>length(fa)                    % we need factorials up to fix(n/2-1)
        fal=length(fa);
        fa(fix(n/2-1),1)=0;                     % enlarge factorial vector
        for i=fal+1:fix(n/2-1)
            fa(i)=i*fa(i-1);                    % create new factorials
        end
    end
    bc(n,n+1)=0;                                % enlarge binomial coefficient array
    mk{n-1,1}=[];                               % enlarge cumulant coefficient array
    mn(n-1,n-1)=0;                              % enlarge cumulant coefficient counts
    for i=n0+1:n
        bc(i,1:i+1)=[1 bc(i-1,1:i)+bc(i-1,2:i+1)];                      % update binomial coefficients
        j=fix((i+1)/2);                                                 % first coefficient row to sum
        nr=1+sum(mn(((j-1:i-3)+(n-1)*(i-j-2:-1:0))));                   % number of terms
        mki=zeros(nr,i+1);                                              % coefficient matrix
        ix=1;
        mki(1,i-1:i+1)=1;                                               % first term always has a coefficient of 1
        for r=j:i-2                                                     % previous coefficients to use
            nk=mn(r-1+(n-1)*(i-r-2));                                   % number of new coefficients for this value of r
            mkk=mk{r-1};                                                % old coefficients for this value of r
            mkik=mkk(1:nk,1:r-1);                                       % extract just the list of powers for each term
            mkik(:,i-r-1)=mkik(:,i-r-1)+1;                              % increment the power of moment  i-r
            mki(ix+1:ix+nk,1:r-1)=mkik;                                 % and save as new terms
            mki(ix+1:ix+nk,i)=mkk(1:nk,r)*bc(i,i-r+1)./mkik(:,i-r-1);   % calculate coefficient for r->m
            rho=sum(mkik,2)-1;                                          % rho is one less than the sum of the moment powers
            mki(ix+1:ix+nk,i+1)= mki(ix+1:ix+nk,i).*fa(rho).*(-1).^rho; % calculate coefficient for m->r
            ix=ix+nk;                                                   % update the number of terms so far
        end
        mki=sortrows(mki);                                              % sort according to the lowest moment that is used
        mn(i-1,1:i-1)=[nr sum(cumprod(mki(:,1:i-2)==0,2),1)];           % update count of terms with lowest moment >= j+1
        mk{i-1}=mki;                                                    % save in persistent cell array
    end
    n0=n;                                       % coefficients are now calculated up to order n
end
% apply scaling if input type is 'c' or 'k'
mu=a*m(1)+b; % calculate new mean
c=m; % initialize output shapes
r=m;
k=m;
m=m(:)'; % now force the input to be a row vector
if any(t=='k')
    tin=3; % set input type
    k(:)=k(:)'.*a.^(1:n);
    k(1)=0; % first cumulant is actually zero
elseif any(t=='r')
    tin=2;
else
    tin=1;
    c(:)=c(:)'.*a.^(1:n);
    c(1)=0;  % first cenral moment is actually zero
end
tout=[(~any(t=='K') && ~any(t=='R')) (nargout>=2 || any(t=='R')) (nargout>=3 || any(t=='K'))]; % outputs required
for il=1:2 % loop through conversion routines twice
    % first convert between moments
    if il==1 % convert unscaled R -> C
        v=[1 m.*a.^(1:n)];
        bb=b-mu;
        doit=tin==2 && (tout(1) || tout(3));
    else  % convert C -> R or unscaled R -> R
        if tin==2 % input type was 'r' (v is OK from previous iteration)
            bb=b;
        else % input type was 'c' or 'k'
            v=[1 c(:)'];
            bb=mu;
        end
        doit=tout(2); % convert if 'R' output required
    end
    if doit
        y=v(2:end);
        if bb~=0 % don't bother if the constant term is zero
            for i=1:n
                y(i)=polyval(bc(i,1:i+1).*v(1:i+1),bb);
            end
        end
        if il==1 % convert unscaled R -> C
            c(:)=y;
        else  % convert C -> R or unscaled R -> R
            r(:)=y;
        end
    end
    % now convert cumulants to/from moments
    if il==1 % convert K -> C
        x=k(:)';
        doit=tin==3 && (tout(1) || tout(2));
    else  % convert C -> K
        x=c(:)';
        doit=(tin<3) && tout(3);
    end
    if doit
        y=x;
        for i=4:n
            mki=mk{i-1}; % get coefficient matrix
            y(i)=mki(:,i-1+il)'*prod(repmat(x(2:i),size(mki,1),1).^mki(:,1:i-1),2); % calculate moment/cumulant (neat but not efficient)
        end
        if il==1 % converted K -> C
            c(:)=y;
        else % converted C -> K
            k(:)=y;
        end
    end
end
c(1)=mu; % restore the means
k(1)=mu;
if any(t=='R')
    c=r;
elseif any(t=='K')
    c=k;
end

author: jtkim-kaist

repo: Speech-enhancement


text:
function welldrawdown
% Well drawdown - comparison for confined / half-confined / unconfined
% aquifers
%    using analytocal solutions
%
%   $Ekkehard Holzbecher  $Date: 2006/04/30 $
%--------------------------------------------------------------------------
K = 1.1e-5;            % hydraulic conductivity
H = 10;                % depth
T = K*H;               % transmissivity
Q = 1.e-4;             % pumping rate
r0 = 0.1;              % well radius
s0 = -1.;               % drawdown in well
c = 0.8e8;             % resistance of half-permeable layer
R = 35;                % maximum radius

r = linspace (r0,R,100);
hc = s0 + (Q/(2*pi*T))*log(r/r0);
hu = -H + s0 + sqrt(H*H + (Q/(pi*K))*log(r/r0));
hh = -(Q/(2*pi*T))*besselk(0,r/sqrt(T*c));

plot (r,hc,r,hh,r,hu);
legend ('confined','half-confined','unconfined')
xlabel ('distance [m]'); ylabel ('drawdown (neg) [m]');
title('Groundwater Drawdown due to Pumping')

author: Sable

repo: mcbench-benchmarks


text:
function [ boolout ] = istiffcapella( filename )
%ISTIFFCAPELLA Test if data is Capella TIFF format
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

tags = read_tiff_tags(filename);
boolout = isfield(tags{1}, 'Software') && strncmpi(tags{1}.Software,'Capella',7);

end

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

author: ngageoint


text:
function varargout = dsift(varargin)
% VL_DSIFT  Dense SIFT
%   [FRAMES,DESCRS] = VL_DSIFT(I) extracts a dense set of SIFT
%   keypoints from image I. I must be of class SINGLE and grayscale.
%   FRAMES is a 2 x NUMKEYPOINTS, each colum storing the center (X,Y)
%   of a keypoint frame (all frames have the same scale and
%   orientation). DESCRS is a 128 x NUMKEYPOINTS matrix with one
%   descriptor per column, in the same format of VL_SIFT().
%
%   VL_DSIFT() does NOT compute a Gaussian scale space of the image
%   I. Instead, the image should be pre-smoothed at the desired scale
%   level, e.b. by using the VL_IMSMOOTH() function.
%
%   The scale of the extracted descriptors is controlled by the option
%   SIZE, i.e. the width in pixels of a spatial bin (recall that a
%   SIFT descriptor is a spatial histogram with 4 x 4 bins).
%
%   The sampling density is controlled by the option STEP, which is
%   the horizontal and vertical displacement of each feature cetner to
%   the next.
%
%   The sampled image area is controlled by the option BOUNDS,
%   defining a rectangle in which features are comptued. A descriptor
%   is included in the rectangle if all the centers of the spatial
%   bins are included. The upper-left descriptor is placed so that the
%   uppler-left spatial bin center is algined with the upper-left
%   corner of the rectangle.
%
%   By default, VL_DSIFT() computes features equivalent to
%   VL_SIFT(). However, the FAST option can be used to turn on an
%   variant of the descriptor (see VLFeat C API documentation for
%   further details) which, while not strictly equivalent, it is much
%   faster.
%
%   VL_DSIFT() accepts the following options:
%
%   Step:: [1]
%     Extracts a SIFT descriptor each STEP pixels.
%
%   Size:: [3]
%     A spatial bin covers SIZE pixels.
%
%   Bounds:: [whole image]
%     Specifies a rectangular area where descriptors should be
%     extracted. The format is [XMIN, YMIN, XMAX, YMAX]. If this
%     option is not specified, the entiere image is used.  The
%     bounding box is clipped to the image boundaries.
%
%   Norm::
%     If specified, adds to the FRAMES ouptut argument a third
%     row containint the descriptor norm, or engergy, before
%     contrast normalization. This information can be used to
%     suppress low contrast descriptors.
%
%   Fast::
%     If specified, use a piecewise-flat, rather than Gaussian,
%     windowing function. While this breaks exact SIFT equivalence,
%     in practice is much faster to compute.
%
%   FloatDescriptors::
%     If specified, the descriptor are returned in floating point
%     rather than integer format.
%
%   Verbose::
%     If specified, be verbose.
%
%   RELATION TO THE SIFT DETECTOR
%
%   In the standard SIFT detector/descriptor, implemented by
%   VL_SIFT(), the size of a spatial bin is related to the keypoint
%   scale by a multiplier, called magnification factor, and denoted
%   MAGNIF. Therefore, the keypoint scale corresponding to the
%   descriptors extracted by VL_DSIFT() is equal to SIZE /
%   MAGNIF. VL_DSIFT() does not use MAGNIF because, by using dense
%   sampling, it avoids detecting keypoints in the first plance.
%
%   VL_DSIFT() does not smooth the image as SIFT does. Therefore, in
%   order to obtain equivalent results, the image should be
%   pre-smoothed approriately. Recall that in SIFT, for a keypoint of
%   scale S, the image is pre-smoothed by a Gaussian of variance S.^2
%   - 1/4 (see VL_SIFT() and VLFeat C API documentation).
%
%   Example::
%     This example produces equivalent SIFT descriptors using
%     VL_DSIFT() and VL_SIFT():
%
%      binSize = 8 ;
%      magnif = 3 ;
%      Is = vl_imsmooth(I, sqrt((binSize/magnif)^2 - .25)) ;
%
%      [f, d] = vl_dsift(Is, 'size', binSize) ;
%      f(3,:) = binSize/magnif ;
%      f(4,:) = 0 ;
%      [f_, d_] = vl_sift(I, 'frames', f) ;
%
%   Remark::
%     The equivalence is never exact due to (i) boundary effects
%     and (ii) the fact that VL_SIFT() downsamples the image to save
%     computation. It is, however, usually very good.
%
%   Remark::
%     In categorization it is often useful to under-smooth the image,
%     comared to standard SIFT, in order to keep the gradients
%     sharp.
%
%   FURTHER DETAILS ON THE GEOMETRY
%
%   As mentioned, the VL_DSIFT() descriptors cover the bounding box
%   specified by BOUNDS = [XMIN YMIN XMAX YMAX]. Thus the top-left bin
%   of the top-left descriptor is placed at (XMIN, YMIN). The next
%   three bins to the right are at XMIN + SIZE, XMIN + 2*SIZE, XMIN +
%   3*SIZE. The X coordiante of the center of the first descriptor is
%   therefore at (XMIN + XMIN + 3*SIZE) / 2 = XMIN + 3/2 * SIZE.  For
%   instance, if XMIN = 1 and SIZE = 3 (default values), the X
%   coordinate of the center of the first descriptor is at 1 + 3/2 * 3
%   = 5.5. For the second descriptor immediately to its right this is
%   5.5 + STEP, and so on.
%
%   See also: VL_SIFT(), VL_HELP().
[varargout{1:nargout}] = vl_dsift(varargin{:});

author: jianxiongxiao

repo: ProfXkit


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Unified Solutions,
%           Solvability/Observability Analysis and Uncertainty Description
%           in A Globally Optimal Framework.
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com

function [XX, ff] = optimize_quat_cov(AA, bb, q, F, cov_left, epsX, epsQuat, solver, verbose)
X = sdpvar(4, 4);
res = AA * vec(X) - bb;
f = res.' * res;
cons = [
    X - epsX * eye(4) >= 0, ...
    vec(F * X * F.' - cov_left) == 0, ...
    q.' * X * q <= epsQuat
    ];
options = sdpsettings('solver', solver, 'verbose', verbose);
if(strcmp(solver, 'sdpa_gmp'))
    options.sdpa_gmp.epsilonDash = 1.0e-35;
    options.sdpa_gmp.precision = 250;
end
optimize(cons, f, options);
XX = value(X);
ff = value(f);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function obj = updateBoxplot(obj, boxIndex)

% y: ...[DONE]
% x0: ...[DONE]
% x: ...[DONE]
% name: ...[DONE]
% boxmean: ...[NOT SUPPORTED IN MATLAB]
% boxpoints: ...[NOT SUPPORTED IN MATLAB]
% jitter: ...[NOT SUPPORTED IN MATLAB]
% pointpos: ...[NOT SUPPORTED IN MATLAB]
% whiskerwidth: ........................[TODO]
% fillcolor: ...[DONE]
% opacity: ---[TODO]
% xaxis: ...[DONE]
% yaxis: ...[DONE]
% showlegend: ...[DONE]
% stream: ...[HANDLED BY PLOTLY STREAM]
% visible: ...[DONE]
% type: ...[DONE]

% MARKER
% color: ...[NA]
% width: ...[NA]
% dash: ...[NA]
% opacity: ...[NA]
% shape: ...[NA]
% smoothing: ...[NA]
% outliercolor: ...[NOT SUPPORTED IN MATLAB]
% outlierwidth: ...[NOT SUPPORTED IN MATLAB]

% LINE
% color: ...[DONE]
% width: ...[DONE]
% dash: ...[DONE]
% opacity: ---[TODO]
% shape: ...[DONE]
% smoothing: ...[NOT SUPPORTED IN MATLAB]

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

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(boxIndex).AssociatedAxis);

%-BOX DATA STRUCTURE-%
box_data = get(obj.State.Plot(boxIndex).Handle);

%-BOX CHILDREN-%
box_child = box_data.Children;

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

%-CONFIRM PROPER BOXPLOT STRUCTURE-%

% check for compact boxplot
isCompact = ~isempty(findobj(obj.State.Plot(boxIndex).Handle,'Tag','Whisker'));

% number of boxplots
if isCompact
    bpcompnum = 6;
    bpnum = length(box_child)/bpcompnum;
    % check for assumed box structure
    if mod(length(box_child), bpcompnum)~=0
        updateAlternativeBoxplot(obj, boxIndex);
        return
    end
else
    bpcompnum = 8;
    bpnum = length(box_child)/bpcompnum;
    % check for assumed box structure
    if mod(length(box_child),bpcompnum)~=0
        updateAlternativeBoxplot(obj, boxIndex);
        return
    end
end

% initialize ydata
ydata = [];

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

%-box groupgap-%
obj.layout.bargroupgap = 1/bpnum;

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

%-box name-%
obj.data{boxIndex}.name = box_data.DisplayName;

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

% iterate through box plot children in reverse order
for bp = bpnum:-1:1

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj,axIndex);

    %-AXIS DATA-%
    eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
    eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

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

    %-box xaxis-%
    obj.data{boxIndex}.xaxis = ['x' num2str(xsource)];

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

    %-box yaxis-%
    obj.data{boxIndex}.yaxis = ['y' num2str(ysource)];

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

    %-box type-%
    obj.data{boxIndex}.type = 'box';

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

    %-box visible-%
    obj.data{boxIndex}.visible = strcmp(box_data.Visible,'on');

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

    %-box fillcolor-%
    obj.data{boxIndex}.fillcolor = 'rgba(0, 0, 0, 0)';

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

    %-box showlegend-%
    leg = get(box_data.Annotation);
    legInfo = get(leg.LegendInformation);

    switch legInfo.IconDisplayStyle
        case 'on'
            showleg = true;
        case 'off'
            showleg = false;
    end

    obj.data{boxIndex}.showlegend = showleg;

    %-boxplot components-%
    Q1 = [];
    Q3 = [];
    median = [];
    outliers = [];
    uwhisker = [];
    lwhisker = [];

    % iterate through boxplot components
    for bpc = 1:bpcompnum

        %get box child data
        box_child_data = get(box_child(bp+bpnum*(bpc-1)));

        %box name
        if strcmp(box_child_data.Type,'text')
            if iscell(box_child_data.String)
                boxname =  box_child_data.String{1};
            else
                boxname =  box_child_data.String;
            end
        end

        % parse boxplot tags
        switch box_child_data.Tag

            %-median-%
            case 'Median'

                median =  box_child_data.YData(1);

                %-upper whisker-%
            case 'Upper Whisker'

                uwhisker = box_child_data.YData(2);

                %-boxplot whisker width-%
                obj.data{boxIndex}.whiskerwidth = 1;

                %-lower whisker-%
            case 'Lower Whisker'

                lwhisker = box_child_data.YData(1);

            case 'Box'

                %-Q1-%
                Q1 = min(box_child_data.YData);

                %-Q3-%
                Q3 = max(box_child_data.YData);

                %-boxplot line style-%

                if isCompact
                    col = 255*box_child_data.Color;
                    obj.data{boxIndex}.fillcolor = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];
                else
                    obj.data{boxIndex}.line = extractLineLine(box_child_data);
                end

                %-outliers-%
            case 'Outliers'

                if ~isnan(box_child_data.YData)
                    %-outlier marker data-%
                    outliers = box_child_data.YData;

                    %-outlier marker style-%
                    obj.data{boxIndex}.marker = extractLineMarker(box_child_data);
                end

                %-compact whiskers-%
            case 'Whisker'

                %-boxplot line style-%
                obj.data{boxIndex}.line = extractLineLine(box_child_data);

                %-boxplot whisker width-%
                obj.data{boxIndex}.whiskerwidth = 0;

                %-whisker data-%
                uwhisker = box_child_data.YData(2);
                lwhisker = box_child_data.YData(1);

                %-compact median-%
            case 'MedianInner'

                median = box_child_data.YData(1);

        end
    end

    %-generate boxplot data-%
    gendata = generateBoxData(outliers, lwhisker, Q1, median, Q3, uwhisker);

    %-boxplot y-data-%
    obj.data{boxIndex}.y(length(ydata)+1:length(ydata)+length(gendata)) = ...
        generateBoxData(outliers, lwhisker, Q1, median, Q3, uwhisker);

    %-boxplot x-data-%
    if (bpnum > 1)
        for n = (length(ydata)+1):(length(ydata)+length(gendata))
            obj.data{boxIndex}.x{n} = boxname;
        end
    end

    %-update ydata-%
    ydata = obj.data{boxIndex}.y;

end

%----------------------------!AXIS UPDATE!--------------------------------%

% take first text object as prototype for axis tick style/layout
text_child = findobj(obj.State.Plot(boxIndex).Handle,'Type','text');

%-STANDARDIZE UNITS-%
fontunits = get(text_child(1),'FontUnits');
set(text_child(1),'FontUnits','points');

%-text data -%
text_data = get(text_child(1));

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

%-xaxis tick font size-%
xaxis.tickfont.size = text_data.FontSize;

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

%-xaxis tick font family-%
xaxis.tickfont.family = matlab2plotlyfont(text_data.FontName);

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

%-xaxis tick font color-%
xaxis.tickfont.color = text_data.Color;

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

%-axis type-%
xaxis.type = 'category';

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

%-show tick labels-%
xaxis.showticklabels = true;

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

%-autorange-%
xaxis.autorange = true;

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

%-set the layout axis field-%
obj.layout = setfield(obj.layout,['xaxis' num2str(xsource)],xaxis);

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

%-REVERT UNITS-%
set(text_child(1),'FontUnits',fontunits);

end

author: plotly


text:
function [ point_num2, edge_num2, face_num2, face_order_max2 ] = ...
  dual_size_3d ( point_num, edge_num, face_num, face_order_max, ...
  point_coord, face_order, face_point )

%*****************************************************************************80
%
%% DUAL_SIZE_3D determines sizes for a dual of a shape in 3D.
%
%  Discussion:
%
%    We don't actually need FACE_POINT as input here.  But since the
%    three arrays occur together everywhere else, it seems unnecessarily
%    user-confusing to vary the usage here!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer EDGE_NUM, the number of edges.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer FACE_ORDER_MAX, the maximum number of vertices per face.
%
%    Input, real POINT_COORD(3,POINT_NUM); POINT_COORD(*,J) is
%    the X, Y and Z coordinates of point J.
%
%    Input, integer FACE_ORDER(FACE_NUM), the number of vertices per face.
%
%    Input, integer FACE_POINT(FACE_ORDER_MAX,FACE_NUM); FACE_POINT(I,J)
%    is the index of the I-th point in the J-th face.  The
%    points are listed in the counter-clockwise direction defined
%    by the outward normal at the face.
%
%    Output, integer POINT_NUM2, the number of points in the dual.
%
%    Output, integer EDGE_NUM2, the number of edges in the dual.
%
%    Output, integer FACE_NUM2, the number of faces in the dual.
%
%    Output, integer FACE_ORDER_MAX2, the maximum number of vertices per face
%    in the dual.
%

%
%  These values are easy to compute:
%
  point_num2 = face_num;
  edge_num2 = edge_num;
  face_num2 = point_num;
%
%  To determine FACE_ORDER_MAX2 is not so easy.
%  You have to construct the FACE_ORDER array for the dual shape.
%  The order of a dual face is the number of edges that the vertex occurs in.
%  But then all we have to do is count how many times each item shows up
%  in the FACE_POINT array.
%
  face_order_max2 = 0;
  face_order2(1:face_num2) = 0;

  for face = 1 : face_num
    for i = 1 : face_order(face)
      face2 = face_point(i,face);
      face_order2(face2) = face_order2(face2) + 1;
    end
  end

  face_order_max2 = max ( face_order2(1:face_num2) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestGetPerspectiveTransform
    %TestGetPerspectiveTransform

    methods (Static)
        function test_1
            ref = [0, 1, 0; -1, 0, 1; 0, 0, 1];
            src = [0.0, 1.0; 1.0, 1.0; 1.0, 0.0; 0.25, 0.25];
            dst = [1.0, 1.0; 1.0, 0.0; 0.0, 0.0; 0.25, 0.75];
            t = cv.getPerspectiveTransform(src, dst);
            validateattributes(t, {'numeric'}, {'size',[3 3]});
            assert(all(abs(t(:)-ref(:)) < 1e-10));
        end

        function test_2
            ref = [0, 1, 0; -1, 0, 1; 0, 0, 1];
            src = num2cell([0.0, 1.0; 1.0, 1.0; 1.0, 0.0; 0.25, 0.25], 2);
            dst = num2cell([1.0, 1.0; 1.0, 0.0; 0.0, 0.0; 0.25, 0.75], 2);
            t = cv.getPerspectiveTransform(src, dst);
            validateattributes(t, {'numeric'}, {'size',[3 3]});
            assert(all(abs(t(:)-ref(:)) < 1e-10));
        end

        function test_error_argnum
            try
                cv.getPerspectiveTransform();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function s = mean(f, dim)
%MEAN   Average or mean value of a SPHEREFUN.
%   MEAN(F) takes the mean in the latitude-direction (default), i.e.,
%          MEAN(F) = 1/pi sum(F).
%
%   MEAN(F, DIM) takes the mean along the direction DIM. If DIM = 1 it is the
%   latitude-direction and if DIM = 2 then it is the longitude-direction.
%
% See also SPHEREFUN/MEAN2, SPHEREFUN/STD2.

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

% Empty check:
if ( isempty( f ) )
    s = chebfun;
    return
end

if ( nargin == 1)
    % Default to the y-direction:
    dim = 1;
end
dom = f.domain;

s = sum( f, dim );
if ( dim == 1 )
    s = s / pi; % Mean in the latitude direction (default)
elseif ( dim == 2 )
    s = s / (2*pi); % Mean in the longitude direction
else
    error('CHEBFUN:SPHEREFUN:mean:dim', ...
        'Mean not in longitude (LAMBDA) or latitude (THETA) direction.')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function view = getFlatCoords(view)
%
% view = getFlatCoords(view)
%
% view: must be a flat view
%
% Loads gLocs2d and gLocs3d coordinates.  Keeps only those voxels
% that correspond to the inplane coordinates.  Sets FLAT
% fields: coordsRight, coordsLeft, grayCoordsRight,
% grayCoordsLeft.
%
% djh, 7/98
%
% djh, 8/4/99.  Use intersectCols instead of intersecting the indices.
% djh, 8/18/99.  loadGLocs now returns curvature.

if ~strcmp(view.viewType,'Flat')
    myErrorDlg('function getFlatCoords only for Flat view.');
end

pathStr=fullfile(viewDir(view),'coords');

if ~check4File(pathStr)
    imSize=[0,0];
    waitHandle = mrvWaitbar(0,'Computing flat coordinates.  Please wait...');
    for h = 1:2
        mrvWaitbar((h-1)/2)

        % Load gLocs2d and gLocs3d
        %
        if h==1
            [gLocs2d,gLocs3d,curvature,leftPath] = loadGLocs('left');
        else
            [gLocs2d,gLocs3d,curvature,rightPath] = loadGLocs('right');
        end

        if isempty(gLocs2d) | isempty(gLocs3d)
            coords{h} = [];
            grayCoords{h} = [];
        else
            % Compute imSize
            %
            imSize = max(imSize,(max(gLocs2d,[],2) - min(gLocs2d,[],2) + 1)');
            imSize = round(imSize);

            % Find gray nodes that are both in the inplanes and included
            % in the unfold.
            % gray.coords are the gray coords that lie in the inplanes.
            % gLocs3d are the gray coords in the unfold.
            hiddenGray = initHiddenGray;
            [grayCoordsTmp,gLocsIndices,coordsIndices] = ...
                intersectCols(gLocs3d,hiddenGray.coords);
            grayCoords{h} = grayCoordsTmp;

            % Flat locations corresponding to those voxels
            %
            coords{h} = gLocs2d(:,gLocsIndices);

            % Warning if there are any NaNs in the coords.
            NaNs = find(isnan(coords{h}(1,:)) | isnan(coords{h}(2,:)));
            if ~isempty(NaNs)
                myWarnDlg(['You have ',int2str(length(NaNs)),' NaNs in your flat coords.  ',...
                        'Those gray matter nodes will not be rendered in the FLAT view.']);
            end
        end
    end
    close(waitHandle)

    % Save to file
    %
    save(pathStr,'coords','grayCoords','imSize','leftPath','rightPath');
end

% Load Flat/coords and fill the fields
%
load(pathStr);
view.coords = coords;
view.grayCoords = grayCoords;
view.leftPath = leftPath;
view.rightPath = rightPath;
view.ui.imSize = imSize;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function rhs_lucas_test ( )

%*****************************************************************************80
%
%% RHS_LUCAS_TEST samples the right hand side at the initial time.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/navier_stokes_2d_exact/rhs_lucas_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2015
%
%  Author:
%
%    John Burkardt
%
  nu = 1.0;
  rho = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'RHS_LUCAS_TEST\n' );
  fprintf ( 1, '  Lucas Bystricky Flow\n' );
  fprintf ( 1, '  Sample the Navier-Stokes right hand sides\n' );
  fprintf ( 1, '  at the initial time T = 0, using the unit square.\n' );
  fprintf ( 1, '  Kinematic viscosity NU = %g\n', nu );
  fprintf ( 1, '  Fluid density RHO = %g\n', rho );

  n = 1000;
  xy_lo = 0.0;
  xy_hi = 1.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, xy_lo, xy_hi, seed );
  [ y, seed ] = r8vec_uniform_ab ( n, xy_lo, xy_hi, seed );
  t = 0.0;

  [ f, g, h ] = rhs_lucas ( nu, rho, n, x, y, t );

  fprintf ( 1, '\n' );
  fprintf ( 1, '           Minimum       Maximum\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F:  %14.6g  %14.6g\n', min ( f ), max ( f ) );
  fprintf ( 1, '  G:  %14.6g  %14.6g\n', min ( g ), max ( g ) );
  fprintf ( 1, '  H:  %14.6g  %14.6g\n', min ( h ), max ( h ) );

  return
end


text:
function [bnet, Unode, Snode, Lnodes, Rnode, Ynode, Lsnode] = ...
    mk_gmux_robot_dbn(nlandmarks, Q, R, init_x, init_V, robot_block, landmark_block)

% Make DBN

% S
% | L1 -------> L1'
% |  | L2 ----------> L2'
% \  | /
%  v v v
%    Ls
%    |
%    v
%    Y
%    ^
%    |
%    R ------->  R'
%    ^
%    |
%    U
%
%
% S is a switch, Ls is a deterministic gmux, Y = Ls-R,
% R(t+1) = R(t) + U(t+1), L(t+1) = L(t)

% number nodes topologically
Snode = 1;
Lnodes = 2:nlandmarks+1;
Lsnode = nlandmarks+2;
Unode = nlandmarks+3;
Rnode = nlandmarks+4;
Ynode = nlandmarks+5;

nnodes = nlandmarks+5;
intra = zeros(nnodes, nnodes);
intra([Snode Lnodes], Lsnode) =1;
intra(Unode,Rnode)=1;
intra([Rnode Lsnode], Ynode)=1;

inter = zeros(nnodes, nnodes);
inter(Rnode, Rnode)=1;
for i=1:nlandmarks
  inter(Lnodes(i), Lnodes(i))=1;
end

Lsz = 2; % (x y) posn of landmark
Rsz = 2; % (x y) posn of robot
Ysz = 2; % relative distance
Usz = 2; % (dx dy) ctrl
Ssz = nlandmarks; % can switch between any landmark

ns = zeros(1,nnodes);
ns(Snode) = Ssz;
ns(Lnodes) = Lsz;
ns(Lsnode) = Lsz;
ns(Ynode) = Ysz;
ns(Rnode) = Rsz;
ns(Ynode) = Usz;
ns(Unode) = Usz;

bnet = mk_dbn(intra, inter, ns, 'discrete', Snode, 'observed', [Snode Ynode Unode]);

bnet.CPD{Snode} = root_CPD(bnet, Snode); % always observed
bnet.CPD{Unode} = root_CPD(bnet, Unode); % always observed
for i=1:nlandmarks
  bi = landmark_block(:,i);
  bnet.CPD{Lnodes(i)} = gaussian_CPD(bnet, Lnodes(i), 'mean', init_x(bi), 'cov', init_V(bi,bi));
end
bi = robot_block;
bnet.CPD{Rnode} = gaussian_CPD(bnet, Rnode, 'mean', init_x(bi), 'cov', init_V(bi,bi), 'weights', eye(2));
bnet.CPD{Lsnode} = gmux_CPD(bnet, Lsnode, 'cov', repmat(zeros(Lsz,Lsz), [1 1 nlandmarks]), ...
			    'weights', repmat(eye(Lsz,Lsz), [1 1 nlandmarks]));
W = [eye(2) -eye(2)]; % Y = Ls - R, where Ls is the lower-numbered parent
bnet.CPD{Ynode} = gaussian_CPD(bnet, Ynode, 'mean', zeros(Ysz,1), 'cov', R, 'weights', W);

% slice 2
eclass = bnet.equiv_class;
W = [eye(2) eye(2)]; % R(t) = R(t-1) + U(t), where R(t-1) is the lower-numbered parent
bnet.CPD{eclass(Rnode,2)} = gaussian_CPD(bnet, Rnode+nnodes, 'mean', zeros(Rsz,1), 'cov', Q, 'weights', W);
for i=1:nlandmarks
  bnet.CPD{eclass(Lnodes(i), 2)} = gaussian_CPD(bnet, Lnodes(i)+nnodes, 'mean', zeros(2,1), ...
						   'cov', zeros(2,2), 'weights', eye(2));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
close all
clear all
clc

%Specify start and end points
start_point = [100;40]
end_point = [55;40]

%Specify external boundaries
external_boundaries = [0,0;60,0;60,45;45,45;45,59;75.5,40;106,59;106,45;91,45;91,0;151,0;151,105;50,105;0,60];

%Let the external boundaries form a closed polygon for graphical
%visualization purposes
external_boundaries_draw = external_boundaries;
external_boundaries_draw(size(external_boundaries,1)+1,:) = external_boundaries_draw(1,:);

hold on
axis equal

%Plot the boundaries, starting points and end points
plot(external_boundaries_draw(:,1), external_boundaries_draw(:,2), 'Color', 'black')
plot(start_point(1), start_point(2), 'X', 'Color', 'green')
plot(start_point(1), start_point(2), 'O', 'Color', 'green')
plot(end_point(1), end_point(2), 'X', 'Color', 'red')
plot(end_point(1), end_point(2), 'O', 'Color', 'red')

tic
waypoint_coordinates = pathfinder(start_point, end_point, external_boundaries)
toc

%Plot the chosen path
plot(waypoint_coordinates(:,1), waypoint_coordinates(:,2), 'Color', 'blue')


text:
function l4_uniform_test ( )

%*****************************************************************************80
%
%% L4_UNIFORM_TEST tests L4_UNIFORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2008
%
%  Author:
%
%    John Burkardt
%

  fprintf ( 1, '\n' );
  fprintf ( 1, 'L4_UNIFORM_TEST\n' );
  fprintf ( 1, '  L4_UNIFORM computes pseudorandom logical values.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ value, seed ] = l4_uniform ( seed );
    fprintf ( 1, '  %8d  %d\n', i, value );
  end

  return
end


text:
%% Script to plot MOP1
% This scripts plot the objective space for MOP1 given in Homework
% Assignment # 5. It used the function eveluate_objective to determine the
% objective values and plots the two dimensional objective space.
hold on
plot3(0,0,0);
axis([0 2 0 2 0 2]);
for i = 1 : 10
    s = sprintf('%d',i);
    %plot3(x(i,13),x(i,14),x(i,15),'String',i);
    %text(x(i,13)+.05,x(i,14)+.05,x(i,14)+.05,s);
    text(x(i,13),x(i,14),x(i,14),s);
end
hold off


text:
% Mathematics Q1344369
% https://math.stackexchange.com/questions/1344369
% Compressive Sensing over the Complex Domain
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     14/04/2018
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

numRows = 64;
numCols = 256;

paramLambda = 4.0;

numIterations = 1000;

cSolversFun = {@(mA, vB, paramLambda, numIterations) SolveLsL1ComplexSubGrad(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexRealSubGrad(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexPgm(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexRealPgm(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexAdmm(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexIrls(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexCd(mA, vB, paramLambda, numIterations); ...
    @(mA, vB, paramLambda, numIterations) SolveLsL1ComplexRealCd(mA, vB, paramLambda, numIterations);};

cMethodString = {['Sub Gradient Method']; ['Sub Gradient Method - Real Domain']; ...
    ['Proximal Gradient Method']; ['Proximal Gradient Method - Real Domain']; ...
    ['ADMM Method']; ['Fixed Point Iteration (IRLS) Method']; ...
    ['Coordinate Descent (CD) Method']; ['Coordinate Descent (CD) Method - Real Domain']};

%% Generate Data

mA = randn([numRows, numCols]) + (1i * randn([numRows, numCols]));
vB = randn([numRows, 1]) + (1i * randn([numRows, 1]));

hCalcErrorNorm = @(mX, vXRef) sum((mX - vXRef) .* conj(mX - vXRef));
hCalcObjFunVal = @(mX) 0.5 * sum(((mA * mX) - vB) .* conj((mA * mX) - vB)) + paramLambda * sum(abs(mX), 1);

%% Solution by CVX

tic();
cvx_begin quiet
    variable vXCvx(numCols) complex
    minimize( (0.5 * sum_square_abs((mA * vXCvx) - vB)) + paramLambda * norm(vXCvx, 1) )
cvx_end
toc();

% disp([' ']);
% disp(['CVX Solution Summary']);
% disp(['The CVX Solver Status - ', cvx_status]);
% disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
% disp(['The Optimal Argument Is Given By - [ ', num2str(vXCvx.'), ' ]']);
% disp([' ']);

%% Method Analysis

numMethods = length(cSolversFun);

vRunTime = zeros([numMethods, 1]);
mErrorNorm = zeros([numIterations, numMethods]);
mObjFunVal = zeros([numIterations, numMethods]);

for ii = 1:numMethods
    hRunTime = tic();
    [vX, mX] = cSolversFun{ii}(mA, vB, paramLambda, numIterations);
    runTime = toc(hRunTime);

    vRunTime(ii)         = runTime;
    mErrorNorm(:, ii)    = hCalcErrorNorm(mX, vXCvx);
    mObjFunVal(:, ii)    = hCalcObjFunVal(mX);
end

%% Display Results

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries = plot([1:numIterations], 10 * log10(mErrorNorm));
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', ['The Error Norm - $ {\left\| {x}^{k} - {x}_{CVX} \right\|}_{2}^{2} $'], ...
    'FontSize', fontSizeTitle, 'Interpreter', 'latex');
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Error Norm [dB]', ...
    'FontSize', fontSizeAxis);
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);
hLegend = ClickableLegend(cMethodString);
set(hLegend, 'FontSize', fontSizeAxis);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries = plot([1:numIterations], mObjFunVal);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(5), 'LineStyle', ':', 'LineWidth', lineWidthThin);
set(get(hAxes, 'Title'), 'String', ['The Objection Function Value - $ \frac{1}{2} {\left\| A {x}^{k} - b \right\|}_{2}^{2} + \lambda {\left\| {x}^{k} \right\|}_{1} $'], ...
    'FontSize', fontSizeTitle, 'Interpreter', 'latex');
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);
hLegend = ClickableLegend(cMethodString);
set(hLegend, 'FontSize', fontSizeAxis);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hBarObj = bar(1:length(vRunTime), vRunTime);
% set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2:end), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', ['L_1 Regularized Least Squares - Methods Run Time'], ...
    'FontSize', fontSizeTitle);
set(hAxes, 'XTickLabel', cMethodString, 'XTickLabelRotation', 45);
set(get(hAxes, 'XLabel'), 'String', 'Method', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Run Time [Sec]', ...
    'FontSize', fontSizeAxis);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function anonStructuresS = defineAnonStructures()
%
% function anonStructuresS = defineAnonStructures()
%
% Defines the anonymous "structures" data structure. The fields containing PHI have
% been left out by default. In order to leave out additional fields, remove
% them from the list below.
%
% The following three operations are allowed per field of the input data structure:
% 1. Keep the value as is: listed as 'keep' in the list below.
% 2. Allow only the specific values: permitted values listed within a cell array.
%    If no match is found, a default anonymous string will be inserted.
% 3. Insert dummy date: listed as 'date' in the list below. Date will be
%    replaced by a dummy value of 11/11/1111.
%
% APA, 1/11/2018

segmentAnonS = struct('points','keep');
contourAnonS = struct('segments',segmentAnonS);

anonStructuresS = struct( ...
    'imageType', {{'STRUCTURE'}}, ...
    'caseNumber', 'keep', ...
    'roiNumber', 'keep', ...
    'structureName', 'keep', ...
    'numberRepresentation', {{'CHARACTER'}}, ...
    'structureFormat', {{'SCAN-BASED'}}, ...
    'numberOfScans', 'keep', ...
    'maximumNumberScans', 'keep', ...
    'maximumPointsPerSegment', 'keep', ...
    'maximumSegmentsPerScan', 'keep', ...
    'unitNumber', 'keep', ...
    'writer', {{'MIM Software Inc.'}}, ...
    'structureColor', 'keep', ...
    'contour', contourAnonS, ...
    'rasterSegments', 'keep', ...
    'DSHPoints', 'keep', ...
    'orientationOfStructure', {{'TRANSVERSE'}}, ...
    'transferProtocol', {{'DICOM'}}, ...
    'visible', 'keep', ...
    'associatedScan', 'keep', ...
    'strUID', 'keep', ...
    'assocScanUID', 'keep', ...
    'rasterized', 'keep' ...
    );

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function kronrod_test02 ( )

%*****************************************************************************80
%
%% KRONROD_TEST02 tests the code for the even case N = 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 August 2010
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1,  'KRONROD_TEST02\n' );
  fprintf ( 1,  '  Request KRONROD to compute the Gauss rule\n' );
  fprintf ( 1,  '  of order 4, and the Kronrod extension of\n' );
  fprintf ( 1,  '  order 4+5=9.\n' );

  tol = 0.000001;

  [ x, w1, w2 ] = kronrod ( n, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  KRONROD returns 3 vectors of length %d\n', n + 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1,  '     I      X               WK              WG\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n + 1
    fprintf ( 1, '  %4d  %14f  %14f  %14f\n', i, x(i), w1(i), w2(i) );
  end

  return
end


text:
function [qrs_pos,sign,en_thres] = qrs_detect2(ecg,varargin)
% QRS detector based on the P&T method. This is an offline implementation
% of the detector.
%
% inputs
%   ecg:            one ecg channel on which to run the detector (required)
%                   in [mV]
%   varargin
%       THRES:      energy threshold of the detector (default: 0.6)
%                   [arbitrary units]
%       REF_PERIOD: refractory period in sec between two R-peaks (default: 0.250)
%                   in [ms]
%       fs:         sampling frequency (default: 1KHz) [Hz]
%       fid_vec:    if some subsegments should not be used for finding the
%                   optimal threshold of the P&Tthen input the indices of
%                   the corresponding points here
%       SIGN_FORCE: force sign of peaks (positive value/negative value).
%                   Particularly usefull if we do window by window detection and want to
%                   unsure the sign of the peaks to be the same accross
%                   windows (which is necessary to build an FECG template)
%       debug:      1: plot to bebug, 0: do not plot
%
% outputs
%   qrs_pos:        indexes of detected peaks (in samples)
%   sign:           sign of the peaks (a pos or neg number)
%   en_thres:       energy threshold used
%
%
%
% Physionet Challenge 2014, version 1.0
% Released under the GNU General Public License
%
% Copyright (C) 2014  Joachim Behar
% Oxford university, Intelligent Patient Monitoring Group
% joachim.behar@gmail.com
%
% Last updates:
%
% 13-09-2014
% Joachim Behar
% - bug on refrac period fixed
% - code a bit more tidy
% - condition added on flatline detection for overall segment (if flatline
% then returns empty matrices rather than some random stuff)
%
% 15-09-2014
% Julien Oster
% Sombrero hat instead of band-pass filter as the prefiltering step.
%
%
% When using this code, please reference the following papers:
%
% [1] Behar Joachim, Jonhson Alistair, Clifford Gari D., Oster Julien A
% Comparison of Single Channel Foetal ECG Extraction Methods.
% Annals of Biomedical Engineering. 42(6), 1340-53. 2014
%
% [2] Johnson Alistair E W, Behar Joachim, Oster Julien and Clifford Gari
% D. R-Peak Estimation using Multimodal Lead Switching.
% Accepted for Computing in Cardiology conference 2014.
%
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

% == managing inputs
REF_PERIOD = 0.250;
THRES = 0.6;
fs = 1000;
fid_vec = [];
SIGN_FORCE = [];
debug = 0;

switch nargin
    case 2
        REF_PERIOD=varargin{1};
    case 3
        REF_PERIOD=varargin{1};
        THRES=varargin{2};
    case 4
        REF_PERIOD=varargin{1};
        THRES=varargin{2};
        fs=varargin{3};
    case 5
        REF_PERIOD=varargin{1};
        THRES=varargin{2};
        fs=varargin{3};
        fid_vec=varargin{4};
    case 6
        REF_PERIOD=varargin{1};
        THRES=varargin{2};
        fs=varargin{3};
        fid_vec=varargin{4};
        SIGN_FORCE=varargin{5};
    case 7
        REF_PERIOD=varargin{1};
        THRES=varargin{2};
        fs=varargin{3};
        fid_vec=varargin{4};
        SIGN_FORCE=varargin{5};
        debug=varargin{6};
    otherwise
        error('qrs_detect: wrong number of input arguments \n');
end

[a b] = size(ecg);
if(a>b); NB_SAMP=a; elseif(b>a); NB_SAMP=b; ecg=ecg'; end;
tm = 1/fs:1/fs:ceil(NB_SAMP/fs);

% == constants
MED_SMOOTH_NB_COEFF = round(fs/100);
INT_NB_COEFF = round(7*fs/256); % length is 7 for fs=256Hz
SEARCH_BACK = 1; % perform search back (FIXME: should be in function param)
MAX_FORCE = []; % if you want to force the energy threshold value (FIXME: should be in function param)
MIN_AMP = 0.1; % if the median of the filtered ECG is inferior to MINAMP then it is likely to be a flatline
               % note the importance of the units here for the ECG (mV)
NB_SAMP = length(ecg); % number of input samples

try
    % == Bandpass filtering for ECG signal
    % this sombrero hat has shown to give slightly better results than a
    % standard band-pass filter. Plot the frequency response to convince
    % yourself of what it does
    b1 = [-7.757327341237223e-05  -2.357742589814283e-04 -6.689305101192819e-04 -0.001770119249103 ...
         -0.004364327211358 -0.010013251577232 -0.021344241245400 -0.042182820580118 -0.077080889653194...
         -0.129740392318591 -0.200064921294891 -0.280328573340852 -0.352139052257134 -0.386867664739069 ...
         -0.351974030208595 -0.223363323458050 0 0.286427448595213 0.574058766243311 ...
         0.788100265785590 0.867325070584078 0.788100265785590 0.574058766243311 0.286427448595213 0 ...
         -0.223363323458050 -0.351974030208595 -0.386867664739069 -0.352139052257134...
         -0.280328573340852 -0.200064921294891 -0.129740392318591 -0.077080889653194 -0.042182820580118 ...
         -0.021344241245400 -0.010013251577232 -0.004364327211358 -0.001770119249103 -6.689305101192819e-04...
         -2.357742589814283e-04 -7.757327341237223e-05];

    b1 = resample(b1,fs,250);
    bpfecg = filtfilt(b1,1,ecg)';

    if (length(find(abs(bpfecg)>MIN_AMP))/NB_SAMP)>0.20
        % if 20% of the samples have an absolute amplitude which is higher
        % than MIN_AMP then we are good to go.

        % == P&T operations
        dffecg = diff(bpfecg');  % (4) differentiate (one datum shorter)
        sqrecg = dffecg.*dffecg; % (5) square ecg
        intecg = filter(ones(1,INT_NB_COEFF),1,sqrecg); % (6) integrate
        mdfint = medfilt1(intecg,MED_SMOOTH_NB_COEFF);  % (7) smooth
        delay  = ceil(INT_NB_COEFF/2);
        mdfint = circshift(mdfint,-delay); % remove filter delay for scanning back through ECG

        % look for some measure of signal quality with signal fid_vec? (FIXME)
        if isempty(fid_vec); mdfintFidel = mdfint; else mdfintFidel(fid_vec>2) = 0; end;

        % == P&T threshold
        if NB_SAMP/fs>90; xs=sort(mdfintFidel(fs:fs*90)); else xs = sort(mdfintFidel(fs:end)); end;

        if isempty(MAX_FORCE)
           if NB_SAMP/fs>10
                ind_xs = ceil(98/100*length(xs));
                en_thres = xs(ind_xs); % if more than ten seconds of ecg then 98% CI
            else
                ind_xs = ceil(99/100*length(xs));
                en_thres = xs(ind_xs); % else 99% CI
            end
        else
           en_thres = MAX_FORCE;
        end

        % build an array of segments to look into
        poss_reg = mdfint>(THRES*en_thres);

        % in case empty because force threshold and crap in the signal
        if isempty(poss_reg); poss_reg(10) = 1; end;

        % == P&T QRS detection & search back
        if SEARCH_BACK
            indAboveThreshold = find(poss_reg); % ind of samples above threshold
            RRv = diff(tm(indAboveThreshold));  % compute RRv
            medRRv = median(RRv(RRv>0.01));
            indMissedBeat = find(RRv>1.5*medRRv); % missed a peak?
            % find interval onto which a beat might have been missed
            indStart = indAboveThreshold(indMissedBeat);
            indEnd = indAboveThreshold(indMissedBeat+1);

            for i=1:length(indStart)
                % look for a peak on this interval by lowering the energy threshold
                poss_reg(indStart(i):indEnd(i)) = mdfint(indStart(i):indEnd(i))>(0.5*THRES*en_thres);
            end
        end

        % find indices into boudaries of each segment
        left  = find(diff([0 poss_reg'])==1);  % remember to zero pad at start
        right = find(diff([poss_reg' 0])==-1); % remember to zero pad at end

        % looking for max/min?
        if SIGN_FORCE
            sign = SIGN_FORCE;
        else
            nb_s = length(left<30*fs);
            loc  = zeros(1,nb_s);
            for j=1:nb_s
                [~,loc(j)] = max(abs(bpfecg(left(j):right(j))));
                loc(j) = loc(j)-1+left(j);
            end
            sign = mean(ecg(loc));  % FIXME: change to median?
        end

        % loop through all possibilities
        compt=1;
        NB_PEAKS = length(left);
        maxval = zeros(1,NB_PEAKS);
        maxloc = zeros(1,NB_PEAKS);
%%%%%%%% add two lines
        max_inv_val = zeros(1,NB_PEAKS);
        max_inv_loc = zeros(1,NB_PEAKS);
        for i=1:NB_PEAKS
            if sign>0
                % if sign is positive then look for positive peaks
                [maxval(compt) maxloc(compt)] = max(ecg(left(i):right(i)));
%%%%%%%%%%%%%%%% add 1 line
                [max_inv_val(compt) max_inv_loc(compt)] = min(ecg(left(i):right(i)));
            else
                % if sign is negative then look for negative peaks
                [maxval(compt) maxloc(compt)] = min(ecg(left(i):right(i)));
%%%%%%%%%%%%%%%% add 1 line
                [max_inv_val(compt) max_inv_loc(compt)] = max(ecg(left(i):right(i)));
            end
            maxloc(compt) = maxloc(compt)-1+left(i); % add offset of present location
%%%%%%%%%%%% add 1 line
            max_inv_loc(compt) = max_inv_loc(compt)-1+left(i); % add offset of present location

            % refractory period - has proved to improve results
            if compt>1
%%%%%%%%%%%%%%%%                if maxloc(compt)-maxloc(compt-1)<fs*REF_PERIOD && abs(maxval(compt))<abs(maxval(compt-1))
                if maxloc(compt)-maxloc(compt-1)<fs*REF_PERIOD && max(abs(maxval(compt)),abs(max_inv_val(compt)))<max(abs(maxval(compt-1)),abs(max_inv_val(compt-1)))
%%%%%%%%%%%%%%%% add 4 lines and modify the next one
                    if abs(maxval(compt-1))<abs(max_inv_val(compt-1))
                        maxloc(compt-1)=max_inv_loc(compt-1);
                        maxval(compt-1)=max_inv_val(compt-1);
                    end
                    maxloc(compt)=[]; maxval(compt)=[];max_inv_loc(compt)=[]; max_inv_val(compt)=[];
%%%%%%%%%%%%%%%%                elseif maxloc(compt)-maxloc(compt-1)<fs*REF_PERIOD && abs(maxval(compt))>=abs(maxval(compt-1))
                elseif maxloc(compt)-maxloc(compt-1)<fs*REF_PERIOD && max(abs(maxval(compt)),abs(max_inv_val(compt)))>=max(abs(maxval(compt-1)),abs(max_inv_val(compt-1)))
%%%%%%%%%%%%%%%% add 4 lines and modify the next one
                    if abs(maxval(compt))<abs(max_inv_val(compt))
                        maxloc(compt)=max_inv_loc(compt);
                        maxval(compt)=max_inv_val(compt);
                    end
                    maxloc(compt-1)=[]; maxval(compt-1)=[];max_inv_loc(compt-1)=[]; max_inv_val(compt-1)=[];
                else
                    compt=compt+1;
                end
            else
                % if first peak then increment
                compt=compt+1;
            end
        end

        qrs_pos = maxloc; % datapoints QRS positions
        R_t = tm(maxloc); % timestamps QRS positions
        R_amp = maxval; % amplitude at QRS positions
        hrv = 60./diff(R_t); % heart rate
    else
        % this is a flat line
        qrs_pos = [];
        R_t = [];
        R_amp = [];
        hrv = [];
        sign = [];
        en_thres = [];
    end
catch ME
    rethrow(ME);
    for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
    qrs_pos = [1 10 20]; sign = 1; en_thres = 0.5;
end

% == plots
if debug
    figure;
    FONTSIZE = 20;
    ax(1) = subplot(4,1,1); plot(tm,ecg); hold on;plot(tm,bpfecg,'r')
        title('raw ECG (blue) and zero-pahse FIR filtered ECG (red)'); ylabel('ECG');
        xlim([0 tm(end)]);  hold off;
    ax(2) = subplot(4,1,2); plot(tm(1:length(mdfint)),mdfint);hold on;
        plot(tm,max(mdfint)*bpfecg/(2*max(bpfecg)),'r',tm(left),mdfint(left),'og',tm(right),mdfint(right),'om');
        title('Integrated ecg with scan boundaries over scaled ECG');
        ylabel('Int ECG'); xlim([0 tm(end)]); hold off;
    ax(3) = subplot(4,1,3); plot(tm,bpfecg,'r');hold on;
        plot(R_t,R_amp,'+k');
        title('ECG with R-peaks (black) and S-points (green) over ECG')
        ylabel('ECG+R+S'); xlim([0 tm(end)]); hold off;
    ax(4) = subplot(4,1,4); plot(R_t(1:length(hrv)),hrv,'r+')
        hold on, title('HR')
        ylabel('RR (s)'); xlim([0 tm(end)]);

    %linkaxes(ax,'x');
    set(gca,'FontSize',FONTSIZE);
    allAxesInFigure = findall(gcf,'type','axes');
    set(allAxesInFigure,'fontSize',FONTSIZE);
end

% NOTES
%   Finding the P&T energy threshold: in order to avoid crash due to local
%   huge bumps, threshold is choosen at 98-99% of amplitude distribution.
%   first sec removed for choosing the thres because of filter init lag.
%
%   Search back: look for missed peaks by lowering the threshold in area where the
%   RR interval variability (RRv) is higher than 1.5*medianRRv
%
%   Sign of the QRS (signForce): look for the mean sign of the R-peak over the
%   first 30sec when looking for max of abs value. Then look for the
%   R-peaks over the whole record that have this given sign. This allows to
%   not alternate between positive and negative detections which might
%   happen in some occasion depending on the ECG morphology. It is also
%   better than forcing to look for a max or min systematically.

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [output] = writeTiffData(imgdata,FileSaving,varargin)
	% Writes out TIF data.
	% Created by Jerome Lecoq in 2012
	% Separate function by Biafra Ahanonu
	% 2015.07.06 [19:29:37]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try

		% FileSaving=fullfile(handles.Path,handles.File);
		% FromFrame=str2num(get(handles.FromFrame,'String'));
		% ToFrame=str2num(get(handles.ToFrame,'String'));
		NumberFrame=size(imgdata,3);
		TifFile = Tiff(FileSaving,'w');
		% We populate the tag on pixel format depending on the current data
		% type
		switch class(imgdata)
		    case 'logical'
		        tagstruct.SampleFormat=1;
		        tagstruct.BitsPerSample=1;
		    case 'uint16'
		        tagstruct.SampleFormat=1;
		        tagstruct.BitsPerSample=16;
		    case 'int16'
		        tagstruct.SampleFormat=2;
		        tagstruct.BitsPerSample=16;
		    case 'single'
		        tagstruct.SampleFormat=3;
		        tagstruct.BitsPerSample=32;
		    case 'double'
		        tagstruct.SampleFormat=3;
		        tagstruct.BitsPerSample=64;
		    case 'uint8'
		        tagstruct.SampleFormat=1;
		        tagstruct.BitsPerSample=8;
		    case 'int8'
		        tagstruct.SampleFormat=2;
		        tagstruct.BitsPerSample=8;
		    case 'uint32'
		        tagstruct.SampleFormat=1;
		        tagstruct.BitsPerSample=32;
		    case 'int32'
		        tagstruct.SampleFormat=2;
		        tagstruct.BitsPerSample=32;
		end

		CompressionScheme = 1;
		% We populate compression scheme
		switch CompressionScheme
		    case 1
		        % No compression
		        tagstruct.Compression=1;
		    case 2
		        % LZW lossless
		        tagstruct.Compression=5;
		    case 3
		        % Packbits lossless
		        tagstruct.Compression=32773;
		    case 4
		        % Deflate lossless
		        tagstruct.Compression=32946;
		end

		tagstruct.ImageLength = size(imgdata,1);
		tagstruct.ImageWidth = size(imgdata,2);
		tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
		tagstruct.SamplesPerPixel = 1;
		tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
		tagstruct.Software = 'MATLAB';

		for i=1:NumberFrame
		    TifFile.setTag(tagstruct);
		    TifFile.write(imgdata(:,:,i));
		    TifFile.writeDirectory();
		end
		TifFile.close();

	catch err
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
clear Models, clear MM

% User-Specified Parameters
% ===============================================================

% Conditions and stimuli
% ---------------------------------------------------------------
GA.conditions = [1 2 3 4 5 6];
	% a vector of integers 1:# of event-related conditions (e.g., [1 2 3 4])
	% to add "rest" intervals, add an extra condition and set its contrast weight
	% to zero in all contrasts.
	% e.g., condition 3 in this design could be passive rest.

GA.freqConditions = [1 1 1 1 1 1] ./ 6;
	% vector of frequencies of each trial type; should sum to 1

GA.scanLength = 36;
	% how long your run is, in seconds.
    % This and the ISI determine how many stimuli are in each design vector

GA.ISI = 1;
	% how long between stimulus presentations? (you can also include "rest" presentations)
	% also the time resolution of stimulus condition function (list of stimuli)
	% designs will be constructed in time units of this resolution

GA.TR = 1;
    % the TR (sampling resolution) of your experiment; time for volume acquisition

GA.doepochs = 0;        % build epochs instead of events
	% use this if you have epochs instead of events; alpha version.  use with caution.
	% if you use this, use the GA function optimizeGA_epochs.m
	% also, the contrast matrix will be different, with 3 periods per trial explicitly modeled
	% this wasn't intended to be a general solution, but a specific optimization for a certain design.

% Hox optimization parameters
% A hox gene is a master gene.  These numbers control the stimulus
% parameters, which can be optimized within the GA rather than
% pre-specified.
% Hox elements in this GA go in the stimlist as the first several elements.
% hox sequence currently codes, in number order:
% ISI TR cuelen cuerest stimlen stimrest resplen resprest (in s)
% the last 6 parameters are for use with doepochs.
% They code epoch lengths for 6 periods within each trial.
% TR does not work right now
% use zeros to use default rather than variable hox parameters
% ---------------------------------------------------------------
GA.numhox = 0;
GA.hoxrange = [];
    % rows are each hox gene's allowable range, cols index hox elements
	% use this if you want to determine the ISI, etc. using the GA
	% also may require significant user input/program modification at this stage.

% Genetic algorithm parameters
% ---------------------------------------------------------------
nmodels = 1;
	% how many runs do you want to optimize?
	% GA runs one separate optimization for each model.

GA.cbalColinPowerWeights = [1 0 0 1];	% 1 = cbal, 2 = eff, 3 = hrf shape, 4 = freq
	% first element: counterbalancing of stimuli
	% second: contrast detection efficiency
	% third: hrf shape estimation efficiency
	% maintenance of input frequencies for each trial type
	% fitness scores for each measure are multiplied by these values
	% before collapsing to a single, final fitness measure for each design.
	% does not have to sum to 1.

GA.numGenerations = 30;
	% how many iterations of the GA to run.
GA.sizeGenerations = 500;
	% how many designs to test  for each generation?  Population size.
GA.maxTime = 30;
	% max time to run in s, or Inf for infinite time
	% The GA stops when either numGenerations or maxTime is reached.

GA.alph = 2.1;
	% "selection pressure": higher is more extreme selection; always pick the best 50% of designs
	%1 is no selection pressure at all, or random selection of designs for recombination
	% selection pressure is like in evolution; refers to whether the best designs are selected to
	% continue on to the next generation.
	% an intermediate value is best, to keep the population heterogeneity high.

GA.plotFlag = 0;
	% plot results after each run.
	% Recommended to leave this off if nmodels > 1, and then plot your final results later.

% Filtering, counterbalancing, and design tolerance
% ---------------------------------------------------------------
% Empty brackets indicate the option is not to be used.

GA.lowerLimit = [];
GA.HPlength = [];
	% high-pass filter length, in s, for analysis; [] for no HP filter
	% Used to filter out noise at lower frequencies than design during data analysis
	% The cutoff you should use depends on how much power in your design is below the cutoff -
	% you want all the power in your design to be above the cutoff.

GA.LPsmooth = [];
	% low-pass filter type for analysis; [] for no LP filter
	% choices are 'hrf' or []

GA.maxOrder = 1;
    % order of counterbalancing to use
    % 1 = each trial type follows each other with equal probabilities, adjusted for base frequencies
    % 2 = one-back and two-back counterbalancing.  3 = 1+2+3 back, etc.

% Hard constraints
% ---------------------------------------------------------------
% Here you can specify tolerances for designs
% If designs do not meet these criteria, they will receive -infinity fitness scores
% Empty brackets indicate the option is not to be used.

GA.NumStimthresh = [];      % maximum number of repeats of any one event type in a row
GA.maxCbalDevthresh = []; 	% maximum acceptable counterbalancing deviation (actual freq - expected freq)
GA.maxFreqDevthresh = .2;   % maximum acceptable deviation from input frequencies (GA.freqConditions)

% Contrast setup
% ---------------------------------------------------------------
% Here you specify contrasts across conditions, for use with the contrast estimation fitness measure
% Each contrasts should be a row in the matrix GA.contrasts
%
GA.contrasts = [1 0 0 0 0 0];
    % there should be one column per condition in your design, not including the intercept.
    %if trans2switch or trans2block = 1, double the number of columns.
    % when using epoch design, three elements per condition for epochs (3 epochs)!

GA.contrastweights = [1];
    % Weighting function for contrasts (rows of GA.contrasts)
    % Contrast efficiencies will be multiplied by these weights before computing overall design fitness
    % If no contrasts are specified, this vector can specify weights for predictors (columns)

% Autocorrelation and special options
% ---------------------------------------------------------------
AutocorrelationFileName = 'myscannerxc';
    % This is the name of a mat file without the .mat extension
    % In the mat file, there should be a variable called myscannerxc
    % which is a row vector containing the autocorrelation fucntion you wish to use.
    % This function is used as the intrinsic noise autocorrelation estimate.
    % Included with the GA scripts are these functions:
    % myscannerxc   (U of M 3T, 7 subjects)
    % hiautocorr    (1/f)
    % hiautocorr3   (another 1/f)
    % noautocorr    (identity; no autocorrelation)

GA.restlength = [];
    % if inserting probe or instruction periods, modeled with periodic rests in the simulation,
    % this is the length of the rest/probe/instr periods to insert, in units of the ISI
GA.restevery = [];
    % This is how many regular events you want to have (in ISIs) between rests.
    % Leave restlength and restevery blank to avoid using these options.
GA.trans2switch = 0;
    % 1 or 0.  This option works if the conditions field is 2 elements long.
    % This option creates new conditions 3 and 4 that occur whenever 1 is followed by 2 or 2 by 1
    % So this option creates trial history-dependent predictors.  1 and 2 are repeats, 3 and 4 are switches
    % Useful when studying habituation, switching between items, etc.
GA.trans2block = 0;
    % 1 or 0.  This option doubles the input length of conditions - so [1 2] is transformed to [1 2 3 4]
    % This option creates alternating ABAB blocks of your input event types (e.g., 1 2) with the new
    % event types (e.g., 3 4).  This is useful when simulating mixed block/event related designs.
    % Trial types 1-2, for example, can be two events within block A, and 3-4 can be the same (or different)
    % events within block B.  The alternation frequency (in ISIs) is specified by GA.restevery.
    % If you use this option, the contrasts field must contain twice as many columns as you have input conditions.
    % for example, if conditions = [1 2 3 4], contrasts = [1 1 1 1 -1 -1 -1 -1] specifies block A vs block B.
    % contrasts = [1 1 1 1 -1 -1 -1 -1; 1 1 -1 -1 1 1 -1 -1] specifies A vs B and events [1 2 5 6] - [3 4 7 8].
GA.dofirst = 0;
    % 1 or 0.  This creates a separate trial type for the first trial following rests in GA.restlength.
    % It only works if you're using trans2switch and trans2block together; unreliable under other circumstances.
GA.nonlinthreshold = [];
    % If the value here is x, predictor heights are clipped (thresholded) at x times the unit HRF height.
    % If empty, no thresholding is performed.
    % This is useful for modeling nonlinearities in the BOLD signal, which is important if your design uses
    % short ISIs (generally, 2 s or less).  The clipping acts as a 'saturation' factor.
    % I like to use 2, because data in our lab suggests the the nonlinearities in the HRF can be
    % roughly approximated by a clipping function with this parameter value.

% ---------------------------------------------------------------
%
% * E N D   U S E R   I N P U T
%
% ---------------------------------------------------------------
eval(['load ' AutocorrelationFileName]);
GA.xc = myscannerxc;

% =============================================
% * vary by parameter - set this in this script
% * if running multiple models, nmodels > 1
% ---------------------------------------------

varyparam = 0;
fieldname = 'alph';		% or 'freqConditions(end), etc.
incrementby = 0.3;
incrementevery = 5;		% every n models
contrastsofinterest = 1;% contrasts or predictors, if no cons specified

% =============================================

if varyparam
	eval(['paramvalue = GA.' fieldname ';'])
	disp(' ');disp('* *********************************************************************')
	disp('*  You have selected to vary a parameter across models in ga_gui_script')
	disp(['*  GA.' fieldname ' starting at ' num2str(paramvalue)])
	disp('* *********************************************************************')
end

for nm = 1:nmodels

   eval(['diary model' num2str(nm) '.log'])
   disp(['Starting Model ' num2str(nm)])
    % freqConditions

   M = optimizeGA(GA);

   Models(:,nm) = M.stimlist;
   MM{nm} = M;
   save GAworkspace

   if varyparam
   	if isfield(M,'consebeta'),
		Fit(1,nm) = mean(M.consebeta(contrastsofinterest));
	else
		Fit(1,nm) = mean(M.sebeta(contrastsofinterest));
	end
   	eval(['Fit(2,nm) = GA.' fieldname ';'])
   	FC(nm,:) = GA.freqConditions;
   	ALPH(nm,:) = GA.alph;
	str = ['if mod(nm,' num2str(incrementevery) ') == 0,GA.' fieldname ' = GA.' fieldname '+' num2str(incrementby) ';,end']
   	eval(str);
	eval(['paramvalue = GA.' fieldname ';'])
	disp('* *********************************************************************')
	disp('*  Varying parameter option on.')
	disp(['*  GA.' fieldname ' is now ' num2str(paramvalue)])
	disp('* *********************************************************************')
   end

   eval(['save model' num2str(nm) '_' M.date ' M'])
   save GAworkspace
   diary off
   fprintf('\n')
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [oz] = qt2oz(qt)
% Convert volume from US liquid quarts to US liquid ounces.
% Chad Greene 2012
oz = qt*32;


text:
function p = linear_search(A,t)
%% linear Search
% This function linear searches target value (t) in array A.
% For this, It sequentially checks each entry of the array until a match is found
% or  or the whole list has been searched.
% If it can find the target returns 1 otherwise 0.

array_length = length(A);
i = 1;
searchTermination = 0;

while searchTermination == 0 && i < array_length+1
    if A(i) == t
        p = 1;
        searchTermination = 1;
        disp('the target is found in the array')
    else
        i = i+1;
    end
end
if i == array_length+1
    p = 0;
    disp('the target is not found in the array')
end
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [overlaps, score, class_ids, avg_cube, avg_cube_per_class] = detector_wiggle_rcnn(filename)
  ld = load(filename);
  scoremaps = ld.scoremaps;

  n = numel(scoremaps);
  all_overlap_scores = cell(n,1);
  for i = 1:n, s = scoremaps(i);
    gt = double(s.gt_box);
    boxes = double(s.boxes);
    invalid = isinf(boxes(:,5)) | isnan(boxes(:,5));
    boxes = boxes(~invalid,:);

    t_overlap = overlap(gt, boxes);
    t_class_ids = repmat(s.class, size(t_overlap));
    all_overlap_scores{i} = cat(2, t_overlap, boxes(:,5), double(t_class_ids));
  end

  all_overlap_scores = cat(1, all_overlap_scores{:});
  overlaps = all_overlap_scores(:,1);
  score = all_overlap_scores(:,2);
  class_ids = all_overlap_scores(:,3);

  avg_cube = [];
  avg_cube_per_class = [];
  if isfield(scoremaps, 'score_map') && ~isempty(scoremaps(1).score_map)
    all_maps = cat(4,scoremaps.score_map);
    avg_cube = nanmean(all_maps, 4);

    classes = unique([scoremaps.class]);
    avg_cube_per_class = cell(1, numel(classes));
    for c = 1:numel(classes)
      t_scoremaps = scoremaps([scoremaps.class] == classes(c));
      avg_cube_per_class{c} = nanmean(cat(4,t_scoremaps.score_map), 4);
    end
  end
end

author: hosang

repo: detection-proposals


text:
function [xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk,dxpdalpha] = project_points_fisheye(X,om,T,f,c,k,alpha)

%project_points2.m
%
%[xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk] = project_points_fisheye(X,om,T,f,c,k,alpha)
%
%Projects a 3D structure onto the image plane of a fisheye camera.
%
%INPUT: X: 3D structure in the world coordinate frame (3xN matrix for N points)
%       (om,T): Rigid motion parameters between world coordinate frame and camera reference frame
%               om: rotation vector (3x1 vector); T: translation vector (3x1 vector)
%       f: camera focal length in units of horizontal and vertical pixel units (2x1 vector)
%       c: principal point location in pixel units (2x1 vector)
%       k: Distortion fisheye coefficients (5x1 vector)
%       alpha: Skew coefficient between x and y pixel (alpha = 0 <=> square pixels)
%
%OUTPUT: xp: Projected pixel coordinates (2xN matrix for N points)
%        dxpdom: Derivative of xp with respect to om ((2N)x3 matrix)
%        dxpdT: Derivative of xp with respect to T ((2N)x3 matrix)
%        dxpdf: Derivative of xp with respect to f ((2N)x2 matrix if f is 2x1, or (2N)x1 matrix is f is a scalar)
%        dxpdc: Derivative of xp with respect to c ((2N)x2 matrix)
%        dxpdk: Derivative of xp with respect to k ((2N)x5 matrix)
%
%Definitions:
%Let P be a point in 3D of coordinates X in the world reference frame (stored in the matrix X)
%The coordinate vector of P in the camera reference frame is: Xc = R*X + T
%where R is the rotation matrix corresponding to the rotation vector om: R = rodrigues(om);
%call x, y and z the 3 coordinates of Xc: x = Xc(1); y = Xc(2); z = Xc(3);
%The pinehole projection coordinates of P is [a;b] where a=x/z and b=y/z.
%call r^2 = a^2 + b^2,
%call theta = atan(r),
%Fisheye distortion -> theta_d = theta * (1 + k(1)*theta^2 + k(2)*theta^4 + k(3)*theta^6 + k(4)*theta^8)
%
%The distorted point coordinates are: xd = [xx;yy] where:
%
%xx = (theta_d / r) * x
%yy = (theta_d / r) * y
%
%Finally, convertion into pixel coordinates: The final pixel coordinates vector xp=[xxp;yyp] where:
%
%xxp = f(1)*(xx + alpha*yy) + c(1)
%yyp = f(2)*yy + c(2)
%
%
%NOTE: About 90 percent of the code takes care fo computing the Jacobian matrices
%
%
%Important function called within that program:
%
%rodrigues.m: Computes the rotation matrix corresponding to a rotation vector
%
%rigid_motion.m: Computes the rigid motion transformation of a given structure

if nargin < 7,
    alpha = 0;
    if nargin < 6,
        k = zeros(4,1);
        if nargin < 5,
            c = zeros(2,1);
            if nargin < 4,
                f = ones(2,1);
                if nargin < 3,
                    T = zeros(3,1);
                    if nargin < 2,
                        om = zeros(3,1);
                        if nargin < 1,
                            error('Need at least a 3D structure to project (in project_points.m)');
                            return;
                        end;
                    end;
                end;
            end;
        end;
    end;
end;

[m,n] = size(X);

if nargout > 1,
    [Y,dYdom,dYdT] = rigid_motion(X,om,T);
else
    Y = rigid_motion(X,om,T);
end;

inv_Z = 1./Y(3,:);

x = (Y(1:2,:) .* (ones(2,1) * inv_Z)) ;

bb = (-x(1,:) .* inv_Z)'*ones(1,3);
cc = (-x(2,:) .* inv_Z)'*ones(1,3);

if nargout > 1,
    dxdom = zeros(2*n,3);
    dxdom(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(1:3:end,:) + bb .* dYdom(3:3:end,:);
    dxdom(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(2:3:end,:) + cc .* dYdom(3:3:end,:);

    dxdT = zeros(2*n,3);
    dxdT(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(1:3:end,:) + bb .* dYdT(3:3:end,:);
    dxdT(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(2:3:end,:) + cc .* dYdT(3:3:end,:);
end;

% Add fisheye distortion:

r2 = x(1,:).^2 + x(2,:).^2;

if nargout > 1,
    dr2dom = 2*((x(1,:)')*ones(1,3)) .* dxdom(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdom(2:2:end,:);
    dr2dT = 2*((x(1,:)')*ones(1,3)) .* dxdT(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdT(2:2:end,:);
end;

% Radial distance:
r = sqrt(r2);
if nargout > 1,
    drdr2 = ones(1,length(r));
    drdr2(r>1e-8) = 1 ./ (2*r(r>1e-8));

    drdom = [ (drdr2').*dr2dom(:,1)   (drdr2').*dr2dom(:,2)  (drdr2').*dr2dom(:,3)  ];
    drdT = [ (drdr2').*dr2dT(:,1)   (drdr2').*dr2dT(:,2)  (drdr2').*dr2dT(:,3)  ];
end;

% Angle of the incoming ray:
theta = atan(r);
if nargout > 1,
    dthetadr = 1 ./ (1 + r2);

    dthetadom = [ (dthetadr').*drdom(:,1)   (dthetadr').*drdom(:,2)  (dthetadr').*drdom(:,3) ];
    dthetadT = [ (dthetadr').*drdT(:,1)   (dthetadr').*drdT(:,2)  (dthetadr').*drdT(:,3) ];
end;

% Add the fisheye distortion:

theta2 = theta.^2;
theta3 = theta2.*theta;
theta4 = theta2.^2;
theta5 = theta4.*theta;
theta6 = theta3.^2;
theta7 = theta6.*theta;
theta8 = theta4.*theta4;
theta9 = theta8.*theta;

% Fisheye distortion -> theta_d = theta * (1 + k(1)*theta2 + k(2)*theta4 + k(3)*theta6 + k(4)*theta8)

theta_d = theta + k(1)*theta3 + k(2)*theta5 + k(3)*theta7 + k(4)*theta9;

if nargout > 1,
    dtheta_ddtheta = 1 + 3*k(1)*theta2 + 5*k(2)*theta4 + 7*k(3)*theta6 + 9*k(4)*theta8;
    dtheta_ddom = [ (dtheta_ddtheta').*dthetadom(:,1)   (dtheta_ddtheta').*dthetadom(:,2)  (dtheta_ddtheta').*dthetadom(:,3) ];
    dtheta_ddT = [ (dtheta_ddtheta').*dthetadT(:,1)   (dtheta_ddtheta').*dthetadT(:,2)  (dtheta_ddtheta').*dthetadT(:,3) ];
    dtheta_ddk = [theta3' theta5' theta7' theta9'];
end;

% ratio:
inv_r = ones(1,length(r));
inv_r(r>1e-8) = 1./r(r>1e-8);

cdist = ones(1,length(r));
cdist(r > 1e-8) = theta_d(r > 1e-8) ./ r(r > 1e-8);

if nargout > 1,
    dcdistdom = [ ((inv_r').*(dtheta_ddom(:,1) - (cdist').*drdom(:,1)))  ((inv_r').*(dtheta_ddom(:,2) - (cdist').*drdom(:,2)))   ((inv_r').*(dtheta_ddom(:,3) - (cdist').*drdom(:,3))) ];
    dcdistdT = [ ((inv_r').*(dtheta_ddT(:,1) - (cdist').*drdT(:,1)))  ((inv_r').*(dtheta_ddT(:,2) - (cdist').*drdT(:,2)))   ((inv_r').*(dtheta_ddT(:,3) - (cdist').*drdT(:,3))) ];
    dcdistdk = [ (inv_r'.*dtheta_ddk(:,1)) (inv_r'.*dtheta_ddk(:,2)) (inv_r'.*dtheta_ddk(:,3))  (inv_r'.*dtheta_ddk(:,4)) ];
end;

xd1 = x .* (ones(2,1)*cdist);

if nargout > 1,
    dxd1dom = zeros(2*n,3);
    dxd1dom(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdom;
    dxd1dom(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdom;
    coeff = (reshape([cdist;cdist],2*n,1)*ones(1,3));
    dxd1dom = dxd1dom + coeff.* dxdom;

    dxd1dT = zeros(2*n,3);
    dxd1dT(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdT;
    dxd1dT(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdT;
    dxd1dT = dxd1dT + coeff.* dxdT;

    dxd1dk = zeros(2*n,4);
    dxd1dk(1:2:end,:) = (x(1,:)'*ones(1,4)) .* dcdistdk;
    dxd1dk(2:2:end,:) = (x(2,:)'*ones(1,4)) .* dcdistdk;
end;

% No tangential distortion:
xd2 = xd1;
if nargout > 1,
    dxd2dom = dxd1dom;
    dxd2dT = dxd1dT;
    dxd2dk = dxd1dk;
end;

% Add Skew:
xd3 = [xd2(1,:) + alpha*xd2(2,:);xd2(2,:)];

% Compute: dxd3dom, dxd3dT, dxd3dk, dxd3dalpha
if nargout > 1,
    dxd3dom = zeros(2*n,3);
    dxd3dom(1:2:2*n,:) = dxd2dom(1:2:2*n,:) + alpha*dxd2dom(2:2:2*n,:);
    dxd3dom(2:2:2*n,:) = dxd2dom(2:2:2*n,:);
    dxd3dT = zeros(2*n,3);
    dxd3dT(1:2:2*n,:) = dxd2dT(1:2:2*n,:) + alpha*dxd2dT(2:2:2*n,:);
    dxd3dT(2:2:2*n,:) = dxd2dT(2:2:2*n,:);
    dxd3dk = zeros(2*n,4);
    dxd3dk(1:2:2*n,:) = dxd2dk(1:2:2*n,:) + alpha*dxd2dk(2:2:2*n,:);
    dxd3dk(2:2:2*n,:) = dxd2dk(2:2:2*n,:);
    dxd3dalpha = zeros(2*n,1);
    dxd3dalpha(1:2:2*n,:) = xd2(2,:)';
end;

% Pixel coordinates:
if length(f)>1,
    xp = xd3 .* (f * ones(1,n))  +  c*ones(1,n);
    if nargout > 1,
        coeff = reshape(f*ones(1,n),2*n,1);
        dxpdom = (coeff*ones(1,3)) .* dxd3dom;
        dxpdT = (coeff*ones(1,3)) .* dxd3dT;
        dxpdk = (coeff*ones(1,4)) .* dxd3dk;
        dxpdalpha = (coeff) .* dxd3dalpha;
        dxpdf = zeros(2*n,2);
        dxpdf(1:2:end,1) = xd3(1,:)';
        dxpdf(2:2:end,2) = xd3(2,:)';
    end;
else
    xp = f * xd3 + c*ones(1,n);
    if nargout > 1,
        dxpdom = f  * dxd3dom;
        dxpdT = f * dxd3dT;
        dxpdk = f  * dxd3dk;
        dxpdalpha = f .* dxd3dalpha;
        dxpdf = xd3(:);
    end;
end;

if nargout > 1,
    dxpdc = zeros(2*n,2);
    dxpdc(1:2:end,1) = ones(n,1);
    dxpdc(2:2:end,2) = ones(n,1);
end;

return;

% Test of the Jacobians:

n = 10;

X = 10*randn(3,n);
om = randn(3,1);
T = [10*randn(2,1);40];
f = 1000*rand(2,1);
c = 1000*randn(2,1);
k = 0.5*randn(4,1);
alpha = 0.01*randn(1,1);

[x,dxdom,dxdT,dxdf,dxdc,dxdk,dxdalpha] = project_points_fisheye(X,om,T,f,c,k,alpha);

% Test on om: not OK

dom = 0.00000000001 * norm(om)*randn(3,1);
om2 = om + dom;

[x2] = project_points_fisheye(X,om2,T,f,c,k,alpha);

x_pred = x + reshape(dxdom * dom,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on T: not OK

dT = 0.0001 * norm(T)*randn(3,1);
T2 = T + dT;

[x2] = project_points_fisheye(X,om,T2,f,c,k,alpha);

x_pred = x + reshape(dxdT * dT,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on f: OK!!

df = 0.001 * norm(f)*randn(2,1);
f2 = f + df;

[x2] = project_points_fisheye(X,om,T,f2,c,k,alpha);

x_pred = x + reshape(dxdf * df,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on c: OK!!

dc = 0.01 * norm(c)*randn(2,1);
c2 = c + dc;

[x2] = project_points_fisheye(X,om,T,f,c2,k,alpha);

x_pred = x + reshape(dxdc * dc,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on k: OK!!

dk = 0.00001 * norm(k)*randn(4,1);
k2 = k + dk;

[x2] = project_points_fisheye(X,om,T,f,c,k2,alpha);

x_pred = x + reshape(dxdk * dk,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on alpha: OK!!

dalpha = 0.001 * norm(k)*randn(1,1);
alpha2 = alpha + dalpha;

[x2] = project_points_fisheye(X,om,T,f,c,k,alpha2);

x_pred = x + reshape(dxdalpha * dalpha,2,n);

norm(x2-x)/norm(x2 - x_pred)

repo: ekfmonoslam


text:
%                     benchmarkTransmissionMatrix.m
%
% Code adapted from Sudarshan Nagesh for the reconstruction of images from
% a transmission matrix.
%
% This script will reconstruct images from phaseless measurements that were
% acquired when imaging through a diffusive media.  The quality of the
% reconstruction will be computed by comparing to a known solution, and
% different algorithms are compared.
%   Three datasets are available with resolutions 16x16, 40x40, and 64x64.
% At each resolution, there are 5 sets of measurements, each corresponding
% to a different image.  The user must select which resolution and dataset
% they want to run on.
%   As the script runs, all reconstructed images will be saved inside the
% 'benchmarkResults' folder.
%
% Note: The empirical dataset must be downloaded and installed into the
% 'data' directory for this to work. See the user guide.
%
%
%%I/O
%  Inputs:
%  imageSize            : Size of image to reconstruct.  Must be {16,40,64}
%  datasetSelection     : Choose which of the sets of measurements to use.
%                           Must be in {1,2,3,4,5}.
%  residualConstant     : Only use rows of the transmission matrix that
%                           had residuals less than this cutoff. Must be
%                           between 0 and 1.  Recommend 0.4.
%  algorithms           : A cell-array of structs that describe algorithms.
%
%
%%References
%  Paper Title:   Coherent inverse scattering via transmission matrices:
%  Efficient phase retrieval algorithms and a public dataset.
%  Authors:       Christopher A. Metzler,  Manoj
%  K. Sharma,  Sudarshan Nagesh,  Richard G. Baraniuk,
%                 Oliver Cossairt,  Ashok Veeraraghavan
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

function benchmarkTransmissionMatrix(imageSize, datasetSelection, residualConstant, algorithms)

    %% Load the transmission matrix, and phaseless measurements
    [A,Xt,b] = loadData(residualConstant, datasetSelection, imageSize);

    % This string contains a text-based table with all the results, and is
    % printed after all calculations are done.
    % Results will also be printed incrementally as the benchmark runs.
    resultsTable = '\nAlgorithm  |  Measurement Error \n';

    %% Loop over each algorithm, and perform reconstruction and evaluation
    for k=1:length(algorithms)
    	opts=algorithms{k};
        fprintf('Running Algorithm: %s\n',opts.algorithm);

        % Call a phase retrieval routine, and record the error of the result
        [reconError, measurementError] = reconstructSignal(A,Xt,b,imageSize,datasetSelection,opts);

        % Record the error in the table
        resultsTable = [resultsTable,...
                        pad(opts.algorithm,16),...
                        sprintf('%0.5f',measurementError),...
                        '\n'];
    end

    %% print the table
    fprintf(resultsTable);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                        Utility functions                        %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load the measurement operator and measurements that correspond to the
% image size and dataset chosen by the user.  This method returns the
% measurement operator 'A' (a dense matrix), a ground truth solution
% 'Xt', and a set of measurements 'Y'.
function [A,Xt,b] = loadData(residualConstant, datasetSelection, imageSize)
    % Create some strings for dealing with different filenames according to the dataset chosen
    switch imageSize
        case {16}
            measurementType = 'AmpSLM_16x16';
        case {40}
            measurementType = 'PhaseSLM_40x40';
        case {64}
            measurementType = 'AmpSLM_64x64';
        otherwise
            error('illegal imageSize: %d. It should be chosen from {16, 40, 64}',imageSize);
    end

    % Make sure the dataset selection is valid
    if imageSize==16 || imageSize==40
        assert(datasetSelection <= 5);
    elseif imageSize==64
        assert(datasetSelection <= 6);
    end

    %% Load the transmission matrix, ground truth image, and measurements
    dataRoot = getFolderPath('data');       % The location of the 'data' folder
    try
        disp(['Loading transmission matrix (A_prVAMP.mat), this may take several minutes']);
        load(strcat(dataRoot,'TransmissionMatrices/Coherent_Data/',measurementType,'/A_prVAMP.mat'));
        disp(['Loading test image (XH_test.mat)']);
        load(strcat(dataRoot,'/TransmissionMatrices/Coherent_Data/',measurementType,'/XH_test.mat'));
        disp(['Loading measurements (YH_squared_test.mat)']);
        load(strcat(dataRoot,'TransmissionMatrices/Coherent_Data/',measurementType,'/YH_squared_test.mat'));
    catch
        error(['You are missing the transmission matrix dataset.  To download it, Go to',...
            'https://rice.app.box.com/v/TransmissionMatrices/,'...
            'and download the "TransmissionMatrices" folder. Unzip the downloaded folder and',...
            'place it in the in "data" folder. See the user guide for details.']);
    end
    % Only use the most accurate rows of the transmission matrix.  This is
    % determined by checking that the residuals found during row calcuation
    % are less than the specificed "residualConstant".
    findPixels = find(residual_vector<residualConstant);
    YH_squared_test = YH_squared_test(:,findPixels);

    A = A(findPixels,:);

    % Unpack the measurement data
    Y = double(YH_squared_test(datasetSelection,:));
    b = sqrt(Y);

    % Unpack the ground-truth image Xt
    % Note:  this image is real-valued
    Xt = double(XH_test(datasetSelection,:))';
    Xt = Xt/max(abs(Xt));
    % reshape original image
    originalImage = imresize(reshape(real(Xt),imageSize,imageSize),4);
    originalImageName = ['TM',num2str(imageSize),'-',num2str(datasetSelection),'-original'];
    imshow(originalImage,[]);  title(originalImageName,'fontsize',16); drawnow;
    % save the image to a folder
    saveBenchmarkImage(originalImage,['TM',num2str(imageSize)],originalImageName);
end

% Reconstruct a signal from a given measurement operator and set of
% measurements.  Compare the results to the ground-truth Xt, and report
% error.
function [reconError, measurementError] = reconstructSignal(A,Xt,Y,imageSize,datasetSelection,opts)
    fprintf('    Reconstructing image...')

    % Convenient variables
    n = imageSize^2;                    % Size of SLM
    imageSizeStr = num2str(imageSize);  % For TM16, this is '16'
    ind = datasetSelection;             % The index of the selected dataset

    % Solve the PR problem
    [X,outs,opts] = solvePhaseRetrieval(A,[],Y(:),n,opts);

    % Compute solution quality by solving a system to get the best affine
    % transform that maps the recovered solution onto the true solution
    Xrec = [X(:), ones(numel(X),1)];
    coeffs = Xrec\Xt(:);
    Xrec = Xrec*coeffs;

    % Reshape images into a square, and save them.
    realIm = imresize(reshape(real(Xrec),imageSize,imageSize),4);
    recoveredImageNameAbs = strcat('TM',imageSizeStr,'-',num2str(ind),'-',opts.algorithm,...
                                num2str(outs.iterationCount),'-real');
    imshow(realIm,[]); title(recoveredImageNameAbs,'fontsize',16); drawnow;
    realIm = realIm-min(realIm(:));
    realIm = realIm/max(realIm(:));
    saveBenchmarkImage(realIm,['TM',imageSizeStr],recoveredImageNameAbs);

    reconError = norm(Xt(:)-Xrec(:))/norm(Xt(:));
    measurementError = norm(abs(A*X)-Y(:))/norm(Y(:));
    %fprintf('    Relative reconstruction error = %s\n',reconError);
    fprintf('    Relative measurement error = %s\n',measurementError);
end

% Find the path for the specified folder by searching the current and
% parent directories
function path = getFolderPath(fname)
    d = dir;
    % check if the folder fname is in current path
    if any(strcmp(fname, {d([d.isdir]).name}))
        path = [fname,'/'];
    else
        % Look for folder fname in parent directory
        d = dir('../');
        if any(strcmp(fname, {d([d.isdir]).name}))
            path = ['../',fname,'/'];
        else
            error(['Unable to find path of folder: ',fname, ...
                '.  Make sure your current directory is the PhasePack root.']);
        end
    end
end

% Save an image to a sub-folder inside 'benchmarkResults'
function saveBenchmarkImage(image,folder,fname)
    bmRoot = getFolderPath('benchmarkResults');
    if ~exist([bmRoot,folder],'dir')
         mkdir([bmRoot,folder]);
    end
    fullName = strcat(bmRoot,folder,'/',fname,'.png');
    fprintf('    Saving reconstructed image: %s\n',fullName);
    try
        imwrite(image,fullName);
    catch
        error(['Unable to save image to file: ', fullName]);
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
% This file is part of TREEQSM.
%
% TREEQSM 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.
%
% TREEQSM 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 TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function plot_models_segmentations(P,cover,segment,cylinder,trunk,triangulation)

% ---------------------------------------------------------------------
% PLOT_MODELS_SEGMENTATION.M   Plots the segmented point clouds and
%                               cylinder/triangulation models
%
% Version 1.1.0
% Latest update     13 July 2020
%
% Copyright (C) 2013-2020 Pasi Raumonen
% ---------------------------------------------------------------------

% Inputs:
% P             Point cloud
% cover         cover-structure array
% segment       segment-structure array
% cylinder      cylinder-structure array
% trunk         point cloud of the trunk
% triangulation triangulation-structure array

% Changes from version 1.0.0 to 1.1.0, 13 July 2020:
% 1) plots now figure 1 and 2 with two subplots; in the first the colors
%    are based on branching order and in the second they are based on
%    branch

%% figure 1: branch-segmented point cloud
% colors denote the branching order and branches
figure(1)
subplot(1,2,1)
plot_branch_segmentation(P,cover,segment,'order')
subplot(1,2,2)
plot_branch_segmentation(P,cover,segment,'branch')

%% figure 2: cylinder model
% colors denote the branching order and branches
Sta = cylinder.start;
P = P-Sta(1,:);
if nargin > 5
    trunk = trunk-Sta(1,:);
    Vert = double(triangulation.vert);
    Vert = Vert-Sta(1,:);
end
Sta = Sta-Sta(1,:);
cylinder.start = Sta;
figure(2)
subplot(1,2,1)
plot_cylinder_model(cylinder,'order',2,10)
subplot(1,2,2)
plot_cylinder_model(cylinder,'branch',2,10)

%% figure 3, segmented point cloud and cylinder model
plot_branch_segmentation(P,cover,segment,'order',3,1)
hold on
plot_cylinder_model(cylinder,'order',3,10,0.7)
hold off

if nargin > 4
    %% figure 4, triangulation model (bottom) and cylinder model (top)
    % of the stem
    Facets = double(triangulation.facet);
    CylInd = triangulation.cylind;
    fvd = triangulation.fvd;
    if max(size(Vert)) > 5
        Bran = cylinder.branch;
        nc = size(Bran,1);
        ind = (1:1:nc)';
        C = ind(Bran == 1);
        n = size(trunk,1);
        I = logical(round(0.55*rand(n,1)));
        figure(4)
        point_cloud_plotting(trunk(I,:),4,3)
        patch('Vertices',Vert,'Faces',Facets,'FaceVertexCData',fvd,...
            'FaceColor','flat')
        alpha(1)
        hold on
        plot_cylinder_model(cylinder,'order',4,20,1,(CylInd:C(end)))
        axis equal
        hold off
    else
        disp('No triangulation model generated!')
    end
end

author: InverseTampere

repo: TreeQSM


text:
function calcselective_search( config)
ssconfig = config.selective_search;

images = dir(config.imageLocation);
images = regexpi({images.name}, '.*jpg|.*jpeg|.*png|.*bmp', 'match');
images = [images{:}];
if(~exist([config.outputLocation '/selective_search'], 'dir'))
	mkdir(config.outputLocation,'/selective_search')
end

for i=1:length(images)
        imname = char(images(i));

        fprintf('Calculating Selective Search Object Proposals for %s\n', imname);

        impath = fullfile(config.imageLocation, imname);
        im=imread(impath);

        proposals=calcselective_searchForIm(im,ssconfig);
        saveFile=[imname '.mat'];
        save([config.outputLocation '/selective_search/' saveFile], 'proposals');
end

author: Cloud-CV

repo: object-proposals


text:
classdef testHorizontalLaminate < TestSequentialLaminateTestedWithNumerics

    methods (Access = public)

        function obj = testHorizontalLaminate()
            obj.compute();
        end

    end

    methods (Access = protected)

        function loadLaminateDirection(obj)
            d = [0 1 0];
            dir = Vector3D;
            dir.setValue(d);
            dir.normalize()
            obj.LaminateDirection = dir;
        end

        function loadFiberDirection(obj)
            d = [1 0 0];
            dir = Vector3D;
            dir.setValue(d);
            dir.normalize()
            obj.FiberDirection = dir;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
msckf5_inf = load('msckf_500_1000_min5_maxInf.mat');
msckf10_20 = load('msckf_500_1000_min10_max20.mat');
msckf10_50 = load('msckf_500_1000_min10_max50.mat');
msckf20_100 = load('msckf_500_1000_min20_max100.mat');
swf10 = load('swf_500_1000_10_dataset3.mat');
swf50 = load('swf_500_1000_50_dataset3.mat');
swf100 = load('swf_500_1000_100_dataset3.mat');
imu = load('imu_500_1000.mat');

%% Calculate Average RMSE (Root-Mean-Squared Error)

% IMU Only RMSE
imu_trans = sqrt(sum(imu.msckf_trans_err.^2, 1)/3);
imu_rot = sqrt(sum(imu.msckf_rot_err.^2, 1)/3);

% MSCKF RMSE
msckf5_inf_trans = sqrt(sum(msckf5_inf.msckf_trans_err.^2, 1)/3);
msckf10_20_trans = sqrt(sum(msckf10_20.msckf_trans_err.^2, 1)/3);
msckf10_50_trans = sqrt(sum(msckf10_50.msckf_trans_err.^2, 1)/3);
msckf20_100_trans = sqrt(sum(msckf20_100.msckf_trans_err.^2, 1)/3);
disp('===============')
msckf5_inf_rot = sqrt(sum(msckf5_inf.msckf_rot_err.^2, 1)/3);
msckf10_20_rot = sqrt(sum(msckf10_20.msckf_rot_err.^2, 1)/3);
msckf10_50_rot = sqrt(sum(msckf10_50.msckf_rot_err.^2, 1)/3);
msckf20_100_rot = sqrt(sum(msckf20_100.msckf_rot_err.^2, 1)/3);

% SWF RMSE
swf10_trans = sqrt(sum(swf10.swf_trans_err.^2, 1)/3);
swf50_trans = sqrt(sum(swf50.swf_trans_err.^2, 1)/3);
swf100_trans = sqrt(sum(swf100.swf_trans_err.^2, 1)/3);

swf10_rot = sqrt(sum(swf10.swf_rot_err.^2, 1)/3);
swf50_rot = sqrt(sum(swf50.swf_rot_err.^2, 1)/3);
swf100_rot = sqrt(sum(swf100.swf_rot_err.^2, 1)/3);

repo: MSCKF