function hh = clabel(cs,varargin)
%CLABEL Contour plot elevation labels.
%   CLABEL(CS,H) adds height labels to the current contour
%   plot.  The labels are rotated and inserted within the contour
%   lines.  CS and H are the contour matrix output and object handle
%   outputs from CONTOUR, CONTOUR3, or CONTOURF.
%
%   CLABEL(CS,H,V) labels just those contour levels given in
%   vector V.  The default action is to label all known contours.
%   The label positions are selected randomly.
%
%   CLABEL(CS,H,'manual') places contour labels at the locations
%   clicked on with a mouse.  Pressing the return key terminates
%   labeling.  Use the space bar to enter contours and the arrow
%   keys to move the crosshair if no mouse is available.
%
%   CLABEL(CS) or CLABEL(CS,V) or CLABEL(CS,'manual') places
%   contour labels as above, except that the labels are drawn as
%   plus signs on the contour with a nearby height value.
%
%   H = CLABEL(...) returns handles to the TEXT (and possibly LINE)
%   objects created.  The UserData property of the TEXT objects contain
%   the height value for each label.
%
%   Text property/value pairs can be appended to the list of 
%   input arguments in the above:
%   CLABEL(clabel options, 'text property', property_value , ... )
%
%   One special property ('label') is also available to specify 
%   the spacing between labels (in points). This defaults to 144, or
%   2 inches. 
%
%   Example
%      subplot(1,3,1), [cs,h] = contour(peaks); clabel(cs,h,'label',72)
%      subplot(1,3,2), cs = contour(peaks); clabel(cs)
%      subplot(1,3,3), [cs,h] = contour(peaks); 
%       clabel(cs,h,'fontsize',15,'color','r','rotation',0)
%
%   See also CONTOUR, CONTOUR3, CONTOURF, LINELABEL.

%   Thanks to R. Pawlowicz (IOS) rich@ios.bc.ca for the algorithm used
%   in 'inline_labels' so that clabel can produce inline labeling.

%   Copyright (c) 1984-96 by The MathWorks, Inc.
%   $Revision: 5.19 $  $Date: 1996/10/25 17:43:58 $

% Modified by R Pawlowicz to allow for text properties as in extcontour code 14/5/97
% 28/10/97 - modified to work in map contouring
%  9/01/98 - improved calculation of gaps for line labels
%   Fix by Eric Firing, efiring@soest.hawaii.edu, 4/97, to
%   make the rotation angles correct when XDir and/or YDir are
%   reverse.

if nargin == 0
    error('Not enough input arguments.')
end
if min(size(cs)) > 2
    error('First input must be a valid contour description matrix.')
end
threeD = IsThreeD(gca);

if nargin == 1,
  h = plus_labels(threeD,cs);
else
  if ~isempty(varargin{1}(1)) & ishandle(varargin{1}(1)) & ...
    (strcmp(get(varargin{1}(1),'type'),'line') | strcmp(get(varargin{1}(1),'type'),'patch')),
    h = inline_labels(cs,varargin{:});
  else
    h = plus_labels(threeD,cs,varargin{:});
  end;
end;


if nargout>0, hh = h; end
if ~ishold, 
  if threeD, view(3), else view(2), end
end

%--------------------------------------------------------------
function H = inline_labels(CS,h,varargin)
%
% Draw the labels along the contours and rotated to match the local slope.
%

% To open up space in the contours, we rely on the order in which
% the handles h are created in CONTOUR3.  If CONTOUR3 changes you
% might need to change the algorithm below.

% Author: R. Pawlowicz IOS rich@ios.bc.ca
%         12/12/94
%         changes - R. Pawlowicz 14/5/97 - small bug in "that ole' matlab magic" fixed,
%                   also another in manual selection of locations.

manual=0; v=[]; inargs=zeros(1,length(varargin));

if nargin>=3 & isstr(varargin{1}) & strcmp(varargin{1},'manual'),
  manual = 1;
  inargs(1)=1;
end

if ~manual & nargin>=3 & ~isstr(varargin{1}),
  v = varargin{1};
  inargs(1)=1;
end;

lab_int=72*2;  % label interval (points)

for k=find(inargs==0),
 if isstr(varargin{k}) & ~isempty(findstr(varargin{k},'lab')),
   inargs([k k+1])=1;
   lab_int=varargin{k+1};
 end;
