{- | Copyright : Copyright (C) 2012 Joachim Breitner License : BSD3 Maintainer : Joachim Breitner Stability : stable Portability: portable -} module Optimisation.CirclePacking ( packCircles -- * Example -- $example ) where import Data.List (sortBy, minimumBy) import Data.Ord (comparing) {- | '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 radius = go . sortBy (comparing radius) where go [] = [] go (c:cs) = let placed = go cs in (c, place c placed) : placed place _ [] = (0,0) place c [(c',(x,y))] = (x + radius c + radius c', y) place c placed = minimumBy (comparing centerDistance) [ p | (c1, c2) <- allPairs placed, touching c1 c2, p <- near c1 c2, all (valid p) placed ] where 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 -- Possible ways to speed this up: -- * Cache radius -- * Remember which circles are touching -- * First sort all points, then take first valid, instead of checking validity -- for all of them. eps :: Double eps = 0.00001 allPairs :: [a] -> [(a,a)] allPairs [] = [] allPairs (x:xs) = [ (x,y) | y <- xs ] ++ allPairs xs {-$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 manges to display it): <> -}