How to evenly distribute points around the perimeter of a rectangle

Question

I'm looking for a way to distribute points along a portion of the perimeter of a rectangle. These points needs to be evenly far away from each other.

I have a rectangular (usually square) bounds, and 2 points (ps and pe) along that perimeter that mark the allowed range for points. Here I've marked the allowed range in red:

I need to place n points along that segment (usually 1-3). These points need to be evenly spaced at a distance d. So the distance between n0..n1 and n1..n2 etc should all be d. The boundary points also count for the purposes of distribution, so the distance between the first and last points and ps/pe should be d as well.

This seemed like a straightforward task going in but I quickly realised the naive method doesn't work here. Taking the length of the segment and dividing by n+1 doesn't factor in corners. For example: n = 1, puts the point too close to pe:

My math is pretty rusty (day job doesn't require much of it usually), but I've tried several different approaches and none have quite worked out. I was able to solve for n = 1 using vectors, by finding the midpoint between ps and pe, finding a perpendicular vector and then intersecting that with the segment, like below. I have no idea how to make this approach work if n is something else though, or even if it can be done.

One final note, if completely even distribution is impracticable then a good enough approximation is fine. Ideally the approximation is off by roughly the same amount throughout the range (instead of say, worse on the edges).

Warman

Interesting problem! Does the range always contain exactly one corner of the rectangle, or can it be also 2 or 3? – Rotberg

It could be more, for example with pe and ps being the SE and SW corners and the allowed segment being all but the S edge. However, due to how these points are used any situation that would include more than 1 corner could probably just use a mirrored 1-corner solution. – Warman

Do you mean d as a geometric distance or Mahalanobis distance? – Moussorgsky

Geometric distance. – Warman

Answer 1

I am going to suggest an algorithm, but since the derivation is a bit mathematically messy, I did not have enough time to think it through carefully and to check it carefully for correctness. Plus I might have included some redundant checks, which if one proves some proper inequalities may become redundant and one may prove the existence of a solution may always exist under reasonable assumptions. I believe the idea is correct, but I might have made some mistakes writing this thing up, so be careful.

Because according to your comment, having only one corner inside the segment along the boundary of the square is enough to solve the rest of the cases due to symmetry, I will focus on the one corner case.

Your polygonal segment with one 90 degree corner is divided into a pair of perpendicular straight line segments, the first one of length l1 and the second of length l2. These two lengths are given to you. You also want to add on the polygonal segment, which is of total length l1 + l2, a given number of n points so that the euclidean straight line distance between any two consecutive points is the same. Call that unknown distance d. When you do that you are going to end up with n1 full segments of unknown length d on l1 and n2 full segments of unknown length d on l2 so that

n1 + n2 = n

In general, you will also end up with an extra segment of length d1 <= d on l1 with one end at the 90 degree corner. Analogously, you will also have an extra segment of length d2 <= d on l2 with one end at the 90 degree corner. Thus, the two segments d1 and d2 share a common end and are perpendicular, so they form a right-angled triangle. According to Pythagoras' theorem, these two segments satisfy the equation

d1^2 + d2^2 = d^2

If we combine all the equations and information up to now, we obtain a system of equations and restrictions which are:

n1*d + d1 = l1
n2*d + d2 = l2
d1^2 + d2^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers

where the variables d, d1, d2, n1, n2 are unknown while l1, l2, n are given. From the first two equations, you can express d1 and d2 and substitute in the third equation:

d1 = l1 - n1*d
d2 = l2 - n2*d
(l1 - n1*d)^2 + (l2 - n2*d)^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers

In the special case, when one wants to add only one point, i.e. n = 1, one has either n1 = n = 1 or n2 = n = 1 depending on whether l1 > l2 or l1 <= l2 respectively. Say l1 > l2. Then n1 = n = 1 and n2 = 0 so

d1 = l1 - d
d2 = l2
(l1 - d)^2 + l2^2 = d^2

Expand the equation, simplify it and solve for d:

l1^2 - 2*l1*d + d^2 + l2^2 = d^2
l1^2 + l2^2 - 2*l1*d = 0
d = (l1^2 + l2^2) / (2*l1)

Next, let us go back to the general case. You have to solve the system

(l1 - n1*d)^2 + (l2 - n2*d)^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers

where the variables d, n1, n2 are unknown while l1, l2, n are given. Expand the first equation:

l1^2  -  2 * l1 * n1 * d  +  n1^2 * d^2  +  l2^2  -  2 * l2 * n2 * d  +  n2^2 * d^2 = d^2
n1 + n2 = n
n1 and n2 are non-negative integers

and group the terms together

(n1^2 + n2^2 - 1) * d^2  - 2 * (l1*n1 + l2*n2) * d  +  (l1^2 + l2^2) = 0
n1 + n2 = n
n1 and n2 are non-negative integers

The first equation is a quadratic equation in d

(n1^2 + n2^2 - 1) * d^2  - 2 * (l1*n1 + l2*n2) * d  +  (l1^2 + l2^2) = 0

By the quadratic formula you expect two solutions (in general, you choose whichever is positive. If both are positive and d < l1 and d < l2, you may have two solutions):

d = ( (l1*n1 + l2*n2) +- sqrt( (l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) ) ) / (n1^2 + n2^2 - 1)
n1 + n2 = n
n1 and n2 are non-negative integers

Now, if you can find appropriate n1 and n2, you can calculate the necessary d using the quadratic formula above. For solutions to exist, the expression in the square root has to be non-negative, so you have the inequality restriciton

d = ( (l1*n1 + l2*n2) +- sqrt( (l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) ) ) / (n1^2 + n2^2 - 1)
(l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) >= 0
n1 + n2 = n
n1 and n2 are non-negative integers

Simplify the inequality espression:

(l1*n1 + l2*n2)^2 - (l1^2 + l2^2)*(n1^2 + n2^2 - 1) = (l1^2 + l2^2) - (l1*n2 - l2*n1)^2

which brings us to the following system

d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
(l1^2 + l2^2) - (l1*n2 - l2*n1)^2 >= 0
n1 + n2 = n
n1 and n2 are non-negative integers

Factorizing the inequality,

d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
(sqrt(l1^2 + l2^2) - l1*n2 + l2*n1) * (sqrt(l1^2 + l2^2) + l1*n2 - l2*n1) >= 0
n1 + n2 = n
n1 and n2 are non-negative integers

So you have two cases for these restrictions:

Case 1:

d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
sqrt(l1^2 + l2^2) - l1*n2 + l2*n1 >= 0 
sqrt(l1^2 + l2^2) + l1*n2 - l2*n1 >= 0
n1 + n2 = n
n1 and n2 are positive integers

or Case 2:

d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
sqrt(l1^2 + l2^2) - l1*n2 + l2*n1 <= 0 
sqrt(l1^2 + l2^2) + l1*n2 - l2*n1 <= 0
n1 + n2 = n
n1 and n2 are positive integers

we focus on case 1 and see that case 2 is not possible. Start by expressing n2 = n - n1, then substitute it in each of the two inequalities and isolate n1 on one side of each inequality. This procedure yields:

Case1:

d = ( (l1*n1 + l2*n2) +- sqrt( (l1^2 + l2^2) - (l1*n2 - l2*n1)^2 ) ) / (n1^2 + n2^2 - 1)
( l1*n - sqrt(l1^2 + l2^2) ) / (l1 + l2) <=  n1  <= ( l1*n + sqrt(l1^2 + l2^2) ) / (l1 + l2) 
n1 + n2 = n
n1 and n2 are positive integers

One can see that case 2 inverts the inequalities, which is impossible because the left side is less than the right one.

So the algorithm could be something like this:

function d = find_d(l1, l2, n)
{
   if n = 1 and l1 > l2 { 
      return d = (l1^2 + l2^2) / (2*l1)
   } else if n = 1 and l1 <= l2 {
      return d = (l1^2 + l2^2) / (2*l2)
   }
   for integer n1 >= 0 starting from floor( ( l1*n - sqrt(l1^2 + l2^2) ) / (l1 + l2) ) to floor( ( l1*n + sqrt(l1^2 + l2^2) ) / (l1 + l2) ) + 1 
   {
      n2 = n - n1
      D = (l1^2 + l2^2) - (l1*n2 - l2*n1)^2
      if D >= 0
      {
         d1 = ( (l1*n1 + l2*n2) - sqrt( D ) ) / (n1^2 + n2^2 - 1)  
         d2 = ( (l1*n1 + l2*n2) + sqrt( D ) ) / (n1^2 + n2^2 - 1) 
         if 0 < d1 < max(l1, l2) {       
            return d1
         } else if  0 < d2 < max(l1, l2) {
            return d2
         } else {
            return "could not find a solution"
         }
      }
   }
}

Spermary

I've gotten this working great for the 1-corner case, seems to be giving accurate results every time. You referenced how I said in one of my comments it can be mirrored to solve for >1 corners. I've just been told the specs were wrong and that isn't actually the case (although it will only ever be 2 corners maximum). I'll mark this as the answer since it does solve the initial question, but any tips on how to expand it to work for that case would be greatly appreciated! – Warman

@Warman I am glad it helps, it is an interesting problem. As you can see one corner already makes things convoluted, I did not try to extend to two corners. Also, I checked the quadratic formulas and some of the inequalities but there might be more rigorous inequalities that rule out some of the cases. I just did not want to go into that too much. After I am done with some work-related tasks, I will try to think about two corners. – Spermary

@Warman I added another answer with a preliminary version of an algorithm for two corners. It's in python and I made it draw some pictures as a sanity check. I have not added the math behind the algorithm, but if you want me to, I can add some explanation over the weekend. – Spermary

Answer 2

This is a preliminary version, so I suggest to approach with some caution. I did not have enough time to check the algorithm whether there might be some close to degenerate cases, for which one may have to add somewhere few short loops with if statements. However, in general, this will probably work almost always. I am just posting a python implementation but when I find a bit more time and if you want me to, I can write down the math behind this algorithm. Some of the ideas from this algorithm can simplify the previous one for one corner.

import numpy as np
import math
import matplotlib.pyplot as plt


def sq_root(x, m, K):
  return math.sqrt(x**2 - (K - m*x)**2)

def f(x, n, L):
  return sq_root(x, n[0], L[0]) + sq_root(x, n[2], L[2]) + n[1]*x - L[1]

def df(x, n, L):
  return ((1-n[0]**2)*x + L[0]*n[0])/sq_root(x, n[0], L[0]) + ((1-n[2]**2)*x + L[2]*n[2])/sq_root(x, n[2], L[2]) + n[1]

#Solving the nonlinear equation for d by using Newton's method:
def solve_f(n, L):
  x = sum(L) / (sum(n) + 2)
  y = f(x, n, L)
  while abs(y) >= 0.0000001:
    x = x - y / df(x, n, L)
    y = f(x, n, L) 
  return x - y / df(x, n, L)

def find_n(L, N):
  x0 = sum(L) / (N + 1)
  # x <= x0
  n = np.array([0,0,0])
  n[0] = math.floor(L[0]/x0)
  n[2] = math.floor(L[2]/x0)
  n[1] = N - n[0] - n[2] - 1
  return n

def find_d(L, N):
  if N==1:
    d2 = (L[2]**2 + L[1]**2 - L[0]**2)/(2*L[1])
    return math.sqrt(L[0]**2 + d2**2), np.array([0,0,0])
  n = find_n(L, N)
  return solve_f(n, L), n

def find_the_points(L, N):
  d, n = find_d(L, N)
  d2 = math.sqrt(d**2 - (L[0]-n[0]*d)**2)
  #d3 = math.sqrt(d**2 - (L[2]-n[2]*d)**2)
  p = np.zeros((sum(n) + 3, 2))
  p[ 0] = np.array([0, L[1]-L[0]])
  p[-1] = np.array([L[1], L[1]-L[2]])
  e_x = np.array([1,0])
  e_y = np.array([0,1])
  corner = np.array([0,L[1]]) 
  for i in range(n[0]):
    p[i+1] = p[0] + (i+1)*d*e_y
  for i in range(n[1]+1):
    p[n[0]+i+1] = corner + (d2 + i*d)*e_x
  for i in range(n[2]):
    p[-(2+i)] = p[-1] + (i+1)*d*e_y
  return p, d, n 

'''
Test example:
'''

# lengths of the three straight segments along the edges of a square of edge_length L2:
L1 = 5
L2 = 7
L3 = 3
L = np.array([L1, L2, L3]) 

N = 7
# N = number of points to be added
# If there are two corners then number of segments aligned with edges of square is N - 1
# total number of equidistant segments is N + 1
# n = n[0], n[1], n[2] represents the number of segments aligned with each 
# striaght segment from the rectangular polyline along square's boundary


points, d, n = find_the_points(L, N)
print(points)
print(d)
print(n)

plt.figure()
plt.plot(points[:,0], points[:,1]) 
for j in range(points.shape[0]):
  plt.plot(points[j,0], points[j,1], 'ro')
axx = plt.gca()
axx.set_aspect('equal')
plt.show() # if you need...

Spermary

@DanSc I spotted a couple of bugs and fixed them. – Spermary

I've incorporated this into my previous code and from what I can tell it's accurate for all the cases I'm testing against right now. I'm sure I can iron out any bugs if I do find them. Thank you for taking the time to make this updated answer, it's very insightful. – Warman

@Warman Thank you for the interesting problem! I also tested a bit more and it works accuratly (after I corrected 'y >= ' with 'abs(y) >= ' in the while loop of function 'solve_f(n, L)', which was kind of important). The only thing that may (or may not) happen is that if the calculations for n[0], n[1], n[2] do not work, one may need to set n[0] = n[0]+1 and adjust the rest of the n's respectively. – Spermary

Answer 3

Pythagoras to the rescue.

Considering the case where n = 2 with one corner¹. Rest of the cases will be likely similar just with more edge cases and possible configurations².

Let us name the length of the vertical segment i and length of the horizontal segment j.
We are looking for a points X, Y, so that distance dist(ps,x) = dist(x,y) = dist(y,pe)

First, let's assume that the corner lies between points X and Y. In that case, we are looking for the solution of this equation: x² = (i - x)² + (j - x)² where i and j are known.

If the abovementioned equation has no solution, it means that the corner has to lie either between ps and X or between Y and pe. I will only cover the case of ps and X since the other one is symetrical:
x² = i² + (j-2x)² is the equation for this case.

For more points and corners, there will just be more possible configurations of points and corners. However since the distance should be equal and the lengths are known, a series of quadratic equations will be most likely enough. The case for 5+ points and three corners will be a wee bit iffy though.

---

¹ The case for n=1 will be nearly the same as the asymmetric variant of the n=2 that I will cover.
²Which will make them rather ugly.

Dejection

Answer 4

The points on any one of the sides are to be equidistant according to the problem. Just guess the distance and use binary search to find the solution to any degree of accuracy. There's a slightly tricky part in determining how many points lie on each rectangle side. (Hopefully it's just "slightly." :)

Sworn