%%%%%
%
%
%%%%%
function [FLIPIND] = accumulatePoints(ipts, opts)
tri = delaunay(ipts(:,1), ipts(:,2),{'QJ','QJ','QJ'});
flipped_tri_ind = find_flipped_triangles(tri, ipts, opts);
FLIPIND = [];
p = [];
if ~isempty(flipped_tri_ind)
    for i=1:length(flipped_tri_ind);
        ind = flipped_tri_ind(i);
        find = getFlippedPoint(opts, tri, ind);
        flippedind = tri(ind, find);
        avind = tri(ind, setdiff(1:3, find));
        %opts(flippedind,:) = sum(opts(avind,:), 1)/2;
        FLIPIND(end+1).point = flippedind;
        FLIPIND(end).neighbours = avind;
    end
end
return;
%%%%%
%
%
%%%%%
function flippedind = getFlippedPoint(opts, tri, ind)
order = [1 2 3; 2 1 3; 3 2 1];
for i=1:size(order,1)
    [x0, x1, x2] = getPosVectors(opts, tri, ind, order(i,:));
    d(i) = getDistance2Line(x0, x1, x2);
end
[val, ind] = min(d);
flippedind = order(ind(1), 1);
return;
%%%%%
%
%
%%%%%
function [x0, x1, x2] = getPosVectors(opts, tri, ind, order)
x0 = [opts(tri(ind,order(1)),:), 0];
x1 = [opts(tri(ind,order(2)),:), 0];
x2 = [opts(tri(ind,order(3)),:), 0];
return;
%%%%%
%
%
%%%%%
function [d] = getDistance2Line(x0, x1, x2)
d = 0;
a = (x2-x1)';
b = (x0-x1)';
if norm(a)>0
    %d = norm(b - (((a'*b)/(a'*a))*a))^2;
    d = norm(cross(x1 - x2, x1 - x0))/norm(x2-x1);
end
return;


