function plot_atlas(atlas, i1, i2, i3)

hold on

for i=1:numel(atlas.charts)
  chart = atlas.charts{i};
  if chart.P.n>0
    plot3(chart.x(i1), chart.x(i2), chart.x(i3),'.', '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(i1,:), pts(i2,:), pts(i3,:), 'r');
        else
          plot3(pts(i1,:), pts(i2,:), pts(i3,:), 'b', 'LineWidth', 2);
        end
      end
    end
  end
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
