dualist: reflex3d2d tweaks
[polytopiary:reflex.git] / reflex3d2d.hs
1 module Main where
2
3 import Control.Monad(forever, when)
4 import Data.IORef(IORef, newIORef, atomicModifyIORef)
5 import Control.Concurrent(forkIO, yield)
6 import System.IO(hGetBuf, withBinaryFile, IOMode(ReadMode), hSetBuffering, BufferMode(LineBuffering), stdin)
7 import Foreign.Marshal.Alloc(allocaBytes)
8 import Foreign.Marshal.Array(peekArray)
9 import Foreign.Storable(sizeOf)
10 import Graphics.UI.GLUT hiding(Point, project)
11
12 import Vector hiding (i, o)
13
14 fr3 :: Int -> Int -> [V3]
15 fr3 p q =
16   let v0 = V3 1 0 0
17       v1 = m1 ^^*^ v0
18       v2 = (m1 ^^*^^ m2) ^^*^ v0
19       m1 = M3 c1 s1 0   (-s1) c1 0   0 0 1
20       m2 = M3 c2 0 s2   0 1 0   (-s2) 0 c2
21       c1 = r1 / r0
22       c2 = r2 / r1
23       s1 = sqrt (1 - c1*c1)
24       s2 = sqrt (1 - c2*c2)
25       r0 = radius 0 [p,q]
26       r1 = radius 1 [p,q]
27       r2 = radius 2 [p,q]
28   in [ norm $ cross3 v1 v2
29      , norm $ cross3 v0 v2 ^* (-1)
30      , norm $ cross3 v0 v1
31      , v0, v1, v2 ]
32
33 delta :: [Int] -> R
34 delta [] = 1
35 delta [p] = let s = sin (pi / fromIntegral p) in s * s
36 delta (p:q:rs) = let c = cos (pi / fromIntegral p) in delta (q:rs) - delta rs * c * c
37
38 radius2 :: Int -> [Int] -> R
39 radius2 0 [] = 1
40 radius2 0 (p:qs) = delta qs / delta (p:qs)
41 radius2 j ps = radius2 0 ps - radius2 0 (take (j-1) ps) 
42
43 radius :: Int -> [Int] -> R
44 radius n ps = sqrt $ radius2 n ps
45
46 dropV :: V v => v -> v -> v
47 dropV u v =
48   let k = u`dot`v / u`dot`u
49   in v ^-^ (k *^ u)
50
51 data Figure3 = Figure3 !V3 !V3 !V3 !V3 !V3 !V3 !V3 !Bool
52
53 instance Approx Figure3 where
54   Figure3 p a b c _ _ _ _ ~=~ Figure3 q e f g _ _ _ _ = and $ zipWith (~=~) [p,a,b,c] [q,e,f,g]
55
56 initialP3 :: [V3] -> V3 -> V3
57 initialP3 [V3 a11 a12 a13, V3 a21 a22 a23, V3 a31 a32 a33, _, _, _] v =
58   let m = M3 a11 a12 a13 a21 a22 a23 a31 a32 a33
59       u = norm $ inv m ^^*^ v
60   in u
61 initialP3 _ _ = error "initialP3"
62
63 initialF3 :: [V3] -> V3 -> Figure3
64 initialF3 [l0,l1,l2,v0,_,v2] p = Figure3 p a b c d e f True
65   where (a, b, c) = (dropV l0 p , dropV l1 p , dropV l2 p)
66         (d, e, f) = (v2 , v0 , (dropV l0 . dropV l2) p)
67 initialF3 _ _ = error "initialF3"
68
69 transformF3 :: M3 -> Figure3 -> Figure3
70 transformF3 m (Figure3 p a b c d e f z) = Figure3 (m ^^*^ p) (m ^^*^ a) (m ^^*^ b) (m ^^*^ c) (m ^^*^ d) (m ^^*^ e) (m ^^*^ f) (if det m > 0 then z else not z)
71
72 transformsF3 :: M3 -> [Figure3] -> [Figure3]
73 transformsF3 m fs = map (m `transformF3`) fs
74
75 polyhedronF :: [M3] -> Figure3 -> [Figure3]
76 polyhedronF ms f = map (`transformF3` f) ms
77
78
79 quality = 6
80
81 drawFigure3 :: (Color3 R, Color3 R, Color3 R) -> (Bool, Bool, Bool) -> Figure3 -> IO ()
82 drawFigure3 (ac,bc,cc) (av,bv,cv) (Figure3 p a b c d e f z) = do
83   let seg p1 p2 = let V2 x1 y1 = project p1
84                       V2 x2 y2 = project p2
85                   in vertex (Vertex2 x1 y1) >> vertex (Vertex2 x2 y2)
86       lin p1 p2 n = let dir = (p2 ^-^ p1 ^/ fromIntegral n) in map (\i -> norm $ p1 ^+^ (fromIntegral i *^ dir)) [0 .. n :: Int]
87       line p1 p2 = let ps = lin p1 p2 quality in sequence_ $ zipWith seg ps (tail ps)
88   when av $ color ac >> line p a >> line p d
89   when bv $ color bc >> line p b >> line p e
90   when cv $ color cc >> line p c >> line p f
91
92 drawFigures3 :: (Color3 R, Color3 R, Color3 R) -> (Bool, Bool, Bool) -> [Figure3] -> IO ()
93 drawFigures3 cs v fs = renderPrimitive Lines $ mapM_ (drawFigure3 cs v) fs
94
95
96 project :: V3 -> V2
97 project p = let (V3 x y z) = norm p in let k = 0.5 / (1 - z) in V2 (k*x) (k*y)
98
99 fillFigure3 :: (Color3 R, Color3 R, Color3 R) -> (Bool, Bool, Bool) -> Figure3 -> IO ()
100 fillFigure3 (ac,bc,cc) (av,bv,cv) (Figure3 p a b c d e f parity) = do
101   when av $ color ac >> triangle' p a d >> triangle' p d b
102   when bv $ color bc >> triangle' p b e >> triangle' p e c
103   when cv $ color cc >> triangle' p c f >> triangle' p f a
104   where
105     triangle' p1 p2 p3 = triangle parity quality p1 p2 p3
106
107 clockwise :: V2 -> V2 -> V2 -> Bool
108 clockwise (V2 ax ay) (V2 bx by) (V2 cx cy) = e1x * e2y - e1y * e2x >= 0
109   where
110     (e1x, e1y) = (ax - bx, ay - by)
111     (e2x, e2y) = (cx - bx, cy - by)
112
113 triangle :: Bool -> Int -> V3 -> V3 -> V3 -> IO ()
114 triangle cw 0 p1 p2 p3 = do
115   let v1@(V2 x1 y1) = project p1
116       v2@(V2 x2 y2) = project p2
117       v3@(V2 x3 y3) = project p3
118   when (clockwise v1 v2 v3 == cw) $ do
119     vertex (Vertex2 x1 y1)
120     vertex (Vertex2 x2 y2)
121     vertex (Vertex2 x3 y3)
122
123 triangle cw n p1 p2 p3 = do
124   let q1 = 0.5 *^ (p2 ^+^ p3)
125       q2 = 0.5 *^ (p3 ^+^ p1)
126       q3 = 0.5 *^ (p1 ^+^ p2)
127   triangle (not cw) (n-1) p1 q2 q3
128   triangle (not cw) (n-1) q1 p2 q3
129   triangle (not cw) (n-1) q1 q2 p3
130   triangle cw (n-1) q1 q2 q3
131
132 fillFigures3 :: (Color3 R, Color3 R, Color3 R) -> (Bool, Bool, Bool) -> [Figure3] -> IO ()
133 fillFigures3 cs v fs = renderPrimitive Triangles $ mapM_ (fillFigure3 cs v) fs
134
135 data State = State { sV :: Size, sVis :: (Bool, Bool, Bool), sP0 :: Point, sP :: Polytope, sFrame :: Int }
136
137 reshape :: IORef State -> Size -> IO ()
138 reshape r v = atomicModifyIORef r $ \s -> (s { sV = v }, ())
139
140 timer :: Timeout -> IORef State -> IO ()
141 timer t r = do
142   addTimerCallback t $ timer t r
143   postRedisplay Nothing
144
145 data Polytope = Polytope [M3] [V3]
146 data Point = Point V3
147
148 red :: Color3 R
149 red = Color3 ((0.114/0.299)**0.2) 0 0
150
151 green :: Color3 R
152 green = Color3 0 ((0.114/0.587)**0.2) 0
153
154 blue :: Color3 R
155 blue  = Color3 0 0 1
156
157 display :: IORef State -> IO ()
158 display ref = do
159   s <- atomicModifyIORef ref $ \x -> (x, x)
160   let vp@(Size vw vh) = sV s
161       w = 2 * fromIntegral vw / fromIntegral vh
162       h = 2
163   viewport $= (Position 0 0, vp)
164   matrixMode $= Projection
165   loadIdentity
166   ortho (-w) (w) (-h) (h) (-1) 1
167   matrixMode $= Modelview 0
168   loadIdentity
169   clear [ColorBuffer, DepthBuffer]
170   depthFunc $= Just Less
171   shadeModel $= Smooth
172   lineSmooth $= Enabled
173   lineWidth $= 2
174   let Point p = sP0 s
175       Polytope syms fr = sP s
176       a1 = fromIntegral (sFrame s) / 100
177       a2 = fromIntegral (sFrame s) / 141
178       a3 = fromIntegral (sFrame s) / 382
179       c1 = cos a1
180       s1 = sin a1
181       c2 = cos a2
182       s2 = sin a2
183       c3 = cos a3
184       s3 = sin a3
185       r1 = M3  c1 s1 0   (-s1) c1 0   0 0 1
186       r2 = M3  c2 0 s2   0 1 0   (-s2) 0 c2
187       r3 = M3  1 0 0   0 c3 s3   0 (-s3) c3
188       r = r1 ^^*^^ r2 ^^*^^ r3
189   fillFigures3 (red, green, blue) (sVis s) . transformsF3 r . polyhedronF syms . initialF3 fr . initialP3 fr $ p
190   flush
191   swapBuffers
192   atomicModifyIORef ref $ \s'' -> (s'' { sFrame = sFrame s'' + 1 }, ())
193
194 main :: IO ()
195 main = do
196   hSetBuffering stdin LineBuffering
197   let v = Size 1024 768
198   initialWindowSize $= v
199   initialDisplayMode $= [RGBAMode, DoubleBuffered, WithDepthBuffer]
200   (_, dims) <- getArgsAndInitialize
201   let [p,q] = map read dims
202       Just (file, size) = lookup (p,q) [((3,3),("3_3.3d", 24))
203                                        ,((4,3),("4_3.3d", 48))
204                                        ,((5,3),("5_3.3d", 120))]
205       n = 3*3 * size
206       b = n * sizeOf (undefined::R)
207   matrices <- fmap toMats $ withBinaryFile file ReadMode $ \h -> do
208    allocaBytes b $ \ptr -> do
209     hGetBuf h ptr b
210     peekArray n ptr
211   createWindow "reflex"
212   sr <- newIORef $ State { sV = v, sVis = (True, True, True), sP0 = Point (V3 1 1 1), sP = Polytope matrices (fr3 p q), sFrame = 500 }
213   forkIO . forever $ commands sr
214   reshapeCallback $= Just (reshape sr)
215   displayCallback $= display sr
216   idleCallback $= Just yield
217   timer 50 sr
218   mainLoop
219
220 toMats :: [R] -> [M3]
221 toMats (a:b:c:d:e:f:g:h:i:zz) = M3 a b c d e f g h i : toMats zz
222 toMats _ = []
223
224 commands :: IORef State -> IO ()
225 commands r = do
226   l <- getLine
227   case map read (words l) of
228     [la,lb,lc,va,vb,vc] -> atomicModifyIORef r $ \s -> (s { sVis = (va/=0,vb/=0,vc/=0), sP0 = Point (V3 la lb lc) }, ())
229     _ -> return ()
230   postRedisplay Nothing