From f860138ad6c6a1a4e8557cff6540db083be89746 Mon Sep 17 00:00:00 2001 From: Frank Staals <991345+noinia@users.noreply.github.com> Date: Sun, 25 Aug 2024 21:02:51 +0200 Subject: [PATCH] Fixed bug related to points in halfspaces (#244) * raytracing stuff + discovered a bug in the hyperplane normal business * fixed one of the testcases since what I watned was too strong * separating onSideTest and verticalSideTest * doctest is now happy * fixing tests :) * adding a progress bar * copied over the camera module from HGeometry 0.X * some extra type annotations; since apparently GHC 9.8 was confused about this * intersecting a ball with a line * halfline ball intersection as well :) --- hgeometry-examples/hgeometry-examples.cabal | 14 ++ hgeometry-examples/raytracer/Main.hs | 102 +++++++++++ hgeometry/data/test-with-ipe/golden/ball.ipe | 156 +++++++++++++++++ hgeometry/hgeometry.cabal | 5 + .../src/HGeometry/Ball/CenterAndRadius.hs | 155 ++++++++++++++--- hgeometry/kernel/src/HGeometry/HalfLine.hs | 21 +++ .../kernel/src/HGeometry/HyperPlane/Class.hs | 63 ++++--- .../src/HGeometry/HyperPlane/NonVertical.hs | 8 +- .../kernel/src/HGeometry/Line/General.hs | 2 + .../src/HGeometry/Line/PointAndVector.hs | 25 ++- .../kernel/src/HGeometry/Triangle/Class.hs | 16 +- hgeometry/kernel/test/HGeometry/BallSpec.hs | 28 ++- .../kernel/test/HGeometry/HyperPlaneSpec.hs | 140 ++++++++------- .../kernel/test/HGeometry/TriangleSpec.hs | 10 +- hgeometry/point/src/HGeometry/Point.hs | 1 + hgeometry/point/src/HGeometry/Point/Class.hs | 12 ++ .../src/HGeometry/ConvexHull/R3/Naive.hs | 2 - hgeometry/src/HGeometry/Graphics/Camera.hs | 161 ++++++++++++++++++ .../LowerEnvelope/Connected/BruteForce.hs | 2 +- .../HGeometry/Plane/LowerEnvelope/Naive.hs | 2 +- hgeometry/test-with-ipe/test/BallSpec.hs | 67 ++++++++ .../test/VoronoiDiagram/VoronoiSpec.hs | 34 ---- hgeometry/test/Graphics/CameraSpec.hs | 64 +++++++ hgeometry/test/LowerEnvelope/RegionsSpec.hs | 2 +- todo.org | 6 +- 25 files changed, 924 insertions(+), 174 deletions(-) create mode 100644 hgeometry-examples/raytracer/Main.hs create mode 100644 hgeometry/data/test-with-ipe/golden/ball.ipe create mode 100644 hgeometry/src/HGeometry/Graphics/Camera.hs create mode 100644 hgeometry/test-with-ipe/test/BallSpec.hs create mode 100644 hgeometry/test/Graphics/CameraSpec.hs diff --git a/hgeometry-examples/hgeometry-examples.cabal b/hgeometry-examples/hgeometry-examples.cabal index aa525b481..361f0e6b8 100644 --- a/hgeometry-examples/hgeometry-examples.cabal +++ b/hgeometry-examples/hgeometry-examples.cabal @@ -215,6 +215,20 @@ executable hgeometry-sampler -- Paths_hgeometry_examples -- Miso.Event.Extra +-------------------------------------------------------------------------------- +-- * Raytracer example + +executable hgeometry-raytracer + import: setup, miso-setup + build-depends: + JuicyPixels >= 3.3.9 && < 4 + , terminal-progress-bar >= 0.4.2 && < 0.5 + hs-source-dirs: raytracer + main-is: Main.hs + other-modules: + Paths_hgeometry_examples + -- Miso.Event.Extra + -------------------------------------------------------------------------------- diff --git a/hgeometry-examples/raytracer/Main.hs b/hgeometry-examples/raytracer/Main.hs new file mode 100644 index 000000000..c72f3af1d --- /dev/null +++ b/hgeometry-examples/raytracer/Main.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE QuasiQuotes #-} +module Main(main) where + +import Codec.Picture +import Control.Lens +import Data.Foldable as F +import Data.Maybe +import Data.Monoid +import Data.Word +import HGeometry.Ball +import HGeometry.Ext +import HGeometry.HalfLine +import HGeometry.Intersection +import HGeometry.Point +import HGeometry.Vector +import Paths_hgeometry_examples +import qualified System.File.OsPath as File +import System.OsPath +import System.ProgressBar + +-------------------------------------------------------------------------------- + +type Color = PixelRGBA8 + +type Scene = [Ball (Point 3 R) :+ Color] + +type R = Double + +type Ray = HalfLine (Point 3 R) + +rayThrough :: Int -> Int -> Ray +rayThrough x y = let q = Point3 x y 0 &coordinates %~ fromIntegral + in HalfLine cameraPos (q .-. cameraPos) + +-- | Intersect to try and get the color +ballColor :: Ray -> Ball (Point 3 R) :+ Color -> Maybe Color +ballColor r (b :+ c) + | r `intersects` b = Just c + | otherwise = Nothing + +shootRay :: Ray -> Scene -> Color +shootRay r = fromMaybe backgroundColor . getFirst . foldMap (First . ballColor r) + +renderPixel :: Vector 2 Int -> Scene -> Int -> Int -> PixelRGBA8 +renderPixel (Vector2 w h) scene x y = shootRay ray scene + where + ray = rayThrough x y + + -- clamped :: Int -> Int -> Word8 + -- clamped x m = fromIntegral $ (255 * x) `div` m + +renderWithProgress :: IO () -> Vector 2 Int -> Scene -> IO (Image PixelRGBA8) +renderWithProgress reportProgress dims@(Vector2 w h) scene = + withImage w h $ \x y -> let !pix = renderPixel dims scene x y + in pix <$ reportProgress + +-------------------------------------------------------------------------------- +-- * Settings + +backgroundColor :: PixelRGBA8 +backgroundColor = PixelRGBA8 maxBound maxBound maxBound 0 + +theScene :: Scene +theScene = [ Ball (Point3 128 128 128) 50 :+ PixelRGBA8 200 0 0 255 + ] + + +cameraPos = Point3 0 0 10000 + + +aspectRatio :: Rational +aspectRatio = 16 / 9 + +outputWidth :: Int +outputWidth = 400 + +outputDimensions :: Vector 2 Int +outputDimensions = Vector2 outputWidth (ceiling $ fromIntegral outputWidth / aspectRatio) + +refreshRate :: Double +refreshRate = 10 + + +viewportDims :: Vector 2 Double +viewportDims = let Vector2 w h = fromIntegral <$> outputDimensions + desiredHeight = 2 + in Vector2 desiredHeight (desiredHeight * (h/w)) + +-------------------------------------------------------------------------------- + +amountOfWork (Vector2 w h) = w * h + +main :: IO () +main = do + let initialProgress = Progress 0 (amountOfWork outputDimensions) () + progressBar <- newProgressBar defStyle refreshRate initialProgress + + imageData <- renderWithProgress (incProgress progressBar 1) + outputDimensions theScene + + let bs = encodePng imageData + File.writeFile [osp|foo.png|] bs diff --git a/hgeometry/data/test-with-ipe/golden/ball.ipe b/hgeometry/data/test-with-ipe/golden/ball.ipe new file mode 100644 index 000000000..97323a475 --- /dev/null +++ b/hgeometry/data/test-with-ipe/golden/ball.ipe @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +0.6 0 0 0.6 0 0 e 0.4 0 0 0.4 0 0 e + + + +0.6 0 0 0.6 0 0 e + + + +0.5 0 0 0.5 0 0 e + +0.6 0 0 0.6 0 0 e 0.4 0 0 0.4 0 0 e + + + +-0.6 -0.6 m 0.6 -0.6 l 0.6 0.6 l -0.6 0.6 l h +-0.4 -0.4 m 0.4 -0.4 l 0.4 0.4 l -0.4 0.4 l h + + +-0.6 -0.6 m 0.6 -0.6 l 0.6 0.6 l -0.6 0.6 l h + + +-0.5 -0.5 m 0.5 -0.5 l 0.5 0.5 l -0.5 0.5 l h + +-0.6 -0.6 m 0.6 -0.6 l 0.6 0.6 l -0.6 0.6 l h +-0.4 -0.4 m 0.4 -0.4 l 0.4 0.4 l -0.4 0.4 l h + + +-0.43 -0.57 m 0.57 0.43 l 0.43 0.57 l -0.57 -0.43 l h + +-0.43 0.57 m 0.57 -0.43 l 0.43 -0.57 l -0.57 0.43 l h + + + +0 0 m -1.0 0.333 l -1.0 -0.333 l h + + +0 0 m -1.0 0.333 l -1.0 -0.333 l h + + +0 0 m -1.0 0.333 l -0.8 0 l -1.0 -0.333 l h + + +0 0 m -1.0 0.333 l -0.8 0 l -1.0 -0.333 l h + + +0 0 m -1.0 0.333 l -1.0 -0.333 l h + + +0 0 m -1.0 0.333 l -0.8 0 l -1.0 -0.333 l h + + +0 0 m -1.0 0.333 l -0.8 0 l -1.0 -0.333 l h + + +-1.0 0.333 m 0 0 l -1.0 -0.333 l + + +0 0 m -1.0 0.333 l -1.0 -0.333 l h +-1 0 m -2.0 0.333 l -2.0 -0.333 l h + + + +0 0 m -1.0 0.333 l -1.0 -0.333 l h +-1 0 m -2.0 0.333 l -2.0 -0.333 l h + + + + + + +-505.0 -1000.0 m +495.0 1000.0 l +-1000.0 20.0 m +1000.0 20.0 l +-1000.0 8.0 m +1000.0 8.0 l +4.0 0.0 m +504.0 -1000.0 l +-5.0 10.0 m +197.0 -1000.0 l +10.0 0.0 0.0 10.0 0.0 0.0 e8.0 0.0 0.0 8.0 20.0 28.0 e16.0 0.0 0.0 16.0 200.0 28.0 e-8.0 -6.0 m +0.0 10.0 l +186.143593539449 20.0 m +213.856406460551 20.0 l +-6.0 8.0 m +6.0 8.0 l +4.0 0.0 m +7.3761226035642204 -6.752245207128441 l +-4.758998912463262 8.794994562316312 m +-1.0102318567675068 -9.948840716162465 l + \ No newline at end of file diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal index 78ba2aa13..baea84bd6 100644 --- a/hgeometry/hgeometry.cabal +++ b/hgeometry/hgeometry.cabal @@ -340,6 +340,8 @@ library exposed-modules: HGeometry + HGeometry.Graphics.Camera + -- HGeometry.DelaunayTriangulation HGeometry.PlaneGraph @@ -715,6 +717,7 @@ test-suite hspec RangeTreeSpec VerticalRayShootingSpec LineSegment.BentheyOttmanNeighbourSpec + Graphics.CameraSpec hs-source-dirs: test build-depends: @@ -751,6 +754,8 @@ test-suite with-ipe-hspec HalfPlane.CommonIntersectionSpec Line.LowerEnvelopeSpec Paths_hgeometry + BallSpec + autogen-modules: Paths_hgeometry hs-source-dirs: test-with-ipe/test diff --git a/hgeometry/kernel/src/HGeometry/Ball/CenterAndRadius.hs b/hgeometry/kernel/src/HGeometry/Ball/CenterAndRadius.hs index 1fb1bff3e..224b2c881 100644 --- a/hgeometry/kernel/src/HGeometry/Ball/CenterAndRadius.hs +++ b/hgeometry/kernel/src/HGeometry/Ball/CenterAndRadius.hs @@ -19,18 +19,24 @@ module HGeometry.Ball.CenterAndRadius , _BallSphere , _DiskCircle - ) where -import Control.Lens -import HGeometry.Ball.Class -import HGeometry.Intersection -import HGeometry.LineSegment --- import qualified HGeometry.Number.Radical as Radical -import HGeometry.Point -import HGeometry.Properties (NumType, Dimension) -import HGeometry.Vector + , IntersectionOf(..) + ) where +import Control.Lens +import HGeometry.Ball.Class +import HGeometry.HalfLine +import HGeometry.HalfSpace +import HGeometry.HyperPlane +import HGeometry.Intersection +import HGeometry.Line.PointAndVector +import HGeometry.LineSegment +import HGeometry.Number.Radical +import HGeometry.Point +import HGeometry.Properties (NumType, Dimension) +import HGeometry.Vector +import Prelude hiding (sqrt) -------------------------------------------------------------------------------- @@ -51,6 +57,7 @@ instance Point_ point (Dimension point) (NumType point) => Ball_ (Ball point) po fromCenterAndSquaredRadius = Ball -------------------------------------------------------------------------------- +-- * Point in ball type instance Intersection (Point d r) (Ball point) = Maybe (Point d r) @@ -67,30 +74,134 @@ instance ( Point_ point d r intersect q b | q `intersects` b = Just q | otherwise = Nothing --- type instance Intersection (Line point) (Ball point) = --- Maybe (IntersectionOf (LineSegment point) (Ball point)) - -data instance IntersectionOf (ClosedLineSegment point) (Ball point) = - Line_x_Ball_Point point - | Line_x_Ball_Segment (ClosedLineSegment point) - -deriving instance (Show point, Show (ClosedLineSegment point)) - => Show (IntersectionOf (ClosedLineSegment point) (Ball point)) -deriving instance (Eq point, Eq (ClosedLineSegment point)) - => Eq (IntersectionOf (ClosedLineSegment point) (Ball point)) +-------------------------------------------------------------------------------- +-- * Testing if a line/ ray /linesegment intersects a ball -type instance Intersection (ClosedLineSegment point) (Ball point) = - Maybe (IntersectionOf (ClosedLineSegment point) (Ball point)) +-- essentially this is all just computing the squared euclidean distance +-- between the object and the center, and testing if it is at most r +instance ( Point_ point d r + , Ord r, Fractional r + , Has_ Metric_ d r + ) => (LinePV d r) `HasIntersectionWith` (Ball point) where + intersects l (Ball c r) = squaredEuclideanDistTo c l <= r instance ( Point_ point d r, Point_ point' d r , Ord r, Fractional r , Has_ Metric_ d r , HasSquaredEuclideanDistance point' + , MkHyperPlaneConstraints d r ) => (ClosedLineSegment point') `HasIntersectionWith` (Ball point) where intersects s (Ball c r) = squaredEuclideanDistTo c s <= r +instance ( Point_ point d r, Point_ point' d r + , Ord r, Fractional r + , Has_ Metric_ d r + , HasSquaredEuclideanDistance point' + , MkHyperPlaneConstraints d r + ) => (HalfLine point') `HasIntersectionWith` (Ball point) where + intersects hl (Ball c r) = squaredEuclideanDistTo c hl <= r + +---------------------------------------- + +type instance Intersection (LinePV d r) (Ball point) = + Maybe (IntersectionOf (LinePV d r) (Ball point)) + +data instance IntersectionOf (LinePV d r) (Ball point) = + Line_x_Ball_Point (Point d r) + | Line_x_Ball_Segment (ClosedLineSegment (Point d r)) + +deriving instance (Show r, Has_ Additive_ d r) => Show (IntersectionOf (LinePV d r) (Ball point)) +deriving instance (Eq r, Eq (Vector d r)) => Eq (IntersectionOf (LinePV d r) (Ball point)) + + +instance ( Point_ point d r + , Ord r, Fractional r, Radical r + , Has_ Metric_ d r + ) => (LinePV d r) `IsIntersectableWith` (Ball point) where + intersect (LinePV p v) (Ball c' r) = case discr `compare` 0 of + LT -> Nothing + EQ -> Just $ Line_x_Ball_Point q0 -- line touches in q + GT -> Just $ Line_x_Ball_Segment (ClosedLineSegment q1 q2) + where + -- main idea: let q = p + \lambda v be an intersection point, we also have + -- squaredEuclideanDist q c' == squaredRadius (=r) this yields some quadratic + -- equation in \lambda, which we just solve using the ABC formula. In particular, we have + -- + -- (sum_i=1^d (p+\lambda v_i - c'_i)^2 = r) + -- (sum_i=1^d (\lambda v_i + (p_i- c'_i))^2 - r = 0 ) + -- (sum_i=1^d ((\lambda v_i)^2 + 2\lambda v_i(p_i- c'_i) + (p_i- c'_i)^2) - r = 0 ) + -- (sum_i=1^d (v_i^2\lambda^2 + 2v_i(p_i- c'_i)\lambda + (p_i- c'_i)^2) - r = 0 ) + -- (lambda^2 sum_i=1^d v_i^2\ + \lambda sum_i=1^d 2v_i(p_i- c'_i) + sum_i=1^d (p_i- c'_i)^2) - r = 0 ) + + + a = v `dot` v -- sum_i v_i^2 + b = 2 * (v `dot` u) -- sum_i v_i(p_i-c_i) + c = (u `dot` u) - r -- sum_i (p_i-c_i)^2 - radius^2 + + u = p .-. (c'^.asPoint) -- helper + + discr = b^2 - 4*a*c + discr' = sqrt discr + da = 2*a + + lambda1 = ((negate discr') - b) / da -- the two solutions + lambda2 = (discr' - b) / da -- + -- note: v must have non-zero length; and thus a (and therefore da) are non-zero. + + q1 = p .+^ (lambda1 *^ v) + q2 = p .+^ (lambda2 *^ v) + + -- if the discr is zer0 there is only one solution: + lambda0 = (negate b) / da + q0 = p .+^ (lambda0 *^ v) + +---------------------------------------- + +type instance Intersection (HalfLine point') (Ball point) = + Maybe (IntersectionOf (LinePV (Dimension point) (NumType point)) (Ball point)) + +instance ( Point_ point d r, Point_ point' d r + , Ord r, Fractional r, Radical r + , Has_ Metric_ d r + , MkHyperPlaneConstraints d r + , HasSquaredEuclideanDistance point' + ) => (HalfLine point') `IsIntersectableWith` (Ball point) where + intersect (HalfLine p v) b = intersect (LinePV p' v) b >>= \case + Line_x_Ball_Point q + | q `intersects` h -> Just $ Line_x_Ball_Point q + | otherwise -> Nothing + Line_x_Ball_Segment seg@(ClosedLineSegment a c) -> + case (a `intersects` h, c `intersects` h) of + (False,False) -> Nothing + (False,True) -> Just $ Line_x_Ball_Segment (ClosedLineSegment p' c) + (True,False) -> Just $ Line_x_Ball_Segment (ClosedLineSegment a p') + (True,True) -> Just $ Line_x_Ball_Segment seg + where + h :: HalfSpace d r + h = HalfSpace Positive (fromPointAndNormal p' v) + p' = p^.asPoint + +type instance Intersection (ClosedLineSegment point') (Ball point) = + Maybe (IntersectionOf (LinePV (Dimension point) (NumType point)) (Ball point)) + +-- data instance IntersectionOf (ClosedLineSegment point) (Ball point) = +-- LineSegment_x_Ball_Point point +-- | LineSegment_x_Ball_Segment (ClosedLineSegment point) + +-- deriving instance (Show point, Show (ClosedLineSegment point)) +-- => Show (IntersectionOf (ClosedLineSegment point) (Ball point)) +-- deriving instance (Eq point, Eq (ClosedLineSegment point)) +-- => Eq (IntersectionOf (ClosedLineSegment point) (Ball point)) + + + + + + + + -------------------------------------------------------------------------------- diff --git a/hgeometry/kernel/src/HGeometry/HalfLine.hs b/hgeometry/kernel/src/HGeometry/HalfLine.hs index 1df07b93b..66a188e12 100644 --- a/hgeometry/kernel/src/HGeometry/HalfLine.hs +++ b/hgeometry/kernel/src/HGeometry/HalfLine.hs @@ -18,6 +18,8 @@ module HGeometry.HalfLine import Control.Lens import GHC.Generics (Generic) import GHC.TypeLits +import HGeometry.HalfSpace +import HGeometry.HyperPlane import HGeometry.Intersection import HGeometry.Interval.Class import HGeometry.Line.Intersection @@ -114,6 +116,25 @@ instance ( Ord r, Fractional r, Point_ point 2 r m = supportingLine hl -- the left side is suposedly the halfplane containing the halfLine +-------------------------------------------------------------------------------- + +instance ( Point_ point d r + , Has_ Metric_ d r + , Ord r, Fractional r + , MkHyperPlaneConstraints d r + ) => HasSquaredEuclideanDistance (HalfLine point) where + pointClosestTo q hl@(HalfLine p v) + | r `intersects` h = r + | otherwise = p' + where + p' = p^.asPoint + r = pointClosestTo q (supportingLine hl) + h = HalfSpace Positive (fromPointAndNormal p' v) :: HalfSpace d r + -- main idea: compute the point closest to the supporting line of the halfline. if this + -- point lies in the halfspace h defined by the ray (e.g. for which the ray is the + -- normal), then we've actually found the point closest to the ray. Otherwise the origin + -- of the ray is simply the closest point. + -------------------------------------------------------------------------------- diff --git a/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs b/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs index a05662ed0..a8b1593f8 100644 --- a/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs +++ b/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs @@ -49,8 +49,9 @@ import Prelude hiding (head,last) -- >>> let myLineAsNV = NonVerticalHyperPlane (Vector2 1 2) :: NonVerticalHyperPlane 2 Double -- >>> let myVerticalLine = HyperPlane (Vector3 5 (-1) 0) :: HyperPlane 2 Double -- --- >>> let myOtherLine = HyperPlane2 4 3 2 :: HyperPlane 2 Double --- >>> let myPlane = HyperPlane3 10 2 3 (-1) :: HyperPlane 3 Double +-- >>> let myOtherLine = HyperPlane2 4 3 2 :: HyperPlane 2 Double +-- >>> let myOtherNVLine = NonVerticalHyperPlane (Vector2 (-3/2 :: Double) (-2)) +-- >>> let myPlane = HyperPlane3 10 2 3 (-1) :: HyperPlane 3 Double -- -- -- 'myLine', 'myLineAgain', and 'myLineAsNV' all represent the line y = 1*x + 2 @@ -129,8 +130,7 @@ class ( NumType hyperPlane ~ r normalVector :: (Num r, Eq r, 1 <= d) => hyperPlane -> Vector d r default normalVector :: (KnownNat d, Num r, Eq r, 1 <= d) => hyperPlane -> Vector d r - normalVector h = let a = suffix $ hyperPlaneEquation h - in if signum (a^.last) == 1 then a else negated a + normalVector h = negated . suffix $ hyperPlaneEquation h {-# INLINE normalVector #-} -- https://en.wikipedia.org/wiki/Normal_(geometry)#Hypersurfaces_in_n-dimensional_space -- states that: if the hyperplane is defined as the solution set of a single linear equation @@ -159,14 +159,8 @@ class ( NumType hyperPlane ~ r q `onHyperPlane` h = (== 0) $ evalHyperPlaneEquation h q {-# INLINE onHyperPlane #-} - -- | Test if a point lies on a hyperplane. For non-vertical hyperplanes, returns whether - -- the point is *above* the hyperplane or not with respect to the d^th - -- dimension. (I.e. if (and only if) q lies above the hyperplane h, then q `onSideTest` - -- h return GT.) - -- - -- For vertical hyperplanes (with respect to dimension d), we return 'LT' when the point - -- is on the "left". - -- + -- | Test if a point lies on a hyperplane. Returns the sign when evaluating the + -- hyperplane equation. -- -- >>> Point2 0 2 `onSideTest` myLineAgain -- EQ @@ -187,23 +181,18 @@ class ( NumType hyperPlane ~ r -- EQ -- -- >>> Point2 0 1 `onSideTest` myOtherLine - -- GT + -- LT -- >>> Point2 0 (-2) `onSideTest` myOtherLine -- EQ -- >>> Point2 1 (-3.5) `onSideTest` myOtherLine -- EQ -- >>> Point2 1 (-4) `onSideTest` myOtherLine - -- LT + -- GT onSideTest :: (Point_ point d r, Ord r, Num r) => point -> hyperPlane -> Ordering - onSideTest q h - | (hyperPlaneEquation h)^.last <= 0 = 0 `compare` evalHyperPlaneEquation h q - | otherwise = evalHyperPlaneEquation h q `compare` 0 + onSideTest q h = 0 `compare` evalHyperPlaneEquation h q {-# INLINE onSideTest #-} - -- if ad is positive, it seems the sign of the comparison should - -- flip i.e. so that we should have evalHyperPlaneEquation h q - -- `compare` 0 instead. This still seems very werid to me. --- | Produce a point that lies on the hyperplane. No gurantees are given abot which point +-- | Produce a point that lies on the hyperplane. No gurantees are given about which point -- -- >>> pointOn myLine -- Point2 (-2.0) 0.0 @@ -324,6 +313,36 @@ class HyperPlane_ hyperPlane d r => NonVerticalHyperPlane_ hyperPlane d r where -- Just (Vector2 (-1.5) (-2.0)) hyperPlaneCoefficients :: Lens' hyperPlane (Vector d r) + -- | Test if a point q lies above a non-vertical hyperplane h (i.e. verticalSideTest q h + -- == GT), on the hyperplane (== EQ), or below (LT). + -- + -- >>> Point2 0 2 `verticalSideTest` myLineAsNV + -- EQ + -- >>> Point2 1 3 `verticalSideTest` myLineAsNV + -- EQ + -- >>> Point2 1 5 `verticalSideTest` myLineAsNV + -- GT + -- >>> Point2 4 5 `verticalSideTest` myLineAsNV + -- LT + -- + -- >>> Point2 0 1 `verticalSideTest` myOtherNVLine + -- GT + -- >>> Point2 0 (-2) `verticalSideTest` myOtherNVLine + -- EQ + -- >>> Point2 1 (-3.5) `verticalSideTest` myOtherNVLine + -- EQ + -- >>> Point2 1 (-4) `verticalSideTest` myOtherNVLine + -- LT + verticalSideTest :: (1 <= d, Point_ point d r, Ord r, Num r) + => point -> hyperPlane -> Ordering + default verticalSideTest :: ( Point_ point d r + , Num r, Ord r + , 1 <= d + , 1 + (d-1) ~ d -- silly silly agian :( + , Has_ Metric_ d r + , Has_ Additive_ (d-1) r + ) => point -> hyperPlane -> Ordering + verticalSideTest q h = (q^.dCoord) `compare` evalAt (projectPoint q :: Point (d-1) r) h -------------------------------------------------------------------------------- -- * Functions on Hyperplanes @@ -374,3 +393,5 @@ instance (NonVerticalHyperPlane_ hyperPlane d r) {-# INLINE evalAt #-} hyperPlaneCoefficients = core.hyperPlaneCoefficients {-# INLINE hyperPlaneCoefficients #-} + verticalSideTest q = verticalSideTest q . view core + {-# INLINE verticalSideTest #-} diff --git a/hgeometry/kernel/src/HGeometry/HyperPlane/NonVertical.hs b/hgeometry/kernel/src/HGeometry/HyperPlane/NonVertical.hs index 4dd795f10..bc5a70989 100644 --- a/hgeometry/kernel/src/HGeometry/HyperPlane/NonVertical.hs +++ b/hgeometry/kernel/src/HGeometry/HyperPlane/NonVertical.hs @@ -113,12 +113,14 @@ asNonVerticalHyperPlane' e -instance ( MkHyperPlaneConstraints d r +instance ( MkHyperPlaneConstraints d r, Has_ Additive_ (d-1) r , 2 <= d ) => HyperPlane_ (NonVerticalHyperPlane d r) d r where + -- normalVector h = let a = suffix $ hyperPlaneEquation h + -- in if signum (a^.last) == 1 then a else negated a -instance ( MkHyperPlaneConstraints d r +instance ( MkHyperPlaneConstraints d r, Has_ Additive_ (d-1) r , Fractional r, Eq r , 2 <= d ) => ConstructableHyperPlane_ (NonVerticalHyperPlane d r) d r where @@ -134,7 +136,7 @@ instance ( MkHyperPlaneConstraints d r {-# INLINE hyperPlaneFromEquation #-} -instance ( MkHyperPlaneConstraints d r, 1 + (d-1) ~ d +instance ( MkHyperPlaneConstraints d r, 1 + (d-1) ~ d, Has_ Additive_ (d-1) r , Num r , 2 <= d ) => NonVerticalHyperPlane_ (NonVerticalHyperPlane d r) d r where diff --git a/hgeometry/kernel/src/HGeometry/Line/General.hs b/hgeometry/kernel/src/HGeometry/Line/General.hs index 4f2b5c107..fb3ab17ab 100644 --- a/hgeometry/kernel/src/HGeometry/Line/General.hs +++ b/hgeometry/kernel/src/HGeometry/Line/General.hs @@ -56,6 +56,8 @@ instance HyperPlane_ (VerticalOrLineEQ r) 2 r where VerticalLineThrough x -> Vector3 1 0 (-x) NonVertical l -> hyperPlaneEquation l onHyperPlane = onLine + + -- FIXME: remove this implementation; just use the default onSideTest q = \case VerticalLineThrough x -> (q^.xCoord) `compare` x NonVertical l -> onSideTest q l diff --git a/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs b/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs index b7673769d..1cea55e91 100644 --- a/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs +++ b/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs @@ -17,9 +17,13 @@ module HGeometry.Line.PointAndVector , HasSupportingLine(..) , fromLinearFunction , toLinearFunction + + , SideTestUpDown(..), OnSideUpDownTest(..) , liesAbove, liesBelow , SideTest(..), onSide + , leftHalfPlane, rightHalfPlane + , bisector , perpendicularTo, isPerpendicularTo @@ -34,12 +38,14 @@ import Data.Ord (comparing) import GHC.Generics (Generic) import GHC.TypeLits import HGeometry.Ext +import HGeometry.Sign import HGeometry.HyperPlane.Class import HGeometry.Intersection import HGeometry.Line.Class import HGeometry.Line.Intersection import HGeometry.Line.LineEQ import HGeometry.Point +import HGeometry.HalfSpace -- import HGeometry.Point.EuclideanDistance -- import HGeometry.Point.Orientation.Degenerate import HGeometry.Properties (NumType, Dimension) @@ -347,9 +353,8 @@ bisector p q = let v = q .-. p h = view asPoint $ p .+^ (v ^/ 2) in perpendicularTo (LinePV h v) --- | Given a line l with anchor point p and vector v, get the line --- perpendicular to l that also goes through p. The resulting line m is --- oriented such that v points into the left halfplane of m. +-- | Given a line l with anchor point p and vector v, get the line m perpendicular to l +-- that also goes through p. The line is oriented *into* the right halfplane of l. -- -- >>> perpendicularTo $ LinePV (Point2 3 4) (Vector2 (-1) 2) -- LinePV (Point2 3 4) (Vector2 2 1) @@ -393,6 +398,20 @@ cmpSlope :: forall r. (Num r, Ord r -------------------------------------------------------------------------------- +-- | Given the oriented line, computes the halfspace left of the line. +leftHalfPlane :: (Num r, Ord r) => LinePV 2 r -> HalfSpaceF (LinePV 2 r) +leftHalfPlane l = HalfSpace sign l + where + sign = let (LinePV p v) = perpendicularTo l + in case onSideTest (p .+^ v) l of + LT -> Positive + _ -> Negative + +-- | Given the oriented line, computes the halfspace right of the line. +rightHalfPlane :: (Num r, Ord r) => LinePV 2 r -> HalfSpaceF (LinePV 2 r) +rightHalfPlane l = let HalfSpace s _ = leftHalfPlane l + in HalfSpace (flipSign s) l + {- -- | Lines are transformable, via line segments instance ( Fractional r diff --git a/hgeometry/kernel/src/HGeometry/Triangle/Class.hs b/hgeometry/kernel/src/HGeometry/Triangle/Class.hs index d90d137ca..7fe53a3c6 100644 --- a/hgeometry/kernel/src/HGeometry/Triangle/Class.hs +++ b/hgeometry/kernel/src/HGeometry/Triangle/Class.hs @@ -17,10 +17,9 @@ module HGeometry.Triangle.Class ) where import Control.Lens -import Data.Ord import HGeometry.HalfSpace -import HGeometry.HyperPlane import HGeometry.Point +import HGeometry.Line.PointAndVector import HGeometry.Properties (NumType, Dimension) import HGeometry.Vector @@ -87,22 +86,19 @@ toCounterClockwiseTriangle t@(Triangle_ a b c) -- -- >>> let t = Triangle origin (Point2 0 (-1)) (Point2 (-1) 0) :: Triangle (Point 2 Int) -- >>> mapM_ print $ intersectingHalfPlanes t --- HalfSpace Negative (HyperPlane [0,0,1]) --- HalfSpace Positive (HyperPlane [-1,-1,-1]) --- HalfSpace Negative (HyperPlane [0,-1,0]) +-- HalfSpace Positive (LinePV (Point2 0 0) (Vector2 (-1) 0)) +-- HalfSpace Positive (LinePV (Point2 (-1) 0) (Vector2 1 (-1))) +-- HalfSpace Positive (LinePV (Point2 0 (-1)) (Vector2 0 1)) intersectingHalfPlanes :: ( Triangle_ triangle point , Point_ point 2 r , Num r, Ord r ) => triangle - -> Vector 3 (HalfSpace 2 r) + -> Vector 3 (HalfSpaceF (LinePV 2 r)) intersectingHalfPlanes (toCounterClockwiseTriangle -> Triangle_ u v w) = Vector3 (leftPlane u v) (leftPlane v w) (leftPlane w u) where - leftPlane p q = let dir = if (p^.xCoord, Down $ p^.yCoord) <= (q^.xCoord, Down $ q^.yCoord) - -- for verical planes, the left side counts as below - then Positive else Negative - in HalfSpace dir . hyperPlaneThrough $ Vector2 p q + leftPlane p q = leftHalfPlane $ LinePV (p^.asPoint) (q .-. p) -- | Given a point q and a triangle, q inside the triangle, get the baricentric diff --git a/hgeometry/kernel/test/HGeometry/BallSpec.hs b/hgeometry/kernel/test/HGeometry/BallSpec.hs index 75449d4a7..725977d38 100644 --- a/hgeometry/kernel/test/HGeometry/BallSpec.hs +++ b/hgeometry/kernel/test/HGeometry/BallSpec.hs @@ -1,6 +1,6 @@ module HGeometry.BallSpec where --- import Control.Lens +import Control.Lens import Control.Monad (forM_) -- import HGeometry.Ext import HGeometry.Number.Real.Rational @@ -8,7 +8,11 @@ import HGeometry.Ball import HGeometry.Intersection import HGeometry.LineSegment import HGeometry.Point --- import HGeometry.Vector +import HGeometry.HalfLine +import HGeometry.Vector +import HGeometry.HyperPlane +import HGeometry.HalfSpace +import HGeometry.Line.PointAndVector import Test.Hspec -- import Test.QuickCheck -- import Test.Util @@ -30,6 +34,26 @@ spec = do (mySeg `intersects` unitCircle @R) `shouldBe` True + it "closest point to ray" $ + let ray = HalfLine (Point3 (0 :: R) (1/2) 10000) (Vector3 0 0 (-1)) + in pointClosestTo (origin :: Point 3 R) ray `shouldBe` (Point3 0 (1/2) 0) + + it "closest point to line" $ + let ray = LinePV (Point3 (0 :: R) (1/2) 10000) (Vector3 0 0 (-1)) + in pointClosestTo (origin :: Point 3 R) ray `shouldBe` (Point3 0 (1/2) 0) + + it "proper halfspace" $ + let p = Point3 0 (1/2) 10000 :: Point 3 R + v = Vector3 0 0 (-1) + h = HalfSpace Positive (fromPointAndNormal p v) :: HalfSpace 3 R + in (h^.boundingHyperPlane.to normalVector) `shouldBe` v + + + it "ball intersects ray" $ + let ray = HalfLine (Point3 (0 :: R) (1/2) 10000) (Vector3 0 0 (-1)) + in (ray `intersects` (unitBall :: Ball (Point 3 R))) `shouldBe` True + + unitCircle :: (Num r) => Circle (Point 2 r) unitCircle = Circle origin 1 diff --git a/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs b/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs index e90515070..b38ec2316 100644 --- a/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs +++ b/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs @@ -3,20 +3,21 @@ {-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-} module HGeometry.HyperPlaneSpec where -import Control.Lens hiding (unsnoc) -import Data.Maybe (isJust) -import GHC.TypeNats -import HGeometry.HyperPlane -import HGeometry.HyperPlane.NonVertical -import HGeometry.Intersection -import HGeometry.Kernel.Instances () -import HGeometry.Line -import HGeometry.Number.Real.Rational -import HGeometry.Point -import HGeometry.Vector -import Test.Hspec -import Test.Hspec.QuickCheck -import Test.QuickCheck ((===), (==>)) +import Control.Lens hiding (unsnoc) +import Data.Maybe (isJust) +import GHC.TypeNats +import HGeometry.HalfSpace +import HGeometry.HyperPlane +import HGeometry.HyperPlane.NonVertical +import HGeometry.Intersection +import HGeometry.Kernel.Instances () +import HGeometry.Line +import HGeometry.Number.Real.Rational +import HGeometry.Point +import HGeometry.Vector +import Test.Hspec +import Test.Hspec.QuickCheck +import Test.QuickCheck ((===), (==>)) -------------------------------------------------------------------------------- @@ -88,34 +89,53 @@ spec = describe "HyperPlane Tests" $ do asNonVerticalHyperPlane (HyperPlane2 0 (-1) 0) `shouldBe` Nothing - prop "onside test non-vertical hyperplane2 (i.e. lines)" $ - \(h :: HyperPlane 2 R) (q :: Point 2 R) -> - isNonVertical h ==> - q `onSideTest` h === (case asNonVerticalHyperPlane h of - Just h' -> onSideTestNonVertical q h' - _ -> error "impossible" - ) - prop "onside test vertical lines" $ - \(x :: R) (q :: Point 2 R) -> - let h = HyperPlane2 x (-1) 0 -- vertical line through x - in q `onSideTest` h === (q^.xCoord) `compare` x - - - prop "onside test non-vertical hyperplane3 (i.e. planes)" $ - \(h :: HyperPlane 3 R) (q :: Point 3 R) -> - isNonVertical h ==> - q `onSideTest` h === (case asNonVerticalHyperPlane h of - Just h' -> onSideTestNonVertical q h' - _ -> error "impossible" - ) - -- prop "onside test vertical lines" $ - -- \(x :: R) (q :: Point 3 R) -> - -- let h = HyperPlane2 x (-1) 0 -- vertical line through x - -- in q `onSideTest` h === (q^.xCoord) `compare` x +-------------------------------------------------------------------------------- + prop "normal vector ok" $ + \(p :: Point 3 R) (n :: Vector 3 R) -> + (quadrance n > 0) ==> + normalVector (fromPointAndNormal p n :: HyperPlane 3 R) `isScalarMultipleOf` n + + prop "normal vector ok for non-vertical" $ + \(p :: Point 3 R) (n :: Vector 3 R) -> + (quadrance n > 0 && n^.zComponent /= 0) ==> + normalVector (fromPointAndNormal p n :: NonVerticalHyperPlane 3 R) `isScalarMultipleOf` n + + + + + + +-- Note that the normal vector we start with and the one we get back from 'normalVector' +--do not have to be identical: Take an arbitrary hyperplane, and consider its two normal +--vectors , and pick the negative one (i.e. pointing into the negative halfspace.; then +--build the plane using this negative normal. We then want the 'normalVector' function to +--return the positive normal vector instead. +-------------------------------------------------------------------------------- + prop "fromPointAnNormal and sideTest consistent for HyperPlane " $ + \(p :: Point 2 R) n -> + allOf components (>0) n ==> + ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: HyperPlane 2 R)) + `shouldBe` GT + prop "fromPointAnNormal and sideTest consistent for NonVerticalHyperPlane " $ + \(p :: Point 2 R) n -> + allOf components (>0) n ==> + ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: NonVerticalHyperPlane 2 R)) + `shouldBe` GT + prop "fromPointAnNormal and sideTest consistent for LineEQ" $ + \(p :: Point 2 R) n -> + allOf components (>0) n ==> + ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: LineEQ R)) + `shouldBe` GT + prop "normalVector and fromPointAndNormal consistent (HyperPlane 2 R)" $ + \(p :: Point 2 R) n -> + allOf components (>0) n ==> + normalVector (fromPointAndNormal p n :: HyperPlane 2 R) `shouldBe` n + +-------------------------------------------------------------------------------- prop "onside test NonVertical Hyperplanes i.e. lines) consitent " $ \(h :: NonVerticalHyperPlane 2 R) (q :: Point 2 R) -> - q `onSideTest` h === onSideTestNonVertical q h + q `verticalSideTest` h === onSideTestNonVertical q h it "onside test NonVertical Hyperplanes i.e. lines) manual " $ do (Point2 3 10) `onSideTestNonVertical` (LineEQ 1 1) `shouldBe` GT @@ -126,7 +146,6 @@ spec = describe "HyperPlane Tests" $ do (Point2 0 0) `onSideTestNonVertical` (LineEQ 2 1) `shouldBe` LT (Point2 0 4) `onSideTestNonVertical` (LineEQ 2 1) `shouldBe` GT - it "on side of vertical line / hyperplane" $ (Point2 0 1 `onSideTest` HyperPlane2 3 (-1) 0) `shouldBe` LT @@ -145,6 +164,7 @@ spec = describe "HyperPlane Tests" $ do it "on side of non-vertical line / hyperplane 5" $ (Point2 0 0 `onSideTest` HyperPlane2 (-1) (-1) (-1)) `shouldBe` GT + prop "pointOn produces a point on the hyperplane (2d)" $ \(h :: HyperPlane 2 R) -> onHyperPlane (pointOn h) h prop "pointOn produces a point on the hyperplane (3d)" $ @@ -154,42 +174,36 @@ spec = describe "HyperPlane Tests" $ do \(l :: LineEQ R) (m :: LineEQ R) -> (l `intersects` m) `shouldBe` (asHyp l `intersects` asHyp m) - - prop "fromPointAnNormal and sideTest consistent for HyperPlane " $ - \(p :: Point 2 R) n -> + prop "fromPointAnNormal and sideTest consistent for HyperPlane in R^3" $ + \(p :: Point 3 R) n -> allOf components (>0) n ==> - ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: HyperPlane 2 R)) + ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: HyperPlane 3 R)) `shouldBe` GT - prop "fromPointAnNormal and sideTest consistent for NonVerticalHyperPlane " $ - \(p :: Point 2 R) n -> - allOf components (>0) n ==> - ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: NonVerticalHyperPlane 2 R)) - `shouldBe` GT - prop "fromPointAnNormal and sideTest consistent for LineEQ" $ - \(p :: Point 2 R) n -> - allOf components (>0) n ==> - ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: LineEQ R)) - `shouldBe` GT - prop "normalVector and fromPointAndNormal consistent (HyperPlane 2 R)" $ - \(p :: Point 2 R) n -> + + prop "fromPointAnNormal and intersects halfSpace consistent in R^3" $ + \(p :: Point 3 R) n -> allOf components (>0) n ==> - normalVector (fromPointAndNormal p n :: HyperPlane 2 R) `shouldBe` n + ((p .+^ n) `intersects` (HalfSpace Positive (fromPointAndNormal p n) + :: HalfSpace 3 R + )) + `shouldBe` True + - prop "nonVertical sidetest means above (NonVHyperplane 2)" $ + prop "nonVertical verticalsidetest means above (NonVHyperplane 2)" $ \(q :: Point 2 R) (h :: NonVerticalHyperPlane 2 R) -> let y = evalAt (projectPoint q) h q' = q&yCoord .~ y+1 - in (q' `onSideTest` h) `shouldBe` GT - prop "nonVertical sidetest means above (NonVHyperplane 3)" $ + in (q' `verticalSideTest` h) `shouldBe` GT + prop "nonVertical verticalsidetest means above (NonVHyperplane 3)" $ \(q :: Point 3 R) (h :: NonVerticalHyperPlane 3 R) -> let z = evalAt (projectPoint q) h q' = q&zCoord .~ z+1 - in (q' `onSideTest` h) `shouldBe` GT - prop "nonVertical sidetest means above LineEQ " $ + in (q' `verticalSideTest` h) `shouldBe` GT + prop "nonVertical verticalsidetest means above LineEQ " $ \(q :: Point 2 R) (h :: LineEQ R) -> let y = evalAt (projectPoint q) h q' = q&yCoord .~ y+1 - in (q' `onSideTest` h) `shouldBe` GT + in (q' `verticalSideTest` h) `shouldBe` GT -- prop "intersect nonvertical conistent" $ diff --git a/hgeometry/kernel/test/HGeometry/TriangleSpec.hs b/hgeometry/kernel/test/HGeometry/TriangleSpec.hs index b7c1e0baf..04ecf50ef 100644 --- a/hgeometry/kernel/test/HGeometry/TriangleSpec.hs +++ b/hgeometry/kernel/test/HGeometry/TriangleSpec.hs @@ -64,14 +64,10 @@ spec = describe "intersection tests" $ do it "intersecting halfspaces" $ let t = Triangle origin (Point2 0 (-1)) (Point2 (-1) 0) :: Triangle (Point 2 Int) in (show <$> intersectingHalfPlanes t) `shouldBe` - (Vector3 "HalfSpace Negative (HyperPlane [0,0,1])" - "HalfSpace Positive (HyperPlane [-1,-1,-1])" - "HalfSpace Negative (HyperPlane [0,-1,0])" + (Vector3 "HalfSpace Positive (LinePV (Point2 0 0) (Vector2 (-1) 0))" + "HalfSpace Positive (LinePV (Point2 (-1) 0) (Vector2 1 (-1)))" + "HalfSpace Positive (LinePV (Point2 0 (-1)) (Vector2 0 1))" ) --- (Vector2 )) - --- undefined undefined undefined --- ) -------------------------------------------------------------------------------- diff --git a/hgeometry/point/src/HGeometry/Point.hs b/hgeometry/point/src/HGeometry/Point.hs index 4f769facd..456a6d43f 100644 --- a/hgeometry/point/src/HGeometry/Point.hs +++ b/hgeometry/point/src/HGeometry/Point.hs @@ -30,6 +30,7 @@ module HGeometry.Point , coord , xCoord, yCoord, zCoord, wCoord + , dCoord , projectPoint diff --git a/hgeometry/point/src/HGeometry/Point/Class.hs b/hgeometry/point/src/HGeometry/Point/Class.hs index 4f482236a..2754052e2 100644 --- a/hgeometry/point/src/HGeometry/Point/Class.hs +++ b/hgeometry/point/src/HGeometry/Point/Class.hs @@ -22,6 +22,7 @@ module HGeometry.Point.Class , pointFromList , coord , xCoord, yCoord, zCoord, wCoord + , dCoord -- , projectPoint -- , PointFor @@ -318,6 +319,17 @@ wCoord :: (4 <= d, Point_ point d r) => IndexedLens' Int point r wCoord = coord @4 {-# INLINABLE wCoord #-} +-- | Shorthand to access the last coordinate +-- +-- >>> (Point2 1 2 :: Point 2 Int) ^. dCoord +-- 2 +-- >>> (Point4 1 2 3 4 :: Point 4 Int) ^. dCoord +-- 4 +-- >>> (Point4 1 2 3 4 :: Point 4 Int) & dCoord %~ (+1) +-- Point4 1 2 3 5 +dCoord :: forall point d r. (1 <= d, Point_ point d r) => IndexedLens' Int point r +dCoord = coord @d + -------------------------------------------------------------------------------- -- | Data types that store points diff --git a/hgeometry/src/HGeometry/ConvexHull/R3/Naive.hs b/hgeometry/src/HGeometry/ConvexHull/R3/Naive.hs index d524213f1..9fa04b8f5 100644 --- a/hgeometry/src/HGeometry/ConvexHull/R3/Naive.hs +++ b/hgeometry/src/HGeometry/ConvexHull/R3/Naive.hs @@ -97,5 +97,3 @@ upperHalfSpaceOf :: (Ord r, Num r, Point_ point 3 r) upperHalfSpaceOf (Triangle p q r) = HalfSpace Positive h where h = hyperPlaneThrough $ Vector3 p q r - -- c = p&zCoord -~ 1 - -- s = if c `onSideTest` h /= GT then Positive else Negative diff --git a/hgeometry/src/HGeometry/Graphics/Camera.hs b/hgeometry/src/HGeometry/Graphics/Camera.hs new file mode 100644 index 000000000..553456603 --- /dev/null +++ b/hgeometry/src/HGeometry/Graphics/Camera.hs @@ -0,0 +1,161 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : HGeometry.Graphics.Camera +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- Description : Data type to represent a camera and some functions for working with it. +-- +-------------------------------------------------------------------------------- +module HGeometry.Graphics.Camera + ( Camera(Camera) + , cameraPosition, rawCameraNormal, rawViewUp + , viewPlaneDepth, nearDist, farDist, screenDimensions + + , cameraNormal, viewUp + + , cameraTransform, worldToView + + , toViewPort, perspectiveProjection, rotateCoordSystem + , flipAxes + + ) where + +import Control.Lens +import HGeometry.Matrix +import HGeometry.Point +import HGeometry.Transformation +import HGeometry.Vector +import HGeometry.Number.Radical + +-------------------------------------------------------------------------------- + +-- | A basic camera data type. The fields stored are: +-- +-- * the camera position, +-- * the raw camera normal, i.e. a unit vector into the center of the screen, +-- * the raw view up vector indicating which side points "upwards" in the scene, +-- * the viewplane depth (i.e. the distance from the camera position to the plane on which we project), +-- * the near distance (everything closer than this is clipped), +-- * the far distance (everything further away than this is clipped), and +-- * the screen dimensions. +-- +data Camera r = Camera { _cameraPosition :: !(Point 3 r) + , _rawCameraNormal :: !(Vector 3 r) + -- ^ unit vector from camera into center of the screen + , _rawViewUp :: !(Vector 3 r) + -- ^ viewUp; assumed to be unit vector + , _viewPlaneDepth :: !r + , _nearDist :: !r + , _farDist :: !r + , _screenDimensions :: !(Vector 2 r) + } deriving (Show,Eq,Ord) + +---------------------------------------- +-- * Field Accessor Lenses + +-- Lemmih: Writing out the lenses by hand so they can be documented. +-- makeLenses ''Camera + +-- | Camera position. +cameraPosition :: Lens' (Camera r) (Point 3 r) +cameraPosition = lens _cameraPosition (\cam p -> cam{_cameraPosition=p}) + +-- | Raw camera normal, i.e. a unit vector into the center of the screen. +rawCameraNormal :: Lens' (Camera r) (Vector 3 r) +rawCameraNormal = lens _rawCameraNormal (\cam r -> cam{_rawCameraNormal=r}) + +-- | Raw view up vector indicating which side points "upwards" in the scene. +rawViewUp :: Lens' (Camera r) (Vector 3 r) +rawViewUp = lens _rawViewUp (\cam r -> cam{_rawViewUp=r}) + +-- | Viewplane depth (i.e. the distance from the camera position to the plane on which we project). +viewPlaneDepth :: Lens' (Camera r) r +viewPlaneDepth = lens _viewPlaneDepth (\cam v -> cam{_viewPlaneDepth=v}) + +-- | Near distance (everything closer than this is clipped). +nearDist :: Lens' (Camera r) r +nearDist = lens _nearDist (\cam n -> cam{_nearDist=n}) + +-- | Far distance (everything further away than this is clipped). +farDist :: Lens' (Camera r) r +farDist = lens _farDist (\cam f -> cam{_farDist=f}) + +-- | Screen dimensions. +screenDimensions :: Lens' (Camera r) (Vector 2 r) +screenDimensions = lens _screenDimensions (\cam d -> cam{_screenDimensions=d}) + + +-------------------------------------------------------------------------------- +-- * Accessor Lenses + +-- | Lens to get and set the Camera normal, makes sure that the vector remains +-- normalized. +cameraNormal :: (Radical r, Fractional r) => Lens' (Camera r) (Vector 3 r) +cameraNormal = lens _rawCameraNormal (\c n -> c { _rawCameraNormal = signorm n} ) + + +-- | Lens to get and set the viewUp vector. Makes sure the vector remains +-- normalized. +viewUp :: (Radical r, Fractional r) => Lens' (Camera r) (Vector 3 r) +viewUp = lens _rawViewUp (\c n -> c { _rawViewUp = signorm n}) + + +-------------------------------------------------------------------------------- +-- * Camera Transformation functions + + +-- | Full transformation that renders the figure +cameraTransform :: Fractional r => Camera r -> Transformation 3 r +cameraTransform c = toViewPort c + |.| perspectiveProjection c + |.| worldToView c + +-- | Translates world coordinates into view coordinates +worldToView :: Fractional r => Camera r -> Transformation 3 r +worldToView c = rotateCoordSystem c |.| translation ((-1) *^ c^.cameraPosition.vector) + +-- | Transformation into viewport coordinates +toViewPort :: Fractional r => Camera r -> Transformation 3 r +toViewPort c = Transformation . Matrix + $ Vector4 (Vector4 (w/2) 0 0 0) + (Vector4 0 (h/2) 0 0) + (Vector4 0 0 (1/2) (1/2)) + (Vector4 0 0 0 1) + where + Vector2 w h = c^.screenDimensions + + +-- | constructs a perspective projection +perspectiveProjection :: Fractional r => Camera r -> Transformation 3 r +perspectiveProjection c = Transformation . Matrix $ + Vector4 (Vector4 (-n/rx) 0 0 0) + (Vector4 0 (-n/ry) 0 0) + (Vector4 0 0 (-(n+f)/(n-f)) (-2*n*f/(n-f))) + (Vector4 0 0 1 0) + where + n = c^.nearDist + f = c^.farDist + Vector2 rx ry = (/2) <$> c^.screenDimensions + +-- | Rotates coordinate system around the camera, such that we look in the negative z +-- direction +rotateCoordSystem :: Num r => Camera r -> Transformation 3 r +rotateCoordSystem c = rotateTo $ Vector3 u v n + where + u = (c^.rawViewUp) `cross` n + v = n `cross` u + n = (-1) *^ c^.rawCameraNormal -- we need the normal from the scene *into* the camera + + + + +-- transformBy' (Transformation m) (Vector3 x y z) = m `mult` (Vector4 x y z (-z)) + +-- | Flips the y and z axis. +flipAxes :: Num r => Transformation 3 r +flipAxes = Transformation . Matrix + $ Vector4 (Vector4 1 0 0 0) + (Vector4 0 0 1 0) + (Vector4 0 1 0 0) + (Vector4 0 0 0 1) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs index f8f20686f..61fe67213 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs @@ -48,7 +48,7 @@ asVertex planes defs = case definers defs of -- | test if v lies below (or on) all the given planes belowAll :: (Plane_ plane r, Ord r, Num r, Foldable f) => Point 3 r -> f plane -> Bool -belowAll v = all (\h -> onSideTest v h /= GT) +belowAll v = all (\h -> verticalSideTest v h /= GT) {-# INLINE belowAll #-} diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs index 5ebd50ced..dffc6b1d8 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs @@ -68,7 +68,7 @@ asVertex hs t@(Three h1 h2 h3) = do v <- intersectionPoint t -- | test if v lies below (or on) all the planes in hs belowAll :: (Plane_ plane r, Ord r, Num r, Foldable f) => Point 3 r -> f plane -> Bool -belowAll v = all (\h -> onSideTest v h /= GT) +belowAll v = all (\h -> verticalSideTest v h /= GT) {-# INLINE belowAll #-} diff --git a/hgeometry/test-with-ipe/test/BallSpec.hs b/hgeometry/test-with-ipe/test/BallSpec.hs new file mode 100644 index 000000000..f5048335e --- /dev/null +++ b/hgeometry/test-with-ipe/test/BallSpec.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE QuasiQuotes #-} +{-# LANGUAGE OverloadedStrings #-} +module BallSpec where + +import Data.Maybe +import Golden +import HGeometry.Ball +import HGeometry.HalfLine +import HGeometry.Intersection +import HGeometry.Line +import HGeometry.Point +import HGeometry.Vector +import Ipe +import Ipe.Color +import System.OsPath +import Test.Hspec +import Test.Hspec.WithTempFile +import Test.QuickCheck.Instances () + +-------------------------------------------------------------------------------- + +type R = Double + +spec :: Spec +spec = describe "ball intersection with line" $ do + goldenWith [osp|data/test-with-ipe/golden/|] + (ipeContentGolden { name = [osp|ball|] }) + (concat + [ map (\l -> iO'' l $ attr SLayer "lines" + ) lines' + , map (\hl -> iO'' hl $ attr SLayer "halfLines" + <> attr SStroke blue + ) halfLines + , map (\b -> iO'' b $ attr SLayer "balls" + ) balls + , map (\case + Line_x_Ball_Point q -> iO'' q $ attr SLayer "LineXBall" + Line_x_Ball_Segment seg -> iO'' seg $ attr SPen (IpePen "fat") + <> attr SLayer "LineXBall" + ) intersections + , map (\case + Line_x_Ball_Point q -> iO'' q $ attr SLayer "HalfLineXBall" + Line_x_Ball_Segment seg -> iO'' seg $ attr SPen (IpePen "fat") + <> attr SLayer "HalfLineXBall" + ) hlIntersections + ]) + +halfLines :: [HalfLine (Point 2 R)] +halfLines = [ HalfLine (Point2 4 0) (Vector2 1 (-2)) + , HalfLine (Point2 (-5) 10) (Vector2 1 (-5)) + ] + +lines' :: [LinePV 2 R] +lines' = [ LinePV (Point2 0 10) (Vector2 1 2) + , LinePV (Point2 (-100) 20) (Vector2 1 0) + , LinePV (Point2 (-100) 8) (Vector2 1 0) + ] + +balls :: [Ball (Point 2 R)] +balls = [ Ball origin 100 + , Ball (Point2 20 28) 64 + , Ball (Point2 200 28) 256 + ] + +intersections = catMaybes [ l `intersect` b | l <- lines', b <- balls ] + +hlIntersections = catMaybes [ l `intersect` b | l <- halfLines, b <- balls ] diff --git a/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs b/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs index 70f223d9b..fc4ae4ed2 100644 --- a/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs +++ b/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs @@ -163,37 +163,3 @@ testIpe inFp outFp = do goldenWith [osp|data/test-with-ipe/VoronoiDiagram/|] (ipeContentGolden { name = outFp }) out - - - - - --- testIpe :: OsPath -> IO [IpeObject R] --- testIpe inFp = do inFp' <- getDataFileName inFp --- (points :: [Point 2 R :+ _]) <- readAllFrom inFp' - --- print $ (Point3 183.02716 93.61106 8869.99979 :: Point 3 R) --- `onSideTest` --- (NonVerticalHyperPlane (Vector3 282 426 (-65250))) - - --- -- mapM_ print points --- -- mapM_ (print . liftPointToPlane . view core) points --- -- let hs = liftPointToPlane . view core <$> points --- -- mapM_ (print . asVertex hs) $ uniqueTriplets hs - --- let vv = voronoiVertices $ (view core) <$> points --- let vd = voronoiDiagram $ (view core) <$> points --- print $ vd - - --- -- print "vertices" --- -- mapM_ print vs --- pure $ [ iO' points --- , iO' vd --- ] <> [ iO'' v $ attr SStroke red | v <- vv ] - - --- -- $ (map iO' points) --- -- -- <> [iO' vd] --- -- <> diff --git a/hgeometry/test/Graphics/CameraSpec.hs b/hgeometry/test/Graphics/CameraSpec.hs new file mode 100644 index 000000000..b78377690 --- /dev/null +++ b/hgeometry/test/Graphics/CameraSpec.hs @@ -0,0 +1,64 @@ +module Graphics.CameraSpec where + +import Control.Lens +import HGeometry.Ext +import HGeometry.Graphics.Camera +import HGeometry.Point +import HGeometry.Transformation +import HGeometry.Triangle +import HGeometry.Vector +import Test.Hspec + +-------------------------------------------------------------------------------- + +spec :: Spec +spec = pure () + + +myCamera :: Camera Rational +myCamera = Camera (Point3 50 0 50) + (Vector3 0 0 (-1)) + (Vector3 0 1 0) + 10 + 15 + 55 + (Vector2 800 600) + + +myCamera1 :: Camera Double +myCamera1 = Camera origin + (Vector3 0 0 (-1)) + (Vector3 0 1 0) + 10 + 10 + 50 -- we can see up to the origin + (Vector2 60 40) + +testProjection :: Camera Double -> [Vector 3 Double] +testProjection c = map (transformBy t) [Vector3 30 30 (-10), Vector3 (30*50/10) 30 (-50)] + where + t = perspectiveProjection c + +myT :: Triangle (Point 3 Rational) +myT = Triangle (Point3 1 1 10) + (Point3 20 1 10) + (Point3 20 30 10) + + + +testToWorld :: Camera Double -> [Vector 3 Double] +testToWorld c = map (transformBy t) [u, v, n, Vector3 80 20 40] + where + u = (c^.rawViewUp) `cross` n + v = n `cross` u + n = (-1) *^ c^.rawCameraNormal -- we need the normal from the scene *into* the camera + t = worldToView c + + +testRotate :: Camera Double -> [Vector 3 Double] +testRotate c = map (transformBy t) [u, v, n] + where + u = (c^.rawViewUp) `cross` n + v = n `cross` u + n = (-1) *^ c^.rawCameraNormal -- we need the normal from the scene *into* the camera + t = rotateCoordSystem c diff --git a/hgeometry/test/LowerEnvelope/RegionsSpec.hs b/hgeometry/test/LowerEnvelope/RegionsSpec.hs index 557bdef16..221ed0f15 100644 --- a/hgeometry/test/LowerEnvelope/RegionsSpec.hs +++ b/hgeometry/test/LowerEnvelope/RegionsSpec.hs @@ -60,7 +60,7 @@ verifyOnPlane h1 h2 h3 = case intersectionPoint (Three h1 h2 h3) of -- | copied from the module belowAll :: (Plane_ plane r, Ord r, Num r, Foldable f) => Point 3 r -> f plane -> Bool -belowAll v = all (\h -> onSideTest v h /= GT) +belowAll v = all (\h -> verticalSideTest v h /= GT) myEnv = bruteForceLowerEnvelope inputs -- myTriEnv = triangulatedLowerEnvelope inputs diff --git a/todo.org b/todo.org index 265687118..751f4ae78 100644 --- a/todo.org +++ b/todo.org @@ -23,7 +23,7 @@ hopefully this will also make the CI story easier ** TODO move LINEQ to HGeometry.Line.NonVertical -* TODO fromPointAndNormal of NonVertical hyperplane +* DONE fromPointAndNormal of NonVertical hyperplane this can't be right: @@ -31,13 +31,11 @@ this can't be right: we should really use the in plane point p -* TODO normal vector business +* DONE normal vector business want: - normalVector points into the positive halfspace -- normalVector (fromPointAnNormal p n) == n - - onSideTest p