{-# 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): <> -}