text:
function [hPatch,api] = imRectRot( varargin )
% Create a draggable, resizable, rotatable rectangle or ellipse.
%
% The 'ellipse' param determines if the displayed object is a rectangle or
% ellipse. The object is identical in both cases, only the display changes.
% The created object may be queried or controlled programatically by using
% the returned api.
%
% The 'rotate' param determines overall behavior of the created object. If
% rotate=0, the resulting object (rect or ellipse) is axis aligned. In
% terms of the graphical interface it is identical to Matlab's imrect (can
% drag by clicking in interior or resize by clicking on edges), although it
% is much less cpu intensive. If rotate>0, the resulting object is
% rotatable. In addition to the interface for a non-rotatable object, four
% control points are present, one at the center of each edge. Three of
% these have color given by the 'color' flag, the last has color 'colorc'.
% The odd colored control point is used to display orientation. Dragging
% this control point or the one opposite changes orientation, dragging the
% remaining two resizes the object symmetrically. Finally, when creating a
% rotatable object, the first drag determines the major axes (height) of
% the object, with the width set to height*rotate (hence the rotate param
% also determines aspect ratio of newly created objects). Using this
% control scheme, an object can be naturally specified with two drags: the
% first is used to draw the major axes, the second to adjust the width.
%
% Position is represented by [x y w h theta], where [x,y] give the top/left
% corner of the rect PRIOR to rotation, [w,h] are the width and height, and
% theta is the angle in degrees. The final rect is given by first placing
% the rect at [x,y], then rotating it by theta around it's center. The
% advantage of this is that if theta=0, the first four elements are
% identical to the standard rect representation. The disadvantage is that
% [x,y] need not lie in the interior of the rect after rotation.
%
% USAGE
%  [h,api] = imRectRot( varargin )
%
% INPUTS
%  varargin   - parameters (struct or name/value pairs)
%   .hParent    - [gca] object parent, typically an axes object
%   .ellipse    - [0] if true display ellipse otherwise display rectangle
%   .rotate     - [1] determines if object is axis aligned
%   .pos        - [] initial pos vector [x y w h theta] or [] or [x y]
%   .lims       - [] rectangle defining valid region for object placement
%   .showLims   - [0] draw rectangle representing lims
%   .color      - ['g'] color for the displayed object
%   .colorc     - ['b'] color for the control point displaying orientation
%   .lw         - [2] 'LineWidth' property for the displayed object
%   .ls         - ['-'] 'LineStyle' property for the displayed object
%   .cross      - [0] if 1 show diagonal, if 2 show cross
%
% OUTPUTS
%  h          - handle used to delete object
%  api        - interface allowing access to created object
%  .getPos()       - get position - returns 5 elt pos
%  .setPos(pos)    - set position (while respecting constraints)
%  .setPosLock(b)  - if lock set (b==true), object cannot change
%  .setSizLock(b)  - if lock set (b==true), object cannot change size
%  .setDrgLock(b)  - if lock set (b==true), object cannot be dragged
%  .setSidLock(lk) - [4x1] set locks for each side (tp/rt/bt/lf)
%  .setPosChnCb(f) - whenever pos changes (even slightly), calls f(pos)
%  .setPosSetCb(f) - whenever pos finished changing, calls f(pos)
%  .uistack(...)   - calls 'uistack( [objectHandles], ... )', see uistack
%  .setStyle(...)  - set line style (ls), width (lw), color and colorc
%
% EXAMPLE - interactively place simple axis aligned rectangle
%  figure(1), imshow peppers.png;
%  [h,api]=imRectRot('rotate',0);
%  api.setPosChnCb( @(pos) disp(num2str(pos)) );
%
% EXAMPLE - create rotatable ellpise that falls inside image
%  figure(1); I=imread('cameraman.tif'); imshow(I); siz=size(I);
%  [h,api]=imRectRot('pos',[60 60 40 40 45],'lims',[1 1 siz(1:2)-2 0],...
%    'showLims',1,'ellipse',1,'rotate',1,'color','w','colorc','y'  );
%  api.setPosSetCb( @(pos) disp(num2str(pos)) );
%
% See also IMRECT, RECTANGLE, PATCH
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.51
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% global variables (shared by all functions below)
[hParent,pos,lims,rotate,crXs,crYs,posChnCb,posSetCb,posLock,sizLock,...
  dgrLock,sidLock,hAx,hFig,hPatch,hBnds,hCntr,hEll] = deal([]);

% intitialize
intitialize( varargin{:} );

