Add source files, commit before reverting pictures to lists
This commit is contained in:
+576
@@ -0,0 +1,576 @@
|
||||
{-# LANGUAGE BangPatterns #-}
|
||||
module Geometry
|
||||
( module Geometry
|
||||
, module Geometry.Data
|
||||
)
|
||||
where
|
||||
import Geometry.Data
|
||||
|
||||
import Data.Function
|
||||
import Data.List
|
||||
import Data.Maybe
|
||||
|
||||
import Control.Applicative
|
||||
|
||||
zeroZ :: Point2 -> Point3
|
||||
zeroZ (x,y) = (x,y,0)
|
||||
|
||||
infixl 6 +.+, -.-
|
||||
infixl 7 *.*
|
||||
|
||||
(+.+) :: Point2 -> Point2 -> Point2
|
||||
{-# INLINE (+.+) #-}
|
||||
(x1, y1) +.+ (x2, y2) =
|
||||
let
|
||||
!x = x1 + x2
|
||||
!y = y1 + y2
|
||||
in (x, y)
|
||||
|
||||
(-.-) :: Point2 -> Point2 -> Point2
|
||||
{-# INLINE (-.-) #-}
|
||||
(x1, y1) -.- (x2, y2) =
|
||||
let
|
||||
!x = x1 - x2
|
||||
!y = y1 - y2
|
||||
in (x, y)
|
||||
|
||||
(*.*) :: Float -> Point2 -> Point2
|
||||
{-# INLINE (*.*) #-}
|
||||
a *.* (x2, y2) =
|
||||
let
|
||||
!x = a * x2
|
||||
!y = a * y2
|
||||
in (x, y)
|
||||
|
||||
infixl 6 +.+.+, -.-.-
|
||||
infixl 7 *.*.*
|
||||
|
||||
(+.+.+) :: Point3 -> Point3 -> Point3
|
||||
{-# INLINE (+.+.+) #-}
|
||||
(x1, y1, z1) +.+.+ (x2, y2, z2) =
|
||||
let
|
||||
!x = x1 + x2
|
||||
!y = y1 + y2
|
||||
!z = z1 + z2
|
||||
in (x, y, z)
|
||||
|
||||
(-.-.-) :: Point3 -> Point3 -> Point3
|
||||
{-# INLINE (-.-.-) #-}
|
||||
(x1, y1, z1) -.-.- (x2, y2, z2) =
|
||||
let
|
||||
!x = x1 - x2
|
||||
!y = y1 - y2
|
||||
!z = z1 - z2
|
||||
in (x, y, z)
|
||||
|
||||
(*.*.*) :: Point3 -> Point3 -> Point3
|
||||
{-# INLINE (*.*.*) #-}
|
||||
(x1, y1, z1) *.*.* (x2, y2, z2) =
|
||||
let
|
||||
!x = x1 * x2
|
||||
!y = y1 * y2
|
||||
!z = z1 * z2
|
||||
in (x, y, z)
|
||||
|
||||
normalizeV :: Point2 -> Point2
|
||||
{-# INLINE normalizeV #-}
|
||||
normalizeV p = (1 / magV p) *.* p
|
||||
|
||||
angleVV :: Point2 -> Point2 -> Float
|
||||
{-# INLINE angleVV #-}
|
||||
angleVV a b = let ma = magV a
|
||||
mb = magV b
|
||||
d = a `dotV` b
|
||||
in acos $ d / (ma * mb)
|
||||
|
||||
dotV :: Point2 -> Point2 -> Float
|
||||
{-# INLINE dotV #-}
|
||||
dotV (x,y) (z,w) = x*z + y*w
|
||||
|
||||
closestPointOnLine :: Point2 -> Point2 -> Point2 -> Point2
|
||||
{-# INLINE closestPointOnLine #-}
|
||||
closestPointOnLine a b p
|
||||
= a +.+ u *.* (b -.- a)
|
||||
where u = closestPointOnLineParam a b p
|
||||
|
||||
closestPointOnLineParam :: Point2 -> Point2 -> Point2 -> Float
|
||||
{-# INLINE closestPointOnLineParam #-}
|
||||
closestPointOnLineParam a b p
|
||||
= (p -.- a) `dotV` (b -.- a) / (b -.- a) `dotV` (b -.- a)
|
||||
|
||||
argV :: Point2 -> Float
|
||||
{-# INLINE argV #-}
|
||||
argV (x,y) = normalizeAngle $ atan2 y x
|
||||
|
||||
detV :: Point2 -> Point2 -> Float
|
||||
{-# INLINE detV #-}
|
||||
detV (x1, y1) (x2, y2)
|
||||
= x1 * y2 - y1 * x2
|
||||
|
||||
-- | Angle in radians, anticlockwise from +ve x-axis.
|
||||
unitVectorAtAngle :: Float -> Point2
|
||||
{-# INLINE unitVectorAtAngle #-}
|
||||
unitVectorAtAngle r
|
||||
= (cos r, sin r)
|
||||
|
||||
-- | Rotate a vector by an angle (in radians). +ve angle is counter-clockwise.
|
||||
rotateV :: Float -> Point2 -> Point2
|
||||
rotateV r (x, y)
|
||||
= ( x * cos r - y * sin r
|
||||
, x * sin r + y * cos r)
|
||||
{-# INLINE rotateV #-}
|
||||
|
||||
-- | Convert degrees to radians
|
||||
degToRad :: Float -> Float
|
||||
degToRad d = d * pi / 180
|
||||
{-# INLINE degToRad #-}
|
||||
|
||||
|
||||
-- | Convert radians to degrees
|
||||
radToDeg :: Float -> Float
|
||||
radToDeg r = r * 180 / pi
|
||||
{-# INLINE radToDeg #-}
|
||||
|
||||
|
||||
-- | Normalize an angle to be between 0 and 2*pi radians
|
||||
normalizeAngle :: Float -> Float
|
||||
normalizeAngle f = f - 2 * pi * floor' (f / (2 * pi))
|
||||
where floor' :: Float -> Float
|
||||
floor' x = fromIntegral (floor x :: Int)
|
||||
{-# INLINE normalizeAngle #-}
|
||||
|
||||
-- the following helper draws a rectangle based on maximal N E S W values
|
||||
rectNESW :: Float -> Float -> Float -> Float -> [Point2]
|
||||
rectNESW a b c d = [(b,a),(b,c),(d,c),(d,a)
|
||||
]
|
||||
rectNSEW n s e w = rectNESW n e s w
|
||||
|
||||
rectNSWE :: Float -> Float -> Float -> Float -> [Point2]
|
||||
rectNSWE n s w e = [ (w,n), (w,s), (e,s), (e,n)]
|
||||
|
||||
-- -- the following filters points in a polygon: supposes the points in the
|
||||
-- polygon are listed in anticlockwise order
|
||||
pointInOrOnPolygon :: Point2 -> [Point2] -> Bool
|
||||
pointInOrOnPolygon p (x:xs) = all (\l -> not (uncurry isRHS l p)) $ zip (x:xs) (xs ++ [x])
|
||||
|
||||
|
||||
pointInPolygon :: Point2 -> [Point2] -> Bool
|
||||
pointInPolygon p [] = False
|
||||
pointInPolygon p (x:xs) = all (\l -> uncurry (errorIsLHS 1) l p) $ zip (x:xs) (xs ++ [x])
|
||||
|
||||
errorPointInPolygon :: Int -> Point2 -> [Point2] -> Bool
|
||||
errorPointInPolygon i p xs | length xs == 1 = error "one point polygon"
|
||||
| length xs == 2 = error "two point polygon"
|
||||
| nub xs == xs = pointInPolygon p xs
|
||||
| otherwise = error $ "errorPointInPolygon "++ show i
|
||||
|
||||
errorNormalizeV :: Int -> Point2 -> Point2
|
||||
errorNormalizeV i (0,0) = error $ "problem with function: errorNormalizeV "++show i
|
||||
errorNormalizeV i p = normalizeV p
|
||||
|
||||
errorAngleVV :: Int -> Point2 -> Point2 -> Float
|
||||
errorAngleVV i (0,0) _ = error $ "problem with function: errorAngleVV "++show i
|
||||
errorAngleVV i _ (0,0) = error $ "problem with function: errorAngleVV "++show i
|
||||
errorAngleVV i p p' = angleVV p p'
|
||||
|
||||
errorIsLHS :: Int -> Point2 -> Point2 -> Point2 -> Bool
|
||||
errorIsLHS i x y | x == y = error $ "problem with function: errorIsLHS "
|
||||
++show i
|
||||
| otherwise = isLHS x y
|
||||
|
||||
errorClosestPointOnLine :: Int -> Point2 -> Point2 -> Point2 -> Point2
|
||||
errorClosestPointOnLine i x y | x == y = error $ "problem with function: errorClosestPointOnLine "
|
||||
++show i
|
||||
| otherwise = closestPointOnLine x y
|
||||
|
||||
errorClosestPointOnLineParam :: Int -> Point2 -> Point2 -> Point2 -> Float
|
||||
errorClosestPointOnLineParam i x y z | x == y = dist x z
|
||||
-- error $ "problem with function: errorClosestPointOnLineParam " ++show i
|
||||
| otherwise = closestPointOnLineParam x y z
|
||||
|
||||
safeNormalizeV :: Point2 -> Point2
|
||||
safeNormalizeV (0,0) = (0,0)
|
||||
safeNormalizeV p = normalizeV p
|
||||
|
||||
-- tests whether a point is on the LHS of a line
|
||||
-- this has been called somewhere with l1 == l2
|
||||
isLHS :: Point2 -> Point2 -> Point2 -> Bool
|
||||
{-# INLINE isLHS #-}
|
||||
isLHS' l1 l2 p | l1 == l2 = False
|
||||
| otherwise = closestPointOnLineParam l1 (l1 +.+ vNormal (l2 -.- l1)) p < 0
|
||||
|
||||
isLHS (x,y) (x',y') (x'',y'')
|
||||
| (x,y) == (x',y') = False
|
||||
| otherwise = a1 * b2 - a2 * b1 > 0
|
||||
where a1 = x' - x
|
||||
a2 = y' - y
|
||||
b1 = x'' - x
|
||||
b2 = y'' - y
|
||||
|
||||
isRHS :: Point2 -> Point2 -> Point2 -> Bool
|
||||
{-# INLINE isRHS #-}
|
||||
isRHS (x,y) (x',y') (x'',y'')
|
||||
| (x,y) == (x',y') = False
|
||||
| otherwise = a1 * b2 - a2 * b1 < 0
|
||||
where a1 = x' - x
|
||||
a2 = y' - y
|
||||
b1 = x'' - x
|
||||
b2 = y'' - y
|
||||
--isRHS l1 l2 p = closestPointOnLineParam l1 (l1 +.+ vNormal (l2 -.- l1)) p > 0
|
||||
|
||||
-- reorders points to be anticlockwise around their center
|
||||
orderPolygon :: [Point2] -> [Point2]
|
||||
orderPolygon [] = []
|
||||
orderPolygon ps = sortBy (compare `on` \p -> argV (p -.- cen)) ps
|
||||
where cen = 1/ fromIntegral (length ps) *.* foldr1 (+.+) ps
|
||||
|
||||
vNormal :: Point2 -> Point2
|
||||
{-# INLINE vNormal #-}
|
||||
vNormal (x,y) = (y,-x)
|
||||
|
||||
vInverse :: Point2 -> Point2
|
||||
vInverse (x,y) = (-x,-y)
|
||||
|
||||
dist :: Point2 -> Point2 -> Float
|
||||
{-# INLINE dist #-}
|
||||
dist p1 p2 = magV (p2 -.- p1)
|
||||
|
||||
normV :: Point2 -> Point2
|
||||
{-# INLINE normV #-}
|
||||
normV (0,0) = (0,0)
|
||||
normV p = (1/magV p ) *.* p
|
||||
|
||||
magV :: Point2 -> Float
|
||||
{-# INLINE magV #-}
|
||||
magV (x,y) = sqrt $ x^2 + y^2
|
||||
|
||||
pHalf :: Point2 -> Point2 -> Point2
|
||||
pHalf a b = 0.5 *.* (a +.+ b)
|
||||
|
||||
linGrad :: (Eq a,Fractional a) => (a,a) -> (a,a) -> Maybe a
|
||||
linGrad (x,y) (a,b) | x-a == 0 = Nothing
|
||||
| otherwise = Just $ (y-b)/(x-a)
|
||||
|
||||
axisInt :: (Eq a,Fractional a) => (a,a) -> (a,a) -> Maybe a
|
||||
axisInt p (a,b) = pure b ^-^ (pure a ^*^ linGrad p (a,b))
|
||||
where (^-^) = liftA2 (-)
|
||||
(^*^) = liftA2 (*)
|
||||
|
||||
-- intersectSegSeg is sometimes broken-- the following fixes at least some of
|
||||
-- the cases
|
||||
-- it is, however, slow
|
||||
myIntersectSegSeg a@(ax,ay) b@(bx,by) c@(cx,cy) d@(dx,dy) = case ratIntersectLineLine a b c d of
|
||||
Nothing -> Nothing
|
||||
Just (x,y) -> if inbetween x && inbetween' y
|
||||
then Just (x,y)
|
||||
else Nothing
|
||||
where inbetween x = ((ax <= x && x <= bx) || (bx <= x && x <= ax)) &&
|
||||
((cx <= x && x <= dx) || (dx <= x && x <= cx))
|
||||
inbetween' y = ((ay <= y && y <= by) || (by <= y && y <= ay)) &&
|
||||
((cy <= y && y <= dy) || (dy <= y && y <= cy))
|
||||
|
||||
crossV :: Point2 -> Point2 -> Float
|
||||
crossV (ax,ay) (bx,by) = ax*by - ay*bx
|
||||
|
||||
myIntersectLineLine :: (Eq a,Fractional a) => (a,a) -> (a,a) -> (a,a) -> (a,a) -> Maybe (a,a)
|
||||
myIntersectLineLine a@(ax,ay) b c@(cx,cy) d
|
||||
| linGrad a b == Nothing = fmap ((,) ax) $ axisInt (c *-* (ax,0)) (d *-* (ax,0))
|
||||
| linGrad c d == Nothing = fmap ((,) cx) $ axisInt (a *-* (cx,0)) (b *-* (cx,0))
|
||||
| otherwise
|
||||
= case linGrad a b ^-^ linGrad c d of
|
||||
Just 0 -> Nothing
|
||||
_ -> liftA2 (,) newx
|
||||
((linGrad a b ^*^ newx) ^+^ axisInt a b)
|
||||
where (^-^) = liftA2 (-)
|
||||
(^+^) = liftA2 (+)
|
||||
(^/^) = liftA2 (/)
|
||||
(^*^) = liftA2 (*)
|
||||
newx = (axisInt c d ^-^ axisInt a b) ^/^ (linGrad a b ^-^ linGrad c d)
|
||||
(*-*) (ax,ay) (bx,by) = (ax-bx,ay-by)
|
||||
|
||||
ratIntersectLineLine :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
|
||||
ratIntersectLineLine a b c d = fmap toNumPoint2 $ myIntersectLineLine (toRatPoint2 a) (toRatPoint2 b) (toRatPoint2 c) (toRatPoint2 d)
|
||||
where toRatPoint2 (x,y) = (toRational x, toRational y)
|
||||
toNumPoint2 (x,y) = (fromRational x, fromRational y)
|
||||
f = toRatPoint2 . roundPoint2
|
||||
|
||||
roundPoint2 :: Point2 -> Point2
|
||||
roundPoint2 (x,y) = (fromIntegral $ round x,fromIntegral $ round y)
|
||||
|
||||
circOnLine' :: Point2 -> Point2 -> Point2 -> Float -> Bool
|
||||
circOnLine' p1 p2 c rad = isJustTrue (fmap (\p -> magV (p -.- c) < rad) y)
|
||||
where y = intersectSegLine' p1 p2 c (c +.+ vNormal (p1 -.- p2))
|
||||
isJustTrue (Just True) = True
|
||||
isJustTrue _ = False
|
||||
|
||||
circOnLine :: Point2 -> Point2 -> Point2 -> Float -> Bool
|
||||
circOnLine p1 p2 c rad = magV (p1 -.- c) <= rad || magV (p2 -.- c) <= rad
|
||||
|| isJustTrue (fmap (\p -> magV (p -.- c) < rad) y)
|
||||
where y = intersectSegLine' p1 p2 c (c +.+ vNormal (p1 -.- p2))
|
||||
isJustTrue (Just True) = True
|
||||
isJustTrue _ = False
|
||||
|
||||
difference x y | x > y = x - y
|
||||
| otherwise = y - x
|
||||
|
||||
reflectIn :: Point2 -> Point2 -> Point2
|
||||
reflectIn line vec = let angle = 2 * angleBetween line vec
|
||||
in rotateV angle vec
|
||||
|
||||
angleBetween v1 v2 = argV v1 - argV v2
|
||||
|
||||
doublePair :: (a,a) -> [(a,a)]
|
||||
doublePair (x,y) = [(x,y),(y,x)]
|
||||
|
||||
polysIntersect :: [Point2] -> [Point2] -> Bool
|
||||
polysIntersect (p:ps) (q:qs)
|
||||
= any isJust $ (\(a,b) (c,d) -> myIntersectSegSeg a b c d) <$> pairs1 <*> pairs2
|
||||
where pairs1 = zip (p:ps) (ps++[p])
|
||||
pairs2 = zip (q:qs) (qs++[q])
|
||||
|
||||
anyPolyssIntersect :: [[Point2]] -> [[Point2]] -> Bool
|
||||
anyPolyssIntersect x y = or $ polysIntersect <$> x <*> y
|
||||
|
||||
nRays :: Int -> [Point2]
|
||||
nRays n = take n $ iterate (rotateV (2*pi/fromIntegral n)) $ (600,0)
|
||||
|
||||
nRaysRad :: Int -> Float -> [Point2]
|
||||
nRaysRad n x = take n $ iterate (rotateV (2*pi/fromIntegral n)) $ (x,0)
|
||||
|
||||
-- angles go from 0 to 2pi, need to work out what is left of another
|
||||
|
||||
isLeftOfA :: Float -> Float -> Bool
|
||||
isLeftOfA angle1 angle2 = (angle1 - angle2 < pi && angle1 > angle2)
|
||||
|| (angle2 - angle1 > pi && angle2 > angle1)
|
||||
|
||||
isLeftOf :: Point2 -> Point2 -> Bool
|
||||
isLeftOf x y = isLeftOfA (argV x) (argV y)
|
||||
|
||||
-- diffAngles has an issue...
|
||||
|
||||
diffAngles :: Float -> Float -> Float
|
||||
diffAngles x y | diff > pi = diffAngles (x - 2*pi) y
|
||||
| diff >= 0 = diff
|
||||
| diff > -pi = -diff
|
||||
| otherwise = diffAngles (x + 2*pi) y
|
||||
where diff = x-y
|
||||
differenceAngles = diffAngles
|
||||
angleDifference = diffAngles
|
||||
|
||||
-- given a triangle where we know the length of a first side,
|
||||
-- the length of a second side, and the angle between the first side and the
|
||||
-- third side, finds the length of the third side
|
||||
-- not this doesn't necessarily find ALL solutions, asin is a map not a function
|
||||
ssaTri :: Float -> Float -> Float -> Float
|
||||
ssaTri ab bc a
|
||||
| sin a == 0 = 0
|
||||
| bc == 0 = ab
|
||||
| otherwise = let c = asin ( (ab * (sin a))/bc)
|
||||
b = pi - (a + c)
|
||||
in sin b * bc / sin a
|
||||
|
||||
-- fix points: we now fix the triangle in the coordinate system, and return a
|
||||
-- third unknown point:
|
||||
-- the point which lies between pa and pc' on a line from b of length bc
|
||||
-- note that there are likely two such points, this seems to return the point
|
||||
-- closer to pc'
|
||||
|
||||
ssaTriPoint :: Point2 -> Point2 -> Point2 -> Float -> Point2
|
||||
ssaTriPoint pa pb pc' bc
|
||||
= let ab = magV (pa -.- pb)
|
||||
a = errorAngleVV 6 (pb -.- pa) (pc' -.- pa)
|
||||
ac = ssaTri ab bc a
|
||||
in pa +.+ (ac *.* errorNormalizeV 47 (pc' -.- pa))
|
||||
|
||||
-- the above SHOULD return a Maybe Point...
|
||||
ssaTriPoint' :: Point2 -> Point2 -> Point2 -> Float -> Maybe Point2
|
||||
ssaTriPoint' pa pb pc' bc
|
||||
| dist pb (closestPointOnSeg pa pc' pb) >= bc
|
||||
= Nothing
|
||||
| otherwise
|
||||
= Just $ ssaTriPoint pa pb pc' bc
|
||||
ssaTriPointCorrect :: Point2 -> Point2 -> Point2 -> Float -> Maybe Point2
|
||||
ssaTriPointCorrect pa pb pc' bc
|
||||
| param <= 1 && param >= 0 = Just p
|
||||
| otherwise = Nothing
|
||||
where p = ssaTriPoint pa pb pc' bc
|
||||
param = closestPointOnLineParam pa pc' p
|
||||
|
||||
|
||||
closestPointOnSeg :: Point2 -> Point2 -> Point2 -> Point2
|
||||
closestPointOnSeg segP1 segP2 p
|
||||
| errorClosestPointOnLineParam 3 segP1 segP2 p <= 0 = segP1
|
||||
| errorClosestPointOnLineParam 4 segP1 segP2 p >= 1 = segP2
|
||||
| otherwise = errorClosestPointOnLine 2 segP1 segP2 p
|
||||
|
||||
pointInCircle :: Point2 -> Float -> Point2 -> Maybe Point2
|
||||
pointInCircle p r c | p == c = Just p
|
||||
| magV (p -.- c) < r = Just p
|
||||
| otherwise = Nothing
|
||||
|
||||
--determines if a moving point intersects with a circle,
|
||||
--if so, returns a point on circle that intersects with the line passing
|
||||
--throught the circle : HOPEFULLY THE CORRECT OF THE TWO!
|
||||
|
||||
collidePointCirc :: Point2 -> Point2 -> Float -> Point2 -> Maybe Point2
|
||||
collidePointCirc p1 p2 rad c = ssaTriPoint' p2 c p1 rad
|
||||
|
||||
-- changes the point to a measure of the distance
|
||||
collidePointCirc' :: Point2 -> Point2 -> Float -> Point2 -> Maybe Float
|
||||
collidePointCirc' p1 p2 rad c = fmap (\x -> magV (x -.- p1))
|
||||
(collidePointCirc p1 p2 rad c)
|
||||
|
||||
--returns both the point and the measure of the distance
|
||||
collidePointCirc'' :: Point2 -> Point2 -> Float -> Point2 -> Maybe (Point2,Float)
|
||||
collidePointCirc'' p1 p2 rad c = (,) <$> collidePointCirc p1 p2 rad c
|
||||
<*> collidePointCirc' p1 p2 rad c
|
||||
|
||||
collidePointCircCorrect :: Point2 -> Point2 -> Float -> Point2 -> Maybe Point2
|
||||
collidePointCircCorrect p1 p2 rad c = ssaTriPointCorrect p2 c p1 rad
|
||||
|
||||
|
||||
-- finds the height of a triangle using herons formula
|
||||
-- the base is the line between the first two points
|
||||
heron :: Point2 -> Point2 -> Point2 -> Float
|
||||
heron x y z | x == y = 0
|
||||
| otherwise = let a = magV $ x -.- y
|
||||
b = magV $ y -.- z
|
||||
c = magV $ z -.- x
|
||||
s = (a+b+c)/2
|
||||
area = sqrt(s*(s-a)*(s-b)*(s-c))
|
||||
in 2*area/a
|
||||
-- multiplies reflection in normal by factor
|
||||
reflectInParam :: Float -> Point2 -> Point2 -> Point2
|
||||
reflectInParam x line vec = let angle = 2 * angleBetween line vec
|
||||
rAng = rotateV angle vec
|
||||
p = x *.* errorClosestPointOnLine 3 (0,0) (vNormal line) rAng
|
||||
in rAng -.- p
|
||||
|
||||
|
||||
reflectIn' :: Point2 -> Point2 -> Point2 -> Point2 -> Point2
|
||||
reflectIn' l1 l2 v1 v2 = v1 +.+ reflectIn (l1 -.- l2) (v2 -.- v1)
|
||||
|
||||
isOnLine :: Point2 -> Point2 -> Point2 -> Bool
|
||||
isOnLine l1 l2 p = errorClosestPointOnLineParam 10 l1 (l1 +.+ vNormal (l2 -.- l1)) p == 0
|
||||
&& errorClosestPointOnLineParam 11 l1 l2 p <= 1
|
||||
&& errorClosestPointOnLineParam 12 l1 l2 p >= 0
|
||||
|
||||
-- the take 5000 here is a hack, otherwise divideLine seems to sometimes
|
||||
-- generate an infinite list, and I don't know why
|
||||
divideLine :: Float -> Point2 -> Point2 -> [Point2]
|
||||
--divideLine x a b = map (\i -> a +.+ (i / (fromIntegral numPoints)) *.* (b -.- a))
|
||||
divideLine x a b = take 5000 $ map (\i -> a +.+ (i / (fromIntegral numPoints) *.* (b -.- a)) )
|
||||
$ map fromIntegral ns
|
||||
where d = dist a b
|
||||
numPoints = max 1 $ ceiling $ d / x
|
||||
ns = [0 .. numPoints]
|
||||
|
||||
|
||||
-- pulled the following from the haskell wiki
|
||||
bresenham :: (Int,Int) -> (Int,Int) -> [(Int,Int)]
|
||||
{-# INLINE bresenham #-}
|
||||
bresenham pa@(xa,ya) pb@(xb,yb) = map maySwitch . unfoldr go $ (x1,y1,0)
|
||||
where
|
||||
steep = abs (yb - ya) > abs (xb - xa)
|
||||
maySwitch = if steep then (\(x,y) -> (y,x)) else id
|
||||
[(x1,y1),(x2,y2)] = sort [maySwitch pa, maySwitch pb]
|
||||
deltax = x2 - x1
|
||||
deltay = abs (y2 - y1)
|
||||
ystep = if y1 < y2 then 1 else -1
|
||||
go (xTemp, yTemp, error)
|
||||
| xTemp > x2 = Nothing
|
||||
| otherwise = Just ((xTemp, yTemp), (xTemp + 1, newY, newError))
|
||||
where
|
||||
tempError = error + deltay
|
||||
(newY, newError) = if (2*tempError) >= deltax
|
||||
then (yTemp+ystep,tempError-deltax)
|
||||
else (yTemp,tempError)
|
||||
|
||||
|
||||
divideCircle :: Float -> Point2 -> Float -> [Point2]
|
||||
divideCircle x cen rad = map (cen +.+) $ nPointsOnCirc n rad
|
||||
where n = ceiling $ rad * 2 * pi / x
|
||||
|
||||
nPointsOnCirc :: Int -> Float -> [Point2]
|
||||
nPointsOnCirc n rad = take n $ iterate (rotateV (2*pi/fromIntegral n)) $ (rad,0)
|
||||
|
||||
lineInPolygon :: Point2 -> Point2 -> [Point2] -> Bool
|
||||
lineInPolygon a b ps = pointInPolygon a ps || pointInPolygon b ps
|
||||
|| any (isJust . uncurry (intersectSegSeg' a b)) pss
|
||||
where pss = zip ps (tail ps ++ [head ps])
|
||||
|
||||
--intersectSegLineFrom :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
|
||||
--intersectSegLineFrom a b c d
|
||||
-- = case intersectSegLine a b c d of
|
||||
-- Just p | closestPointOnLineParam c d p >= 0 -> Just p
|
||||
-- | otherwise -> Nothing
|
||||
-- Nothing -> Nothing
|
||||
|
||||
intersectSegSeg' :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
|
||||
{-# INLINE intersectSegSeg' #-}
|
||||
intersectSegSeg' (x1,y1) (x2,y2) (x3,y3) (x4,y4)
|
||||
-- | t' < 0 || u' < 0 || t' > den || u' > den || den == 0
|
||||
-- = Nothing
|
||||
| den == 0 = Nothing
|
||||
| den > 0 && (t' < 0 || u' < 0 || t' > den || u' > den)
|
||||
= Nothing
|
||||
| den < 0 && (t' > 0 || u' > 0 || t' < den || u' < den)
|
||||
= Nothing
|
||||
-- | den > 0 && (t' < 0 || t' > den)
|
||||
-- = Nothing
|
||||
-- | den < 0 && (t' > 0 || t' < 0-den)
|
||||
-- = Nothing
|
||||
| otherwise = Just (x1 + (x2-x1)*t'/den, y1 + (y2-y1)*t'/den)
|
||||
where
|
||||
den = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
|
||||
t' = (x1-x3)*(y3-y4) - (y1-y3)*(x3-x4)
|
||||
u' = (y1-y2)*(x1-x3) - (x1-x2)*(y1-y3)
|
||||
|
||||
intersectSegLineFrom' :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
|
||||
{-# INLINE intersectSegLineFrom' #-}
|
||||
intersectSegLineFrom' (x1,y1) (x2,y2) (x3,y3) (x4,y4)
|
||||
-- | t' < 0 || u' < 0 || t' > den || u' > den || den == 0
|
||||
-- = Nothing
|
||||
| den == 0 = Nothing
|
||||
| den > 0 && ( t' < 0 || u' < 0 || t' > den )
|
||||
= Nothing
|
||||
| den < 0 && ( t' > 0 || u' > 0 || t' < den )
|
||||
= Nothing
|
||||
-- | den > 0 && (t' < 0 || t' > den)
|
||||
-- = Nothing
|
||||
-- | den < 0 && (t' > 0 || t' < 0-den)
|
||||
-- = Nothing
|
||||
| otherwise = Just (x1 + (x2-x1)*t'/den, y1 + (y2-y1)*t'/den)
|
||||
where
|
||||
den = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
|
||||
t' = (x1-x3)*(y3-y4) - (y1-y3)*(x3-x4)
|
||||
u' = (y1-y2)*(x1-x3) - (x1-x2)*(y1-y3)
|
||||
|
||||
intersectSegLine' :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
|
||||
{-# INLINE intersectSegLine' #-}
|
||||
intersectSegLine' (x1,y1) (x2,y2) (x3,y3) (x4,y4)
|
||||
-- | t' < 0 || u' < 0 || t' > den || u' > den || den == 0
|
||||
-- = Nothing
|
||||
| den == 0 = Nothing
|
||||
| den > 0 && (t' < 0 || t' > den)
|
||||
= Nothing
|
||||
| den < 0 && (t' > 0 || t' < den)
|
||||
= Nothing
|
||||
-- | den > 0 && (t' < 0 || t' > den)
|
||||
-- = Nothing
|
||||
-- | den < 0 && (t' > 0 || t' < 0-den)
|
||||
-- = Nothing
|
||||
| otherwise = Just (x1 + (x2-x1)*t'/den, y1 + (y2-y1)*t'/den)
|
||||
where
|
||||
den = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
|
||||
t' = (x1-x3)*(y3-y4) - (y1-y3)*(x3-x4)
|
||||
u' = (y1-y2)*(x1-x3) - (x1-x2)*(y1-y3)
|
||||
|
||||
intersectLineLine' :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
|
||||
{-# INLINE intersectLineLine' #-}
|
||||
intersectLineLine' (x1,y1) (x2,y2) (x3,y3) (x4,y4)
|
||||
| den == 0 = Nothing
|
||||
| otherwise = Just (x1 + (x2-x1)*t'/den, y1 + (y2-y1)*t'/den)
|
||||
where
|
||||
den = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
|
||||
t' = (x1-x3)*(y3-y4) - (y1-y3)*(x3-x4)
|
||||
Reference in New Issue
Block a user