{-# LANGUAGE CPP #-}
{- |
Copyright : Copyright (C) 2012 Joachim Breitner
License : BSD3
Maintainer : Joachim Breitner
Stability : stable
Portability: portable
-}
module Optimisation.CirclePacking (
packCircles
-- * Example
-- $example
) where
#if defined(__GLASGOW_HASKELL__ )
import Data.List (sortBy, find)
import Data.Ord (comparing)
import Data.Maybe (fromMaybe)
#else
import Language.Fay.Prelude
-- Not provided in Fay set, it seems
comparing :: (b -> Double) -> b -> b -> Ordering
comparing p x y = compare (p x) (p y)
fromMaybe :: a -> Maybe a -> a
fromMaybe d x = case x of {Nothing -> d;Just v -> v}
#endif
type Radius = Double
type Circle a = (Double, a)
type Coordinate = (Double, Double)
type PlacedCircle a = ((Double, a), Coordinate)
type TouchingCircles a = [(PlacedCircle a, PlacedCircle a)]
{- | 'packCircles' takes a list of circles and a function that yields the
radius of the circle.
It returns a list of all circles, in unspecified order, together with
coordinates such that they do not overlap but sit as tight as possible,
filling a large circle.
Finding the optimal solution to this is NP hard, so only
heuristics are feasible. This particular
implementation is neither very good nor very fast,
compared to the state of the art in research. Nevertheless
it is simple to use and gives visually acceptable results.
The heuristics begins by placing the largest circle first, and the
next-to-largest next to it. From then on it adds circles by considering all
points where the circle to be added would touch two circles but overlap with
none, and picks the one that is closest to the center of mass of the current
placements.
-}
packCircles :: (a -> Double) -> [a] -> [(a, (Double, Double))]
packCircles radiusFunction =
-- Forget the cached radius
map (\t -> case t of (x,p) -> (snd x, p)) .
-- Make sure we center up the circles
centerUp .
-- Forget the pairs of touching circles
fst .
-- Run the main algorithm
go .
-- Look at large circles last
sortBy (comparing radius) .
-- Cache the radius
map (\x -> (radiusFunction x, x))
-- Just for a nicer name
radius :: Circle a -> Radius
radius = fst
-- Place the tail, then try to place the head
go :: [Circle a] -> ([PlacedCircle a], TouchingCircles a)
go [] = ([],[])
go (c:cs) = case go cs of
(placed, pairs) -> case place c placed pairs of
(cp, newPairs) -> (cp : placed, newPairs ++ pairs)
place :: Circle a -> [PlacedCircle a] -> TouchingCircles a ->
(PlacedCircle a, TouchingCircles a)
place c [] _ = ((c, (0,0)), [])
place c [cp'@(c',(x,y))] _ =
let cp = (c, (x + radius c + radius c', y))
in (cp, [(cp, cp')])
place c placed pairs = (cp, newPairs)
where
newPairs = [ (cp, cp') | cp' <- placed, touching cp cp' ]
cp = (c, p)
p = fromMaybe (error "packCircles: The end of the real plane has been reached?") $
find (\p' -> all (valid p') placed) $
sortBy (comparing centerDistance)
[ p' | (c1, c2) <- pairs,
touching c1 c2,
p' <- near c1 c2
]
centerDistance (x,y) = sqrt ((centerx - x)**2 + (centery - y)**2)
centerx = sum [x * (radius c')**2 | (c',(x,_)) <- placed] / area
centery = sum [y * (radius c')**2 | (c',(_,y)) <- placed] / area
area = sum [(radius c')**2 | (c',_) <- placed]
valid (x1,y1) (c2,(x2,y2)) =
sqrt ((x2 - x1)**2 + (y2-y1)**2) >= (radius c + radius c2) * (1-eps)
touching (c1,(x1,y1)) (c2, (x2, y2)) =
sqrt ((x2 - x1)**2 + (y2-y1)**2) <= (radius c1 + radius c2) * (1+eps)
near (c1,(x1,y1)) (c2, (x2, y2)) = [(c1x,c1y), (c2x,c2y)]
where
base = sqrt ((x2 - x1)**2 + (y2 - y1)**2)
lat1 = radius c1 + radius c
lat2 = radius c2 + radius c
--From http://stackoverflow.com/a/11356687/946226
ad_length = (base**2 + lat1**2 - lat2**2)/(2 * base)
h = sqrt (abs (lat1**2 - ad_length**2))
dx = x1 + ad_length * (x2 - x1)/base
dy = y1 + ad_length * (y2 - y1)/base
c1x = dx + h * (y2 - y1) / base
c1y = dy - h * (x2 - x1) / base
c2x = dx - h * (y2 - y1) / base
c2y = dy + h * (x2 - x1) / base
centerUp :: [PlacedCircle a] -> [PlacedCircle a]
centerUp placed = map (\(o,(x,y)) -> (o, (x-centerx, y-centery))) placed
where
centerx = sum [x * (radius c')**2 | (c',(x,_)) <- placed] / area
centery = sum [y * (radius c')**2 | (c',(_,y)) <- placed] / area
area = sum [(radius c')**2 | (c',_) <- placed]
eps :: Double
eps = 0.00001
{-$example
The following code demonstrates how one can use 'packCircles' together with
the diagrams library:
>import Diagrams.Prelude
>import Diagrams.Backend.SVG.CmdLine
>
>import Optimisation.CirclePacking
>
>colorize = zipWith fc $
> cycle [red,blue,yellow,magenta,cyan,bisque,firebrick,indigo]
>
>objects = colorize $
> [ circle r | r <- [0.1,0.2..1.6] ] ++
> [ hexagon r | r <- [0.1,0.2..0.7] ] ++
> [ decagon r | r <- [0.1,0.2..0.7] ]
>
>-- Just a approximation, diagram objects do not have an exact radius
>radiusApproximation o = maximum [ radius (e (CircleFrac alpha)) o | alpha <- [0,0.1..1.0]]
>
>main = defaultMain $
> position $ map (\(o,(x,y)) -> (p2 (x,y),o)) $
> packCircles radiusApproximation objects
This generates the following SVG file (if your browser manages to display it):
<>
-}