dph-examples: pluecker line / triangle projection, also a very stupid way to draw...
[packages/dph.git] / dph-examples / examples / spectral / Pluecker / Common.hs
1 module Common where
2
3 -- Type synonyms so they work with Data.Vector.Unboxed etc.
4
5 -- A 3-vector *or* points
6 type Vec3 = (Double, Double, Double)
7 -- Line or ray between two points in space
8 -- These are both points, so direction of (a,b) is the vector (b-a)
9 type Line = (Vec3, Vec3)
10
11 -- Pl\"{u}cker coordinates of a line.
12 -- To convert line (a,b) into pluecker, (b-a, a X b)
13 type Pluecker = (Vec3,Vec3)
14
15 -- Triangles are cool because they're obviously always planar.
16 type Triangle = (Vec3, Vec3, Vec3)
17
18 -- Normal and distance of normal from (0,0,0)
19 type Plane = (Vec3,Double)
20
21
22 {-# INLINE dot #-}
23 dot :: Vec3 -> Vec3 -> Double
24 dot (u,v,w) (x,y,z) = u*x + v*y + w*z
25
26 {-# INLINE cross #-}
27 cross :: Vec3 -> Vec3 -> Vec3
28 cross (u,v,w) (x,y,z) = ( v * z - w * y
29 , w * x - u * z
30 , u * y - v * x)
31
32 {-# INLINE over1 #-}
33 over1 :: (Double -> Double) -> Vec3 -> Vec3
34 over1 f (x,y,z) = (f x, f y, f z)
35
36 {-# INLINE over2 #-}
37 over2 :: (Double -> Double -> Double) -> Vec3 -> Vec3 -> Vec3
38 over2 f (u,v,w) (x,y,z) = (f u x, f v y, f w z)
39
40 {-# INLINE mag #-}
41 mag a = sqrt (a `dot` a)
42
43 {-# INLINE norm #-}
44 norm a = over1 (/m) a
45 where m = mag a
46
47 {-# INLINE vsub #-}
48 vadd = over2 (+)
49 vsub = over2 (-)
50
51 -- Convert a line into pluecker coordinates
52 {-# INLINE plueckerOfLine #-}
53 plueckerOfLine :: Line -> Pluecker
54 plueckerOfLine (p,q) = (q `vsub` p, q `cross` p)
55
56 -- Find intersection of a line and plucker, if exists.
57 -- Otherwise
58 {-# INLINE projectPluecker2 #-}
59 projectPluecker2 :: Pluecker -> Pluecker -> Double
60 projectPluecker2 (u1,v1) (u2,v2) = (u1 `dot` v2) + (u2 `dot` v1)
61
62 -- Check whether pluecker line intersects given triangle
63 {-# INLINE inside #-}
64 inside :: Pluecker -> Triangle -> Bool
65 inside p (a,b,c)
66 = let pr = projectPluecker2 p . plueckerOfLine
67 ab = pr (a,b) < 0
68 bc = pr (b,c) < 0
69 ca = pr (c,a) < 0
70 in (ab && bc && ca) || not (ab || bc || ca)
71
72
73 -- Get plane from triangle, eg for projection
74 {-# INLINE planeOfTriangle #-}
75 planeOfTriangle :: Triangle -> Plane
76 planeOfTriangle (a,b,c)
77 = let n = (a `vsub` b) `cross` (c `vsub` b)
78 n' = norm n
79 d = negate (n' `dot` b)
80 in (n', d)
81
82 -- should this return Maybe Double?
83 {-# INLINE lineOnPlane #-}
84 lineOnPlane :: Line -> Plane -> Double
85 lineOnPlane (p,q) (n,d)
86 = let d1 = (n `dot` p) + d
87 d2 = (n `dot` q) + d
88 in if d1 == d2
89 then 1e100 -- meh. big. fail.
90 else d1 / (d1 - d2)
91
92 -- Project line onto triangle's plane
93 -- Disregard whether intersection is actually inside triangle
94 -- Intersection point is scale from line
95 {-# INLINE lineOnTriangle #-}
96 lineOnTriangle :: Line -> Triangle -> Double
97 lineOnTriangle p t
98 = lineOnPlane p (planeOfTriangle t)