In order to use the following procedures with an arbitrary set of points in the plane, the convex hull of the points should first be computed.
/* * Minimum Rectangle * The minimum rectangle of a set of points is that rectangle of smallest * area containing all the points */ /* minimum_rectangle(Points, Rectangle) is true if Points are the */ /* consecutive vertices going clockwise of a convex polygon with the first */ /* vertex repeated as the last vertex, and Rectangle is the rectangle of */ /* minimal area enclosing the Points. */ /* This uses a "double rotating caliper" algorithm. */ /* e.g. minimum_rectangle([p(3,1),p(2,0),p(1,0),p(0,1),p(2,5),p(3,1)], Rect) */ /* gives Rect=rect(p(2,5),p(4.4,3.8),p(2.2,-0.6),p(-0.2,0.6)). */ minimum_rectangle(Ps, Rectangle):- Ps=[P1,P2|Ps1], square(P1, P2, P3, P4), /* Initialize calipers */ farthest( Ps, void, [P2|Ps1], P4, P1, -1.0E9, QsW), farthest(QsW, void, [P2|Ps1], P1, P2, -1.0E9, QsN), farthest(QsN, void, [P2|Ps1], P2, P3, -1.0E9, QsE), minimum_rectangle_1(Ps, Ps, QsW, QsN, QsE, 1.0E9, void, Rectangle). minimum_rectangle_1([_], _, _, _, _, _, Rect, Rect):-!. minimum_rectangle_1([P1,P2|Ps], PsCopy, QsW, QsN, QsE, Area, Rect0, Rect):- square(P1, P2, P3, P4), /* Rotate calipers */ farthest(QsW, void, PsCopy, P4, P1, -1.0E9, QsW1), farthest(QsN, void, PsCopy, P1, P2, -1.0E9, QsN1), farthest(QsE, void, PsCopy, P2, P3, -1.0E9, QsE1), minimum_rectangle_2([P1,P2|Ps], PsCopy, QsW1, QsN1, QsE1, Area, Rect0, Rect). minimum_rectangle_2([P1,P2|Ps], PsCopy, QsW, QsN, QsE, Area, _, Rect):- QsW=[Pw|_], QsN=[Pn|_], QsE=[Pe|_], perpendicular_intersection(P1, P2, Pw, Psw), perpendicular_intersection(P2, P1, Pe, Pse), square(P1, Psw, P3, _), perpendicular_intersection(Psw, P3, Pn, Pnw), area(Pse, Pnw, Psw, Area1), Area1 < Area, !, square(Pse, P2, _, P4), perpendicular_intersection(Pse, P4, Pn, Pne), Rect1=rect(Psw,Pnw,Pne,Pse), minimum_rectangle_1([P2|Ps], PsCopy, QsW, QsN, QsE, Area1, Rect1, Rect). minimum_rectangle_2([_,P2|Ps], PsCopy, QsW, QsN, QsE, Area, Rect0, Rect):- minimum_rectangle_1([P2|Ps], PsCopy, QsW, QsN, QsE, Area, Rect0, Rect). /* farthest(Qs, void, PsCopy, P1, P2, 0, [Q|Qs1]) is true if Q is the vertex */ /* in Qs which is farthest to the right of the line through P1 and P2, and */ /* Qs1 are the remaining vertices in Qs. PsCopy is a copy of the original */ /* list of vertices of the convex polygon, and allows wraparound from the */ /* last vertex to the first vertex. */ farthest([], Qprev, PsCopy, P1, P2, Area, Qs1):-!, PsCopy=[_|_], /* Fail if only two points in polygon, for example */ farthest(PsCopy, Qprev, [], P1, P2, Area, Qs1). farthest([Q|Qs], _, PsCopy, P1, P2, Area, Qs1):- area(P2, P1, Q, Area1), /* Anticlockwise */ Area1 >= Area, !, farthest(Qs, Q, PsCopy, P1, P2, Area1, Qs1). farthest(Qs, Qprev, _, _, _, _, [Qprev|Qs]). /* area(Pa, Pb, Pc, Area) is true if Area is twice the signed area of the */ /* triangle determined by Pa, Pb and Pc. The area is positive if Pa, Pb */ /* and Pc are oriented anti-clockwise, negative if clockwise, and zero if */ /* the points are collinear. */ area(p(Xa,Ya), p(Xb,Yb), p(Xc,Yc), Area):- Area is (Xb-Xa) * (Yc-Ya) - (Xc-Xa) * (Yb-Ya). /* square(P1, P2, P3, P4) is true if the points P1, P2, P3 and P4 are the */ /* vertices going clockwise round a square. */ /* e.g. square(p(11,4), p(5,0), p(1,6), p(7,10)). */ square(p(X1,Y1), p(X2,Y2), p(X3,Y3), p(X4,Y4)):- Y2_Y1 is Y2 - Y1, X1_X2 is X1 - X2, X4 is X1 + Y2_Y1, Y4 is Y1 + X1_X2, X3 is X2 + Y2_Y1, Y3 is Y2 + X1_X2. /* perpendicular_intersection(P1, P2, P3, P4) is true if P4 is the */ /* intersection of the line through P1 and P2 with the perpendicular from */ /* P3 to this line. */ /* e.g. perpendicular_intersection(p(0,0), p(0,4), p(2,6), p(0,6)). */ perpendicular_intersection(_, P, P, P):-!. perpendicular_intersection(P1, P2, p(X3,Y3), p(X4,Y4)):- line_through_points(P1, P2, A, B, C), Temp1 is A*A + B*B, Temp1 =\= 0.0, Temp2 is B*X3 - A*Y3, X4 is (B*Temp2 - A*C) / Temp1, Y4 is (-A*Temp2 - B*C) / Temp1. /* line_through_points(P1, P2, A, B, C) is true if Ax + By + C = 0 is the */ /* equation of the line through P1 and P2. */ /* e.g. line_through_points(p(3,8), p(6,12), 4, -3, 12). */ line_through_points(p(X1,Y1), p(X2,Y2), A, B, C):- A is Y2 - Y1, B is X1 - X2, C is - A*X1 - B*Y1.