June 3, 2022

Intersecting rectangles

Main problem and motivation»

At work, we had the task of calculating the area of a partially overlapping rectangular image. Since I was busy with a larger overhaul of the frontend and the problem was mostly in the backend, my colleague was assigned to work on it.

However, I was fascinated by the problem and decided to work on it in my spare time. It turns out that this is a “solved problem,” and a good description of the algorithm can be found in Sinan Oz’s blog post intersection polygon of convex polygons, which I’ll briefly summarize here for context.

His main algorithm consists of the following steps and works not only with rectangles, but with all convex polygons:

Create an empty polygon point list
Add all points of polygon 1 that are within polygon 2
Add all points of polygon 2 that are within polygon 1
Add all intersection points of the sides of polygon 1 with the sides of polygon 2
Sort the point list counter-clockwise

To calculate the size of the resulting polygon, we used the Shoelace formula.

Sub problems»

I was fascinated by the solution, the clear divison into sub-problems, which make the main algoritm look much simpler. So I decided to create interactive visualizations with Svelte of the steps as well as the result, while translating the C Code to Javascript. My colleague was using PHP for the Backend code and had to come up with relevant unit tests to make sure, everything’s working. With Svelte, I had immediate visual feedback during implementation. Also I replaced the for loops with map/reduce wherever possible.

Is a point inside a polygon?»

To figure out if a point is inside a given convex polygon (or rectangle), Sinan used a Ray Casting Algorithm, most probably based on the ‘crossing number algorithm’ of M. Shimrat (1962) and mentions its relation to the Jordan Curve Theorem.

The used adaption checks, if a semi-infinite horizontal ray from the point to check intersects an odd number of times with the polygon. It’s really little code when using reduce.

You can drag and drop either the points or the corners of the four-sided figure to play around.

const isSegmentIntersectingRay = (p, q, r) =>
  (q.y > p.y) != (r.y > p.y) &&
  (p.x < (r.x - q.x) * (p.y - q.y) / (r.y - q.y) + q.x)
;

export const isPointInsidePoly = (p, points) =>
  points.reduce(
    (oddNumOfXings, q, i) =>
      isSegmentIntersectingRay(p, q, points[(i + 1) % points.length]) !== oddNumOfXings
    ,
    false
  )
;

Find the Intersection Point of two lines.»

Well this is something everyone remembers from school, right? Given are two lines and you want to know its intersection.

But because we are interested in the intersection point within certain line segments, we also need to check if the point is inside these bounds. It has to be inside the bounding boxes of both line segments.

const isInLineBounds = (q, r, p) =>
  isLessOrEqual(Math.min(q.x, r.x), p.x) &&
  isMoreOrEqual(Math.max(q.x, r.x), p.x) &&
  isLessOrEqual(Math.min(q.y, r.y), p.y) &&
  isMoreOrEqual(Math.max(q.y, r.y), p.y)
;

const calculateLineCoefficients = (p, q) => {
  const a = q.y - p.y;
  const b = p.x - q.x;
  const c = a * p.x + b * p.y;
  return { a, b, c };
};

export const getIntersectionPoint = (fa, fb, ga, gb) => {
  const lcf = calculateLineCoefficients(fa, fb);
  const lcg = calculateLineCoefficients(ga, gb);

  const det = lcf.a * lcg.b - lcg.a * lcf.b;
  if (det === 0) return null; // lines are parallel

  const ip = point(
    (lcg.b * lcf.c - lcf.b * lcg.c) / det,
    (lcf.a * lcg.c - lcg.a * lcf.c) / det
  );

  if (isInLineBounds(fa, fb, ip) && isInLineBounds(ga, gb, ip)) return ip;

  return null; // intersection is at least out of one bounding box
};

Get all Intersection Points of a line
with a Rectangle.»

I very much agree with Sinan, that it is so easy to just check the line segment with each edge of the rectangle, even easier with reduce().

const getIntersectionPoints = (fa, fb, points) =>
  points.reduce((ips, p, i) => {
    const q = points[(i + 1) % points.length];
    const ip = getIntersectionPoint(fa, fb, p,q);
    return ip !== null ? [...ips, ip] : ips;
  }, [])
;

Order Point List clockwise.»

To make it possible to draw the intersection polygon we need to order the points clockwise. For this reason Sinan first determines the mid point between all polygon points, and then sorts all of them according to the angle of the lines between the mid point and each polygon point.

You can drag and drop all the points in the visualization and see the mid point move. If you move some points past others, you might notice the order number changing accordingly.

const orderBy = (arr, orderCallback) =>
  arr.sort((a, b) => orderCallback(a) - orderCallback(b))
;

export const getPolygonMidPoint = (points) => {
  const sum = points.reduce(plus, point(0,0));
  return point(
    sum.x/points.length,
    sum.y/points.length
  );
};

export const orderClockwise = (points) => {
  const m = getPolygonMidPoint(points);
  return orderBy(points, (p) => Math.atan2(p.y - m.y, p.x - m.x));
};

Add points without duplicates.»

I’ve also translated Sinans addPoints method using reduce(). It also eliminates all duplicte points.

const addPoint = (points, p) =>
  points.find((q) => isEqual(q.x, p.x) && isEqual(q.y, p.y)) ? points : [...points, p]
;

const addPoints = (points, ps) =>
  ps.reduce(addPoint, points)
;

Calculate Area of a Polygon.»

Now to a function not part of Sinans Blog post. We need to calculate the area of a Polygon (Rectangle). We use the so called Shoelace formula also known as Gauss’s area formula. Here’s a short version using, your guessed right, reduce().

export const calcPolygonArea = (points) =>
  Math.abs(points.reduce((area, p, i) => {
    const q = points[(i + 1) % points.length];
    return area + (p.x * q.y * 0.5) - (q.x * p.y * 0.5);
  }, 0))
;

Everything together»

I did really enjoy the clean code and speaking function used in that algorithm, so Sinan should be happy. For my part, I hope the reader of my blog post did enjoy the Svelte based, interactive Illustrations.

Using all the part together, we can find the new polygon from two rectangles.

export const getIntersectionOfPolygons = (poly1, poly2) => {
  let clippedCorners = [];

  // add  the corners of poly1 which are inside poly2
  clippedCorners = poly1.reduce((acc, p) => {
    return isPointInsidePoly(p, poly2) ? addPoint(acc, p) : acc;
  }, clippedCorners);

  // add the corners of poly2 which are inside poly1
  clippedCorners = poly2.reduce((acc, p) => {
    return isPointInsidePoly(p, poly1) ? addPoint(acc, p) : acc;
  }, clippedCorners);

  // add  the intersection points
  clippedCorners = poly1.reduce((acc, p, i) => {
    const q = poly1[(i + 1) % poly1.length];
    const ips = getIntersectionPoints(p, q, poly2);
    return addPoints(acc, ips);
  }, clippedCorners);


  return orderClockwise(clippedCorners);
};

Simple Svelte visualization.»

In this visualization, two four-sided figures can overlap each other.

You can move each corner by drag and drop.

However, since the shapes can no longer be reliably considered as rectangles or convex polygons, the calculation of the intersection area may be less accurate.

This was easily put together with my existing components (DragPoint) and SVG path elements.

Forced rectangles»

The next visualization offers only two rectangles. I built a Rectangle component which offers three handles for manipulation: one for translation, one for rotation and one for scaling.

Afterword»

I admit, that most of this content already existed in the mentioned blog post, but it was fun to do the Svelte components to bring it to life.

I hope you enjoyed reading and feel free to provide feedback or suggestions.