diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Changelog + +## [0.1] - 2018-03-20 + +[0.1]: https://github.com/dzhus/csg/tree/0.1 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2012-2018, Dmitry Dzhus + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Dmitry Dzhus nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,140 @@ +# CSG: constructive solid geometry library + +[![Travis CI build status](https://travis-ci.org/dzhus/csg.svg)](https://travis-ci.org/dzhus/csg) +[![Hackage](https://img.shields.io/hackage/v/csg.svg?colorB=5e5184&style=flat)](https://hackage.haskell.org/package/csg) +[![Hackage deps](https://img.shields.io/hackage-deps/v/csg.svg)](http://packdeps.haskellers.com/feed?needle=csg) + +CSG is a [constructive solid geometry][csg-wiki] library with support +for ray casting. CSG allows you to define a complex solid as a +composition of primitives. It also provides functions to perform ray +casting (find an intersection of a ray and the defined solid) or test +whether a point belongs to the solid (for Monte Carlo volume +calculation). + +```haskell +-- "Data.CSG" uses 'Vec3' to represent vectors and points: +>>> let p1 = fromXYZ (5, -6.5, -5) +>>> toXYZ (origin :: Point) +(0.0,0.0,0.0) + +-- Define some solids: +>>> let s = sphere origin 5.0 +>>> let b = cuboid (fromXYZ (-1, -1, -1)) (fromXYZ (1, 1, 1)) + +-- Test if a point is inside the solid: +>>> origin `inside` (s `intersect` b) +True + +>>> origin `inside` (s `subtract` b) +False + +-- Find the distance to the next intersection of a ray with a solid, along with the +-- surface normal: +>>> let axis = fromXYZ (1, 2, 10) +>>> let solid = cylinder origin axis 2.0 `intersect` sphere origin 3.5 +>>> let ray = Ray (p1, origin <-> p1) +>>> ray `cast` solid +Just (HitPoint 0.7422558525331708 (Just (CVec3 0.7155468474912454 (-0.6952955216188516) 6.750441957464598e-2))) + +-- Load a solid definition from a file: +>>> import Data.CSG.Parser +>>> Right solid2 <- parseGeometryFile "examples/reentry.geo" +>>> ray `cast` solid2 +Just (HitPoint 10.877824491509912 (Just (CVec3 (-0.5690708596937849) 0.7397921176019203 0.3589790793088691))) +``` + +Please consult the [Hackage page for csg][hackage-doc] +for full documentation. + +By default `csg` is built using `CVec3` from [simple-vec3][] to +represent vectors and points, which according to benchmarks shows +better performance with Unboxed and Storable vectors. Build `csg` with +`triples` flag to use `(Double, Double, Double)` instead which may be +a more convenient programmatic interface that needs no +`fromXYZ`/`toXYZ`. + +See [alternatives](#alternatives) too. + +## csg-raycaster + +The package also includes `csg-raycaster` executable, which is a +simple interactive GUI for the ray casting algorithm. + +`csg-raycaster` takes a geometry defintion file as input. See +[`cube.geo`](examples/cube.geo): + +``` +solid box = orthobrick (-150, -150, -150; 150, 150, 150); + +solid rounded = sphere (0, 0, 0; 200); + +solid roundedbox = rounded and box; + +solid cylinder1 = cylinder (-160, 0, 0; 160, 0, 0; 100); +solid cylinder2 = cylinder (0, -160, 0; 0, 160, 0; 100); +solid cylinder3 = cylinder (0, 0, -160; 0, 0, 160; 100); + +solid cross = cylinder1 or cylinder2 or cylinder3; + +solid cutout = not cross; + +solid top = roundedbox and cutout; + +tlo top; +``` + +Please consult the [Hackage page for Data.CSG.Parser][parser-doc] for +full format specification. + +`csg-raycaster` may be run as + +``` +csg-raycaster cube.geo +``` + +Run as `csg-raycaster --help` to see all options. + +![csg-raycaster demo](csg-raycaster.gif) + +When run without a file argument, `csg-raycaster` will try to display +an arbitrary CSG solid. + +In the GUI window the following controls are supported: + +| Input | Function | +|---------------------------|--------------------------------------| +| Left mouse button + drag | Rotate | +| Right mouse button + drag | Pan | +| Mouse wheel up | Zoom in | +| Mouse wheel down | Zoom out | +| `r` | Reset zoom level and camera position | + +## Alternatives + +csg library performs no surface interpolation when doing ray casting. +Instead, we only solve ray-surface intersection equation numerically. + +There're other Haskell libraries for CSG: + +- [implicit][]: + + - Offers a much richer operation set + + - Uses function representation for CSG solids + + - If `implicit` had ray-casting support in early 2012 then I + probably wouldn't write `csg`. + +- [mecha][]: + + - Only provides types and functions to define solids and export + definitions to external formats + + - No support for ray casting + +[csg-wiki]: https://en.wikipedia.org/wiki/Constructive_solid_geometry +[hackage-doc]: http://hackage.haskell.org/package/csg/docs/Data-CSG.html +[implicit]: https://hackage.haskell.org/package/implicit +[mecha]: https://hackage.haskell.org/package/mecha +[parser-doc]: http://hackage.haskell.org/package/csg/docs/Data-CSG-Parser.html +[simple-vec3]: https://hackage.haskell.org/package/simple-vec3 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmark/Benchmark.hs b/benchmark/Benchmark.hs new file mode 100644 --- /dev/null +++ b/benchmark/Benchmark.hs @@ -0,0 +1,117 @@ +{-# LANGUAGE BangPatterns #-} + +{-| + +Single-threaded benchmark for CSG operations. + +-} + +import Criterion.Main +import Data.Vector.Unboxed as V hiding ((++)) +import qualified Data.Strict.Maybe as S + +import Data.CSG hiding (distance) + +-- | Pixels in meter at unit distance. +resolution :: Double +resolution = 50.0 + + +-- | Build cartesian axes from yaw and pitch with 0 roll. Angles are +-- in radians. +buildCartesian :: Double -> Double -> (Vec3, Vec3, Vec3) +buildCartesian yaw pitch = (u, v, w) + where u = fromXYZ (cos yaw * cos pitch, sin yaw * cos pitch, sin pitch) + v = fromXYZ (- (sin yaw), cos yaw, 0) + w = u >< v +{-# INLINE buildCartesian #-} + + +-- | Generate initial points of rays within a square region on a plane +-- parallel to YZ. +generateRayPoints :: Int + -- ^ Total ray count. + -> Double + -- ^ Distance of a plane from the origin. + -> V.Vector Vec3 +generateRayPoints rayCount distance = + let + !(n, sX, sY) = buildCartesian 0 0 + !p = n .^ (-distance) + -- Find a dimension of a square viewport used to generate + -- rays. The dimension will most closely fit the required + -- particle count number without exceeding it + dim :: Int + !dim = floor (sqrt $ fromIntegral rayCount :: Double) + !halfDim = dim `div` 2 + -- This differs from Raycaster module by a (fromIntegral dim) + -- term because range for ray index values are integers (as + -- opposed to (-1,1) used by gloss) + !scale = fromIntegral halfDim * distance / + (resolution * fromIntegral dim) + -- Initial point of I-th ray on a plane, row-major. 0 ray + -- starts at (-halfDim, -halfDim). + ithRay i = + p <+> (sX .^ (scale * rx)) + <+> (sY .^ (scale * ry)) + where + (y, x) = i `divMod` dim + rx = fromIntegral $ x - halfDim + ry = fromIntegral $ y - halfDim + {-# INLINE ithRay #-} + in + V.generate (dim * dim) ithRay + + +-- | Test raycasting performance for a solid. +test :: V.Vector Vec3 + -- ^ Initial points of test rays. Test rays are directed along + -- the Ox axis. + -> Solid + -> V.Vector Bool +test rayPoints solid = + let + !(n, _, _) = buildCartesian 0 0 + posToTrace pos = + -- Only head of trace gets evaluated + case cast (Ray (pos, n)) solid of + S.Nothing -> False + _ -> True + in + V.map posToTrace rayPoints + + +rays :: Int +rays = 1000 * 1000 + +dist :: Double +dist = 100 + + +solid1 :: Solid +solid1 = plane origin (fromXYZ (0, -0.5, 1)) + `intersect` + sphere origin 2.5 + + +solid2 :: Solid +solid2 = cylinder origin (fromXYZ (1, 0, 0)) 4 + `intersect` + cylinder origin (fromXYZ (0, 1, 0)) 4 + `intersect` + cylinder origin (fromXYZ (0, 0, 1)) 4 + + +main :: IO () +main = defaultMain + [ bgroup "Misc" + [ bench ("Test rays generation (" ++ show rays ++ ")" ) $ + whnf (uncurry generateRayPoints) (rays, dist) + ] + , bgroup "Raycasting" + [ bench "solid1" $ whnf (uncurry test) (rs, solid1) + , bench "solid2" $ whnf (uncurry test) (rs, solid2) + ] + ] + where + !rs = generateRayPoints rays dist diff --git a/csg.cabal b/csg.cabal new file mode 100644 --- /dev/null +++ b/csg.cabal @@ -0,0 +1,123 @@ +-- This file has been generated from package.yaml by hpack version 0.20.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: f2c1fd60fb41426a33ebda6ab5c083364c788ebc2458967a297173ddb3acb0b0 + +name: csg +version: 0.1 +synopsis: Analytical CSG (Constructive Solid Geometry) library +category: Graphics +homepage: https://github.com/dzhus/csg#readme +bug-reports: https://github.com/dzhus/csg/issues +author: Dmitry Dzhus +maintainer: dima@dzhus.org +license: BSD3 +license-file: LICENSE +build-type: Simple +cabal-version: >= 1.10 + +extra-source-files: + CHANGELOG.md + examples/cube.geo + examples/reentry.geo + README.md + +source-repository head + type: git + location: https://github.com/dzhus/csg + +flag triples + description: Use triples of Doubles to represent vectors (slower with vector library arrays) + manual: True + default: False + +library + hs-source-dirs: + src + ghc-options: -Wall -Wcompat -O2 + build-depends: + QuickCheck + , attoparsec + , base <5 + , bytestring + , containers + , simple-vec3 >=0.4 + , strict + , transformers + if flag(triples) + cpp-options: -DWITH_TRIPLES + exposed-modules: + Data.CSG + Data.CSG.Parser + other-modules: + Paths_csg + default-language: Haskell2010 + +executable csg-raycaster + main-is: raycaster.hs + hs-source-dirs: + exe + ghc-options: -Wall -Wcompat -O2 -threaded -rtsopts "-with-rtsopts=-N" + build-depends: + QuickCheck + , base <5 + , csg + , gloss + , gloss-raster + , simple-vec3 >=0.4 + , strict + , system-filepath + , turtle + other-modules: + Paths_csg + default-language: Haskell2010 + +test-suite csg-doctests + type: exitcode-stdio-1.0 + main-is: doctest-driver.hs + hs-source-dirs: + tests + ghc-options: -Wall -Wcompat -O2 -threaded + build-depends: + base <5 + , doctest + , doctest-discover >=0.1.0.8 + , simple-vec3 >=0.4 + other-modules: + Main + Paths_csg + default-language: Haskell2010 + +test-suite csg-tests + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: + tests + ghc-options: -Wall -Wcompat -O2 + build-depends: + base <5 + , bytestring + , csg + , simple-vec3 >=0.4 + , tasty + , tasty-hunit + , tasty-quickcheck + other-modules: + Paths_csg + default-language: Haskell2010 + +benchmark csg-benchmark + type: exitcode-stdio-1.0 + main-is: benchmark/Benchmark.hs + ghc-options: -Wall -Wcompat -O2 + build-depends: + base <5 + , criterion + , csg + , simple-vec3 >=0.4 + , strict + , vector + other-modules: + Paths_csg + default-language: Haskell2010 diff --git a/examples/cube.geo b/examples/cube.geo new file mode 100644 --- /dev/null +++ b/examples/cube.geo @@ -0,0 +1,17 @@ +solid box = orthobrick (-150, -150, -150; 150, 150, 150); + +solid rounded = sphere (0, 0, 0; 200); + +solid roundedbox = rounded and box; + +solid cylinder1 = cylinder (-160, 0, 0; 160, 0, 0; 100); +solid cylinder2 = cylinder (0, -160, 0; 0, 160, 0; 100); +solid cylinder3 = cylinder (0, 0, -160; 0, 0, 160; 100); + +solid cross = cylinder1 or cylinder2 or cylinder3; + +solid cutout = not cross; + +solid top = roundedbox and cutout; + +tlo top; diff --git a/examples/reentry.geo b/examples/reentry.geo new file mode 100644 --- /dev/null +++ b/examples/reentry.geo @@ -0,0 +1,13 @@ +# https://en.wikipedia.org/wiki/Atmospheric_entry#Blunt_body_entry_vehicles + +solid body = cone (0, 0, 0; 100; 0, 0, 130; 50); +solid head = cylinder (0, 0, 130; 0, 0, 145; 46); + +solid rounding = sphere (0, 0, 149.6629; 180); +solid cutoff = plane (0, 0, 0; 0, 0, 10); + +solid butt = rounding and cutoff; + +solid all = body or butt or head; + +tlo all; diff --git a/exe/raycaster.hs b/exe/raycaster.hs new file mode 100644 --- /dev/null +++ b/exe/raycaster.hs @@ -0,0 +1,230 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE OverloadedStrings #-} + +{-| + +Standalone raycaster for CSG objects powered by Gloss. + +|-} + +import Prelude hiding (FilePath) + +import Control.Applicative +import qualified Data.Strict.Maybe as S + +import GHC.Float +import Data.String +import Data.Monoid +import Data.Version + +import Graphics.Gloss.Data.Color +import Graphics.Gloss.Data.Display +import qualified Graphics.Gloss.Data.Point as G +import qualified Graphics.Gloss.Interface.Pure.Game as G +import Graphics.Gloss.Raster.Field hiding (Point) + +import Test.QuickCheck hiding ((><)) +import Turtle.Options + +import Filesystem.Path.CurrentOS + +import Data.CSG +import Data.CSG.Parser + +import Paths_csg + +data InteractionMode = None | Rotate | Pan + + +-- | World state with observation point parameters and event helpers. +-- Roll is always 0. +data World = World { dist :: Double + -- ^ Distance to origin. + , pitch :: Double + , yaw :: Double + -- ^ Yaw of camera as if it was in origin. + , target :: Point + -- ^ Where camera looks at. + , holdPoint :: Maybe (Float, Float) + -- ^ Point where mouse button was held down. + , mode :: InteractionMode + } + + +-- | Command line options for caster. +data Options = Options + { geoFile :: Maybe FilePath + , width :: Int + , height :: Int + , pixels :: Int + , brightRGBA :: (Float, Float, Float, Float) + -- ^ Color for bright surfaces parallel to view plane (RGB) + , darkRGBA :: (Float, Float, Float, Float) + -- ^ Color for dark surfaces perpendicular to view plane (RGB) + } + + +optParser :: Parser Options +optParser = Options + <$> optional (argPath "geo-file" "Geometry definition file") + <*> (optInt "width" 'w' "Window width" <|> pure 500) + <*> (optInt "height" 'h' "Window height" <|> pure 500) + <*> (optInt "pixels" 'p' "Pixels per ray (set to 1 for top quality)" + <|> pure 1) + <*> (optRead "bright-color" 'b' + (fromString $ + "Color for bright surfaces (parallel to the view plane) " <> + "as (R, G, B, A) with each component between 0.0 and 1.0") + <|> pure (0.9, 0.9, 0.9, 1)) + <*> (optRead "dark-color" 'd' "Color for dark surfaces" + <|> pure (0, 0, 0, 1)) + + +-- | The factor between distance and view port width/height. +scaleFactor :: Double +scaleFactor = 225.0 +{-# INLINE scaleFactor #-} + + +initialDistance :: Double +initialDistance = 225 + + +-- | Initial world. +start :: World +start = World initialDistance 0 0 origin Nothing None + + +uncurry4 :: (a -> b -> c -> d -> e) -> (a, b, c, d) -> e +uncurry4 f (a, b, c, d) = f a b c d + + +-- | Scale deltas between hold and release coordinates by this number. +dragFactor :: Double +dragFactor = pi / 180 + + +-- | Change distance by this amount per one mouse wheel step. +zoomFactor :: Double +zoomFactor = 0.1 + + +-- | Handle mouse drag to change pitch & yaw and mouse wheel to zoom. +handleEvents :: G.Event -> World -> World +handleEvents e w = + case e of + G.EventKey (G.MouseButton G.LeftButton) G.Down _ c -> + w{holdPoint = Just c, mode = Rotate} + G.EventKey (G.MouseButton G.RightButton) G.Down _ c -> + w{holdPoint = Just c, mode = Pan} + G.EventKey (G.MouseButton _) G.Up _ _ -> + w{holdPoint = Nothing} + G.EventKey (G.MouseButton G.WheelDown) _ _ _ -> + w{dist = dist w + zoomFactor} + G.EventKey (G.MouseButton G.WheelUp) _ _ _ -> + w{dist = dist w - zoomFactor} + G.EventKey (G.Char 'r') G.Down _ _ -> + w{ target = origin + , yaw = 0 + , pitch = 0 + , dist = initialDistance + } + G.EventMotion p@(x, y) -> + case holdPoint w of + Nothing -> w + Just (u, v) -> + let + xdelta = float2Double (x - u) * dragFactor + ydelta = float2Double (y - v) * dragFactor + in + case mode w of + Rotate -> w{ holdPoint = Just p + , yaw = yaw w - xdelta + , pitch = pitch w + ydelta + } + Pan -> + let + !(_, sX, sY) = buildCartesian (yaw w) (pitch w) + in + w{ holdPoint = Just p + , target = target w <+> (sX .^ xdelta) <-> (sY .^ ydelta) + } + _ -> w + _ -> w + + +-- | Build cartesian axes from yaw and pitch with 0 roll. Angles are +-- in radians. +buildCartesian :: Double -> Double -> (Vec3, Vec3, Vec3) +buildCartesian yaw pitch = (u, v, w) + where u = fromXYZ (cos yaw * cos pitch, sin yaw * cos pitch, sin pitch) + v = fromXYZ (- (sin yaw), cos yaw, 0) + w = u >< v +{-# INLINE buildCartesian #-} + + +programName :: String +programName = "csg-raycaster " ++ showVersion version + + +casterField :: Int + -- ^ Window width. + -> Int + -- ^ Window height. + -> Int + -- ^ Pixels per point. + -> Solid + -- ^ Solid to show. + -> Color + -- ^ Bright color. + -> Color + -- ^ Dark color. + -> IO () +casterField width height pixels solid bright' dark' = + let + display = InWindow programName (width, height) (0, 0) + makePixel :: World -> G.Point -> Color + !wS = fromIntegral (width `div` 2) + !hS = fromIntegral (height `div` 2) + makePixel !w (x, y) = + let + !d = dist w + !wScale = -(wS * d / scaleFactor) + !hScale = (hS * d / scaleFactor) + !(n, sX, sY) = buildCartesian (yaw w) (pitch w) + !p = n .^ (-d) <+> target w + ray :: Ray + !ray = Ray (p + <+> (sX .^ (float2Double x * wScale)) + <+> (sY .^ (float2Double y * hScale)), n) + in + case ray `cast` solid of + S.Just (HitPoint _ (S.Just hn)) -> + mixColors factor (1 - factor) bright' dark' + where + factor = abs $ double2Float $ invert n .* hn + _ -> white + {-# INLINE makePixel #-} + in + playField display (pixels, pixels) 5 start makePixel + handleEvents + (flip const) + + +-- | Read solid def and program arguments, run the actual caster on +-- success. +main :: IO () +main = do + Options{..} <- options (fromString programName) optParser + solid <- + case geoFile of + Nothing -> Right <$> generate arbitrary + Just fp -> parseGeometryFile (encodeString fp) + case solid of + Right b -> do + putStrLn $ "Rendering " <> show b + casterField width height pixels b + (uncurry4 makeColor brightRGBA) + (uncurry4 makeColor darkRGBA) + Left e -> error $ "Problem when reading solid definition: " ++ e diff --git a/src/Data/CSG.hs b/src/Data/CSG.hs new file mode 100644 --- /dev/null +++ b/src/Data/CSG.hs @@ -0,0 +1,650 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} + +{-| + +Types and routines for constructive solid geometry. + +This module provides constructors for complex solids as well as +membership predicates and routines to compute intersections of such +solids with a ray. + +-} + +module Data.CSG + ( -- * Examples + -- $examples + -- * Solids + Solid + -- ** Primitives + , plane + , sphere + , cylinder + , cone + -- ** Complex solids + , cuboid + , coneFrustum + , cylinderFrustum + + -- ** Operations + , intersect + , unite + , complement + , subtract + + -- * Ray casting + , Point + , Vec3 + , Ray(..) + , HitPoint(..) + , HitSegment + , Trace + , trace + , cast + + -- * Membership + , inside + + , module V3 + ) + +where + +import Prelude hiding (Just, Nothing, Maybe, subtract) + +import Data.Strict.Maybe +import Data.Strict.Tuple +import Test.QuickCheck (Arbitrary(..), frequency, sized) + +import Data.Vec3 hiding (Vec3, Matrix) +import qualified Data.Vec3 as V3 + + +-- $examples +-- +-- "Data.CSG" uses 'Vec3' to represent vectors and points: +-- +-- >>> let p1 = fromXYZ (5, -6.5, -5) +-- >>> toXYZ (origin :: Point) +-- (0.0,0.0,0.0) +-- +-- Define some solids: +-- +-- >>> let s = sphere origin 5.0 +-- >>> let b = cuboid (fromXYZ (-1, -1, -1)) (fromXYZ (1, 1, 1)) +-- +-- See "Data.CSG.Parser" for a non-programmatic way to define solids. +-- +-- Test if a point is inside the solid: +-- +-- >>> origin `inside` (s `intersect` b) +-- True +-- +-- >>> origin `inside` (s `subtract` b) +-- False +-- +-- Find the distance to the next intersection of a ray with a solid, along with the +-- surface normal: +-- +-- >>> let axis = fromXYZ (1, 2, 10) +-- >>> let solid = cylinder origin axis 2.0 `intersect` sphere origin 3.5 +-- >>> let ray = Ray (p1, origin <-> p1) +-- >>> ray `cast` solid +-- Just (HitPoint 0.7422558525331708 (Just (CVec3 0.7155468474912454 (-0.6952955216188516) 6.750441957464598e-2))) +-- +-- Load a solid definition from a file: +-- +-- >>> import Data.CSG.Parser +-- >>> Right solid2 <- parseGeometryFile "examples/reentry.geo" +-- >>> ray `cast` solid2 +-- Just (HitPoint 10.877824491509912 (Just (CVec3 (-0.5690708596937849) 0.7397921176019203 0.3589790793088691))) + +#ifdef WITH_TRIPLES +type Vec3 = TVec3 +#else +-- | We use 'CVec3' as a simple replacement for @(Double, Double, +-- Double)@. 'CVec3' implements a contiguous storage scheme for +-- Unboxed and Storable vectors which shows better performance. +-- Compile this package with @triples@ flag and run benchmarks to see +-- the difference. +type Vec3 = CVec3 +#endif +type Point = Vec3 +type Matrix = V3.Matrix Vec3 + + +-- | A ray described by the equation @p(t) = p_0 + v * t@ with an +-- initial point @p_0@ and a direction @v@. Substituting a specific +-- time @t'@ in the equation yields a position of a point @p(t')@ on +-- the ray. For negative values of @t'@, position precedes the initial +-- point. +newtype Ray = Ray (Point, Vec3) + + +-- | A point at which a ray intersects a surface, given as a distance +-- from the ray's initial point and an outward normal to the surface +-- at the hit point. If hit is in infinity, then normal is 'Nothing'. +-- If the hit occures on the same line but precedes the initial point +-- of the ray, the distance is negative. +-- +-- Note that this datatype is strict only on first argument: we do not +-- compare normals when combining traces and thus do not force +-- calculation of normals. +data HitPoint = HitPoint !Double (Maybe Vec3) + deriving (Eq, Show) + + +instance Ord HitPoint where + compare (HitPoint t1 _) (HitPoint t2 _) = compare t1 t2 + + +-- | A segment of ray inside a solid. +type HitSegment = (Pair HitPoint HitPoint) + + +-- | Trace of a ray/line on a solid is a list of segments +-- corresponding to the portions of the ray inside the solid. +-- +-- > O - ray +-- > \ +-- > \ +-- > +------------ +-- > ---/ * \--- +-- > -/ * \- +-- > / * \ +-- > ( * - trace ) +-- > solid -\ * / +-- > -\ * /- +-- > ---\ * /--- +-- > --------+---- +-- > \ +-- > \ +-- > _\/ +-- > \ +-- +-- Each 'HitSegment' is defined by a pair of 'HitPoint's on the ray +-- line. +-- +-- Ends of segments or intervals are calculated by intersecting the +-- ray and the surface of the primitive. This is done with 'trace', +-- which substitutes the equation of ray @p(t) = p_o + v * t@ into the +-- equation which defines the surface and solves it for @t@. +-- +-- Hit points may in lie in infinity. For example, because a ray +-- intersects a plane only once, a half-space primitive defined by +-- this plane results in a half-interval trace of a ray: +-- +-- > / +-- > / +-- > / +-- > O========================+*****************> +-- > | / | +-- > ray / goes to infinity +-- > / +-- > / +-- > / +-- > / - surface of half-space +-- +-- If the solid is a composition, traces from primitives are then +-- combined according to operations used to define the solid (union, +-- intersection or complement). +-- +-- Although only convex primitives are used in the current +-- implementation, operations on solids may result in concave solids, +-- which is why trace is defined as a list of segments. +-- +-- In this example, solid is a sphere with a cutout: +-- +-- > ------------- +-- > ----/ \---- +-- > -/ \- +-- > -/ \- +-- > -/ ----------- \- +-- > / --/ \-- \ +-- > / -/ \- \ +-- > / / \ \ +-- > / / \ \ +-- > | hs1 | | hs2 | +-- > - - -+*******+- - - - - - O============+*******+=========> +-- > | | | | | +-- > \ \ ray / / +-- > \ \ / / +-- > \ -\ /- / +-- > \ --\ /-- / +-- > -\ ----------- /- +-- > -\ /- +-- > -\ /- +-- > ----\ /---- +-- > ------------- +-- +-- Here, the full trace contains two segments: @hs1@ and @hs2@. +-- Moreover, 'trace' treats ray as a line with a designated point on +-- it, in reference to which distances to hit points are calculated. +-- This means that @hs1@ will have negative distances from the initial +-- point as that segment precedes it. +type Trace = [HitSegment] + + +-- | IEEE positive infinity. +infinityP :: Double +infinityP = (/) 1 0 + + +-- | Negative infinity. +infinityN :: Double +infinityN = -infinityP + + +-- | Hit in negative infinity. +hitN :: HitPoint +hitN = HitPoint infinityN Nothing + + +-- | Hit in positive infinity. +hitP :: HitPoint +hitP = HitPoint infinityP Nothing + + +-- | CSG solid is a recursive composition of primitive objects or other +-- solids. +data Solid = Plane !Vec3 !Double + -- ^ Half-space defined by a unit outward normal and a + -- distance of boundary plane from the origin. + | Sphere !Vec3 !Double + -- ^ Sphere defined by a center point and a radius. + | Cylinder !Vec3 !Point !Double + -- ^ Infinite circular cylinder defined by a normalized axis + -- vector, a point on axis and a radius. + | Cone !Vec3 !Point !Double !Matrix !Double !Double + -- ^ Cone defined by an inward axis direction, a vertex and + -- a cosine to the angle h between the axis and the + -- generatrix. + -- + -- Additionally, a transformation matrix $n * n - cos^2 h$, + -- tangent of angle and odelta are stored for intersection + -- calculations. + | Union !Solid !Solid + | Intersection !Solid !Solid + | Complement !Solid + deriving (Eq, Show) + + +instance Arbitrary Solid where + -- There's got to be a nicer way to write this. + -- + -- We can't use generic-random here because Solid fields must be + -- populated by smart constructors. Generating them independently + -- will break internal assumptions (such as normals being unit + -- vectors or the way matrix field for Cone is populated). + arbitrary = sized $ \n -> + frequency $ + [ (4, sphere <$> arbitrary <*> arbitrary) + , (4, cuboid <$> arbitrary <*> arbitrary) + , (3, cylinder <$> arbitrary <*> arbitrary <*> arbitrary) + , (3, cone <$> arbitrary <*> arbitrary <*> arbitrary) + , (3, cylinderFrustum <$> arbitrary <*> arbitrary <*> arbitrary) + , (3, coneFrustum <$> arbitrary <*> arbitrary) + ] ++ + -- Recurse + if n == 0 then [] else + [ (8, unite <$> arbitrary <*> arbitrary) + , (3, intersect <$> arbitrary <*> arbitrary) + , (3, complement <$> arbitrary) + ] + +-- | A half-space defined by an arbitary point on the boundary plane +-- and an outward normal (not necessarily a unit vector). +plane :: Point -> Vec3 -> Solid +plane p n = Plane nn (p .* nn) + where + nn = normalize n + + +-- | A sphere defined by a center point and a radius. +sphere :: Vec3 -> Double -> Solid +sphere = Sphere + + +-- | A rectangular cuboid with faces parallel to axes, defined by two +-- opposite vertices. +cuboid :: Point -> Point -> Solid +cuboid p1 p2 = + plane p1' (fromXYZ (1, 0, 0)) + `intersect` + plane p1' (fromXYZ (0, 1, 0)) + `intersect` + plane p1' (fromXYZ (0, 0, 1)) + `intersect` + plane p2' (fromXYZ (-1, 0, 0)) + `intersect` + plane p2' (fromXYZ (0, -1, 0)) + `intersect` + plane p2' (fromXYZ (0, 0, -1)) + where + (x1, y1, z1) = toXYZ p1 + (x2, y2, z2) = toXYZ p2 + p2' = fromXYZ (min x1 x2, min y1 y2, min z1 z2) + p1' = fromXYZ (max x1 x2, max y1 y2, max z1 z2) + + +-- | An infinite circular cylinder defined by two arbitary points on +-- axis and a radius. +cylinder :: Point -> Point -> Double -> Solid +cylinder p1 p2 = Cylinder (normalize $ p2 <-> p1) p1 + + +-- | A finite right circular cylinder defined by two points on its top +-- and bottom and a radius. +cylinderFrustum :: Point -> Point -> Double -> Solid +cylinderFrustum pb pt r = + plane pt axis + `intersect` + plane pb (invert axis) + `intersect` + cylinder pb pt r + where + axis = pt <-> pb + + +-- | An infinite right circular cone defined by an outward axis +-- vector, an apex point and an angle between the generatrix and the +-- axis (in degrees, less than 90). +cone :: Vec3 -> Point -> Double -> Solid +cone a o h = + let + rads = h * pi / 180 + h' = cos rads + n = normalize $ invert a + gamma = diag (-h' * h') + m = addM (n `vxv` n) gamma + ta = tan rads + odelta = n .* o + in + Cone n o h' m ta odelta + + +-- | A conical frustum defined by two points on its axis with radii at +-- that points. One of radii may be zero (in which case one of frustum +-- ends will be the apex). +coneFrustum :: (Point, Double) -> (Point, Double) -> Solid +coneFrustum (p1, r1) (p2, r2) = + let + -- Direction from pb to pt is towards apex. Corresponding + -- radii are rb > rt. + (pb, rb, pt, rt) = if r1 < r2 + then (p2, r2, p1, r1) + else (p1, r1, p2, r2) + -- Cone axis and frustum height + gap = pt <-> pb + height = norm gap + axis = normalize gap + -- Calculate distance from pt to apex. + dist = if rt == 0 + then 0 + else height / (rb / rt - 1) + apex = pt <+> (axis .^ dist) + -- Angle between generatrix and axis + degs = atan (rb / (dist + norm (pt <-> pb))) * (180 / pi) + in + plane pt axis + `intersect` + plane pb (invert axis) + `intersect` + cone axis apex degs + + +-- | Intersection of two solids. +intersect :: Solid -> Solid -> Solid +intersect !b1 !b2 = Intersection b1 b2 + + +-- | Union of two solids. +unite :: Solid -> Solid -> Solid +unite !b1 !b2 = Union b1 b2 + + +-- | Complement to a solid (normals flipped). +complement :: Solid -> Solid +complement !b = Complement b + + +-- | Subtract a solid from another. +subtract :: Solid -> Solid -> Solid +subtract !b1 !b2 = intersect b1 $ complement b2 + + +-- | Trace of a ray on a solid. +trace :: Solid -> Ray -> Trace +{-# INLINE trace #-} + +trace b@(Plane n d) (Ray (pos, v)) = + let + !f = -(n .* v) + in + if f == 0 + then + -- If ray is parallel to plane and is inside, then trace is + -- the whole timeline. + [hitN :!: hitP | inside pos b] + else + let + !t = (pos .* n - d) / f + in + if f > 0 + then [HitPoint t (Just n) :!: hitP] + else [hitN :!: HitPoint t (Just n)] + +trace (Sphere c r) (Ray (pos, v)) = + let + !d = pos <-> c + !roots = solveq (v .* v) (v .* d * 2) (d .* d - r * r) + normal !u = normalize (u <-> c) + in + case roots of + Nothing -> [] + Just (t1 :!: t2) -> + [HitPoint t1 (Just $ normal $ moveBy pos v t1) :!: + HitPoint t2 (Just $ normal $ moveBy pos v t2)] + +trace (Cylinder n c r) (Ray (pos, v)) = + let + d = (pos <-> c) >< n + e = v >< n + roots = solveq (e .* e) (d .* e * 2) (d .* d - r * r) + normal u = normalize $ h <-> (n .^ (h .* n)) + where h = u <-> c + in + case roots of + Nothing -> [] + Just (t1 :!: t2) -> + [HitPoint t1 (Just $ normal $ moveBy pos v t1) :!: + HitPoint t2 (Just $ normal $ moveBy pos v t2)] + +trace (Cone n c _ m ta odelta) (Ray (pos, v)) = + let + delta = pos <-> c + c2 = dotM v v m + c1 = dotM v delta m + c0 = dotM delta delta m + roots = solveq c2 (2 * c1) c0 + normal !u = normalize $ nx .^ (1 / ta) <-> ny .^ ta + where h = u <-> c + -- Component of h parallel to cone axis + ny' = n .^ (n .* h) + ny = normalize ny' + -- Perpendicular component + nx = normalize $ h <-> ny' + in + case roots of + Nothing -> [] + Just (t1 :!: t2) -> + let + pos1 = moveBy pos v t1 + pos2 = moveBy pos v t2 + in + case ((pos1 .* n - odelta) > 0, (pos2 .* n - odelta) > 0) of + (True, True) -> [HitPoint t1 (Just $ normal pos1) :!: + HitPoint t2 (Just $ normal pos2)] + (True, False) -> [hitN :!: + HitPoint t1 (Just $ normal pos1)] + (False, True) -> [HitPoint t2 (Just $ normal pos2) :!: + hitP] + (False, False) -> [] + +trace (Intersection b1 b2) !p = + intersectTraces tr1 tr2 + where + tr1 = trace b1 p + tr2 = trace b2 p + +trace (Union b1 b2) !p = + uniteTraces tr1 tr2 + where + tr1 = trace b1 p + tr2 = trace b2 p + +trace (Complement b) !p = + complementTrace $ trace b p + + +-- | Find the next point where a ray hits a solid, if any. +-- +-- Here we consider only future intersections: the 'HitPoint' is +-- guaranteed to have non-negative distance (unlike when using +-- 'trace'). +-- +-- This means that if the ray starts inside the solid the only way to +-- tell that from 'cast' result is to compare it's direction and the +-- surface normal at the hit point. +cast :: Ray -> Solid -> Maybe HitPoint +cast r b = + case intersectTraces onlyFutureHits (trace b r) of + (:!:) hp@(HitPoint _ (Just _)) _ : _ -> Just hp + (:!:) (HitPoint _ Nothing) hp@(HitPoint _ (Just _)) : _ -> Just hp + _ -> Nothing + where + onlyFutureHits = [HitPoint 0 Nothing :!: HitPoint infinityP Nothing] + + +-- | Union of two traces. +uniteTraces :: Trace -> Trace -> Trace +uniteTraces u [] = u +uniteTraces u (v:t2) = + uniteTraces (unite1 u v) t2 + where + merge :: HitSegment -> HitSegment -> HitSegment + merge (a1 :!: b1) (a2 :!: b2) = min a1 a2 :!: max b1 b2 + {-# INLINE merge #-} + unite1 :: Trace -> HitSegment -> Trace + unite1 [] hs = [hs] + unite1 t@(hs1@(a1 :!: b1):tr') hs2@(a2 :!: b2) + | b1 < a2 = hs1:unite1 tr' hs2 + | a1 > b2 = hs2:t + | otherwise = unite1 tr' (merge hs1 hs2) + {-# INLINE unite1 #-} +{-# INLINE uniteTraces #-} + + +-- | Intersection of two traces. +intersectTraces :: Trace -> Trace -> Trace +intersectTraces tr1 tr2 = + let + -- Overlap two overlapping segments + overlap :: HitSegment -> HitSegment -> HitSegment + overlap (a1 :!: b1) (a2 :!: b2) = max a1 a2 :!: min b1 b2 + {-# INLINE overlap #-} + in + case tr2 of + [] -> [] + (hs2@(a2 :!: b2):tr2') -> + case tr1 of + [] -> [] + (hs1@(a1 :!: b1):tr1') | b1 < a2 -> intersectTraces tr1' tr2 + | b2 < a1 -> intersectTraces tr1 tr2' + | otherwise -> overlap hs1 hs2:intersectTraces tr1' tr2 +{-# INLINE intersectTraces #-} + + +-- | Complement to a trace (normals flipped). +complementTrace :: Trace -> Trace +complementTrace ((sp@(HitPoint ts _) :!: ep):xs) = + start ++ complementTrace' ep xs + where + flipNormals :: HitSegment -> HitSegment + flipNormals (HitPoint t1 n1 :!: HitPoint t2 n2) = + HitPoint t1 (invert <$> n1) :!: HitPoint t2 (invert <$> n2) + {-# INLINE flipNormals #-} + -- Start from infinity if first hitpoint is finite + start = if isInfinite ts + then [] + else [flipNormals $ hitN :!: sp] + complementTrace' :: HitPoint -> Trace -> Trace + complementTrace' c ((a :!: b):tr) = + -- Bridge between the last point of the previous segment and + -- the first point of the next one. + flipNormals (c :!: a):complementTrace' b tr + complementTrace' a@(HitPoint t _) [] = + -- End in infinity if last hitpoint is finite + [flipNormals (a :!: hitP) | not (isInfinite t)] +complementTrace [] = [hitN :!: hitP] +{-# INLINE complementTrace #-} + + +-- | True if the point is in inside the solid. +inside :: Point -> Solid -> Bool +{-# INLINE inside #-} + +inside !pos (Plane n d) = (pos .* n - d) < 0 + +inside !pos (Sphere c r) = norm (pos <-> c) < r + +inside !pos (Cylinder n c r) = + norm (h <-> (n .^ (h .* n))) < r + where + h = pos <-> c + +inside !pos (Cone n c a _ _ _) = + n .* normalize (pos <-> c) > a + +inside !p (Intersection b1 b2) = inside p b1 && inside p b2 + +inside !p (Union b1 b2) = inside p b1 || inside p b2 + +inside !p (Complement b) = not $ inside p b + + +-- | Move point by velocity vector for given time and return new +-- position. +moveBy :: Point + -- ^ Current position. + -> Vec3 + -- ^ Velocity. + -> Double + -- ^ Time step. + -> Point +moveBy !p !v !t = p <+> (v .^ t) +{-# INLINE moveBy #-} + + +-- | Solve quadratic equation @ax^2 + bx + c = 0@. +-- +-- If less than two roots exist, Nothing is returned. +solveq :: Double + -- ^ a + -> Double + -- ^ b + -> Double + -- ^ c + -> Maybe (Pair Double Double) +solveq !a !b !c + | d > 0 = Just $ min r1 r2 :!: max r1 r2 + | otherwise = Nothing + where + d = b * b - 4 * a * c + q = sqrt d + t = 2 * a + r = - b / t + s = q / t + r1 = r - s + r2 = r + s +{-# INLINE solveq #-} diff --git a/src/Data/CSG/Parser.hs b/src/Data/CSG/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Data/CSG/Parser.hs @@ -0,0 +1,316 @@ +{-# LANGUAGE OverloadedStrings #-} + +{-| + +Parser for CSG solid definition format. The format uses text files and +is inspired by NETGEN 4.x @.geo@ format. + +Each definition may contain several solid definitions and ends with +the top level object declaration. Right hand side of solid +equations may reference other solids to allow composing of complex +solids. + +> # comment +> +> # define several primitives +> solid b1 = sphere (0, 0, 0; 5); +> solid p1 = plane (0, 0, 0; 1, 0, 0); +> +> # define a composition +> solid comp = b1 and p1; +> +> # make it the top level object +> tlo comp; + +Statements must end with a semicolon (newlines are optional). +Whitespace is ignored. + +Multiple-solid compositions are __right-associative__, so @b1 and b2 +or b3@ really means @b1 and (b2 or b3)@. Keep simpler objects on the +left and when in doubt stick to combining two solids at a time. + +Top-level object line must reference a previously defined solid. + +Syntax for primitives follows the signatures of 'CSG' constructors +for 'CSG.plane' and 'CSG.sphere', but differs for cylinder and +cone, as this module provides access only to frustums +('CSG.cylinderFrustum' and 'CSG.coneFrustum'). + +[Half-space] @plane (px, py, pz; nx, ny, nz)@, where @(px, py, pz)@ +is a point on a plane which defines the half-space and @(nx, ny, +nz)@ is a normal to the plane (outward to the half-space), not +necessarily a unit vector. + +[Brick] @orthobrick (ax, ay, az; bx, by, bz)@, where @(ax, ay, az)@ is +a vertex with minimum coordinates and @(bx, by, bz)@ is a vertex with +maximum coordinates. + +[Sphere] @sphere (cx, cy, cz; r)@, where @(cx, cy, cz)@ is a +central point of a sphere and @r@ is radius. + +[Right circular cylinder] @cylinder (p1x, p1y, p1z; p2x, p2y, p2z; +r)@ where @(p1x, p1y, p1z)@ and @(p2x, p2y, p2z)@ are bottom and +top points on axis and @r@ is radius. + +[Right circular conical frustum] @cone (p1x, p1y, p1z; r1; p2x, +p2y, p2z; r2)@ where @(p1x, p1y, p1z)@ and @(p2x, p2y, p2z)@ are +bottom and top points on cone axis and @r1@, @r2@ are the +corresponding radii. + +-} + +module Data.CSG.Parser + ( parseGeometry + , parseGeometryFile + ) + +where + +import Prelude as P + +import Control.Applicative +import qualified Control.Exception as E +import Control.Monad.Trans.Class +import Control.Monad.Trans.State.Strict + +import Data.Attoparsec.ByteString.Char8 +import Data.ByteString.Char8 as B + +import qualified Data.Map as M + +import Data.Vec3 hiding (Vec3, Matrix) + +import qualified Data.CSG as CSG + + +-- | Transformer which adds a lookup table to a monad. +type TableT a k v = StateT (M.Map k v) a + + +-- | Add an entry to the lookup table. +addEntry :: (Ord k, Monad a) => k -> v -> TableT a k v () +addEntry key value = fmap (M.insert key value) get >>= put + + +-- | Lookup entry in the table. +getEntry :: (Ord k, Monad a) => k -> TableT a k v (Maybe v) +getEntry key = fmap (M.lookup key) get + + +-- | Parser with a lookup table. +type CSGParser = TableT Parser String CSG.Solid + + +lp :: Parser Char +lp = char '(' + + +rp :: Parser Char +rp = char ')' + + +eq :: Parser Char +eq = char '=' + + +cancer :: Parser Char +cancer = char ';' + + +comma :: Parser Char +comma = char ',' + + +-- | Read three comma-separated doubles into point. +-- +-- > ::= ',' ',' +triple :: Parser CSG.Point +triple = fmap fromXYZ $ + (,,) <$> double + <*> + (skipSpace *> comma *> skipSpace *> + double + <* skipSpace <* comma <* skipSpace) + <*> + double + + +keywords :: [String] +keywords = [ "solid" + , "tlo" + , "orthobrick" + , "plane" + , "sphere" + , "cylinder" + , "cone" + ] + + +-- | Read variable name or fail if it's a keyword. +varName :: CSGParser String +varName = do + k <- lift $ many1 (letter_ascii <|> digit) + if k `P.notElem` keywords + then return k + else fail ("Unexpected keyword when reading a solid name: " ++ k) + + +-- | Look up a solid in the table by its name or fail if it's not +-- defined yet. +readName :: CSGParser CSG.Solid +readName = do + k <- varName + v <- getEntry k + case v of + Just b -> return b + _ -> fail $ "Undefined solid: " ++ k + + +-- > ::= +-- > 'plane (' ';' ')' +plane :: Parser CSG.Solid +plane = CSG.plane <$> + (string "plane" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> triple <* skipSpace <* rp) + + +-- > ::= +-- > 'orthobrick (' ';' ')' +orthobrick :: Parser CSG.Solid +orthobrick = CSG.cuboid <$> + (string "orthobrick" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> triple <* skipSpace <* rp) + + +-- > ::= +-- > 'sphere (' ';' ')' +sphere :: Parser CSG.Solid +sphere = CSG.sphere <$> + (string "sphere" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp) + + +-- > ::= +-- > 'cylinder (' ';' ';' ')' +cylinder :: Parser CSG.Solid +cylinder = CSG.cylinderFrustum <$> + (string "cylinder" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp) + + +-- > ::= +-- > 'cone (' ';' ';' ';' ')' +cone :: Parser CSG.Solid +cone = CSG.coneFrustum <$> + ((,) <$> + (string "cone" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> double)) <*> + ((,) <$> + (skipSpace *> cancer *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp)) + + +primitive :: Parser CSG.Solid +primitive = plane <|> orthobrick <|> sphere <|> cylinder <|> cone + + +-- > ::= 'not' +complement :: CSGParser CSG.Solid +complement = CSG.complement <$> (lift (string "not" *> skipSpace) *> solid) + + +-- > ::= 'or' +union :: CSGParser CSG.Solid +union = binary "or" CSG.unite + + +-- > ::= 'and' +intersection :: CSGParser CSG.Solid +intersection = binary "and" CSG.intersect + + +-- | Parse binary operation on two bodies with given composition +-- operators. +-- +-- Note that due to the way 'binary' and 'solid' combinators recurse +-- into each other multi-solid composition chains are +-- __right-associative__. However, this also means that if we keep +-- simpler solids on the left then ray casting routines will have a +-- chance to work faster and terminate earlier. +binary :: ByteString -> (CSG.Solid -> CSG.Solid -> CSG.Solid) -> CSGParser CSG.Solid +binary op compose = do + b1 <- uncomposedSolid + lift (skipSpace *> string op *> skipSpace) + b2 <- solid + return $ compose b1 b2 + + +-- | Read a stamement which adds a new solid entry to the lookup +-- table. +-- +-- > ::= +-- > 'solid' '=' ';' +statement :: CSGParser () +statement = do + lift $ skipSpace *> string "solid" *> skipSpace + k <- varName + lift $ skipSpace <* eq <* skipSpace + v <- solid <* lift (cancer *> skipSpace) + addEntry k v + + +-- | Expression is either a primitive, a reference to previously +-- defined solid or an operation on expressions. +-- +-- > ::= | | | | +solid :: CSGParser CSG.Solid +solid = union <|> intersection <|> complement <|> uncomposedSolid + + +-- | Used to terminate left branch of binary compositions. +-- +-- > ::= | +uncomposedSolid :: CSGParser CSG.Solid +uncomposedSolid = lift primitive <|> readName + + +-- | Top-level object declaration. +-- +-- > ::= 'tlo' ';' +topLevel :: CSGParser CSG.Solid +topLevel = lift (string "tlo" *> skipSpace) *> + readName + <* lift (cancer <* skipSpace) + + +-- | Read one-line comment starting with hash sign. +comment :: Parser () +comment = char '#' >> manyTill anyChar endOfLine >> return () + + +-- | Read sequence of statements which define solids, and finally read +-- top level object definition. +-- +-- > ::= | | +geoFile :: CSGParser CSG.Solid +geoFile = many1 (lift comment <|> statement) *> topLevel + + +-- | Read solid definition. If parsing fails, return error message as a +-- string. +parseGeometry :: ByteString -> Either String CSG.Solid +parseGeometry input = + case parseOnly (runStateT geoFile M.empty) input of + Right (b, _) -> Right b + Left msg -> Left msg + + +-- | Read solid definition from a file. If parsing fails, return error +-- message as a string. +parseGeometryFile :: FilePath -> IO (Either String CSG.Solid) +parseGeometryFile file = do + res <- E.try $ B.readFile file + return $ case res of + Right d -> parseGeometry d + Left e -> Left $ show (e :: E.IOException) diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -0,0 +1,34 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE ScopedTypeVariables #-} + +import Data.ByteString.Char8 as B +import Test.Tasty +import Test.Tasty.QuickCheck +import Test.Tasty.HUnit + +import Data.CSG +import Data.CSG.Parser + +tests :: [TestTree] +tests = + [ testCase "Parsing cube.geo" $ do + f <- B.readFile "examples/cube.geo" + let box = cuboid (fromXYZ (-150, -150, -150)) (fromXYZ (150, 150, 150)) + rounded = sphere origin 200 `intersect` box + cyl1 = cylinderFrustum + (fromXYZ (-160, 0, 0)) (fromXYZ (160, 0, 0)) 100 + cyl2 = cylinderFrustum + (fromXYZ (0, -160, 0)) (fromXYZ (0, 160, 0)) 100 + cyl3 = cylinderFrustum + (fromXYZ (0, 0, -160)) (fromXYZ (0, 0, 160)) 100 + cross = cyl1 `unite` (cyl2 `unite` cyl3) + cutout = complement cross + top = rounded `intersect` cutout + parseGeometry f @=? Right top + , testProperty "CSG complement membership" + (\(b :: Solid) (p :: Point) -> + p `inside` b == not (p `inside` complement b)) + ] + +main :: IO () +main = defaultMain $ testGroup "Tests" tests diff --git a/tests/doctest-driver.hs b/tests/doctest-driver.hs new file mode 100644 --- /dev/null +++ b/tests/doctest-driver.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF doctest-discover #-}