function [cout,hout] = tricontour(p,t,Hn,N)

% Contouring for functions defined on triangular meshes
%
%   TRICONTOUR(p,t,F,N)
%
% Draws contours of the surface F, where F is defined on the triangulation
% [p,t]. These inputs define the xy co-ordinates of the points and their
% connectivity:
%
%   P = [x1,y1; x2,y2; etc],            - xy co-ordinates of nodes in the 
%                                         triangulation
%   T = [n11,n12,n13; n21,n23,n23; etc] - node numbers in each triangle
%
% The last input N defines the contouring levels. There are several
% options:
%
%   N scalar - N number of equally spaced contours will be drawn
%   N vector - Draws contours at the levels specified in N
%
% A special call with a two element N where both elements are equal draws a
% single contour at that level.
%
%   [C,H] = TRICONTOUR(...)
%
% This syntax can be used to pass the contour matrix C and the contour
% handels H to clabel by adding clabel(c,h) or clabel(c) after the call to
% TRICONTOUR.
%
% TRICONTOUR can also return 3D contours similar to CONTOUR3 by adding
% view(3) after the call to TRICONTOUR.
%
% Type "contourdemo" for some examples.
%
% See also, CONTOUR, CLABEL

% This function does NOT interpolate back onto a Cartesian grid, but
% instead uses the triangulation directly.
%
% If your going to use this inside a loop with the same [p,t] a good
% modification is to make the connectivity "mkcon" once outside the loop
% because usually about 50% of the time is spent in "mkcon".
%
% Darren Engwirda - 2005 (d_engwirda@hotmail.com)
% Updated 15/05/2006


% I/O checking
if nargin~=4
    error('Incorrect number of inputs')
end
if nargout>2
    error('Incorrect number of outputs')
end

% Error checking
if (size(p,2)~=2) || (size(t,2)~=3) || (size(Hn,2)~=1)
    error('Incorrect input dimensions')
end
if size(p,1)~=size(Hn,1)
    error('F and p must be the same length')
end
if (max(t(:))>size(p,1)) || (min(t(:))<=0)
    error('t is not a valid triangulation of p')
end
if (size(N,1)>1) && (size(N,2)>1)
    error('N cannot be a matrix')
end

% Make mesh connectivity data structures (edge based pointers)
[e,eINt,e2t] = mkcon(p,t);

numt = size(t,1);       % Num triangles
nume = size(e,1);       % Num edges



%==========================================================================
%                Quadratic interpolation to centroids
%==========================================================================

% Nodes
t1 = t(:,1); t2 = t(:,2); t3 = t(:,3);

% FORM FEM GRADIENTS
% Evaluate centroidal gradients (piecewise-linear interpolants)
x23 = p(t2,1)-p(t3,1);  y23 = p(t2,2)-p(t3,2);
x21 = p(t2,1)-p(t1,1);  y21 = p(t2,2)-p(t1,2);

% Centroidal values
Htx = (y23.*Hn(t1) + (y21-y23).*Hn(t2) - y21.*Hn(t3)) ./ (x23.*y21-x21.*y23);
Hty = (x23.*Hn(t1) + (x21-x23).*Hn(t2) - x21.*Hn(t3)) ./ (y23.*x21-y21.*x23);

% Form nodal gradients.
% Take the average of the neighbouring centroidal values
Hnx = 0*Hn; Hny = Hnx; count = Hnx;
for k = 1:numt
    % Nodes
    n1 = t1(k); n2 = t2(k); n3 = t3(k);
    % Current values
    Hx = Htx(k); Hy = Hty(k);
    % Average to n1
    Hnx(n1)   = Hnx(n1)+Hx;
    Hny(n1)   = Hny(n1)+Hy;
    count(n1) = count(n1)+1;
    % Average to n2
    Hnx(n2)   = Hnx(n2)+Hx;
    Hny(n2)   = Hny(n2)+Hy;
    count(n2) = count(n2)+1;
    % Average to n3
    Hnx(n3)   = Hnx(n3)+Hx;
    Hny(n3)   = Hny(n3)+Hy;
    count(n3) = count(n3)+1;
end
Hnx = Hnx./count;
Hny = Hny./count;

% Centroids [x,y]
pt = (p(t1,:)+p(t2,:)+p(t3,:))/3;

% Take unweighted average of the linear extrapolation from nodes to centroids
Ht = ( Hn(t1) + (pt(:,1)-p(t1,1)).*Hnx(t1) + (pt(:,2)-p(t1,2)).*Hny(t1) + ...
       Hn(t2) + (pt(:,1)-p(t2,1)).*Hnx(t2) + (pt(:,2)-p(t2,2)).*Hny(t2) + ...
       Hn(t3) + (pt(:,1)-p(t3,1)).*Hnx(t3) + (pt(:,2)-p(t3,2)).*Hny(t3) )/3;



