Resize and reposition a disk so its bounding box accommodates other objects

Question

This is a continuation of this post.

I have a circle or disk with some rectangles located on its circumference which have arbitrary widths, heights, and angles (by angle I mean the line connecting center of the circle to center of the rectangle). An example angle is shown for one of the rectangles with pink lines.

The center coordinate of the circle is (centerX, centerY) and its radius is r.
The widths and heights of the rectangles and angles of the lines connecting circle center to each rectangle center are available too.

How can I find a new radius r and a new centerX and CenterY for the disk so that it has maximum area and the previous bounding box of the disk (dashed green line) can fully accommodate all the rectangles?

A relaxed version of the problem

If the above problem is too hard to solve, then, how can I just modify the circle radius r (its center stays the same) so that its previous bounding box can contain all the rectangles?

Answer 1

If an axis-aligned box is centered at $(x, y)$, its width being $2 w$, height $2 h$, the distance from the box to origin is $d$, $$d = \begin{cases} \sqrt{c_x^2 + c_y^2}, & \text{if } c_x \gt 0 \text{ and } c_y \gt 0 \\ \max(0,\, c_x,\, c_y), & \text{otherwise} \\ \end{cases} \tag{1a}\label{G1a}$$ where $$c_x = \lvert x \rvert - w , \quad c_y = \lvert y \rvert - h \tag{1b}\label{G1b}$$ are the distance between the box and the $y$ axis, and the distance between the box and the $x$ axis, respectively.

If we use $$r = \sqrt{x^2 + y^2}, \quad a_x = \frac{\lvert x \rvert}{r}, \quad a_y = \frac{\lvert y \rvert}{r} \tag{1c}\label{G1c}$$ then we can solve $\eqref{G1a}$ for $r$ via $$\begin{aligned} b &= w a_x + h a_y \\ c &= h a_x - w a_y \\ r_{xy} &= b + \sqrt{d^2 - c^2} \\ r_{x} &= \frac{d + w}{a_x} \\ r_{y} &= \frac{d + h}{a_y} \\ r &= \begin{cases} r_{xy}, & d^2 \ge c^2 \text{ and } r_{xy} a_x \ge w \text{ and } r_{xy} a_y \ge h \\ r_{x}, & d^2 \ge c^2 \text{ and } r_{xy} a_x \ge w \text{ and } r_{xy} a_y \lt h \\ r_{y}, & d^2 \ge c^2 \text{ and } r_{xy} a_x \lt w \text{ and } r_{xy} a_y \ge h \\ 0, & d^2 \ge c^2 \text{ and } r_{xy} a_x \lt w \text{ and } r_{xy} a_y \lt h \\ r_{x}, & d^2 \lt c^2 \text{ and } c \gt 0 \\ r_{y}, & d^2 \lt c^2 \text{ and } c \lt 0 \\ \end{cases} \\ \end{aligned} \tag{1d}\label{G1d}$$
Here is an example Python 3 implementation, rect.py. You can import the Rect type to use in your own code via from rect import Rect, run pydoc3 rect to see the documentation of the Rect class it implements, or run the file using e.g. python3 -B rect.py to test it with randomly generated rectangles. It is licensed under the Creative Commons Zero v1.0 Universal license; so feel free to reuse as you wish.

# SPDX-License-Identifier: CC0-1.0

from math import pi, atan2, tan, sqrt, sin, cos, pi

