Thread Subject: The largest Triangle that can fit in convex hull

Hi
Could any one tell me how to find the largest triangle can fit into a
convex hull.
Thanks in advance.

Before anyone can attempt to answer that question, you have to define what
you mean by "large".
Do you mean the triangle that has the largest diameter, the triangle that
has the largest area, or what ?
Also, do you need an accurate solution or are you willing to accept an
approximation?

Damian

"Dodo" <lovelygirl_2220@hotmail.com> wrote in message
news:b9bb3ca1-66f4-4c5f-9553-cd970bfb7eb1@v4g2000vba.googlegroups.com...
> Hi
> Could any one tell me how to find the largest triangle can fit into a
> convex hull.
> Thanks in advance.

On May 4, 9:18 pm, "Damian Sheehy" <Damian.She...@mathworks.com>
wrote:
> Before anyone can attempt to answer that question, you have to define wha=
t
> you mean by "large".
> Do you mean thetrianglethat has the largest diameter, thetrianglethat
> has the largest area, or what ?
> Also, do you need an accurate solution or are you willing to accept an
> approximation?
>
> Damian
>
> "Dodo" <lovelygirl_2...@hotmail.com> wrote in message
>
> news:b9bb3ca1-66f4-4c5f-9553-cd970bfb7eb1@v4g2000vba.googlegroups.com...
>
> > Hi
> > Could any one tell me how to find the largesttrianglecanfitintoa
> >convexhull.
> > Thanks in advance.

Hi Damian
i want the largest-area triangle
i want an accurate solution but if there is not you can tell me an
approximate one and i will test it
Thanks in advance

Two observations:

- If T is a triangle that maximizes the area, we can always find a triangle with the same area where the corners are coincides with the nodes of the conver hull (since triangle area must be monotonic when corner moves along an edge).
- Given two arbitrary nodes of the convex hull where two corners of the triangles coincide with. Take this as the constraint, its is easy to find the third corner (that maximizes the area) in constant time.

So with this two observations, we can do this method: just generate a list of all the couples of nodes of the convex hull, then compute the list of respective third corner that maximizes the area, and compare the areas, take the maximum, and you are done. The complexity is O(N^2), where N is number of points in the convex hull.

Does it sound OK for you?

Bruno

Dodo <lovelygirl_2220@hotmail.com> wrote in message <a3eecd7c-3fb1-4f8f-ba6c-0e39057b0279@j12g2000vbl.googlegroups.com>...
> On May 4, 9:18?pm, "Damian Sheehy" <Damian.She...@mathworks.com>
> wrote:
> > Before anyone can attempt to answer that question, you have to define wha=
> t
> > you mean by "large".
> > Do you mean thetrianglethat has the largest diameter, thetrianglethat
> > has the largest area, or what ?
> > Also, do you need an accurate solution or are you willing to accept an
> > approximation?
> >
> > Damian
> >
> > "Dodo" <lovelygirl_2...@hotmail.com> wrote in message
> >
> > news:b9bb3ca1-66f4-4c5f-9553-cd970bfb7eb1@v4g2000vba.googlegroups.com...
> >
> > > Hi
> > > Could any one tell me how to find the largesttrianglecanfitintoa
> > >convexhull.
> > > Thanks in advance.
>
> Hi Damian
> i want the largest-area triangle
> i want an accurate solution but if there is not you can tell me an
> approximate one and i will test it
> Thanks in advance

My gut says that this is just a combinatorial problem.
The nice thing is it is easily vectorized, almost trivially
so, and it is only O(n^3) in the number of vertices
of the convex hull. Since convex hulls tend to have
a limited number of points, I don't see much of a
problem.

Given the N vertices of the convex hull, pick any
subset of three of them. There are nchoosek(N,3)
such combinations. Each subset of vertices defines
a triangle. Pick the one with the largest area.

V = randn(100,2);
edges = convhulln(V)
edges =
    49 2
     2 68
    68 27
    27 14
    14 44
    44 95
    95 84
    84 28
    28 49