end;

varargin(find(inargs))=[]; 

if strcmp(get(h(1),'type'),'patch') & ~strcmp(get(h(1),'facecolor'),'none'),
  isfilled = 1;
else
  isfilled = 0;
end

%% EF 4/97
if (strcmp(get(gca, 'XDir'), 'reverse')), XDir = -1; else XDir = 1; end
if (strcmp(get(gca, 'YDir'), 'reverse')), YDir = -1; else YDir = 1; end
%%

% Compute scaling to make sure printed output looks OK. We have to go via
% the figure's 'paperposition', rather than the the absolute units of the
% axes 'position' since those would be absolute only if we kept the 'units'
% property in some absolute units (like 'points') rather than the default
% 'normalized'.

UN=get(gca,'units');
if (UN(1:3)=='nor'),
  UN=get(gcf,'paperunits');
  set(gcf,'paperunits','points');
  PA=get(gcf,'paperposition');
  set(gcf,'paperunits',UN);
  PA=PA.*[get(gca,'position')];
else
  set(gca,'units','points');
  PA=get(gca,'pos');
  set(gca,'units',UN); 
end;  

% Find beginning of all lines

lCS=size(CS,2);

if ~isempty(get(gca,'children')),
  XL=get(gca,'xlim');
  YL=get(gca,'ylim');
else
  iL=[];
  k=1;
  XL=[Inf -Inf];
  YL=[Inf -Inf];
  while (k<lCS),
    x=CS(1,k+(1:CS(2,k)));
    y=CS(2,k+(1:CS(2,k)));
    XL=[ min([XL(1),x]) max([XL(2),x]) ];
    YL=[ min([YL(1),y]) max([YL(2),y]) ]; 
    iL=[iL k];
    k=k+CS(2,k)+1;
  end;
  set(gca,'xlim',XL,'ylim',YL);
end;

Aspx=PA(3)/diff(XL);  % To convert data coordinates to paper (we need
                      % to do this
Aspy=PA(4)/diff(YL);  % to get the gaps for text the correct size)

H=[];

% Set up a dummy text object from which you can get text extent info
H1=text(XL(1),YL(1),'dummyarg','units','points','visible','off',varargin{:});

% Decompose contour data structure if manual mode.

