dph-examples: pluecker. work around bug, possibly vectoriser? "not" behaving strangely.
[packages/dph.git] / dph-examples / examples / spectral / Pluecker / CommonVectorised.hs
1 {-# OPTIONS -fvectorise #-}
2 module CommonVectorised where
3 import Data.Array.Parallel.Prelude.Bool
4 import Data.Array.Parallel.Prelude.Double as D hiding (pi)
5 import qualified Prelude as P
6
7 -- Type synonyms so they work with Data.Vector.Unboxed etc.
8
9 -- A 3-vector *or* points
10 type Vec3 = (Double, Double, Double)
11 -- Line or ray between two points in space
12 -- These are both points, so direction of (a,b) is the vector (b-a)
13 type Line = (Vec3, Vec3)
14
15 -- Pl\"{u}cker coordinates of a line.
16 -- To convert line (a,b) into pluecker, (b-a, a X b)
17 type Pluecker = (Vec3,Vec3)
18
19 -- Triangles are cool because they're obviously always planar.
20 type Triangle = (Vec3, Vec3, Vec3)
21
22 -- Normal and distance of normal from (0,0,0)
23 type Plane = (Vec3,Double)
24
25 {-# INLINE dot #-}
26 dot :: Vec3 -> Vec3 -> Double
27 dot (u,v,w) (x,y,z) = u*x + v*y + w*z
28
29 {-# INLINE cross #-}
30 cross :: Vec3 -> Vec3 -> Vec3
31 cross (u,v,w) (x,y,z) = ( v * z - w * y
32 , w * x - u * z
33 , u * y - v * x)
34
35 {-# INLINE over1 #-}
36 over1 :: (Double -> Double) -> Vec3 -> Vec3
37 over1 f (x,y,z) = (f x, f y, f z)
38
39 {-# INLINE over2 #-}
40 over2 :: (Double -> Double -> Double) -> Vec3 -> Vec3 -> Vec3
41 over2 f (u,v,w) (x,y,z) = (f u x, f v y, f w z)
42
43 {-# INLINE mag #-}
44 mag a = sqrt (a `dot` a)
45
46 {-# INLINE norm #-}
47 norm a = over1 (/m) a
48 where m = mag a
49
50 {-# INLINE vsub #-}
51 vadd = over2 (+)
52 vsub = over2 (-)
53
54 -- Convert a line into pluecker coordinates
55 {-# INLINE plueckerOfLine #-}
56 plueckerOfLine :: Line -> Pluecker
57 plueckerOfLine (p,q) = (q `vsub` p, q `cross` p)
58
59 -- Find intersection of a line and plucker, if exists.
60 -- Otherwise
61 {-# INLINE projectPluecker2 #-}
62 projectPluecker2 :: Pluecker -> Pluecker -> Double
63 projectPluecker2 (u1,v1) (u2,v2) = (u1 `dot` v2) + (u2 `dot` v1)
64
65 -- Check whether pluecker line intersects given triangle
66 {-# INLINE inside #-}
67 inside :: Pluecker -> Triangle -> Bool
68 inside p (a,b,c)
69 = let pr l = projectPluecker2 p (plueckerOfLine l)
70 ab = pr (a,b) < 0
71 bc = pr (b,c) < 0
72 ca = pr (c,a) < 0
73 in (ab && bc && ca) || not' (ab || bc || ca)
74 where
75 -- TODO ISSUE WTF
76 -- Vectoriser bug?
77 -- Before, I was just using "not (ab || ...)". It compiled fine but when running it I got this error.
78 -- MainGloss: Data.Array.Parallel.Unlifted.Sequential.Flat.combine2ByTag: tags length /= sum of args length (first = 50000; second = 148577)
79 not' True = False
80 not' False = True
81
82
83 -- Get plane from triangle, eg for projection
84 {-# INLINE planeOfTriangle #-}
85 planeOfTriangle :: Triangle -> Plane
86 planeOfTriangle (a,b,c)
87 = let n = (a `vsub` b) `cross` (c `vsub` b)
88 n' = norm n
89 d = negate (n' `dot` b)
90 in (n', d)
91
92 -- should this return Maybe Double?
93 {-# INLINE lineOnPlane #-}
94 lineOnPlane :: Line -> Plane -> Double
95 lineOnPlane (p,q) (n,d)
96 = let d1 = (n `dot` p) + d
97 d2 = (n `dot` q) + d
98 in if d1 == d2
99 then 1e100 -- meh. big. fail.
100 else d1 / (d1 - d2)
101
102 -- Project line onto triangle's plane
103 -- Disregard whether intersection is actually inside triangle
104 -- Intersection point is scale from line
105 {-# INLINE lineOnTriangle #-}
106 lineOnTriangle :: Line -> Triangle -> Double
107 lineOnTriangle p t
108 = lineOnPlane p (planeOfTriangle t)
109