hlist = unique(edges(:))
hlist =
     2
    14
    27
    28
    44
    49
    68
    84
    95

There are exactly 84 distinct triangles that can
be created from these points.

nchoosek(9,3)
ans =
    84

tri = nchoosek(hlist,3);
A = V(tri(:,1),:);
B = V(tri(:,2),:);
C = V(tri(:,3),:);
nt = size(tri,1);

areas = sum(abs(cross([B-A,zeros(nt,1)],[C-A,zeros(nt,1)])/2),2);

That triangle with maximum area:

[maxarea,ind] = max(areas);
maxarea =
       7.5423

ind =
     6

The vertices that contributed to that triangle:
tri(ind,:)
ans =
     2 14 84

The vertices that define the triangle:
V(tri(ind,:),:)
ans =
      -1.6656 -2.2023
       2.1832 0.055744
     -0.99209 2.1122

HTH,
John

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gtvifq$dur$1@fred.mathworks.com>...
> Two observations:
>
> - If T is a triangle that maximizes the area, we can always find a triangle with the same area where the corners are coincides with the nodes of the conver hull (since triangle area must be monotonic when corner moves along an edge).
> - Given two arbitrary nodes of the convex hull where two corners of the triangles coincide with. Take this as the constraint, its is easy to find the third corner (that maximizes the area) in constant time.
>
> So with this two observations, we can do this method: just generate a list of all the couples of nodes of the convex hull, then compute the list of respective third corner that maximizes the area, and compare the areas, take the maximum, and you are done. The complexity is O(N^2), where N is number of points in the convex hull.
>
> Does it sound OK for you?
>
> Bruno

I'll claim that this is an O(N^3) operation, as
pointed out in my response.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gtvjn5$979$1@fred.mathworks.com>...

>
> I'll claim that this is an O(N^3) operation, as
> pointed out in my response.
>

John, I'll claim that by book keeping correctly the couple of vertexes, the third one can be found in constant time (no need to go the full round the hull to find it). Think like a basic move of simplex method. ;-)

Do I miss something?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gtvkuh$443$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gtvjn5$979$1@fred.mathworks.com>...
>
> >
> > I'll claim that this is an O(N^3) operation, as
> > pointed out in my response.
> >
>
> John, I'll claim that by book keeping correctly the couple of vertexes, the third one can be found in constant time (no need to go the full round the hull to find it). Think like a basic move of simplex method. ;-)
>
> Do I miss something?

