{-# LANGUAGE DataKinds, TypeFamilies, FlexibleContexts, FlexibleInstances, PolyKinds #-}
{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses #-}
module Polynomial.Hypersurface where
import Polynomial.Prelude
import Polynomial.Monomial
import Geometry.ConvexHull2
import qualified Data.Map.Strict as MS
import qualified Data.Sized as DS
import Data.Maybe
import Debug.Trace
import Data.List
import Geometry.ConvexHull3 (Point3D)
import Geometry.Polytope
type Polygon = [Point2D]
type Normals = [Point2D]
mapTermPoint :: (IsMonomialOrder ord, Ord k, Integral k)
=> Polynomial k ord n -> MS.Map Point2D (Monomial ord n, k)
mapTermPoint poly = MS.fromList $ zip points terms
where
terms = MS.toList $ getTerms poly
monExps = DS.toList . getMonomial
toPoints (mon,coef) = let [a,b] = monExps mon in (a, b)
points = map toPoints terms
computeIntersection :: (Integral k) => (Monomial ord n, k) -> (Monomial ord n, k) -> (Monomial ord n, k) -> Point2D
computeIntersection (mon1, c1) (mon2, c2) (mon3, c3) = (x, y)
where
[a,b] = DS.toList $ getMonomial mon1
[d,e] = DS.toList $ getMonomial mon2
[g,h] = DS.toList $ getMonomial mon3
c = fromIntegral c1
f = fromIntegral c2
i = fromIntegral c3
y = ((f-c)*(d-g) - (i-f)*(a-d)) `div` ((b-e)*(d-g)-(e-h)*(a-d))
x = ((f-c)-(b-e)*y) `div` (a-d)
findFanNVertex :: (Integral k) => MS.Map Point2D (Monomial ord n, k) -> Polygon -> (Point2D,Normals)
findFanNVertex mapPointMon points = (computeIntersection mon1 mon2 mon3, sort [inormal1,inormal2,inormal3])
where
[pp1,pp2,pp3] = sort points
mon1 = fromJust $ MS.lookup pp3 mapPointMon
mon2 = fromJust $ MS.lookup pp1 mapPointMon
mon3 = fromJust $ MS.lookup pp2 mapPointMon
inormal1 = innerNormal pp1 pp2 pp3
inormal2 = innerNormal pp2 pp3 pp1
inormal3 = innerNormal pp3 pp1 pp2
findPolygonNVertex :: (Integral k) => MS.Map Point2D (Monomial ord n, k) -> Polygon -> (Polygon, Point2D)
findPolygonNVertex mapPointMon points = ( [pp1,pp2,pp3], computeIntersection mon1 mon2 mon3)
where
[pp1,pp2,pp3] = sort points
mon1 = fromJust $ MS.lookup pp3 mapPointMon
mon2 = fromJust $ MS.lookup pp1 mapPointMon
mon3 = fromJust $ MS.lookup pp2 mapPointMon
innerNormal :: Point2D -> Point2D -> Point2D -> Point2D
innerNormal a@(x1,y1) b@(x2,y2) c@(x3,y3)
| dot > 0 = simplify nab
| dot < 0 = simplify (y1-y2,x2-x1)
where
ab = (x2-x1,y2-y1)
ac = (x3-x1,y3-y1)
nab = (y2-y1,x1-x2)
dot = (y2-y1)*(x3-x1) + (x1-x2)*(y3-y1)
simplify (0, q) = (0, div q (abs q))
simplify (q, 0) = (div q (abs q), 0)
simplify (p, q) = let g = gcd p q in (div p g, div q g)
innerNormals :: Point2D -> Point2D -> Point2D -> Normals
innerNormals a@(x1,y1) b@(x2,y2) c@(x3,y3) = map simplify [inner1, inner2, inner3]
where
inner1 = innerNormal a b c
inner2 = innerNormal b c a
inner3 = innerNormal c a b
simplify (0, q) = (0, div q (abs q))
simplify (q, 0) = (div q (abs q), 0)
simplify (p, q) = let g = gcd p q in (div p g, div q g)
verticesNormals :: (IsMonomialOrder ord, Ord k, Integral k) => Polynomial k ord n -> MS.Map Point2D Normals
verticesNormals poly = MS.fromList $ map (findFanNVertex polyMap) triangles
where
polyMap = mapTermPoint poly
triangles = subdivision poly
neighborTriangles :: [Polygon] -> MS.Map Polygon [Polygon] -> MS.Map Polygon [Polygon]
neighborTriangles [] mapContainer = mapContainer
neighborTriangles (p:ps) mapContainer
| null ps && MS.null mapContainer = MS.fromList [(p, [])]
| otherwise = MS.map (map sort) $ neighborTriangles ps (foldr (lookAndInsert p) mapContainer ps)
where
lookAndInsert p1 p2 acc = if length (p1\\p2) == (length p1) - 2 then
MS.insertWith (++) p2 [p1] $ MS.insertWith (++) p1 [p2] acc
else acc
pointsWithTriangles :: (IsMonomialOrder ord, Ord k, Integral k) => Polynomial k ord n -> MS.Map Polygon Point2D
pointsWithTriangles poly = MS.fromList $ map (findPolygonNVertex polyMap) triangles
where
polyMap = mapTermPoint poly
triangles = subdivision poly
polygonCenter :: MS.Map Polygon Point2D -> Polygon -> Point2D
polygonCenter pointPolygonMap polygon = fromJust $ MS.lookup polygon pointPolygonMap
convertMap :: MS.Map Polygon [Polygon] -> MS.Map Polygon Point2D -> MS.Map Point2D [Point2D]
convertMap map1 map2 = MS.fromList $ map fromPolygons $ MS.toList map1
where
fromPolygons (polygon, polygons) = (polygonCenter map2 polygon, map (polygonCenter map2) polygons)
computeEdges :: MS.Map Point2D [Point2D] -> MS.Map Point2D Normals -> [(Point2D, Point2D)]
computeEdges map1 map2 = concatMap getEdges pointsWithNormals
where
listMap1 = MS.toList map1
listMap2 = MS.toList map2
getNormals point2D = (point2D, fromJust $ MS.lookup point2D map2)
attachNormals = map (\(e,l) -> (getNormals e, map getNormals l))
getEdges ((p,n), []) = map (\normal-> (p, p+ (10 >*< normal))) n
getEdges ((p,n), (p1,n1):ps) = let ((point, newNormals),edge) = analizeNormals (p,n) (p1,n1) in
edge:(getEdges ((point, newNormals), ps))
pointsWithNormals = attachNormals listMap1
isInverse :: Point2D -> Point2D -> Point2D -> Point2D ->Bool
isInverse (x1,y1) (nx1, ny1) (x2,y2) (nx2, ny2)
| x1 == x2 = nx1 == 0 && nx2 == 0 && ny1*ny2 < 0
| y1 == y2 = ny1 == 0 && ny2 == 0 && nx1*nx2 < 0
| nx1 == 0 = False
| div (y2-y1) (x2-x1) == div ny1 nx1 = True
| otherwise = False
analizeNormals :: (Point2D, Normals) -> (Point2D, Normals) -> ((Point2D,Normals),(Point2D, Point2D))
analizeNormals (p1, n1) (p2, n2) = ((p1, newNormals), (p1,p2))
where
isThereTwin p1 normal p2 normals = any (isInverse p1 normal p2) normals
newNormals = foldr (\normal acc -> if isThereTwin p1 normal p2 n2 then acc else normal:acc) [] n1
hypersurface :: (IsMonomialOrder ord, Ord k, Integral k) => Polynomial k ord n -> [(Point2D, Point2D)]
hypersurface poly = nub $ computeEdges (convertMap neighbors pointTriangles) pointNormals
where
pointNormals = verticesNormals poly
neighbors = neighborTriangles (map sort $ subdivision poly) MS.empty
pointTriangles = pointsWithTriangles poly