function plot_atlas(atlas, dim)
% PLOT_ATLAS   Skeletal visualization of one- and two-dimensional atlases
%
% PLOT_ATLAS(ATLAS, DIM)
%
% ATLAS - Manifold atlas.
% DIM   - Embedding dimension.

% This file is part of the atlas_kd toolbox, copyright (C) Michael
% Henderson, Frank Schilder, Harry Dankowicz, Erika Fotsch, of the package
% COCO (http://sourceforge.net/projects/cocotools).

hold on
switch dim
  case 2
    for i=1:numel(atlas.charts)
      chart = atlas.charts{i};
      if chart.P.n>0
        plot(chart.x(1), chart.x(2), '.', 'MarkerSize', 1);
        pts = repmat(chart.x, [1,2]) + chart.TS*[ chart.P.v{1} chart.P.v{2} ];
        if isempty(chart.bv)
          plot(pts(1,:), pts(2,:), 'ro-');
        else
          plot(pts(1,:), pts(2,:), 'bo-');
        end
      end
    end
  case 3
    for i=1:numel(atlas.charts)
      chart = atlas.charts{i};
      if chart.P.n>0
        plot3(chart.x(1), chart.x(2), chart.x(3),'.', 'MarkerSize', 1);
      end
      for j=1:chart.P.n
        for k=j+1:chart.P.n
          inter = intersectSets(chart.P.indices{j}, chart.P.indices{k});
          if numel(inter) >= 1
            pts = repmat(chart.x, [1,2]) + chart.TS*[ chart.P.v{j} chart.P.v{k} ];
            if isempty(chart.bv)
              plot3(pts(1,:), pts(2,:), pts(3,:), 'r');
            else
              plot3(pts(1,:), pts(2,:), pts(3,:), 'b', 'LineWidth', 2);
            end
          end
        end
      end
    end
  otherwise
end
box on
hold off

end

function C = intersectSets(A, B)

%  Returns the intersection of two SORTED index sets

nA = numel(A);
nB = numel(B);
C  = [];

iA = 1;
iB = 1;
while iA<=nA && iB<=nB
  if A(iA)==B(iB)
    C  = [C A(iA)]; %#ok<AGROW>
    iA = iA+1;
    iB = iB+1;
  elseif A(iA)<B(iB)
    iA = iA+1;
  else
    iB = iB+1;
  end
end

end
