------------------------- --RASTERIZER ------------------ {- each point satisfies the following relationship, leaving us to find a, b, c as scalar multiples of the 'northern' vector, 'eastwards' vector and 'view' vector, respectively. / n_x e_x v_x \ / a \ / dx \ | n_y e_y v_y | | b | = | dy | \ n_z e_z v_z / \ c / \ dz / -} import LSystem import List --vectors are arrays to allow for easy manipulation and interaction with matrices type Point3D = (Float, Float, Float) type Point2D = (Float, Float) type Vector3D = [Float] --3x type Vector2D = [Float] --2x type Vector = [Float] --Nx type Matrix = [[Float]] data Object3D = Object3D { points3D :: [Point3D] , color3D :: Pen } deriving (Show, Eq) data Object2D = Object2D { points2D :: [Point2D] , distance2D :: Float --preserve some idea of original 3D object's distance from view point , color2D :: Pen } deriving (Show, Eq) --HELPER VECTOR/MATRIX/POINT FUNCTIONS --inverts a 3x3 matrix invMat :: Matrix -> Matrix invMat mat = divdet [ [ i * e - h * f, h * c - i * b, f * b - e * c ] , [ g * f - i * d, i * a - g * c, d * c - f * a ] , [ h * d - g * e, g * b - h * a, e * a - d * b ] ] where [ [ a, b, c ] , [ d, e, f ] , [ g, h, i ] ] = mat divdet m = map (\line -> map (/det) line) m det = a * (e * i - f * h) - d * (i * b - h * c) + g * (f * b - e * c) --generates a matrix with x, y, z mappings of vectors v1, v2, v3 respectively genMat :: Vector -> Vector -> Vector -> Matrix genMat v1 v2 v3 = [ [a, b, c] | (a, b, c) <- zip3 v1 v2 v3 ] --dots two vectors (any dimension) dotProduct :: Vector -> Vector -> Float dotProduct vec1 vec2 = sum [ a * b | (a, b) <- zip vec1 vec2 ] --crosses two vectors (3rd dimension) crossProduct :: Vector3D -> Vector3D -> Vector3D crossProduct [a, b, c] [x, y, z] = [ c * y - b * z, a * z - c * x, b * x - a * y ] --multiplies matrix mat against vector vec, as a vector function fMat :: Matrix -> Vector -> Vector fMat mat vec = [ dotProduct line vec | line <- mat ] --length of a 3D vector vecLen :: Vector3D -> Float vecLen [x, y, z] = sqrt (x ^ 2 + y ^ 2 + z ^ 2) --translates a 3D vector into a unit vector of the same direction vecUnit :: Vector -> Vector vecUnit v = map (/ vecLen v) v --distance between two points (max 3rd dimension) lineLen :: Point3D -> Point3D -> Float lineLen (a, b, c) (x, y, z) = sqrt ((x - a) ^ 2 + (y - b) ^ 2 + (z - c) ^ 2) lineLen2D :: Point2D -> Point2D -> Float lineLen2D (a, b) (x, y) = lineLen (a, b, 0) (x, y, 0) --RASTERIZER --maps a collection of points in 3D space onto a 2D plane based on the parameters passed mapPoints :: Point3D -> Vector3D -> Vector3D -> Object3D -> Object2D mapPoints viewPoint viewVec northVec (Object3D {points3D = points, color3D = color}) = Object2D { points2D = [ (x' * zZoom / z', y' * zZoom / z') | (x, y, z) <- points, let [x', y', z'] = fMat mat [x - xp, y - yp, z - zp] ] , distance2D = average $ map (lineLen viewPoint) points --this is minimum for reasons of rendering order , color2D = color } where mat = invMat (genMat eastVec' northVec' viewVec') --unit vector equivelents of the three view-defining vector parameters northVec' = vecUnit northVec viewVec' = vecUnit viewVec eastVec' = vecUnit $ crossProduct viewVec northVec (xp, yp, zp) = viewPoint zZoom = 10 --this can pretty much be arbitrary thanks to the auto-zoom feature of LSystem.h average xs = (sum xs) / (fromIntegral $ length xs) --a viewpoint facing towards the origin with a given radius and rotation --rotZ = rotation around z-axis; rotUpDown = 'pitch' relative to xy plane inwardView :: Float -> (Float, Float) -> Object3D -> Object2D inwardView radius (rotZ, rotUpDown) obj = mapPoints (x, y, z) [-x, -y, -z] [-sinZ * sinU, -cosZ * sinU, cosU] obj where sinZ = sin rotZ cosZ = cos rotZ sinU = sin rotUpDown cosU = cos rotUpDown x = radius * sinZ * cosU y = radius * cosZ * cosU z = radius * sinU --TRACE -> TURTLE TRANSLATOR --sorting function, based off callback --Ord does not play a role due to the presence of a callback being implicit testamony to 'orderability' qSort :: (a -> a -> Int) -> [a] -> [a] qSort callback [] = [] qSort callback [x] = [x] qSort callback (x:xs) = qSort callback [ y | y <- xs, callback y x <= 0 ] ++ [x] ++ qSort callback [ y | y <- xs, callback y x > 0 ] --retrieve aspects of an Object2D structure getPoints :: Object2D -> [Point2D] getPoints (Object2D {points2D = points}) = points getColor :: Float -> Object2D -> Pen getColor scaleFactor (Object2D {distance2D = scale, color2D = color}) = let (Colour r g b) = color in Colour (r * scale') (g * scale') (b * scale') where scale' = scale / scaleFactor getDistance :: Object2D -> Float getDistance (Object2D {distance2D = distance}) = distance --helper function for points2Turtle; writes the turtle commands to move from point (a, b) to (x, y) --preserves orientation, which is assumed to be northwards pointPair2Turtle :: (Point2D, Point2D) -> Command pointPair2Turtle (p1, p2) = Turn delta :#: Go len :#: Turn (-delta) where delta = makeAngle p1 p2 * 180 / pi - 90 len = lineLen2D p1 p2 --translates a trace to turtle commands points2Turtle :: [Point2D] -> Command points2Turtle points = joinTurtle . optimizeTurtle . splitTurtle . joinTurtle $ map pointPair2Turtle $ zip points (tail points) --translates an array of 2D traces into a single string of turtle graphics commands objects2Turtle :: Float -> [Object2D] -> Command objects2Turtle scaleFactor os = GrabPen Inkless :#: (joinTurtle $ map doObject os) --(map mergePair $ zip os (tail os)) -- ++ [points2Turtle . getPoints $ last os] --the 'glue' function where doObject obj = Branch ( (pointPair2Turtle ((0, 0), head $ getPoints obj)) :#: (GrabPen $ getColor scaleFactor obj) :#: (points2Turtle $ getPoints obj) ) {-mergePair (obj1, obj2) = points2Turtle (getPoints obj1) :#: GrabPen Inkless :#: pointPair2Turtle (last $ getPoints obj1, head $ getPoints obj2) :#: GrabPen (getColor scaleFactor obj2)-} --translates an array of 3D traces mapped with inwardView into turtle commands showObjects :: Float -> (Float, Float) -> [Object3D] -> Command showObjects radius rot objs = --GrabPen (getColor scaleFactor $ head objs2D) :#: objects2Turtle scaleFactor objs2D objects2Turtle scaleFactor objs2D --we need to now only map the objects to a 2D plane but sort the 2D objects according to distance from viewpoint, --so as to properly render them where objs2D = sortObjs $ map (inwardView radius rot) objs sortObjs os = qSort sortCallback os sortCallback o1 o2 = if len1 == len2 then 0 else if len1 > len2 then -1 else 1 where (len1, len2) = (getDistance o1, getDistance o2) scaleFactor = maximum $ map getDistance objs2D --bla joinTurtle :: [Command] -> Command joinTurtle cs = foldr (:#:) Sit cs splitTurtle :: Command -> [Command] splitTurtle (a :#: b) = splitTurtle a ++ splitTurtle b splitTurtle Sit = [] splitTurtle c = [c] --deals only with lists of commands optimizeTurtle :: [Command] -> [Command] optimizeTurtle c | optc == c = c | otherwise = optimizeTurtle optc where optc = optimizeTurtle' c --the work-turtle optimizeTurtle' (Go x : Go y : cs) = optimizeTurtle' (Go (x + y) : cs) optimizeTurtle' (Turn x : Turn y : cs) = optimizeTurtle' (Turn (x + y) : cs) optimizeTurtle' (Sit : cs) = cs optimizeTurtle' (Go 0 : cs) = cs optimizeTurtle' (Turn 0 : cs) = cs optimizeTurtle' (c:cs) = c : optimizeTurtle' cs optimizeTurtle' [] = [] --MENGER SPONGE GENERATION --p = center of cube of side length n --this needs a bit more attention because it must represent the uninterrupted trace of a cube, preferrably with the least amount --of retracing possible cube :: Point3D -> Float -> [Point3D] cube (px, py, pz) n = [ --two full sides (bx, by, bz), (bx + n, by, bz), (bx + n, by, bz + n), (bx, by, bz + n), (bx, by, bz), (bx, by + n, bz), (bx + n, by + n, bz), (bx + n, by + n, bz + n), (bx, by + n, bz + n), (bx, by + n, bz), --the unattended edges + some retracing (bx + n, by + n, bz), (bx + n, by, bz), (bx + n, by, bz + n), (bx + n, by + n, bz + n), (bx, by + n, bz + n), (bx, by, bz + n) ] where (bx, by, bz) = (px - n', py - n', pz - n') n' = n / 2 --creates a 2D lambda in plane parallel to x-z plane lambda :: Point3D -> Float -> [Point3D] lambda (px, py, pz) n = [ (px - n', py, pz + n'), (px + n', py, pz - n'), (px, py, pz), (px - n', py, pz - n') ] where n' = n / 2 --purdy culurs spongePalette :: [Pen] --spongePalette = [Colour 0.5 0.5 1, Colour 1 0.5 0.5, Colour 1 1.5 1.5] spongePalette = [Colour 0.5 0.5 1, Colour 1 0.8 0.8, Colour 1.3 1.0 1.3] --it's an lsystem at heart... sponge :: Int -> Float -> [Object3D] sponge i n = f i (0, 0, 0) where f 0 p = [makeObject 0 $ cube p n'] -- , makeObject 0 $ lambda' p n'] where n' = n / 3 ^ (i + 1) f (x + 1) p@(px, py, pz) = concat [ f x (px + a, py + b, pz + c) | a <- [-n', 0, n'], b <- [-n', 0, n'], c <- [-n', 0, n'], --the absolute and planar centers must not be covered (a /= 0 && b /= 0 && c == 0) || (a /= 0 && b == 0 && c /= 0) || (a == 0 && b /= 0 && c /= 0) || (a /= 0 && b /= 0 && c /= 0) ] ++ [makeObject (x + 1) $ cube p n'] -- , makeObject (x + 1) $ lambda' p n'] where n' = n / 3 ^ (i - x) makeObject ind points = Object3D { points3D = points, color3D = spongePalette !! (i - ind) } lambda' (px, py, pz) n = lambda (px, py + n / 2, pz) (n / 2) --------------------------------------- --TIME TO ROLL OUR OWN REGION CLIPPING SYSTEM! TEH WOOTZAWRS! ------------------------------ --yes, this is all rather ugly. sorry about that. type Polygon = [Point2D] data LineEq = LineEq Float Float Float deriving (Show) data Intersect a = NoIntersect | Intersect a | InfiniteIntersect deriving (Show, Eq) --uses an intersect algorithm; regrettably not my own idea --however I pardon myself for this heresy by the consideration that I did make my own algorithm/code for this, but unfortunately --is potentially less efficient for concave polygons --and it's longer --and it dosn't fit in as well with the other stuff pointInPolygon :: Point2D -> Polygon -> Bool pointInPolygon (x, y) points = mod (length $ filter (\line -> let (a, b) = lineIntersect' (LineEq 0 1 y) line in a < x) $ map makeLine relavantSides) 2 == 1 where sides = zip points (tail points ++ [head points]) relavantSides = filter (\((_, y1), (_, y2)) -> y `between` (y1, y2)) sides makeAngle :: Point2D -> Point2D -> Float makeAngle (a, b) (x, y) | dx == 0 = if dy == 0 then 0 else if dy > 0 then pi / 2 else 3 * pi / 2 | dx < 0 = pi + tn | dy < 0 = 2 * pi + tn | otherwise = tn where dx = x - a dy = y - b tn = atan (dy / dx) --traces out the outer-most points into a convex polygon convexPolygon :: Polygon -> Polygon convexPolygon points = p2CP' pStart pStart $ nub points where pStart = foldr1 (\(x, y) (a, b) -> if y < b || (y == b && x < a) then (x, y) else (a, b)) points p2CP' pOrig pCur [] = [pCur] p2CP' pOrig pCur pList | pNext == pOrig = [pCur] | otherwise = pCur : p2CP' pOrig pNext (pList \\ [pNext]) --find the next point which generates the line with the least angle with this where (_, pNext) = foldr1 foldrCmp assocAngles assocAngles = map (\p -> (positiveAngle $ makeAngle pCur p - makeAngle pOrig pCur, p)) (pList \\ [pCur]) positiveAngle x | x < 0 = 2 * pi + x | otherwise = x --essentially a maximum function which preserves extra information foldrCmp v1@(a1, p1) v2@(a2, p2) | a1 < a2 = v1 | a1 > a2 = v2 | otherwise = if lineLen2D pCur p1 < lineLen2D pCur p2 then v1 else v2 convexObject :: Object2D -> Object2D convexObject (Object2D { points2D = points, color2D = color, distance2D = distance }) = Object2D { points2D = convexPolygon points, color2D = color, distance2D = distance } linesBetween :: [a] -> [(a, a)] linesBetween xs = zip xs (tail xs ++ [head xs]) between :: Float -> (Float, Float) -> Bool a `between` (x, y) | x > y = x > a && a > y --you have NO IDEA how much trouble >= <= caused me until I found them... | otherwise = y > a && a > x pointBetween :: Point2D -> (Point2D, Point2D) -> Bool pointBetween (x, y) ((a, b), (c, d)) = x `between` (a, c) && y `between` (b, d) --returns a list (side index, point intersect) --optimize? linePolygonIntersects :: (Point2D, Point2D) -> Polygon -> [(Int, Point2D)] linePolygonIntersects ppair points = filter (\(id, p) -> pointBetween p ppair && pointBetween p (lines !! id)) intersects where intersects = [ (i, p) | (i, int) <- zip [0..] $ map ((lineIntersect qline).makeLine) lines, int /= NoIntersect && int /= InfiniteIntersect, let Intersect p = int ] lines = linesBetween points qline = makeLine ppair --possible superfluous... --but whatever makeLine :: (Point2D, Point2D) -> LineEq makeLine ((a, b), (x, y)) | dx == 0 = LineEq dy dx x | otherwise = LineEq dy dx (y - dy / dx * x) where dx = x - a dy = y - b --for cases which we are guaranteed to get an intersection point lineIntersect' :: LineEq -> LineEq -> Point2D lineIntersect' lineA lineB = case lineIntersect lineA lineB of Intersect p -> p _ -> (0, 0) filterIntersects :: [Intersect Point2D] -> [Point2D] filterIntersects is = [ p | i <- is, i /= NoIntersect && i /= InfiniteIntersect, let Intersect p = i ] lineIntersect :: LineEq -> LineEq -> Intersect Point2D lineIntersect (LineEq dy1 dx1 c1) (LineEq dy2 dx2 c2) --are any of the lines vertical? | dx1 == 0 && dx2 == 0 = NoIntersect | dx1 == 0 = Intersect (c1, dy2 / dx2 * c1 + c2) | dx2 == 0 = Intersect (c2, dy1 / dx1 * c2 + c1) --we're free to divide by dx --unless of course we're chuck norris, in which case we were never restricted by silly divide-by-0 cases anyhow | dy1 / dx1 == dy2 / dx2 = InfiniteIntersect | otherwise = Intersect (x, y) where x = (c1 - c2) / (dy2 / dx2 - dy1 / dx1) y = dy1 / dx1 * x + c1 --bounding rectangles polygonBoundaries :: Polygon -> (Point2D, Point2D) polygonBoundaries poly = foldr f ((px, py), (px, py)) poly where (px, py) = head poly f (x, y) ((xx, xy), (nx, ny)) = ( ( max x xx, max y xy ), ( min x nx, min y ny ) ) boundariesIntersect :: (Point2D, Point2D) -> (Point2D, Point2D) -> Bool boundariesIntersect a@((pxx, pxy), (pnx, pny)) b@((qxx, qxy), (qnx, qny)) = pointBetween (pxx, pxy) b || pointBetween (pxx, pny) b || pointBetween (pnx, pxy) b || pointBetween (pnx, pny) b || pointBetween (qxx, qxy) a || pointBetween (qxx, qny) a || pointBetween (qnx, qxy) a || pointBetween (qnx, qny) a --bow down before my mighty naive homemade algorithms! --which are buggy showObjectsClipped :: Float -> [Polygon] -> [Object2D] -> Command showObjectsClipped scaleFactor regions [] = Sit showObjectsClipped scaleFactor regions (o:os) = if startState == 0 && any (nestedPolygon opoints) relRegions then Sit else Branch ( GrabPen Inkless :#: (pointPair2Turtle ((0, 0), (head opoints))) :#: GrabPen (transState startState) :#: (showObject relRegions startState sides) ) :#: (showObjectsClipped scaleFactor regions os) where sides = linesBetween opoints opoints = getPoints o --restrict our search of regions only to those which are of concern bounds = polygonBoundaries opoints relRegions = filter ((boundariesIntersect bounds).polygonBoundaries) regions --is our first point hidden? startState = if or [ pointInPolygon (head opoints) region | region <- relRegions ] then 0 else 1 --translates a given line state into an actual color command transState ind = states !! (mod ind $ length states) where states = [Inkless, getColor scaleFactor o] --progresses along each line of an object creating the turtle commands to properly render them --renders each object 'around' the specified regions showObject regions curState [] = Sit showObject regions curState ((p1, p2):ps) = Turn delta :#: (joinTurtle $ map (\((cmd, p1), (_, p2)) -> cmd :#: Go (lineLen2D p1 p2)) $ zip points (tail points)) :#: Turn (-delta) :#: (showObject regions (curState + length poi') ps) where delta = makeAngle p1 p2 * 180 / pi - 90 --each point has an associated command to switch between colored/inkless as necessary cmds = [ GrabPen $ transState c | c <- [curState + 1, curState + 2..] ] poi = concat [ [ p | (_, p) <- linePolygonIntersects (p1, p2) region ] | region <- regions ] poi' = zip cmds $ qSort (\pa pb -> if lineLen2D p1 pa < lineLen2D p1 pb then -1 else 1) poi points = [(Sit, p1)] ++ poi' ++ [(Sit, p2)] --analagous to object2Turtle, except with clipping objects2Turtle' :: Float -> [Polygon] -> [Object2D] -> Command objects2Turtle' scaleFactor regions [] = Sit objects2Turtle' scaleFactor regions (o:os) = showObjectsClipped scaleFactor regions [o] :#: objects2Turtle' scaleFactor regions' os where regions' = mergePolygonCollection points regions points = convexPolygon (getPoints o) -- ^ what he said showObjects' :: Float -> (Float, Float) -> [Object3D] -> Command showObjects' radius rot objs = objects2Turtle' scaleFactor [] objs2D where objs2D = sortObjs $ map (convexObject.(inwardView radius rot)) objs sortObjs os = qSort sortCallback os sortCallback o1 o2 = if len1 == len2 then 0 else if len1 > len2 then 1 else -1 where (len1, len2) = (getDistance o1, getDistance o2) scaleFactor = maximum $ map getDistance objs2D --merges a polygon into a collection of polygons mergePolygonCollection :: Polygon -> [Polygon] -> [Polygon] mergePolygonCollection p ps = ps' where relps = filter ((boundariesIntersect (polygonBoundaries p)).polygonBoundaries) ps mzps = zip relps $ map (\q -> mergePolygons q p) relps stay = [ r | (r, r') <- mzps, r == r' ] change = [ r | (r, r') <- mzps, r /= r' ] ps' = ( if null mzps then [p] --no relavant ones, outside all of them else if null change --nothings changes then ( if all (nestedPolygon p) relps then [] else [p] ) else [ foldl1 mergePolygons (p:change) ] ) ++ stay ++ (ps \\ relps) --is p inside q? nestedPolygon p q = and [ pointInPolygon point q | point <- p ] --is poly traced out in the clockwise (-1) or counterclockwise (1) direction? polygonDirection :: Polygon -> Int polygonDirection poly = if makeAngle (poly !! ind) (poly !! nind) < makeAngle (poly !! ind) (poly !! pind) then -1 else 1 where lowest = foldl1 (\(a, b) (x, y) -> if y < b || y == b && (a < x) then (x, y) else (a, b)) poly ind = head [ i | (i, p) <- zip [0..] poly, p == lowest ] nind = mod (ind + 1) (length poly) pind = mod (ind - 1 + length poly) (length poly) --traces the outer-most union of two polygons --todo: deal with spaces in between --hotodoit: fudge each 'seperate' space as a thinly (read: zero-width) connected portion of a concave polygon --hotodoit2.0: blindly ignore the possibility of this occuring, probably resulting the bug which I'm -- tearing my hair out trying to find. in fact, now that I think about it, I'm quite certain it is: in -- larger sponges, such situations /do/ arise, and in some cases it will trace out the /inner/ space -- rather than the outer one, depending on where it happens to begin mergePolygons :: Polygon -> Polygon -> Polygon mergePolygons polyA polyB | null firstpoints = polyB | otherwise = mergePolygons' 0 ([], []) (points', linesA) (polyB, linesB) where --our starting point must obviously land outside the other polygon firstpoints = dropWhile (\p -> pointInPolygon p polyB) polyA points = firstpoints ++ takeWhile (\p -> pointInPolygon p polyB) polyA --we must make sure that in jumping from one to another we keep the same angular direction points' = if polygonDirection polyA == polygonDirection polyB then points else reverse (tail points ++ [head points]) linesA = linesBetween points' linesB = linesBetween polyB mergePolygons' curA (doneA, doneB) a@(polyA, linesA) b@(polyB, linesB) | curA `elem` doneA = [] | (not . null) poi = (polyA !! curA) : ip : mergePolygons' (incId bid polyB) (doneB, curA : doneA) b a --swap | otherwise = (polyA !! curA) : mergePolygons' (incId curA polyA) (curA : doneA, doneB) a b where (bid, ip) = head poi poi = qSort (\(_, a) (_, b) -> if lineLen2D curp a < lineLen2D curp b then -1 else 1) $ linePolygonIntersects curline polyB (curp, _) = curline curline = linesA !! curA incId id poly = mod (id + 1) $ length poly --note to self: send to w.b.heijltjes@sms.ed.ac.uk main = display $ showObjects 30 (pi / 10, pi / 14) $ sponge 2 27 main' = display $ showObjects' 30 (pi / 10, pi / 14) $ sponge 2 27