Find the bounding box of a circle segment.

Question

You have three points. A,B and C. They define a circle segment that starts at A, goes through B and ends at C. Find the smallest bounding box that encompases the circle segment.

Here is a picture: https://docs.google.com/drawings/d/14YwCO0UeMzu-rTLmqULg7HDU5XPYZvSBQdyfu0l71Fs/edit

I started creating thee equatios for the circle:

(x_0-k)^2 + (y_0-h)^2 == r^2
(x_1-k)^2 + (y_1-h)^2 == r^2
(x_2-k)^2 + (y_2-h)^2 == r^2

Here A = (x_0,y_0) and so on. k and h is the origin of the circle.

Solving this by hand quickly showed ineffective so I turned to the computer for CAS. It manages to solve the equations but the answer is really verbose. Solving for r, eliminating k and h gives an answer with several hundred terms. And this answer only gives the whole circle. There is still the problem of deciding what part of the circle we want to create a bbox around. Is there any smarter way to do this?

Any help appreciated.

Answer 1

You may find a detailed algorithm at the Drexel Math Forum, due to Vladimir Zajic. It is a messy but elementary calculation. You will first have to compute the center of the circle containing the arc, and the start and end angles, a calculation that can be found all over the web, e.g., at the MathWorld page on Circle.

Answer 2

I managed to do it. The first phase is to calculate the circle centre and the angles, just as Joseph said. Then you have to calculate the angles between the vector (1,0) and the (x_n-cx, y_n-cy) where n is [0,1,2]. The first and last angle together with the circle define a sector. Now, calculate the four compass points on the circle. North, East, South and West. A key insight is that the only points that can intersect with the bbox are the arc end, start and theese four "compass ponits." Filter out those of these six points, start, end, north, east, south, west that are outside the sector defined previously. Afterwards, take the extreme points from each point and create a bbox from them.

EDIT: I ended up using a CAS system for solving the equations giving cx and cy since it's messy and long, but conceptually easy math.

Here is the Haskell code that does this:

import Data.Fixed

type Point = [Double]
type BBox = (Point, Point)
type Angle = Double


angle :: Point -> Point -> Angle
angle [cx,cy] [x,y] = atan2 (x-cx) (y-cy)

-- insideSector takes three angles, the start, some intermideate and the end angle.
-- It takes the centre of the circle and a point. It returns true
-- if the point is inside the sector described by the angles and the centre.
insideSector :: Angle -> Angle -> Angle -> Point -> Point -> Bool
insideSector a b c centre point =
  let ang = angle centre point
  in  insideDomain (2*pi) a b c ang

-- D is a circular domain starting at 0 and ending at m.
-- A subdomain E is given by the values btween start and end, that goes trough point
-- This function checks if the point check is inside or not E
insideDomain :: Double -> Double -> Double -> Double -> Double -> Bool
insideDomain m start point end check =
  let span = mod' (end - start) m
      dist = mod' (point - start) m
      dist2 = mod' (check - start) m
  in (dist <= span) == (dist2 <= span)


-- To understand calculation of bboxes for arc segments we'll
-- have to take a look at geometry.
-- First, and arc is defined by three points, A, B and C.
-- The arc starts at A, goes trough B and ends in C.
-- In a 2D world there exists only one circle that intersects
-- every combination of 3 points.
-- To find this circle, we do some math.
--
-- (x_0-cx)^2 + (y_0-cy)^2 == r^2 (1)
-- (x_1-cx)^2 + (y_1-cy)^2 == r^2 (2)
-- (x_2-cx)^2 + (y_2-cy)^2 == r^2 (3)
--
-- Here (x_0,y_0) = A and so on.
--
-- Theese equations can be solved for cx and cy
-- eliminating r and cy, and r and cx respectively.
-- This was done with a Computer Algebra System
-- to reduse human error and speed up the process.
--
-- We now have an equations of the centre of the circle
-- given three coordinates.
--
-- Next, a key insight is that maximum six points of an arc
-- can tangent the bounding box. Namely the beginning and end
-- of the arc and the four points furthest to each compass direction
-- (north, east, south, west).
-- 
-- The four 'compass points' are calculated given the centre and
-- the radius. The start and beginnig are already supplied.
--
-- Next, we calculate the angle between the vector (1,0) and
-- each point. Afterwards we check which of those points which
-- are inside the sector given by the circle centre, start and
-- end points. All points outside this sector is filtered out.
-- Finally, a bounding box for the remaining points are calculated.

--calcBbox :: Point -> Point -> Point -> [Point]
calcBbox a b c
  | a == b = undefined
  | otherwise =
  let [x0, y0] = a
      [x1, y1] = b
      [x2, y2] = c
      -- (cx,cy) centre of circle
      cx = -1/2*((y0 - y1)*y2^2 - y0*y1^2 - (x1^2 - x2^2)*y0 + (x0^2 - x2^2 + y0^2)*y1 - (x0^2 - x1^2 + y0^2 - y1^2)*y2)/((x1 - x2)*y0 - (x0 - x2)*y1 + (x0 - x1)*y2)
      cy = 1/2*((x1 - x2)*y0^2 - (x0 - x2)*y1^2 + (x0 - x1)*x2^2 + (x0 - x1)*y2^2 + x0^2*x1 - x0*x1^2 - (x0^2 - x1^2)*x2)/((x1 - x2)*y0 - (x0 - x2)*y1 + (x0 - x1)*y2)
      centre = [cx,cy]
      r = sqrt ((x0-cx)^2 + (y0-cy)^2)
      -- an angle between the vector (1,0) and (xn,yn)
      a0 = angle centre a
      a1 = angle centre b
      a2 = angle centre c
      -- the only points that can touch the bounding box are the start, end
      -- and the compass points on the circle: north, east, south, west.
      criticalPoints = a:c:[[cx+r,cy],[cx-r,cy],[cx,cy+r],[cx,cy-r]]
      -- we now have to check which points are inside the sector of the
      -- arc
      validPoints = filter (insideSector a0 a1 a2 [cx,cy]) criticalPoints
  in bboxFromPoints validPoints

bboxFromPoints ([a,b]:rest) =
    foldr (\[x,y] [[minx,miny],[maxx,maxy]] -> [[minimum [x,minx],minimum [y,miny]]
                                               ,[maximum [x,maxx],maximum [y,maxx]]])
          [[a,b],[a,b]] rest


metaBbox :: [BBox] -> BBox
metaBbox =
    foldr1 (\([a,b],[c,d]) ([e,f],[g,h]) -> ([minimum [a,e],minimum [b,f]]
                                            ,[maximum [c,g],maximum [d,h]]))