% DEAL WITH CONTOURING LEVELS
if length(N)==1
    lev = linspace(max(Ht),min(Ht),N+1);
    num = N;
else
    if (length(N)==2) && (N(1)==N(2))
        lev = N(1);
        num = 1;
    else
        lev = sort(N);
        num = length(N);
        lev = lev(num:-1:1);
    end
end

% MAIN LOOP
c   = [];
h   = [];
in  = false(numt,1);
vec = 1:numt;
old = in;
for v = 1:num       % Loop over contouring levels
    
    % Find centroid values >= current level
    i     = vec(Ht>=lev(v));
    i     = i(~old(i));         % Don't need to check triangles from higher levels
    in(i) = true;
    
    % Locate boundary edges in group
    bnd  = [i; i; i];       % Just to alloc
    next = 1;
    for k = 1:length(i)
        ct    = i(k);
        count = 0;
        for q = 1:3     % Loop through edges in ct
            ce = eINt(ct,q);
            if ~in(e2t(ce,1)) || ((e2t(ce,2)>0)&&~in(e2t(ce,2)))    
                bnd(next) = ce;     % Found bnd edge
                next      = next+1;
            else
                count = count+1;    % Count number of non-bnd edges in ct
            end
        end
        if count==3                 % If 3 non-bnd edges ct must be in middle of group
            old(ct) = true;         % & doesn't need to be checked for the next level
        end
    end
    numb = next-1; bnd(next:end) = [];
    
    % Skip to next lev if empty
    if numb==0
        continue
    end
    
    % Place nodes approximately on contours by interpolating across bnd
    % edges    
    t1  = e2t(bnd,1);
    t2  = e2t(bnd,2);
    ok  = t2>0;
    
    % Get two points for interpolation. Always use t1 centroid and 
    % use t2 centroid for internal edges and bnd midpoint for boundary 
    % edges
    
    % 1st point is always t1 centroid
    H1 = Ht(t1);                                                % Centroid value
    p1 = ( p(t(t1,1),:)+p(t(t1,2),:)+p(t(t1,3),:) )/3;          % Centroid [x,y]
    
    % 2nd point is either t2 centroid or bnd edge midpoint
    i1        = t2(ok);                                         % Temp indexing
    i2        = bnd(~ok);
    H2        = H1;
    H2(ok)    = Ht(i1);                                         % Centroid values internally
    H2(~ok)   = ( Hn(e(i2,1))+Hn(e(i2,2)) )/2;                  % Edge values at boundary
    p2        = p1;
    p2(ok,:)  = ( p(t(i1,1),:)+p(t(i1,2),:)+p(t(i1,3),:) )/3;   % Centroid [x,y] internally
    p2(~ok,:) = ( p(e(i2,1),:)+p(e(i2,2),:) )/2;                % Edge [x,y] at boundary
    
    % Linear interpolation
    r     = (lev(v)-H1)./(H2-H1);
    penew = p1 + [r,r].*(p2-p1);
    
    % Do a temp connection between adjusted node & endpoint nodes in
    % ce so that the connectivity between neighbouring adjusted nodes
    % can be determined
    vecb    = (1:numb)';
    m       = 2*vecb-1;
    c1      = 0*m;
    c2      = 0*m;
    c1(m)   = e(bnd,1);
    c1(m+1) = e(bnd,2);
    c2(m)   = vecb;
    c2(m+1) = vecb;
    
    % Sort connectivity to place connected edges in sucessive rows
    [c1,i] = sort(c1); c2 = c2(i);
    
    % Connect adjacent adjusted nodes
    k    = 1;
    next = 1;
    while k<(2*numb)
        if c1(k)==c1(k+1)
            c1(next) = c2(k);
            c2(next) = c2(k+1);
            next     = next+1;
            k        = k+2;         % Skip over connected edge
        else
            k = k+1;                % Node has only 1 connection - will be picked up above
        end
    end
    ncc          = next-1; 
    c1(next:end) = []; 
    c2(next:end) = [];
    
    
    % Plot the contours
    % If an output is required, extra sorting of the
    % contours is necessary for CLABEL to work.   
    if nargout>0
        
        % Form connectivity for the contour, connecting 
        % its edges (rows in cc) with its vertices.
        ndx = repmat(1,nume,1);
        n2e = 0*penew;
        for k = 1:ncc
            % Vertices
            n1 = c1(k); n2 = c2(k);
            % Connectivity
            n2e(n1,ndx(n1)) = k; ndx(n1) = ndx(n1)+1;
            n2e(n2,ndx(n2)) = k; ndx(n2) = ndx(n2)+1;
        end
        bndn = n2e(:,2)==0;         % Boundary nodes
        bnde = bndn(c1)|bndn(c2);   % Boundary edges
        
        % Alloc some space
        tmpv = repmat(0,1,ncc);
        
        % Loop through the points at the current contour level (lev(v))
        % Try to assemble the CS data structure introduced in "contours.m"
        % so that clabel will work. Assemble CS by "walking" around each 
        % subcontour segment contiguously.
        ce    = 1;
        start = ce;
        next  = 2;
        cn    = c2(1);
        flag  = false(ncc,1);        
        x     = tmpv; x(1) = penew(c1(ce),1);
        y     = tmpv; y(1) = penew(c1(ce),2);
        for k = 1:ncc
            
            % Checked this edge
            flag(ce) = true;
            
            % Add vertices to patch data
            x(next) = penew(cn,1);
            y(next) = penew(cn,2);
            next    = next+1;
            
            % Find edge (that is not ce) joined to cn
            if ce==n2e(cn,1)
                ce = n2e(cn,2);
            else
                ce = n2e(cn,1);
            end
            
            % Check the new edge
            if (ce==0)||(ce==start)||(flag(ce))     
               
                % Plot current subcontour as a patch and save handles
                x   = x(1:next-1);
                y   = y(1:next-1);
                z   = repmat(lev(v),1,next);
                h   = [h; patch('Xdata',[x,NaN],'Ydata',[y,NaN],'Zdata',z, ...
                                'Cdata',z,'facecolor','none','edgecolor','flat')]; hold on      
                
                % Update the CS data structure as per "contours.m"
                % so that clabel works
                c = horzcat(c,[lev(v), x; next-1, y]);
                
                if all(flag)    % No more points at lev(v)
                    break
                else            % More points, but need to start a new subcontour
                    
                    % Find the unflagged edges
                    edges = find(~flag);
                    ce    = edges(1);
                    % Try to select a boundary edge so that we are 
                    % not repeatedly running into the boundary
                    for i = 1:length(edges)
                        if bnde(edges(i))
                            ce = edges(i); break
                        end
                    end
                    % Reset counters
                    start = ce;
                    next  = 2;
                    % Get the non bnd node in ce
                    if bndn(c2(ce))
                        cn = c1(ce);
                        % New patch vectors
                        x = tmpv; x(1) = penew(c2(ce),1);
                        y = tmpv; y(1) = penew(c2(ce),2);
                    else
                        cn = c2(ce);
                        % New patch vectors
                        x = tmpv; x(1) = penew(c1(ce),1);
                        y = tmpv; y(1) = penew(c1(ce),2);
                    end                    
                    
                end
            
            else                            
                % Find node (that is not cn) in ce
                if cn==c1(ce)
                    cn = c2(ce);
                else
                    cn = c1(ce);
                end
            end
            
        end
        
    else        % Just plot the contours as is, this is faster...
        
        z = repmat(lev(v),2,ncc);
        
        patch('Xdata',[penew(c1,1),penew(c2,1)]', ...
              'Ydata',[penew(c1,2),penew(c2,2)]', ...
              'Zdata',z,'Cdata',z,'facecolor','none','edgecolor','flat'); 
          
        hold on
        
    end
    
