text:
classdef imtool3DROI_rect < imtool3DROI

    properties
        fixedAspectRatio = false;
    end

    properties (SetAccess = protected, GetAccess = protected)
        position        %defines the center of the box and its width/height [cx cy width height]
        tbuff           %amount of space (in pixels) to place the text above the ROI
    end

    properties (Dependent = true, Access = public)
        aspectRatio
    end

    methods

        %constructor
        function ROI = imtool3DROI_rect(varargin)

            switch nargin
                case 0  %use the current figure

                    %let the user draw the ROI
                    h = imrect;

                    %get the parent axes
                    ha = get(h,'Parent');

                    %get the handle of the image
                    hi = imhandles(ha);
                    if length(hi)>1
                        for i=1:length(hi)
                            if ndims(get(hi(i),'CData'))<3
                                imageHandle = hi(i);
                            end
                        end
                    else
                        imageHandle = hi;
                    end

                    %get the position
                    pos = getPosition(h);
                    position = [pos(1)+pos(3)/2 pos(2)+pos(4)/2 pos(3) pos(4)];

                    %delete the imroi object
                    delete(h);
                case 1 %user inputs only the handle to the image
                    imageHandle = varargin{1};
                    parent = get(imageHandle,'Parent');
                    h = imrect(parent);
                    pos = getPosition(h);
                    position = [pos(1)+pos(3)/2 pos(2)+pos(4)/2 pos(3) pos(4)];
                    delete(h);
                case 2 %user inputs both the parent handle and a position
                    imageHandle = varargin{1};
                    position = varargin{2};
                case 3
                    imageHandle = varargin{1};
                    position = varargin{2};
                    if isempty(position)
                        parent = get(imageHandle,'Parent');
                        h = imrect(parent);
                        pos = getPosition(h);
                        position = [pos(1)+pos(3)/2 pos(2)+pos(4)/2 pos(3) pos(4)];
                        delete(h);
                    end
                    tool = varargin{3};
            end

            %get the parent axis handle
            parent = get(imageHandle,'Parent');

            %find the top left corner of the box
            pos = [position(1)-position(3)/2 position(2)-position(4)/2 position(3) position(4)];

            %Draw the rectangle at the desired spot
            graphicsHandles(1) = rectangle('Position',pos,'Parent',parent,'EdgeColor','r','LineWidth',1.5,'Curvature',[0 0]);

            %Draw a cross at the center of the box and sqaures on the sides
            nextPlot = get(parent,'NextPlot');  %make sure the new graphics don't delete the old ones
            set(parent,'NextPlot','add');
            graphicsHandles(2) = plot(position(1),position(2),'+r','MarkerSize',12,'Parent',parent); %middle cross
            graphicsHandles(3) = plot(pos(1),pos(2),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %top left corner
            graphicsHandles(4) = plot(pos(1)+pos(3),pos(2),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %top right corner
            graphicsHandles(5) = plot(pos(1),pos(2)+pos(4),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %bottom left corner
            graphicsHandles(6) = plot(pos(1)+pos(3),pos(2)+pos(4),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %bottom right corner
            graphicsHandles(7) = plot(pos(1),position(2),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %left
            graphicsHandles(8) = plot(pos(1)+pos(3),position(2),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %right
            graphicsHandles(9) = plot(position(1),pos(2),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %top
            graphicsHandles(10) = plot(position(1),pos(2)+pos(4),'sr','MarkerSize',8,'Parent',parent,'MarkerFaceColor','r'); %bottom
            set(graphicsHandles,'Clipping','off')
            set(parent,'NextPlot',nextPlot);

            %Define the context menu options (i.e., what happens when you
            %right click on the ROI)
            menuLabels = {'Export stats','Fix Aspect Ratio','Hide Text','Delete','poly2mask'};
            if ~exist('tool','var') || isempty(tool), menuLabels(end) = []; tool = []; end
            menuFunction = @contextMenuCallback;

            %create the ROI object from the superclass
            ROI@imtool3DROI(imageHandle,graphicsHandles,menuLabels,menuFunction, tool);

            %Create the text box
            I = get(ROI.imageHandle,'CData');
            ROI.tbuff = .02*size(I,1);
            ROI.textHandle = text(pos(1),pos(2)-ROI.tbuff,'text','Parent',parent,'Color','w','FontSize',10,'EdgeColor','w','BackgroundColor','k','HorizontalAlignment','Left','VerticalAlignment','bottom','Clipping','on');

            %Set the position property of the ROI
            ROI.position = position;

            %Set the button down functions of the graphics
            for i=1:length(graphicsHandles)
                fun = @(hObject,evnt) ButtonDownFunction(hObject,evnt,ROI,i); set(graphicsHandles(i),'ButtonDownFcn',fun);
            end

            %add a listener for changes in the image. This automatically
            %updates the ROI text when the image changes
            ROI.listenerHandle = addlistener(ROI.imageHandle,'CData','PostSet',@ROI.handlePropEvents);

            %update the text
            newPosition(ROI,position)

        end

        function position = getPosition(ROI)
            position = ROI.position;
        end

        function newPosition(ROI,position, notifoff)

            %Adjust the position if the aspect ratio should be fixed
            if ROI.fixedAspectRatio

            end

            %set the position property of the ROI
            ROI.position = position;

            %find the top left corner of the box
            pos = [position(1)-position(3)/2 position(2)-position(4)/2 position(3) position(4)];

            %get the graphics handles
            graphicsHandles = ROI.graphicsHandles;

            %set the new position of the rectangle and other graphics
            %objects
            set(graphicsHandles(1),'Position',pos);
            set(graphicsHandles(2),'Xdata',position(1),'Ydata',position(2));
            set(graphicsHandles(3),'Xdata',pos(1),'Ydata',pos(2));
            set(graphicsHandles(4),'Xdata',pos(1)+pos(3),'Ydata',pos(2));
            set(graphicsHandles(5),'Xdata',pos(1),'Ydata',pos(2)+pos(4));
            set(graphicsHandles(6),'Xdata',pos(1)+pos(3),'Ydata',pos(2)+pos(4));
            set(graphicsHandles(7),'Xdata',pos(1),'Ydata',position(2));
            set(graphicsHandles(8),'Xdata',pos(1)+pos(3),'Ydata',position(2));
            set(graphicsHandles(9),'Xdata',position(1),'Ydata',pos(2));
            set(graphicsHandles(10),'Xdata',position(1),'Ydata',pos(2)+pos(4));

            %get the ROI measurements
            stats = getMeasurements(ROI);

            %set the textbox
            str = {['Mean: ' num2str(stats.mean,'%+.2f')], ['STD:     ' num2str(stats.STD,'%.2f')]};
            set(ROI.textHandle,'String',str,'Position',[pos(1) pos(2)-ROI.tbuff]);

            %notify a new position
            if ~exist('notifoff','var') || ~notifoff
                notify(ROI,'newROIPosition');
            end

        end

        function newPositionSameSize(ROI,pos)
            position = ROI.position;
            position(1:2)=pos;
            newPosition(ROI,position);
        end

        function [x, y] = getPoly(ROI)
            %get the position
            position = ROI.position;

            %find the top left corner of the box
            pos = [position(1)-position(3)/2 position(2)-position(4)/2 position(3) position(4)];

            %make the polygon
            x = [pos(1) pos(1)+pos(3) pos(1)+pos(3) pos(1) pos(1)];
            y = [pos(2) pos(2) pos(2)+pos(4) pos(2)+pos(4) pos(2)];
        end

        function stats = getMeasurements(ROI)
            [x, y] = getPoly(ROI);
            im = double(get(ROI.imageHandle,'CData'));

            m = size(im,1);
            n = size(im,2);
            %Scale the polygon to match the size of the displayed image (in
            %the case that the displayed image is being upsampled to match
            %the screen resolution).
            x = x*n/ROI.imageHandle.XData(2);
            y = y*m/ROI.imageHandle.YData(2);

            mask = poly2mask(x,y,m,n);

            im = im(mask);

            stats.mean = mean(im);
            stats.STD = std(im);
            stats.min = min(im);
            stats.max = max(im);
            stats.mask = mask;
            stats.position = ROI.position;

        end

        function handlePropEvents(ROI,src,evnt)
            position = getPosition(ROI);
            newPosition(ROI,position);
        end

        function set.fixedAspectRatio(ROI,fixedAspectRatio)
            %Set the fixed aspect ratio property
            ROI.fixedAspectRatio = fixedAspectRatio;
            %update the context menu
            for i=1:length(ROI.menuHandles)
                if strcmp(get(ROI.menuHandles(i),'Label'),'Fix Aspect Ratio')
                    if fixedAspectRatio
                        set(ROI.menuHandles(i),'Check','on')
                    else
                        set(ROI.menuHandles(i),'Check','off')
                    end
                end
            end
        end

        function aspectRatio = get.aspectRatio(ROI)
            position = ROI.position;
            aspectRatio = position(3)/position(4);
        end

        function BB = getBoundingBox(ROI)
            %BB = [rowMin rowMax; colMin ColMax];
            lims = size(get(ROI.imageHandle,'CData'));

            BB=zeros(2);

            rowMin = floor(ROI.position(2) - ROI.position(4)/2);
            if rowMin < 1
                rowMin = 1;
            end
            colMin = floor(ROI.position(1) - ROI.position(3)/2);
            if colMin < 1
                colMin = 1;
            end
            rowMax = floor(ROI.position(2) + ROI.position(4)/2 - 1);
            if rowMax > lims(1)
                rowMax = lims(1);
            end
            colMax = floor(ROI.position(1) + ROI.position(3)/2 - 1);
            if colMax > lims(2)
                colMax = lims(2);
            end
            BB = [rowMin rowMax; colMin colMax];

        end

    end

end

function ButtonDownFunction(hObject,evnt,ROI,n)

%get the parent figure handle
fig = ROI.figureHandle;

%get the type of click
click = get(fig,'SelectionType');

if strcmp(click,'normal')
    %get the current button motion and button up functions of the figure
    WBMF_old = get(fig,'WindowButtonMotionFcn');
    WBUF_old = get(fig,'WindowButtonUpFcn');

    %set the new window button motion function and button up function of the figure
    fun = @(src,evnt) ButtonMotionFunction(src,evnt,ROI,n);
    fun2=@(src,evnt)  ButtonUpFunction(src,evnt,ROI,WBMF_old,WBUF_old);
    set(fig,'WindowButtonMotionFcn',fun,'WindowButtonUpFcn',fun2);
end
end

function ButtonMotionFunction(src,evnt,ROI,n)
cp = get(ROI.axesHandle,'CurrentPoint'); cp=[cp(1,1) cp(1,2)];

position = getPosition(ROI);

switch n
    case 2                                          %middle cross
        position(1) = cp(1); position(2) = cp(2);

    case 3                                          %top left corner

        %find the bottom edge
        bottom = position(2)+position(4)/2;
        %get the new height
        height = bottom - cp(2);
        if height>1
            cy = cp(2)+height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)+position(3)/2-height*ROI.aspectRatio;
            end
        end

        %find the right edge
        right = position(1)+position(3)/2;
        %find the new width
        width = right-cp(1);
        if width>1
            %find the new center
            cx = cp(1)+width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 4                                          %top right corner
        %find the bottom edge
        bottom = position(2)+position(4)/2;
        %get the new height
        height = bottom - cp(2);
        if height>1
            cy = cp(2)+height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)-position(3)/2+height*ROI.aspectRatio;
            end
        end

        %find the left edge
        left = position(1)-position(3)/2;
        %find the new width
        width = cp(1) - left;
        if width>1
            cx = cp(1)-width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 5                                          %bottom left corner
        %find the top edge
        top = position(2)-position(4)/2;
        %get the new height
        height = cp(2) - top;
        if height>1
            cy = cp(2)-height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)+position(3)/2+-height*ROI.aspectRatio;
            end
        end

        %find the right edge
        right = position(1)+position(3)/2;
        %find the new width
        width = right-cp(1);
        if width>1
            %find the new center
            cx = cp(1)+width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 6                                          %bottom right corner
        %find the top edge
        top = position(2)-position(4)/2;
        %get the new height
        height = cp(2) - top;
        if height>1
            cy = cp(2)-height/2;
            position(2) = cy;
            position(4) = height;
            %Adjust cp(1) if you want to fix the aspect ratio
            if ROI.fixedAspectRatio
                cp(1) = position(1)-position(3)/2+height*ROI.aspectRatio;
            end
        end

        %find the left edge
        left = position(1)-position(3)/2;
        %find the new width
        width = cp(1) - left;
        if width>1
            cx = cp(1)-width/2;
            position(1) = cx;
            position(3) = width;
        end

    case 7                                          %left
        %find the right edge
        right = position(1)+position(3)/2;
        %find the new width
        width = right-cp(1);
        if width>1
            %find the new center
            cx = cp(1)+width/2;
            position(1) = cx;
            position(3) = width;
        end

        if ROI.fixedAspectRatio
            position(4) = width/ROI.aspectRatio;
        end

    case 8                                          %right
        %find the left edge
        left = position(1)-position(3)/2;
        %find the new width
        width = cp(1) - left;
        if width>1
            cx = cp(1)-width/2;
            position(1) = cx;
            position(3) = width;
        end

        if ROI.fixedAspectRatio
            position(4) = width/ROI.aspectRatio;
        end

    case 9                                          %top
        %find the bottom edge
        bottom = position(2)+position(4)/2;
        %get the new height
        height = bottom - cp(2);
        if height>1
            cy = cp(2)+height/2;
            position(2) = cy;
            position(4) = height;
        end

        if ROI.fixedAspectRatio
            position(3) = height*ROI.aspectRatio;
        end

    case 10                                         %bottom
        %find the top edge
        top = position(2)-position(4)/2;
        %get the new height
        height = cp(2) - top;
        if height>1
            cy = cp(2)-height/2;
            position(2) = cy;
            position(4) = height;
        end

        if ROI.fixedAspectRatio
            position(3) = height*ROI.aspectRatio;
        end

end

newPosition(ROI,position);

end

function ButtonUpFunction(src,evnt,ROI,WBMF_old,WBUF_old)
fig = ROI.figureHandle;
set(fig,'WindowButtonMotionFcn',WBMF_old,'WindowButtonUpFcn',WBUF_old);
end

function contextMenuCallback(source,callbackdata,ROI, tool)

switch get(source,'Label')
    case 'Delete'
        delete(ROI);
    case 'Export stats'
        stats = getMeasurements(ROI);
        name = inputdlg('Enter variable name');
        if isempty(name); return; end
        name=name{1};
        assignin('base', name, stats)
    case 'Fix Aspect Ratio'
        switch get(source,'Checked')
            case 'off'
                ROI.fixedAspectRatio=true;
            case 'on'
                ROI.fixedAspectRatio=false;
        end
    case 'Hide Text'
        switch get(source,'Checked')
            case 'off'
                set(source,'Checked','on');
                ROI.textVisible = false;
            case 'on'
                set(source,'Checked','off');
                ROI.textVisible = true;
        end
    case 'poly2mask'
        mask = tool.getCurrentMaskSlice;

        [x, y] = getPoly(ROI);

        m = size(mask,1);
        n = size(mask,2);
        %Scale the polygon to match the size of the displayed image (in
        %the case that the displayed image is being upsampled to match
        %the screen resolution).
        x = x*n/ROI.imageHandle.XData(2);
        y = y*m/ROI.imageHandle.YData(2);

        masknew = poly2mask(x,y,m,n);
        combine = true;
        tool.setCurrentMaskSlice(masknew,combine);
        notify(tool,'maskChanged')
end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [F,V]=patchNanFix(F,V)

numPoints=size(V,1); %Get original number of entries in vertex list

logicValid=~any(isnan(V),2); %Logic for valid vertices
logicValid_F=all(logicValid(F),2); %Logic vor valid faces

F=F(logicValid_F,:); %Get valid faces
V=V(logicValid,:); %Get valid vertices

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function z = times( x, y, oper )

%   Disciplined convex programming information for TIMES:
%      In general, disciplined convex programs must obey the "no-product
%      rule" which states that two non-constant expressions cannot be
%      multiplied together. Under this rule, the following combinations
%      are allowed:
%         {constant} .* {non-constant}  {non-constant} .* {constant}
%      Furthermore, if the non-constant expression is nonlinear, then
%      the constant term must be real.
%
%      A lone exception to the no-product rule is made for quadratic
%      forms: two affine expressions may be multiplied together if the
%      result is convex or concave. For example, the construction
%         variable x(n)
%         x.*x  <= 1;
%      would be permitted because each element of x.*x is convex.
%
%   Disciplined geometric programming information for TIMES:
%      Both terms in a multiplication must have the same log-curvature,
%      so the following products are permitted:
%         {log-convex} .* {log-convex}  {log-concave} .* {log-concave}
%         {log-affine} .* {log-affine}
%      Note that log-affine expressions are both log-convex and
%      log-concave.
%
%   For vectors, matrices, and arrays, these rules are verified
%   indepdently for each element.

persistent P
if isempty( P ),
    P.map = cvx_remap( ...
        ... % Invalid combinations: negative of log-concave, gp term
        { { 'negative' },   { 'l_concave_', 'g_valid' } }, ...
        { { 'l_concave_', 'g_valid' }, { 'negative' }   }, ...
        ... % Left-multiply by constant
        { { 'constant' },   { 'affine' }     }, ...
        { { 'real' },       { 'valid'  }     }, ...
        ... % Right-multiply by constant
        { { 'affine' },     { 'constant' }   }, ...
        { { 'valid' },      { 'real'  }      }, ...
        ... % Geometric
        { { 'l_convex' } }, ...
        { { 'l_concave' } }, ...
        ... % Potential quadratic form
        { { 'affine', 'p_convex', 'n_concave' } }, ...
        [ 0, 0, 1, 1, 2, 2, 3, 3, 4 ] );
    P.funcs = { @times_l, @times_r, @times_g, @times_q };
    P.constant = [];
end
if nargin < 3, oper = '.*'; end
P.name = oper;
z = cvx_binary_op( P, x, y );

function z = times_l( x, y )
% constant .* something
z = bsxfun( @times, cvx_constant( x ).', cvx_basis( y ) );
z = cvx( size(z,2), z );

function z = times_r( x, y )
% constant .* something
z = bsxfun( @times, cvx_basis( x ), cvx_constant( y ).' );
z = cvx( size(z,2), z );

function z = times_g( x, y )
% geom .* geom
z = exp( plus_nc( log( x ), log( y ) ) );

function z = times_q( x, y )
% affine .* affine
% p_convex .* p_convex
% n_concave .* n_concave
xA = x.basis_;
yA = y.basis_;
mm = max( size( xA, 1 ), size( yA, 1 ) );
xA( end + 1 : mm, : ) = 0;
yA( end + 1 : mm, : ) = 0;
xB  = xA( 2 : end, : );
yB  = yA( 2 : end, : );
cyB = conj( yB );
alpha = sum( yB .* cyB, 1 );
alpha = sum( bsxfun( @times, xB, yB ), 1 ) ./ max( alpha, realmin );
if nnz( xB - bsxfun( @times, alpha, cyB ) > 2 * eps * sqrt( conj( xB ) .* xB ) ),
    cvx_throw( 'Disciplined convex programming error:\n    Invalid quadratic form(s): not a square.' );
elseif any( abs(imag(alpha)) > 2 * eps * abs(real(alpha)) ),
    cvx_throw( 'Disciplined convex programming error:\n    Invalid quadratic form(s): not real.' );
else
    if isreal( xB ),
        z = square( y );
    else
        z = square( abs( y ) );
    end
    z.basis_ = bsxfun( @times, alpha, z.basis_ );
    offset = xA(1,:) - alpha .* conj( yA(1,:) );
    if any( offset ),
        z = z + offset.' .* y;
    end
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.tx for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function varargout = euclideanMST(points)
%EUCLIDEANMST Build euclidean minimal spanning tree of a set of points.
%
%   EDGES = euclideanMST(POINTS)
%   POINTS is a [NxP] array, N being the number of points and P being the
%   dimension.
%   Result EDGES is a [Mx2] array, containing indices of each vertex for
%   each edges.
%
%   [EDGES DISTS] = euclideanMST(POINTS)
%   Also returns the lengths of edges computed by MST algorithm.
%
%   Algorithm first computes Delaunay triangulation of the set of points,
%   then computes euclidean length of each edge of triangulation, and
%   finally uses prim algorithm to simplify the graph.
%
%   Example
%     % choose random points in the plane and display their Euclidean MST
%     pts = rand(50, 2)*100;
%     edges = euclideanMST(pts);
%     drawGraph(pts, edges)
%
%   See also
%   prim_mst, distancePoints, delaunayn
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2007-07-27, using Matlab 7.4.0.287 (R2007a)
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

% dimension
D   = size(points, 2);
Df  = factorial(D);

% compute all couples of vertices in unit triangle, tetrahedron, or n-dim
% simplex
subs = zeros(Df, 2);
k = 1;
for i = 1:D
    for j = i+1:D+1
        subs(k, 1) = i;
        subs(k, 2) = j;
        k = k + 1;
    end
end

% compute delaunay triangulation in D dimensions
tri = delaunayn(points);
Nt  = size(tri, 1);

% compute all possible edges
edges = zeros(Nt*Df, 2);
for t = 1:Nt
    for i = 1:Df
        edges((t-1)*Df+i, 1) = tri(t, subs(i, 1));
        edges((t-1)*Df+i, 2) = tri(t, subs(i, 2));
    end
end

% simplify edges
edges = unique(sort(edges, 2), 'rows');

% compute euclidean length of each edge
val = zeros(size(edges, 1), 1);
for i = 1:size(edges,1)
    val(i) = distancePoints(points(edges(i,1), :), points(edges(i,2), :));
end

% compute MST of created graph
[edges2, vals2] = prim_mst(edges, val);

% process output arguments
if nargout == 1
    varargout{1} = edges2;
elseif nargout==2
    varargout{1} = edges2;
    varargout{2} = vals2;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
clear all, close all, clc

% x = ones(500,1);
% t = .1*(1:500)';
% x = cos(t*2*pi);
% y = fft(x);
%
% N = length(x);
%
% plot(real(y),'k')
% hold on
% plot(real(DFT*x),'r--')

%% slow and steady
clear all, close all, clc
N = 256;
w = exp(-i*2*pi/N);
for i=1:N
    for j=1:N
        DFT(i,j) = w^((i-1)*(j-1));
    end
end

imagesc(real(DFT))
colormap jet
axis off
axis equal
box on
colorbar
set(gcf,'Position',[100 100 1000 900])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/DFTMatrixreal');
% print('-dpng', '-loose', '../figures/DFTMatrixreal');

%%

imagesc(imag(DFT))
colormap jet
axis off
axis equal
box on
colorbar
set(gcf,'Position',[100 100 1000 900])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/DFTMatriximag');
% print('-dpng', '-loose', '../figures/DFTMatriximag');
%% really fast
[I,J] = meshgrid(1:N,1:N);
DFT = w.^((I-1).*(J-1));

author: dynamicslab


text:
function tc_setAxisBounds(tc,newBounds);
%
% tc_setAxisBounds(tc,[newBounds]):
%
% Set the axis bounds for the time Course UI plots.
%
% If not provided,  puts up a dialog.
%
%
% ras 03/05.
if ieNotDefined('tc')
    tc = get(gcf,'UserData');
end

if ieNotDefined('newBounds')
    % put up a dialog
    def = {num2str(axis)};
    prompt{1} = 'New Axis Bounds [xmin xmax ymin ymax]: ';
	dlgTitle='Set Axis Bounds';
	answer=inputdlg(prompt,dlgTitle,1,def);
    newBounds = str2num(answer{1});
elseif isequal(newBounds,'auto')
    % set empty, so that will always autoscale
    newBounds = [];
end

tc.params.axisBounds = newBounds;

set(tc.ui.fig,'UserData',tc);
timeCourseUI;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http ://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

function A = buildDerivXMatrix(grid_size)

% d/dx for every entry in the first slice except the last column.
m = grid_size(1) * (grid_size(2) - 1);
n = grid_size(1) * grid_size(2);
e = ones(m, 1);
d_dx = spdiags([-e, e], [0, grid_size(1)], m, n);

A = sparse(0, 0);

for v = 1:grid_size(5)
    for u = 1:grid_size(4)
        for k = 1:grid_size(3)
            A = blkdiag(A, d_dx);
        end
    end
end

author: mahmoudnafifi


text:
%% Copyright (C) 2014, 2015, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

function z = numeric_array_to_sym(A)
%private helper for sym ctor
%   convert an array to syms, currently on 1D, 2D.

  [n, m] = size(A);

  if (n == 0 || m == 0)
    cmd = { sprintf('return sp.Matrix(%d, %d, []),', n, m) };
    z = pycall_sympy__ (cmd);
    return
  end

  Ac = cell(n,1);
  for i=1:n
    % we want all sym creation to go through the ctor.
    Ac{i} = cell(m,1);
    for j=1:m
      Ac{i}{j} = sym(A(i,j));
    end
  end

  %Ac = {{x 2}; {3 4}; {8 9}};

  d = size(A);
  if (length(d) > 2)
    error('conversion not supported for arrays of dim > 2');
  end

  cmd = { 'L = _ins[0]'
          'M = sp.Matrix(L)'
          'return M,' };
  z = pycall_sympy__ (cmd, Ac);

repo: octsympy


text:
classdef K3Supervisor < simiam.controller.Supervisor
%% SUPERVISOR switches between controllers and handles their inputs/outputs.
%
% Properties:
%   current_controller      - Currently selected controller
%   controllers             - List of available controllers
%   goal_points             - Set of goal points
%   goal_index              - Pointer to current goal point
%   v                       - Robot velocity
%
% Methods:
%   execute - Selects and executes the current controller.

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties
    %% PROPERTIES

        states
        eventsd

        current_state

        prev_ticks          % Previous tick count on the left and right wheels
        v
        goal

        goal_prev
        d_stop
        d_at_obs
        d_unsafe
        d_prog
        p

        direction

        v_gtg
        v_ao
        v_fw

%         is_init
    end

    methods
    %% METHODS

        function obj = K3Supervisor()
        %% SUPERVISOR Constructor
            obj = obj@simiam.controller.Supervisor();

            % initialize the controllers
            obj.controllers{1} = simiam.controller.AvoidObstacles();
            obj.controllers{2} = simiam.controller.GoToGoal();
            obj.controllers{3} = simiam.controller.GoToAngle();
            obj.controllers{4} = simiam.controller.AOandGTG();
            obj.controllers{5} = simiam.controller.Stop();
            obj.controllers{6} = simiam.controller.FollowWall();
            obj.controllers{7} = simiam.controller.SlidingMode();

            % set the initial controller
            obj.current_controller = obj.controllers{4};
            obj.current_state = 4;

            % generate the set of states
            for i = 1:length(obj.controllers)
                obj.states{i} = struct('state', obj.controllers{i}.type, ...
                                       'controller', obj.controllers{i});
            end

            % define the set of eventsd
            obj.eventsd{1} = struct('event', 'at_obstacle', ...
                                   'callback', @at_obstacle);

            obj.eventsd{2} = struct('event', 'at_goal', ...
                                   'callback', @at_goal);

            obj.eventsd{3} = struct('event', 'obstacle_cleared', ...
                                    'callback', @obstacle_cleared);

            obj.eventsd{4} = struct('event', 'unsafe', ...
                                    'callback', @unsafe);

            obj.eventsd{5} = struct('event', 'progress_made', ...
                                    'callback', @progress_made);

            obj.eventsd{6} = struct('event', 'sliding_left', ...
                                    'callback', @sliding_left);

            obj.eventsd{7} = struct('event', 'sliding_right', ...
                                    'callback', @sliding_right);

            obj.prev_ticks = struct('left', 0, 'right', 0);

            obj.v               = 0.15;

            obj.goal            = [-1;1];
            obj.goal_prev       = obj.goal;
            obj.d_stop          = 0.05;
            obj.d_at_obs        = 0.10;
            obj.d_unsafe        = 0.03;

            obj.d_prog = 10;

            obj.p = simiam.util.Plotter();
            obj.current_controller.p = obj.p;
        end

        function execute(obj, dt)
        %% EXECUTE Selects and executes the current controller.
        %   execute(obj, robot) will select a controller from the list of
        %   available controllers and execute it.
        %
        %   See also controller/execute

%             if(~obj.is_init)
%                 hold(obj.robot.parent, 'on');
%                 obj.v_gtg = plot(obj.robot.parent, [0 0], [0 0], 'b-');
%                 obj.v_ao = plot(obj.robot.parent, [0 0], [0 0], 'r-');
%                 obj.v_fw = plot(obj.robot.parent, [0 0], [0 0], 'g-');
%                 obj.is_init = true;
%             end

            inputs = obj.controllers{4}.inputs;
            inputs.x_g = obj.goal(1);
            inputs.y_g = obj.goal(2);
            inputs.v = obj.v;

%             if(~all(obj.goal==obj.goal_prev))
% %                 [x,y,theta] = obj.state_estimate.unpack();
% %                 if(norm(obj.goal-[x;y])<norm(obj.goal_prev-[x;y]))
%                     obj.set_progress_point();
% %                 end
%                 obj.goal_prev = obj.goal;
%                 obj.switch_to_state('go_to_goal');
%             end
%
            if (obj.check_event('at_goal'))
                obj.switch_to_state('stop');

%                 [x,y,theta] = obj.state_estimate.unpack();
%                 fprintf('stopped at (%0.3f,%0.3f)\n', x, y);
            else
                obj.switch_to_state('ao_and_gtg');
            end
%             elseif(obj.check_event('unsafe'))
%                 obj.switch_to_state('avoid_obstacles');
%             else
%                 if (obj.is_in_state('go_to_goal'))
%                     if(obj.check_event('at_obstacle') && obj.check_event('sliding_left'))
%                         obj.direction = 'left';
% %                         fprintf('now following to the left\n');
%                         obj.switch_to_state('follow_wall');
%                         obj.set_progress_point();
%                     elseif(obj.check_event('at_obstacle') && obj.check_event('sliding_right'))
%                         obj.direction = 'right';
% %                         fprintf('now following to the right\n');
%                         obj.switch_to_state('follow_wall');
%                         obj.set_progress_point();
%                     end
%                 elseif (obj.is_in_state('follow_wall') && strcmp(obj.direction,'left'))
%                     if(obj.check_event('progress_made') && ~obj.check_event('sliding_left'))
%                         obj.switch_to_state('go_to_goal');
%                     end
%                 elseif (obj.is_in_state('follow_wall') && strcmp(obj.direction, 'right'))
%                     if(obj.check_event('progress_made') && ~obj.check_event('sliding_right'))
%                         obj.switch_to_state('go_to_goal');
%                     end
%                 elseif (obj.is_in_state('avoid_obstacles'))
%                     if(obj.check_event('obstacle_cleared'))
% %                         if(obj.check_event('sliding_left'))
% %                             obj.direction = 'left';
% %                             obj.switch_to_state('follow_wall');
% %                         elseif(obj.check_event('sliding_right'))
% %                             obj.direction = 'right';
% %                             obj.switch_to_state('follow_wall');
% %                         else
%                             obj.switch_to_state('go_to_goal');
% %                         end
%                     end
%                 end
%
%             end

%             elseif(obj.check_event('unsafe'))
%                 obj.switch_to_state('avoid_obstacles');
%             else
%                 if (obj.is_in_state('go_to_goal'))
%                     if(obj.check_event('at_obstacle') && obj.check_event('sliding_left'))
%                         obj.direction = 'left';
% %                         fprintf('now following to the left\n');
%                         obj.switch_to_state('follow_wall');
%                         obj.set_progress_point();
%                     elseif(obj.check_event('at_obstacle') && obj.check_event('sliding_right'))
%                         obj.direction = 'right';
% %                         fprintf('now following to the right\n');
%                         obj.switch_to_state('follow_wall');
%                         obj.set_progress_point();
%                     end
%                 elseif (obj.is_in_state('follow_wall') && strcmp(obj.direction,'left'))
%                     if(obj.check_event('progress_made') && ~obj.check_event('sliding_left'))
%                         obj.switch_to_state('go_to_goal');
%                     end
%                 elseif (obj.is_in_state('follow_wall') && strcmp(obj.direction, 'right'))
%                     if(obj.check_event('progress_made') && ~obj.check_event('sliding_right'))
%                         obj.switch_to_state('go_to_goal');
%                     end
%                 elseif (obj.is_in_state('avoid_obstacles'))
%                     if(obj.check_event('obstacle_cleared'))
% %                         if(obj.check_event('sliding_left'))
% %                             obj.direction = 'left';
% %                             obj.switch_to_state('follow_wall');
% %                         elseif(obj.check_event('sliding_right'))
% %                             obj.direction = 'right';
% %                             obj.switch_to_state('follow_wall');
% %                         else
%                             obj.switch_to_state('go_to_goal');
% %                         end
%                     end
%                 end

%             end

%             inputs.direction = obj.direction;

            outputs = obj.current_controller.execute(obj.robot, obj.state_estimate, inputs, dt);

            [vel_r, vel_l] = obj.robot.dynamics.uni_to_diff(outputs.v, outputs.w);
            obj.robot.set_wheel_speeds(vel_r, vel_l);

%             fprintf('(v,w) = (%0.3f,%0.3f)\n', outputs.v, outputs.w);

            obj.update_odometry();
%             [x, y, theta] = obj.state_estimate.unpack();
%             fprintf('current_pose: (%0.3f,%0.3f,%0.3f)\n', x, y, theta);
        end

        function set_progress_point(obj)
            [x, y, theta] = obj.state_estimate.unpack();
            obj.d_prog = norm([x-obj.goal(1);y-obj.goal(2)]);
        end

        %% Events %%

        function rc = sliding_left(obj, state, robot)
            inputs = obj.controllers{7}.inputs;
            inputs.x_g = obj.goal(1);
            inputs.y_g = obj.goal(2);
            inputs.direction = 'left';

            obj.controllers{7}.execute(obj.robot, obj.state_estimate, inputs, 0);

            u_gtg = obj.controllers{7}.u_gtg;
            u_ao = obj.controllers{7}.u_ao;
            u_fw = obj.controllers{7}.u_fw;

            %% START CODE BLOCK %%
            sigma = [u_gtg u_ao]\u_fw;
            %% END CODE BLOCK %%

            rc = false;
            if sigma(1) > 0 && sigma(2) > 0
%                 fprintf('now sliding left\n');
                rc = true;
            end
        end

        function rc = sliding_right(obj, state, robot)
            inputs = obj.controllers{7}.inputs;
            inputs.x_g = obj.goal(1);
            inputs.y_g = obj.goal(2);
            inputs.direction = 'right';

            obj.controllers{7}.execute(obj.robot, obj.state_estimate, inputs, 0);

            u_gtg = obj.controllers{7}.u_gtg;
            u_ao = obj.controllers{7}.u_ao;
            u_fw = obj.controllers{7}.u_fw;

            %% START CODE BLOCK
            sigma = [u_gtg u_ao]\u_fw;
            %% END CODE BLOCK

            rc = false;
            if sigma(1) > 0 && sigma(2) > 0
%                 fprintf('now sliding right\n');
                rc = true;
            end
        end

        function rc = at_obstacle(obj, state, robot)
            ir_distances = obj.robot.get_ir_distances();
            rc = false;                                     % Assume initially that the robot is clear of obstacle

            % Loop through and test the sensors (only the front set)
            if any(ir_distances(2:7) < obj.d_at_obs)
                rc = true;
            end
        end

        function rc = unsafe(obj, state, robot)
            ir_distances = obj.robot.get_ir_distances();
            rc = false;             % Assume initially that the robot is clear of obstacle

            % Loop through and test the sensors (only the front set)
            if any(ir_distances(2:7) < obj.d_unsafe)
                    rc = true;
            end
        end

        function rc = at_goal(obj, state, robot)
            [x,y,theta] = obj.state_estimate.unpack();
            rc = false;

            % Test distance from goal
            if norm([x - obj.goal(1); y - obj.goal(2)]) < obj.d_stop
                rc = true;
            end
        end

        function rc = obstacle_cleared(obj, state, robot)
            ir_distances = obj.robot.get_ir_distances();
            rc = false;              % Assume initially that the robot is clear of obstacle

            % Loop through and test the sensors (only front set)
            if all(ir_distances(2:7) > obj.d_at_obs)
                rc = true;
            end
        end

        function rc = progress_made(obj, state, robot)

            % Check for any progress
            [x, y, theta] = obj.state_estimate.unpack();
            epsilon = 0.1;

            %% START CODE BLOCK %%
            rc = false;
            if norm([x-obj.goal(1); y-obj.goal(2)]) < (obj.d_prog - epsilon)
                rc = true;          % progress has been made
            end
            %% END CODE BLOCK %%

        end

        %% State machine support functions

        function set_current_controller(obj, ctrl)
            % save plots
            obj.current_controller = ctrl;
            obj.p.switch_2d_ref();
            obj.current_controller.p = obj.p;
        end

        function rc = is_in_state(obj, name)
            rc = strcmp(name, obj.states{obj.current_state}.state);
        end

        function switch_to_state(obj, name)

            if(~obj.is_in_state(name))
                for i=1:numel(obj.states)
                    if(strcmp(obj.states{i}.state, name))
                        obj.set_current_controller(obj.states{i}.controller);
                        obj.current_state = i;
                        fprintf('switching to state %s\n', name);
                        return;
                    end
                end
            else
%                 fprintf('already in state %s\n', name);
                return
            end

            fprintf('no state exists with name %s\n', name);
        end

        function rc = check_event(obj, name)
           for i=1:numel(obj.eventsd)
               if(strcmp(obj.eventsd{i}.event, name))
                   rc = obj.eventsd{i}.callback(obj, obj.states{obj.current_state}, obj.robot);
                   return;
               end
           end

           % return code (rc)
           fprintf('no event exists with name %s\n', name);
           rc = false;
        end

        %% Odometry

        function update_odometry(obj)
        %% UPDATE_ODOMETRY Approximates the location of the robot.
        %   obj = obj.update_odometry(robot) should be called from the
        %   execute function every iteration. The location of the robot is
        %   updated based on the difference to the previous wheel encoder
        %   ticks. This is only an approximation.
        %
        %   state_estimate is updated with the new location and the
        %   measured wheel encoder tick counts are stored in prev_ticks.

            % Get wheel encoder ticks from the robot
            right_ticks = obj.robot.encoders(1).ticks;
            left_ticks = obj.robot.encoders(2).ticks;

            % Recal the previous wheel encoder ticks
            prev_right_ticks = obj.prev_ticks.right;
            prev_left_ticks = obj.prev_ticks.left;

            % Previous estimate
            [x, y, theta] = obj.state_estimate.unpack();

            % Compute odometry here
            R = obj.robot.wheel_radius;
            L = obj.robot.wheel_base_length;
            m_per_tick = (2*pi*R)/obj.robot.encoders(1).ticks_per_rev;

            d_right = (right_ticks-prev_right_ticks)*m_per_tick;
            d_left = (left_ticks-prev_left_ticks)*m_per_tick;
            d_center = (d_right + d_left)/2;

            x_dt = d_center*cos(theta);
            y_dt = d_center*sin(theta);
            theta_dt = (d_right - d_left)/L;

            theta_new = theta + theta_dt;
            x_new = x + x_dt;
            y_new = y + y_dt;
%             fprintf('Estimated pose (x,y,theta): (%0.3g,%0.3g,%0.3g)\n', x_new, y_new, theta_new);

            % Save the wheel encoder ticks for the next estimate
            obj.prev_ticks.right = right_ticks;
            obj.prev_ticks.left = left_ticks;

            % Update your estimate of (x,y,theta)
            obj.state_estimate.set_pose([x_new, y_new, theta_new]);
        end
    end
end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function reportResults(status)
global logger
global testidx
global procStreamStyle

logger.Write('\n');

for iTest = 1:length(status)
    if status(iTest)~=0
        logger.Write(sprintf('Unit test #%d FAILED.\n', iTest));
    else
        logger.Write(sprintf('Unit test #%d PASSED.\n', iTest));
    end
end
logger.Write('\n');

testidx=[];
procStreamStyle=[];

author: BUNPC


text:
function [v]=tetVolMeanEst(F,V)

% function [v]=tetVolMeanEst(F,V)
% ------------------------------------------------------------------------
%
% This function calculates the volume of an ideal regular tetrahedron with
% edge lengths (all equal) that match the mean edge lengths occuring for
% the input surface defined by F (faces) and V (vertices).

%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/10/17
%------------------------------------------------------------------------
%%

[edgeLengths]=patchEdgeLengths(F,V);
edgeLengthsMean=mean(edgeLengths);
meanProposedVolume=edgeLengthsMean^3./(6*sqrt(2)); %For a regular tetrahedron
v=meanProposedVolume;

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


text:
function signal = logCompression(signal, a, normalise)
%LOGCOMPRESSION   Log compress an input signal.
%
% DESCRIPTION:
%       logCompression compresses the input signal using the expression
%       signal = log10(1 + a*signal)./log10(1 + a)
%
% USAGE:
%       signal = logCompression(signal, a)
%       signal = logCompression(signal, a, normalise)
%
% INPUTS:
%       signal      - input signal
%       a           - compression factor
%
% OPTIONAL INPUTS
%       normalise   - Boolean controlling whether the maximum of the input
%                     signal is normalised to unity before compression
%                     (default = false). If set to true, the original
%                     magnitude is restored after compression.
%
% OUTPUTS:
%       signal      - log compressed signal
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 24th February 2011
%       last update - 24th February 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% check for optional normalise input
if nargin == 2
    normalise = false;
end

% compress signal
if normalise
    mx = max(signal(:));
    signal = mx*(log10(1 + a*signal./mx)./log10(1 + a));
else
    signal = log10(1 + a*signal)./log10(1 + a);
end

repo: Aero-matlab


text:
function s = plus(f, g)
%+   Addition of DELTAFUN objects.
%   F + G adds F and G, where F and G may be DELTAFUN objects or scalars.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%

% If one of the arguments is empty:
if ( isempty(f) || isempty(g) )
    % Return with an empty DELTAFUN:
    s = deltafun;
    return
end

% One of the arguments (i.e., f or g) is necessarily a DELTAFUN object.
% (Otherwise, this overloaded plus would not have been called!). Ensure it is f:
if ( ~isa(f, 'deltafun') )
    s = plus(g, f);
    return
end

% Cast g to a DELTAFUN if it is a double or a different kind of FUN:
if ( isa(g, 'double') )
    g = deltafun(0*f.funPart + g);
elseif ( isa(g, 'fun') )
    g = deltafun(g);
end

% Add the funParts:
funPart = f.funPart + g.funPart;

% Add the delta functions:
[deltaMag, deltaLoc] = deltafun.mergeDeltas(f.deltaMag, f.deltaLoc, ...
    g.deltaMag, g.deltaLoc);

% Assemble the output DELTAFUN:
s = deltafun(funPart, struct('deltaMag', deltaMag, 'deltaLoc', deltaLoc));

% Simplify:
s = simplifyDeltas(s);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_ft_connectivitysimulation

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_connectivitysimulation

cfg = [];
cfg.method      = 'ar';
cfg.nsignal     = 32;
cfg.ntrials     = 100;
cfg.triallength = 2;
cfg.fsample     = 500;

% test 'ar' method
cfg.params = zeros(cfg.nsignal,cfg.nsignal,cfg.ntrials);
cfg.params(14,21,42) = 0.5;
cfg.noisecov = 0.5*eye(cfg.nsignal);
dataout = ft_connectivitysimulation(cfg);

return
%% test linear_mix method => wait for issue to be solved
cfg.mix = 0.5*ones(32,2);
cfg.delay = 50*ones(32,2);

dataout = ft_connectivitysimulation(cfg);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

function data = getSimulatorBlockUserData(modelHandle)

blockName = find_system(modelHandle,...
    'SearchDepth', 1, ...
    'MaskType', 'Simulator');

if numel(blockName) == 0
    error('Cannot find Simulator Block at top level of model');
elseif numel(blockName) > 1
    error('Multiple Simulator Blocks');
end

data = get_param(blockName, 'UserData');
if isempty(data)
   error('Simulator Block is not initialized');
end


text:
function c = tapas_ehgf_binary_pu_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF)
% for binary inputs in the *presence* of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_hgf_config.
%
% This file refers to UNCERTAIN inputs (Eqs 45-47 in Mathys et al., (2011));
% for inputs without uncertainty, refer to tapas_hgf_binary_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_binary_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.rho        row vector of rhos (representing drift; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%         est.p_prc.al         scalar alpha (perceptual uncertainty)
%         est.p_prc.eta0       scalar eta0 (mean of first input category)
%         est.p_prc.eta1       scalar eta1 (mean of second input category)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_hgf_binary_pu_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - When analyzing a new dataset, take your inputs u and use
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'ehgf_binary_pu';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          0,      0];

% Rhos
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero to turn off drift.
c.rhomu = [NaN, 0, 0];
c.rhosa = [NaN, 0, 0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -3,  -6];
c.omsa = [NaN, 4^2, 4^2];

% Alpha
% Format: scalar.
c.logalmu = log(0.5);
c.logalsa = 1;

% Eta0
% Format: scalar.
c.eta0mu = 0;
c.eta0sa = 0;

% Eta1
% Format: scalar.
c.eta1mu = 1;
c.eta1sa = 0;

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.rhomu,...
    c.logkamu,...
    c.ommu,...
    c.logalmu,...
    c.eta0mu,...
    c.eta1mu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.rhosa,...
    c.logkasa,...
    c.omsa,...
    c.logalsa,...
    c.eta0sa,...
    c.eta1sa,...
         ];

% Check whether we have the right number of priors
expectedLength = 3*c.n_levels+2*(c.n_levels-1)+4;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_ehgf_binary_pu;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_ehgf_binary_pu_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function [d] = spm_kl_normal (m_q,c_q,m_p,c_p)
% KL divergence between two multivariate normal densities
% FORMAT [d] = spm_kl_normal (m_q,c_q,m_p,c_p)
%
% KL (Q||P) = <log Q/P> where avg is wrt Q
%
% between two Normal densities Q and P
%
% m_q, c_q    Mean and covariance of first Normal density
% m_p, c_p    Mean and covariance of second Normal density
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_kl_normal.m 2696 2009-02-05 20:29:48Z guillaume $

d=length(m_q);
m_q=m_q(:);
m_p=m_p(:);

Term1=0.5*spm_logdet(c_p)-0.5*spm_logdet(c_q);

inv_c_p=inv(c_p);
Term2=0.5*trace(inv_c_p*c_q)+0.5*(m_q-m_p)'*inv_c_p*(m_q-m_p);

d=Term1+Term2-0.5*d;

author: spm


text:
function c=ref_dstiv(f)
%REF_DSTIV  Reference Discrete Sine Transform type IV
%   Usage:  c=ref_dstiv(f);
%
%

L=size(f,1);
W=size(f,2);

% Create weights.
w=sqrt(2/L);

% Create transform matrix.
F=zeros(L);

for m=0:L-1
  for n=0:L-1
    F(m+1,n+1)=w*sin(pi*(n+.5)*(m+.5)/L);
  end;
end;

% Compute coefficients.
c=F*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function values = coeffs2vals(coeffs)
%COEFFS2VALS   Convert Fourier coefficients to values at N equally spaced
%points between [-1 1), where N is the number of coefficients.
%   V = COEFFS2VALS(C) returns the values of the trignometric polynomial
%   as follows:
%   If N is odd
%       F(x) = C(1)*z^(-(N-1)/2) + C(2)*z^(-(N-1)/2-1) + ... + C(N)*z^((N-1)/2)
%   If N is even
%       F(x) = C(1)*z^(-N/2) + C(2)*z^(-N/2+1) + ... + C(N)*z^(N/2-1)
%   where z = exp(1i*pi*x) and -1 <= x <= 1.
%
%   If the input C is an (N+1)xM matrix then V = COEFFS2VALS(C) returns the
%   (N+1)xM matrix of values V such that V(i,j) is the ith value of the
%   trignometric polynomial corresponding to the jth column.
%
% See also VALS2COEFFS, TRIGPTS.

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

% *Note about symmetry*.  Some of the code below is designed to
% enforce two symmetries whose failure might disturb users:
% COEFFS exactly real ==> VALUES exactly hermitian
% COEFFS exactly imaginary ==> VALUES exactly skew-hermitian
% This is necessary because the MATLAB FFT code does not
% exactly preserve symmetries.

% Get the length of the input:
n = size(coeffs, 1);

% Trivial case (constant or empty):
if ( n <= 1 )
    values = coeffs;
    return
end

% The coefficients are for interpolation defined on [-pi,pi), but the FFT
% works for values on [0,2*pi). To fix the coefficients for this we just need to
% assign c_k = (-1)^k c_k, with k=-(N-1)/2:(N-1)/2 for N odd, and
% k = -N/2:N/2-1 for N even.
if ( mod(n, 2) )
    even_odd_fix = (-1).^(-(n-1)/2:(n-1)/2).';
else
    even_odd_fix = (-1).^((-n/2):(n/2-1)).';
end
coeffs = bsxfun(@times, coeffs, even_odd_fix);

% test for symmetry
isHerm = max(abs(imag(coeffs)),[],1) == 0;
isSkew = max(abs(real(coeffs)),[],1) == 0;

% Shift the coefficients properly.
values = ifft(ifftshift(n*coeffs, 1), [], 1);

% correct if symmetric
vals = [values;values(1,:)];
hermvals = (vals+flipud(conj(vals)))/2;
skewvals = (vals-flipud(conj(vals)))/2;
values(:,isHerm) = hermvals(1:end-1,isHerm);
values(:,isSkew) = skewvals(1:end-1,isSkew);

end


text:
% timefdetails() - details of the timef() function for time/frequency analysis
%                  of multiple epochs of single-channel event-related data.
%
% Global Description:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% timef() performs normalized time/frequency averaging using either
% FFT-, wavelet-, or multitaper DFT estimates. The wavelets are N-cycle
% Hanning-windowed sinusoids. (Note: To substitute for hanning() windowing
% gauss() or other windowing, replace the timef.m reference to hanning()).
%
% By default, the two image panels of the output plot show, respectively, the
% event-related spectral perturbation (ERSP) and inter-trial coherence (ITC)
% of the input data.
%
% The ERSP:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The ERSP (S. Makeig, Electroencephalogr Clin Neurophysiol 86:283-93, 1993)
% shows mean log event-loced deviations from epoch-mean (or baseline-mean) power
% at each frequency. If bootstrap statistics are computed (via the 'alpha'
% probability input parameter), (time, freq) points with non-significant
% differences from 0 are colored green in the image (but not in the 'ersp'
% output variable - use the output 'erspboot' variable to re-mask the 'ersp'
% output if desired). The baseline mean spectrum removed from each epoch
% is available as output parameter "powbase". Note that log(power) differences
% are equivalent to log(power ratios) between baseline and observed spectra,
% meaning the implicit ERSP model is one of (multiplicative) amplitude modulation
% of the EEG/MEG spectrum by e.g. subcortical and/or intra-cortical influences.
%
% In the default view, the thin bottom panel below the upper (ERSP) image shows
% the ERSP envelope (the most positive and most negative values at each output
% time point). The thin left panel shows the mean (or baseline) log spectrum
% (blue trace). When bootstrap statistics are computed (via the "alpha" argument),
% the left panel (green trace) also shows the bootstrap significant levels (+/-)
% at each frequency.
%
% The ITC:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (Inter-trial Coherence, cf. Tallon-Baudry et al., "Phase-locking factor")
% The lower panel shows the degree of tendency for spectral phase at each
% (time, freq) point to repeat across the data epochs. If bootstrap statistics
% are computed (as per the 'alpha' input parameter), non-significant points
% are colored green (but again not in the 'itc' output variable - use the
% itcboot output to re-mask if desired in later plotting).
%
% The lower thin panel shows the time-domain average (ERP) of the input data
% (blue) plus a zero-line (green). The average (ERP) is created principally by
% partial phase resetting of the EEG as measured by the ITC. (While phase resetting
% dominates, event-related spectral power changes (as measured by the ERSP)
% may also play a minor role). The thin left panel shows the frequency-mean ITC
% (blue trace) and, if bootstrap statistics are computed, the ITC significance
% limits at each frequency (green trace).
%
% ITC Math Derivation:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% By definition, linear coherence is
%       R=mean(Fxy)/sqrt(mean(abs(Fxx))*mean(abs(Fyy)));
% where Fxy is the cross-spectrum (FxFy*) and Fxx and Fyy the autospectra of
% processes x and y.  We define the phase coherence to be
%       R=mean(Fxy/(abs(Fxx))*abs(Fyy));  % mean of individually normed Fxy
% To derive the ITC, we consider y to be a stimulus-locked process
% such that, at each time and frequency, angle(Fyy) = 0 and abs(Fyy) = 1.
% Thus Pxy = Pxx, and the (complex) inter-trial phase coherence between x and
% the constant stimulus-locked process y is
%       ITC=mean(Pxx/abs(Pxx));
%
% USAGE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   >> [ersp,itc,powbase,times,freqs,erspboot,itcboot]  ...
%                = timef(data,frames,tlimits,srate,cycles,...
%                              'key1',value1,'key2',value2, ... );
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE:
% * Left-click on subplots to view and zoom in separate windows (uses axcopy()).
%
% Required Inputs:
%   data        = Single-channel data vector (1,frames*ntrials) (required)
%
%   frames      = Frames per trial                        {750}
%  Here, a frame is a data point (one channel at one time point) and the data
%  is assumed to be composed of concatenated epochs of the same length (frames).
%
%   tlimits     = [mintime maxtime] (ms) Epoch time limits {[-1000 2000]}
%  These should be the starting and ending times of the input epochs.
%
%   srate       = data sampling rate (Hz)                 {250}
%  This is sample rate per channel (i.e., data frames per second).
%
%   cycles      = >0 -> Number of cycles in each analysis wavelet
%                 =0 -> Use FFTs (with constant window length) {0}
%                 If [wavecycles factor] -> wavelet cycles increase with frequency
%                 beginning at wavecyles (0<factor<1; factor=1 -> no increase,
%                 standard wavelets; factor=0 -> fixed epoch length, as in FFT.
%                 OR multitaper decomposition (with 'mtaper').
%  Here, the user chooses either to use the FFT method (fixed window size
%  for all frequencies) or the wavelet DFT method (the data window length
%  depends inversely on the frequency, e.g. '3' means that the data
%  windows will each be three cycles wide (at each frequency). A higher
%  number here (e.g., '5') will narrow the frequency band and widen the
%  time window.
%
% Optional Inter-Irial Coherence Type:
%   'type'      = ['coher'|'phasecoher'] Compute either linear coherence
%                  ('coher') or phase coherence ('phasecoher') also known
%                  as the phase coupling factor           {'phasecoher'}.
%
% Optional Detrending:
%   'detret'    = ['on'|'off'], Detrend data in time.       {'off'}
%   'detrep'    = ['on'|'off'], Detrend data across trials (at each time point
%                 compute the linear trend across the set of ordered trials and
%                 remove it; not that this also subtract the ERP) {'off'}
%
% Optional FFT/DFT Parameters:
%   'winsize'   = If cycles==0: data subwindow length (fastest is 2^n < frames);
%                 If cycles >0: *longest* window length to use. This determines
%                 the lowest output frequency  {default: ~frames/8}
%  When cycles>0, winsize determines the lowest computed frequency. For example,
%  with srate=100 and cycles=3, a winsize of 100 means that 3 cycles must fit
%  within a 1-sec ( 100-sample) window. So, the lowest output frequency is 3 Hz.
%
%  When cycles=0, winsize is the length of data in each FFT window. This may be
%  extended with zeroes (to give more output frequencies) using padratio (below).
%
%  'timesout'  = Number of output times (int<frames-winframes) {200}
%  The number of FFTs or wavelet DFTs computed and plotted.
%
%   'padratio'  = FFT-length/winframes (2^k)                    {2}
%                  Multiplies the number of output frequencies by
%                  dividing their spacing. When cycles==0, frequency
%                  spacing is (low_freq/padratio).
%  This factor multiplies the number of output frequencies. In the FFT method
%  (cycles=0), this is done by zero-padding each analysis window. In the wavelet
%  DFT method (cycles>0), this gives the number of frequencies per Hz.
%
%   'maxfreq'   = Maximum frequency (Hz) to plot (& to output, if cycles>0)
%                  If cycles==0, all FFT frequencies are output. {50}
%   'baseline'  = Spectral baseline end-time (in ms). Use NaN for no baseline
%                 removal{0}
%   'powbase'   = Baseline spectrum to log-subtract. 'baseline' parameter is
%                 ignored if this parameter is used {def|NaN->from data}
%  This is useful only when you want to use a known baseline spectrum (e.g. from
%  another condition) instead of using the actual mean baseline spectrum of the data.
%  Otherwise, leave this out or specify as 'NaN' (not a number).
%
% Optional Multitaper Parameters:
%   'mtaper'    = If [N W], performs multitaper decomposition.
%                  (N is the time resolution and W the frequency resolution;
%                  maximum taper number is 2NW-1). Overwrites 'winsize' and
%                  'padratio'.
%                  If [N W K], forces the use of K Slepian tapers (if possible).
%                  Phase is calculated using standard methods.
%                  The use of mutitaper with wavelets (cycles>0) is not
%                  recommended (as multiwavelets are not implemented).
%                  Uses Matlab functions DPSS, PMTM.   {no multitaper}
%
% Optional Bootstrap Parameters:
%   'alpha'     = If non-0, compute two-tailed bootstrap significance prob.
%                  level. Show non-signif. output values as green.   {0}
%  This optional parameter lengthens the computation time but gives bootstrap
%  estimates of which ERSP and ITC values are significantly different from 0,
%  by setting to 0 (green) all non-significant values in the ERSP and ITC images.
%  Normal values for alpha are 0 ([], or none) -> no bootstrap computation, or
%  0.01 (which should allow about 1% of random images to appear "significant"
%
%   'naccu'     = Number of bootstrap replications to accumulate     {200}
%   'baseboot'  = Bootstrap baseline subtract (0 -> use 'baseline';
%                                                  1 -> use whole trial) {0}
% Optional Scalp Map Plotting Parameters:
%   'topovec'   = Scalp topography (map) to plot                     {none}
%   'elocs'     = Electrode location file for scalp map   {no default}
%                     File should be ascii in format of  >> topoplot example
%  This is an optional map-plotting feature. Given an input map vector
%  (one weight at each channel, and a electrode location file, timef() plots
%  a topoplot()-style 2-d scalp map on the left side of the figure. See
%  >> topoplot example % for the format of the electrode location file.
%
% Other Optional Plotting Parameters:
%   'vert'      = [vector of ms times] -> plot vertical dashed lines  {0 only}
%  Use this to add extra vertical dashed lines at significant epoch times.
%  Time 0 is marked by default.
%
%   'plotersp'  = ['on'|'off'] Plot power spectral perturbations    {'on'}
%   'plotitc'   = ['on'|'off'] Plot inter trial coherence            {'on'}
%   'title'     = Optional figure title                              {none}
%
%   'pboot'     = Bootstrap power limits (e.g., from timef())   {from data}
%   'rboot'     = Bootstrap ITC limits (e.g., from timef())     {from data}
%  These are useful if you want to apply significance limits from another condition
%  to new data. {default|NaN, compute from data}
%
%   'linewidth' = Line width for 'marktimes' traces (thick=2, thin=1) {2}
%   'axesfont'  = Axes text font size                                {10}
%   'titlefont' = Title text font size                               {8}
%
% Outputs:
%        ersp   = Matrix (nfreqs,timesout) of log spectral diffs. from baseline (dB)
%        itc    = Matrix (nfreqs,timesout) of inter-trial phase coherence (range: [0 1])
%   Note that when cycles=0, nfreqs is total number of FFT frequencies, which
%   typically include frequencies higher than maxfreq. When cycles>0, *no* extra
%   (higher) frequencies are computed.
%
%      powbase  = Baseline power spectrum (removed to compute the ERSP)
%        times  = Vector of output times (subwindow centers) (in ms).
%        freqs  = Vector of frequency bin centers (in Hz).
%
%     erspboot  = Matrix (2,nfreqs) of [lower;upper] ERSP significance diffs.
%      itcboot  = Matrix (2,nfreqs) of [lower;upper] ITC thresholds (abs., not diffs)
%  Note that the itcboot lower bound is practically meaningless.
%
%  Plot description:
%    Assuming both 'plotersp' and 'plotitc' options are 'on' (= default). The upper panel
%    presents the data ERSP (Event-Related Spectral Perturbation) in dB, with mean baseline
%    spectral activity (in dB) subtracted. Use "'baseline', NaN" to prevent timef() from
%    removing the baseline. The lower panel presents the data ITC (Inter-Trial Coherence).
%    Click on any plot axes to pop up a new window (using 'axcopy()')
%    -- Upper left marginal panel presents the mean spectrum during the baseline period
%       (blue), and when significance is set, the significance threshold at each frequency
%       (dotted green-black trace).
%    -- The marginal panel under the ERSP image shows the maximum (green) and minimum
%       (blue) ERSP values relative to baseline power at each frequency.
%    -- The lower left marginal panel shows mean ITC across the imaged time range (blue),
%       and when significance is set, the significance threshold (dotted green-black).
%    -- The marginal panel under the ITC image shows the ERP (which is produced by ITC
%       across the data spectral pass band).
%
% Author: Sigurd Enghoff, Arnaud Delorme & Scott Makeig
%          CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
% See also: crossf() - event-related cross-spectral coherence between two input
%                      time series.
%
% History:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% timef() was coded by Sigurd Enghoff and Scott Makeig at The Salk
% Institute, La Jolla CA in August, 1998, using methods developed
% in Makeig, 1993. Arno Delorme added the multitaper option, recoded
% the function to use 'keyword','parameter' argument pairs, and added the
% 'type' argument with advice from Joern Anemueller at SCCN/Institute for
% Neural Computation, UCSD in early 2002.

% Copyright (C) 8/01/00 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 01-25-02 reformated help & license -ad

author: PatternRecognition

repo: OpenBMI


text:
function [ vX, mX ] = SolveLsL1Admm( vX, mA, vY, paramLambda, sSolverParams )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveProxTvAdmm( vY, mD, paramLambda, numIterations )
% Solves the L1 Norm Regualrized Least Squares using Alternating
% Direction Method of Multipliers (ADMM) Method.
% Basically solves the problem given by:
% $$ \arg \min_{ x \in \mathbb{R}^{n} } \frac{1}{2} {\left\| A x - y \right|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} $$
% Input:
%   - vX                -   Optimization Vector.
%                           The vector to be optimized. Initialization of
%                           the iterative process.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vY                -   Measurements Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - paramRho          -   The Rho Parameter.
%                           Sets the weight of the equality constraint.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia ADMM - https://en.wikipedia.org/wiki/Augmented_Lagrangian_method#Alternating_direction_method_of_multipliers.
% Remarks:
%   1.  Using vanilla ADMM with no optimization of the parameter or
%       smoothing.
%   2.  For high values of `paramLambda` it will require many iterations.
% Known Issues:
%   1.  C
% TODO:
%   1.  D
% Release Notes:
%   -   1.0.000     05/08/2021  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    vX (:, 1) {mustBeFloat, mustBeReal}
    mA (:, :) {mustBeFloat, mustBeReal}
    vY (:, 1) {mustBeFloat, mustBeReal}
    paramLambda (1, 1) {mustBeFloat, mustBeReal, mustBePositive}
    sSolverParams.paramRho (1, 1) {mustBeNumeric, mustBeReal, mustBePositive} = 5
    sSolverParams.numIterations (1, 1) {mustBeNumeric, mustBeReal, mustBePositive, mustBeInteger} = 100
end

MAT_TYPE_SKINNY = 1;
MAT_TYPE_FAT    = 2;

if(size(mA, 1) >= size(mA, 2))
    matType = MAT_TYPE_SKINNY;
else
    matType = MAT_TYPE_FAT;
end

paramRho        = sSolverParams.paramRho;
numIterations   = sSolverParams.numIterations;

mX = zeros(size(mA, 2), numIterations);

% Caching facotirzation
mC = MatrixFactorize(mA, paramRho, matType);
vAy = mA.' * vY;

vZ = ProxL1(vX, paramLambda / paramRho);
vU = vX - vZ;

mX(:, 1) = vX;

for ii = 2:numIterations

    vQ = vAy + (paramRho * (vZ - vU));

    % Matrix Inversion Lemma
    switch(matType)
        case(MAT_TYPE_SKINNY)
            vX = mC \ vQ;
        case(MAT_TYPE_FAT)
            vX = (vQ / paramRho) - ((mA.' * (mC \ (mA * vQ))) / (paramRho * paramRho));
    end

    vZ = ProxL1(vX + vU, paramLambda / paramRho);
    vU = vU + vX - vZ;

    mX(:, ii) = vX;

end

end

function [ mC ] = MatrixFactorize( mA, paramRho, matType )

MAT_TYPE_SKINNY = 1;
MAT_TYPE_FAT    = 2;

switch(matType)
    case(MAT_TYPE_SKINNY)
        mI = speye(size(mA, 2));
        mC = decomposition((mA.' * mA) + (paramRho * mI), 'chol');
    case(MAT_TYPE_FAT)
        mI = speye(size(mA, 1));
        mC = decomposition(mI + ((1 / paramRho) * (mA * mA.')), 'chol');
end

end

function [ vX ] = ProxL1( vX, lambdaFactor )

% Soft Thresholding
vX = max(vX - lambdaFactor, 0) + min(vX + lambdaFactor, 0);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%VGG_KR_FROM_P Extract K, R from camera matrix.
%
%    [K,R,t] = VGG_KR_FROM_P(P [,noscale]) finds K, R, t such that P = K*R*[eye(3) -t].
%    It is det(R)==1.
%    K is scaled so that K(3,3)==1 and K(1,1)>0. Optional parameter noscale prevents this.
%
%    Works also generally for any P of size N-by-(N+1).
%    Works also for P of size N-by-N, then t is not computed.

% Author: Andrew Fitzgibbon <awf@robots.ox.ac.uk>
% Modified by werner.
% Date: 15 May 98

function [K, R, t] = vgg_KR_from_P(P, noscale)

N = size(P,1);
H = P(:,1:N);

[K,R] = vgg_rq(H);

if nargin < 2
  K = K / K(N,N);
  if K(1,1) < 0
    D = diag([-1 -1 ones(1,N-2)]);
    K = K * D;
    R = D * R;

  %  test = K*R;
  %  vgg_assert0(test/test(1,1) - H/H(1,1), 1e-07)
  end
end

if nargout > 2
  t = -P(:,1:N)\P(:,end);
end

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
incl = {'ucm', 'COM', 'utils', ...
	'semantics', 'segmentation', ...
	'ucmGT', 'classify', 'sceneClassification', 'allBenchmarks', ...
	'external', 'external/vlfeat/toolbox/', 'external/BSR/grouping/lib', 'external/BSR/bench/benchmarks'};
for i = 1:length(incl)
	addpath(incl{i});
end
vl_setup();
fprintf('Set up and ready to use...!\n');

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function [SRE,LRE,IV,RLV,RP,LIRE,HIRE,LISRE,HISRE,LILRE,HILRE] = getSizeParams(structNum,numLevels,planC)
%function getSizeParams(structNum)
%
%This function returns zone size based features for structure structNum.
%
%APA, 03/09/2015

if numel(structNum) == 1
    if ~exist('planC')
        global planC
    end
    indexS = planC{end};
    scanNum                             = getStructureAssociatedScan(structNum,planC);
    [rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
    [mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
    scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));
    SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
    [minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
    volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
    volToEval(volToEval==0)             = NaN;
else
    volToEval = structNum;
end

% Quantize using percentiles
%numLevels = 16;
qtlImage = zeros(size(volToEval));
minVal = min(volToEval(:));
maxVal = max(volToEval(:));
levelsV = linspace(minVal,maxVal,numLevels+1);
%qtlsV = 0:0.1:1;
qtlsV = 1:(numLevels+1);
qltCount = 1;
for qtl = 1:length(qtlsV)-1
    %lowQtl = quantile(volToEval(:),qtlsV(qltCount));
    %hiQtl = quantile(volToEval(:),qtlsV(qltCount+1));
    lowQtl = levelsV(qltCount);
    hiQtl = levelsV(qltCount+1);
    qtlImage(volToEval > lowQtl & volToEval <= hiQtl) = qltCount;
    qltCount = qltCount + 1;
end

% Compute Zone-Size matrix
clear global run_length_matrix_all image_property
global run_length_matrix_all image_property
image_property.num_voxels = sum(~isnan(volToEval(:)));
run_length_matrix_all = zeros(max(qtlImage(:)), max(size(qtlImage)));
for idx_intensity = 1:max(qtlImage(:))
    mat_isintensity = (qtlImage==idx_intensity);
    mat_connection = bwlabeln(mat_isintensity, 26); % allow the max connectivity
    for idx_group = 1:max(mat_connection(:))
        if size(run_length_matrix_all,2)< length(find(mat_connection== idx_group))
            run_length_matrix_all(idx_intensity, length(find(mat_connection== idx_group))) = 1;
        else
            run_length_matrix_all(idx_intensity, length(find(mat_connection== idx_group))) = run_length_matrix_all(idx_intensity, length(find(mat_connection== idx_group))) +1;
        end
    end
end

% Compute scalar features
SRE = run_length_SRE();
LRE = run_length_LRE();
IV = run_length_IV();
RLV = run_length_RLV();
RP = run_length_RP();
LIRE = run_length_LIRE();
HIRE = run_length_HIRE();
LISRE = run_length_LISRE();
HISRE = run_length_HISRE();
LILRE = run_length_LILRE();
HILRE = run_length_HILRE();

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [C_smooth]=patchSmoothFaceMeasure(varargin)

% function [C_smooth]=patchSmoothFaceMeasure(F,V,C,smoothPar)

%% Parse input
switch nargin
    case 3
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        smoothPar=[];
    case 4
        F=varargin{1};
        V=varargin{2};
        C=varargin{3};
        smoothPar=varargin{4};
end

smoothParDefault.lambda=0.5;
smoothParDefault.n=1;
smoothParDefault.faceFaceConnectivity=[];
smoothPar=structComplete(smoothPar,smoothParDefault,1);

%% Get connectivity array
if isempty(smoothPar.faceFaceConnectivity)
    [connectivityStruct]=patchConnectivity(F,V);
    faceFaceConnectivity=connectivityStruct.face.face;
end

%%

numSmoothIterations=smoothPar.n;
lambdaSmooth=smoothPar.lambda;

%%

nDims=size(C,2); %Number of dimensions
logicValid=faceFaceConnectivity>0;
C_smooth=C;
C_smooth_step=C;
for qIter=1:numSmoothIterations
    %Loop for all dimensions
    for qDim=1:1:nDims
        Xp=NaN(size(C,1),size(faceFaceConnectivity,2));
        Xp(logicValid)=C_smooth(faceFaceConnectivity(logicValid),qDim);
        Xp=gnanmean(Xp,2);
        C_smooth_step(:,qDim)=Xp;
    end
    C_smooth=((1-lambdaSmooth).*C_smooth)+(lambdaSmooth.*C_smooth_step);
end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode


text:
function below = isBelowPlane(point, varargin)
%ISBELOWPLANE Test whether a point is below or above a plane.
%
%   BELOW = isBelowPlane(POINT, PLANE)
%   where POINT is given as coordinate row vector [XP YP ZP], and PLANE is
%   given as a row containing initial point and 2 direction vectors,
%   return TRUE if POINT lie below PLANE.
%
%   Example
%   isBelowPlane([1 1 1], createPlane([1 2 3], [1 1 1]))
%   ans =
%       1
%   isBelowPlane([3 3 3], createPlane([1 2 3], [1 1 1]))
%   ans =
%       0
%
%   See also
%   planes3d, points3d, linePosition3d, planePosition

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2007-01-05
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

if length(varargin)==1
    plane = varargin{1};
elseif length(varargin)==2
    plane = createPlane(varargin{1}, varargin{2});
end

% ensure same dimension for parameters
if size(point, 1)==1
    point = repmat(point, [size(plane, 1) 1]);
end
if size(plane, 1)==1
    plane = repmat(plane, [size(point, 1) 1]);
end

% compute position of point projected on 3D line corresponding to plane
% normal, and returns true for points locatd below the plane (pos<=0).
below = linePosition3d(point, [plane(:, 1:3) planeNormal(plane)]) <= 0;

author: mattools


text:
function res = invert_nn_dw(net, ref, opts)

% setup param
[net,opts,x] = invert_nn_pre(net, ref, opts);
x0_size = cat(2,net.normalization.imageSize,opts.numRepeats);
x_momentum = zeros(x0_size, 'single') ;
load('x0_sigma.mat', 'x0_sigma');
y0 = repmat(ref, [1 1 1 opts.numRepeats]);
switch opts.task
    case {0,2}
        y0_sigma = norm(y0(:)) ;
    case {1}
        y0_sigma = 1;
end

%% --------------------------------------------------------------------
%%                                                 Perform optimisation
%% --------------------------------------------------------------------

% Run forward propogation once on the modified network before we
% begin backprop iterations - this is to exploit an optimization in
% vl_simplenn

% recored results
output = {} ;
prevlr = 0 ;

% Iterate until maxNumIterations to optimize the objective
% and generate the reconstuction
regu_out=[];
if ~isempty(opts.regu)
    regu_out=cell(1,numel(opts.regu));
end

res=[];
% visualize dropout mask
E=zeros(5,opts.maxNumIterations);
for t=1:opts.maxNumIterations

  % 1. Effectively does both forward and backward passes
  res = U_cnn_fb(net,res,x,opts.mask);

 % The current best feature we could generate
  switch net.cnn_mode
  case 'matconvnet'; y = res(end-1).x ;
  case 'caffe'; y = res.y ;
  end
  dr = zeros(size(x),'single'); % The derivative

  % 2. update the mask
  if ~isempty(opts.regu)
    for j=1:numel(opts.regu)
        regu_out{j} = res(opts.regu{j}{1}).x;
    end
  end

  % 3. calcuate gradient
  % 3.1 Data term
  E(1,t) = opts.lambdaD*res(end).x/(y0_sigma^2);

  % 3.2 Cost and derivative for TV\beta norm
  if opts.lambdaTV > 0
    tmp_x =x;
      if ~isempty(opts.im_w)
        tmp_x=reshape(reshape(x,[],3)*opts.im_w',size(x));
      end
    [r_,dr_] = tv(tmp_x,opts.TVbeta,opts.xmask) ;
    E(2,t) = opts.lambdaTV/2 * r_ ;
    dr = dr + opts.lambdaTV/2 * dr_ ;
  end
  % 3.2 Cost and derivative of L\beta norm
  if opts.lambdaL2 > 0
     tmp_x=x;
      if ~isempty(opts.im_w)
        tmp_x=reshape(reshape(x,[],3)*opts.im_w',size(x));
      end
    if ~isempty(opts.xmask)
        tmp_x = bsxfun(@times,x,opts.xmask);
    end
    if ~isempty(opts.regu_c)
      tmp_x = reshape(bsxfun(@minus,reshape(tmp_x,[],3),opts.regu_c),size(x));
    end
    r_ = sum(tmp_x(:).^opts.beta) ;
    dr_ = opts.beta * tmp_x.^(opts.beta-1) ;

    E(3,t) = opts.lambdaL2/2 * r_ ;
    dr = dr + opts.lambdaL2/2 * dr_ ;
  end

  E(2:3,t) = E(2:3,t) / (x0_sigma^2) ;
  E(4,t) = sum(E(1:3,t)) ;
  lr = opts.learningRate(min(t, numel(opts.learningRate))) ;

  % when the learning rate changes suddently, it is not
  % possible for the gradient to crrect the momentum properly
  % causing the algorithm to overshoot for several iterations
  if lr ~= prevlr
    fprintf('switching learning rate (%f to %f) and resetting momentum\n', ...
      prevlr, lr) ;
    x_momentum = 0 * x_momentum ;
    prevlr = lr ;
  end

  % x_momentum combines the current gradient and the previous gradients
  % with decay (opts.momentum)
  dr_d = opts.lambdaD*res(1).dzdx;
  %keyboard
  % [max(dr_d(:)) max(dr(:)/(x0_sigma^2)) max(opts.lambdaD*res(1).dzdx(:))]
  dsp_p=x0_sigma^2/y0_sigma^2;
  x_momentum = opts.momentum * x_momentum ...
      - lr * dr ...
      - (lr *  (x0_sigma^2/y0_sigma^2)*dr_d);

    % gradient clipping
    tmp_max =max(abs(x_momentum(:)));
    if opts.grad_ran(1)>=0
        if tmp_max>opts.grad_ran(2)
            x_momentum=x_momentum/tmp_max*opts.grad_ran(2);
        elseif tmp_max<opts.grad_ran(1)
            x_momentum=x_momentum/tmp_max*opts.grad_ran(1);
        end
    end
    tmp_max =max(abs(x_momentum(:)));

    if opts.dsp>=0
      fprintf('iter:%05d: max_grad=%3.4g, err_data=%5.4g, err_TV=%5.4g, err_L2=%5.4g, err_all=%5.4g;\n', t, tmp_max,E(1,t),E(2,t),E(3,t), E(4,t)) ;
  end

  % This is the main update step (we are updating the the variable
  % along the gradient
    if ~isempty(opts.xmask)
        x_momentum = bsxfun(@times,x_momentum,opts.xmask);
    end
    x = x + x_momentum ;

  % truncate output
  if opts.do_thres
      tmp_m = opts.denormalize(x);
      tmp_m(tmp_m<0)=0;
      tmp_m(tmp_m>255)=255;
      x = opts.normalize(tmp_m);
  end

  %% -----------------------------------------------------------------------
  %% Plots - Generate several plots to keep track of our progress
  %% -----------------------------------------------------------------------

  if opts.dsp>0 && mod(t-1,opts.dsp)==0
    figure(1) ; clf ;
    output{end+1} = opts.denormalize(x) ;
    subplot(3,2,[1 3]) ;
    if opts.numRepeats > 1
      vl_imarraysc(output{end}) ;
    else
      imagesc(uint8(output{end})) ;
    end
    axis image ;

    subplot(3,2,6) ;
    if ~isempty(opts.xmask)
        hist(reshape(x(opts.xmask==1),[],1),100);
    else
        hist(x(:),100) ;
    end
    grid on ;
    title('histogram of x') ;

    subplot(3,2,5) ;
    plot(E') ;
    h = legend('recon', 'tv_reg', 'l2_reg', 'tot') ;
    set(h,'color','none') ; grid on ;
    title(sprintf('iter:%d \\lambda_{tv}:%g \\lambda_{l2}:%g rate:%g obj:%s', ...
                  t, opts.lambdaTV, opts.lambdaL2, lr, opts.objective)) ;

   drawnow ;

  end % end if(mod(t-1,25) == 0)
end % end loop over maxNumIterations
if opts.dsp<=0
    output{end+1} = opts.denormalize(x) ;
end

% Compute the features optained using feedforward on the computed inverse
res_nn = res;

clear res;
%res.input = NaN;
res.output = output ;
res.energy = E ;
%res.y0 = y0 ;
%res.y = res_nn(end-1).x ;
%res.opts = opts ;
%res.err = res_nn(end).x / y0_sigma^2 ;

% --------------------------------------------------------------------
function [e, dx] = tv(x,beta,mask)
% dw: increase only boundary weight doesn't help
% which equivalently push the boundary inside by 1-pix
% --------------------------------------------------------------------
if(~exist('beta', 'var'))
  beta = 1; % the power to which the TV norm is raized
end
if(~exist('mask', 'var'));mask=[];end

d1 = x(:,[2:end end],:,:) - x ;
d2 = x([2:end end],:,:,:) - x ;
if ~isempty(mask)
    new_m = imdilate(mask,strel('disk',1))>0;
    d1 = bsxfun(@times,d1, new_m);
    d2 = bsxfun(@times,d2, new_m);
end
v = sqrt(d1.*d1 + d2.*d2).^beta ;
e = sum(sum(sum(sum(v)))) ;
if nargout > 1
  d1_ = (max(v, 1e-5).^(2*(beta/2-1)/beta)) .* d1;
  d2_ = (max(v, 1e-5).^(2*(beta/2-1)/beta)) .* d2;
  % take grad from both end
  %  x_{i,j}-x_{i,j+1}
  %  x_{i,j-1}-x_{i,j}
  d11 = d1_(:,[1 1:end-1],:,:) - d1_ ;
  d22 = d2_([1 1:end-1],:,:,:) - d2_ ;
  d11(:,1,:,:) = - d1_(:,1,:,:) ;
  d22(1,:,:,:) = - d2_(1,:,:,:) ;
  dx = beta*(d11 + d22);
if ~isempty(mask)
    dx = bsxfun(@times,dx,mask);
end
  if(any(isnan(dx)))
  end
end

% --------------------------------------------------------------------
function test_tv()
% --------------------------------------------------------------------
x = randn(5,6,1,1) ;
[e,dr] = tv(x,6) ;
vl_testder(@(x) tv(x,6), x, 1, dr, 1e-3) ;

author: donglaiw

repo: mNeuron

save_path: github-repos/MATLAB/donglaiw-mNeuron