csg (empty) → 0.1
raw patch · 13 files changed
+1678/−0 lines, 13 filesdep +QuickCheckdep +attoparsecdep +basesetup-changed
Dependencies added: QuickCheck, attoparsec, base, bytestring, containers, criterion, csg, doctest, doctest-discover, gloss, gloss-raster, simple-vec3, strict, system-filepath, tasty, tasty-hunit, tasty-quickcheck, transformers, turtle, vector
Files
- CHANGELOG.md +5/−0
- LICENSE +30/−0
- README.md +140/−0
- Setup.hs +2/−0
- benchmark/Benchmark.hs +117/−0
- csg.cabal +123/−0
- examples/cube.geo +17/−0
- examples/reentry.geo +13/−0
- exe/raycaster.hs +230/−0
- src/Data/CSG.hs +650/−0
- src/Data/CSG/Parser.hs +316/−0
- tests/Main.hs +34/−0
- tests/doctest-driver.hs +1/−0
+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Changelog++## [0.1] - 2018-03-20++[0.1]: https://github.com/dzhus/csg/tree/0.1
+ LICENSE view
@@ -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.
+ README.md view
@@ -0,0 +1,140 @@+# CSG: constructive solid geometry library++[](https://travis-ci.org/dzhus/csg)+[](https://hackage.haskell.org/package/csg)+[](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.++++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
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmark/Benchmark.hs view
@@ -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
+ csg.cabal view
@@ -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
+ examples/cube.geo view
@@ -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;
+ examples/reentry.geo view
@@ -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;
+ exe/raycaster.hs view
@@ -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
+ src/Data/CSG.hs view
@@ -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 #-}
+ src/Data/CSG/Parser.hs view
@@ -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> ::= <double> ',' <double> ',' <double>+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 (' <triple> ';' <triple> ')'+plane :: Parser CSG.Solid+plane = CSG.plane <$>+ (string "plane" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> triple <* skipSpace <* rp)+++-- > <orthobrick> ::=+-- > 'orthobrick (' <triple> ';' <triple> ')'+orthobrick :: Parser CSG.Solid+orthobrick = CSG.cuboid <$>+ (string "orthobrick" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> triple <* skipSpace <* rp)+++-- > <sphere> ::=+-- > 'sphere (' <triple> ';' <double> ')'+sphere :: Parser CSG.Solid+sphere = CSG.sphere <$>+ (string "sphere" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp)+++-- > <cylinder> ::=+-- > 'cylinder (' <triple> ';' <triple> ';' <double> ')'+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 (' <triple> ';' <double> ';' <triple> ';' <double> ')'+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+++-- > <complement> ::= 'not' <solid>+complement :: CSGParser CSG.Solid+complement = CSG.complement <$> (lift (string "not" *> skipSpace) *> solid)+++-- > <union> ::= <uncomposed-solid> 'or' <solid>+union :: CSGParser CSG.Solid+union = binary "or" CSG.unite+++-- > <intersection> ::= <uncomposed-solid> 'and' <solid>+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.+--+-- > <statement> ::=+-- > 'solid' <varname> '=' <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> ::= <union> | <intersection> | <complement> | <primitive> | <reference>+solid :: CSGParser CSG.Solid+solid = union <|> intersection <|> complement <|> uncomposedSolid+++-- | Used to terminate left branch of binary compositions.+--+-- > <uncomposed-solid> ::= <primitive> | <reference>+uncomposedSolid :: CSGParser CSG.Solid+uncomposedSolid = lift primitive <|> readName+++-- | Top-level object declaration.+--+-- > <tlo> ::= 'tlo' <solid> ';'+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> ::= <statement> <geoFile> | <comment> <geoFile> | <tlo>+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)
+ tests/Main.hs view
@@ -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
+ tests/doctest-driver.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF doctest-discover #-}