implicit-0.1.0: Graphics/Implicit/Export/MarchingSquares.hs
-- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
-- Copyright 2014 2015 2016, Julia Longtin (julial@turinglace.com)
-- Released under the GNU AGPLV3+, see LICENSE
module Graphics.Implicit.Export.MarchingSquares (getContour) where
import Prelude(Int, Bool(True, False), ceiling, fromIntegral, (/), (+), (-), filter, map, ($), (*), (/=), (<=), (>), (.), splitAt, div, unzip, length, (++), (<), (++), head, concat, not, null, (||), Eq, Int, fst, snd)
import Graphics.Implicit.Export.Render.HandlePolylines (reducePolyline)
import Graphics.Implicit.Definitions (ℝ2, Polyline, Obj2, (⋯/), (⋯*))
-- FIXME: commented out for now, parallelism is not properly implemented.
-- import Control.Parallel.Strategies (using, parList, rdeepseq)
import Data.VectorSpace ((^-^), (^+^))
both :: (a -> b) -> (a,a) -> (b,b)
both f (x,y) = (f x, f y)
-- | getContour gets a polyline describe the edge of your 2D
-- object. It's really the only function in this file you need
-- to care about from an external perspective.
getContour :: ℝ2 -> ℝ2 -> ℝ2 -> Obj2 -> [Polyline]
getContour p1 p2 d obj =
let
-- How many steps will we take on each axis?
n :: (Int, Int)
n = (ceiling) `both` ((p2 ^-^ p1) ⋯/ d)
nx = fst n
ny = snd n
-- Divide it up and compute the polylines
gridPos :: (Int,Int) -> (Int,Int) -> ℝ2
gridPos (nx',ny') (mx,my) =
let
p :: ℝ2
p = ( fromIntegral mx / fromIntegral nx'
, fromIntegral my / fromIntegral ny')
in
p1 ^+^ (p2 ^-^ p1) ⋯* p
linesOnGrid :: [[[Polyline]]]
linesOnGrid = [[getSquareLineSegs
(gridPos n (mx,my))
(gridPos n (mx+1,my+1))
obj
| mx <- [0.. nx-1] ] | my <- [0..ny-1] ]
-- Cleanup, cleanup, everybody cleanup!
-- (We connect multilines, delete redundant vertices on them, etc)
multilines = (filter polylineNotNull) $ (map reducePolyline) $ orderLinesDC $ linesOnGrid
in
multilines
-- FIXME: Commented out, not used?
{-
getContour2 :: ℝ2 -> ℝ2 -> ℝ2 -> Obj2 -> [Polyline]
getContour2 p1@(x1, y1) p2@(x2, y2) d obj =
let
-- How many steps will we take on each axis?
n@(nx,ny) = (fromIntegral . ceiling) `both` ((p2 ^-^ p1) ⋯/ d)
-- Grid mapping funcs
fromGrid (mx, my) = let p = (mx/nx, my/ny)
in (p1 ^+^ (p2 ^-^ p1) ⋯/ p)
toGrid (x,y) = (floor $ nx*(x-x1)/(x2-x1), floor $ ny*(y-y1)/(y2-y1))
-- Evaluate obj on a grid, in parallel.
valsOnGrid :: [[ℝ]]
valsOnGrid = [[ obj (fromGrid (mx, my)) | mx <- [0.. nx-1] ] | my <- [0..ny-1] ]
`using` parList rdeepseq
-- A faster version of the obj. Sort of like memoization, but done in advance, in parallel.
preEvaledObj p = valsOnGrid !! my !! mx where (mx,my) = toGrid p
-- Divide it up and compute the polylines
linesOnGrid :: [[[Polyline]]]
linesOnGrid = [[getSquareLineSegs (fromGrid (mx, my)) (fromGrid (mx+1, my+1)) preEvaledObj
| mx <- [0.. nx-1] ] | my <- [0..ny-1] ]
-- Cleanup, cleanup, everybody cleanup!
-- (We connect multilines, delete redundant vertices on them, etc)
multilines = (filter polylineNotNull) $ (map reducePolyline) $ orderLinesDC $ linesOnGrid
in
multilines
-}
-- | This function gives line segments to divide negative interior
-- regions and positive exterior ones inside a square, based on its
-- values at its vertices.
-- It is based on the linearly-interpolated marching squares algorithm.
getSquareLineSegs :: ℝ2 -> ℝ2 -> Obj2 -> [Polyline]
getSquareLineSegs (x1, y1) (x2, y2) obj =
let
(x,y) = (x1, y1)
-- Let's evlauate obj at a few points...
x1y1 = obj (x1, y1)
x2y1 = obj (x2, y1)
x1y2 = obj (x1, y2)
x2y2 = obj (x2, y2)
c = obj ((x1+x2)/2, (y1+y2)/2)
dx = x2 - x1
dy = y2 - y1
-- linearly interpolated midpoints on the relevant axis
-- midy2
-- _________*__________
-- | |
-- | |
-- | |
--midx1* * midx2
-- | |
-- | |
-- | |
-- ---------*----------
-- midy1
midx1 = (x, y + dy*x1y1/(x1y1-x1y2))
midx2 = (x + dx, y + dy*x2y1/(x2y1-x2y2))
midy1 = (x + dx*x1y1/(x1y1-x2y1), y )
midy2 = (x + dx*x1y2/(x1y2-x2y2), y + dy)
notPointLine :: Eq a => [a] -> Bool
notPointLine (p1:p2:[]) = p1 /= p2
notPointLine ([]) = False
notPointLine ([_]) = False
notPointLine (_ : (_ : (_ : _))) = False
in filter (notPointLine) $ case (x1y2 <= 0, x2y2 <= 0,
x1y1 <= 0, x2y1 <= 0) of
-- Yes, there's some symetries that could reduce the amount of code...
-- But I don't think they're worth exploiting...
(True, True,
True, True) -> []
(False, False,
False, False) -> []
(True, True,
False, False) -> [[midx1, midx2]]
(False, False,
True, True) -> [[midx1, midx2]]
(False, True,
False, True) -> [[midy1, midy2]]
(True, False,
True, False) -> [[midy1, midy2]]
(True, False,
False, False) -> [[midx1, midy2]]
(False, True,
True, True) -> [[midx1, midy2]]
(True, True,
False, True) -> [[midx1, midy1]]
(False, False,
True, False) -> [[midx1, midy1]]
(True, True,
True, False) -> [[midx2, midy1]]
(False, False,
False, True) -> [[midx2, midy1]]
(True, False,
True, True) -> [[midx2, midy2]]
(False, True,
False, False) -> [[midx2, midy2]]
(True, False,
False, True) -> if c > 0
then [[midx1, midy2], [midx2, midy1]]
else [[midx1, midy1], [midx2, midy2]]
(False, True,
True, False) -> if c <= 0
then [[midx1, midy2], [midx2, midy1]]
else [[midx1, midy1], [midx2, midy2]]
-- $ Functions for cleaning up the polylines
-- Many have multiple implementations as efficiency experiments.
-- At some point, we'll get rid of the redundant ones....
{-
orderLines :: [Polyline] -> [Polyline]
orderLines [] = []
orderLines (present:remaining) =
let
findNext ((p3:ps):segs) = if p3 == last present then (Just (p3:ps), segs) else
if last ps == last present then (Just (reverse $ p3:ps), segs) else
case findNext segs of (res1,res2) -> (res1,(p3:ps):res2)
findNext [] = (Nothing, [])
in
case findNext remaining of
(Nothing, _) -> present:(orderLines remaining)
(Just match, others) -> orderLines $ (present ++ tail match): others
-}
orderLinesDC :: [[[Polyline]]] -> [Polyline]
orderLinesDC segs =
let
halve :: [a] -> ([a], [a])
halve l = splitAt (div (length l) 2) l
splitOrder segs' = case (\(x,y) -> (halve x, halve y)) . unzip . map (halve) $ segs' of
((a,b),(c,d)) -> orderLinesDC a ++ orderLinesDC b ++ orderLinesDC c ++ orderLinesDC d
in
if (length segs < 5 || length (head segs) < 5 ) then concat $ concat segs else
splitOrder segs
{-
orderLinesP :: [[[Polyline]]] -> [Polyline]
orderLinesP segs =
let
halve l = splitAt (div (length l) 2) l
splitOrder segs = case (\(x,y) -> (halve x, halve y)) $ unzip $ map (halve) segs of
((a,b),(c,d)) -> orderLinesDC a ++ orderLinesDC b ++ orderLinesDC c ++ orderLinesDC d
-- force is frome real world haskell
force xs = go xs `pseq` ()
where go (_:xs) = go xs
go [] = 1
in
if (length segs < 5 || length (head segs) < 5 ) then concat $ concat segs else
case (\(x,y) -> (halve x, halve y)) $ unzip $ map (halve) segs of
((a,b),(c,d)) -> orderLines $
let
a' = orderLinesP a
b' = orderLinesP b
c' = orderLinesP c
d' = orderLinesP d
in (force a' `par` force b' `par` force c' `par` force d') `pseq`
(a' ++ b' ++ c' ++ d')
-}
polylineNotNull :: [a] -> Bool
polylineNotNull (_:l) = not (null l)
polylineNotNull [] = False