image-triangles/src/Triangles.hs

{-# LANGUAGE TupleSections, TypeFamilies, MultiParamTypeClasses, TemplateHaskell, TypeSynonymInstances, FlexibleInstances #-}
{-# LANGUAGE InstanceSigs #-}
-- module Triangles
-- ( getRandomPixel
-- , getRandomTriangle
-- , getPointsInTriangle
-- , getTriangleAverageRGB
-- ) where

module Triangles where

import System.Random
import Data.Ratio
import qualified Debug.Trace
import qualified Data.Colour.SRGB.Linear as C
import Control.Parallel.Strategies
import qualified Data.Colour as C
import Data.Colour.SRGB.Linear (Colour)
import Data.Vector.Generic.Base (Vector)
import Data.Vector.Generic.Mutable (MVector)
import Data.List
import Data.Maybe
import qualified Data.Vector.Unboxed as Vec
import Data.Fixed
import Data.Vector.Unboxed.Deriving
import Diagrams.TwoD
import Diagrams.Prelude
import qualified Data.Map as M
import qualified Data.Set as S
import Control.Arrow
import qualified Data.List as L
import qualified Data.Function as F
import qualified Diagrams.TwoD.Path.IntersectionExtras as I
import Diagrams.Trail (trailPoints)
import qualified Graphics.Image as H
import Debug.Trace (traceShow)
import qualified Data.Colour.Names as C
import qualified Debug.Trace as D
import Diagrams.TwoD.Segment.Bernstein (listToBernstein)
-- import qualified Linear.Affine as L
type Image_ = Vec.Vector Pixel_
type Pixel_ = Colour Double
-- type Point = (Double, Double)
-- type Triangle = (Point, Point, Point)

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)

fromSRGBTuple :: (Double, Double, Double) -> Pixel_
fromSRGBTuple (r, g, b) = C.rgb r g b

derivingUnbox "Pixel_"
    [t| Pixel_ -> (Double, Double, Double)|]
    [| toSRGBTuple |]
    [| fromSRGBTuple |]

-- from -0.05 to 1.05 so there aren't missing/elongated triangles at the edges
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) - 0.1

combinations :: (Ord b, Floating b, NFData b) => [P2 b] -> [(P2 b, P2 b)]
combinations xs =
  sortOn (abs . uncurry distanceA)
  . S.toList . S.fromList
  . filter (uncurry (/=))
  . sortOn (abs . uncurry distanceA)
  . concatMap (\(x:xs) -> map (x,) xs)
  . init . tails $ xs
  where
    edgeLengthThreshold = 10

toPlanarGraph :: forall n . (NFData n, Floating n, Ord n) => [P2 n] -> [(Point V2 n, Point V2 n)]
toPlanarGraph points =
  removeIntersections .
  combinations $ points
  where
    numPoints = length points

    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

    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

-- data VoroniRegion = VoroniRegion {
--   center :: P2 Double,
--   neighbors :: [P2 Double]
--                                  }

-- voroniDiagramCorners :: Point V2 Double -> [Point V2 Double] -> [Point V2 Double]
-- voroniDiagramCorners :: Point (V c) (N c) -> [Point V2 Double] -> c
-- voroniDiagramCorners :: forall t. (N t ~ Double, V t ~ V2, TrailLike t) => Point V2 Double -> [Point V2 Double] -> t
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


sortOnAngle center = sortOn (normalizeAngle . signedAngleBetweenDirs xDir . dirBetween center)

-- voroniDiagramCorners :: forall n . (Floating n, Real n, Show n) => Point V2 n -> [P2 n] -> [P2 n]
voroniDiagramCorners center midpoints -- = midpoints
  = sortOnAngle center . filter isValidMidpoint . concat $ [intersectPointsT l0 l1 | l0 <- tangentTrails, l1 <- tangentTrails]
  where

    lessThanQuarterTurn candidate = candidate <= (10001 / 40000) @@ turn || candidate >= (29999 / 40000) @@ turn
      -- where
      --   candidate' = (candidate ^. turn)  - (fromIntegral . floor $ (candidate ^. turn))
    -- candidateVertecies = sortOn (normalizeAngle . signedAngleBetweenDirs xDir . dirBetween center) . concat $ [ intersectPointsT (tangentTrail x) (tangentTrail y)  | x <- midpoints, y <- midpoints ]

    tangentTrails = map tangentTrail midpoints
    -- validMidpoints = sortOn (normalizeAngle . signedAngleBetweenDirs xDir . dirBetween center) . filter isValidMidpoint $ midpoints

    appendHead (x:xs) = xs ++ [x]

    isValidMidpoint candidate = all isNonObtuseMidpoint . filter (/= candidate) $ midpoints
      where
        isNonObtuseMidpoint m = -- traceShow ((dirBetween m center, dirBetween m candidate), (candidate, m)) .
          lessThanQuarterTurn . normalizeAngle $ angleBetweenDirs (dirBetween m center) (dirBetween m candidate)

    tangentTrail midpoint = fromVertices $ [midpoint .-^ tangentVec, midpoint .+^ tangentVec]
      where
        -- implicitly uses the unit vector as an infinitely long vector
        tangentVec = scale 8 .
          fromDirection . rotateBy (1 / 4) $ dirBetween midpoint center


