MarchingSquaresFill.hs

-- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
-- Released under the GNU GPL, see LICENSE

module Graphics.Implicit.Export.MarchingSquaresFill (getContourMesh) where

import Graphics.Implicit.Definitions
import Control.Parallel (par, pseq)

-- | 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.

getContourMesh :: ℝ2 -> ℝ2 -> ℝ2 -> Obj2 -> [(ℝ2,ℝ2,ℝ2)]
getContourMesh (x1, y1) (x2, y2) (dx, dy) obj = 
	let
		-- How many steps will we take on each axis?
		nx = fromIntegral $ ceiling $ (x2 - x1) / dx
		ny = fromIntegral $ ceiling $ (y2 - y1) / dy
		-- Divide it up and compute the polylines
		trisOnGrid :: [[[(ℝ2,ℝ2,ℝ2)]]]
		trisOnGrid = [[getSquareTriangles
		           (x1 + (x2 - x1)*mx/nx,     y1 + (y2 - y1)*my/ny)
		           (x1 + (x2 - x1)*(mx+1)/nx, y1 + (y2 - y1)*(my+1)/ny)
		           obj
		     | mx <- [0.. nx-1] ] | my <- [0..ny-1] ]
		triangles = concat $ concat trisOnGrid
	in
		triangles
		

-- | This function gives line segmensts to divde 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.

getSquareTriangles :: ℝ2 -> ℝ2 -> Obj2 -> [(ℝ2,ℝ2,ℝ2)]
getSquareTriangles (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)

		square a b c d = [(a,b,c), (a,c,d)]

	in 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)  -> square (x1,y1) (x2,y1) (x2,y2) (x1,y2)
		(False, False,
		 False, False) -> []
		(True,  True, 
		 False, False) -> square midx1 midx2 (x2,y2) (x1,y2) 
		(False, False,
		 True,  True)  -> square (x1,y1) (x2,y1) midx2 midx1 
		(False, True, 
		 False, True)  -> square midy1 (x2,y1) (x2,y2) midy2
		(True,  False,
		 True,  False) -> square (x1,y1) midy1 midy2 (x1,y2)
		(True,  False,
		 False, False) -> [((x1,y2), midx1, midy2)]
		(False, True, 
		 True,  True)  -> 
			[(midx1, (x1,y1), midy2), ((x1,y1), (x2,y1), midy2), (midy2, (x2,y1), (x2,y2))]
		(True,  True, 
		 False, True)  -> 
			[((x1,y2), midx1, (x2,y2)), (midx1, midy1, (x2,y2)), ((x2,y2), midy1, (x2,y1))] 
		(False, False,
		 True,  False) -> [(midx1, (x1,y1), midy1)]
		(True,  True, 
		 True,  False) -> 
			[(midy1,midx2,(x2,y2)), ((x2,y2), (x1,y2), midy1), (midy1, (x1,y2), (x1,y1))]
		(False, False,
		 False, True)  -> [(midx2, midy1, (x2,y1))]
		(True,  False,
		 True,  True)  -> 
			[(midy2, (x2,y1), midx2), ((x2,y1), midy2, (x1,y1)), ((x1,y1), midy2, (x1,y2))]
		(False, True, 
		 False, False) -> [(midx2, (x2,y2), midy2)]
		(True,  False,
		 False, True)  -> if c > 0
			then [((x1,y2), midx1, midy2), ((x2,y1), midy1, midx2)]
			else [] --[[midx1, midy1], [midx2, midy2]]
		(False, True, 
		 True,  False) -> if c <= 0
			then [] --[[midx1, midy2], [midx2, midy1]]
			else [((x1,y1), midy1, midx1), ((x2,y2), midx2, midy2)] --[[midx1, midy1], [midx2, midy2]]