end

% Assign outputs if needed
if nargout>0
    cout = c;
    hout = h;
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [e,eINt,e2t] = mkcon(p,t)

numt = size(t,1);
vect = 1:numt;

% DETERMINE UNIQUE EDGES IN MESH
 
e       = [t(:,[1,2]); t(:,[2,3]); t(:,[3,1])];             % Edges - not unique
vec     = (1:size(e,1))';                                   % List of edge numbers
[e,j,j] = unique(sort(e,2),'rows');                         % Unique edges
vec     = vec(j);                                           % Unique edge numbers
eINt    = [vec(vect), vec(vect+numt), vec(vect+2*numt)];    % Unique edges in each triangle

% DETERMINE EDGE TO TRIANGLE CONNECTIVITY

% Each row has two entries corresponding to the triangle numbers
% associated with each edge. Boundary edges have one entry = 0.
nume = size(e,1);
e2t  = repmat(0,nume,2);
ndx  = repmat(1,nume,1);
for k = 1:numt
    % Edge in kth triangle
    e1 = eINt(k,1); e2 = eINt(k,2); e3 = eINt(k,3);
    % Edge 1
    e2t(e1,ndx(e1)) = k; ndx(e1) = ndx(e1)+1;
    % Edge 2
    e2t(e2,ndx(e2)) = k; ndx(e2) = ndx(e2)+1;
    % Edge 3
    e2t(e3,ndx(e3)) = k; ndx(e3) = ndx(e3)+1;
end

return