% create api
api = struct('getPos',@getPos, 'setPos',@setPos, 'uistack',@uistack1, ...
  'setPosChnCb',@setPosChnCb, 'setPosSetCb',@setPosSetCb, ...
  'setPosLock',@setPosLock, 'setSizLock',@setSizLock, ...
  'setDrgLock',@setDrgLock, 'setSidLock',@setSidLock, ...
  'setStyle',@setStyle );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function intitialize( varargin )
    % get default arguments
    dfs={'hParent',gca, 'ellipse',0, 'rotate',1, 'pos',[], 'lims',[], ...
      'showLims',0,'color','g','colorc','b','lw',2,'ls','-','cross',0};
    [hParent,ellipse,rotate,pos,lims,showLims,color,colorc,lw,ls,cross]=...
      getPrmDflt(varargin,dfs,1);
    if(numel(pos)==4), pos=[pos(:)' 0]; end
    if(numel(lims)==4), lims=[lims(:)' 0]; end

    % get figure and axes handles
    hAx = ancestor(hParent,'axes'); assert(~isempty(hAx));
    hFig = ancestor(hAx,'figure'); assert(~isempty(hFig));
    set( hFig, 'CurrentAxes', hAx );

    % optionally display limits
    if(showLims && ~isempty(lims)), [~,xs,ys]=rectToCorners(lims);
      for j=1:4, ids=mod([j-1 j],4)+1; line(xs(ids),ys(ids)); end
    end

    % create objects for rectangle display/interface
    if( ellipse ), linePrp={'color',[.7 .7 .7],'LineWidth',1};
    else linePrp={'color',color,'LineWidth',lw,'LineStyle',ls}; end
    if(isempty(pos)); vis='off'; else vis='on'; end;
    assert(cross==0 || cross==1 || cross==2);
    for i=1:4+cross, hBnds(i)=line(linePrp{:},'Visible',vis); end

    % create transparent patch to capture clicks in object interior
    hPatch=patch('FaceColor','none','EdgeColor','none');

    % create objects for ellipse display/interface
    if( ellipse )
      ellPrp={'color',color,'LineWidth',lw,'LineStyle',ls};
      ts=linspace(-180,180,50); crXs=cosd(ts); crYs=sind(ts);
      hold on; hEll=plot(crXs,crYs,ellPrp{:},'Visible',vis); hold off;
    end

    % create objects for rotation display/interface
    if( rotate )
      circPrp={'Curvature',[1 1],'FaceColor',color,'EdgeColor',color};
      for i=1:4, hCntr(i)=rectangle(circPrp{:},'Visible',vis); end
      set(hCntr(1),'LineWidth',lw,'FaceColor',colorc);
    end

    % set callbacks on all objects
    hs=[hPatch hBnds hCntr hEll];
    set(hs,'ButtonDownFcn',{@btnDwn,0},'DeleteFcn',@deleteFcn);

    % set or query initial position
    posLock=0; sizLock=0; dgrLock=0; sidLock=[0 0 0 0];
    if(length(pos)==5), setPos(pos,1); else
      btnDwn([],[],-1); waitfor(hFig,'WindowButtonUpFcn','');
    end
  end

  function [pc,rs,R] = rectInfo( pos0 )
    % return rectangle center, radii, and rotation matrix
    t=pos0(5); c=cosd(t); s=sind(t); R=[c -s; s c];
    rs=pos0(3:4)/2; pc=pos0(1:2)+rs;
  end

  function [pts,xs,ys,pc] = rectToCorners( pos0 )
    % return 4 rect corners in real world coords
    [pc,rs,R]=rectInfo( pos0 );
    x0=-rs(1); x1=rs(1); y0=-rs(2); y1=rs(2);
    pts=[x0 y0; x1 y0; x1 y1; x0 y1]*R'+pc(ones(4,1),:);
    xs=pts(:,1); ys=pts(:,2);
  end

  function pos1 = cornersToRect( pos0, x0, y0, x1, y1 )
    % compute pos from 4 corners given in rect coords
    [pc,~,R] = rectInfo( pos0 );
    p0=[x0 y0]*R'+pc; p1=[x1 y1]*R'+pc;
    pc=(p1+p0)/2; p0=(p0-pc)*R; p1=(p1-pc)*R;
    pos1 = [pc+p0 p1-p0 pos0(5)];
  end

  function setPos( posNew, ignoreLims )
    if(isempty(hBnds) || isempty(hPatch)); return; end; L=lims;
    % if corners fall outside lims don't use posNew
    if( ~isempty(L) && (nargin<2 || ~ignoreLims) )
      if( posNew(5)==lims(5) )
        % lims and posNew at same orient (do everything in rect coords)
        [pc,rs,R]=rectInfo(posNew); p0=-rs; p1=rs;
        pts=(rectToCorners(L)-pc(ones(4,1),:))*R;
        L0=[pts(1,1) pts(1,2)]; L1=[pts(2,1) pts(3,2)];
        p0=min(max(p0,L0),L1); p1=max(min(p1,L1),p0);
        posNew=cornersToRect(posNew,p0(1),p0(2),p1(1),p1(2));
      else
        % lims and posNew at diff orient (do everything in lims coords)
        [pc,rs,R]=rectInfo(L); p0=-rs; p1=rs;
        p=(rectToCorners(posNew)-pc(ones(4,1),:))*R; xs=p(:,1); ys=p(:,2);
        valid=[xs>=p0(1) xs<=p1(1) ys>=p0(2) ys<=p1(2)];
        if(~all(valid(:))), return; end
      end
    end
    % create invisible patch to captures key presses
    pos=posNew; [~,xs,ys]=rectToCorners(pos);
    vert=[xs ys ones(4,1)]; face=1:4;
    set(hPatch,'Faces',face,'Vertices',vert);
    % draw ellipse
    if(~isempty(hEll)), [pc,rs]=rectInfo(posNew); th=pos(5);
      xsEll = rs(1)*crXs*cosd(-th)+rs(2)*crYs*sind(-th)+pc(1);
      ysEll = rs(2)*crYs*cosd(-th)-rs(1)*crXs*sind(-th)+pc(2);
      set(hEll,'XData',xsEll,'YData',ysEll);
    end
    % draw rectangle boundaries and control circles
    r=axisUnitsPerCentimeter();
    r=max( r*.1, min(mean(pos(3:4))/10,r) );
    for i=1:length(hBnds), ids=mod([i-1 i],4)+1;
      if(i==5), ids=[1 3]; elseif(i==6), ids=[2 4]; end
      set(hBnds(i),'Xdata',xs(ids),'Ydata',ys(ids));
      if(rotate && i<=4), x=mean(xs(ids)); y=mean(ys(ids));
        set(hCntr(i),'Position',[x-r y-r 2*r 2*r]);
      end
    end
  end

  function pnt = getCurPnt()
    pnt=get(hAx,'CurrentPoint'); pnt=pnt([1 3]);
  end

  function d = axisUnitsPerCentimeter()
    % approximage absolute axes size
    units=get(hAx,'units'); set(hAx,'units','centimeters');
    cm=get(hAx,'Position'); cm=max(cm(3:4));
    set(hAx,'units',units);
    % axes size in axes units
    xLim=get(gca,'xLim'); yLim=get(gca,'yLim');
    sizInner = max([xLim(2)-xLim(1),yLim(2)-yLim(1)]);
    % units per centimeter
    d = sizInner/cm;
  end

  function [side,cursor,flag] = getSide( pnt0 )
    % get pnt in coordintates of center of rectangle
    [pc,rs,R]=rectInfo(pos); pnt=(pnt0-pc)*R; th=pos(5);
    % side(i) -1=lf/tp bnd; 0=interior; +1:rt/bt bnd; 2=center
    t = axisUnitsPerCentimeter(); side=[0 0]; flag=0;
    t = max( t*.1, min(mean(pos(3:4))/5,t) );
    for i=1:2
      ti = min(t,rs(i)/3);
      if( abs(pnt(i))<ti && rotate ), side(i)=2; % near center
      elseif( pnt(i)<-rs(i)+ti ), side(i)=-1; % near lf/tp boundary
      elseif( pnt(i)>rs(i)-ti), side(i)=1; % near rt/bt boundary
      end
    end
    % flag: -1: none, 0=resize; 1=drag; 2=rotate; 3=symmetric-resize
    if((sidLock(4)&&side(1)==-1)||(sidLock(2)&&side(1)==1)), side(1)=0; end
    if((sidLock(1)&&side(2)==-1)||(sidLock(3)&&side(2)==1)), side(2)=0; end
    if(any(side==0) || all(side==2)), side(side==2)=0; end
    if(side(1)==2), flag=2; cursor='crosshair'; return; end
    if(sizLock), side=[0 0]; end
    if(all(side==0) && dgrLock), cursor=''; flag=-1; return; end
    if(all(side==0)), flag=1; side=pnt0; cursor='fleur'; return; end
    if(side(2)==2), flag=3; side(2)=0; end
    % select cursor based on position
    cs={'bottom','botr','right','topr','top','topl','left','botl'};
    a=mod(round((atan2(side(1),side(2))/pi-th/180)*4),8)+1; cursor=cs{a};
  end

  function btnDwn( h, evnt, flag ) %#ok<INUSL>
    if(isempty(hBnds) || isempty(hPatch) || posLock), return; end
    if( flag==-1 ) % create new rectangle
      if(isempty(pos)), anchor=ginput(1); else anchor=pos(1:2); end
      pos=[anchor 1 1 0]; setPos(pos);
      set([hBnds hCntr hEll],'Visible','on');
      if(rotate), cursor='crosshair'; flag=4;
      else cursor='botr'; flag=0; end
    else % resize, rotate or drag rectangle
      pnt=getCurPnt(); [anchor,cursor,flag]=getSide( pnt );
      if(flag==-1), return; end
    end
    set( hFig, 'Pointer', cursor );
    set( hFig, 'WindowButtonMotionFcn',{@drag,flag,anchor,pos} );
    set( hFig, 'WindowButtonUpFcn', @stopDrag );
  end

  function drag( h, evnt, flag, anchor, pos0 ) %#ok<INUSL>
    if(isempty(hBnds) || isempty(hPatch)); return; end;
    pnt = getCurPnt();
    if( flag==1 ) % shift rectangle by del=pnt-anchor
      setPos( [pos0(1:2)+(pnt-anchor) pos0(3:5)] );
    elseif( flag==0 || flag==3 ) % resize in rectangle coordinate frame
      [pc,rs,R]=rectInfo(pos0); p0=-rs; p1=rs; pnt=(pnt-pc)*R;
      if( flag==3 ), p0(1)=-pnt(1); p1(1)=pnt(1); else
        for i=1:2, if(anchor(i)>0), p1(i)=pnt(i); end; end
        for i=1:2, if(anchor(i)<0), p0(i)=pnt(i); end; end
      end; p0a=min(p0,p1); p1=max(p0,p1); p0=p0a;
      setPos(cornersToRect(pos0,p0(1),p0(2),p1(1),p1(2)));
    elseif( flag==2 || flag==4 ) % rotate rectangle
      [~,xs,ys]=rectToCorners(pos0);
      if(anchor(2)==-1), ids=[3 4]; else ids=[1 2]; end
      p0=mean([xs(ids) ys(ids)]); pc=(p0+pnt)/2;
      if(anchor(2)==-1), d=pnt-p0; else d=p0-pnt; end
      h=sqrt(sum(d.^2)); th=atan2(d(1),-d(2))/pi*180;
      if(flag==2), w=pos(3); else w=h*rotate; end
      setPos([pc(1)-w/2 pc(2)-h/2 w h th]);
    end
    if(~isempty(posChnCb)); posChnCb(pos); end;
    drawnow; % update display
  end

  function stopDrag( h, evnt ) %#ok<INUSD>
    set( hFig, 'WindowButtonMotionFcn','');
    set( hFig, 'WindowButtonUpFcn','');
    set( hFig, 'Pointer', 'arrow' );
    if(~isempty(posSetCb)); posSetCb(pos); end;
  end

  function deleteFcn( h, evnt ) %#ok<INUSD>
    [posChnCb,posSetCb]=deal([]); hs=[hPatch hBnds hCntr hEll];
    for i=1:length(hs), if(ishandle(hs(i))); delete(hs(i)); end; end
    [hBnds,hPatch,hCntr,hEll]=deal([]);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function pos0 = getPos(), pos0=pos; end

  function setPosChnCb( f ), posChnCb=f; end

  function setPosSetCb( f ), posSetCb=f; end

  function setPosLock( b ), posLock=b; end

  function setSizLock( b ), sizLock=b; end

  function setDrgLock( b ), dgrLock=b; end

  function setSidLock( lk )
    sidLock=lk; vis={'off','on'};
    for i=1:4, set(hCntr(i),'Visible',vis{2-lk(i)}); end
  end

  function uistack1( varargin )
    if(isempty(hBnds) || isempty(hPatch)); return; end;
    uistack( [hBnds hPatch hCntr hEll], varargin{:} );
  end

  function setStyle( ls, lw, color, colorc )
    if(isempty(hEll)), h=hBnds; else h=hEll; end
    set(h,'LineStyle',ls,'LineWidth',lw,'color',color);
    if(~rotate), return; end; if(nargin<4), colorc='b'; end
    set(hCntr,'FaceColor',color,'EdgeColor',color);
    set(hCntr(1),'LineWidth',lw,'FaceColor',colorc);
  end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function goldstein_price_test ( )

%*****************************************************************************80
%
%% GOLDSTEIN_PRICE_TEST works with the Goldstein-Price function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GOLDSTEIN_PRICE_TEST:\n' );
  fprintf ( 1, '  Test COORDINATE_SEARCH with the Goldstein-Price function.\n' );
  n = 2;

  x = [ -0.5, 0.25 ];
  r8vec_print ( n, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', goldstein_price ( x ) );

  flag = 0;
  x = coordinate_search ( x, @goldstein_price, flag );
  r8vec_print ( n, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g\n', goldstein_price ( x ) );
%
%  Repeat with more difficult start.
%
  x = [ -4.0, 5.0 ];
  r8vec_print ( n, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', goldstein_price ( x ) );

  flag = 0;
  x = coordinate_search ( x, @goldstein_price, flag );
  r8vec_print ( n, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g\n', goldstein_price ( x ) );
%
%  Demonstrate correct minimizer.
%
  x = [ 0.0, -1.0 ];
  r8vec_print ( n, x, '  Correct minimizer X*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', goldstein_price ( x ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = gtrack(newVarName,titleFmt)
% default format for printing coordinates in title
if nargin<2,	titleFmt = '%3.5f';		end

% === mouse states ===
global mouseButton; % 'down' or 'up'

% track mouseMovement
global mouseMotion; % 'moving' or 'stationary'
global mouseMoveDirection; % 'N', 'S', 'E', 'W'; 'NW', 'NE', 'SW', 'SE', 'stationary'
global whichMouseButton; % 'left' or 'right'
global currentLocation; % x, y coordinates
global previousLocation; % x, y coordinates

previousLocation = 'unassigned';

% get current figure event functions
currFcn = get(gcf, 'windowbuttonmotionfcn');
currFcn2 = get(gcf, 'windowbuttondownfcn');
currTitle = get(get(gca, 'Title'), 'String');

% add data to figure handles
handles = guidata(gca);
if (isfield(handles,'ID') & handles.ID==1)
	disp('gtrack is already active.');
	return;
else
	handles.ID = 1;
end
handles.currFcn = currFcn;
handles.currFcn2 = currFcn2;
handles.currTitle = currTitle;
handles.theState = uisuspend(gcf);
guidata(gca, handles);

% set event functions
set(gcf,'Pointer','crosshair');
set(gcf, 'windowbuttonmotionfcn', @gtrack_OnMouseMove);
set(gcf, 'windowbuttondownfcn', @gtrack_OnMouseDown);
set(gcf, 'windowbuttonupfcn', @gtrack_OnMouseUp);

% declare variables
xInd = 0;
yInd = 0;
clickData = [];

% set output mode
if nargout,
	uiMode = 'uiwait';		% use UIWAIT and return to clickData
	uiwait;
elseif nargin && isvarname(newVarName)	% dont'e use UIWAIT and assign in caller
	uiMode = 'nowait';				% workspace to variable newVarName
else
	uiMode = 'noreturn';	% dont't use UIWAIT and don't return results (print only)
end

%% --- nested functions ---------------------------------------------------

%% mouse move callback
function gtrack_OnMouseMove(src,evnt)

% get mouse position
pt = get(gca, 'CurrentPoint');
xInd = pt(1, 1);
yInd = pt(1, 2);
currentLocation = [xInd yInd];

if (~strcmp(previousLocation, 'unassigned'))
    changesInVectoralDistance = currentLocation - previousLocation;

    xDelta = changesInVectoralDistance(1, 1);
    yDelta = changesInVectoralDistance(1, 2);

    % degreesInMovement = NaN, if there is no movement
    degreesInMovement = atan(yDelta/xDelta)*(180/pi);

    if (xDelta == 0) && (yDelta == 0)
        mouseMoveDirection = 'stationary';
        mouseMotion = 'stationary';
    elseif (xDelta == 0) && (yDelta > 0)
        mouseMoveDirection = 'N';
        mouseMotion = 'moving';
    elseif (xDelta == 0) && (yDelta < 0)
        mouseMoveDirection = 'S';
        mouseMotion = 'moving';
    elseif (xDelta > 0) && (yDelta == 0)
        mouseMoveDirection = 'E';
        mouseMotion = 'moving';
    elseif (xDelta < 0) && (yDelta == 0)
        mouseMoveDirection = 'W';
        mouseMotion = 'moving';
        degreesInMovement = degreesInMovement + 180;
    elseif (xDelta > 0) && (yDelta > 0)
        mouseMoveDirection = 'NE';
        mouseMotion = 'moving';
    elseif (xDelta < 0) && (yDelta > 0)
        mouseMoveDirection = 'NW';
        mouseMotion = 'moving';
        degreesInMovement = degreesInMovement + 180;
    elseif (xDelta < 0) && (yDelta < 0)
        mouseMoveDirection = 'SW';
        mouseMotion = 'moving';
        degreesInMovement = degreesInMovement + 180;
    elseif (xDelta > 0) && (yDelta < 0)
        mouseMoveDirection = 'SE';
        mouseMotion = 'moving';
    end

    distanceMoved = sqrt((xDelta^2)+(yDelta^2));
    polarCoordinates = [distanceMoved degreesInMovement];

    fprintf('Mouse moved in the direction of %s at %s degrees. Distance: %s\n', ...
        mouseMoveDirection, num2str(degreesInMovement), num2str(distanceMoved));
end

% check if its within axes limits
xLim = get(gca, 'XLim');
yLim = get(gca, 'YLim');
if xInd < xLim(1) | xInd > xLim(2)
	title('Out of X limit');
	return;
end
if yInd < yLim(1) | yInd > yLim(2)
	title('Out of Y limit');
	return;
end

% update figure title
try
	title(['X = ' num2str(xInd,titleFmt) ', Y = ' num2str(yInd,titleFmt)]);
% possibility of wrong format strings...
catch
	gtrack_Off()
	error('TRACK: Error printing coordinates. Check that you used a valid format string.')
end

% update the location
previousLocation = currentLocation;

end

%% mouse click callback
function gtrack_OnMouseDown(src,evnt)
whichMouseButton = get(gcf,'SelectionType');
if strcmp(whichMouseButton,'alt')
	% gtrack_Off % if left button, terminate
    %return
    whichMouseButton = 'right';
end

if strcmp(whichMouseButton,'normal')
	% gtrack_Off
	%return
    whichMouseButton = 'left';
end

% else add click to clickData
clickData(end+1).x = xInd;
clickData(end).y = yInd;
fprintf('%s Mouse Button Down at X = %f   Y = %f\n',whichMouseButton, xInd, yInd); % update mouse position

if (strcmp(mouseMotion, 'moving'))
    fprintf('Mouse is dragged.\n\n');
end

end

function gtrack_OnMouseUp(src,evnt)
whichMouseButton = get(gcf,'SelectionType');
if strcmp(whichMouseButton,'alt')
	% gtrack_Off % if left button, terminate
	%return
    whichMouseButton = 'right';
end

if strcmp(whichMouseButton,'normal')
	% gtrack_Off
	%return
    whichMouseButton = 'left';
end

% else add click to clickData
clickData(end+1).x = xInd;
clickData(end).y = yInd;
fprintf('%s Mouse Button Up at X = %f   Y = %f\n',whichMouseButton, xInd, yInd); % update mouse position

end

%% terminate callback
function gtrack_Off(src,evnt)

% restore default figure properties
handles = guidata(gca);
set(gcf, 'windowbuttonmotionfcn', handles.currFcn);
set(gcf, 'windowbuttondownfcn', handles.currFcn2);
set(gcf,'Pointer','arrow');
title(handles.currTitle);
uirestore(handles.theState);
handles.ID=0;
guidata(gca,handles);

% if there are outputs to assign do so
switch uiMode
	case 'uiwait'	% data return as output argument (clickData)
		varargout{1} = clickData;
		uiresume,
	case 'nowait'	% data assigned in base workspace as new variable
		assignin('base',newVarName,clickData);
		fprintf('Variable %s assigned with click data.\n',newVarName);
	case 'noreturn'
					% nothing to return
end

end

%% --- end nested functions -----------------------------------------------

end % end everything

% improved version of GTRACK
%
% GTRACK Track mouse position and show coordinates in figure title.
%
% 	GTRACK Activates GTRACK. Once it is active the mouse position is
% 	constantly tracked and printed on the figure title. A left-click will
% 	print the coordinates in the command line and store them. Clicking the
% 	mouse right button deactivates GTRACK.
%
% USAGE
% 	gtrack() tracks the mouse and prints coordinates in the command line.
%
% 	clickData = gtrack() will return the click positions in clickData using
% 	UIWAIT. Matlab will be in wait mode until the user finishes clicking.
%
% 	gtrack('newVar') tracks the mouse and creates a new variable in the
% 	base workspace called 'newVar' with the click coordinates. This mode
% 	does not use UIWAIT.
%
%	gtrack([],titleFormat) uses titleFormat as the format string for
%	printing the mouse coordinates in the title.
%

author: Sable

repo: mcbench-benchmarks


text:
function same=strcomp(str1,str2)
%A local strcmp routine.   Used.  Not sure why. Left over from 1996.
%
%  same=strcomp(str1,str2);
%
%Returns 1 if str1 == str2
%Returns 0 otherwise.
%
%4/12/96 gmb	Why did I have to write this?
%

if length(str1)~=length(str2)
	same=0;
	return
end

for i=1:length(str1)
	if (str1(i)~=str2(i))
		same=0;
		return
	end
end

same=1;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% std_spec() - Returns the ICA component spectra for a dataset. Updates the EEG structure
%              in the Matlab environment and in the .set file as well. Saves the spectra
%              in a file.
% Usage:
%           >> [spec freqs] = std_spec(EEG, 'key', 'val', ...);
%
%              Computes the mean spectra of the activites of specified components of the
%              supplied dataset. The spectra are saved in a Matlab file. If such a file
%              already exists, loads the spectral information from this file.
%              Options (below) specify which components to use, and the desired frequency
%              range. There is also an option to specify other spectopo() input variables
%              (see >> help spectopo for details).
%
%              Returns the removed mean spectra of the selected ICA components in the
%              requested frequency range. If the spectra were computed previously but a
%              different frequency range is selected, there is an overwrite option.
%              so. The function will load previously computed log spectra, if any, and
%              will remove the mean from the requested frequency range. The frequencies
%              vector is also returned.
% Inputs:
%   EEG - a loaded epoched EEG dataset structure.
%
% Optional inputs:
%   'components' - [numeric vector] components of the EEG structure for which
%                  activation ERPs will be computed. Note that because
%                  computation of component spectra is relatively fast, all
%                  components spectra are computed and saved. Only selected
%                  component are returned by the function to Matlab
%                  {default|[] -> all}
%   'channels'   - [cell array] channels of the EEG structure for which
%                  activation spectrum will be computed. Note that because
%                  computation of spectrum is relatively fast, all channels
%                  spectrum are computed and saved. Only selected channels
%                  are returned by the function to Matlab
%                  {default|[] -> none}
%   'specmode'   - ['psd'|'fft'] method to compute spectral
%                  decomposition. 'psd' uses the spectopo function. 'fft'
%                  uses a simple fft on each trial.
%   'epochlim'   - [min max] for FFT on continuous data, extract data
%                  epochs with specific epoch limits in seconds (see also
%                  'epochrecur' below). Default is [0 1].
%   'epochrecur' - [float] for FFT on continuous data, set the automatic
%                  epoch extraction recurence interval (default is 1 second).
%   'timerange'  - [min max] use data within a specific time range before
%                  computing the data spectrum. For instance, for evoked
%                  data trials, it is recommended to use the baseline time
%                  period.
%   'freqrange'  - [minhz maxhz] frequency range (in Hz) within which to
%                  return the spectrum {default|[]: [0 sample rate/2]}.
%   'recompute'  - ['on'|'off'] force recomputing ERP file even if it is
%                  already on disk.
%
% Other optional spectral parameters:
%   All optional parameters to the spectopo function may be provided to this function
%   as well.
%
% Outputs:
%   spec      - the mean spectra (in dB) of the requested ICA components in the selected
%               frequency range (with the mean of each spectrum removed).
%   freqs     - a vector of frequencies at which the spectra have been computed.
%
% Files output or overwritten for ICA:
%               [dataset_filename].icaspec,   % raw spectrum of ICA components
%               [dataset_filename].icaspecm   % spectrum with the mean baseline removed
% Files output or overwritten for data:
%               [dataset_filename].datspec,
%               [dataset_filename].datspecm
%
% See also  spectopo(), std_erp(), std_ersp(), std_map(), std_preclust()
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, January, 2005

% Defunct:      0 -> if frequency range is different from saved spectra, ask via a
%                    pop-up window whether to keep existing spectra or to overwrite them.

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2004, arno@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

function [X, f, overwrt] = std_spec(EEG, varargin)

overwrt = 1; % deprecated
if nargin < 1
    help std_spec;
    return;
end;

% decode inputs
% -------------
if ~isempty(varargin)
    if ~isstr(varargin{1})
        varargin = { varargin{:} [] [] };
        if all(varargin{1} > 0)
            options = { 'components' varargin{1} 'freqrange' varargin{2} };
        else
            options = { 'channels' -varargin{1} 'freqrange' varargin{2} };
        end;
    else
        options = varargin;
    end;
else
    options = varargin;
end;

[g spec_opt] = finputcheck(options, { 'components' 'integer' []         [];
                                      'channels'   'cell'    {}         {};
                                      'timerange'  'float'   []         [];
                                      'specmode'   'string'  {'fft' 'psd' 'pmtm' 'pburg'} 'psd';
                                      'recompute'  'string'  { 'on' 'off' } 'off';
                                      'savetrials' 'string'  { 'on' 'off' } 'off';
                                      'epochlim'   'real'    []         [0 1];
                                      'epochrecur' 'real'    []         1;
                                      'rmcomps'    'cell'    []         cell(1,length(EEG));
                                      'nw'         'float'   []         4;
                                      'fileout'    'string'  []         '';
                                      'burgorder'  'integer' []         20;
                                      'interp'     'struct'  { }        struct([]);
                                      'nfft'       'integer' []         [];
                                      'freqrange'  'real'    []         [] }, 'std_spec', 'ignore');
if isstr(g), error(g); end;
if isfield(EEG,'icaweights')
   numc = size(EEG(1).icaweights,1);
else
   error('EEG.icaweights not found');
end
if isempty(g.components)
    g.components = 1:numc;
end

EEG_etc = [];

% filename
% --------
if isempty(g.fileout), g.fileout = fullfile(EEG(1).filepath, EEG(1).filename(1:end-4)); end;
if ~isempty(g.channels)
    filename = [ g.fileout '.datspec'];
    prefix = 'chan';
else
    filename = [ g.fileout '.icaspec'];
    prefix = 'comp';
end;

% SPEC information found in datasets
% ---------------------------------
if exist(filename) & strcmpi(g.recompute, 'off')

    fprintf('File "%s" found on disk, no need to recompute\n', filename);
    setinfo.filebase = g.fileout;
    if strcmpi(prefix, 'comp')
        [X tmp f] = std_readfile(setinfo, 'components', g.components, 'freqlimits', g.freqrange, 'measure', 'spec');
    else
        [X tmp f] = std_readfile(setinfo, 'channels', g.channels, 'freqlimits', g.freqrange, 'measure', 'spec');
    end;
    if ~isempty(X), return; end;
end

oritrials = EEG.trials;
if ~strcmpi(g.specmode, 'psd')
    if EEG(1).trials == 1,
        EEG = eeg_checkset(EEG, 'loaddata');
        EEG = eeg_regepochs(EEG, g.epochrecur, g.epochlim);
        g.trialindices = { [1:EEG(1).trials] };
        disp('Warning: continuous data, extracting 1-second epochs');
    end;
end;

% No SPEC information found
% ------------------------
options = {};
if ~isempty(g.rmcomps), options = { options{:} 'rmcomps' g.rmcomps }; end;
if ~isempty(g.interp),  options = { options{:} 'interp' g.interp }; end;
X          = [];
boundaries = [];
for dat = 1:length(EEG)
    if strcmpi(prefix, 'comp')
        tmpdata = eeg_getdatact(EEG(dat), 'component', [1:size(EEG(dat).icaweights,1)], 'trialindices', g.trialindices{dat} );
    else
        EEG(dat).data = eeg_getdatact(EEG(dat), 'channel', [1:EEG(dat).nbchan], 'rmcomps', g.rmcomps{dat}, 'trialindices', g.trialindices{dat});
        EEG(dat).trials = size(EEG(dat).data,3);
        EEG(dat).event  = [];
        EEG(dat).epoch  = [];
        if ~isempty(g.interp),
            EEG(dat) = eeg_interp(EEG(dat), g.interp, 'spherical');
        end;
        tmpdata = EEG(dat).data;
    end;
    if all([ EEG.trials ] > 1)
        if size(X,2) ~= size(tmpdata,2) && size(X,3) ~= 1, error('Datasets to be concatenated do not have the same number of time points'); end;
        if isempty(X),
            X = tmpdata;
        else
            if size(X,1) ~= size(tmpdata,1), error('Datasets to be concatenated do not have the same number of channels'); end;
            X(:,:,end+1:end+size(tmpdata,3)) = tmpdata; % concatenating trials
        end;
    else
        % get boundaries for continuous data
        if ~isempty(EEG(dat).event) && isfield(EEG(dat).event, 'type') && ischar(EEG(dat).event(1).type)
             tmpevent = EEG(dat).event;
             tmpbound = strmatch('boundary', lower({ tmpevent.type }));
             if ~isempty(tmpbound)
                 boundaries = [boundaries size(X,2) [ tmpevent(tmpbound).latency ]-0.5+size(X,2) ];
             end;
        else
        end;
        if isempty(X),
            X = tmpdata;
        else
            if size(X,1) ~= size(tmpdata,1), error('Datasets to be concatenated do not have the same number of channels or components'); end;
            X(:,end+1:end+size(tmpdata,2)) = tmpdata;
        end;
    end;
end;
if ~isempty(boundaries), boundaries = [boundaries size(X,2)]; end;

% get specific time range for epoched and continuous data
if ~isempty(g.timerange)
    if oritrials > 1
        timebef  = find(EEG(1).times >= g.timerange(1) & EEG(1).times < g.timerange(2) );
        X        = X(:,timebef,:);
        EEG(1).pnts = length(timebef);
    else
        disp('warning: ''timerange'' option cannot be used with continuous data');
    end;
end;

% compute spectral decomposition
if strcmpi(g.specmode, 'psd')
    [X, f] = spectopo(X, size(X,2), EEG(1).srate, 'plot', 'off', 'boundaries', boundaries, 'nfft', g.nfft, spec_opt{:});
    if strcmpi(g.savetrials, 'on')
        disp('Cannot save trials using ''psd'' specmode option');
    end;
elseif strcmpi(g.specmode, 'pmtm')
    if all([ EEG.trials ] == 1) && ~isempty(boundaries), disp('Warning: multitaper does not take into account boundaries in continuous data'); end;
    fprintf('Computing multitaper:');
    for cind = 1:size(X,1)
        fprintf('.');
        for tind = 1:size(X,3)
            [tmpdat f] = pmtm(X(cind,:,tind), g.nw, g.nfft, EEG.srate);
            if cind == 1 && tind == 1
                X2 = zeros(size(X,1), length(tmpdat), size(X,3));
            end;
            X2(cind,:,tind) = tmpdat;
        end;
    end;
    fprintf('\n');
    X     = 10*log10(X2);
    if strcmpi(g.savetrials, 'off')
        X = mean(X,3);
    end;
elseif strcmpi(g.specmode, 'pburg')
    if all([ EEG.trials ] == 1) && ~isempty(boundaries), disp('Warning: pburg does not take into account boundaries in continuous data'); end;
    for cind = 1:size(X,1)
        fprintf('.');
        for tind = 1:size(X,3)
            [tmpdat f] = pburg(X(cind,:,tind), g.burgorder, g.nfft, EEG.srate);
            if cind == 1 && tind == 1
                X2 = zeros(size(X,1), length(tmpdat), size(X,3));
            end;
            X2(cind,:,tind) = tmpdat;
        end;
    end;
    fprintf('\n');
    if strcmpi(g.savetrials, 'off')
        X = mean(X,3);
    end;
else % fft mode
    if all([ EEG.trials ] == 1) && ~isempty(boundaries), disp('Warning: fft does not take into account boundaries in continuous data'); end;
    tmp   = fft(X, g.nfft, 2);
    f     = linspace(0, EEG(1).srate/2, floor(size(tmp,2)/2));
    f     = f(2:end); % remove DC (match the output of PSD)
    tmp   = tmp(:,2:floor(size(tmp,2)/2),:);
    X     = 10*log10(abs(tmp).^2);
    if strcmpi(g.savetrials, 'off')
        X = mean(X,3);
    end;
end;

% Save SPECs in file (all components or channels)
% ----------------------------------
options = { options{:} spec_opt{:} 'timerange' g.timerange 'nfft' g.nfft 'specmode' g.specmode };
if strcmpi(prefix, 'comp')
    savetofile( filename, f, X, 'comp', 1:size(X,1), options);
else
    if ~isempty(g.interp)
        savetofile( filename, f, X, 'chan', 1:size(X,1), options, { g.interp.labels });
    else
        tmpchanlocs = EEG(1).chanlocs;
        savetofile( filename, f, X, 'chan', 1:size(X,1), options, { tmpchanlocs.labels });
    end;
end;
return;

% -------------------------------------
% saving SPEC information to Matlab file
% -------------------------------------
function savetofile(filename, f, X, prefix, comps, params, labels);

    disp([ 'Saving SPECTRAL file ''' filename '''' ]);
    allspec = [];
    for k = 1:length(comps)
        allspec = setfield( allspec, [ prefix int2str(comps(k)) ], squeeze(X(k,:,:)));
    end;
    if nargin > 6
        allspec.labels = labels;
    end;
    allspec.freqs      = f;
    allspec.parameters = params;
    allspec.datatype   = 'SPECTRUM';
    allspec.average_spec = mean(X,1);
    std_savedat(filename, allspec);

author: PatternRecognition

repo: OpenBMI


text:
filename='ImprovedBridge_hexahedra';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volume'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 10;
Vfrac_final = 0.15;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function startIdx=findPatternInArray(pattern,array,maxTimes)
%%FINDPATTERNINARRAY Find all or a maximum number of occurrences of a
%                    pattern in an array of values.
%
%INPUTS: pattern A 1Xm or mX1 array holding elements of a pattern to be
%                found.
%        array   A 1Xn or nX1 array in which occurrences of the pattern are
%                to be found.
%       maxTimes An optional parameter specifying the maximum number of
%                occurrences of the pattern to find in array. If omitted or
%                an empty matrix is passed, then all occurrences of pattern
%                in array will be found.
%
%OUTPUTS: startIdx A numFoundX1 array of the atarting indices of the
%                occurrences of pattern in array, given in increasing
%                order.
%
%The algorithm used is a linear time method given in [1]. It is the
%algorithm of Section 2 with the modification to the next array of Section
%4 to find all occurrences of the pattern. Not all of the suggestions for
%improving efficiency of Section 3 are used, because they would require
%lengthening the input arrays and also would require defining certain
%values the inputs could not take.
%
%EXAMPLE:
% array='baabbabbaabaabbaabbaabbabaa';
% pattern='aabba';
% startIdx=findPatternInArray(pattern,array)
%One finds the string at startIdx=[2;12;16;20]. Note that some of the
%solutions overlap.
%
%REFERENCES:
%[1] D. E. Knuth, J. H. Morris Jr., and V. R. Pratty, "Fast pattern
%    matching in strings," SIAM Journal on Computing, vol. 6, no. 2, pp.
%    323-350, Jun. 1977.
%
%August 2015, David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

m=length(pattern);
n=length(array);

if(m>n)
    %If the pattern is longer than the array, then a match is impossible.
    startIdx=[];
    return;
end

if(nargin<3||isempty(maxTimes))
   maxTimes=n-m+1;%The maximum
else
   maxTimes=min([n-m+1,maxTimes]);
end

%Allocate space for the maximum possible number of matches.
startIdx=zeros(maxTimes,1);
numFound=0;

%First, compute the elements of the next table. We are computing them for a
%length m+1 table where the last element does not match anything else so
%that we can get a resume index to use when finding multiple matches.
next=zeros(m+1,1);
j=1;
t=0;
next(1)=0;
while(j<m)
   while((t>0)&&pattern(j)~=pattern(t))
       t=next(t);
   end
   t=t+1;
   j=j+1;
   if(pattern(j)==pattern(t))
       next(j)=next(j);
   else
        next(j)=t;
   end
end

%Add in the value for j=m;
while((t>0)&&pattern(j)~=pattern(t))
   t=next(t);
end
t=t+1;
j=j+1;
next(j)=t;%This sets next(m+1);

%Next, run the actual algorithm using the next table.
j=1;
k=1;

%Consider the special case of j=1 --skip forward to the first possible
%match.
while(array(k)~=pattern(1))
    k=k+1;
    if(k>n)%Nothing matches the first character of the pattern.
        startIdx=[];
        return;
    end
end

while(1)
    while((j<=m)&&(k<=n))
        while((j>0)&&(array(k)~=pattern(j)))
           j=next(j);
        end
        k=k+1;
        j=j+1;
    end

    if(j<=m)%All of the matches have been found.
        startIdx=startIdx(1:numFound);
        return;
    end

    %With j>m, the leftmost match has been found in positions k-m through
    %k-1.
    numFound=numFound+1;
    startIdx(numFound)=k-m;
    j=next(m+1);

    if(numFound==maxTimes)
        startIdx=startIdx(1:numFound);
        return;
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [errorcode, varargout] = common_size (varargin)

if (nargin < 2)
  error ('common_size: only makes sense if nargin >= 2');
end

% Find scalar args.
nscal = cellfun (@numel, varargin) ~= 1;

i = find (nscal, 1);

if (isempty (i))
  errorcode = 0;
  varargout = varargin;
else
  match = cellfun (@size_equal, varargin, repmat(varargin(i),size(varargin)));
  if (any (nscal & ~match))
    errorcode = 1;
    varargout = varargin;
  else
    errorcode = 0;
    if (nargout > 1)
      scal = ~nscal;
      varargout = varargin;
      dims = size (varargin{i});
      for s = find(scal)
        varargout{s} = repmat(varargin{s}, dims);
      end
    end
  end
end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function opt=varargin2struct(varargin)
%
% opt=varargin2struct('param1',value1,'param2',value2,...)
%   or
% opt=varargin2struct(...,optstruct,...)
%
% convert a series of input parameters into a structure
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% date: 2012/12/22
%
% input:
%      'param', value: the input parameters should be pairs of a string and a value
%       optstruct: if a parameter is a struct, the fields will be merged to the output struct
%
% output:
%      opt: a struct where opt.param1=value1, opt.param2=value2 ...
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

len=length(varargin);
opt=struct;
if(len==0) return; end
i=1;
while(i<=len)
    if(isstruct(varargin{i}))
        opt=mergestruct(opt,varargin{i});
    elseif(ischar(varargin{i}) && i<len)
        opt.(lower(varargin{i}))=varargin{i+1};
        i=i+1;
    else
        error('input must be in the form of ...,''name'',value,... pairs or structs');
    end
    i=i+1;
end


text:
function batch_data=batch_do_data_crop(train_opts, imdb, work_info_epoch, batch_data, work_info_batch, work_info)

data_crop_config=train_opts.data_crop_config;

if work_info.ref.run_eva || ~data_crop_config.do_crop
    batch_data.crop_info=[];
    return;
end

img_data=batch_data.img_data;
mask_data=batch_data.label_data;
img_size=size(img_data);
img_size=img_size(1:2);

if any(size(mask_data)~=img_size)
    mask_data=imresize(mask_data, img_size, 'nearest');
end

min_edge_size=min(img_size);
need_changed_flags=mod(min_edge_size,2)~=0;
min_edge_size(need_changed_flags)=min_edge_size(need_changed_flags)-1;

crop_box_size=data_crop_config.crop_box_size;
crop_box_size=min(min_edge_size, crop_box_size);

gen_crop_point_type=data_crop_config.gen_crop_point_type;

start_point=[];

if strcmp(gen_crop_point_type, 'random')
    start_point=gen_crop_point_random(data_crop_config, imdb, work_info_batch, mask_data, crop_box_size);

elseif strcmp(gen_crop_point_type, 'class_sample')

    start_point=gen_crop_point_class_sample(data_crop_config, imdb, work_info_batch, mask_data, crop_box_size);
    if isempty(start_point)
        start_point=gen_crop_point_random(data_crop_config, imdb, work_info_batch, mask_data, crop_box_size);
    end

else

	error('gen_crop_point_type not support!');
end

assert(~isempty(start_point));

stop_point1=start_point(1)+crop_box_size-1;
if stop_point1>img_size(1)
    stop_point1=img_size(1);
    start_point(1)=stop_point1-crop_box_size+1;
end

stop_point2=start_point(2)+crop_box_size-1;
if stop_point2>img_size(2)
    stop_point2=img_size(2);
    start_point(2)=stop_point2-crop_box_size+1;
end

row_idxes=start_point(1):stop_point1;
col_idxes=start_point(2):stop_point2;

crop_img_data=img_data(row_idxes, col_idxes, :);
crop_mask_data=mask_data(row_idxes, col_idxes);

crop_info=[];
crop_info.img_size=img_size;
crop_info.row_idxes=row_idxes;
crop_info.col_idxes=col_idxes;
crop_info.crop_mask_data=crop_mask_data;
crop_info.crop_img_data=crop_img_data;
crop_info.mask_data=mask_data;
crop_info.img_data=img_data;
crop_info.do_crop=true;

batch_data.crop_info=crop_info;

batch_data.img_data=crop_img_data;
batch_data.label_data=crop_mask_data;

crop_img_size=size(crop_img_data);
crop_img_size=crop_img_size(1:2);
batch_data.img_size=crop_img_size;

end

function start_point=gen_crop_point_random(data_crop_config, imdb, work_info_batch, mask_data, crop_box_size)

img_size=size(mask_data);

crop_box_step_ratio=data_crop_config.crop_box_step_ratio;
if ~isempty(crop_box_step_ratio)
    step_size=round(crop_box_size*crop_box_step_ratio);
    step_size=max(step_size, 1);
else
    step_size=1;
end

max_range=img_size-crop_box_size+1;

can_row_points=[1:step_size:max_range(1) max_range(1)];
can_col_points=[1:step_size:max_range(2) max_range(2)];
start_point=[my_random_sample(can_row_points, 1) my_random_sample(can_col_points, 1)];

end

function start_point=gen_crop_point_class_sample(data_crop_config, imdb, work_info_batch, mask_data, crop_box_size)

if ~isfield(imdb.ref, 'crop_cache_info') || isempty(imdb.ref.crop_cache_info)
    imdb.ref.crop_cache_info=gen_crop_cache_info(work_info_batch.ref.train_opts, imdb.ref.ds_info);
end

crop_cache_info=imdb.ref.crop_cache_info;
task_idx=work_info_batch.ref.task_idxes;
class_idxes=crop_cache_info.class_idxes_imgs{task_idx};

start_point=[];

if isempty(class_idxes)
    disp('################# class_sample crop: empty class_idxes!');
    return;
end

one_class_idx=my_random_sample(class_idxes, 1);
valid_flags=mask_data==one_class_idx;

valid_idxes=find(valid_flags);

if isempty(valid_idxes)
    disp('################# class_sample crop: no points with the selected class label are found!');
    return;
end

tmp_point_idx=my_random_sample(valid_idxes, 1);
[start_point(1), start_point(2)]=ind2sub(size(valid_flags), tmp_point_idx);

assert(~isempty(start_point));

start_point=round(start_point-crop_box_size./2);
start_point=max(start_point, 1);

end

function crop_cache_info=gen_crop_cache_info(train_opts, ds_info)

disp('gen_crop_cache_info...');
class_idxes_imgs=ds_info.class_idxes_imgs;

void_class_idxes=ds_info.class_info.void_class_idxes;
for img_idx=1:length(class_idxes_imgs)
    class_idxes=class_idxes_imgs{img_idx};
    assert(max(class_idxes)<256);
    if ~isempty(void_class_idxes)
        class_idxes=class_idxes(class_idxes~=void_class_idxes);
    end
    class_idxes_imgs{img_idx}=uint8(class_idxes);
end

crop_cache_info=[];
crop_cache_info.class_idxes_imgs=class_idxes_imgs;

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function RegOpticalFlow(handles, clipBox)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS;

indexS = planC{end};

    [originF, spacingF, centerF] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));

    [originM, spacingM, centerM] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationMovDataset));

    method = get(handles.para.optFlowMethod, 'value');
    RegLevels = str2double(get(handles.para.optFlowRegistLevel, 'string'));

    output = cell(1, 16);

    if isempty(clipBox)
        FixedImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
        MovingImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;
    else
        FixedImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        MovingImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
    end

    xyScale =  1; zScale = 1; DesampleTime = '0';
    %downsample the input datasets
    if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))

        tic;
        xyScale =  2; zScale = 2;
        spacingF = [spacingF(1)*xyScale spacingF(2)*xyScale spacingF(3)*zScale];
        spacingM = [spacingM(1)*xyScale spacingM(2)*xyScale spacingM(3)*zScale];
        disp('downSampling ...');

%         FImg = imdesample3d(FixedImg,xyScale,zScale);
%         MImg = imdesample3d(MovingImg,xyScale,zScale);
        FImg=GPReduce2(FixedImg,2,0);
        MImg=GPReduce2(MovingImg,2,0);

        toc;
        DesampleTime = num2str(ceil(toc/60));
    else
        FImg = FixedImg;
        MImg = MovingImg;
    end

%     call registration method
    switch method
        case 1
            methodIndex = 1;
        case 2
            methodIndex = 2;
        case 3
            methodIndex = 9;
        case 4
            methodIndex = 13;
        case 5
            methodIndex = 15;
        case 6
            methodIndex = 16;
        case 7
            methodIndex = 17;
        case 8
            methodIndex = 18;

    end

    tic;
    RegTime = '0';

    try
        [mvy,mvx,mvz,im]=multigrid_nogui6(          methodIndex, ... %method,
                                                    MImg, ...%Moving Image,
                                                    FImg, ...%Fixed Image,
                                                    spacingF,   ...%ratio,
                                                    RegLevels);    %steps,

        maxv = max(max(MImg(:)),max(FImg(:)));
        im = im * single(maxv);
        im = cast(im, class(FImg));

        toc;
        RegTime = num2str(ceil(toc/60));
    catch
        disp('Out of Memory, Registration failed, please restart Matlab!');
        return;
    end

% if downsampled, upsample the deformation fields.
    imIsDownsampled = 0;
    if ~(xyScale == 1 && zScale == 1)
        try
            [mvy,mvx,mvz] = recalculate_mvs(mvy,mvx,mvz);
            im = move3dimage(single(MovingImg),mvy,mvx,mvz); %deform the image with upsampled df fields.
        catch
            disp('Out of Memory, can not upscaling the motion field, system will save the downsampled field!');
            imIsDownsampled = 1;
        end
    end

% save the base image for analysis
    if ~isempty(clipBox) || imIsDownsampled

        planC{indexS.scan}(end+1) = planC{indexS.scan}(stateS.imageRegistrationBaseDataset);
        planC{indexS.scan}(end).scanArraySuperior = [];
        planC{indexS.scan}(end).scanArrayInferior = [];

        if (imIsDownsampled)
            bImg = FImg;
        else
            bImg = FixedImg;
        end
        planC{indexS.scan}(end).scanArray = int16(bImg);

        % setup scan info
        [xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));
        if ~isempty(clipBox)
            xRange = [xV(clipBox(1)) xV(clipBox(3))];
            yRange = [yV(clipBox(4)) yV(clipBox(2))];
        else
            xRange = [xV(1) xV(end)];
            yRange = [yV(end) yV(1)];
        end

        for i = 1:length(planC{indexS.scan}(end).scanInfo)
            planC{indexS.scan}(end).scanInfo(i).sizeOfDimension1 = size(bImg,1);
            planC{indexS.scan}(end).scanInfo(i).sizeOfDimension2 = size(bImg,2);
            %planC{indexS.scan}(end).scanInfo(i).CTOffset = -1*min(bImg(:));
            planC{indexS.scan}(end).scanInfo(i).xOffset = (xRange(1)+xRange(2))/2;
            planC{indexS.scan}(end).scanInfo(i).yOffset = (yRange(1)+yRange(2))/2;
        end
    end

% save the deformed image and deformation fields
    scanSetM = stateS.imageRegistrationMovDataset;
    scanSetF = stateS.imageRegistrationBaseDataset;

    planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
    planC{indexS.scan}(end).scanArraySuperior = [];
    planC{indexS.scan}(end).scanArrayInferior = [];
    planC{indexS.scan}(end).scanArray = int16(im);
    planC{indexS.scan}(end).deformInfo.BaseUID = planC{indexS.scan}(scanSetF).scanUID;
    planC{indexS.scan}(end).deformInfo.MoveUID = planC{indexS.scan}(scanSetM).scanUID;
    planC{indexS.scan}(end).deformInfo.DFX = mvx;
    planC{indexS.scan}(end).deformInfo.DFY = mvy;
    planC{indexS.scan}(end).deformInfo.DFZ = mvz;

    planC{indexS.scan}(end).scanUID = createUID('CT');

    planC{indexS.scan}(end).deformInfo.clipBox = clipBox;
    if ~(xyScale == 1 && zScale == 1)
        planC{indexS.scan}(end).deformInfo.downSampled = 'true';
    else
        planC{indexS.scan}(end).deformInfo.downSampled = 'false';
    end
% setup scan info
    [xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(scanSetF));
    xRange = [xV(clipBox(1)) xV(clipBox(3))];
    yRange = [yV(clipBox(4)) yV(clipBox(2))];
    for i = 1:length(planC{indexS.scan}(end).scanInfo)
        planC{indexS.scan}(end).scanInfo(i).sizeOfDimension1 = size(im,1);
        planC{indexS.scan}(end).scanInfo(i).sizeOfDimension2 = size(im,2);
        %planC{indexS.scan}(end).scanInfo(i).CTOffset = -1*min(im(:));
        planC{indexS.scan}(end).scanInfo(i).xOffset = (xRange(1)+xRange(2))/2;
        planC{indexS.scan}(end).scanInfo(i).yOffset = (yRange(1)+yRange(2))/2;
    end

% refresh view
    controlFrame('refresh');
    sliceCallBack('refresh');
    close('Registration Setup');

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function jacobi_eigenvalue_test03 ( )

%*****************************************************************************80
%
%% JACOBI_EIGENVALUE_TEST03 uses a 5x5 test matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'JACOBI_EIGENVALUE_TEST03\n' );
  fprintf ( 1, '  For a symmetric matrix A,\n' );
  fprintf ( 1, '  JACOBI_EIGENVALUE computes the eigenvalues D\n' );
  fprintf ( 1, '  and eigenvectors V so that A * V = D * V.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use the discretized second derivative matrix.\n' );

  n = 5;

  a = zeros ( n, n );
  for i = 1 : n
    a(i,i) = -2.0;
  end
  for i = 1 : n - 1
    a(i,i+1) = +1.0;
    a(i+1,i) = +1.0;
  end

  r8mat_print ( n, n, a, '  Input matrix A:' );

  it_max = 100;

  [ v, d, it_num, rot_num ] = jacobi_eigenvalue ( n, a, it_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', it_num );
  fprintf ( 1, '  Number of rotations  = %d\n', rot_num );

  r8vec_print ( n, d, '  Eigenvalues D:' );

  r8mat_print ( n, n, v, '  Eigenvector matrix V:' );
%
%  Compute eigentest.
%
  error_frobenius = r8mat_is_eigen_right ( n, n, a, v, d );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm error in eigensystem A*V-D*V = %g\n', ...
    error_frobenius );

  return
end


text:
function p_xy = cvt_circle ( r, np, p_xy, p_type )

%*****************************************************************************80
%
%% CVT_CIRCLE applies the CVT algorithm to points on a circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, integer NP, the number of points.
%
%    Input, real P_XY(2,NP), the point coordinates.
%
%    Input, integer P_TYPE(NP),
%    1, the point is constrained to the circle.
%    2, the point is constrainted to the disk.
%
%    Output, real P_XY(2,NP), the updated point coordinates.
%
  for i = 1 : 30
    p_xy = cvt_circle_step ( r, np, p_xy, p_type );
  end

  return
end


text:
function [g,dgdx,dgdp] = g_conv0(x,P,u,in)
% this function evaluates a linear convolution of the inputs
% function [gx,dgdx,dgdp] = g_conv0(x,P,u,in)
% This function evaluates the following multiple convolution model:
%   y(t) = cst + sum_k sum_tau w(k,tau)*u(k,t-tau)
% where u can be of any dimension (k=1,...,K).
% Note that the convolution (Volterra) kernels are constructed from the
% parameters' vector P. The size of the parameters' vector P is determined
% by the maximum lag (tau) considered.
% IN:
%   - x: [useless]
%   - P: (K*maxLag+1)x1 vector of kernel parameters.
%   - u: (K*nt)x1 vectorized input to the system
%   - in: [useless]
% OUT:
%   - g: the predicted system's output
%   - dgdx: [useless]
%   - dgdp: the gradient of the system's output w.r.t. kernel parameters
% SEE ALSO: g_convSig
if ~ isfield (in, 'K')
    in.K = 0;
end

try
    dgdp = in.dgdp;
catch
    nt = size(u,1)/in.dim.nu;
    dgdp = zeros(size(P,1),nt);
    for i=1:in.dim.nu
        ui = u((i-1)*nt+1:i*nt);
        for j=1:nt
            if j<=in.dim.n_t
                dgdp((i-1)*in.dim.n_t+j,:) = circshift(ui',[0,j-1]);
                dgdp((i-1)*in.dim.n_t+j,1:j-1) = 0;
            end
        end
    end
    dgdp(end,:) = 1;
end
g = dgdp'*P + in.K;
dgdx = [];

author: MBB-team

repo: VBA-toolbox


text:
function [ly] = nm2ly(nm)
% Convert length from nanometers to light years.
% Chad A. Greene 2012
ly = nm*1.057000834025e-25;


text:
function theta = polyhedronNormalAngle(varargin)
%POLYHEDRONNORMALANGLE Compute normal angle at a vertex of a 3D polyhedron.
%
%   THETA = polyhedraNormalAngle(NODES, EDGES, FACES, IND);
%   THETA = polyhedraNormalAngle(NODES, FACES, IND);
%   where NODES is a set of 3D points, and FACES a set of faces, whose
%   elements are indices to NODES array, compute the normal angle at the
%   vertex whose index is given by IND.
%
%   THETA = polyhedraNormalAngle(GRAPH, IND);
%   Uses a graph structure. GRAPH should contain at least fields : 'nodes'
%   and 'faces'.
%
%   Example :
%   % create a simple (irregular) tetrahedra
%   nodes = [0 0 0;1 0 0;0 1 0;0 0 1];
%   faces = [1 2 3;1 2 4;1 3 4;2 3 4];
%   % compute normal angle at each vertex
%   theta = polyhedronNormalAngle(nodes, faces, 1:size(nodes, 1));
%   % sum of normal angles should be equal to 4*pi :
%   sum(theta)
%
%   To-do
%       Works only for polyhedra with convex faces!!!
%
%   See also
%   polyhedra, polygon3dNormalAngle

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2005-11-30
% Copyright 2005-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

if length(varargin)==4
    nodes = varargin{1};
    faces = varargin{3};
    ind   = varargin{4};

elseif length(varargin)==3
    nodes = varargin{1};
    faces = varargin{2};
    ind   = varargin{3};

elseif length(varargin)==2
    graph = varargin{1};
    nodes = graph.nodes;
    faces = graph.faces;
    ind   = varargin{2};
else
    error('wrong number of arguments');
end

% number of angles to compute
na = length(ind);

theta = zeros(na, 1);
for i=1:na

    thetaf = [];

    % find faces containing given vertex,
    % and compute normal angle at each face containing vertex
    if iscell(faces)
        for j=1:length(faces)
            if ismember(ind(i), faces{j})
                % create 3D polygon
                face = nodes(faces{j}, :);

                % index of point in polygon
                indp = find(faces{j}==i);

                % compute normal angle of vertex
                thetaf = [thetaf polygon3dNormalAngle(face, indp)]; %#ok<AGROW>
            end
        end
    else
        indf = find(sum(ismember(faces, ind(i)), 2));

        thetaf = zeros(length(indf), 1);
        for j=1:length(indf)
            ind2 = faces(indf(j), :);
            face = nodes(ind2, :);
            indp = find(ind2==ind(i));
            thetaf(j) = pi - polygon3dNormalAngle(face, indp);
        end
    end

    % compute normal angle of polyhedron, by use of angle defect formula
    theta(i) = 2*pi - sum(thetaf);

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function OUT = VARmakelags(DATA,lag)
% =======================================================================
% Builds a matrix with lagged values of DATA, i.e. if DATA = [x1 x2],
% VARmakelags(DATA,1) yields OUT = [x1 x2 x1(-1) x2(-1)]. Serves as an
% input to VARmakexy
% =======================================================================
% OUT = VARmakelags(DATA,lag)
% -----------------------------------------------------------------------
% INPUT
%   - DATA: matrix containing the original data
%   - lag: lag order
% -----------------------------------------------------------------------
% OUTPUT
%   - OUT: matrix of lagged values
% -----------------------------------------------------------------------
% EXAMPLE
%   x = [1 2; 3 4; 5 6; 7 8; 9 10];
%   OUT = VARmakelags(x,2)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

% Get dimesion of DATA
[nobs, ~]= size(DATA);

% Create the lagged matrix
OUT = [];
for jj=0:lag-1
    OUT = [DATA(jj+1:nobs-lag+jj,:), OUT];
end

% Finally, save the non-lagged values...
aux = DATA(lag+1:end,:);

%... and append to the lagged matrix
OUT = [aux OUT];

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [hrf, fit, e, param] = Fit_sFIR(tc,TR,Run,T,mode)
% function [hrf, fit, e, param] = Fit_sFIR(tc,TR,Runs,T,mode)
%
% Fits FIR and smooth FIR model
%
% INPUTS:
%
% tc    - time course
% TR    - time resolution
% Runs  - expermental design
% T     - length of estimated HRF
% mode  - FIR or smooth FIR
%   options:
%       0 - standard FIR
%       1 - smooth FIR
%
% OUTPUTS:
%
% hrf   - estimated hemodynamic response function
% fit   - estimated time course
% e     - residual time course
% param - estimated amplitude, height and width
%
% Created by Martin Lindquist on 10/02/09
% Last edited: 05/26/10 (ML)

numstim = length(Run);
len = length(Run{1});

Runs = zeros(len,numstim);
for i=1:numstim,
    Runs(:,i) = Run{i};
end;

[DX] = tor_make_deconv_mtx3(Runs,T,1);

% DX2 = DX(:,1:T);
% num = T;

if mode == 1

    C=(1:T)'*(ones(1,T));
    h = sqrt(1/(7/TR));                       % 7 seconds smoothing - ref. Goutte

    v = 0.1;
    sig = 1;

    R = v*exp(-h/2*(C-C').^2);
    RI = inv(R);
    MRI = zeros(numstim*T+1);
    for i=1:numstim,
        MRI(((i-1)*T+1):(i*T),((i-1)*T+1):(i*T)) = RI;
    end;

    b = inv(DX'*DX+sig^2*MRI)*DX'*tc;
    fit = DX*b;
    e = tc - DX*b;

elseif mode == 0

    b = pinv(DX)*tc;
    fit = DX*b;
    e = tc - DX*b;

end

hrf =zeros(T,numstim);
param = zeros(3,numstim);

for i=1:numstim,
    hrf(:,i) = b(((i-1)*T+1):(i*T))';
    param(:,i) = get_parameters2(hrf(:,i),T);
end;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function SDP = rmlindep_sdpt3(SDP)
%% remove linearly dependent constraints
[Atnew,bnew,~,~,~,~,~] = ...
              checkdepconstr(SDP.blk,SDP.At,SDP.b,zeros(length(SDP.b),1),1);
SDP.At  = Atnew;
SDP.b   = bnew;
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function [dat] = ft_preproc_rectify(dat)

% FT_PREPROC_RECTIFY rectifies the data, i.e. converts all samples with a
% negative value into the similar magnitude positive value
%
% Use as
%   [dat] = ft_preproc_rectify(dat)
% where
%   dat        data matrix (Nchans X Ntime)
%
% See also PREPROC

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

dat = abs(dat);

author: spm


text:
function tests = TwistTest
  tests = functiontests(localfunctions);

  clc
end

function twist2d_test(tc)
    %%2D twists

    % check basics work
    s = [1 2 3];
    tw = Twist(s);
    verifyEqual(tc, tw.S, s', 'absTol', 1e-6);
    verifyEqual(tc, tw.v', s(1:2), 'absTol', 1e-6);
    verifyEqual(tc, tw.w, s(3), 'absTol', 1e-6);
    verifyEqual(tc, tw.se, [skew(s(3)) [s(1:2)]'; 0 0 0], 'absTol', 1e-6);
end

function operator2d_test(tc)
    % check overloaded *

    tw = Twist([1 2 3]);
    tw2 = Twist([4 6 5]);

    both = tw * tw2;

    Tboth = tw.T * tw2.T;
    verifyEqual(tc, both.T, Tboth, 'absTol', 1e-6);

    % check rotational twist
    tw = Twist('R', [1 2]);
    verifyEqual(tc, tw.S, [2 -1 1]', 'absTol', 1e-6);

    % check prismatic twist
    tw = Twist('P', [2 3]);
    verifyEqual(tc, tw.S, [unit([2 3]) 0]', 'absTol', 1e-6);
    tw = Twist('T', [2 3]);
    verifyEqual(tc, tw.S, [unit([2 3]) 0]', 'absTol', 1e-6);

    % check twist from SE(2)
    tw = Twist( trot2(0) );
    verifyEqual(tc, tw.S, [0 0 0]', 'absTol', 1e-6);
    tw = Twist( trot2(pi/2) );
    verifyEqual(tc, tw.S, [0 0 pi/2]', 'absTol', 1e-6);
    tw = Twist( SE2(1,2,0) );
    verifyEqual(tc, tw.S, [1 2 0]', 'absTol', 1e-6);
    tw = Twist( SE2(1,2,pi/2) );
    verifyEqual(tc, tw.S, [3*pi/4 pi/4 pi/2]', 'absTol', 1e-6);

    % test expm and T
    verifyEqual(tc, tw.T, trexp2(tw.S), 'absTol', 1e-6);
    verifyEqual(tc, tw.exp.double, trexp2(tw.S), 'absTol', 1e-6);

    tw = Twist('R', [1 2]);
    verifyEqual(tc, tw.T(pi/2), [0 -1 3; 1 0 1; 0 0 1], 'absTol', 1e-6);
end

function operator3d_test(tc)
    %% 3D twists

    % check basics work
    s = [1 2 3 4 5 6];
    tw = Twist(s);
    verifyEqual(tc, tw.S, s', 'absTol', 1e-6);
    verifyEqual(tc, tw.v', s(1:3), 'absTol', 1e-6);
    verifyEqual(tc, tw.w', s(4:6), 'absTol', 1e-6);
    verifyEqual(tc, tw.se, [skew(s(4:6)) [s(1:3)]'; 0 0 0 0], 'absTol', 1e-6);

    % check overloaded *
    s2 = [4 6 5 7 9 8];
    tw2 = Twist(s2);

    both = tw * tw2;

    Tboth = tw.T * tw2.T;
    verifyEqual(tc, double(both.T), double(Tboth), 'absTol', 1e-6);

    % check rotational twist
    tw = Twist('R', [1 2 3], [0 0 0]);
    verifyEqual(tc, tw.S, [0 0 0 unit([1 2 3])]', 'absTol', 1e-6);

    % check prismatic twist
    tw = Twist('P', [1 2 3]);
    verifyEqual(tc, tw.S, [unit([1 2 3]) 0 0 0 ]', 'absTol', 1e-6);
    tw2 = Twist('T', [1 2 3]);
    verifyEqual(tc, tw.S, tw2.S, 'absTol', 1e-6);

    % check twist from SE(3)
    tw = Twist( trotx(0) );
    verifyEqual(tc, tw.S, [0 0 0  0 0 0]', 'absTol', 1e-6);
    tw = Twist( trotx(pi/2) );
    verifyEqual(tc, tw.S, [0 0 0  pi/2 0 0]', 'absTol', 1e-6);
    tw = Twist( troty(pi/2) );
    verifyEqual(tc, tw.S, [0 0 0  0 pi/2 0]', 'absTol', 1e-6);
    tw = Twist( trotz(pi/2) );
    verifyEqual(tc, tw.S, [0 0 0  0 0 pi/2]', 'absTol', 1e-6);

    tw = Twist( transl([1 2 3]) );
    verifyEqual(tc, tw.S, [1 2 3  0 0 0]', 'absTol', 1e-6);
    tw = Twist( transl([1 2 3])*troty(pi/2) );
    verifyEqual(tc, tw.S, [-pi/2 2 pi  0 pi/2 0]', 'absTol', 1e-6);

    % test expm and T
    verifyEqual(tc, tw.T, SE3.exp(tw.S).T);
    verifyEqual(tc, tw.exp, SE3.exp(tw.S));

    tw = Twist('R', [1 0 0], [0 0 0]);
    verifyEqual(tc, tw.T(pi/2), SE3.Rx(pi/2).T, 'absTol', 1e-6);
    tw = Twist('R', [0 1 0], [0 0 0]);
    verifyEqual(tc, tw.T(pi/2), SE3.Ry(pi/2).T, 'absTol', 1e-6);
    tw = Twist('R', [0 0 1], [0 0 0]);
    verifyEqual(tc, tw.T(pi/2), SE3.Rz(pi/2).T, 'absTol', 1e-6);
end

function char_test(tc)
    % 2d

    tw = Twist([1 2 3]);

    s = char(tw);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 1);
    s = char([tw tw tw]);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 3);
    s = char([tw tw tw]');
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 3);

    % 3d
    tw = Twist([4 6 5 7 9 8]);
    s = char(tw);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 1);
    s = char([tw tw tw]);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 3);
    s = char([tw tw tw]');
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 3);
end

author: petercorke

repo: spatialmath-matlab


text:
% Copyright 2017 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
%
%     https://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.

% Returns set of parameters
% set light_mode = true to run the code in a sub optimal but faster mode
% set light_mode = false to obtain the results reported in the RED paper

function params = GetGaussianDeblurSDParams(light_mode)

% regularization factor
params.lambda = 0.01;
% number of outer iterations
if light_mode
    params.outer_iters = 500;
else
    params.outer_iters = 1500;
end
% level of noise assumed in the regularization-denoiser
params.effective_sigma = 4.1;

return

author: google

repo: RED

save_path: github-repos/MATLAB/google-RED


text:
function [ NRI , pval ] = NetReclassificationImprovement( pred_old , pred_new , outcome )
% Net Reclassification Index (NRI)
% Compare classification of an old vs new classification technique
% Usage:
%   [ NRI , pval ] = NetReclassificationImprovement( pred_old , pred_new , outcome )
%
% http://www.epibiostat.ucsf.edu/courses/RoadmapK12/SIGS/Pencina.pdf
% Statist. Med. (in press) (www.interscience.wiley.com) DOI: 10.1002/sim.2929
% Evaluating the added predictive ability of a new marker: From area under the ROC curve to reclassi�cation and beyond
% Michael J. Pencina Ralph B. D'agostino Sr Ralph B. D'Agostino Jr and Ramachandran S. Vasan
%
% (c) Louis Mayaud, 2011 (louis.mayaud@gmail.com)
% Please reference :
% Mayaud, Louis, et al. "Dynamic Data During Hypotensive Episode Improves
%        Mortality Predictions Among Patients With Sepsis and Hypotension*."
%        Critical care medicine 41.4 (2013): 954-962.

% Remove NaNs
Idx = (isnan(pred_old) & isnan(pred_new));
pred_old(Idx) = [];
pred_new(Idx) = [];

% Need to define the following probabilities
PEventUp = sum(outcome==1 & pred_old==0 & pred_new==1)/sum(outcome==1);
PEventDown = sum(outcome==1 & pred_old==1 & pred_new==0)/sum(outcome==1);
PNoneventUp= sum(outcome==0 & pred_old==0 & pred_new==1)/sum(outcome==0);
PNoneventDown= sum(outcome==0 & pred_old==1 & pred_new==0)/sum(outcome==0);

NRI = (PEventUp - PEventDown) - (PNoneventUp - PNoneventDown) ;
z = abs(NRI)/sqrt( (PEventUp + PEventDown)/sum(outcome==1) + (PNoneventUp + PNoneventDown)/sum(outcome==0)  );
[~,pval] = ztest(z,0,1);

% Implemented from Mc Nemar 1947, as detailed in
% Biometrics, 66, 1185-1191, 2010 Dec.
% Westfall et al.
% N01 = sum(pred_old==0 & pred_new==1);
% N10 = sum(pred_old==1 & pred_new==0);
% Nd = N01 + N10 ;
%
% pval = binocdf(N01,Nd,0.5);


text:
%   Distribution code Version 1.0 -- 02/31/2020 by Wei Liu Copyright 2020
%
%   The code is created based on the method described in the following paper
%   [1] "Real-time Image Smoothing via Iterative Least Squares", Wei Liu, Pingping Zhang,
%        Xiaolin Huang, Jie Yang, Chunhua Shen and Ian Reid, ACM Transactions on Graphics,
%        presented at SIGGRAPH 2020.
%
%   The code and the algorithm are for non-comercial use only.

%  ---------------------- Input------------------------
%  F:              input image, can be gray image or RGB color image
%  lambda:     \lambda in Eq.(1), control smoothing strength
%  gamma:       the \gamma in the Welsch's penalty in Eq. (18)
%  iter:           iteration number of the ILS

%  ---------------------- Output------------------------
%  U:             smoothed image

function U =ILS_Welsch_GPU(F, lambda, gamma, iter)

F = gpuArray(single(F)); % 'single' precision is very important to reduce the computational cost

c = 2;

[N, M, D] = size(F);
sizeI2D = [N, M];

otfFx = psf2otf_Dx_GPU(sizeI2D); % equal to otfFx = psf2otf(fx, sizeI2D) where fx = [1, -1];
otfFy = psf2otf_Dy_GPU(sizeI2D); % equal to otfFy = psf2otf(fy, sizeI2D) where fy = [1; -1];

Denormin = abs(otfFx).^2 + abs(otfFy ).^2;
Denormin = repmat(Denormin, [1, 1, D]);
Denormin = 1 + 0.5 * c * lambda * Denormin;

U = F;  % smoothed image

Normin1 = fft2(U);

for k = 1: iter

    % Intermediate variables \mu update, in x-axis and y-axis direction
    u_h = [diff(U,1,2), U(:,1,:) - U(:,end,:)];
    u_v = [diff(U,1,1); U(1,:,:) - U(end,:,:)];

    mu_h = c .* u_h - 2 .* u_h .* exp(- u_h .* u_h / (2 * gamma^2));
    mu_v = c .* u_v - 2 .* u_v .* exp(- u_v .* u_v / (2 * gamma^2));

    % Update the smoothed image U
    Normin2_h = [mu_h(:,end,:) - mu_h(:, 1,:), - diff(mu_h,1,2)];
    Normin2_v = [mu_v(end,:,:) - mu_v(1, :,:); - diff(mu_v,1,1)];

    FU = (Normin1 + 0.5 * lambda * (fft2(Normin2_h + Normin2_v))) ./ Denormin;
    U = real(ifft2(FU));

    Normin1 = FU;  % This helps to further enlarge the smoothing strength

end

U = gather(U);

author: wliusjtu


text:
function ResampledFiberGroup=dtiResampleFiberGroup(fg, value, flag)
%ResampledFiberGroup=resampleFibers(fiberGroup, numNodes, 'N')
%ResampledFiberGroup=resampleFibers(fiberGroup, stepSize, 'L')
%Resamples the fibers so that each has number of nodes = value (if flag=='N')
%or to fibers of step size = value (flag=='L')

%ER 02/2008 SCSNL
warning off;

if ~exist('flag', 'var') || isempty(flag)
flag='N';
end

ResampledFiberGroup=fg;

numFibers=size(fg.fibers,1);

for fiberIndex=1:numFibers
ResampledFiberGroup.fibers{fiberIndex}=dtiFiberResample(fg.fibers{fiberIndex}, value, flag);
end
 warning on;


text:
function [ddtData] = get_all_from_ddt(filename)

ddtData  =[];
[dirstr, name, ext] = fileparts(filename);
ddtData.FileName = [name ext];

[ddtData.raw.nch, ddtData.raw.npoints, ddtData.raw.freq, ddtData.raw.d] = ddt(filename);
[ddtData.v.nch, ddtData.v.npoints, ddtData.v.freq, ddtData.v.d] = ddt_v(filename);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [x, w] = trigpts(n)
%TRIGPTS   Equispaced points in [-1, 1).
%   TRIGPTS(N) returns N equispaced points in [-1, 1).
%
%   [X, W] = TRIGPTS(N) returns also a row vector of the weights for
%   the trapezoidal rule.
%
% See also CHEBPTS, LEGPTS, JACPTS, LAGPTS, and HERMPTS.

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

% Special case (no points).
if ( n <= 0 )
    x = [];
    w = [];
    return

end

x = linspace(-1, 1, n+1).';
x(end) = [];

% Quadrature weights:
if ( nargout > 1 )
    w = trigtech.quadwts(n);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [au] = ft2au(ft)
% Convert length from feet to astronomical units.
% Chad A. Greene 2012
au = ft*2.03746215499e-12;

author: Sable

repo: mcbench-benchmarks


text:
function I = argmax(X, DIM)
%ARGMAX   Find the index of the maximum value of a matrix or vector along
% the dimension specified (DIM).
%
%   Inputs:
%       X                       = Input vector or matrix.
%       DIM                     = Dimension along which to computer maximum argument.
%
%   Outputs:
%       I                       = Index of the maximum value.
%
% Daniel McDuff, Ethan Blackford, January 2019
% Copyright (c)
% Licensed under the MIT License and the RAIL AI License.

[~,I] = max(X,[],DIM);

author: danmcduff

repo: iphys-toolbox


text:
%% t_mrdViewFibers
%
% Load and view a set of fibers stored in a PDB (v3) file using mrMesh.
%
% (c) Stanford VISTA Team

%% Make sure that vistadata is on your path
% vistaDataPath;

%% dtiStart for scripts

dataDir = fullfile(mrvDataRootPath,'diffusion','sampleData');
dt6Name = fullfile(dataDir,'dti40','dt6.mat');
% Initialize mrDiffusion
[dtiFig, dtiH]= mrDiffusion('on',dt6Name);
% At this point you can use dtiSet/Get on dtiH

if isunix, disp('Start mrMeshSrv.exe')
else       mrmStart
end

%% Load fiber group structure
difDir = fullfile(mrvDataRootPath,'diffusion','sampleData','dti40');
chdir(difDir)

fgName = fullfile(mrvDataRootPath,'diffusion','sampleData','fibers','leftArcuate.pdb');
% The mtr<> function name is old from metrotrac days.  It will probably change some
% day.
fg = mtrImportFibers(fgName);

%% Show fibers in mrMesh

% Attach the fiber group to the handles
dtiH = dtiSet(dtiH,'add fiber group',fg);

set(dtiH.cbUseMrMesh, 'Value',1);     % Use mrMesh for update
set(dtiH.cbShowFibers,'Value',1);     % Show loaded fibers
set(dtiH.cbShowMatlab3d,'Value',0);   % No 3d Matlab window
set(dtiH.popupBackground,'Value',2);  % Mean diffusivity

guidata(dtiFig,dtiH);  % Refresh the Matlab window handles.

mrmCloseWindow(dtiH.mrMesh.id,dtiH.mrMesh.host);

showMeshWindow = 1;
dtiH = dtiRefreshFigure(dtiH,showMeshWindow);

%% The code below pulls out key routines from dtiRefreshFigure

% Closing the window and then running this is much faster than removing the
% previous actors.
mrmCloseWindow(dtiH.mrMesh.id,dtiH.mrMesh.host);

% This code, which is not normally used, permits a faster reload of the
% mrMesh window. The code below here could be placed in a routine like
%
%   id = dtiMeshView(dtiH,varargin);
%
set(dtiH.popupBackground,'Value',1);  % Choose type of background (1-4)
set(dtiH.rbSagittal,'Value',1);       % Choose which planes
set(dtiH.rbCoronal, 'Value',0);
set(dtiH.rbAxial,   'Value',1);

% With mrMesh - This started with the code from dtiRefreshFogure/dtiFiberUI
% Needs anat, anatXform
[xSliceRgb,ySliceRgb,zSliceRgb,anat,anatXform, ...
    mmPerVoxel,xform,xSliceAxes,ySliceAxes,zSliceAxes] = ...
    dtiGetCurSlices(dtiH);

curPosition = dtiGet(dtiH,'curpos');

% Should be:
% anatXform = dtiGet(dtiH,'anatXform');

[zIm] = dtiGetSlice(anatXform, anat, 3, curPosition(3), [], dtiH.interpType);
[yIm] = dtiGetSlice(anatXform, anat, 2, curPosition(2), [], dtiH.interpType);
[xIm] = dtiGetSlice(anatXform, anat, 1, curPosition(1), [], dtiH.interpType);
% figure; imagesc(zIm); axis image; colormap(gray)

% This is a little slow.
[xIm,yIm,zIm] = dtiMrMeshSelectImages(dtiH,xIm,yIm,zIm);
origin = dtiMrMeshOrigin(dtiH);
dtiH = dtiMrMesh3AxisImage(dtiH, origin, xIm, yIm, zIm);

%%  Manipulating the mesh view

% Need to permute the rotation matrix, I think, for various cases.
% I haven't worked that out.

msh = dtiGet(dtiH,'mesh'); % If you have a mrDiffusion guidata in dtiH
mrmRotateCamera(msh.id,'front',1); pause(1)
mrmRotateCamera(msh.id,'back',1);  pause(1)
mrmRotateCamera(msh.id,'top',1);   pause(1)
mrmRotateCamera(msh.id,'bottom',1);pause(1)
mrmRotateCamera(msh.id,'right',1); pause(1)
mrmRotateCamera(msh.id,'left',1);

%% End after here

author: vistalab

repo: vistasoft


text:
function order = strong_elim_order(G, node_sizes, partial_order)
% STRONG_ELIM_ORDER Find an elimination order to produce a strongly triangulated graph.
% order = strong_elim_order(moral_graph, node_sizes, partial_order)
%
% partial_order(i,j)=1 if we must marginalize i *after* j
% (so i will be nearer the strong root).
% e.g., if j is a decision node and i is its information set:
%   we cannot maximize j if we have marginalized out some of i
% e.g., if j is a continuous child and i is its discrete parent:
%   we want to integrate out the cts nodes before the discrete ones,
%   so that the marginal is strong.
%
% For details, see
% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
% - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical
%   association models", JASA 87(420):1098--1108, 1992.
%
% On p369 of the Jensen paper, they state "the reverse of the elimination order must be some
% extension of [the partial order] to a total order".
% We make no attempt to find the best such total ordering, in the sense of minimizing the weight
% of the resulting cliques.

% Example from the Jensen paper:
% Let us number the nodes in Fig 1 from top to bottom, left to right,
% so a=1,b=2,D1=3,c=4,...,l=14,j=15,k=16.
% The elimination ordering they propose on p370 is [14 15 16 11 12 1 4 5 10 8 13 9 7 6 3 2];

if 0
  total_order = topological_sort(partial_order);
  order = total_order(end:-1:1); % no attempt to find an optimal constrained ordering!
  return;
end

% The following implementation is due to Ilya Shpitser and seems to give wrong
% results on cg1

n = length(G);
MG = G; % copy the original graph
uneliminated = ones(1,n);
order = zeros(1,n);

for i=1:n
  roots = [];
  k = 1;
  for j=1:n
    if sum(partial_order(j,:)) == 0
      roots(k) = j;
      k = k + 1;
    end
  end
  U = find(uneliminated);
  valid = myintersect(U, roots);
  % Choose the best node from the set of valid candidates
  score1 = zeros(1,length(valid));
  score2 = zeros(1,length(valid));
  for j=1:length(valid)
    k = valid(j);
    ns = myintersect(neighbors(G, k), U);
    l = length(ns);
    M = MG(ns,ns);
    score1(j) = l^2 - sum(M(:)); % num. added edges
    score2(j) = prod(node_sizes([k ns])); % weight of clique
  end
  j1s = find(score1==min(score1));
  j = j1s(argmin(score2(j1s)));
  k = valid(j);
  uneliminated(k) = 0;
  order(i) = k;
  ns = myintersect(neighbors(G, k), U);
  if ~isempty(ns)
    G(ns,ns) = 1;
    G = setdiag(G,0);
  end
  partial_order(:,k) = 0;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [newMap] = colorSubsystemCD(map, model, subsystem, color, areaWidth)
% Color and increase areaWidth of every reaction in a specific subsystem
%
% USAGE:
%
%   [newMap] = colorSubsystemCD(map, model, subsystem, color, areaWidth);
%
% INPUTS:
%   map:            File from CellDesigner parsed to MATLAB format
%   model:          COBRA model structure
%   subsystem:      Name of a subsystem as a String
%
% OPTIONAL INPUTS:
%   color:          Color desired for reactions in CAPITALS
%   areaWidth:          Width desired for reactions
%
% OUTPUT:
%   newMap          MATLAB structure of map with reaction modifications
%
% .. Authors:
%       - A.Danielsdottir 17/07/2017 LCSB. Belval. Luxembourg
%       - N.Sompairac - Institut Curie, Paris, 11/10/2017.

    if nargin < 5
        areaWidth = 8;
    end
    if nargin < 4
        color = 'RED';
    end

    newMap = map;
    rxnList = model.rxns(ismember([model.subSystems{:}]', subsystem));
    colors = createColorsMap;

    index = find(ismember(newMap.rxnName, rxnList));
    for j = index'
        newMap.rxnColor{j, 1} = colors(color);
        newMap.rxnWidth{j, 1} = areaWidth;
    end

    % Use the existence of reactant lines to check if the map has the
    % complete structure, and if so change also secondary lines.
    if any(strcmp('rxnReactantLineColor', fieldnames(newMap))) == 1
        for j = index'
            if ~isempty(newMap.rxnReactantLineColor{j})
                for k = 1:length(newMap.rxnReactantLineColor{j})
                    newMap.rxnReactantLineColor{j, 1}{k, 1} = colors(color);
                    newMap.rxnReactantLineWidth{j, 1}{k, 1} = areaWidth;
                end
            end
            if ~isempty(newMap.rxnProductLineColor{j})
                for m = 1:1:length(newMap.rxnProductLineColor{j})
                    newMap.rxnProductLineColor{j, 1}{m, 1} = colors(color);
                    newMap.rxnProductLineWidth{j, 1}{m, 1} = areaWidth;
                end
            end
        end
    end

end

author: opencobra

repo: cobratoolbox