module Triangles where

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

type Image_ = Vec.Vector Pixel_
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)

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

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

borderSize = 0.05

-- 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 * 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 15 . sortOn (abs . uncurry distanceA) . 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
    . 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 []

  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

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 :: 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

  tangentTrails = map tangentTrail midpoints

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

  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 8
        . 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

  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)

  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 unp2 . S.toList $ pts

blendEqually colors = C.affineCombo (zip (repeat fraction) colors) $ C.blue
 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)

  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

  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)

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 (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 (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)

data LineMXB = LineMXB
  { m :: Rational
  , b :: Rational
  , startX :: Int
  , endX :: Int
  }
  deriving (Show, Ord, Eq)