I don't see it. I'll claim that for any pair of vertices
in the hull, you must check each of the other N-2
vertices. (Actually, not quite, since you can exclude
the triangles you have already seen. In essence, you
must check all nchoosek(N,3) possible triangles. The
nice thing is nchoosek gives you the complete list
of triangles to check.

This is simple to do in a vectorized operation, as I
showed, but it is still O(N^3) in the size of the
convex hull.

It sounds like you think that one of the edges of
the convex hull must be an edge of the maximal
area triangle. As it turns out, the example I ran
in my response did not have that happen. The
maximal area triangle that was found had a
considerably larger area than if I had only looked
at triangles with an edge that was conincident
with an edge of the convex hull.

The maximal area was 7.5423, yet by looking
only at those triangles which shared an edge
of the convex hull, 4.7898 was the maximal area.

John

Dodo <lovelygirl_2220@hotmail.com> wrote in message <b9bb3ca1-66f4-4c5f-9553-cd970bfb7eb1@v4g2000vba.googlegroups.com>...
> Hi
> Could any one tell me how to find the largest triangle can fit into a
> convex hull.
> Thanks in advance.

Convex hull of what? All the responders are assuming that it is the convex hull of a list of vertices (i.e. a polytope), but you haven't confirmed this.

> It sounds like you think that one of the edges of
> the convex hull must be an edge of the maximal
> area triangle.

No, I don't assume such thing.

Bruno

"Matt " <xys@whatever.com>
>
> Convex hull of what?

No need for precision. A convex hull is generally defined as {x : A*x <= b }. That's enough to work with. Of course, it must be bounded so as the maximizing triangle exists.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gtvoba$hrd$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com>
> >
> > Convex hull of what?
>
> No need for precision. A convex hull is generally defined as {x : A*x <= b }. That's enough to work with. Of course, it must be bounded so as the maximizing triangle exists.
-------


The most general definition of a convex hull of a region S is the smallest convex region enclosing S. In particular, the convex hull need not be polyhedral and need not have vertices, e.g. the convex hull of an ellipse, or some more complicated curve.

Your posts seem to be assuming that the convex hull will have vertices...

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gtvo3p$2mj$1@fred.mathworks.com>...
>
> > It sounds like you think that one of the edges of
> > the convex hull must be an edge of the maximal
> > area triangle.
>
> No, I don't assume such thing.
>
> Bruno

Then there are n*(n-1)/2 pairs of vertices.

Finding the largest area triangle requires that
of the remaining vertices, you find that one
which is farthest from the edge defined by that
pair of vertices. This requires a dot product,
and a max operation. But the fact is, it will
require (n-2) such operations to find the largest
triangle given that edge asa base. This is not a
constant time operation, independent of the
number of vertices. It is an O(n) operation. So
the total computation must be O(n^3). How do
you claim to get around this?

John

> The most general definition of a convex hull of a region S is the smallest convex region enclosing S. In particular, the convex hull need not be polyhedral and need not have vertices, e.g. the convex hull of an ellipse, or some more complicated curve.
>

The example you gave can be still considered as intersection of half planes Matt. In the definition the matrix A and b can have infinity number of rows (A is just a linear operator between two vectorial spaces, no need to be R^m to R^n).

I assume the computer does not have infinity memory, of A and b is finite matrix. And yes, if they are finite, vertexes are indeed defined when A has finite size.

Bruno

> How do
> you claim to get around this?
>

Like this:

% Data
x=randn(1,10000);
y=randn(1,10000);
%
k=convhull(x,y);
X=[x(k); y(k)];

% Engine
% size
n=size(X,2);
disp(n)

% wrap around index
fwrap = @(i) mod(i-1,n)+1;

% Select two vertexes
IJ = nchoosek(1:n,2);
I=IJ(:,1); J=IJ(:,2);

% Normal
normal = [X(2,J)-X(2,I); X(1,I)-X(1,J)];
% length
l = sqrt(sum(normal.^2,1));
normal = normal ./ repmat(l,2,1);

isplit = mat2cell(1:size(IJ,1),1,n-1:-1:1);
AREAMAX = -Inf;
for k=1:length(isplit)
    idx = isplit{k};
    X1 = X(:,I(idx(1)));
    for j=1:length(idx)
        lj = l(idx(j));
        nj = normal(:,idx(j));
        if j==1
            p = nj'*X;
            [dummy imin]=min(p);
            [dummy imax]=max(p);
        else
            % This loop is O(1) complexity in average.
            while nj'*X(:,fwrap(imin+1)) < nj'*X(:,imin)
                imin = fwrap(imin+1);
            end
            while nj'*X(:,fwrap(imax+1)) > nj'*X(:,imax)
                imax = fwrap(imax+1);
            end
        end
        
        % Compute area
        Tmin = 0.5*lj*abs(nj'*(X(:,imin)-X1));
        Tmax = 0.5*lj*abs(nj'*(X(:,imax)-X1));
        % Keep track the largest
        if Tmin>AREAMAX
           AREAMAX = Tmin;
           IdxTriMax = [I(idx(j)) J(idx(j)) imin];
        end
        if Tmax>AREAMAX
           AREAMAX = Tmax;
           IdxTriMax = [I(idx(j)) J(idx(j)) imax];
        end
    end % for loop on j
end % for loop on split

disp(AREAMAX)
disp(IdxTriMax)

clf;
plot(X(1,[1:end 1]),X(2,[1:end 1]),'k-',...
     X(1,IdxTriMax([1:end 1])),X(2,IdxTriMax([1:end 1])),'ro-');
 axis equal

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gtvt1a$leg$1@fred.mathworks.com>...
> > The most general definition of a convex hull of a region S is the smallest convex region enclosing S. In particular, the convex hull need not be polyhedral and need not have vertices, e.g. the convex hull of an ellipse, or some more complicated curve.
> >
>
> The example you gave can be still considered as intersection of half planes Matt. In the definition the matrix A and b can have infinity number of rows (A is just a linear operator between two vectorial spaces, no need to be R^m to R^n).
>
> I assume the computer does not have infinity memory, of A and b is finite matrix. And yes, if they are finite, vertexes are indeed defined when A has finite size.

So you would approximate the convex hull with a polyhedron and then apply your vertex searching method? Is it clear how many vertices and where the vertices need to be placed to get a sufficiently good approximation? And couldn't the required number of vertices be prohibitvely large if your method is O(N^2) or higher?

"Matt " <xys@whatever.com> wrote in message <gu039t$m0c$1@fred.mathworks.com>...

>
> So you would approximate the convex hull with a polyhedron and then apply your vertex searching method?

No, *I* don't (OP probably can tell you yes or no). I assumed this is an input. When OP tells "convex hull", I assume he can provide either (A & b matrices) or alternatively the set of vertexes by other computation mean. These are to me the basic boiled-down representations of convex set in R^2 (on a computer). I'm not awared of other method, and if they are the only one that exist or known, then this should be implicitly understood as the representation of the geometry object.

To take a parallel, when people tells: "I would like to solve a linear equation with Matlab", I would assume they already have some how a matrix and right hand side in matlab arrays.

> Is it clear how many vertices and where the vertices need to be placed to get a sufficiently good approximation?

It's probably a valid and very important question question. But this question falls outside the scope. To take a parallel with the above example, we usually do not need to ask people: "how do you build your matrix?" in order to advise them "use backslash". This should be classified as discretization method (think the parallel with solving a PDE, which has infinity of number of dimensions), it's entirely an art by itself, and there is probably no best universal answer.

Bruno

I have cleaned up my code. I have performed few experimental runs to check the time follows O(N^2) law. It looks like it is. Though it's quite long: for 1000 vertexes, it takes about 4 full minutes, where as for hull of Gaussian random data points, a few mini seconds maximum is needed because the hull has very few vertexes (#20).

I would say do not bother with complex algorithm if the number points is less than 100. The code as it is is certainly not optimized for anything.

% Data
caseid = 2;
switch caseid
    case 1,
        m = 1e5;
        x=randn(1,m)*3;
        y=randn(1,m);
    case 2,
        % LONG, few minutes
        theta=linspace(0,2*pi,1e3);
        theta(end)=[];
        x=cos(theta);
        y=sin(theta);
end
%
k=convhull(x,y);
X=[x(k); y(k)];
clear x y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Engine
% size
n=size(X,2);

% All two-vertices combo
[J I] = find(tril(ones(n),-1));

% Normal vector to all couple of points
normal = [X(2,J)-X(2,I); X(1,I)-X(1,J)].';

% wrap around index
fwrap = @(i) mod(i-1,n)+1;

% Split with the same j
isplit = mat2cell(1:size(I,1),1,n-1:-1:1);
AREAMAX = -Inf;
innerlooptime=zeros(1,length(isplit));
for k=1:length(isplit)
    if n>20
        disp(k);
    end
    idx = isplit{k};
    DX = bsxfun(@minus, X, X(:,I(idx(1))));
    firstj = true;
    tic
    for j=idx
        nj = normal(j,:);
        if firstj
            p = nj*DX;
            [dummy imin]=min(p);
            [dummy imax]=max(p);
            firstjúlse;
        else
            % This loop is O(1) complexity in average.
            while nj*DX(:,fwrap(imin+1)) < nj*DX(:,imin)
                imin = fwrap(imin+1);
            end
            while nj*DX(:,fwrap(imax+1)) > nj*DX(:,imax)
                imax = fwrap(imax+1);
            end
        end
        % Compute area
        Tmin = 0.5*abs(nj*DX(:,imin));
        Tmax = 0.5*abs(nj*DX(:,imax));
        % Keep track the largest
        if Tmin>AREAMAX
           AREAMAX = Tmin;
           IdxTriMax = [I(j) J(j) imin];
        end
        if Tmax>AREAMAX
           AREAMAX = Tmax;
           IdxTriMax = [I(j) J(j) imax];
        end
    end % for loop on j
    innerlooptime(k)=toc;
end % for loop on split

clear normal I J

% Display
fprintf('Dimensions = %d\n', n)
fprintf('AREAMAX = %g\n', AREAMAX)
fprintf('IdxTriMax @ [%d %d %d]\n', IdxTriMax)

%%%%%%%% Graphic
clf;
subplot(1,2,1);
plot(n-1:-1:1,innerlooptime);
xlabel('size');
ylabel('Inner loop time [s]');
subplot(1,2,2);
plot(X(1,[1:end 1]),X(2,[1:end 1]),'k-',...
     X(1,IdxTriMax([1:end 1])),X(2,IdxTriMax([1:end 1])),'ro-');
axis equal

% Bruno

I have optimized the code using profiler, put it in a mfile function and now it takes about 18s for 1000 vertexes points. ;-)

function result = tmaxt(caseid)

if nargin<1
    caseid = 1;
end


% Data
switch caseid
    case 1,
        m = 1e5;
        x=randn(1,m)*3;
        y=randn(1,m);
    case 2,
        % LONG, 18 seconds
        theta=linspace(0,2*pi,1e3);
        theta(end)=[];
        x=cos(theta);
        y=sin(theta);
end
%
k=convhull(x,y);
X=[x(k); y(k)];
clear x y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Engine
% size
n=size(X,2);

% All two-vertices combo
[J I] = find(tril(ones(n),-1));

% Normal vector to all couple of points
normal = [X(2,J)-X(2,I); X(1,I)-X(1,J)].';

% Split with the same j
isplit = mat2cell(1:size(I,1),1,n-1:-1:1);
AREAMAX = -Inf;
innerlooptime=zeros(1,length(isplit));
% array used to move to the next point by indexing it
nextwrap = [2:n 1];
for k=1:length(isplit)
% if n>50
% disp(k);
% end
    idx = isplit{k};
    DX = bsxfun(@minus, X, X(:,I(idx(1))));
    firstj = true;
    tic
    for j=idx
        nj = normal(j,:);
        if firstj
            p = nj*DX;
            [Tmin imin] = min(p);
            [Tmax imax] = max(p);
            firstjúlse;
        else
            % Compute area
            Tmin = nj*DX(:,imin);
            Tmax = nj*DX(:,imax);
            % This loop is O(1) complexity in average.
            while true
                next = nextwrap(imin);
                Tnext = nj*DX(:,next);
                if Tnext < Tmin
                    Tmin = Tnext;
                    imin = next;
                else
                    break
                end
            end
            while true
                next = nextwrap(imax);
                Tnext = nj*DX(:,next);
                if Tnext > Tmax
                    Tmax = Tnext;
                    imax = next;
                else
                    break
                end
            end
        end

        % Keep track the largest
        if abs(Tmin)>AREAMAX
           AREAMAX = abs(Tmin);
           IdxTriMax = [I(j) J(j) imin];
        end
        if abs(Tmax)>AREAMAX
           AREAMAX = abs(Tmax);
           IdxTriMax = [I(j) J(j) imax];
        end
    end % for loop on j
    innerlooptime(k)=toc;
end % for loop on split

AREAMAX = AREAMAX/2;

clear normal I J

% Display
fprintf('Dimensions = %d\n', n)
fprintf('AREAMAX = %g\n', AREAMAX)
fprintf('IdxTriMax @ [%d %d %d]\n', IdxTriMax)

result = struct('X',X,...
             'AREAMAX',AREAMAX, ...
             'IdxTriMax', IdxTriMax, ...
             'innerlooptime', innerlooptime);

%%%%%%%% Graphic
clf;
subplot(1,2,1);
plot(n-1:-1:1,innerlooptime);
xlabel('size');
ylabel('Inner loop time [s]');
subplot(1,2,2);
plot(X(1,[1:end 1]),X(2,[1:end 1]),'k-',...
     X(1,IdxTriMax([1:end 1])),X(2,IdxTriMax([1:end 1])),'ro-');
axis equal

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gu0lbe$g1v$1@fred.mathworks.com>...

> No, *I* don't (OP probably can tell you yes or no). I assumed this is an input. When OP tells "convex hull", I assume he can provide either (A & b matrices) or alternatively the set of vertexes by other computation mean.
-------------

Then, as I said initially, you are effectively assuming the given shape is a polytope. If the OP told us the hull was that of an elliptic boundary curve, for example, you would use a very different method to find the max-area triangle

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gu0lbe$g1v$1@fred.mathworks.com>...
 These are to me the basic boiled-down representations of convex set in R^2 (on a computer). I'm not awared of other method, and if they are the only one that exist or known, then this should be implicitly understood as the representation of the geometry object.
-----------------------

Just to be clear, the OP could be talking about a convex hull bounded by curves, e.g.

{(x,y): x.^2<=y<=1-x^2}

In this case, you would not look for an approximation {x: Ax<=b} to describe this convex hull. You would use fmincon to find the max-area triangle.

My code now runs about 12.5s for 1000 points; and 50s for 2000 points. ;-)

I stop to post the code from now on. If interested please email me.

Bruno

On May 8, 9:41 pm, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> My code now runs about 12.5s for 1000 points; and 50s for 2000 points. ;-=
)
>
> I stop to post the code from now on. If interested please email me.
>
> Bruno

