text:
function [ncst,Area,k]=mu_coast(optn,varargin)
% MU_COAST Add a coastline to a given map.
%         MU_COAST draw a coastline as either filled patches (slow) or
%         lines (fast) on a given projection. It uses a coastline database with
%         a resolution of about 1/4 degree.
%
%         It is not meant to be called directly by the user, instead it contains
%         code common to various functions in the m_map directory.
%
%         Calling sequence is MU_COAST(option,arguments) where
%
%         Option string
%           c,l,i,h,f :  Accesses various GSHHS databases. Next argument is
%                        GSHHS filename.
%
%           u[ser] :   Accesses user-specified coastline file (a mat-file of
%                      data saved from a previous MU_COAST call). Next argument
%                      is filename
%
%           v[ector] : Uses input vector of data. Next argument is the
%                     data in the form of a nx2 matrix of [longitude latitude].
%                     Patches must have first and last points the same. In a vector,
%                     different patches can be separated by NaN.
%
%           d[efault] :  Accesses default coastline.
%
%         The arguments given above are then followed by (optional) arguments
%         specifying lines or patches, in the form:
%
%          optional arguments:  <line property/value pairs>, or
%                               'line',<line property/value pairs>.
%                               'patch',<patch property/value pairs>.
%                               'speckle',width,density,<line property/value pairs>.
%
%         If no or one output arguments are specified then the coastline is drawn, with
%         patch handles returned.
%         If 3 output arguments are specified in the calling sequence, then no drawing
%         is done. This can be used to save subsampled coastlines for future use
%         with the 'u' option, for fast drawing of multiple instances of a particular
%         coastal region.
%
%
%
%         See also M_PROJ, M_GRID

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997
%
% Notes: 15/June/98 - fixed some obscure problems that occured sometimes
%                     when using conic projections with large extents that
%                     crossed the 180-deg longitude line.
%        31/Aug/98  - added "f" gshhs support (Thanks to RAMO)
%        17/June/99 - 'line' option as per manual (thanks to Brian Farrelly)
%        3/Aug/01   - kludge fix for lines to South Pole in Antarctica (response
%                     to Matt King).
%        30/May/02  - fix to get gshhs to work in Antarctica.
%        15/Dec/05  - speckle additions
%        21/Mar/06  - handling of gshhs v1.3 (developed from suggestions by
%                     Martin Borgh)
%        26/Nov/07 - changed 'finite' to 'isfinite' after warnings
%        26/Sep/14 - added hierarchy flag to borders and rivers
%        13/Nov/14 - suggested matlab2014b graphics fix
%         8/Mar/17 - conic wrap-around not quite working right in a weird case when
%                    'rectbox' was on.
%           Nov/17 - redid a lot of the gshhs-related stuff
%
% This software is provided "as is" without warranty of any kind. But
% its mine, so you can't sell it.

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

% m_coasts.mat contains 3 variables:
% ncst: a Nx2 matrix of [LONG LAT] line segments, each of which form
%       a closed contour (i.e. first and last points are the SAME),
%       separated by NaN
%  k=[find(isnan(ncst(:,1)))];
% Area: a vector giving the areas of the different patches. Both ncst
%     and Area should be ordered with biggest regions first. Area can
%     be computed as follows:
%
%  Area=zeros(length(k)-1,1);
%  for i=1:length(k)-1,
%    x=ncst([k(i)+1:(k(i+1)-1) k(i)+1],1);
%    y=ncst([k(i)+1:(k(i+1)-1) k(i)+1],2);
%    nl=length(x);
%    Area(i)=sum( diff(x).*(y(1:nl-1)+y(2:nl))/2 );
%  end;
%
%     Area should be >0 for land, and <0 for lakes and inland seas.

switch optn(1)
  case {'c','l','i','h','f'}
    [ncst,k,Area]=get_coasts(optn,varargin{1});
    varargin(1)=[];
  case  'u'
    eval(['load ' varargin{1} ' -mat']);
    varargin(1)=[];
  case 'v'
    ncst=[NaN NaN;varargin{1};NaN NaN];
    varargin(1)=[];
    k=[find(isnan(ncst(:,1)))];  % Get k
    Area=ones(size(k));          % Make dummy Area vector (all positive).
    otherwise
    currloc=mfilename('fullpath');  % Octave has problems finding the file
    load([currloc(1:end-8) 'm_coasts.mat']);

end

% If all we wanted was to extract a sub-coastline, return.
if nargout==3
  return;
end

% Handle wrap-arounds (not needed for azimuthal and oblique projections)

switch MAP_PROJECTION.routine
 case {'mp_cyl','mp_conic','mp_tmerc'}
  if MAP_VAR_LIST.longs(2)<-180
   ncst(:,1)=ncst(:,1)-360;
  elseif MAP_VAR_LIST.longs(1)>180
   ncst(:,1)=ncst(:,1)+360;
  elseif MAP_VAR_LIST.longs(1)<-180
   Area=[Area;Area];
   k=[k;k(2:end)+k(end)-1];
   ncst=[ncst;[ncst(2:end,1)-360 ncst(2:end,2)]];
   % This is kinda kludgey - but sometimes adding all these extra points
   % causes wrap-around in the conic projection, so we want to limit the
   % longitudes to the range needed. However, we don't just clip them to
   % min long because that can cause problems in trying to decide which way
   % curves are oriented when doing the fill algorithm below. So instead
   % I sort of crunch the scale, preserving topology.
   %
   % 12/Sep/2006 - in the gsshs_crude database we have a lot of long skinny
   % things which interact badly with this - so I offset the scrunch 2 degrees
   % away from the bdy
   %
   % Mar/2017 - OK, so somebody used a rectbox 'on' which gave a really wide range of
   % lat/long, and it turns out I need to scale the *other* side as well!
   nn=ncst(:,1)>MAP_VAR_LIST.longs(2)+2;                                          % added 2017
   ncst(nn,1)=(ncst(nn,1)-MAP_VAR_LIST.longs(2)-2)/100+MAP_VAR_LIST.longs(2)+2;   %  "
   nn=ncst(:,1)<MAP_VAR_LIST.longs(1)-2;
   ncst(nn,1)=(ncst(nn,1)-MAP_VAR_LIST.longs(1)+2)/100+MAP_VAR_LIST.longs(1)-2;
  elseif MAP_VAR_LIST.longs(2)>180
   Area=[Area;Area];
   k=[k;k(2:end)+k(end)-1];
   ncst=[ncst;[ncst(2:end,1)+360 ncst(2:end,2)]];
   % Ditto.
   nn=ncst(:,1)>MAP_VAR_LIST.longs(2)+2;
   ncst(nn,1)=(ncst(nn,1)-MAP_VAR_LIST.longs(2)-2)/100+MAP_VAR_LIST.longs(2)+2;
   nn=ncst(:,1)<MAP_VAR_LIST.longs(1)-2;                                          % added 2017
   ncst(nn,1)=(ncst(nn,1)-MAP_VAR_LIST.longs(1)+2)/100+MAP_VAR_LIST.longs(1)-2;   % "
 %%%  disp('beep');
  end
end

if ~isempty(varargin)
  if strcmp(varargin(1),'patch')
    optn='patch';
  end
  if strcmp(varargin(1),'speckle')
    optn='speckle';
  end
  if strcmp(varargin(1),'line')
    optn='line';
    varargin=varargin(2:end); % ensure 'line' does not get passed to line
  end
else
  optn='line';
end

switch optn
 case {'patch','speckle'}

  switch MAP_VAR_LIST.rectbox
    case 'on'
      xl=MAP_VAR_LIST.xlims;
      yl=MAP_VAR_LIST.ylims;
      [X,Y]=m_ll2xy(ncst(:,1),ncst(:,2),'clip','on');
      oncearound=4;
    case 'off'
      xl=MAP_VAR_LIST.longs;
      yl=MAP_VAR_LIST.lats;
      X=ncst(:,1);
      Y=ncst(:,2);
      [X,Y]=mu_util('clip','on',X,xl(1),X<xl(1),Y);
      [X,Y]=mu_util('clip','on',X,xl(2),X>xl(2),Y);
      [Y,X]=mu_util('clip','on',Y,yl(1),Y<yl(1),X);
      [Y,X]=mu_util('clip','on',Y,yl(2),Y>yl(2),X);
      oncearound=4;
    case 'circle'
      rl=MAP_VAR_LIST.rhomax;
      [X,Y]=m_ll2xy(ncst(:,1),ncst(:,2),'clip','on');
      oncearound=2*pi;
  end

  if ~MAP_PROJECTION.newgraphics
     p_hand=zeros(length(k)-1,1); % Patch handles
  else
     p_hand=gobjects(length(k)-1,1); % Patch handles
  end

  for i=1:length(k)-1
    x=X(k(i)+1:k(i+1)-1);
    fk=isfinite(x);
    if any(fk)
      y=Y(k(i)+1:k(i+1)-1);XX=x;YY=y;
%% if i>921, disp('pause 1'); pause; end;
      nx=length(x);
      if Area(i)<0, x=flipud(x);y=flipud(y); fk=flipud(fk); end
%clf
%line(x,y,'color','m');
      st=find(diff(fk)==1)+1;
      ed=find(diff(fk)==-1);
%length(x),
%ed,
%st
      if length(st)<length(ed) || isempty(st), st=[ 1;st]; end
      if length(ed)<length(st),                ed=[ed;nx]; end
%ed
%st
      if  ed(1)<st(1)
        if c_edge(x(1),y(1))==9999
          x=x([(ed(1)+1:end) (1:ed(1))]);
          y=y([(ed(1)+1:end) (1:ed(1))]);
          fk=isfinite(x);
          st=find(diff(fk)==1)+1;
          ed=[find(diff(fk)==-1);nx];
          if length(st)<length(ed), st=[1;st]; end
        else
          ed=[ed;nx];
          st=[1;st];
        end
      end
%ed
%st
      % Get rid of 2-point curves (since often these are out-of-range lines with
      % both endpoints outside the region of interest)
      k2=(ed-st)<3;
      ed(k2)=[]; st(k2)=[];
%%[XX,YY]=m_ll2xy(x(st),y(st),'clip','off');
%line(x,y,'color','r','linestyle','-');
%line(x(st),y(st),'marker','o','color','r','linestyle','none');
%line(x(ed),y(ed),'marker','o','color','g','linestyle','none');
      edge1=c_edge(x(st),y(st));
      edge2=c_edge(x(ed),y(ed));
%-edge1*180/pi
%-edge2*180/pi
      mi=1;
      while ~isempty(st)
        if mi==1
          xx=x(st(1):ed(1));
          yy=y(st(1):ed(1));
        end
        estart=edge2(1);
        s_edge=edge1;
        s_edge(s_edge<estart)=s_edge(s_edge<estart)+oncearound;