if manual
   disp(' '), disp('    Please wait a moment...')
   x = []; y = []; ilist = []; klist = []; plist = [];
   ii = 0; k = 0; n = 0;
   while (1)
      k = k + 1;
      ii = ii + n + 1; if ii > lCS, break, end
      c = CS(1,ii); n = CS(2,ii); nn = 2 .* n -1;
      xtemp = zeros(nn, 1); ytemp = zeros(nn, 1);
      xtemp(1:2:nn) = CS(1, ii+1:ii+n);
      xtemp(2:2:nn) = (xtemp(1:2:nn-2) + xtemp(3:2:nn)) ./ 2;
      ytemp(1:2:nn) = CS(2, ii+1:ii+n);
      ytemp(2:2:nn) = (ytemp(1:2:nn-2) + ytemp(3:2:nn)) ./ 2;
      x = [x; xtemp]; y = [y; ytemp];   % Keep these.
      ilist = [ilist; ii(ones(nn,1))];
      klist = [klist; k(ones(nn,1))];
      plist = [plist; (1:.5:n)'];
   end
   ax = axis; 
   xmin = ax(1); xmax = ax(2); ymin = ax(3); ymax = ax(4);
   xrange = xmax - xmin; yrange = ymax - ymin;
   xylist = (x .* yrange + sqrt(-1) .* y .* xrange);
   view(2)
   disp(' ');
   disp('   Carefully select contours for labeling.')
   disp('   When done, press RETURN while the Graph window is the active window.')
end

% Get labels all at once to get the length of the longest string.  
% This allows us to call extent only once, thus speeding up this routine
if ~manual, 
  labels = getlabels(CS);
  % Get the size of the label
  set(H1,'string',repmat('9',1,size(labels,2)),'visible','on',varargin{:})
  EX=get(H1,'extent'); set(H1,'visible','off')
  len_lab=EX(3)/2;
end

ii=1; k = 0;
levels = []; positions = [];
while (ii<lCS),
  if manual
    [xx, yy, button] = ginput(1);
    if isempty(button) | isequal(button,13), break, end
    if xx < xmin | xx > xmax, break, end
    if yy < ymin | yy > ymax, break, end
    xy = xx .* yrange + sqrt(-1) .* yy .* xrange;
    dist = abs(xylist - xy);
    f = find(dist == min(dist));
    if ~isempty(f)
      f = f(1);
      ii  = ilist(f);
      k   = klist(f);
      p   = floor(plist(f));
    end
  else
    k = k+1;
  end

  if ~isfilled & k>length(h), error('Not enough contour handles.'); end

  l=CS(2,ii);
  x=CS(1,ii+(1:l));
  y=CS(2,ii+(1:l));

  lvl=CS(1,ii);

  if manual
    lab=num2str(lvl);
    % Get the size of the label
    set(H1,'string',lab,'visible','on',varargin{:})
    EX=get(H1,'extent'); set(H1,'visible','off')
    len_lab=EX(3)/2;
  else
    %RP - get rid of all blanks in label
    lab=labels(k,labels(k,:)~=' ');
    %RP - scale label length by string size instead of a fixed length
    len_lab=EX(3)/2*length(lab)/size(labels,2);
  end
  
  % RP28/10/97 - Contouring sometimes returns x vectors with NaN in them - we want to handle
  % this case!
  sx=x*Aspx;
  sy=y*Aspy;
  d=[0 sqrt(diff(sx).^2 +diff(sy).^2) ];d(isnan(d))=0;
  d=cumsum(d);
  if ~manual
    psn=[max(len_lab,lab_int+lab_int*(rand(1)-.5)):lab_int:d(l)-len_lab];
  else
    psn = min(max(max(d(p),d(2)+eps*d(2)),d(1)+len_lab),d(end)-len_lab);
%%RP    levels = [levels;lvl];
    levels = [levels;ii];
    positions = [positions;psn];
%%RP    kk = find(levels==lvl)
    kk = find(levels==ii);
    psn = positions(kk)';
  end
  lp=size(psn,2);

  if (lp>0) & isfinite(lvl)  & ...
     (isempty(v) | any(abs(lvl-v)/max(eps+abs(v)) < .00001)),
 
    Ic=sum( d(ones(1,lp),:)' < psn(ones(1,l),:) );
    Il=sum( d(ones(1,lp),:)' <= psn(ones(1,l),:)-len_lab );
    Ir=sum( d(ones(1,lp),:)' < psn(ones(1,l),:)+len_lab );
 
    Ir = max(0,min(Ir,length(d)-1));
    Il = max(0,min(Il,length(d)-1));
    Ic = max(0,min(Ic,length(d)-1));

    % Endpoints of text in data coordinates
    wl=(d(Il+1)-psn+len_lab)./(d(Il+1)-d(Il));
    wr=(psn-len_lab-d(Il)  )./(d(Il+1)-d(Il));
    xl=x(Il).*wl+x(Il+1).*wr;
    yl=y(Il).*wl+y(Il+1).*wr;
   
    wl=(d(Ir+1)-psn-len_lab)./(d(Ir+1)-d(Ir));
    wr=(psn+len_lab-d(Ir)  )./(d(Ir+1)-d(Ir));
    xr=x(Ir).*wl+x(Ir+1).*wr;
    yr=y(Ir).*wl+y(Ir+1).*wr;
   
    trot=atan2( (yr-yl)*YDir*Aspy, (xr-xl)*XDir*Aspx )*180/pi; %% EF 4/97
    backang=abs(trot)>90;
    trot(backang)=trot(backang)+180;
    
    % Text location in data coordinates 

    wl=(d(Ic+1)-psn)./(d(Ic+1)-d(Ic));
    wr=(psn-d(Ic)  )./(d(Ic+1)-d(Ic));    
    xc=x(Ic).*wl+x(Ic+1).*wr;
    yc=y(Ic).*wl+y(Ic+1).*wr;

    % Shift label over a little if in a curvy area
    shiftfrac=.5;
    
    xc=xc*(1-shiftfrac)+(xr+xl)/2*shiftfrac;
    yc=yc*(1-shiftfrac)+(yr+yl)/2*shiftfrac;
    
    % Remove data points under the label...
    % First, find endpoint locations as distances along lines
  
    dr=d(Ir)+sqrt( ((xr-x(Ir))*Aspx).^2 + ((yr-y(Ir))*Aspy).^2 );
    dl=d(Il)+sqrt( ((xl-x(Il))*Aspx).^2 + ((yl-y(Il))*Aspy).^2 );
  
    % Now, remove the data points in those gaps using that
    % ole' Matlab magic. We use the sparse array stuff instead of
    % something like:  
    %        f1=zeros(1,l); f1(Il)=ones(1,lp);
    % because the sparse functions will sum into repeated indices,
    % rather than just take the last accessed element - compare
    %   x=[0 0 0]; x([2 2])=[1 1]
    % with
    %   x=full(sparse([1 1],[2 2],[1 1],1,3))
    % (bug fix in original code 18/7/95 - RP)

    f1=full(sparse(ones(1,lp),Il,ones(1,lp),1,l));
    f2=full(sparse(ones(1,lp),Ir,ones(1,lp),1,l));
    irem=find(cumsum(f1)-cumsum(f2))+1;
    x(irem)=[];
    y(irem)=[];
    d(irem)=[];
    l=l-size(irem,2);
    
    % Put the points in the correct order...
    
    xf=[x(1:l),xl,repmat(NaN,size(xc)),xr];
    yf=[y(1:l),yl,yc,yr];

    [df,If]=sort([d(1:l),dl,psn,dr]);
  
    % ...and draw.
    %
    % Here's where we assume the order of the h(k).  
    %

    z = get(h(k),'zdata');
    if ~isfilled, % Only modify lines or patches if unfilled
      set(h(k),'xdata',[xf(If) NaN],'ydata',[yf(If) NaN])

      % Handle contour3 case (z won't be empty).
      if ~isempty(z), z = repmat(z(1),size(get(h(k),'xdata'))); end

      if strcmp(get(h(k),'type'),'patch')
        set(h(k),'cdata',lvl+[0*xf(If) 0])
      end
    end

    for jj=1:lp,
      % Handle contour3 case (z won't be empty).
      if ~isempty(z),
        H = [H;text(xc(jj),yc(jj),z(1),lab,'rotation',trot(jj), ...
             'verticalAlignment','middle','horizontalAlignment','center',...
             'clipping','on','userdata',lvl,varargin{:})];
      else
        H = [H;text(xc(jj),yc(jj),lab,'rotation',trot(jj), ...
             'verticalAlignment','middle','horizontalAlignment','center',...
             'clipping','on','userdata',lvl,varargin{:})];       
      end
    end;
  else
    if ~isfilled, % Only modify lines or patches if unfilled
      %
      % Here's another place where we assume the order of the h(k)
       %
      set(h(k),'xdata',[x NaN],'ydata',[y NaN], ...
               'zdata',[0*x 0])
      if strcmp(get(h(k),'type'),'patch')
         set(h(k),'cdata',lvl+[0*x 0])
      end
    end
  end;
  
  if ~manual
    ii=ii+1+CS(2,ii);
  end
end;
  
% delete dummy string
delete(H1);
%-------------------------------------------------------

%-------------------------------------------------------
function h = plus_labels(threeD,cs,varargin);
%
% Draw the labels as plus symbols next to text (v4 compatible)
%

%    RP - 14/5/97
%    Clay M. Thompson 6-7-96
%    Charles R. Denham, MathWorks, 1988, 1989, 1990.
cax = gca;
manual = 0;
choice = 0;


if nargin > 2
  if isstr(varargin{1}),
    if strcmp(varargin{1}, 'manual')
      varargin(1)=[];
      manual=1;
    end;
  else
    choice = 1;
    v = sort(varargin{1}(:));
    varargin(1)=[];
  end
end

[mcs, ncs] = size(cs);

% Find range of levels.
k = 1; i = 1;
while k <= ncs
   levels(i) = cs(1,k);
   i = i + 1;
   k = k + cs(2,k) + 1;
end
cmin = min(levels);
cmax = max(levels);
crange = max(abs(levels));
cdelta = abs(diff(levels)); 
cdelta = min(cdelta(cdelta > eps))/max(eps,crange); % Minimum significant change
if isempty(cdelta), cdelta = 0; end

% Decompose contour data structure if manual mode.

if manual
   disp(' '), disp('    Please wait a moment...')
   x = []; y = []; clist = []; k = 0; n = 0;
   while (1)
      k = k + n + 1; if k > ncs, break, end
      c = cs(1,k); n = cs(2,k); nn = 2 .* n -1;
      xtemp = zeros(nn, 1); ytemp = zeros(nn, 1);
      xtemp(1:2:nn) = cs(1, k+1:k+n);
      xtemp(2:2:nn) = (xtemp(1:2:nn-2) + xtemp(3:2:nn)) ./ 2;
      ytemp(1:2:nn) = cs(2, k+1:k+n);
      ytemp(2:2:nn) = (ytemp(1:2:nn-2) + ytemp(3:2:nn)) ./ 2;
      x = [x; xtemp]; y = [y; ytemp];   % Keep these.
      clist = [clist; c .* ones(2*n-1, 1)];
   end
   ax = axis; 
   xmin = ax(1); xmax = ax(2); ymin = ax(3); ymax = ax(4);
   xrange = xmax - xmin; yrange = ymax - ymin;
   xylist = (x .* yrange + sqrt(-1) .* y .* xrange);
   view(2)
   disp(' ');
   disp('   Carefully select contours for labeling.')
   disp('   When done, press RETURN while the Graph window is the active window.')
end

k = 0; n = 0; flip = 0; h = [];

while (1)

% Use GINPUT and select nearest point if manual.

   if manual
      [xx, yy, button] = ginput(1);
      if isempty(button) | isequal(button,13), break, end
      if xx < xmin | xx > xmax, break, end
      if yy < ymin | yy > ymax, break, end
      xy = xx .* yrange + sqrt(-1) .* yy .* xrange;
      dist = abs(xylist - xy);
      f = find(dist == min(dist));
      if length(f) > 0
         f = f(1); xx = x(f); yy = y(f); c = clist(f);
         okay = 1;
        else
         okay = 0;
      end
   end

% Select a labeling point randomly if not manual.

   if ~manual
      k = k + n + 1; if k > ncs, break, end
      c = cs(1, k); n = cs(2, k);
      if choice
         f = find(abs(c-v)/max(eps+abs(v)) < .00001);
         okay = length(f) > 0;
      else
         okay = 1;
      end
      if okay
         r = rand;
         j = fix(r.* (n - 1)) + 1;
         if flip, j = n - j; end
         flip = ~flip;
         if n == 1    % if there is only one point
           xx = cs(1, j+k); yy = cs(2, j+k);
         else
           x1 = cs(1, j+k); y1 = cs(2, j+k);
           x2 = cs(1, j+k+1); y2 = cs(2, j+k+1);
           xx = (x1 + x2) ./ 2; yy = (y1 + y2) ./ 2;  % Test was here; removed.
         end
      end
   end

% Label the point.

   if okay
      % Set tiny labels to zero.
      if abs(c) <= 10*eps*crange, c = 0; end
      % Determine format string number of digits
      if cdelta > 0, 
        ndigits = max(3,ceil(-log10(cdelta)));
      else
        ndigits = 3;
      end
      s = num2str(c,ndigits);
      hl = line('xdata',xx,'ydata',yy,'marker','+','erasemode','none');
      ht = text(xx, yy, s, 'verticalalignment', 'bottom', ...
             'horizontalalignment', 'left','erasemode','none', ...
             'clipping','on','userdata',c,varargin{:});
      if threeD, 
        set(hl,'zdata',c);
        set(ht,'position',[xx yy c]);
      end
      h = [h;hl];
      h = [h;ht];
   end
end
%-------------------------------------------------------

%-------------------------------------------------------
function labels = getlabels(CS)
%GETLABELS Get contour labels
v = []; i =1;
while i < size(CS,2), 
  v = [v,CS(1,i)];
  i = i+CS(2,i)+1;
end
labels = num2str(v');

%---------------------------------------------------
function threeD = IsThreeD(cax)
%ISTHREED  True for a contour3 plot
hp = findobj(cax,'type','patch');
if isempty(hp), hp = findobj(gca,'type','line'); end
if ~isempty(hp),
  % Assume a contour3 plot if z data not empty
  threeD = ~isempty(get(hp(1),'zdata'));
else
  threeD = 0;
end