class Rect:
    __slots__ = ( "_w", "_h", "_a", "_d", "_x", "_y" )

    def __init__(self, width=0.0, height=0.0, distance=0.0, angle=0.0):
        w = 0.5 * abs(float(width))
        h = 0.5 * abs(float(height))
        a = float(angle) * pi / 180.0
        d = abs(float(distance))

        self._w = w         # Half-width
        self._h = h         # Half-height
        self._a = a         # Angle in radians
        self._d = d         # Distance from perimeter to origin
        self._set_x_y_d(d)

    def _set_x_y_d(self, distance):
        """Internal: Recalculates self._x and self._y based on perimeter to origin distance"""
        d = abs(float(distance))
        w = self._w
        h = self._h

        dx = cos(self._a)
        dy = sin(self._a)
        ax = abs(dx)
        ay = abs(dy)

        b = w*ax + h*ay
        c = h*ax - w*ay

        dd = d*d
        cc = c*c

        if dd >= cc:
            r = b + sqrt(dd - cc)
            cx = r*ax - w
            cy = r*ay - h
            if cx > 0 and cy > 0:
                # r is correct
                pass
            elif cx > 0:
                r = (d + w) / ax
            elif cy > 0:
                r = (d + h) / ay
            else:
                r = d + sqrt(w*w + h*h)
        elif c > 0:
            r = (d + w) / ax
        else:
            r = (d + h) / ay

        self._d = d
        self._x = r*dx
        self._y = r*dy

    def _set_x_y_r(self, distance):
        """Internal: Recalculates self._x and self._y and self._d based on center to center distance"""
        r  = abs(float(distance))
        dx = cos(self._a)
        dy = sin(self._a)
        cx = r * abs(dx) - self._w
        cy = r * abs(dy) - self._h
        if cx >= 0 and cy >= 0:
            d = sqrt(cx*cx + cy*cy)
        elif cx > 0:
            d = cx
        elif cy > 0:
            d = cy
        else:
            r = 0
            d = 0

        self._d = d
        self._x = r*dx
        self._y = r*dy

    @property
    def copy(self):
        """Yields an exact duplicate of this rectangle object"""
        r = Rect()
        r._w = self._w
        r._h = self._h
        r._a = self._a
        r._d = self._d
        r._x = self._x
        r._y = self._y
        return r

    @property
    def width(self):
        """Rectangle size along X axis (width)"""
        return 2 * self._w

    @width.setter
    def width(self, width):
        self._w = 0.5 * abs(float(width))
        self._set_x_y_d(self._d)

    @property
    def height(self):
        """Rectangle size along Y axis (height)"""
        return 2 * self._h

    @height.setter
    def height(self, height):
        self._h = 0.5 * abs(float(height))
        self._set_x_y_d(self._d)

    @property
    def distance(self):
        """Rectangle distance from center of disk"""
        return self._d

    @distance.setter
    def distance(self, distance):
        self._set_x_y_d(float(distance))

    @property
    def radius(self):
        """Distance from center of rectangle to center of disk"""
        return sqrt(self._x*self._x + self._y*self._y)

    @radius.setter
    def radius(self, radius):
        self._set_x_y_r(float(radius))

    @property
    def angle(self):
        """Rectangle direction, angle in degrees"""
        return self._a * 180.0 / pi

    @angle.setter
    def angle(self, angle):
        self._a = float(angle) * pi / 180.0
        self._set_x_y_d(self._d)

    @property
    def center(self):
        """Rectangle center as a tuple"""
        return (self._x, self._y)

    @center.setter
    def center(self, p):
        x = float(p[0])
        y = float(p[1])
        self._a = atan2(y, x)
        self._set_x_y_r(sqrt(x*x + y*y))

    @property
    def min(self):
        """Rectangle corner with the smallest coordinates as a tuple"""
        return (self._x-self._w, self._y-self._h)

    @property
    def max(self):
        """Rectangle corner with the largest coordinates as a tuple"""
        return (self._x+self._w, self._y+self._h)

    @property
    def aabb(self):
        """Axis-aligned bounding box (with respect to its disk)"""
        return (self._x-self._w, self._y-self._h, self._x+self._w, self._y+self._h)

    @property
    def x(self):
        """Center X coordinate (with respect to its disk)"""
        return self._x

    @property
    def y(self):
        """Center Y coordinate (with respect to its disk)"""
        return self._y

    @property
    def xmin(self):
        """Minimum X coordinate (left edge)"""
        return self._x - self._w

    @property
    def ymin(self):
        """Minimum Y coordinate (top edge)"""
        return self._y - self._h

    @property
    def xmax(self):
        """Maximum X coordinate (right edge)"""
        return self._x + self._w

    @property
    def ymax(self):
        """Maximum Y coordinate (bottom edge)"""
        return self._y + self._h


