Minimum Rectangle

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.

LPA Index     Home Page

Valid HTML 4.0 Strict