extractloops.m

function loops = extractloops(edges)
%
% loops=extractloops(edges)
%
% extract individual loop or polyline segment from a collection of edges
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2007/11/21
%
% input:
%    edges:  two column matrix recording the starting/ending
%             points of all edge segments
%
% output:
%    loops:  output, a single vector separated by NaN, each segment
%             is a 3D polyline or loop consisted of node IDs
%
% example:
%    edges=[1 2;2 3;1 4;3 4;7 3;1 9;5 6;6 7;10 9; 8 10;1 8;9 3;11 11;11 12];
%    loops=extractloops(edges)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

loops = [];
edges(edges(:, 1) == edges(:, 2), :) = []; % remove degenerated edges
loops = [loops, edges(1, :)];
loophead = edges(1, 1);
loopend = edges(1, end);
edges(1, :) = [];

while (~isempty(edges))
    idx = [find(edges(:, 1) == loopend)', find(edges(:, 2) == loopend)'];
    if (length(idx) > 1) % when a node with multiple connection found
        idx = idx(1);   % take the first connection and continue
    end
    if (isempty(idx)) % when an open-line segment gets to one end
        % when both open ends are found
        if (isempty([find(edges(:, 1) == loophead)', find(edges(:, 2) == loophead)']))
            loops = [loops, nan];
            loops = [loops, edges(1, :)];
            loophead = edges(1, 1);
            loopend = edges(1, end);
            edges(1, :) = [];
        else % only the first open end is found, flip and trace the other
            [loophead, loopend] = deal(loopend, loophead);
            lp = fliplr(loops);
            seg = find(isnan(lp), 1);
            if (isempty(seg))
                loops = lp;
            else
                loops = [loops(1:end - seg(1) + 1) lp(1:seg(1) - 1)];
            end
        end
        continue
    end
    if (length(idx) == 1) % tracing along a single line thread
        idx = idx(1);
        ed = edges(idx, :);
        ed(ed == loopend) = [];
        newend = ed(1);
        if (newend == loophead)  % when a loop is found
            loops = [loops loophead nan];
            edges(idx, :) = [];
            if (size(edges, 1) == 0)
                break
            end
            loops = [loops, edges(1, :)];
            loophead = edges(1, 1);
            loopend = edges(1, end);
            edges(1, :) = [];
            continue
        else
            loops = [loops, newend];
        end
        loopend = newend;
        edges(idx, :) = [];
    end
end