%s_edge,estart
        [md,mi]=min(s_edge-estart);
        switch MAP_VAR_LIST.rectbox
          case {'on','off'}
            for e=floor(estart):floor(s_edge(mi))
              if e==floor(s_edge(mi)), xe=x(st([mi mi])); ye=y(st([mi mi]));
              else,  xe=xl; ye=yl; end
              switch rem(e,4)
                case 0
                  xx=[xx; xx(end*ones(10,1))];
                  yy=[yy; yy(end)+(ye(2)-yy(end))*[1:10]'/10 ];
                case 1
                  xx=[xx; xx(end)+(xe(2)-xx(end))*[1:10]'/10 ];
                  yy=[yy; yy(end*ones(10,1))];
                case 2
                  xx=[xx; xx(end*ones(10,1))];
                  yy=[yy; yy(end)+(ye(1)-yy(end))*[1:10]'/10;];
                case 3
                  xx=[xx; xx(end)+(xe(1)-xx(end))*[1:10]'/10 ];
                  yy=[yy; yy(end*ones(10,1))];
              end
            end
          case 'circle'
            if estart~=9999
%s_edge(mi),estart
              xx=[xx; rl*cos(-[(estart:.1:s_edge(mi))]')];
              yy=[yy; rl*sin(-[(estart:.1:s_edge(mi))]')];
            end
        end
        if mi==1 % joined back to start
    	   if strcmp(optn,'patch')
               if strcmp(MAP_VAR_LIST.rectbox,'off')
                 [xx,yy]=m_ll2xy(xx,yy,'clip','off');
               end
               if Area(i)<0
                  p_hand(i)=patch(xx,yy,varargin{2:end},'facecolor',get(gca,'color'));
                  if ishandle(p_hand(i)), set(p_hand(i),'tag',[ get(p_hand(i),'tag') '_lake']);end
              else
                  p_hand(i)=patch(xx,yy,varargin{2:end});
                  if ishandle(p_hand(i)),set(p_hand(i),'tag',[ get(p_hand(i),'tag') '_land']);end
               end

	     else  % speckle
            if ~strcmp(MAP_VAR_LIST.rectbox,'off')  % If we were clipping in
	           [xx,yy]=m_xy2ll(xx,yy);                % screen coords go back
            end
            if Area(i)>0
               p_hand(i)=m_hatch(xx,yy,'speckle',varargin{2:end});
               if ishandle(p_hand(i)),set(p_hand(i),'tag',[get(p_hand(i),'tag') '_lake']); end
            else
	           p_hand(i)=m_hatch(xx,yy,'outspeckle',varargin{2:end});
               if ishandle(p_hand(i)),set(p_hand(i),'tag',[ get(p_hand(i),'tag') '_land']); end
            end
         end
%%%if i>921, disp(['paused-2 ' int2str(i)]);pause; end;
         ed(1)=[];st(1)=[];edge2(1)=[];edge1(1)=[];
       else
         xx=[xx;x(st(mi):ed(mi))];
         yy=[yy;y(st(mi):ed(mi))];
         ed(1)=ed(mi);st(mi)=[];ed(mi)=[];
         edge2(1)=edge2(mi);edge2(mi)=[];edge1(mi)=[];
       end
%%if i>921, disp(['paused-2 ' int2str(i)]);pause; end;
     end
    end
   end

   otherwise

  % This handles the odd points required at the south pole by any Antarctic
  % coastline by setting them to NaN (for lines only)
  ii=ncst(:,2)<=-89.9;
  if any(ii), ncst(ii,:)=NaN; end

  [X,Y]=m_ll2xy(ncst(:,1),ncst(:,2),'clip','on');

  % Get rid of 2-point lines (these are probably clipped lines spanning the window)
  fk=isfinite(X);
  st=find(diff(fk)==1)+1;
  ed=find(diff(fk)==-1);
  k=find((ed-st)==1);
  X(st(k))=NaN;

  p_hand=line(X,Y,varargin{:});

end

ncst=p_hand;

end
%-----------------------------------------------------------------------
function edg=c_edge(x,y)
% C_EDGE tests if a point is on the edge or not. If it is, it is
%        assigned a value representing it's position oon the perimeter
%        in the clockwise direction. For x/y or lat/long boxes, these
%        values are
%           0 -> 1 on left edge
%           1 -> 2 on top
%           2 -> 3 on right edge
%           3 -> 4 on bottom
%        For circular boxes, these values are the -ve of the angle
%        from center.

global MAP_VAR_LIST

switch MAP_VAR_LIST.rectbox
  case 'on'
    xl=MAP_VAR_LIST.xlims;
    yl=MAP_VAR_LIST.ylims;
  case 'off'
    xl=MAP_VAR_LIST.longs;
    yl=MAP_VAR_LIST.lats;
  case 'circle'
    rl2=MAP_VAR_LIST.rhomax^2;
end

edg=9999+zeros(length(x),1);
tol=1e-10;

switch MAP_VAR_LIST.rectbox
  case {'on','off'}
    i=abs(x-xl(1))<tol;
    edg(i)=(y(i)-yl(1))/diff(yl);

    i=abs(x-xl(2))<tol;
    edg(i)=2+(yl(2)-y(i))/diff(yl);

    i=abs(y-yl(1))<tol;
    edg(i)=3+(xl(2)-x(i))/diff(xl);

    i=abs(y-yl(2))<tol;
    edg(i)=1+(x(i)-xl(1))/diff(xl);

  case 'circle'
    i=abs(x.^2 + y.^2 - rl2)<tol;
    edg(i)=-atan2(y(i),x(i));   % use -1*angle so that numeric values
                                % increase in CW direction
end

end

function [ncst,k,Area]=get_coasts(optn,file)
%
%  GET_COASTS  Loads various GSHHS coastline databases and does some preliminary
%              processing to get things into the form desired by the patch-filling
%              algorithm.
%
% Changes; 3/Sep/98 - RP: better decimation, fixed bug in limit checking.
%          Nov/17  - mostly rewritten

global MAP_PROJECTION MAP_VAR_LIST

llim=rem(MAP_VAR_LIST.longs(1)+360,360)*1e6;
rlim=rem(MAP_VAR_LIST.longs(2)+360,360)*1e6;
tlim=MAP_VAR_LIST.lats(2)*1e6;
blim=MAP_VAR_LIST.lats(1)*1e6;

m.rlim=rem(MAP_VAR_LIST.longs(2)+360+180,360)-180;
m.llim=rem(MAP_VAR_LIST.longs(1)+360+180,360)-180;
m.tlim=MAP_VAR_LIST.lats(2);
m.blim=MAP_VAR_LIST.lats(1);

% decfac is for decimation of areas outside the lat/long bdys.
% Sizes updated for gshhs v1.3, checked for v2.3.6
% Sizes updated for v1.10, and for river/border databases.
switch optn(1)
  case 'f'  % 'full' (undecimated) database
    ncst=NaN(10810000,2);Area=zeros(188611,1);k=ones(188612,1);
    decfac=12500;
  case 'h'
    ncst=NaN(2080000,2);Area=zeros(153545,1);k=ones(153546,1);
   decfac=2500;
  case 'i'
    ncst=NaN(493096,2);Area=zeros(41529,1);k=ones(41530,1);
    decfac=500;
  case 'l'
    ncst=NaN(124871,2);Area=zeros(20776,1);k=ones(27524,1);
    decfac=100;
  case 'c'
    ncst=NaN(14403,2);Area=zeros(1766,1);k=ones(1767,1);
    decfac=20;
end

flaglim=9;
if length(optn)>=2
  flaglim=str2num(optn(2));
end

fid=fopen(file,'r','ieee-be');

if fid==-1
  warning(['Coastline file ' file ...
          ' not found \n(Have you installed it? See the M_Map User''s Guide for details)' ...
          '\n ---Using default coastline instead']);
  load m_coasts
  return
end

Area2=Area;

% Read the File header

[cnt,g]=get_gheader(fid);

l=0;
while cnt>0

   C=fread(fid,g.N*2,'int32'); % Read all points in the current segment.

   % For versions > 12 there are 2 Antarctics - ice-line and grounding line,
   % also they changed the limits from 0-360 to -180 to 180.
   if g.ver>14
     switch g.level
       case 6
         g.level=flaglim+1; % ice line - don't show
       case 5
         g.level=1;         % grounding line - use this.
     end
   else
     %For Antarctica the lime limits are 0 to 360 (exactly), thus c==0 and the
     %line is not chosen for (e.g. a conic projection of part of Antarctica)
     % Fix 30may/02
     if g.extentE==360e6, g.extentE=g.extentE-1; end
   end

   a=rlim>llim;  % Map limits cross longitude jump? (a==1 is no)
   b=g.greenwich<65536; % Cross boundary? (b==1 if no).
   c=llim<rem(g.extentE+360e6,360e6);
   d=rlim>rem(g.extentW+360e6,360e6);
   e=tlim>g.extentS & blim<g.extentN;

   % This test checks whether the lat/long box containing the line overlaps that of
   % the map. There are various cases to consider, depending on whether map limits
   % and/or the line limits cross the longitude jump or not.

   if     e  ...
      && (    (  a&&( (b&&c&&d) || (~b&&(c||d) ) ) )    ...
           || ( ~a&&(    ~b     || ( b&&(c||d) ) ) )  ) ...
      && g.level<=flaglim

     l=l+1;
     x=C(1:2:end)*1e-6;
     y=C(2:2:end)*1e-6;

     %  line(x,y,'color','r');pause;

     %  make things continuous (join edges that cut across 0-meridian)

     dx=diff(x);
     x=x-360*cumsum([x(1)>180;(dx>356) - (dx<-356)]);

     % Antarctic is a special case - extend contour to make nice closed polygon
     % that doesn't surround the pole.
     if g.ver<15
       % Range from 0-360
       if abs(x(1))<1 && abs(y(1)+68.9)<1
         y=[-89.9;-78.4;y(x<=-180);y(x>-180);   -78.4;-89.9*ones(18,1)];
         x=[  180; 180 ;x(x<=-180)+360;x(x>-180);-180; [-180:20:160]'];
       end
     else
       % new version is from -180 to 180
       if abs(x(1))==180 && y(1)<-76
         y=[-89.9;-78.4;y;   -78.4;-89.9*ones(19,1)];
         x=[ 180; 180 ;x;    -180; [-180:20:180]'];
       end
     end

     %plot(x,y);pause;

     % First and last point should be the same IF THIS IS A POLYGON
     % if the Area=0 then this is a line, and don't add points!

     if g.area>0

       if x(end)~=x(1) || y(end)~=y(1)   % First and last points should be the same
             x=[x;x(1)];y=[y;y(1)];
       end

       % get correct curve orientation for patch-fill algorithm.

       Area2(l)=sum( diff(x).*(y(1:(end-1))+y(2:end))/2 ); % Area "on the page"
       Area(l)=g.area/10;                                  % Area "on the globe"

       if rem(g.level,2)==0  % Make lakes (2) and islands (1) differently oriented
         Area(l)=-abs(Area(l));
         if Area2(l)>0, x=x(end:-1:1);y=y(end:-1:1); end
       else
         if Area2(l)<0, x=x(end:-1:1);y=y(end:-1:1); end
       end

     else
       % Later on 2 point lines are clipped so we want to avoid that
       if length(x)==2
          x=[x(1);mean(x);x(2)];
          y=[y(1);mean(y);y(2)];
       end
     end

     % Here we try to reduce the number of points and clip them
     % close to map limits (if they are too far away this can cause
     % problems with some projections as they wrap around into the
     % wrong place)

     [cx,cy]=clip_to_lims(x,y,m,decfac);

     k(l+1)=k(l)+length(cx)+1;
     ncst(k(l)+1:k(l+1),:)=[cx,cy;NaN NaN];

     % This is a little tricky...the filling algorithm expects data to be in the
     % range -180 to 180 deg long. However, there are some land parts that cut across
     % this divide so they appear at +190 but not -170. This causes problems later...
     % so as a kludge I replicate some of the problematic features at 190-360=-170.
     % I cut out the problematic points only, make sure the first and last
     % points are the same, and then add this curve to ncst
     if max(x)>180
        ii=find(x>179);ii=[ii;ii(1)];
        [cx,cy]=clip_to_lims(x(ii)-360,y(ii),m,decfac);
        l=l+1;
        k(l+1)=k(l)+length(cx)+1;
        ncst(k(l)+1:k(l+1),:)=[cx,cy;NaN NaN];
     end

     %plot(ncst(:,1),ncst(:,2));drawnow;
   end
  [cnt,g]=get_gheader(fid);
end

fclose(fid);

%plot(ncst(:,1),ncst(:,2));pause;clf;

ncst((k(l+1)+1):end,:)=[];  % get rid of unused part of data matrices
Area((l+1):end)=[];
k((l+2):end)=[];
%%%fprintf('Size ncst: %d, Area: %d, k %d\n',length(ncst),length(Area),length(k))

end

function [x,y]=clip_to_lims(x,y,m,decfac)
% Look for points outside the lat/long boundaries, and then decimate them
   % by a factor of about 'decfac' (don't get rid of them completely because that
   % can sometimes cause problems when polygon edges cross curved map edges).

   tol=.2;

   % Do y limits, then x so we can keep corner points.

   nn=(y>m.tlim+tol) | (y<m.blim-tol);
     % keep one extra point when crossing limits, also the beginning/end point.
   nn=logical(nn-min(1,([0;diff(nn)]>0)+([diff(nn);0]<0)));
   nn([1 end])=0;
     % decimate vigorously
   nn=nn & rem(1:length(nn),decfac)'~=0;
   x(nn)=[];y(nn)=[];

   if m.rlim>m.llim  % no wraparound
       % sections of line outside lat/long limits
     nn=(x>m.rlim+tol | x<m.llim-tol) & y<m.tlim & y>m.blim;
    else            % wraparound case
     nn=(x>m.rlim+tol & x<m.llim-tol ) & y<m.tlim & y>m.blim;
   end
   nn=logical(nn-min(1,([0;diff(nn)]>0)+([diff(nn);0]<0)));nn([1 end])=0;
   nn=nn & rem(1:length(nn),decfac)'~=0;
   x(nn)=[];y(nn)=[];

   % Move all points "near" to map boundaries.
   % I'm not sure about the wisdom of this - it might be better to clip
   % to the boundaries instead of moving. Hmmm.

   y(y>m.tlim+tol)=m.tlim+tol;
   y(y<m.blim-tol)=m.blim-tol;
   if m.rlim>m.llim   % Only clip long bdys if I can tell I'm on the right
                     % or left (i.e. not in wraparound case)
     x(x>m.rlim+tol)=m.rlim+tol;
     x(x<m.llim-tol)=m.llim-tol;
   end

end

function [cnt,g]=get_gheader(fid)
% Reads the gshhs file header
%
% A bit of code added because header format changed with version 1.3.
%
% 17/Sep/2008 - added material to handle latest GSHHS version.
%
% For version 1.1 this is the header ( 9*4 = 36 bytes long)
%
%int id;				/* Unique polygon id number, starting at 0 */
%int n;				/* Number of points in this polygon */
%int level;			/* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
%int west, east, south, north;	/* min/max extent in micro-degrees */
%int area;			/* Area of polygon in 1/10 km^2 */
%short int greenwich;		/* Greenwich is 1 if Greenwich is crossed */
%short int source;		/* 0 = CIA WDBII, 1 = WVS */
%
% For version 1.3 of GMT format was changed to this ( 10*4 = 40 bytes long)
%
%int id;				/* Unique polygon id number, starting at 0 */
%int n;				/* Number of points in this polygon */
%int level;			/* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
%int west, east, south, north;	/* min/max extent in micro-degrees */
%int area;			/* Area of polygon in 1/10 km^2 */
%int version;			/* Polygon version, set to 3
%short int greenwich;		/* Greenwich is 1 if Greenwich is crossed */
%short int source;		/* 0 = CIA WDBII, 1 = WVS */
%
% For version 1.4, we have (8*4 = 32 bytes long)
%
%int id;				/* Unique polygon id number, starting at 0 */
%int n;				/* Number of points in this polygon */
%int flag;			/* level + version << 8 + greenwich << 16 + source << 24
%int west, east, south, north;	/* min/max extent in micro-degrees */
%int area;			/* Area of polygon in 1/10 km^2 */
%
%Here, level, version, greenwhich, and source are
%level:		1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
%version:	Set to 4 for GSHHS version 1.4
%greenwich:	1 if Greenwich is crossed
%source:		0 = CIA WDBII, 1 = WVS
%
% For version 2.0 it all changed again, we have (11*4 = 44 bytes)
%
%	int id;		/* Unique polygon id number, starting at 0 */
%	int n;		/* Number of points in this polygon */
%	int flag;	/* = level + version << 8 + greenwich << 16 + source << 24 + river << 25 */
%	/* flag contains 5 items, as follows:
%	 * low byte:	level = flag & 255: Values: 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
%                                             For Antarctic starting with 2.3.0 ice-front=5 and grounding line=6
%                                           for border database: 1=country, 2=state/province
%	 * 2nd byte:	version = (flag >> 8) & 255: Values: Should be 7 for GSHHS release 7 (i.e., version 2.0)
%                                                                     12 for version 2.2
%                                                                     15 for 2.3.6
%	 * 3rd byte:	greenwich = (flag >> 16) & 1: Values: Greenwich is 1 if Greenwich is crossed
%	 * 4th byte:	source = (flag >> 24) & 1: Values: 0 = CIA WDBII, 1 = WVS
%	 * 4th byte:	river = (flag >> 25) & 1: Values: 0 = not set, 1 = river-lake and level = 2
%	 */
%	int west, east, south, north;	/* min/max extent in micro-degrees */
%	int area;	/* Area of polygon in 1/10 km^2 */
%	int area_full;	/* Area of original full-resolution polygon in 1/10 km^2 */
%	int container;	/* Id of container polygon that encloses this polygon (-1 if none) */
%	int ancestor;	/* Id of ancestor polygon in the full resolution set that was the source of this polygon (-1 if none)

% Now, in the calling code I have to use A(2),A(3),A(5-7), A(8), A(9) from original.

[A,cnt]=fread(fid,8,'int32');

if cnt<8  % This gets triggered by the EOF
  g=[];
  return;
end

g.ver=bitand(bitshift(A(3),-8),255);

if g.ver==0  % then its an old version, fake it to look new

  % This works for version 1.2, but not 1.3.

  [A2,cnt2]=fread(fid,1,'int32');
  A=[A;A2];

  if (cnt+cnt2)==9 && A(9)==3  % we have version 1.3, this would be one of 0,1,65535,65536 in v 1.2
    % Read one more byte
    A2=fread(fid,1,'int32');
    % This is the easiest way not to break existing code.
    A(9)=A2;   % one of 0,1,65536,65537
  end

elseif g.ver>=7  % After v2.0 some more bytes around for original area and container/ancerstor IDs

  A2=fread(fid,3,'int32');

end

% a newest versions
g.id=A(1);
g.N=A(2);
g.level=bitand(A(3),255);
g.greenwich=bitand(bitshift(A(3),-16),255)*65536;
g.source=bitand(bitshift(A(3),-24),255);
A(3)=g.level;
A(9)=g.greenwich;

g.extentW=A(4);
g.extentE=A(5);
g.extentS=A(6);
g.extentN=A(7);
g.area=A(8);

end

author: goGPS-Project


text:
function u = cosh(a)
%COSH         Slope hyperbolic cosine cosh(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  u = a;

  u.r = cosh(a.r);
  u.s = slopeconvexconcave('cosh','sinh(%)',a,1);

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function neuron = run_cnmfe_matlab(nam)

    %% make a sources object to store results
    neuron = Sources2D();
    nam = neuron.select_data(nam);

    %% don't ask about anything
    choose_params = false; % manually choose parameters

    %% parameters
    % -------------------------    COMPUTATION    -------------------------  %
    pars_envs = struct('memory_size_to_use', 20, ...   % GB, memory space you allow to use in MATLAB
        'memory_size_per_patch', 5, ...   % GB, space for loading data within one patch
        'patch_dims', [64, 64]);  %GB, patch size

    % -------------------------      SPATIAL      -------------------------  %
    gSig = 3;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
    gSiz = 10;          % pixel, neuron diameter
    ssub = 1;           % spatial downsampling factor
    with_dendrites = true;   % with dendrites or not
    if with_dendrites
        % determine the search locations by dilating the current neuron shapes
        updateA_search_method = 'dilate';  %#ok<UNRCH>
        updateA_bSiz = 5;
        updateA_dist = neuron.options.dist;
    else
        % determine the search locations by selecting a round area
        updateA_search_method = 'ellipse'; %#ok<UNRCH>
        updateA_dist = 5;
        updateA_bSiz = neuron.options.dist;
    end
    spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'
    spatial_algorithm = 'hals_thresh';

    % -------------------------      TEMPORAL     -------------------------  %
    Fs = 30;             % frame rate
    tsub = 1;           % temporal downsampling factor
    deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
        'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
        'smin', -1, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level (was -2)
        'optimize_pars', true, ...  % optimize AR coefficients
        'optimize_b', true, ...% optimize the baseline);
        'max_tau', 100);    % maximum decay time (unit: frame);

    nk = 3;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
    % when changed, try some integers smaller than total_frame/(Fs*30)
    detrend_method = 'spline';  % compute the local minimum as an estimation of trend.

    % -------------------------     BACKGROUND    -------------------------  %
    bg_model = 'ring';  % model of the background {'ring', 'svd'(default), 'nmf'}
    nb = 1;             % number of background sources for each patch (only be used in SVD and NMF model)
    bg_neuron_factor = 1.4;
    ring_radius = round(bg_neuron_factor * gSiz);  % when the ring model used, it is the radius of the ring used in the background model.
    %otherwise, it's just the width of the overlapping area
    num_neighbors = 50; % number of neighbors for each neuron

    % -------------------------      MERGING      -------------------------  %
    show_merge = false;  % if true, manually verify the merging step
    merge_thr = 0.65;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
    method_dist = 'max';   % method for computing neuron distances {'mean', 'max'}
    dmin = 4;       % minimum distances between two neurons. it is used together with merge_thr
    dmin_only = 1;  % merge neurons if their distances are smaller than dmin_only.
    merge_thr_spatial = [0.8, 0.4, -inf];  % merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)

    % -------------------------  INITIALIZATION   -------------------------  %
    K = [];             % maximum number of neurons per patch. when K=[], take as many as possible.
    min_corr = 0.4;     % minimum local correlation for a seeding pixel
    min_pnr = 5;       % minimum peak-to-noise ratio for a seeding pixel
    min_pixel = gSig^2;      % minimum number of nonzero pixels for each neuron
    bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
    frame_range = [];   % when [], uses all frames
    save_initialization = false;    % save the initialization procedure as a video.
    use_parallel = true;    % use parallel computation for parallel computing
    show_init = true;   % show initialization results
    center_psf = true;  % set the value as true when the background fluctuation is large (usually 1p data)
    % set the value as false when the background fluctuation is small (2p)

    % -------------------------  Residual   -------------------------  %
    min_corr_res = 0.7;
    min_pnr_res = 6;
    seed_method_res = 'auto';  % method for initializing neurons from the residual
    update_sn = true;

    % ----------------------  WITH MANUAL INTERVENTION  --------------------  %
    with_manual_intervention = false;

    % -------------------------  FINAL RESULTS   -------------------------  %
    save_demixed = true;    % save the demixed file or not
    kt = 3;                 % frame intervals

    % -------------------------    UPDATE ALL    -------------------------  %
    neuron.updateParams('gSig', gSig, ...       % -------- spatial --------
        'gSiz', gSiz, ...
        'ring_radius', ring_radius, ...
        'ssub', ssub, ...
        'search_method', updateA_search_method, ...
        'bSiz', updateA_bSiz, ...
        'dist', updateA_bSiz, ...
        'spatial_constraints', spatial_constraints, ...
        'spatial_algorithm', spatial_algorithm, ...
        'tsub', tsub, ...                       % -------- temporal --------
        'deconv_options', deconv_options, ...
        'nk', nk, ...
        'detrend_method', detrend_method, ...
        'background_model', bg_model, ...       % -------- background --------
        'nb', nb, ...
        'ring_radius', ring_radius, ...
        'num_neighbors', num_neighbors, ...
        'merge_thr', merge_thr, ...             % -------- merging ---------
        'dmin', dmin, ...
        'method_dist', method_dist, ...
        'min_corr', min_corr, ...               % ----- initialization -----
        'min_pnr', min_pnr, ...
        'min_pixel', min_pixel, ...
        'bd', bd, ...
        'center_psf', center_psf);
    neuron.Fs = Fs;

    %% distribute data and be ready to run source extraction
    neuron.getReady(pars_envs);

    %% initialize neurons from the video data within a selected temporal range
    if choose_params
        % change parameters for optimized initialization
        [gSig, gSiz, ring_radius, min_corr, min_pnr] = neuron.set_parameters();
    end

    [center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);
    neuron.compactSpatial();
    if show_init
        figure();
        ax_init= axes();
        imagesc(Cn, [0, 1]); colormap gray;
        hold on;
        plot(center(:, 2), center(:, 1), '.r', 'markersize', 10);
    end

    %% estimate the background components
    neuron.update_background_parallel(use_parallel);
    neuron_init = neuron.copy();

    %%  merge neurons and update spatial/temporal components
    neuron.merge_neurons_dist_corr(show_merge);
    neuron.merge_high_corr(show_merge, merge_thr_spatial);

    %% update spatial components

    % %% pick neurons from the residual
    % [center_res, Cn_res, PNR_res] =neuron.initComponents_residual_parallel([], save_initialization, use_parallel, min_corr_res, min_pnr_res, seed_method_res);
    % if show_init
    %     axes(ax_init);
    %     plot(center_res(:, 2), center_res(:, 1), '.g', 'markersize', 10);
    % end
    % neuron_init_res = neuron.copy();
    if size(neuron.A, 2) == 0
        disp('NO NEURONS FOUND WITH CURRENT INITIALIZATION SETTINGS')
        return
    end

    %% udpate spatial&temporal components, delete false positives and merge neurons
    % update spatial
    if update_sn
        neuron.update_spatial_parallel(use_parallel, true);
        udpate_sn = false;
    else
        neuron.update_spatial_parallel(use_parallel);
    end
    % merge neurons based on correlations
    neuron.merge_high_corr(show_merge, merge_thr_spatial);

    for m=1:2
        % update temporal
        neuron.update_temporal_parallel(use_parallel);

        % delete bad neurons
        neuron.remove_false_positives();

        % merge neurons based on temporal correlation + distances
        neuron.merge_neurons_dist_corr(show_merge);
    end

    %% add a manual intervention and run the whole procedure for a second time
    neuron.options.spatial_algorithm = 'nnls';
    if with_manual_intervention
        show_merge = true;
        neuron.orderROIs('snr');   % order neurons in different ways {'snr', 'decay_time', 'mean', 'circularity'}
        neuron.viewNeurons([], neuron.C_raw);

        % merge closeby neurons
        neuron.merge_close_neighbors(true, dmin_only);

        % delete neurons
        tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
        ids = find(tags>0);
        if ~isempty(ids)
            neuron.viewNeurons(ids, neuron.C_raw);
        end
    end
    %% run more iterations
    neuron.update_background_parallel(use_parallel);
    neuron.update_spatial_parallel(use_parallel);
    neuron.update_temporal_parallel(use_parallel);

    K = size(neuron.A,2);
    tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
    neuron.remove_false_positives();
    neuron.merge_neurons_dist_corr(show_merge);
    neuron.merge_high_corr(show_merge, merge_thr_spatial);

    if K~=size(neuron.A,2)
        neuron.update_spatial_parallel(use_parallel);
        neuron.update_temporal_parallel(use_parallel);
        neuron.remove_false_positives();
    end

    %% save the workspace for future analysis
    neuron.orderROIs('snr');
    cnmfe_path = neuron.save_workspace();

    %% show neuron contours
    Coor = neuron.show_contours(0.6);

    %% create a video for displaying the
    %amp_ac = 140;
    %range_ac = 5+[0, amp_ac];
    %multi_factor = 10;
    %range_Y = 1300+[0, amp_ac*multi_factor];
    %avi_filename = neuron.show_demixed_video(save_demixed, kt, [], amp_ac, range_ac, range_Y, multi_factor);

    %% save neurons shapes
    neuron.save_neurons();

end

author: zhoupc


text:
function TwoDimEllipsoid(Location,Square_Dispersion,Scale,PlotEigVectors,PlotSquare)
% this function computes the location-dispersion ellipsoid
% see "Risk and Asset Allocation"-Springer (2005), by A. Meucci

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the ellipsoid in the r plane, solution to  ((R-Location)' * Dispersion^-1 * (R-Location) ) = Scale^2
[EigenVectors,EigenValues] = eig(Square_Dispersion);
EigenValues=diag(EigenValues);
Centered_Ellipse=[];
Angle = [0 : pi/500 : 2*pi];
NumSteps=length(Angle);
for i=1:NumSteps
    y=[cos(Angle(i))                                   % normalized variables (parametric representation of the ellipsoid)
        sin(Angle(i))];
    Centered_Ellipse=[Centered_Ellipse EigenVectors*diag(sqrt(EigenValues))*y];
end
R= Location*ones(1,NumSteps) + Scale*Centered_Ellipse;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%draw plots

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the ellipsoid
hold on
h=plot(R(1,:),R(2,:));
set(h,'color','r','linewidth',2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot a rectangle centered in Location with semisides of lengths Dispersion(1) and Dispersion(2), respectively
if PlotSquare
    Dispersion=sqrt(diag(Square_Dispersion));
    Vertex_LowRight_A=Location(1)+Scale*Dispersion(1); Vertex_LowRight_B=Location(2)-Scale*Dispersion(2);
    Vertex_LowLeft_A=Location(1)-Scale*Dispersion(1); Vertex_LowLeft_B=Location(2)-Scale*Dispersion(2);
    Vertex_UpRight_A=Location(1)+Scale*Dispersion(1); Vertex_UpRight_B=Location(2)+Scale*Dispersion(2);
    Vertex_UpLeft_A=Location(1)-Scale*Dispersion(1); Vertex_UpLeft_B=Location(2)+Scale*Dispersion(2);

    Square=[Vertex_LowRight_A            Vertex_LowRight_B
        Vertex_LowLeft_A             Vertex_LowLeft_B
        Vertex_UpLeft_A              Vertex_UpLeft_B
        Vertex_UpRight_A             Vertex_UpRight_B
        Vertex_LowRight_A            Vertex_LowRight_B];
        hold on;
        h=plot(Square(:,1),Square(:,2));
        set(h,'color','r','linewidth',2)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot eigenvectors in the r plane (centered in Location) of length the
% square root of the eigenvalues (rescaled)
if PlotEigVectors
    L_1=Scale*sqrt(EigenValues(1));
    L_2=Scale*sqrt(EigenValues(2));

    % deal with reflection: matlab chooses the wrong one
    Sign= sign(EigenVectors(1,1));
    Start_A=Location(1);                               % eigenvector 1
    End_A= Location(1) + Sign*(EigenVectors(1,1)) * L_1;
    Start_B=Location(2);
    End_B= Location(2) + Sign*(EigenVectors(1,2)) * L_1;
    hold on
    h=plot([Start_A End_A],[Start_B End_B]);
    set(h,'color','r','linewidth',2)
    axis equal;

    Start_A=Location(1);                               % eigenvector 2
    End_A= Location(1) + (EigenVectors(2,1)* L_2);
    Start_B=Location(2);
    End_B= Location(2) + (EigenVectors(2,2)* L_2);
    hold on;
    h=plot([Start_A End_A],[Start_B End_B]);
    set(h,'color','r','linewidth',2)
end

grid on

author: Sable

repo: mcbench-benchmarks


text:
% change in max anglesub
function [data,units] = compute_danglesub(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  if trx(fly).nframes <= 1 ,
    data{i} = [];
  else
    data{i} = diff(trx(fly).anglesub,1,2) ./ trx(fly).dt;
  end
end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
%                            [delta,h,alpha] = iswnbr(vSQR,thetaSQR)
% ISWNBR  Checks feasibility w.r.t. wide region/neighborhood of Sturm-Zhang.
%   vTAR:= (1-alpha)*max(h,v) projection v onto theta-central region
%   delta = (sqrt(n)/theta) * norm(vTAR - v) / norm(v)
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

function [delta,h,alpha] = iswnbr(w,thetaSQR) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

% ----------------------------------------
% r = n/thetaSQR
% hSQR = sumwNT/(r-|T|),  hubSQR = sumwNT/(r-|T| - |Q|)
% sumdifv = h*|T| - sumvT   (sumvT = sum(v_T), growing )
% sumdifw = hSQR*|T| - sumwT
% alpha = sumdifv / (r*h)
% deltaSQR = r * ( 2*alpha-alpha^2 - (1-alpha)^2 * sumdifw/gap )
% WE UPDATE sumdifv AND sumdifw IN A STABLE WAY
% ----------------------------------------
n = length(w); gap = sum(w); %#ok
sumwNT = gap;
r = n / thetaSQR;
cardT = 0; wQ = []; sumdifv = 0; sumdifw = 0;
cardQ = n;
hSQR = sumwNT / (r - cardT);  hubSQR = sumwNT / (r-(n-1));
for j = 1:n
    wj = w(j);
    if wj >= hubSQR              % wj >= hubSQR ==> not in T
        cardQ = cardQ - 1;
        hubSQR = sumwNT / (r-cardT-cardQ);
    elseif wj < hSQR             % wj < hSQR ==> in T
        cardT = cardT + 1;
        cardQ = cardQ - 1;
        hubSQR = (1-wj/sumwNT) * hubSQR;
        sumwNT = sumwNT - wj;
        oldhSQR = hSQR;
        hSQR = sumwNT / (r - cardT);
        sumdifw = sumdifw + (oldhSQR-wj) + cardT * (hSQR-oldhSQR);
        sumdifv = sumdifv + (sqrt(oldhSQR)-sqrt(wj)) + ...
            cardT * (sqrt(hSQR)-sqrt(oldhSQR));
    else                    % Inconclusive: j in Q
        wQ = [wQ;wj];
    end % if
end % for
% ----------------------------------------
% The same treatment for the Q set, but we
% sort the (presumably short) wQ first.
% ----------------------------------------
if ~isempty(wQ)
    sort(wQ);
    STOP = 0; j = 1;
    while ~STOP
        wj = wQ(j);
        if wj >= hSQR
            STOP = 1;
        else
            cardT = cardT + 1;
            sumwNT = sumwNT - wj;
            oldhSQR = hSQR;
            hSQR = sumwNT / (r - cardT);
            sumdifw = sumdifw + (oldhSQR-wj) + cardT * (hSQR-oldhSQR);
            sumdifv = sumdifv + (sqrt(oldhSQR)-sqrt(wj)) + ...
                cardT * (sqrt(hSQR)-sqrt(oldhSQR));
            j = j+1;
            if j > length(wQ)
                STOP = 1;
            end
        end
    end
end % treatment Q
% ----------------------------------------
% alpha = sumdifv/(r*h)
% deltaSQR = r * ( 2*alpha-alpha^2 - (1-alpha)^2 * sumdifw/gap )
%  (THE ABOVE DIFFERENCE SHOULD NOT BE NUMERICALLY DANGEROUS,
%    SINCE alpha IS *SIGNIF* BIGGER THAN sumdifw/gap )
% ----------------------------------------
h = sqrt(hSQR);
alpha = sumdifv/ (r*h);
deltaSQR = alpha*(2-alpha) - (1-alpha)^2 * sumdifw/gap;
delta = sqrt(r*deltaSQR);

repo: REST


text:
function z = ldivide( x, y )

%Disciplined convex programming information for LDIVIDE:
%   For DCP purposes, the LDIVIDE division operator X.\Y is equivalent to
%   (1./X).*Y. The left-hand term must be constant and non-zero; and if the
%   right-hand term is nonlinear, t constant must also be real.
%
%Disciplined geometric programming information for LDIVIDE:
%   Terms in a left divide must have opposite log-curvature, so the
%   following products are permitted:
%      {log-convex} .\ {log-concave}  {log-concave} .\ {log-convex}
%      {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.

z = times( x, y, '.\' );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function stirling2_test ( )

%*****************************************************************************80
%
%% STIRLING2_TEST tests STIRLING2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  m = 8;
  n = 8;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'STIRLING2_TEST\n' );
  fprintf ( 1, '  STIRLING2: Stirling numbers of second kind.\n' );
  fprintf ( 1, '  Get rows 1 through %d\n', m );
  fprintf ( 1, '\n' );

  s2 = stirling2 ( m, n );

  for i = 1 : m
    fprintf ( 1, '  %4d', i );
    for j = 1 : n
      fprintf ( 1, '  %6d', s2(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function[lstms,all_h_t,all_c_t]=Forward(batch,parameter,isTraining)%Forward
    N=size(batch.Word,1);
    zeroState=zeroMatrix([parameter.hidden,N]);
    if isTraining==1
        T=batch.MaxLen;
    else
        T=batch.MaxLenSource;
    end
    all_h_t=cell(parameter.layer_num,T);
    all_c_t=cell(parameter.layer_num,T);
    lstms = cell(parameter.layer_num,T);

    for ll=1:parameter.layer_num
        for tt=1:T
            all_h_t{ll,tt}=zeroMatrix([parameter.hidden,N]);
            all_c_t{ll,tt}=zeroMatrix([parameter.hidden,N]);
        end
    end
    for t=1:T
        for ll=1:parameter.layer_num
            if t<batch.MaxLenSource+1;
                W=parameter.W_S{ll};
            else
                W=parameter.W_T{ll};
            end
            if t==1
                h_t_1=zeroState;
                c_t_1 =zeroState;
            else
                c_t_1 = all_c_t{ll, t-1};
                h_t_1 = all_h_t{ll, t-1};
            end
            if ll==1
                x_t=parameter.vect(:,batch.Word(:,t));
            else
                x_t=all_h_t{ll-1,t};
            end
            x_t(:,batch.Delete{t})=0;
            h_t_1(:,batch.Delete{t})=0;
            c_t_1(:,batch.Delete{t})=0;
            [lstms{ll, t},all_h_t{ll, t},all_c_t{ll, t}]=lstmUnit(W,parameter,x_t,h_t_1,c_t_1,ll,t,isTraining);%LSTM unit calculation
        end
    end
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
% Select a region of interest for a given dataset at multiple resolution
% levels.
function imageDataset = selectROI(imageDataset, factors, ROI)

    for binningFactorIdx = 1:length(factors)

        s = factors(binningFactorIdx);

        if factors(binningFactorIdx) == 1
            % Extract the given ROI form the ground truth data.
            imageDataset.groundTruth = imageDataset.groundTruth(ROI(2):(ROI(2) + ROI(4) - 1), ROI(1):ROI(1) + ROI(3) - 1, :);
        else
            frames = imageDataset.(['Bin', num2str( factors(binningFactorIdx) )]);
            % Re-calculate ROI for the current binning factor.
            width = ceil(ROI(3) / s);
            hight = ceil(ROI(4) / s);
            x0 = ceil(ROI(1) / s);
            y0 = ceil(ROI(2) / s);
            imageDataset.(['Bin', num2str( factors(binningFactorIdx) )]) = frames(y0:(y0 + hight - 1), x0:(x0 + width - 1), :);
        end

    end

author: thomas-koehler

repo: SupER


text:
%% Copyright 2013 The MathWorks, Inc.
% Read input image
I = imread('visionteam.jpg');
figure,imshow(I);

%% Detect upright people
peopleDetector = vision.PeopleDetector;
[bboxes, scores] = step(peopleDetector,I);
I_people = insertObjectAnnotation(I,'rectangle',bboxes,scores);
figure, imshow(I_people);

%% Try the other model
peopleDetector = vision.PeopleDetector('ClassificationModel','UprightPeople_128x64');
peopleDetector.WindowStride = [4 4];
peopleDetector.MinSize = [256 128];
[bboxes, scores] = step(peopleDetector,I);
I_people = insertObjectAnnotation(I,'rectangle',bboxes,scores);
figure, imshow(I_people);

%% Use PeopleDetector with video
peopleDetector = vision.PeopleDetector;
video = vision.VideoFileReader('viptrain.avi');
viewer = vision.VideoPlayer;
while ~isDone(video)
    image = step(video);
    [bboxes, scores] = step(peopleDetector,image);
    I_people = insertObjectAnnotation(image,'rectangle',bboxes,scores);
    step(viewer,I_people);
end

%% Detect Faces in the image
% Create a detector object
faceDetector = vision.CascadeObjectDetector('FrontalFaceCART');

% Detect faces
bbox = step(faceDetector, I);

% Draw boxes around detected faces and display results
shapeInserter = vision.ShapeInserter('BorderColor','Custom','CustomBorderColor',[255 255 0]);
I_faces = step(shapeInserter, I, int32(bbox));
imshow(I_faces);

%% Use new object annotation function
I_faces = insertObjectAnnotation(I,'rectangle',bbox,[1:size(bbox,1)]);
imshow(I_faces);

%% Detect Upper Bodies in the image
% Create a detector object
bodyDetector = vision.CascadeObjectDetector('UpperBody');
bodyDetector.MinSize = [60 60];
bodyDetector.ScaleFactor = 1.05;

bbox_body = step(bodyDetector, I);

% Draw bounding boxes
I_body = insertObjectAnnotation(I,'rectangle',bbox_body,1:size(bbox_body,1));
figure, imshow(I_body);

%% Remove false detections by looking for face
bbox_face = zeros(size(bbox_body));
for i=1:length(bbox_body)
    Icrop = imcrop(I,bbox_body(i,:));
    bbox = step(faceDetector,Icrop);
    if ~isempty(bbox)
        bbox_face(i,:) = bbox + [bbox_body(i,1:2)-1 0 0];
    end
end

I_faces2 = insertObjectAnnotation(I,'rectangle',bbox_face,1:size(bbox_face,1));
figure, imshow(I_faces2);

%% Extract one face to use
Icrop = imcrop(I,bbox_body(1,:));
figure;imshow(Icrop);
bbox = step(faceDetector,Icrop);
hold on;rectangle('Position',bbox,'EdgeColor','y');

%% See how much we can rotate the image with the face still being detected
x = 10;
Irotate = imrotate(Icrop,x);
imshow(Irotate);
bbox = step(faceDetector, Irotate);
if bbox > 0
    hold on;rectangle('Position',bbox,'EdgeColor','y'); hold off;
end


text:
% SP_EVAL: Compute the value or the derivatives of a function, given by its degrees of freedom, at a given set of points.
%
%   [eu, F] = sp_eval (u, space, geometry, pts, [options]);
%   [eu, F] = sp_eval (u, space, geometry, npts, [options]);
%
% INPUT:
%
%     u:           vector of dof weights
%     space:       object defining the discrete space (see sp_scalar)
%     geometry:    geometry structure (see geo_load)
%     pts:         cell array with coordinates of points along each parametric direction
%     npts:        number of points along each parametric direction
%     options:     cell array with the fields to plot
%                   accepted options are 'value' (default), 'gradient' and 'laplacian'
%
% OUTPUT:
%
%     eu: cell-array with the fields evaluated at the given points
%     F:  grid points in the physical domain, that is, the mapped points
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2012, 2014, 2015, 2018 Rafael Vazquez
%
%    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/>.

function [eu, F] = sp_eval (u, space, geometry, npts, options)

  if (numel (u) ~= space.ndof)
    error ('The number of degrees of freedom of the vector and the space do not match')
  end

  if (nargin < 5)
    options = {'value'};
  end
  if (~iscell (options))
    options = {options};
  end
  nopts = numel (options);

  ndim = numel (space.knots);

  endpoints = zeros (2, ndim);
  if (isfield (geometry, 'nurbs'))
    nurbs = geometry.nurbs;
    if (ndim == 1)
      nurbs.knots = {nurbs.knots};
    end
    for idim=1:ndim
      endpoints(:,idim) = nurbs.knots{idim}([nurbs.order(idim), end-nurbs.order(idim)+1]);
    end
    clear nurbs
  elseif (isfield (struct(space), 'knots'))
    degree = space.degree;
    for idim=1:ndim
      endpoints(:,idim) = space.knots{idim}([degree(idim)+1, end-degree(idim)]);
    end
  else
    endpoints(2,:) = 1;
  end

  if (iscell (npts))
    pts = npts;
    npts = cellfun (@numel, pts);
  elseif (isvector (npts))
    if (numel (npts) == 1)
      npts = npts * ones (1,ndim);
    end
    for idim = 1:ndim
      pts{idim} = linspace (endpoints(1,idim), endpoints(2,idim), npts(idim));
    end
  end

  for jj = 1:ndim
    pts{jj} = pts{jj}(:)';
    if (numel (pts{jj}) > 1)
      brk{jj} = [endpoints(1,jj), pts{jj}(1:end-1) + diff(pts{jj})/2, endpoints(2,jj)];
    else
      brk{jj} = endpoints(:,jj).';
    end
  end

  msh = msh_cartesian (brk, pts, [], geometry, 'boundary', false);
  sp  = space.constructor (msh);

  value = false; grad = false; laplacian = false; hessian = false;

  for iopt = 1:nopts
    switch (lower (options{iopt}))
      case 'value'
        eu{iopt} = zeros (msh.nqn, msh.nel);
        eunum{iopt} = {1:msh.nqn};
        eusize{iopt} = npts;
        value = true;

      case 'gradient'
        eu{iopt} = zeros (msh.rdim, msh.nqn, msh.nel);
        eunum{iopt} = {1:msh.rdim, 1:msh.nqn};
        eusize{iopt} = [msh.rdim, npts];
        grad = true;

      case 'laplacian'
        eu{iopt} = zeros (msh.nqn, msh.nel);
        eunum{iopt} = {1:msh.nqn};
        eusize{iopt} = npts;
        laplacian = true;

      case 'hessian'
        eu{iopt} = zeros (msh.rdim, msh.rdim, msh.nqn, msh.nel);
        eunum{iopt} = {1:msh.rdim, 1:msh.rdim, 1:msh.nqn};
        eusize{iopt} = [msh.rdim, msh.rdim, npts];
        hessian = true;
    end
  end

  F = zeros (msh.rdim, msh.nqn, msh.nel);

  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    sp_col  = sp_evaluate_col (sp, msh_col, 'value', value, 'gradient', grad, ...
          'laplacian', laplacian, 'hessian', hessian);

    eu_aux = sp_eval_msh (u, sp_col, msh_col, options);

    F(:,:,msh_col.elem_list) = msh_col.geo_map;
    for iopt = 1:nopts
      eu{iopt}(eunum{iopt}{:},msh_col.elem_list) = eu_aux{iopt};
    end
  end

  F = reshape (F, [msh.rdim, npts]);
  for iopt = 1:nopts
    eu{iopt} = reshape (eu{iopt}, [eusize{iopt}, 1]); % The extra 1 makes things work also in 1D
  end

  if (nopts == 1)
    eu = eu{1};
  end
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function con = spm_cfg_con
% SPM Configuration file for contrast specification
%__________________________________________________________________________
% Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_con.m 6952 2016-11-25 16:03:13Z guillaume $

%--------------------------------------------------------------------------
% spmmat Select SPM.mat
%--------------------------------------------------------------------------
spmmat         = cfg_files;
spmmat.tag     = 'spmmat';
spmmat.name    = 'Select SPM.mat';
spmmat.help    = {'Select SPM.mat file for contrasts.'};
spmmat.filter  = 'mat';
spmmat.ufilter = '^SPM\.mat$';
spmmat.num     = [1 1];

%--------------------------------------------------------------------------
% name Name
%--------------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Name of contrast.'};
name.strtype = 's';
name.num     = [1 Inf];

%--------------------------------------------------------------------------
% convec T-contrast weights vector
%--------------------------------------------------------------------------
convec         = cfg_entry;
convec.tag     = 'weights';
convec.name    = 'Weights vector';
convec.help    = {
    'Enter T-contrast weights vector.'
    'This is done similarly to the contrast manager. A 1 x n vector should be entered for T-contrasts.'
    'Contrast weight vectors will be padded with zeros to the correct length.'
    };
convec.strtype = 'r';
convec.num     = [1 Inf];

%--------------------------------------------------------------------------
% sessrep Replicate over sessions
%--------------------------------------------------------------------------
sessrep         = cfg_menu;
sessrep.tag     = 'sessrep';
sessrep.name    = 'Replicate over sessions';
sessrep.val = {'none'};
sessrep.help    = {
    'If there are multiple sessions with identical conditions, one might want to specify contrasts which are identical over sessions. This can be done automatically based on the contrast spec for one session.'
    'Contrasts can be either replicated (thus testing average effects over sessions) or created per session. In both cases, zero padding up to the length of each session and the block effects is done automatically. In addition, weights of replicated contrasts can be scaled by the number of sessions. This allows to use the same contrast manager batch for fMRI analyses with a variable number of sessions. The scaled contrasts can then be compared in a 2nd level model without a need for further adjustment of effect sizes.'
    };
sessrep.labels = {
    'Don''t replicate'
    'Replicate'
    'Replicate&Scale'
    'Create per session'
    'Both: Replicate + Create per session'
    'Both: Replicate&Scale + Create per session'
    }';
sessrep.values = {
    'none'
    'repl'
    'replsc'
    'sess'
    'both'
    'bothsc'
    }';

%==========================================================================
% tcon T-contrast
%==========================================================================
tcon      = cfg_branch;
tcon.tag  = 'tcon';
tcon.name = 'T-contrast';
tcon.val  = {name convec sessrep};
tcon.help = {
    '* Simple one-dimensional contrasts for an SPM{T}'
    ''
    'A simple contrast for an SPM{T} tests the null hypothesis c''B=0 against the one-sided alternative c''B>0, where c is a column vector. '
    ''
    '    Note that throughout SPM, the transpose of the contrast weights is used for display and input. That is, you''ll enter and visualise c''. For an SPM{T} this will be a row vector.'
    ''
    'For example, if you have a design in which the first two columns of the design matrix correspond to the effects for "baseline" and "active" conditions respectively, then a contrast with weights c''=[-1,+1,0,...] (with zero weights for any other parameters) tests the hypothesis that there is no "activation" (the parameters for both conditions are the same), against the alternative that there is some activation (i.e. the parameter for the "active" condition is greater than that for the "baseline" condition). The resulting SPM{T} (created by spm_getSPM.m) is a statistic image, with voxel values the value of the t-statistic for the specified contrast at that location. Areas of the SPM{T} with high voxel values indicate evidence for "activation". To look for areas of relative "de-activation", the inverse contrast could be used c''=[+1,-1,0,...].'
    ''
    'Similarly, if you have a design where the third column in the design matrix is a covariate, then the corresponding parameter is essentially a regression slope, and a contrast with weights c''=[0,0,1,0,...] (with zero weights for all parameters but the third) tests the hypothesis of zero regression slope, against the alternative of a positive slope. This is equivalent to a test no correlation, against the alternative of positive correlation. If there are other terms in the model beyond a constant term and the covariate, then this correlation is apartial correlation, the correlation between the data Y and the covariate, after accounting for the other effects.'
    }';

%--------------------------------------------------------------------------
% name Name
%--------------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Name of contrast.'};
name.strtype = 's';
name.num     = [1 Inf];

%--------------------------------------------------------------------------
% convec F-contrast weights matrix
%--------------------------------------------------------------------------
convec         = cfg_entry;
convec.tag     = 'weights';
convec.name    = 'Weights matrix';
convec.help    = {
    'Enter F-contrast weights matrix.'
    'This is done similarly to the contrast manager.'
    'Contrast weight matrices will be padded with zeros to the correct length.'
    };
convec.strtype = 'r';
convec.num     = [Inf Inf];

%==========================================================================
% fcon F-contrast
%==========================================================================
fcon      = cfg_branch;
fcon.tag  = 'fcon';
fcon.name = 'F-contrast';
fcon.val  = {name convec sessrep};
fcon.help = {
    '* Linear constraining matrices for an SPM{F}'
    ''
    'The null hypothesis c''B=0 can be thought of as a (linear) constraint on the full model under consideration, yielding a reduced model. Taken from the viewpoint of two designs, with the full model an extension of the reduced model, the null hypothesis is that the additional terms in the full model are redundent.'
    ''
    'Statistical inference proceeds by comparing the additional variance explained by full design over and above the reduced design to the error variance (of the full design), an "Extra Sum-of-Squares" approach yielding an F-statistic for each voxel, whence an SPM{F}.'
    ''
    'This is useful in a number of situations:'
    ''
    '* Two sided tests'
    ''
    'The simplest use of F-contrasts is to effect a two-sided test of a simple linear contrast c''B, where c is a column vector. The SPM{F} is the square of the corresponding SPM{T}. High values of the SPM{F} therefore indicate evidence against the null hypothesis c''B=0 in favour of the two-sided alternative c''B~=0.'
    ''
    '* General linear hypotheses'
    ''
    'Where the contrast weights is a matrix, the rows of the (transposed) contrast weights matrix c'' must define contrasts in their own right, and the test is effectively simultaneously testing the null hypotheses associated with the individual component contrasts with weights defined in the rows. The null hypothesis is still c''B=0, but since c is a matrix, 0 here is a zero vector rather than a scalar zero, asserting that under the null hypothesis all the component hypotheses are true.'
    ''
    'For example: Suppose you have a language study with 3 word categories (A,B & C), and would like to test whether there is any difference at all between the three levels of the "word category" factor.'
    ''
    'The design matrix might look something like:'
    ''
    '         [ 1 0 0 ..]'
    '         [ : : : ..]'
    '         [ 1 0 0 ..]'
    '         [ 0 1 0 ..]'
    '    X =  [ : : : ..]'
    '         [ 0 1 0 ..]'
    '         [ 0 0 1 ..]'
    '         [ : : : ..]'
    '         [ 0 0 1 ..]'
    '         [ 0 0 0 ..]'
    '         [ : : : ..]'
    ''
    ' ...with the three levels of the "word category" factor modelled in the  first three columns of the design matrix.'
    ''
    'The matrix of contrast weights will look like:'
    ''
    ' c'' = [1 -1  0 ...;'
    '       0  1 -1 ...]'
    ''
    'Reading the contrasts weights in each row of c'', we see that row 1 states that category A elicits the same response as category B, row 2 that category B elicits the same response as category C, and hence together than categories A, B & C all elicit the same response.'
    ''
    'The alternative hypothesis is simply that the three levels are not all the same, i.e. that there is some difference in the paraeters for the three levels of the factor: The first and the second categories produce different brain responses, OR the second and third categories, or both.'
    ''
    'In other words, under the null hypothesis (the categories produce the same brain responses), the model reduces to one in which the three level "word category" factor can be replaced by a single "word" effect, since there is no difference in the parameters for each category. The corresponding design matrix would have the first three columns replaced by a single column that is the sum (across rows) of the first three columns in the design matric above, modelling the brain response to a word, whatever is the category. The F-contrast above is in fact testing the hypothesis that this reduced design doesn''t account for significantly less variance than the full design with an effect for each word category.'
    ''
    'Another way of seeing that, is to consider a reparameterisation of the model, where the first column models effects common to all three categories, with the second and third columns modelling the differences between the three conditions, for example:'
    ''
    '         [ 1  1  0 ..]'
    '         [ :  :  : ..]'
    '         [ 1  1  0 ..]'
    '         [ 1  0  1 ..]'
    '    X =  [ :  :  : ..]'
    '         [ 1  0  1 ..]'
    '         [ 1 -1 -1 ..]'
    '         [ :  :  : ..]'
    '         [ 1 -1 -1 ..]'
    '         [ 0  0  0 ..]'
    '         [ :  :  : ..]'
    ''
    'In this case, an equivalent F contrast is of the form'
    ' c'' = [ 0 1 0 ...;'
    '        0 0 1 ...]'
    'and would be exactly equivalent to the previous contrast applied to the previous design. In this latter formulation, you are asking whewher the two columns modelling the "interaction space" account for a significant amount of variation (variance) of the data. Here the component contrasts in the rows of c'' are simply specifying that the parameters for the corresponding rows are are zero, and it is clear that the F-test is comparing this full model with a reduced model in which the second and third columns of X are omitted.'
    ''
    '    Note the difference between the following two F-contrasts:'
    '         c'' = [ 0 1 0 ...;     (1)'
    '                0 0 1 ...]'
    '     and'
    '         c'' = [ 0 1 1 ...]     (2)'
    ''
    '    The first is an F-contrast, testing whether either of the parameters for the effects modelled in the 2nd & 3rd columns of the design matrix are significantly different from zero. Under the null hypothesis c''B=0, the first contrast imposes a two-dimensional constraint on the design. The second contrast tests whether the SUM of the parameters for the 2nd & 3rd columns is significantly different from zero. Under the null hypothesis c''B=0, this second contrast only imposes a one dimensional constraint on the design.'
    ''
    '    An example of the difference between the two is that the first contrast would be sensitive to the situation where the 2nd & 3rd parameters were +a and -a, for some constant a, wheras the second contrast would not detect this, since the parameters sum to zero.'
    ''
    'The test for an effect of the factor "word category" is an F-test with 3-1=2 "dimensions", or degrees of freedom.'
    ''
    '* Testing the significance of effects modelled by multiple columns'
    ''
    'A conceptially similar situation arises when one wonders whether a set of coufound effects are explaining any variance in the data. One important advantage of testing the with F contrasts rather than one by one using SPM{T}''s is the following. Say you have two covariates that you would like to know whether they can "predict" the brain responses, and these two are correlated (even a small correlation would be important in this instance). Testing one and then the other may lead you to conclude that there is no effect. However, testing with an F test the two covariates may very well show a not suspected effect. This is because by testing one covariate after the other, one never tests for what is COMMON to these covariates (see Andrade et al, Ambiguous results in functional neuroimaging, NeuroImage, 1999).'
    ''
    ''
    'More generally, F-tests reflect the usual analysis of variance, while t-tests are traditionally post hoc tests, useful to see in which direction is an effect going (positive or negative). The introduction of F-tests can also be viewed as a first means to do model selection.'
    ''
    ''
    'Technically speaking, an F-contrast defines a number of directions (as many as the rank of the contrast) in the space spanned by the column vectors of the design matrix. These directions are simply given by X*c if the vectors of X are orthogonal, if not, the space define by c is a bit more complex and takes care of the correlation within the design matrix. In essence, an F-contrast is defining a reduced model by imposing some linear constraints (that have to be estimable, see below) on the parameters estimates. Sometimes, this reduced model is simply made of a subset of the column of the original design matrix but generally, it is defined by a combination of those columns. (see spm_FcUtil for what (I hope) is an efficient handling of F-contrats computation).'
    }';

%--------------------------------------------------------------------------
% name Name
%--------------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Name of contrast.'};
name.strtype = 's';
name.num     = [1 Inf];

%--------------------------------------------------------------------------
% conweight Contrast weight
%--------------------------------------------------------------------------
conweight         = cfg_entry;
conweight.tag     = 'conweight';
conweight.name    = 'Contrast weight';
conweight.help    = {'The contrast weight for the selected column.'};
conweight.strtype = 'r';
conweight.num     = [1 1];

%--------------------------------------------------------------------------
% colcond Condition #
%--------------------------------------------------------------------------
colcond         = cfg_entry;
colcond.tag     = 'colcond';
colcond.name    = 'Condition #';
colcond.help    = {'Select which condition function set is to be contrasted.'};
colcond.strtype = 'n';
colcond.num     = [1 1];

%--------------------------------------------------------------------------
% colbf Basis function #
%--------------------------------------------------------------------------
colbf         = cfg_entry;
colbf.tag     = 'colbf';
colbf.name    = 'Basis function #';
colbf.help    = {'Select which basis function from the basis function set is to be contrasted.'};
colbf.strtype = 'n';
colbf.num     = [1 1];

%--------------------------------------------------------------------------
% colmod Parametric modulation #
%--------------------------------------------------------------------------
colmod         = cfg_entry;
colmod.tag     = 'colmod';
colmod.name    = 'Parametric modulation #';
colmod.help    = {'Select which parametric modulation is to be contrasted. If there is no time/parametric modulation, enter "1". If there are both time and parametric modulations, then time modulation comes before parametric modulation.'};
colmod.strtype = 'n';
colmod.num     = [1 1];

%--------------------------------------------------------------------------
% colmodord Parametric modulation order
%--------------------------------------------------------------------------
colmodord         = cfg_entry;
colmodord.tag     = 'colmodord';
colmodord.name    = 'Parametric modulation order';
colmodord.help    = {
    'Order of parametric modulation to be contrasted. '
    ''
    '0 - the basis function itself, 1 - 1st order mod etc'
    }';
colmodord.strtype = 'w';
colmodord.num     = [1 1];

%--------------------------------------------------------------------------
% colconds Contrast entry
%--------------------------------------------------------------------------
colconds      = cfg_branch;
colconds.tag  = 'colconds';
colconds.name = 'Contrast entry';
colconds.val  = {conweight colcond colbf colmod colmodord};
colconds.help = {''};

%--------------------------------------------------------------------------
% generic T contrast for conditions
%--------------------------------------------------------------------------
generic        = cfg_repeat;
generic.tag    = 'generic';
generic.name   = 'T contrast for conditions';
generic.help   = {'Assemble your contrast column by column.'};
generic.values = {colconds};
generic.val    = {colconds};
generic.num    = [1 Inf];

%--------------------------------------------------------------------------
% colreg T contrast for extra regressors
%--------------------------------------------------------------------------
colreg         = cfg_entry;
colreg.tag     = 'colreg';
colreg.name    = 'T contrast for extra regressors';
colreg.help    = {'Enter T contrast vector for extra regressors.'};
colreg.strtype = 'r';
colreg.num     = [1 Inf];

%--------------------------------------------------------------------------
% coltype Contrast columns
%--------------------------------------------------------------------------
coltype        = cfg_choice;
coltype.tag    = 'coltype';
coltype.name   = 'Contrast columns';
coltype.val    = {generic };
coltype.help   = {'Contrasts can be specified either over conditions or over extra regressors.'};
coltype.values = {generic colreg };

%--------------------------------------------------------------------------
% sessions Session(s)
%--------------------------------------------------------------------------
sessions         = cfg_entry;
sessions.tag     = 'sessions';
sessions.name    = 'Session(s)';
sessions.help    = {'Enter session number(s) for which this contrast should be created. If more than one session number is specified, the contrast will be an average contrast over the specified conditions or regressors from these sessions.'};
sessions.strtype = 'n';
sessions.num     = [1 Inf];

%==========================================================================
% tconsess T-contrast (cond/sess based)
%==========================================================================
tconsess         = cfg_branch;
tconsess.tag     = 'tconsess';
tconsess.name    = 'T-contrast (cond/sess based)';
tconsess.val     = {name coltype sessions};
tconsess.help    = {
    'Define a contrast in terms of conditions or regressors instead of columns of the design matrix. This allows to create contrasts automatically even if some columns are not always present (e.g. parametric modulations).'
    ''
    'Each contrast column can be addressed by specifying'
    '* session number'
    '* condition number'
    '* basis function number'
    '* parametric modulation number and'
    '* parametric modulation order.'
    ''
    'If the design is specified without time or parametric modulation, SPM creates a "pseudo-modulation" with order zero. To put a contrast weight on a basis function one therefore has to enter "1" for parametric modulation number and "0" for parametric modulation order.'
    ''
    'Time and parametric modulations are not distinguished internally. If time modulation is present, it will be parametric modulation "1", and additional parametric modulations will be numbered starting with "2".'
    ''
    '* Simple one-dimensional contrasts for an SPM{T}'
    ''
    'A simple contrast for an SPM{T} tests the null hypothesis c''B=0 against the one-sided alternative c''B>0, where c is a column vector. '
    ''
    '    Note that throughout SPM, the transpose of the contrast weights is used for display and input. That is, you''ll enter and visualise c''. For an SPM{T} this will be a row vector.'
    ''
    'For example, if you have a design in which the first two columns of the design matrix correspond to the effects for "baseline" and "active" conditions respectively, then a contrast with weights c''=[-1,+1,0,...] (with zero weights for any other parameters) tests the hypothesis that there is no "activation" (the parameters for both conditions are the same), against the alternative that there is some activation (i.e. the parameter for the "active" condition is greater than that for the "baseline" condition). The resulting SPM{T} (created by spm_getSPM.m) is a statistic image, with voxel values the value of the t-statistic for the specified contrast at that location. Areas of the SPM{T} with high voxel values indicate evidence for "activation". To look for areas of relative "de-activation", the inverse contrast could be used c''=[+1,-1,0,...].'
    ''
    'Similarly, if you have a design where the third column in the design matrix is a covariate, then the corresponding parameter is essentially a regression slope, and a contrast with weights c''=[0,0,1,0,...] (with zero weights for all parameters but the third) tests the hypothesis of zero regression slope, against the alternative of a positive slope. This is equivalent to a test no correlation, against the alternative of positive correlation. If there are other terms in the model beyond a constant term and the covariate, then this correlation is apartial correlation, the correlation between the data Y and the covariate, after accounting for the other effects.'
    }';

%--------------------------------------------------------------------------
% consess Contrast Sessions
%--------------------------------------------------------------------------
consess        = cfg_repeat;
consess.tag    = 'consess';
consess.name   = 'Contrast Sessions';
consess.help   = {
    'For general linear model Y = XB + E with data Y, desgin matrix X, parameter vector B, and (independent) errors E, a contrast is a linear combination of the parameters c''B. Usually c is a column vector, defining a simple contrast of the parameters, assessed via an SPM{T}. More generally, c can be a matrix (a linear constraining matrix), defining an "F-contrast" assessed via an SPM{F}.'
    ''
    'The vector/matrix c contains the contrast weights. It is this contrast weights vector/matrix that must be specified to define the contrast. The null hypothesis is that the linear combination c''B is zero. The order of the parameters in the parameter (column) vector B, and hence the order to which parameters are referenced in the contrast weights vector c, is determined by the construction of the design matrix.'
    ''
    'There are two types of contrast in SPM: simple contrasts for SPM{T}, and "F-contrasts" for SPM{F}.'
    ''
    'For a thorough theoretical treatment, see the Human Brain Function book and the statistical literature referenced therein.'
    ''
    ''
    '* Non-orthogonal designs'
    ''
    'Note that parameters zero-weighted in the contrast are still included in the model. This is particularly important if the design is not orthogonal (i.e. the columns of the design matrix are not orthogonal). In effect, the significance of the contrast is assessed *after* accounting for the other effects in the design matrix. Thus, if two covariates are correlated, testing the significance of the parameter associated with one will only test for the part that is not present in the second covariate. This is a general point that is also true for F-contrasts. See Andrade et al, Ambiguous results in functional neuroimaging, NeuroImage, 1999, for a full description of the effect of non othogonal design testing.'
    ''
    ''
    '* Estimability'
    ''
    'The contrast c''B is estimated by c''b, where b are the parameter estimates given by b=pinv(X)*Y.'
    ''
    'However, if a design is rank-deficient (i.e. the columns of the design matrix are not linearly independent), then the parameters are not unique, and not all linear combinations of the parameter are valid contrasts, since contrasts must be uniquely estimable.'
    ''
    'A weights vector defines a valid contrast if and only if it can be constructed as a linear combination of the rows of the design matrix. That is c'' (the transposed contrast vector - a row vector) is in the row-space of the design matrix.'
    ''
    'Usually, a valid contrast will have weights that sum to zero over the levels of a factor (such as condition).'
    ''
    'A simple example is a simple two condition design including a constant, with design matrix'
    ''
    '          [ 1 0 1 ]'
    '          [ : : : ]'
    '     X =  [ 1 0 1 ]'
    '          [ 0 1 1 ]'
    '          [ : : : ]'
    '          [ 0 1 1 ]'
    ''
    'The first column corresponds to condition 1, the second to condition 2, and the third to a constant (mean) term. Although there are three columns to the design matrix, the design only has two degrees of freedom, since any one column can be derived from the other two (for instance, the third column is the sum of the first two). There is no unique set of parameters for this model, since for any set of parameters adding a constant to the two condition effects and subtracting it from the constant effect yields another set of viable parameters. However, the difference between the two condition effects is uniquely estimated, so c''=[-1,+1,0] does define a contrast.'
    ''
    'If a parameter is estimable, then the weights vector with a single "1" corresponding to that parameter (and zero elsewhere) defines a valid contrast.'
    ''
    ''
    '* Multiple comparisons'
    ''
    'Note that SPM implements no corrections to account for you looking at multiple contrasts.'
    ''
    'If you are interested in a set of hypotheses that together define a consistent question, then you should account for this when assessing the individual contrasts. A simple Bonferroni approach would assess N simultaneous contrasts at significance level alpha/N, where alpha is the chosen significance level (usually 0.05).'
    ''
    'For two sided t-tests using SPM{T}s, the significance level should be halved. When considering both SPM{T}s produced by a contrast and it''s inverse (the contrast with negative weights), to effect a two-sided test to look for both "increases" and "decreases", you should review each SPM{T} at at level 0.05/2 rather than 0.05. (Or consider an F-contrast!)'
    ''
    ''
    '* Contrast images and ESS images'
    ''
    'For a simple contrast, SPM (spm_getSPM.m) writes a contrast image: con_????.{img,nii}, with voxel values c''b. (The ???? in the image names are replaced with the contrast number.) These contrast images (for appropriate contrasts) are suitable summary images of an effect at this level, and can be used as input at a higher level when effecting a random effects analysis.'
    ''
    'For an F-contrast, SPM (spm_getSPM.m) writes the Extra Sum-of-Squares (the difference in the residual sums of squares for the full and reduced model) as ess_????.{img,nii}. (Note that the ess_????.{img,nii} and SPM{T,F}_????.{img,nii} images are not suitable input for a higher level analysis.)'
    }';
consess.values = {tcon fcon tconsess};
consess.num    = [0 Inf];

%--------------------------------------------------------------------------
% delete Delete existing contrasts
%--------------------------------------------------------------------------
delete        = cfg_menu;
delete.tag    = 'delete';
delete.name   = 'Delete existing contrasts';
delete.help   = {'Delete existing contrasts.'};
delete.labels = {'Yes', 'No'};
delete.values = {1 0};
delete.val    = {0};

%--------------------------------------------------------------------------
% con Contrast Manager
%--------------------------------------------------------------------------
con      = cfg_exbranch;
con.tag  = 'con';
con.name = 'Contrast Manager';
con.val  = {spmmat consess delete};
con.help = {'Specify T and F contrasts.'};
con.prog = @spm_run_con;
con.vout = @vout_stats;

%==========================================================================
function dep = vout_stats(job)
dep(1)            = cfg_dep;
dep(1).sname      = 'SPM.mat File';
dep(1).src_output = substruct('.','spmmat');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
%dep(2)            = cfg_dep;
%dep(2).sname      = 'SPM Variable';
%dep(2).src_output = substruct('.','spmvar');
%dep(2).tgt_spec   = cfg_findspec({{'strtype','e'}});
dep(2)            = cfg_dep;
dep(2).sname      = 'All Con Images';
dep(2).src_output = substruct('.','con');
dep(2).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
dep(3)            = cfg_dep;
dep(3).sname      = 'All Stats Images';
dep(3).src_output = substruct('.','spm');
dep(3).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});

author: spm


text:
%LOGLC Trainable logistic linear classifier
%
%   W = LOGLC(A)
%   W = A*LOGLC
%
% INPUT
%   A   Dataset
%
% OUTPUT
%   W   Logistic linear classifier
%
% DESCRIPTION
% Computation of the linear classifier for the dataset A by maximizing the
% likelihood criterion using the logistic (sigmoid) function.
% This routine becomes very slow for feature sizes above 1000.
%
% REFERENCES
% A. Webb, Statistical Pattern Recognition, John Wiley & Sons, New York, 2002.
% J. A. Anderson, Logistic discrimination, in: P. R. Krishnaiah and L. N.
% Kanal (eds.), Handbook of Statistics 2: Classification, Pattern Recognition
% and Reduction of Dimensionality, North Holland, Amsterdam, 1982, 169--191.
%
%  SEE ALSO (PRTools Guide)
%  MAPPINGS, DATASETS, LDC, FISHERC

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: loglc.m,v 1.8 2010/02/08 15:31:48 duin Exp $

function W = loglc(a)

		% No input data, return an untrained classifier.
	if (nargin == 0) | (isempty(a))
		W = prmapping(mfilename);
		W = setname(W,'Logistic');
		return
	end

	islabtype(a,'crisp');
	isvaldfile(a,1,2); % at least 2 object per class, 2 classes
	a = testdatasize(a);
	fid = []; % Progress messages default destination
	[m,k,c] = getsize(a);
	nlab = getnlab(a);
	prior = getprior(a);
	a = setprior(a,prior);
	if (c > 2)
		% Compute C classifiers: each class against all others.
		W = mclassc(a,prmapping(mfilename));
	else
		v = scalem(a,'variance');
		a = a*v;
		accuracy = 0.0001;		% An accuracy for the ML loop.
		x = [a,ones(m,1)];
		% A standard trick to set the labels to +1 for the first class
		% (NLAB=1) and to -1 for the second one (NLAB=2). Then, each
		% object vector is multiplied by its new label +/-1.
		x(find(nlab==2),:) = -x(find(nlab==2),:);
		x = +x;
		alf = sum(nlab==2)/sum(nlab==1);
		weights = zeros(1,k+1);

		% Maximize the likelihood L to find WEIGHTS
		L = -inf; Lnew = -realmax;
		prwaitbar(100,'loglc: Optimizing log likelihoods',k > 100)
		d0 = 10*log(10);
		while (abs(Lnew - L) > accuracy)
			prwaitbar(100, 100-100*(log(abs(Lnew-L))-log(accuracy))/d0);
			pax = ones(m,1) ./ (1 + exp(-x*weights'));	% Estimate of P(class +1|x).
			pbx = 1 - pax;                          		% Estimate of P(class -1|x).
			L = Lnew; Lnew = sum(log(pax+realmin));     % Update likelihood.
			p2x = sqrt(pax.*pbx);
			y = x .* p2x(:,ones(1,k+1));
			%size(y'*y)
			weights = pbx' * x * prpinv(y'*y) + weights;
		end
		prwaitbar(0);

		% Define LOGLC by an affine (linear) mapping based
		% on the [K x 1] matrix R and the offset w0.
		w0 = weights(k+1) + log(alf*prior(1)/prior(2));
		R  = weights(1:k)';
		%DXD: for a two-class classifier we have to supply two posterior
		%probabilities:
		%W  = v*affine(R,w0,a,getlablist(a));  % wrong
		W  = v*affine([R -R],[w0 -w0],a,getlablist(a));
		W  = setout_conv(W,1);
	end
	W = setname(W,'Logistic');

	return

	function [weights,L,Lnew,len2] = optimw(x,weights,L,Lnew,accuracy,len2,fid)
		% this function is never called
		[m,k] = size(x);
		while (abs(Lnew - L) > accuracy)
			pax = ones(m,1) ./ (1 + exp(-x*weights'));	% Estimate of P(class +1|x).
			pbx = 1 - pax;                          		% Estimate of P(class -1|x).
			L = Lnew; Lnew = sum(log(pax+realmin));     % Update likelihood.
			p2x = sqrt(pax.*pbx);
			y = x .* p2x(:,ones(1,k));
			weights = pbx' * x * prpinv(y'*y) + weights;
		end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
I = im2double(imread('fruit.bmp'));
I_d = Yang2010(I);
figure(1)
subplot(1, 3, 1), imshow(I), xlabel('(a)', 'FontName', 'Helvetica', 'FontSize', 8, 'VerticalAlignment', 'bottom')
subplot(1, 3, 2), imshow(I_d), xlabel('(b)', 'FontName', 'Helvetica', 'FontSize', 8, 'VerticalAlignment', 'bottom')
subplot(1, 3, 3), imshow(I-I_d), xlabel('(c)', 'FontName', 'Helvetica', 'FontSize', 8, 'VerticalAlignment', 'bottom')
set(gcf, 'PaperUnits', 'inches')
set(gcf, 'PaperPosition', [0, 0, 5.35, 3])
print('example.png', '-dpng', '-r300')
% convert figure1.png -trim figure1.png

author: vitorsr

repo: SIHR

save_path: github-repos/MATLAB/vitorsr-SIHR


text:
function varargout = spy(X,variablesonly)
%SPY (overloaded)

if nargin < 2
    variablesonly = 0;
end
 if isa(X,'blkvar')
    X = sdpvar(X);
 end

if variablesonly
    Z = X.basis(:,2:end);
else
    Z = X.basis;
end

Z = reshape(sum(abs(Z),2),X.dim(1),X.dim(2));
Z = Z~=0;
if nargout==0
    spy(Z)
else
    varargout{1}=Z;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [x,value] = steepestDescent(fun,x,varargin)
% find maximum with steepest descent
%
% Input
%  fun - @S2Fun, @SO3Fun
%  x   - initial seed @vector3d, @rotation
%
% Output
%  x     - local maxima
%  value - function values
%
% Options
%  maxIter    -
%  resolution -

% prepare stepsize computation by Armijo search
maxIter = get_option(varargin,'maxIter',30);
res = get_option(varargin,'resolution',0.05*degree);
omega = 1.25.^(-30:1:10) * degree;
omega(omega<res) = [];
omega = [0,omega];

% use fixed number of iterations
for k = 1:maxIter

  % compute gradient - element of tangential space
  g = normalize(fun.grad(x));

  % prepare for linesearch
  line_x = exp(repmat(x(:),1,length(omega)),g(:) * omega);

  % evaluate along lines
  line_v = fun.eval(line_x);

  % find maxima along line
  [value,id] = max(line_v,[],2);

  % update nodes
  x = line_x(sub2ind(size(line_x),(1:length(x)).',id));

  % if there was no change -> abbort
  if all(id == 1), break; end
  %fprintf('.')
end

% [o2,v1,v2] = unique(ori)
% v = accumarray(v2,1)
%id = v>5;
%o2 = o2(id);
%v = v(id)

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function opennii(nii)
imtool3D_nii_3planes(nii)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% EX_STOKES_SQUARE_TH: solve the Stokes problem in the unit square with generalized Taylor-Hood elements.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_square.txt';

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:4;
problem_data.nmnn_sides = [];

% Physical parameters
problem_data.viscosity = @(x, y) ones (size (x));

% Force term
fx = @(x, y) (6 * x + y .* cos(x .* y) + 2 * cos(y) .* sin(x));
fy = @(x, y) (x .* cos(x .* y) - 2 * cos(x) .* sin(y));
problem_data.f  = @(x, y) cat(1, ...
                reshape (fx (x,y), [1, size(x)]), ...
                reshape (fy (x,y), [1, size(x)]));

% Boundary terms
uxex = @(x, y) (sin(x) .* cos(y));
uyex = @(x, y) (-sin(y) .* cos(x));

problem_data.h = @(x, y, iside) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

% Exact solution, to compute the errors
problem_data.velex = @(x, y) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

problem_data.pressex = @(x, y) (3 * x.^2 + sin(x .* y) - 1.239811742000564725943866);

problem_data.gradvelex = @test_stokes_square_bc_graduex;

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.element_name = 'th';   % Element type for discretization
method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
method_data.nsub         = [10 10];  % Number of subdivisions
method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = ...
                       solve_stokes (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) COMPARISON WITH EXACT SOLUTION
error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex)
[error_h1_v, error_l2_v] = ...
   sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex)

% 4.2) EXPORT TO PARAVIEW
output_file = 'SQUARE_TH_Deg3_Reg2_Sub10';

fprintf ('The result is saved in the files %s \n and %s \n \n', ...
           [output_file '_vel'], [output_file '_press']);
vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
sp_to_vtk (press, space_p, geometry, vtk_pts, [output_file '_press'], 'press')
sp_to_vtk (vel,   space_v, geometry, vtk_pts, [output_file '_vel'  ], {'velocity', 'divergence'}, {'value', 'divergence'})

% 4.3) PLOT IN MATLAB
[eu, F] = sp_eval (vel, space_v, geometry, vtk_pts);
[X,  Y] = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));

figure()
subplot (1,2,2)
eu2 = problem_data.velex (X, Y);
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
axis equal
title('Exact solution')
subplot (1,2,1)
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
axis equal
title('Computed solution')

%!demo
%! ex_stokes_square_th

%!test
%! problem_data.geo_name = 'geo_square.txt';
%! problem_data.drchlt_sides = 1:4;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y) ones (size (x));
%! fx = @(x, y) (6 * x + y .* cos(x .* y) + 2 * cos(y) .* sin(x));
%! fy = @(x, y) (x .* cos(x .* y) - 2 * cos(x) .* sin(y));
%! problem_data.f  = @(x, y) cat(1, ...
%!                 reshape (fx (x,y), [1, size(x)]), ...
%!                 reshape (fy (x,y), [1, size(x)]));
%! uxex = @(x, y) (sin(x) .* cos(y));
%! uyex = @(x, y) (-sin(y) .* cos(x));
%! problem_data.h = @(x, y, iside) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.velex = @(x, y) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.pressex = @(x, y) (3 * x.^2 + sin(x .* y) - 1.239811742000564725943866);
%! problem_data.gradvelex = @test_stokes_square_bc_graduex;
%! method_data.element_name = 'th';   % Element type for discretization
%! method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
%! method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
%! method_data.nsub         = [10 10];  % Number of subdivisions
%! method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule
%! [geometry, msh, space_v, vel, space_p, press] = ...
%!                        solve_stokes (problem_data, method_data);
%! error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex);
%! [error_h1_v, error_l2_v] = ...
%!    sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex);
%! assert (msh.nel, 100)
%! assert (space_p.ndof, 169)
%! assert (space_v.ndof, 1058)
%! assert (error_l2_p, 1.95560715994218e-08, 1e-15)
%! assert (error_h1_v, 8.98425390938543e-08, 1e-15)
%! assert (error_l2_v, 1.37608469955346e-09, 1e-15)


text:
function dcentral = dcentral_pair(trx,larva1,larva2)

% initialize
dcentral = nan(1,trx(larva1).nframes);

% get start and end frames of overlap
t0 = max(trx(larva1).firstframe,trx(larva2).firstframe);
t1 = min(trx(larva1).endframe,trx(larva2).endframe);

% no overlap
if t1 < t0,
  return;
end

% indices for these frames
i0 = t0 + trx(larva1).off;
i1 = t1 + trx(larva1).off;
j0 = t0 + trx(larva2).off;
j1 = t1 + trx(larva2).off;

% centroid distance
dx = trx(larva2).xcentral_mm(j0:j1)-trx(larva1).xcentral_mm(i0:i1);
dy = trx(larva2).ycentral_mm(j0:j1)-trx(larva1).ycentral_mm(i0:i1);
z = sqrt(dx.^2 + dy.^2);
dcentral(i0:i1) = z;