if __name__ == '__main__':
    from random import Random

    u = Random().uniform

    rerrmax = 0.0
    derrmax = 0.0

    print("Testing Rect() radius/distance correctness.  Press Ctrl+C to stop.")

    n = 0
    while True:
        try:
            w = u(0.1, 1.0)
            h = u(0.1, 1.0)
            d = u(0.0, 1.0)
            a = u(0.0, 360.0)

            rect = Rect(w, h, d, a)
            r = rect.radius
            rect.radius = r
            R = rect.radius
            D = rect.distance

            rerr = abs(R - r)
            derr = abs(D - d)
            if rerr > rerrmax or derr > derrmax:
                rerrmax = max(rerrmax, rerr)
                derrmax = max(derrmax, derr)
                print("%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s" % ("w = %.7f" % w, "h = %.7f" % h, "a = %6.2f" % a, "d = %.7f" % d, "D = %.7f" % D, "D-d = %+.7g" % (D-d), "r = %.7f" % r, "R = %.7f" % R, "R-r = %+.7g" % (R-r)))

        except KeyboardInterrupt:
            print("After %d tests, maximum r error is %.16f and maximum d error is %.16f" % (n, rerrmax, derrmax))
            break

        n += 1
        if (n & 65535) == 0:
            print("After %d tests, maximum r error is %.16f and maximum d error is %.16f" % (n, rerrmax, derrmax))

Another file, disk.py, implements the disk (with associated rectangles) logic:

# SPDX-License-Identifier: CC0-1.0

from rect import Rect

class Disk:
    slots = ('_x', '_y', '_d', '_r')

    def __init__(self, *rects, x=0.0, y=0.0, radius=0.0):
        self._x = float(x)
        self._y = float(y)
        self._d = abs(float(radius))
        self._r = []

        if len(rects) > 0:
            for r in rects:
                if isinstance(r, Rect):
                    # Choose minimum distance among rectangles
                    if self._d == 0 or self._d > r.distance:
                        self._d = r.distance
                    self._r.append(r)
                else:
                    raise TypeError("Expected a Rect, got a %s." % str(type(r)))

    @property
    def x(self):
        """Center of the disk, X coordinate"""
        return self._x

    @x.setter
    def x(self, x):
        self._x = float(x)

    @property
    def y(self):
        """Center of the disk, Y coordinate"""
        return self._y

    @y.setter
    def y(self, y):
        self._y = float(y)

    @property
    def center(self):
        """Center of the disk"""
        return (self._x, self._y)

    @center.setter
    def center(self, c):
        self._x = float(c[0])
        self._y = float(c[1])

    @property
    def aabb(self):
        """Axis-aligned bounding box"""
        abb = (-self._d, -self._d, self._d, self._d)
        for r in self._r:
            rbb = r.aabb
            abb = (min(abb[0], rbb[0]), min(abb[1], rbb[1]), max(abb[2], rbb[2]), max(abb[3], rbb[3]))
        return (abb[0]+self._x, abb[1]+self._y, abb[2]+self._x, abb[3]+self._y)

    @property
    def size(self):
        """Width and height of the axis-aligned bounding box"""
        bb = self.aabb
        return (bb[2]-bb[0], bb[3]-bb[1])

    @property
    def width(self):
        """Width of the axis-aligned bounding box"""
        bb = self.aabb
        return bb[2] - bb[0]

    @property
    def height(self):
        """Height of the axis-aligned bounding box"""
        bb = self.aabb
        return bb[3] - bb[1]

    @property
    def copy(self):
        """Independent copy of the disk and related rectangles"""
        disk = Disk()
        disk._x = self._x
        disk._y = self._y
        disk._d = self._d
        disk._r = [ r.copy for r in self._r ]
        return disk

    def size_at_radius(self, radius):
        """Axis-aligned bounding box size, if all rectangles moved to a specific distance"""
        d = abs(float(radius))
        dbb = (-d, -d, d, d)
        for r in self._r:
            c = r.copy
            c.distance = d
            rbb = c.aabb
            dbb = (min(dbb[0], rbb[0]), min(dbb[1], rbb[1]), max(dbb[2], rbb[2]), max(dbb[3], rbb[3]))
        return (dbb[2]-dbb[0], dbb[3]-dbb[1])

    def copy_at(self, radius, x=None, y=None):
        """Independent copy of the disk and rectangles, with rectangles moved to a specific distance"""
        disk = Disk()
        d = abs(float(radius))
        disk._d = d
        if x is None:
            disk._x = self._x
        else:
            disk._x = float(x)
        if y is None:
            disk._y = self._y
        else:
            disk._y = float(y)

        for r in self._r:
            c = r.copy
            c.distance = d
            disk._r.append(c)

        return disk

    @property
    def radius(self):
        """Disk radius"""
        return self._d

    @radius.setter
    def radius(self, radius=-1.0):
        if radius >= 0.0:
            d = float(radius)
        else:
            d = self._d
        for r in self._r:
            r.distance = d

    @property
    def rects(self):
        """Immutable list (tuple) of rectangles"""
        return tuple(self._r)

    def resized_to(self, width, height):
        size = ( abs(float(width)), abs(float(height)) )

        if len(self._r) < 1:
            disk = Disk()
            disk._x = self._x
            disk._y = self._y
            disk._d = 0.5*min(w, h)
            return disk
        else:
            d = self._r[0].distance
            dmin = d
            dmax = d
            for r in self._r:
                d = r.distance
                dmin = min(dmin, d)
                dmax = max(dmax, d)

        if dmin < 0:
            dmin = 0
        if dmax <= dmin:
            dmax = dmin + 1.0

        minsize = self.size_at_radius(dmin)
        maxsize = self.size_at_radius(dmax)

        while (dmin > 0) and (minsize[0] > size[0] or minsize[1] > size[1]):
            dmin = 0.5 * dmin
            minsize = self.size_at_radius(dmin)

        if (minsize[0] > size[0] or minsize[1] > size[1]):
            return self.copy_at(radius=0)

        while maxsize[0] < size[0] or maxsize[1] < size[1]:
            dmax = 2.0 * dmax
            maxsize = self.size_at_radius(dmax)

        while True:
            d = 0.5 * dmin + 0.5 * dmax
            if d == dmin or d == dmax:
                break

            dsize = self.size_at_radius(d)
            if dsize[0] > size[0] or dsize[1] > size[1]:
                dmax = d
            elif dsize[0] < size[0] or dsize[1] < size[1]:
                dmin = d
            else: break

        return self.copy_at(radius=d)