-- findVoroniDiagram :: (Ord (v n), Additive v, Fractional n) => [(Point v n, Point v n)] -> [[Point v n]]
-- findVoroniDiagram :: (Ord n, Ord (v n), Metric v, Floating n, R1 v, Real n, r) => [(Point v n, Point v n)] -> [[Point v n]]
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

    -- I'm not sure this part works
    pointBetween p0 p1 = p0 .+^ ((p1 .-. p0) ^/ 2)

-- filterNotInVoroniRegion :: RealFloat n => Point V2 n -> [Point V2 n] -> [Point V2 n] -> [Point V2 n]
-- filterNotInVoroniRegion center midpoints = id -- filter allowed
--   where
--     allowed point
--       = all ((< 0.25) . abs . (^. turn))
--       . zipWith signedAngleBetweenDirs midpointAngles
--       . map (dirBetween point) $ midpoints

--     midpointAngles = map (dirBetween center) midpoints


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)

    adjacencyMap = adjacencyMapOf edges

adjacencyMapOf edges = M.fromListWith S.union . map (second S.singleton) $ (edges ++ edgesReversed)
  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 (\(y, x) -> (x,y)) . map unp2 . S.toList $ pts

blendEqually colors = C.affineCombo (zip (repeat fraction) colors) $ C.blue
    where
        fraction = 1.0 / (fromIntegral . length $ colors)

-- voroniRegionAverageColor :: Image_ -> (Int, Int) -> (P2 Int, [P2 Int]) -> C.Colour Double
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)

    scaleToUnitCoords :: P2 Int -> P2 Double
    scaleToUnitCoords p = p2 ((fromIntegral x' / (fromIntegral $ p ^. _x)), fromIntegral y' / (fromIntegral $ p ^. _y))
  -- = blendEqually
  -- . mapMaybe (index' . unp2)
  -- . filter withinShape'
  -- $ candidatePoints
  -- where

  --   withinShape' :: P2 Int -> Bool
  --   withinShape' = withinShape center verticies . scaleToUnitCoords

  --   candidatePoints = [p2 (x, y) | x <- [minX .. maxX], y <- [minY .. maxY]]

  --   maxX = fst . unp2 $ maximumBy (compare `F.on` fst . unp2) verticies'
  --   minX = fst . unp2 $ minimumBy (compare `F.on` fst . unp2) verticies'
  --   maxY = snd . unp2 $ maximumBy (compare `F.on` snd . unp2) verticies'
  --   minY = snd . unp2 $ minimumBy (compare `F.on` snd . unp2) verticies'

  --   verticies' = map scaleToImageCoords verticies

  --   index' :: (Int, Int) -> Maybe Pixel_
  --   index' (y, x)
  --       | y >= y' = Nothing
  --       | x >= x' = Nothing
  --       | y < 0   = Nothing
  --       | x < 0   = Nothing
  --       | otherwise = image Vec.!? ((y * x') + x)


getColorsInTriangle :: Image_ -> (Int, Int) -> S.Set (P2 Int)  -> [C.Colour Double]
getColorsInTriangle image (x', y') triangle = pixels
    where

        pixels :: [Pixel_]
        pixels = mapMaybe index' points

        points :: [(Int, Int)]
        points = getPointsInTriangle image triangle

        index' :: (Int, Int) -> Maybe Pixel_
        index' (y, x)
            | 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)

        rasterLine x = map (\y -> (y, x)) $ range' (yAt l1 x) (yAt l2 x)

noSingletons :: [a] -> [a]
noSingletons [x] = []
noSingletons l   = l

range' :: Int -> Int -> [Int]
range' a b = [(min a b) .. (max a b)]

yAt :: LineMXB -> Int -> Int
yAt (LineMXB {m = m, b = b}) x = round $ (m * (fromIntegral x)) + b

makeLine :: (Int, Int) -> (Int, Int) -> LineMXB
makeLine (y1, x1) (y2, x2) = 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)