function [cout,H,CS] = contourf(varargin)
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      


%elseif (nin == 2),
[z,nv] = deal(varargin{1:2});
[mz,nz] = size(z);
x = 1:nz;
y = (1:mz)';

%levels
lims = [min(varargin{1}(:)),max(varargin{1}(:)),  min(varargin{2}(:)),max(varargin{2}(:))];


i = find(isfinite(z));
minz = min(z(i));
maxz = max(z(i));

% Don't fill contours below the lowest level specified in nv.
% To fill all contours, specify a value of nv lower than the
% minimum of the surface. 
draw_min=0;
if any(nv <= minz),
  draw_min=1;
end

% Get the unique levels
nv = sort([minz nv(:)']);
zi = [1, find(diff(nv))+1];
nv = nv(zi);

% Surround the matrix by a very low region to get closed contours, and
% replace any NaN with low numbers as well.

zz=[ repmat(NaN,1,nz+2) ; repmat(NaN,mz,1) z repmat(NaN,mz,1) ; repmat(NaN,1,nz+2)];
kk=find(isnan(zz(:)));
zz(kk)=minz-1e4*(maxz-minz)+zeros(size(kk));

xx = [2*x(:,1)-x(:,2), x, 2*x(:,nz)-x(:,nz-1)];
yy = [2*y(1,:)-y(2,:); y; 2*y(mz,:)-y(mz-1,:)];

[CS,msg]=contours(xx,yy,zz,nv);


% Find the indices of the curves in the c matrix, and get the
% area of closed curves in order to draw patches correctly. 
ii = 1;
ncurves = 0;
I = [];
Area=[];
while (ii < size(CS,2)),
  nl=CS(2,ii);
  ncurves = ncurves + 1;
  I(ncurves) = ii;
  xp=CS(1,ii+(1:nl));  % First patch
  yp=CS(2,ii+(1:nl));
  Area(ncurves)=sum( diff(xp).*(yp(1:nl-1)+yp(2:nl))/2 );
  ii = ii + nl + 1;
end

newplot;
if ~ishold,
  view(2);
  set(gca,'box','on');
  set(gca,'xlim',lims(1:2),'ylim',lims(3:4))
end

% Plot patches in order of decreasing size. This makes sure that
% all the levels get drawn, not matter if we are going up a hill or
% down into a hole. When going down we shift levels though, you can
% tell whether we are going up or down by checking the sign of the
% area (since curves are oriented so that the high side is always
% the same side). Lowest curve is largest and encloses higher data
% always.

H=[];
[FA,IA]=sort(-abs(Area));
if ~isstr(get(gca,'color')),
  bg = get(gca,'color');
else
  bg = get(gcf,'color');
end
if isempty(col)
  edgec = get(gcf,'defaultsurfaceedgecolor');
else
  edgec = col;
end
if isempty(lin)
  edgestyle = get(gcf,'defaultpatchlinestyle');
else
  edgestyle = lin;
end

% Tolerance for edge comparison
xtol = 0.1*(lims(2)-lims(1))/size(z,2);
ytol = 0.1*(lims(4)-lims(3))/size(z,1);

if nargout>0
  cout = [];
end
for jj=IA,
  nl=CS(2,I(jj));
  lev=CS(1,I(jj));
  if (lev ~= minz | draw_min ),
    xp=CS(1,I(jj)+(1:nl));  
    yp=CS(2,I(jj)+(1:nl));
    clev = lev;           % color for filled region above this level
    if (sign(Area(jj)) ~=sign(Area(IA(1))) ),
      kk=find(nv==lev);
      kk0 = 1 + sum(nv<=minz) * (~draw_min);
      if (kk > kk0)
        clev=nv(kk-1);    % in valley, use color for lower level
      elseif (kk == kk0)
        clev=NaN;
      else 
        clev=NaN;         % missing data section
        lev=NaN;
      end
    end

    if (isfinite(clev)),
      H=[H;patch(xp,yp,clev,'facecolor','flat','edgecolor',edgec, ...
              'linestyle',edgestyle,'userdata',lev)];
    else
      H=[H;patch(xp,yp,clev,'facecolor',bg,'edgecolor',edgec, ...
              'linestyle',edgestyle,'userdata',CS(1,I(jj)))];
    end
    
    if nargout>0
      % Ignore contours that lie along a boundary
      
      % Get +1 along lower boundary, -1 along upper, 0 in middle
      tx = (abs(xp - lims(1)) < xtol ) - (abs(xp - lims(2)) < xtol);
      ty = (abs(yp - lims(3)) < ytol ) - (abs(yp - lims(4)) < ytol);

      % Locate points with a boundary contour segment leading up to them
      bcf = find((tx & [0 ~diff(tx)]) | (ty & [0 ~diff(ty)]));

      if (~isempty(bcf))
         % Get a logical vector that has 0 inserted before each such location
         wuns = logical(ones(1,length(xp) + length(bcf)));
         wuns(bcf + (0:(length(bcf)-1))) = 0;

         % Create new arrays so that NaN breaks each boundary contour segment
         xp1 = NaN * wuns;
         yp1 = xp1;
         xp1(wuns) = xp;
         yp1(wuns) = yp;

         % Remove unnecessary elements
         if (length(xp1) > 2)
            % Blank out segments consisting of a single point
            tx = ([1 isnan(xp1(1:end-1))] & [isnan(xp1(2:end)) 1]);
            xp1(tx) = NaN;
            
            % Remove consecutive NaNs or NaNs on either end
            tx = isnan(xp1) & [isnan(xp1(2:end)) 1];
            xp1(tx) = [];
            yp1(tx) = [];
            if (length(xp1)>2 & isnan(xp1(1)))
               xp1 = xp1(2:end);
               yp1 = yp1(2:end);
            end
            
            % No empty contours allowed
            if (length(xp1) == 0)
               xp1 = NaN;
               yp1 = NaN;
            end
         end

         % Update the contour segments and their length
         xp = xp1;
         yp = yp1;
         nl = length(xp);
      end

      cout = [cout,[lev xp;nl yp]];
    end
  end
end

numPatches = length(H);
if numPatches>1
  for i=1:numPatches
    set(H(i), 'faceoffsetfactor', 0, 'faceoffsetbias', (1e-3)+(numPatches-i)/(numPatches-1)/30); 
  end
end