if __name__ == '__main__':
    def describe(disk, name):
        print("%s:" % name)
        print("  Width: %.3f" % disk.width)
        print(" Height: %.3f" % disk.height)
        print(" Radius: %.3f" % disk.radius)
        for r in disk.rects:
            print("   Rect: (%.3f, %.3f) - (%.3f, %.3f), distance %.3f" % (*r.aabb, r.distance))

    disk = Disk( Rect(2.0, 1.0,  4.0,  15),
                 Rect(4.0, 2.0,  4.0,  85),
                 Rect(1.0, 1.0,  4.0, 140),
                 Rect(2.0, 3.0,  4.0, 300) )
    describe(disk, "Initial disk")

    other = disk.resized_to(8.0, 6.0)
    describe(other, "If bounding box set to 8.0 by 6.0")

    other = disk.resized_to(30.0, 24.0)
    describe(other, "If bounding box set to 30.0 by 24.0")

If run, it constructs a disk of radius $4$ with four rectangles around it. It then uses the resized_to() method of the Disk class to resize the disk so that the resulting bounding box matches the specified width and/or height.

The size_at_radius() method computes the bounding box size, if the disk was resized to the specified radius and all rectangles moved correspondingly. The copy_at() method returns the same as a new Disk object, without modifying the original.

Both functions use the Rect.distance setter method to compute the new location of each rectangle, with the calculations in Rect._set_x_y_d in the first file; which in turn is a straightforward implementation of $\eqref{G1d}$.

The resized_to() method uses the size_at_radius() method to bracket the radius range, then a binary search to find the radius that yields the best bounding box size.

I have not proven that the disk bounding box size is a smooth enough function for a binary search to be reasonable (converge to either the global extremum, or to a local one so near to the global one that the difference is within rounding error), but the behaviour of $\eqref{G1a}$ makes me believe so.

The code is not production quality; it probably contains typos and other bugs (none intentionally!), so anyone intending to use this should really test their code well with different datasets before relying on it. However, it should clearly show how $\eqref{G1a}$ is inverted (and running rect.py verifies this with lots of random rectangles), and how a binary search can be used to find a disk radius that yields a desired bounding box size.

It may not behave sensibly when the origin (center of disk) is within a rectangle, though.

At minimum, consider limiting the binary search to meaningful differences in disk radiuses. As it is, the binary search only ends if a desired bounding box size is found, or when the limit of precision is reached.

The bracketing part extends the minimum and maximum radius, until it covers the desired bounding box size. Note that when the maximum has to be extended, we could set the minimum radius to the old maximum radius, so that we move the range rather than scale it; and vice versa for when extending the minimum. So, there are definitely lots of optimization opportunities.
It is just example code, after all.