can you send the last code that you had optimize and tell me what
exactly you have done it
Thanks

lovelygirl2220@gmail.com wrote in message <fd9acebb-3282-4649-9c70-9df0076bbc6e@v4g2000vba.googlegroups.com>...

>
> can you send the last code that you had optimize and tell me what
> exactly you have done it

On your way...

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gtvifq$dur$1@fred.mathworks.com>...
> Two observations:
>
> - If T is a triangle that maximizes the area, we can always find a triangle with the same area where the corners are coincides with the nodes of the conver hull (since triangle area must be monotonic when corner moves along an edge).
> - Given two arbitrary nodes of the convex hull where two corners of the triangles coincide with. Take this as the constraint, its is easy to find the third corner (that maximizes the area) in constant time.
>
> So with this two observations, we can do this method: just generate a list of all the couples of nodes of the convex hull, then compute the list of respective third corner that maximizes the area, and compare the areas, take the maximum, and you are done. The complexity is O(N^2), where N is number of points in the convex hull.
>
> Does it sound OK for you?
>
> Bruno


*********************************************************
Can u please elaborate on your two observations. I did not understand the algorithm. How the third vertex can be found in constant time?

Let T is the maximized triangle, AB is the base of T. C is the opposite triangle corner. C belongs to the convexhull. If C is on the edge of the hull, then the edge must be parallel to AB (otherwise the derivative of the area of T with respect to C is not zero). So I can move C to one of the corner of the hull edge, and the area of the triangle remains unchanged. So I can select C is a corner. Thus I can select A, B, and C.

Let A is one fixed corner, B is the second corner (varying), and C(B) is the corner such that

Area (ABC) = max(ABD) for all D in the hull.

1. When B changes monotonically on the circle (A fixes), C changes monotonically (using convex property).
2. When B runs on the whole circle, C runs less than the whole circle.

Therefore to find C with respect to B, we can scans B on the circle step-by-step; and when B move 1 step we update the maximizing C by looking forward and sequentially from its current position (prop 1). Because of (pro 2) the number of substep to find C in average is 1 or less. That means C can be found in constant time.

End of the proof.

Bruno