diff --git a/src/Main.hs b/src/Main.hs index f60d2a9..dc0fa8f 100644 --- a/src/Main.hs +++ b/src/Main.hs @@ -3,7 +3,9 @@ module Main where import Control.Arrow import qualified Control.Monad as M import qualified Control.Monad.Parallel as MP +import Control.Monad.Zip (MonadZip (mzipWith)) import Control.Parallel.Strategies +import qualified Data.Bifunctor as Bi import qualified Data.Colour as C import qualified Data.Colour.Names as CN import Data.Colour.RGBSpace (uncurryRGB) @@ -15,6 +17,7 @@ import qualified Data.Maybe as My import qualified Data.Set as S import qualified Data.Vector.Unboxed as Vec import Debug.Trace +import qualified Debug.Trace as D import qualified Debug.Trace as DT import qualified Debug.Trace as T import qualified Diagrams as DP @@ -26,6 +29,7 @@ import GHC.Generics import Graphics.Image as Img hiding (map, zipWith) import qualified Graphics.Image.ColorSpace as G import qualified Graphics.Image.Interface as Int +import qualified MinDistanceSample as MDS import Options.Generic import qualified System.Environment as Env import System.Random @@ -33,10 +37,6 @@ import System.Random.Internal import System.Random.SplitMix import Triangles (getTriangleAverageRGB) import qualified Triangles as Tri -import qualified MinDistanceSample as MDS -import Control.Monad.Zip (MonadZip(mzipWith)) -import qualified Data.Bifunctor as Bi -import qualified Debug.Trace as D -- CL.rgb might be the wrong fn... tosRGB' :: (Ord b, Floating b) => Pixel G.RGB b -> CL.Colour b @@ -50,78 +50,77 @@ corners = (,) <$> [0, 1] <*> [0, 1] shapeCircumference :: [Point V2 Double] -> Double shapeCircumference = Data.List.sum . map D.norm . loopOffsets . fromVertices -genImage :: Image VU G.RGB Double -> Double -> StdGen -> QDiagram SVG V2 Double Any -genImage image minDistance gen = - scaleX widthHeightRatio - . reflectY - . rectEnvelope (mkP2 0 0) (1 ^& 1) - . mconcat - . map drawVoroniRegion - . sortOn shapeCircumference - . withStrategy (parListChunk 50 rdeepseq) - . map (uncurry Tri.voroniDiagramCorners) - $ voroni - where - drawVoroniRegion shape = - lw 0 - . fillColor (Tri.voroniRegionAverageColor img' dimensions shape) - . strokeLocLoop - . fromVertices - $ shape +genImage :: Image VU G.RGB Double -> Double -> StdGen -> QDiagram SVG V2 Double Any +genImage image minDistance gen = + scaleX widthHeightRatio + . reflectY + . rectEnvelope (mkP2 0 0) (1 ^& 1) + . mconcat + . map drawVoroniRegion + . sortOn shapeCircumference + . withStrategy (parListChunk 50 rdeepseq) + . map (uncurry Tri.voroniDiagramCorners) + $ voroni + where + drawVoroniRegion shape = + lw 0 + . fillColor (Tri.voroniRegionAverageColor img' dimensions shape) + . strokeLocLoop + . fromVertices + $ shape - widthHeightRatio :: Double - widthHeightRatio = (fromIntegral . fst $ dimensions) / (fromIntegral . snd $ dimensions) + widthHeightRatio :: Double + widthHeightRatio = (fromIntegral . fst $ dimensions) / (fromIntegral . snd $ dimensions) - img' = convImage image - dimensions = uncurry (flip (,)) . Img.dims $ image - dimensionsVec = fromIntegral <$> uncurry V2 dimensions + img' = convImage image + dimensions = uncurry (flip (,)) . Img.dims $ image + dimensionsVec = fromIntegral <$> uncurry V2 dimensions - singleVoroni = last voroni + singleVoroni = last voroni - visualizeGraph :: QDiagram SVG V2 Double Any - visualizeGraph = - lc red - . lw 1 - . position - . map (\(x0, x1) -> (x0,) . strokeLine $ (x0 ~~ x1)) - . Tri.toPlanarGraph - $ corners' + visualizeGraph :: QDiagram SVG V2 Double Any + visualizeGraph = + lc red + . lw 1 + . position + . map (\(x0, x1) -> (x0,) . strokeLine $ (x0 ~~ x1)) + . Tri.toPlanarGraph + $ corners' - voroni = - Tri.findVoroniDiagram - . Tri.toPlanarGraph - $ corners' + voroni = + Tri.findVoroniDiagram + . Tri.toPlanarGraph + $ corners' + averageSideSize = (fromIntegral (uncurry (+) dimensions)) / 2 - averageSideSize = (fromIntegral (uncurry (+) dimensions)) / 2 + padding = (/) 10 . (*) widthHeightRatio <$> V2 1 1 - padding = (/) 10 . (*) widthHeightRatio <$> V2 1 1 - - corners' :: [P2 Double] - corners' = map (p2 . Bi.first (/ widthHeightRatio) . unp2 . (.-^ ((/ 2) <$> padding))) . MDS.randomPoints ((mkP2 widthHeightRatio 1) .+^ padding) ( minDistance * widthHeightRatio) $ gen + corners' :: [P2 Double] + corners' = map (p2 . Bi.first (/ widthHeightRatio) . unp2 . (.-^ ((/ 2) <$> padding))) . MDS.randomPoints ((mkP2 widthHeightRatio 1) .+^ padding) (minDistance * widthHeightRatio) $ gen deriving instance Generic (CL.RGB a) deriving instance NFData a => NFData (CL.RGB a) toDimensionVector :: (Int.BaseArray arr cs e, Fractional n) => Image arr cs e -> SizeSpec V2 n toDimensionVector image = - Diagrams.Prelude.dims $ - p2 (fromIntegral $ cols image, fromIntegral $ rows image) .-. p2 (0.0, 0.0) + Diagrams.Prelude.dims $ + p2 (fromIntegral $ cols image, fromIntegral $ rows image) .-. p2 (0.0, 0.0) data CLIOptions = CLIOptions - { input :: FilePath - , output :: FilePath - , minDistance :: Double - } - deriving (Generic) + { input :: FilePath + , output :: FilePath + , minDistance :: Double + } + deriving (Generic) instance ParseRecord CLIOptions main :: IO () main = do - CLIOptions{..} <- getRecord "image options" - gen' <- getStdGen -- for consistency, swap with something like: pure . mkStdGen $ 2344 - print gen' - image <- Img.readImageRGB VU input - let dimVector = toDimensionVector image - renderSVG output dimVector (genImage image minDistance gen') + CLIOptions{..} <- getRecord "image options" + gen' <- getStdGen -- for consistency, swap with something like: pure . mkStdGen $ 2344 + print gen' + image <- Img.readImageRGB VU input + let dimVector = toDimensionVector image + renderSVG output dimVector (genImage image minDistance gen') diff --git a/src/MinDistanceSample.hs b/src/MinDistanceSample.hs index dbc1039..c85d780 100644 --- a/src/MinDistanceSample.hs +++ b/src/MinDistanceSample.hs @@ -3,14 +3,14 @@ module MinDistanceSample where import qualified Control.Monad as M import qualified Data.Array as A import qualified Data.Bifunctor as B +import qualified Data.Ix as Ix import qualified Data.List as L import qualified Data.List.NonEmpty as NE +import qualified Data.Map.Strict as M import qualified Data.Maybe as My import qualified Debug.Trace as D import Diagrams.Prelude import System.Random.Stateful -import qualified Data.Map.Strict as M -import qualified Data.Ix as Ix k :: Int k = 10 @@ -68,6 +68,8 @@ randomPoints dims minDistance gen' = runStateGen_ gen' randomPointsM isValidPoint :: Point V2 Double -> Bool isValidPoint p = - (Ix.inRange gridBounds (floor <$> p)) && (all ((>= 1) . abs . norm . (p .-.)) - . My.mapMaybe ((grid M.!?) . fmap floor . (p .-^)) - $ unitVectorsAround) + (Ix.inRange gridBounds (floor <$> p)) + && ( all ((>= 1) . abs . norm . (p .-.)) + . My.mapMaybe ((grid M.!?) . fmap floor . (p .-^)) + $ unitVectorsAround + ) diff --git a/src/Triangles.hs b/src/Triangles.hs index 594013f..869833c 100644 --- a/src/Triangles.hs +++ b/src/Triangles.hs @@ -34,193 +34,190 @@ type Pixel_ = Colour Double toSRGBTuple :: Pixel_ -> (Double, Double, Double) toSRGBTuple = srgb' . C.toRGB - where - srgb' (C.RGB{C.channelRed = red, C.channelGreen = green, C.channelBlue = blue}) = (red, green, blue) + where + srgb' (C.RGB{C.channelRed = red, C.channelGreen = green, C.channelBlue = blue}) = (red, green, blue) fromSRGBTuple :: (Double, Double, Double) -> Pixel_ fromSRGBTuple (r, g, b) = C.rgb r g b derivingUnbox - "Pixel_" - [t|Pixel_ -> (Double, Double, Double)|] - [|toSRGBTuple|] - [|fromSRGBTuple|] - - + "Pixel_" + [t|Pixel_ -> (Double, Double, Double)|] + [|toSRGBTuple|] + [|fromSRGBTuple|] -- from -0.05 to 1.05 so there aren't missing/elongated triangles at the edges borderSize = 0.05 randomPoints :: StdGen -> [(Double, Double)] randomPoints = map (bimap toZeroToOneTuple toZeroToOneTuple) . randomRs ((0 :: Word, 0 :: Word), (maxBound, maxBound)) - where - toZeroToOneTuple :: Word -> Double - toZeroToOneTuple x = ((fromIntegral x / (fromIntegral (maxBound :: Word))) * (1 + (2 * borderSize))) - borderSize + where + toZeroToOneTuple :: Word -> Double + toZeroToOneTuple x = ((fromIntegral x / (fromIntegral (maxBound :: Word))) * (1 + (2 * borderSize))) - borderSize combinations :: (Ord b, Floating b, NFData b) => [P2 b] -> [(P2 b, P2 b)] combinations xs = - sortOn (abs . uncurry distanceA) - . S.toList -- deduplicate - . S.fromList - . filter (uncurry (/=)) - . concat - . withStrategy (parListChunk 50 rdeepseq) - . map (\(x : xs) -> take 10 . sortOn (abs . uncurry distanceA) . map (x,) $ xs) - . init -- last output of tails is empty list - . tails - $ xs - where - - xsLen = length xs + sortOn (abs . uncurry distanceA) + . S.toList -- deduplicate + . S.fromList + . filter (uncurry (/=)) + . concat + . withStrategy (parListChunk 50 rdeepseq) + . map (\(x : xs) -> take 10 . sortOn (abs . uncurry distanceA) . map (x,) $ xs) + . init -- last output of tails is empty list + . tails + $ xs + where + xsLen = length xs toPlanarGraph :: forall n. (NFData n, Floating n, Ord n) => [P2 n] -> [(Point V2 n, Point V2 n)] toPlanarGraph points = - removeIntersections - . sortOn (abs . uncurry distanceA) - . combinations - $ points - where - numPoints = length points + removeIntersections + . sortOn (abs . uncurry distanceA) + . combinations + $ points + where + numPoints = length points - removeIntersections :: [(Point V2 n, Point V2 n)] -> [(Point V2 n, Point V2 n)] - removeIntersections = foldl' addIfNoIntersection [] + removeIntersections :: [(Point V2 n, Point V2 n)] -> [(Point V2 n, Point V2 n)] + removeIntersections = foldl' addIfNoIntersection [] - addIfNoIntersection xs x - | all (noIntersection x) xs = x : xs - | otherwise = xs + addIfNoIntersection xs x + | all (noIntersection x) xs = x : xs + | otherwise = xs - noIntersection l1 l2 = sharedEndPoint || (null $ intersectPointsT (uncurry toLocatedTrail l1) (uncurry toLocatedTrail l2)) - where - sharedEndPoint = (< 4) . length . nub $ [fst l1, snd l1, fst l2, snd l2] + noIntersection l1 l2 = sharedEndPoint || (null $ intersectPointsT (uncurry toLocatedTrail l1) (uncurry toLocatedTrail l2)) + where + sharedEndPoint = (< 4) . length . nub $ [fst l1, snd l1, fst l2, snd l2] toLocatedTrail :: TrailLike a => Point (V a) (N a) -> Point (V a) (N a) -> Located a toLocatedTrail p1 p2 = fromVertices [p1, p2] `at` p1 withinShape :: RealFloat v => Point V2 v -> [Point V2 v] -> Point V2 v -> Bool withinShape pointInShape verticies candidate = all ((< quarterTurn) . fmap abs . uncurry signedAngleBetweenDirs) $ zip (shapeDirections pointInShape) (shapeDirections candidate) - where - shapeDirections p = map (dirBetween p) verticies + where + shapeDirections p = map (dirBetween p) verticies sortOnAngle center = sortOn (normalizeAngle . signedAngleBetweenDirs xDir . dirBetween center) voroniDiagramCorners :: RealFloat n => Point V2 n -> [Point V2 n] -> [Point V2 n] voroniDiagramCorners center midpoints = - sortOnAngle center - . filter isValidMidpoint - . concat - $ [intersectPointsT l0 l1 | l0 <- tangentTrails, l1 <- tangentTrails] - where - lessThanQuarterTurn candidate = candidate <= (10001 / 40000) @@ turn || candidate >= (29999 / 40000) @@ turn + sortOnAngle center + . filter isValidMidpoint + . concat + $ [intersectPointsT l0 l1 | l0 <- tangentTrails, l1 <- tangentTrails] + where + lessThanQuarterTurn candidate = candidate <= (10001 / 40000) @@ turn || candidate >= (29999 / 40000) @@ turn - tangentTrails = map tangentTrail midpoints + tangentTrails = map tangentTrail midpoints - appendHead (x : xs) = xs ++ [x] + appendHead (x : xs) = xs ++ [x] - isValidMidpoint candidate = all isNonObtuseMidpoint . filter (/= candidate) $ midpoints - where - isNonObtuseMidpoint m = - lessThanQuarterTurn . normalizeAngle $ angleBetweenDirs (dirBetween m center) (dirBetween m candidate) + isValidMidpoint candidate = all isNonObtuseMidpoint . filter (/= candidate) $ midpoints + where + isNonObtuseMidpoint m = + lessThanQuarterTurn . normalizeAngle $ angleBetweenDirs (dirBetween m center) (dirBetween m candidate) - tangentTrail midpoint = fromVertices [midpoint .-^ tangentVec, midpoint .+^ tangentVec] - where - -- implicitly uses the unit vector * 8 as an infinitely long vector - tangentVec = - scale 2 - . fromDirection - . rotateBy (1 / 4) - $ dirBetween midpoint center + tangentTrail midpoint = fromVertices [midpoint .-^ tangentVec, midpoint .+^ tangentVec] + where + -- implicitly uses the unit vector * 8 as an infinitely long vector + tangentVec = + scale 2 + . fromDirection + . rotateBy (1 / 4) + $ dirBetween midpoint center findVoroniDiagram :: RealFloat n => [(Point V2 n, Point V2 n)] -> [(Point V2 n, [Point V2 n])] findVoroniDiagram edges = - M.toList - . M.mapWithKey - ( \key -> - L.sortOn (normalizeAngle . signedAngleBetweenDirs xDir . dirBetween key) - . map (pointBetween key) - . S.toList - ) - $ adjacencyMap - where - adjacencyMap = adjacencyMapOf edges + M.toList + . M.mapWithKey + ( \key -> + L.sortOn (normalizeAngle . signedAngleBetweenDirs xDir . dirBetween key) + . map (pointBetween key) + . S.toList + ) + $ adjacencyMap + where + adjacencyMap = adjacencyMapOf edges - pointBetween p0 p1 = p0 .+^ ((p1 .-. p0) ^/ 2) + pointBetween p0 p1 = p0 .+^ ((p1 .-. p0) ^/ 2) findTriangles :: Ord b => [(b, b)] -> S.Set (S.Set b) findTriangles edges = S.unions . S.map threeCyclesOf . M.keysSet $ adjacencyMap - where - threeCyclesOf node = - S.unions - . S.map (\x -> S.map (\y -> S.fromList [node, x, y]) $ S.delete node . S.intersection originalNodeNeighbors . (M.!) adjacencyMap $ x) - $ originalNodeNeighbors - where - originalNodeNeighbors = fromMaybe S.empty (adjacencyMap M.!? node) + where + threeCyclesOf node = + S.unions + . S.map (\x -> S.map (\y -> S.fromList [node, x, y]) $ S.delete node . S.intersection originalNodeNeighbors . (M.!) adjacencyMap $ x) + $ originalNodeNeighbors + where + originalNodeNeighbors = fromMaybe S.empty (adjacencyMap M.!? node) - adjacencyMap = adjacencyMapOf edges + adjacencyMap = adjacencyMapOf edges adjacencyMapOf edges = M.fromListWith S.union . map (second S.singleton) $ (edges ++ edgesReversed) - where - edgesReversed = map (\(a, b) -> (b, a)) edges + where + edgesReversed = map (\(a, b) -> (b, a)) edges triangleAdjacencyMap :: Ord b => S.Set (S.Set b) -> M.Map b (S.Set (S.Set b)) triangleAdjacencyMap s = M.fromListWith S.union . concatMap (\s' -> map (,S.singleton s') . S.toList $ s') $ S.toList s getPointsInTriangle :: p -> S.Set (P2 Int) -> [(Int, Int)] getPointsInTriangle image pts = - S.toList . S.unions . map S.fromList $ - [ ptsBtween (makeLine p1 p3) (makeLine p1 p2) - , ptsBtween (makeLine p1 p3) (makeLine p2 p3) - , ptsBtween (makeLine p1 p2) (makeLine p2 p3) - ] - where - [p1, p2, p3] = sortOn fst . map unp2 . S.toList $ pts + S.toList . S.unions . map S.fromList $ + [ ptsBtween (makeLine p1 p3) (makeLine p1 p2) + , ptsBtween (makeLine p1 p3) (makeLine p2 p3) + , ptsBtween (makeLine p1 p2) (makeLine p2 p3) + ] + where + [p1, p2, p3] = sortOn fst . map unp2 . S.toList $ pts blendEqually colors = C.affineCombo (zip (repeat fraction) colors) $ C.blue - where - fraction = 1.0 / (fromIntegral . length $ colors) + where + fraction = 1.0 / (fromIntegral . length $ colors) voroniRegionAverageColor image (x', y') verticies = - blendEqually - . concatMap (getColorsInTriangle image (x', y')) - . filter ((== 3) . S.size) - . map (S.fromList . take 3) - . tails - . L.nub - . map scaleToImageCoords - $ verticies - where - scaleToImageCoords :: P2 Double -> P2 Int - scaleToImageCoords p = fmap round $ p2 ((fromIntegral x' * p ^. _x), fromIntegral y' * p ^. _y) + blendEqually + . concatMap (getColorsInTriangle image (x', y')) + . filter ((== 3) . S.size) + . map (S.fromList . take 3) + . tails + . L.nub + . map scaleToImageCoords + $ verticies + where + scaleToImageCoords :: P2 Double -> P2 Int + scaleToImageCoords p = fmap round $ p2 ((fromIntegral x' * p ^. _x), fromIntegral y' * p ^. _y) - scaleToUnitCoords :: P2 Int -> P2 Double - scaleToUnitCoords p = p2 ((fromIntegral x' / (fromIntegral $ p ^. _x)), fromIntegral y' / (fromIntegral $ p ^. _y)) + scaleToUnitCoords :: P2 Int -> P2 Double + scaleToUnitCoords p = p2 ((fromIntegral x' / (fromIntegral $ p ^. _x)), fromIntegral y' / (fromIntegral $ p ^. _y)) getColorsInTriangle :: Image_ -> (Int, Int) -> S.Set (P2 Int) -> [C.Colour Double] getColorsInTriangle image (x', y') triangle = pixels - where - pixels :: [Pixel_] - pixels = mapMaybe index' points + where + pixels :: [Pixel_] + pixels = mapMaybe index' points - points :: [(Int, Int)] - points = getPointsInTriangle image triangle + points :: [(Int, Int)] + points = getPointsInTriangle image triangle - index' :: (Int, Int) -> Maybe Pixel_ - index' (x, y) - | y >= y' = Nothing - | x >= x' = Nothing - | y < 0 = Nothing - | x < 0 = Nothing - | otherwise = image Vec.!? ((y * x') + x) + index' :: (Int, Int) -> Maybe Pixel_ + index' (x, y) + | y >= y' = Nothing + | x >= x' = Nothing + | y < 0 = Nothing + | x < 0 = Nothing + | otherwise = image Vec.!? ((y * x') + x) getTriangleAverageRGB :: Image_ -> (Int, Int) -> S.Set (P2 Int) -> C.Colour Double getTriangleAverageRGB image (x', y') triangle = blendEqually $ getColorsInTriangle image (x', y') triangle ptsBtween :: LineMXB -> LineMXB -> [(Int, Int)] ptsBtween l1 l2 = concatMap rasterLine . noSingletons $ [startingX .. endingX] - where - startingX = max (startX l1) (startX l2) - endingX = min (endX l1) (endX l2) + where + startingX = max (startX l1) (startX l2) + endingX = min (endX l1) (endX l2) - rasterLine x = map (x,) $ range' (yAt l1 x) (yAt l2 x) + rasterLine x = map (x,) $ range' (yAt l1 x) (yAt l2 x) noSingletons :: [a] -> [a] noSingletons [x] = [] @@ -234,22 +231,22 @@ yAt (LineMXB{m = m, b = b}) x = round $ (m * (fromIntegral x)) + b makeLine :: (Int, Int) -> (Int, Int) -> LineMXB makeLine (x1, y1) (x2, y2) = - LineMXB - { m = slope - , b = (fromIntegral y1) - (slope * (fromIntegral x1)) - , startX = min x1 x2 - , endX = max x1 x2 - } - where - slope = - if x1 /= x2 - then (fromIntegral $ y1 - y2) % (fromIntegral $ x1 - x2) - else fromIntegral . ceiling $ ((10.0 :: Double) ** 100.0) + LineMXB + { m = slope + , b = (fromIntegral y1) - (slope * (fromIntegral x1)) + , startX = min x1 x2 + , endX = max x1 x2 + } + where + slope = + if x1 /= x2 + then (fromIntegral $ y1 - y2) % (fromIntegral $ x1 - x2) + else fromIntegral . ceiling $ ((10.0 :: Double) ** 100.0) data LineMXB = LineMXB - { m :: Rational - , b :: Rational - , startX :: Int - , endX :: Int - } - deriving (Show, Ord, Eq) + { m :: Rational + , b :: Rational + , startX :: Int + , endX :: Int + } + deriving (Show, Ord, Eq)