Thу post is the answer to http://www.mapleprimes.com/questions/200355-Clever-Cutter-3?reply=reply

This message contains two procedures . The first procedure called **LargestIncircle** symbolically looks for the largest circle inscribed in a convex polygon . Polygon must be specified as a list of lists of its vertices . If the polygon does not have parallel sides , the problem has a unique solution , otherwise may be several solutions . It is easy to prove that there is always a maximum circle that touches three sides of the polygon . Procedure code based on this property of the optimal solution.

The second procedure called **MinCoveringCircle** symbolically seeking smallest circle that encloses a given set of points. The set of points can be anyone not necessarily being vertices of a convex polygon. Procedure code based on the following geometric properties of the optimal solution , which can easily be proved :

1) The minimum covering circle is unique.

2) The minimum covering circle of a set *Set* can be determined by at most three points in *Set* which lie on the boundary of the circle. If it is determined by only two points, then the line segment joining those two points must be a diameter of the minimum circle. If it is determined by three points, then the triangle consisting of those three points is not obtuse.

Code of the 1st procedure:

**restart;**

**LargestIncircle:=proc(L::listlist)**

**local n, x, y, P, S, T, Eqs, IsInside, OC, Pt, E, R, t, Sol, P0, R1, Circle;**

**uses combinat, geometry;**

** **

**n:=nops(L); _EnvHorizontalName := x; _EnvVerticalName := y;**

**P:=[seq(point(A||i, L[i]), i=1..n)];**

**if nops(convexhull(P))<n then error "L should be a convex polygon" fi;**

**S:={seq(line(l||i,[A||i, A||(i+1)]), i=1..n-1), line(l||n,[A||n, A||1])};**

**Eqs:=map(lhs,{seq(Equation(l||i), i=1..n)});**

**T:=choose({seq(l||i, i=1..n)}, 3);**

** **

**IsInside:=proc(Point, OC, Eqs)**

**convert([seq(is(eval(Eqs[i], {x=Point[1], y=Point[2]})*eval(Eqs[i], {x=OC[1], y=OC[2]})>0), i=1..n)], `and`);**

**end proc;**

** **

**OC:=[add(L[i,1], i=1..n)/n, add(L[i,2], i=1..n)/n];**** **

**R:=0; point(Pt,[x,y]);**

** for t in T do**

** solve([distance(Pt,t[1])=distance(Pt,t[2]), distance(Pt,t[2])=distance(Pt,t[3])],[x,y]);**

** Sol:=map(a->[rhs(a[1]), rhs(a[2])], %);**

** P0:=op(select(b->IsInside(b, OC, Eqs), Sol)); if P0<>NULL then R1:=distance(point(E,P0), t[1]);**

** if convert([seq(is(R1<=distance(E, i)), i=S minus t)], `and`) and is(R1>R) then R:=R1; Circle:=[P0, R] fi; fi;**

** od;**** **

**Circle;**** **

**end proc:**

Code the 2nd procedure:

**MinCoveringCircle:=proc(Set::{set, list}(list))**

**local S, T, t, d, A, B, C, P, R, R1, O1, Coord, L;**

**uses combinat, geometry;**** **

**if type(Set, list) then S:=convert(Set, set) else S:=Set fi;**

**T:=choose(S, 2);**

** for t in T do**

** d:=simplify(distance(point(A, t[1]), point(B, t[2]))); midpoint(C, A, B);**

** if convert([seq(is(distance(point(P, p), C)<=d/2), p in S minus {t[1], t[2]})], `and`) then return [coordinates(C), d/2] fi;**

** od;**

**T:=choose(S, 3);**

**R:=infinity;**

** for t in T do**

** point(A, t[1]); point(B, t[2]); point(C, t[3]);**

** if is(distance(A, B)^2<=distance(B, C)^2+distance(A, C)^2 and distance(A,C)^2<=distance(B, C)^2+distance(A, B)^2 and distance(B, C)^2 <= distance(A, C)^2 + distance(A, B)^2) then**

** circle(c, [A,B,C],'centername'=O1); R1:=radius(c); Coord:=coordinates(O1);**

** if convert([seq(is(distance(point(P, p), O1)<=R1), p in S minus {t[1], t[2], t[3]})], `and`) and is(R1<R) then R:=R1; L:=[Coord, R] fi; fi;**

** od;**** **

**L;**** **

**end proc: **

Examples of using the first procedure:

**L:=[[0,2],[1,4],[4,4],[5,1],[3,0]]:**

**C:=LargestIncircle(L);**

**A:=plottools[disk](op(C), color=green):**

**B:=plottools[polygon](L, style=line, thickness=2):**

**plots[display](A,B, scaling=constrained);**

The procedure can be used to find the largest circle inscribed in a convex shape, bounded by the curves, if these curves are approximated by broken lines:

**L:=[seq([k,k^2], k=-1..2, 0.1)]:**

**C:=LargestIncircle(L);**

**A:=plottools[disk](op(C), color=green):**

**B:=plottools[polygon](L, style=line, thickness=2):**

**plots[display](A,B, scaling=constrained);**

** **

Examples of using the second procedure:

**L:=[[0,2],[1,4],[4,4],[5,1],[3,0]]:**

**C:=MinCoveringCircle(L);**

**A:=plottools[circle](op(C), color=red, thickness=2):**

**B:=plottools[polygon](L, color=cyan):**

**plots[display](A,B);**

The smallest circle covering 30 random points:

> **roll:=rand(1..10):**

**L:=[seq([roll(), roll()], i=1..30)]:**

**C:=MinCoveringCircle(L);**

**A:=plottools[circle](op(C), color=red, thickness=2):**

**B:=seq(plottools[disk](L[i],0.07,color=blue), i=1..nops(L)):**

**C1:=plottools[disk](C[1], 0.07, color=red):**

**T:=plots[textplot]([C[1,1]-0.2, C[1,2]-0.2, "C"], font=[TIMES,ROMAN,16]):**

**plots[display](A, B, C1, T);**