diff --git a/hgeometry-combinatorial/src/HGeometry/Sequence/Alternating.hs b/hgeometry-combinatorial/src/HGeometry/Sequence/Alternating.hs index 00c0f55d9..e6ae3327d 100644 --- a/hgeometry-combinatorial/src/HGeometry/Sequence/Alternating.hs +++ b/hgeometry-combinatorial/src/HGeometry/Sequence/Alternating.hs @@ -13,6 +13,7 @@ module HGeometry.Sequence.Alternating ( Alternating(..) , fromNonEmptyWith , mapF + , firstWithNeighbors , withNeighbours , mergeAlternating , insertBreakPoints @@ -21,6 +22,7 @@ module HGeometry.Sequence.Alternating , consElemWith , unconsAlt , snocElemWith + , separators ) where import Control.DeepSeq @@ -194,3 +196,19 @@ snocElemWith f (Alternating x0 xs) y = Alternating x0 $ view (re _Snoc) (xs, (s, s = case xs^?_last of Nothing -> f x0 y Just (_,x) -> f x y + + +-------------------------------------------------------------------------------- + +-- | Map over the separators of the alterating together with the neighbours +firstWithNeighbors :: Traversable f + => (a -> sep -> a -> sep') + -> Alternating f sep a -> Alternating f sep' a +firstWithNeighbors f (Alternating x0 xs) = Alternating x0 xs' + where + (_, xs') = List.mapAccumL (\x (sep,y) -> (y, (f x sep y, y))) x0 xs + + +-- | Get the separators out of the alternating +separators :: Functor f => Alternating f sep a -> f sep +separators (Alternating _ xs) = fmap fst xs diff --git a/hgeometry-examples/skia/Model.hs b/hgeometry-examples/skia/Model.hs index 006766371..76418bbcd 100644 --- a/hgeometry-examples/skia/Model.hs +++ b/hgeometry-examples/skia/Model.hs @@ -28,6 +28,7 @@ import Color import Control.Lens hiding (view, element) import Data.Default.Class import qualified Data.IntMap as IntMap +import qualified Data.Set as Set import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import Data.Maybe (fromJust) @@ -94,7 +95,7 @@ data Model = Model { _canvas :: SkiaCanvas.Canvas R , _planeGraphs :: IntMap.IntMap (PlaneGraph' R) - , _diagram :: Maybe [Point 2 R] + , _diagram :: Maybe (Set.Set (Point 2 R)) , __layers :: Layers , _stroke :: !StrokeFill , _fill :: !StrokeFill diff --git a/hgeometry-examples/voronoiDiagram/Main.hs b/hgeometry-examples/voronoiDiagram/Main.hs index 1bc3e19d1..92a0c5215 100644 --- a/hgeometry-examples/voronoiDiagram/Main.hs +++ b/hgeometry-examples/voronoiDiagram/Main.hs @@ -3,11 +3,12 @@ module Main(main) where import Control.Lens hiding (view, element) +import Data.Foldable (toList) import qualified Data.IntMap as IntMap import qualified Data.List.NonEmpty as NonEmpty import qualified Data.Map as Map +import qualified Data.Set as Set import GHC.TypeNats -import HGeometry.VoronoiDiagram import HGeometry.Ext import HGeometry.Miso.OrphanInstances () import HGeometry.Miso.Svg @@ -15,6 +16,7 @@ import HGeometry.Miso.Svg.Canvas (Canvas, blankCanvas, mouseCoordinate import qualified HGeometry.Miso.Svg.Canvas as Canvas import HGeometry.Number.Real.Rational import HGeometry.Point +import HGeometry.VoronoiDiagram import qualified Language.Javascript.JSaddle.Warp as JSaddle import Miso import Miso.String (MisoString,ToMisoString(..), ms) @@ -26,7 +28,7 @@ type R = RealNumber 5 data Model = Model { _canvas :: Canvas R , _points :: IntMap.IntMap (Point 2 R) - , _diagram :: Maybe [Point 2 R] + , _diagram :: Maybe (Set.Set (Point 2 R)) } deriving (Eq) makeLenses ''Model @@ -96,7 +98,7 @@ viewModel m = div_ [ ] canvasBody = [ g_ [] [ draw v [ fill_ "red" ] ] - | v <- m^..diagram.traverse.traverse ] + | v <- m^..diagram.traverse.to toList.traverse ] <> [ g_ [] [ draw p [ fill_ "black" ] , textAt p [] (ms i) diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate1_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate1_out.ipe index 43321fae7..321205bb4 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate1_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate1_out.ipe @@ -143,20 +143,20 @@ -1096.000000000000 -904.000000000000 m -1096.000000000000 1096.000000000000 l +-904.000000000000 1096.000000000000 m +-904.000000000000 -904.000000000000 l 96.000000000000 96.000000000000 l h -1096.000000000000 1096.000000000000 m --904.000000000000 1096.000000000000 l +-904.000000000000 -904.000000000000 m +1096.000000000000 -904.000000000000 l 96.000000000000 96.000000000000 l h --904.000000000000 -904.000000000000 m -1096.000000000000 -904.000000000000 l +1096.000000000000 1096.000000000000 m +-904.000000000000 1096.000000000000 l 96.000000000000 96.000000000000 l h --904.000000000000 1096.000000000000 m --904.000000000000 -904.000000000000 l +1096.000000000000 -904.000000000000 m +1096.000000000000 1096.000000000000 l 96.000000000000 96.000000000000 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate2_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate2_out.ipe index 40f54ddbd..edf26ea35 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate2_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate2_out.ipe @@ -143,46 +143,46 @@ -1128.000225351379 -286.210430677666 m -1128.000566767257 542.208941532476 l -128.000566767257 128.000566767257 l -128.000019914361 128.000019914361 l +-872.000225351380 542.213019517769 m +-872.000048077500 -286.214576669489 l +127.999951922500 128.000048077500 l +127.999774648620 128.000225351379 l +h +-872.000048077500 -286.214576669489 m +-871.999951922500 -2286.207418497661 l 128.000048077500 127.999951922500 l -128.000225351379 127.999774648620 l +128.000019914361 128.000019914361 l +127.999951922500 128.000048077500 l h -1128.000566767257 542.208941532476 m -1128.000566767257 2542.244365113483 l -128.000566767257 128.000566767257 l +-872.000225351380 2542.233354185735 m +-872.000225351380 542.213019517769 l +127.999774648620 128.000225351379 l h -1128.000225351379 -2286.218265169573 m -1128.000225351379 -286.210430677666 l +-871.999951922500 -2286.207418497661 m +1128.000225351379 -2286.218265169573 l 128.000225351379 127.999774648620 l +128.000048077500 127.999951922500 l h -1128.000566767257 2542.244365113483 m +1128.000566767257 2542.244365113483 m -872.000225351380 2542.233354185735 l 127.999774648620 128.000225351379 l 127.999951922500 128.000048077500 l 128.000019914361 128.000019914361 l 128.000566767257 128.000566767257 l h --871.999951922500 -2286.207418497661 m -1128.000225351379 -2286.218265169573 l +1128.000225351379 -2286.218265169573 m +1128.000225351379 -286.210430677666 l 128.000225351379 127.999774648620 l -128.000048077500 127.999951922500 l h --872.000225351380 2542.233354185735 m --872.000225351380 542.213019517769 l -127.999774648620 128.000225351379 l +1128.000566767257 542.208941532476 m +1128.000566767257 2542.244365113483 l +128.000566767257 128.000566767257 l h --872.000048077500 -286.214576669489 m --871.999951922500 -2286.207418497661 l -128.000048077500 127.999951922500 l +1128.000225351379 -286.210430677666 m +1128.000566767257 542.208941532476 l +128.000566767257 128.000566767257 l 128.000019914361 128.000019914361 l -127.999951922500 128.000048077500 l -h --872.000225351380 542.213019517769 m --872.000048077500 -286.214576669489 l -127.999951922500 128.000048077500 l -127.999774648620 128.000225351379 l +128.000048077500 127.999951922500 l +128.000225351379 127.999774648620 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate_out.ipe index 16fbaaf23..60b0e6c2b 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/degenerate_out.ipe @@ -143,20 +143,20 @@ -1096.000000000000 96.000000000000 m -96.000000000000 1096.000000000000 l +-904.000000000000 96.000000000000 m +96.000000000000 -904.000000000000 l 96.000000000000 96.000000000000 l h -96.000000000000 -904.000000000000 m -1096.000000000000 96.000000000000 l +96.000000000000 1096.000000000000 m +-904.000000000000 96.000000000000 l 96.000000000000 96.000000000000 l h -96.000000000000 1096.000000000000 m --904.000000000000 96.000000000000 l +96.000000000000 -904.000000000000 m +1096.000000000000 96.000000000000 l 96.000000000000 96.000000000000 l h --904.000000000000 96.000000000000 m -96.000000000000 -904.000000000000 l +1096.000000000000 96.000000000000 m +96.000000000000 1096.000000000000 l 96.000000000000 96.000000000000 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/foo_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/foo_out.ipe index 552094146..99d145a49 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/foo_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/foo_out.ipe @@ -143,48 +143,47 @@ -1464.975609756097 -2346.439024390244 m -1481.313432835820 290.812603648424 l -481.313432835820 179.701492537313 l -464.975609756097 153.560975609756 l +-889.000000000000 398.250000000000 m +-888.000000000000 -28312.000000000000 l +112.000000000000 -3312.000000000000 l +112.000000000000 146.000000000000 l +111.000000000000 148.250000000000 l h -1481.313432835820 290.812603648424 m -1427.534883720930 552.883720930232 l -427.534883720930 302.883720930232 l -425.254901960784 277.803921568627 l -481.313432835820 179.701492537313 l +-892.444444444445 705.777777777777 m +-889.000000000000 398.250000000000 l +111.000000000000 148.250000000000 l +154.285714285714 224.000000000000 l +107.555555555555 305.777777777777 l h -1427.534883720930 552.883720930232 m -1375.292035398230 3120.053097345132 l -375.292035398230 370.053097345132 l -427.534883720930 302.883720930232 l +-758.666666666667 2872.666666666666 m +-892.444444444445 705.777777777777 l +107.555555555555 305.777777777777 l +241.333333333333 372.666666666666 l h --888.000000000000 -28312.000000000000 m -1464.975609756097 -2346.439024390244 l -464.975609756097 153.560975609756 l -388.075471698113 120.603773584905 l -248.000000000000 -1000.000000000000 l +248.000000000000 -1000.000000000000 m +248.000000000000 112.000000000000 l +112.000000000000 146.000000000000 l 112.000000000000 -3312.000000000000 l h -425.254901960784 277.803921568627 m -341.913043478260 211.130434782608 l -325.473684210526 172.771929824561 l -388.075471698113 120.603773584905 l -464.975609756097 153.560975609756 l -481.313432835820 179.701492537313 l +272.500000000000 161.000000000000 m +220.000000000000 224.000000000000 l +154.285714285714 224.000000000000 l +111.000000000000 148.250000000000 l +112.000000000000 146.000000000000 l +248.000000000000 112.000000000000 l h -375.292035398230 370.053097345132 m +268.571428571428 321.142857142857 m 271.917525773195 351.257731958762 l -268.571428571428 321.142857142857 l -341.913043478260 211.130434782608 l -425.254901960784 277.803921568627 l -427.534883720930 302.883720930232 l +241.333333333333 372.666666666666 l +107.555555555555 305.777777777777 l +154.285714285714 224.000000000000 l +220.000000000000 224.000000000000 l h -388.075471698113 120.603773584905 m -325.473684210526 172.771929824561 l +268.571428571428 321.142857142857 m +220.000000000000 224.000000000000 l 272.500000000000 161.000000000000 l -248.000000000000 112.000000000000 l -248.000000000000 -1000.000000000000 l +325.473684210526 172.771929824561 l +341.913043478260 211.130434782608 l h 1375.292035398230 3120.053097345132 m -758.666666666667 2872.666666666666 l @@ -192,46 +191,47 @@ h 271.917525773195 351.257731958762 l 375.292035398230 370.053097345132 l h -268.571428571428 321.142857142857 m -220.000000000000 224.000000000000 l -272.500000000000 161.000000000000 l +388.075471698113 120.603773584905 m 325.473684210526 172.771929824561 l -341.913043478260 211.130434782608 l +272.500000000000 161.000000000000 l +248.000000000000 112.000000000000 l +248.000000000000 -1000.000000000000 l h -268.571428571428 321.142857142857 m +375.292035398230 370.053097345132 m 271.917525773195 351.257731958762 l -241.333333333333 372.666666666666 l -107.555555555555 305.777777777777 l -154.285714285714 224.000000000000 l -220.000000000000 224.000000000000 l +268.571428571428 321.142857142857 l +341.913043478260 211.130434782608 l +425.254901960784 277.803921568627 l +427.534883720930 302.883720930232 l h -272.500000000000 161.000000000000 m -220.000000000000 224.000000000000 l -154.285714285714 224.000000000000 l -111.000000000000 148.250000000000 l -112.000000000000 146.000000000000 l -248.000000000000 112.000000000000 l +425.254901960784 277.803921568627 m +341.913043478260 211.130434782608 l +325.473684210526 172.771929824561 l +388.075471698113 120.603773584905 l +464.975609756097 153.560975609756 l +481.313432835820 179.701492537313 l h -248.000000000000 -1000.000000000000 m -248.000000000000 112.000000000000 l -112.000000000000 146.000000000000 l +-888.000000000000 -28312.000000000000 m +1464.975609756097 -2346.439024390244 l +464.975609756097 153.560975609756 l +388.075471698113 120.603773584905 l +248.000000000000 -1000.000000000000 l 112.000000000000 -3312.000000000000 l h --758.666666666667 2872.666666666666 m --892.444444444445 705.777777777777 l -107.555555555555 305.777777777777 l -241.333333333333 372.666666666666 l +1427.534883720930 552.883720930232 m +1375.292035398230 3120.053097345132 l +375.292035398230 370.053097345132 l +427.534883720930 302.883720930232 l h --892.444444444445 705.777777777777 m --889.000000000000 398.250000000000 l -111.000000000000 148.250000000000 l -154.285714285714 224.000000000000 l -107.555555555555 305.777777777777 l +1481.313432835820 290.812603648424 m +1427.534883720930 552.883720930232 l +427.534883720930 302.883720930232 l +425.254901960784 277.803921568627 l +481.313432835820 179.701492537313 l h --889.000000000000 398.250000000000 m --888.000000000000 -28312.000000000000 l -112.000000000000 -3312.000000000000 l -112.000000000000 146.000000000000 l -111.000000000000 148.250000000000 l +1464.975609756097 -2346.439024390244 m +1481.313432835820 290.812603648424 l +481.313432835820 179.701492537313 l +464.975609756097 153.560975609756 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple1_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple1_out.ipe index 283197295..1d33ec548 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple1_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple1_out.ipe @@ -143,21 +143,21 @@ --821.763977197983 -15225.338851129139 m -1274.778112187036 187.289342337422 l -274.778112187036 175.728648695803 l +-916.167328699107 1463.621598088741 m +-907.769230769231 -584.076923076924 l +92.230769230769 165.923076923076 l +83.832671300893 217.244786494538 l +h +-907.769230769231 -584.076923076924 m +-821.763977197983 -15225.338851129139 l 178.236022802017 74.661148870861 l +128.529411764705 153.823529411764 l +92.230769230769 165.923076923076 l h -1274.778112187036 187.289342337422 m -1185.530696065715 68765.102680501513 l +1185.530696065715 68765.102680501513 m +-916.167328699107 1463.621598088741 l +83.832671300893 217.244786494538 l 185.530696065715 265.102680501513 l -206.675427069645 216.340341655716 l -274.778112187036 175.728648695803 l -h -274.778112187036 175.728648695803 m -206.675427069645 216.340341655716 l -128.529411764705 153.823529411764 l -178.236022802017 74.661148870861 l h 206.675427069645 216.340341655716 m 185.530696065715 265.102680501513 l @@ -165,20 +165,20 @@ h 92.230769230769 165.923076923076 l 128.529411764705 153.823529411764 l h -1185.530696065715 68765.102680501513 m --916.167328699107 1463.621598088741 l -83.832671300893 217.244786494538 l +274.778112187036 175.728648695803 m +206.675427069645 216.340341655716 l +128.529411764705 153.823529411764 l +178.236022802017 74.661148870861 l +h +1274.778112187036 187.289342337422 m +1185.530696065715 68765.102680501513 l 185.530696065715 265.102680501513 l +206.675427069645 216.340341655716 l +274.778112187036 175.728648695803 l h --907.769230769231 -584.076923076924 m --821.763977197983 -15225.338851129139 l +-821.763977197983 -15225.338851129139 m +1274.778112187036 187.289342337422 l +274.778112187036 175.728648695803 l 178.236022802017 74.661148870861 l -128.529411764705 153.823529411764 l -92.230769230769 165.923076923076 l -h --916.167328699107 1463.621598088741 m --907.769230769231 -584.076923076924 l -92.230769230769 165.923076923076 l -83.832671300893 217.244786494538 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple_out.ipe index 35900beb3..748895846 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simple_out.ipe @@ -143,28 +143,28 @@ --725.221887812964 -871.146351304197 m -1274.778112187036 187.289342337422 l -274.778112187036 175.728648695803 l -h -1274.778112187036 187.289342337422 m -1185.530696065715 68765.102680501513 l +1185.530696065715 68765.102680501513 m +-814.469303934285 -205.485554792605 l 185.530696065715 265.102680501513 l +h +-814.469303934285 -205.485554792605 m +-793.324572930355 -583.659658344284 l 206.675427069645 216.340341655716 l -274.778112187036 175.728648695803 l +185.530696065715 265.102680501513 l h -793.324572930355 -583.659658344284 m -725.221887812964 -871.146351304197 l 274.778112187036 175.728648695803 l 206.675427069645 216.340341655716 l h --814.469303934285 -205.485554792605 m --793.324572930355 -583.659658344284 l -206.675427069645 216.340341655716 l +1274.778112187036 187.289342337422 m +1185.530696065715 68765.102680501513 l 185.530696065715 265.102680501513 l +206.675427069645 216.340341655716 l +274.778112187036 175.728648695803 l h -1185.530696065715 68765.102680501513 m --814.469303934285 -205.485554792605 l -185.530696065715 265.102680501513 l +-725.221887812964 -871.146351304197 m +1274.778112187036 187.289342337422 l +274.778112187036 175.728648695803 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simpler_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simpler_out.ipe index 43ecb975f..6af0d6660 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simpler_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simpler_out.ipe @@ -143,22 +143,22 @@ -1092.230769230769 -167.410256410257 m -1083.832671300893 687.833021788655 l -83.832671300893 217.244786494538 l +-916.167328699107 1463.621598088741 m +-907.769230769231 -584.076923076924 l 92.230769230769 165.923076923076 l -h -1083.832671300893 687.833021788655 m --916.167328699107 1463.621598088741 l 83.832671300893 217.244786494538 l h --907.769230769231 -584.076923076924 m +-907.769230769231 -584.076923076924 m 1092.230769230769 -167.410256410257 l 92.230769230769 165.923076923076 l h --916.167328699107 1463.621598088741 m --907.769230769231 -584.076923076924 l -92.230769230769 165.923076923076 l +1083.832671300893 687.833021788655 m +-916.167328699107 1463.621598088741 l +83.832671300893 217.244786494538 l +h +1092.230769230769 -167.410256410257 m +1083.832671300893 687.833021788655 l 83.832671300893 217.244786494538 l +92.230769230769 165.923076923076 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simplest_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simplest_out.ipe index 78a638b8a..a5fa08cd4 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simplest_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/simplest_out.ipe @@ -143,16 +143,16 @@ -1112.695099818511 96.422840015398 m --887.304900181489 1427.648136458086 l +-887.304900181489 1427.648136458086 m +-887.304900181489 -568.728675136117 l 112.695099818511 181.271324863883 l h -887.304900181489 -568.728675136117 m 1112.695099818511 96.422840015398 l 112.695099818511 181.271324863883 l h --887.304900181489 1427.648136458086 m --887.304900181489 -568.728675136117 l +1112.695099818511 96.422840015398 m +-887.304900181489 1427.648136458086 l 112.695099818511 181.271324863883 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/trivial_out.ipe b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/trivial_out.ipe index 8da9b0e9c..57d5e549b 100644 --- a/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/trivial_out.ipe +++ b/hgeometry/data/test-with-ipe/Plane/LowerEnvelope/trivial_out.ipe @@ -143,16 +143,16 @@ -1337.545454545454 236.181818181818 m --662.454545454546 18636.181818181818 l +-662.454545454546 18636.181818181818 m +-662.454545454546 -919.373737373738 l 337.545454545454 636.181818181818 l h -662.454545454546 -919.373737373738 m 1337.545454545454 236.181818181818 l 337.545454545454 636.181818181818 l h --662.454545454546 18636.181818181818 m --662.454545454546 -919.373737373738 l +1337.545454545454 236.181818181818 m +-662.454545454546 18636.181818181818 l 337.545454545454 636.181818181818 l h \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/foo_out.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/foo_out.ipe index b0688f298..99d145a49 100644 --- a/hgeometry/data/test-with-ipe/VoronoiDiagram/foo_out.ipe +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/foo_out.ipe @@ -1,5 +1,15 @@ - + + + + + + + + + + + @@ -133,72 +143,95 @@ -464.975609756097 153.560975609756 m -1000.000000000000 -1184.000000000000 l -481.313432835820 179.701492537313 m -1000.000000000000 237.333333333333 l -427.534883720930 302.883720930232 m -1000.000000000000 446.000000000000 l -112.000000000000 -3312.000000000000 m -111.960000000000 -3313.000000000000 l -375.292035398230 370.053097345132 m -604.363636363636 1000.000000000000 l -241.333333333333 372.666666666666 m --9.600000000000 1000.000000000000 l -107.555555555555 305.777777777777 m --1000.000000000000 748.800000000000 l -111.000000000000 148.250000000000 m --1000.000000000000 426.000000000000 l -107.555555555555 305.777777777777 m +-889.000000000000 398.250000000000 m +-888.000000000000 -28312.000000000000 l +112.000000000000 -3312.000000000000 l +112.000000000000 146.000000000000 l +111.000000000000 148.250000000000 l +h +-892.444444444445 705.777777777777 m +-889.000000000000 398.250000000000 l +111.000000000000 148.250000000000 l 154.285714285714 224.000000000000 l -107.555555555555 305.777777777777 m +107.555555555555 305.777777777777 l +h +-758.666666666667 2872.666666666666 m +-892.444444444445 705.777777777777 l +107.555555555555 305.777777777777 l 241.333333333333 372.666666666666 l -111.000000000000 148.250000000000 m -154.285714285714 224.000000000000 l -111.000000000000 148.250000000000 m +h +248.000000000000 -1000.000000000000 m +248.000000000000 112.000000000000 l 112.000000000000 146.000000000000 l -112.000000000000 -3312.000000000000 m +112.000000000000 -3312.000000000000 l +h +272.500000000000 161.000000000000 m +220.000000000000 224.000000000000 l +154.285714285714 224.000000000000 l +111.000000000000 148.250000000000 l 112.000000000000 146.000000000000 l -112.000000000000 -3312.000000000000 m -248.000000000000 -1000.000000000000 l -112.000000000000 146.000000000000 m 248.000000000000 112.000000000000 l -154.285714285714 224.000000000000 m -220.000000000000 224.000000000000 l -220.000000000000 224.000000000000 m -272.500000000000 161.000000000000 l -220.000000000000 224.000000000000 m -268.571428571428 321.142857142857 l -241.333333333333 372.666666666666 m +h +268.571428571428 321.142857142857 m 271.917525773195 351.257731958762 l -248.000000000000 -1000.000000000000 m -388.075471698113 120.603773584905 l -248.000000000000 -1000.000000000000 m -248.000000000000 112.000000000000 l -248.000000000000 112.000000000000 m +241.333333333333 372.666666666666 l +107.555555555555 305.777777777777 l +154.285714285714 224.000000000000 l +220.000000000000 224.000000000000 l +h +268.571428571428 321.142857142857 m +220.000000000000 224.000000000000 l 272.500000000000 161.000000000000 l -268.571428571428 321.142857142857 m +325.473684210526 172.771929824561 l 341.913043478260 211.130434782608 l -268.571428571428 321.142857142857 m +h +1375.292035398230 3120.053097345132 m +-758.666666666667 2872.666666666666 l +241.333333333333 372.666666666666 l 271.917525773195 351.257731958762 l -271.917525773195 351.257731958762 m 375.292035398230 370.053097345132 l -272.500000000000 161.000000000000 m +h +388.075471698113 120.603773584905 m 325.473684210526 172.771929824561 l -325.473684210526 172.771929824561 m +272.500000000000 161.000000000000 l +248.000000000000 112.000000000000 l +248.000000000000 -1000.000000000000 l +h +375.292035398230 370.053097345132 m +271.917525773195 351.257731958762 l +268.571428571428 321.142857142857 l 341.913043478260 211.130434782608 l -325.473684210526 172.771929824561 m -388.075471698113 120.603773584905 l -341.913043478260 211.130434782608 m 425.254901960784 277.803921568627 l -375.292035398230 370.053097345132 m 427.534883720930 302.883720930232 l -388.075471698113 120.603773584905 m +h +425.254901960784 277.803921568627 m +341.913043478260 211.130434782608 l +325.473684210526 172.771929824561 l +388.075471698113 120.603773584905 l 464.975609756097 153.560975609756 l -425.254901960784 277.803921568627 m +481.313432835820 179.701492537313 l +h +-888.000000000000 -28312.000000000000 m +1464.975609756097 -2346.439024390244 l +464.975609756097 153.560975609756 l +388.075471698113 120.603773584905 l +248.000000000000 -1000.000000000000 l +112.000000000000 -3312.000000000000 l +h +1427.534883720930 552.883720930232 m +1375.292035398230 3120.053097345132 l +375.292035398230 370.053097345132 l +427.534883720930 302.883720930232 l +h +1481.313432835820 290.812603648424 m +1427.534883720930 552.883720930232 l 427.534883720930 302.883720930232 l -425.254901960784 277.803921568627 m +425.254901960784 277.803921568627 l 481.313432835820 179.701492537313 l -464.975609756097 153.560975609756 m +h +1464.975609756097 -2346.439024390244 m +1481.313432835820 290.812603648424 l 481.313432835820 179.701492537313 l - \ No newline at end of file +464.975609756097 153.560975609756 l +h + \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/simple1_out.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/simple1_out.ipe index 9fe5b9d1a..1d33ec548 100644 --- a/hgeometry/data/test-with-ipe/VoronoiDiagram/simple1_out.ipe +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/simple1_out.ipe @@ -1,5 +1,15 @@ - + + + + + + + + + + + @@ -133,30 +143,42 @@ -178.236022802017 74.661148870861 m -107.996732026143 -1000.000000000000 l -274.778112187036 175.728648695803 m -1000.000000000000 184.112716763005 l -185.530696065715 265.102680501513 m -196.259124087591 1000.000000000000 l -92.230769230769 165.923076923076 m --1000.000000000000 -653.250000000000 l -83.832671300893 217.244786494538 m --544.191860465117 1000.000000000000 l -83.832671300893 217.244786494538 m +-916.167328699107 1463.621598088741 m +-907.769230769231 -584.076923076924 l 92.230769230769 165.923076923076 l -83.832671300893 217.244786494538 m +83.832671300893 217.244786494538 l +h +-907.769230769231 -584.076923076924 m +-821.763977197983 -15225.338851129139 l +178.236022802017 74.661148870861 l +128.529411764705 153.823529411764 l +92.230769230769 165.923076923076 l +h +1185.530696065715 68765.102680501513 m +-916.167328699107 1463.621598088741 l +83.832671300893 217.244786494538 l 185.530696065715 265.102680501513 l -92.230769230769 165.923076923076 m +h +206.675427069645 216.340341655716 m +185.530696065715 265.102680501513 l +83.832671300893 217.244786494538 l +92.230769230769 165.923076923076 l +128.529411764705 153.823529411764 l +h +274.778112187036 175.728648695803 m +206.675427069645 216.340341655716 l 128.529411764705 153.823529411764 l -128.529411764705 153.823529411764 m 178.236022802017 74.661148870861 l -128.529411764705 153.823529411764 m +h +1274.778112187036 187.289342337422 m +1185.530696065715 68765.102680501513 l +185.530696065715 265.102680501513 l 206.675427069645 216.340341655716 l -178.236022802017 74.661148870861 m 274.778112187036 175.728648695803 l -185.530696065715 265.102680501513 m -206.675427069645 216.340341655716 l -206.675427069645 216.340341655716 m +h +-821.763977197983 -15225.338851129139 m +1274.778112187036 187.289342337422 l 274.778112187036 175.728648695803 l - \ No newline at end of file +178.236022802017 74.661148870861 l +h + \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/simple_out.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/simple_out.ipe index 0b3eb2fb1..748895846 100644 --- a/hgeometry/data/test-with-ipe/VoronoiDiagram/simple_out.ipe +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/simple_out.ipe @@ -1,5 +1,15 @@ - + + + + + + + + + + + @@ -133,18 +143,28 @@ -274.778112187036 175.728648695803 m --848.305970149254 -1000.000000000000 l -274.778112187036 175.728648695803 m -1000.000000000000 184.112716763005 l -206.675427069645 216.340341655716 m --1000.000000000000 -749.000000000000 l -185.530696065715 265.102680501513 m --1000.000000000000 -292.794117647059 l -185.530696065715 265.102680501513 m -196.259124087591 1000.000000000000 l -185.530696065715 265.102680501513 m +1185.530696065715 68765.102680501513 m +-814.469303934285 -205.485554792605 l +185.530696065715 265.102680501513 l +h +-814.469303934285 -205.485554792605 m +-793.324572930355 -583.659658344284 l 206.675427069645 216.340341655716 l -206.675427069645 216.340341655716 m +185.530696065715 265.102680501513 l +h +-793.324572930355 -583.659658344284 m +-725.221887812964 -871.146351304197 l 274.778112187036 175.728648695803 l - \ No newline at end of file +206.675427069645 216.340341655716 l +h +1274.778112187036 187.289342337422 m +1185.530696065715 68765.102680501513 l +185.530696065715 265.102680501513 l +206.675427069645 216.340341655716 l +274.778112187036 175.728648695803 l +h +-725.221887812964 -871.146351304197 m +1274.778112187036 187.289342337422 l +274.778112187036 175.728648695803 l +h + \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/simpler_out.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/simpler_out.ipe index 2836da311..6af0d6660 100644 --- a/hgeometry/data/test-with-ipe/VoronoiDiagram/simpler_out.ipe +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/simpler_out.ipe @@ -1,5 +1,15 @@ - + + + + + + + + + + + @@ -133,14 +143,22 @@ -92.230769230769 165.923076923076 m -1000.000000000000 -136.666666666667 l -83.832671300893 217.244786494538 m -1000.000000000000 648.382352941176 l -92.230769230769 165.923076923076 m --1000.000000000000 -653.250000000000 l -83.832671300893 217.244786494538 m --544.191860465117 1000.000000000000 l -83.832671300893 217.244786494538 m +-916.167328699107 1463.621598088741 m +-907.769230769231 -584.076923076924 l 92.230769230769 165.923076923076 l - \ No newline at end of file +83.832671300893 217.244786494538 l +h +-907.769230769231 -584.076923076924 m +1092.230769230769 -167.410256410257 l +92.230769230769 165.923076923076 l +h +1083.832671300893 687.833021788655 m +-916.167328699107 1463.621598088741 l +83.832671300893 217.244786494538 l +h +1092.230769230769 -167.410256410257 m +1083.832671300893 687.833021788655 l +83.832671300893 217.244786494538 l +92.230769230769 165.923076923076 l +h + \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/simplest_out.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/simplest_out.ipe index 17cfb79cf..a5fa08cd4 100644 --- a/hgeometry/data/test-with-ipe/VoronoiDiagram/simplest_out.ipe +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/simplest_out.ipe @@ -1,5 +1,15 @@ - + + + + + + + + + + + @@ -133,10 +143,16 @@ -112.695099818511 181.271324863883 m -1000.000000000000 105.984848484848 l -112.695099818511 181.271324863883 m --1000.000000000000 -653.250000000000 l -112.695099818511 181.271324863883 m --544.191860465117 1000.000000000000 l - \ No newline at end of file +-887.304900181489 1427.648136458086 m +-887.304900181489 -568.728675136117 l +112.695099818511 181.271324863883 l +h +-887.304900181489 -568.728675136117 m +1112.695099818511 96.422840015398 l +112.695099818511 181.271324863883 l +h +1112.695099818511 96.422840015398 m +-887.304900181489 1427.648136458086 l +112.695099818511 181.271324863883 l +h + \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/trivialVoronoi.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/trivialVoronoi.ipe new file mode 100644 index 000000000..5169117ac --- /dev/null +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/trivialVoronoi.ipe @@ -0,0 +1,148 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +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 + + + + + + +-995.000000000000 1005.000000000000 m +5.000000000000 -995.000000000000 l +5.000000000000 5.000000000000 l +h +5.000000000000 -995.000000000000 m +1005.000000000000 5.000000000000 l +5.000000000000 5.000000000000 l +h +1005.000000000000 5.000000000000 m +-995.000000000000 1005.000000000000 l +5.000000000000 5.000000000000 l +h + \ No newline at end of file diff --git a/hgeometry/data/test-with-ipe/VoronoiDiagram/trivial_out.ipe b/hgeometry/data/test-with-ipe/VoronoiDiagram/trivial_out.ipe index 73de28fdc..57d5e549b 100644 --- a/hgeometry/data/test-with-ipe/VoronoiDiagram/trivial_out.ipe +++ b/hgeometry/data/test-with-ipe/VoronoiDiagram/trivial_out.ipe @@ -1,5 +1,15 @@ - + + + + + + + + + + + @@ -133,10 +143,16 @@ -337.545454545454 636.181818181818 m -1000.000000000000 371.200000000000 l -337.545454545454 636.181818181818 m --714.285714285715 -1000.000000000000 l -337.545454545454 636.181818181818 m -317.333333333333 1000.000000000000 l - \ No newline at end of file +-662.454545454546 18636.181818181818 m +-662.454545454546 -919.373737373738 l +337.545454545454 636.181818181818 l +h +-662.454545454546 -919.373737373738 m +1337.545454545454 236.181818181818 l +337.545454545454 636.181818181818 l +h +1337.545454545454 236.181818181818 m +-662.454545454546 18636.181818181818 l +337.545454545454 636.181818181818 l +h + \ No newline at end of file diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal index 896293ae5..90d84bfe0 100644 --- a/hgeometry/hgeometry.cabal +++ b/hgeometry/hgeometry.cabal @@ -76,6 +76,7 @@ common all-setup , HsYAML >= 0.2 && < 1 , semialign >= 1.3 && < 1.4 , these >= 1.0.1 && < 1.3 + , subcategories >= 0.2 && < 0.3 , ghc-typelits-natnormalise >= 0.7.7 && < 1 , ghc-typelits-knownnat >= 0.7.6 && < 1 @@ -389,8 +390,8 @@ library HGeometry.ConvexHull.Melkman - HGeometry.ConvexHull.R3.Naive - HGeometry.ConvexHull.R3.Naive.Dual + -- HGeometry.ConvexHull.R3.Naive + -- HGeometry.ConvexHull.R3.Naive.Dual HGeometry.BezierSpline @@ -413,12 +414,7 @@ library HGeometry.Plane.LowerEnvelope HGeometry.Plane.LowerEnvelope.Naive - HGeometry.Plane.LowerEnvelope.VertexForm - HGeometry.Plane.LowerEnvelope.AdjListForm HGeometry.Plane.LowerEnvelope.Connected - - -- newer attempt: - HGeometry.Plane.LowerEnvelope.ConnectedNew HGeometry.Plane.LowerEnvelope.Connected.Graph @@ -452,22 +448,17 @@ library HGeometry.LineSegment.Intersection.Types HGeometry.Plane.LowerEnvelope.Type - HGeometry.Plane.LowerEnvelope.Connected.Type - HGeometry.Plane.LowerEnvelope.Connected.FromVertexForm -- HGeometry.Plane.LowerEnvelope.Sample - - --- files corresponding to the newer attempt + HGeometry.Plane.LowerEnvelope.Connected.Type HGeometry.Plane.LowerEnvelope.Connected.Primitives + HGeometry.Plane.LowerEnvelope.Connected.VertexForm HGeometry.Plane.LowerEnvelope.Connected.Regions HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap HGeometry.Plane.LowerEnvelope.Connected.BruteForce HGeometry.Plane.LowerEnvelope.Connected.Separator - HGeometry.Plane.LowerEnvelope.DivideAndConquer - HGeometry.Plane.LowerEnvelope.EpsApproximation - - -- HGeometry.Plane.LowerEnvelope.AtMostThree - -- HGeometry.Plane.LowerEnvelope.Triangulate + -- HGeometry.Plane.LowerEnvelope.DivideAndConquer + -- HGeometry.Plane.LowerEnvelope.EpsApproximation HGeometry.VoronoiDiagram.ViaLowerEnvelope @@ -710,7 +701,6 @@ test-suite hspec ConvexHull.MelkmanSpec ClosestPair.ClosestPairSpec IntervalTreeSpec - LowerEnvelope.NaiveSpec LowerEnvelope.RegionsSpec Polygon.Convex.ConvexSpec Polygon.TriangulateSpec diff --git a/hgeometry/ipe/src/Ipe/IpeOut.hs b/hgeometry/ipe/src/Ipe/IpeOut.hs index e6340d5f9..9e5a42eab 100644 --- a/hgeometry/ipe/src/Ipe/IpeOut.hs +++ b/hgeometry/ipe/src/Ipe/IpeOut.hs @@ -33,6 +33,7 @@ import HGeometry.Foldable.Util import HGeometry.HalfLine import HGeometry.Intersection import HGeometry.Line +import HGeometry.Line.General import HGeometry.LineSegment import HGeometry.Number.Radical import HGeometry.Point @@ -158,6 +159,16 @@ instance (Fractional r, Ord r, Show r) => HasDefaultIpeOut (LinePV 2 r) where type DefaultIpeOut (LinePV 2 r) = Path defIO = ipeLine +instance (Fractional r, Ord r, Show r) => HasDefaultIpeOut (LineEQ r) where + type DefaultIpeOut (LineEQ r) = Path + defIO (LineEQ a b) = ipeLine $ fromLinearFunction a b + +instance (Fractional r, Ord r, Show r) => HasDefaultIpeOut (VerticalOrLineEQ r) where + type DefaultIpeOut (VerticalOrLineEQ r) = Path + defIO = \case + VerticalLineThrough x -> ipeLine $ verticalLine x + NonVertical l -> defIO l + instance (Fractional r, Ord r, Point_ point 2 r,Show r, Show point) => HasDefaultIpeOut (HalfLine point) where type DefaultIpeOut (HalfLine point) = Path defIO = ipeHalfLine diff --git a/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs b/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs index d4849a332..9402781e0 100644 --- a/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs +++ b/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs @@ -36,7 +36,7 @@ type UpperHull point = LowerEnvelope (NonVerticalHyperPlane 3 (NumType point) :+ -- | Computes the upper hull of a set of points in R^3 -- --- O(n^4) +-- \(O(n^4)\) upperHull :: ( Point_ point 3 r , Ord r, Fractional r , Foldable1 f, Functor f @@ -48,7 +48,7 @@ upperHull = lowerEnvelope . fmap (\p -> dualHyperPlane p :+ p) type Facet point = [point] --- | Outputs the facet sof the upper hull. +-- | Outputs the facets of the upper hull. facets :: UpperHull point -> [Facet point] facets = \case ParallelStrips _ -> error "facets: parallel strips; no bounded facets" diff --git a/hgeometry/src/HGeometry/Line/LowerEnvelope.hs b/hgeometry/src/HGeometry/Line/LowerEnvelope.hs index 94030d28a..044b4a9ef 100644 --- a/hgeometry/src/HGeometry/Line/LowerEnvelope.hs +++ b/hgeometry/src/HGeometry/Line/LowerEnvelope.hs @@ -42,6 +42,7 @@ import HGeometry.Sequence.Alternating -- | The lower envelope of a set of lines newtype LowerEnvelopeF f vertex line = LowerEnvelope (Alternating f vertex line) + deriving newtype (Functor, Bifunctor, Foldable) -- | A lower envelope, where the data structure is a vector. type LowerEnvelope = LowerEnvelopeF Vector.Vector @@ -61,11 +62,7 @@ deriving instance ( Eq line, Eq (f (vertex, line)) deriving instance ( Ord line, Ord (f (vertex, line)) ) => Ord (LowerEnvelopeF f vertex line) -instance Functor f => Functor (LowerEnvelopeF f r) where - fmap f = bimap id f -instance Functor f => Bifunctor (LowerEnvelopeF f) where - bimap f g (LowerEnvelope alt) = LowerEnvelope $ bimap f g alt -- -- | change the functor type of the Lower envelope diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope.hs index 1048240ef..45b8f3dc7 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope.hs @@ -9,13 +9,12 @@ -- -------------------------------------------------------------------------------- module HGeometry.Plane.LowerEnvelope - ( module HGeometry.Plane.LowerEnvelope.Naive - , module HGeometry.Plane.LowerEnvelope.Type - , module HGeometry.Plane.LowerEnvelope.AdjListForm - -- , module HGeometry.Plane.LowerEnvelope.Naive + ( module HGeometry.Plane.LowerEnvelope.Naive + , module HGeometry.Plane.LowerEnvelope.Type + , module HGeometry.Plane.LowerEnvelope.Connected ) where -import HGeometry.Plane.LowerEnvelope.AdjListForm +import HGeometry.Plane.LowerEnvelope.Connected import HGeometry.Plane.LowerEnvelope.Naive import HGeometry.Plane.LowerEnvelope.Type diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/AdjListForm.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/AdjListForm.hs deleted file mode 100644 index 2a1451a6a..000000000 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/AdjListForm.hs +++ /dev/null @@ -1,93 +0,0 @@ -{-# LANGUAGE UndecidableInstances #-} --------------------------------------------------------------------------------- --- | --- Module : HGeometry.Plane.LowerEnvelope.AdjListForm --- Copyright : (C) Frank Staals --- License : see the LICENSE file --- Maintainer : Frank Staals --- --- A Representation of the Lower envelope of planes in Adjacency-list --- form. --- --------------------------------------------------------------------------------- -module HGeometry.Plane.LowerEnvelope.AdjListForm - ( LowerEnvelope(..) - , LowerEnvelope'(LowerEnvelope) - , ParallelPlane - , theUnboundedVertex, boundedVertices - - , singleton - , fromVertexForm - - , BoundedVertexF(Vertex) - , location, definers, location2 - - , UnboundedVertex(UnboundedVertex) - , unboundedVertexId - , HasUnboundedEdges(..) - - , EdgeGeometry - , projectedEdgeGeometries, projectedEdgeGeometry - ) where - - --------------------------------------------------------------------------------- - -import Control.Lens -import qualified Data.Foldable as F -import Data.Foldable1 -import Data.List.NonEmpty (NonEmpty(..)) -import qualified Data.Set as Set -import HGeometry.HyperPlane.NonVertical -import HGeometry.Plane.LowerEnvelope.Connected -import qualified HGeometry.Plane.LowerEnvelope.VertexForm as VertexForm -import HGeometry.Properties -import HGeometry.Vector.NonEmpty.Util () - --------------------------------------------------------------------------------- --- * Data type defining a lower envelope - --- | The lower enevelope of planes in R^3. (Or rather, its minimization diagram) -data LowerEnvelope plane = - ParallelStrips !(Set.Set (ParallelPlane plane)) - | ConnectedEnvelope !(LowerEnvelope' plane) - -deriving instance (Show plane, Show (NumType plane)) => Show (LowerEnvelope plane) -deriving instance (Eq plane, Eq (NumType plane)) => Eq (LowerEnvelope plane) - -type instance NumType (LowerEnvelope plane) = NumType plane -type instance Dimension (LowerEnvelope plane) = 3 - --- | Just a newtype around plane, to be used to model parallel strips in the Lower envelope. -newtype ParallelPlane plane = - ParallelPlane plane deriving (Show,Eq) - -instance Wrapped (ParallelPlane plane) where - type Unwrapped (ParallelPlane plane) = plane - _Wrapped' = coerced - --- instance Rewrapped (ParallelPlane plane) plane - --------------------------------------------------------------------------------- - --- | Given a Lower envelope in vertex form, construct the AdjacencyList representation out --- of it. --- --- \(O(n\log n)\) -fromVertexForm :: forall f plane r. ( Plane_ plane r, Ord plane, Ord r, Fractional r - , Show plane, Show r - , Foldable1 f - ) - => f plane -> VertexForm.VertexForm plane -> LowerEnvelope plane -fromVertexForm _hs lEnv - | VertexForm.hasVertices lEnv = ConnectedEnvelope $ fromVertexForm' lEnv - | otherwise = ParallelStrips $ Set.fromDistinctAscList (F.toList strips) - where - strips :: NonEmpty (ParallelPlane plane) - strips = error "TODO" - -- coerce - -- $ divideAndConquer1With (mergeAndDiscardBy cmpPlanes) NonEmpty.singleton hs - --- withBisectors :: Plane_ plane r --- => f plane -> Alternating.Alternating Set.Set plane (IntersectionLine r) --- withBisectors = undefined diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/AtMostThree.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/AtMostThree.hs deleted file mode 100644 index aebefd7d8..000000000 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/AtMostThree.hs +++ /dev/null @@ -1,52 +0,0 @@ -module HGeometry.LowerEnvelope.AtMostThree - ( AtMostThree(..) - , commonElems - ) where - -import qualified Data.Set as Set -import HGeometry.Vector -import Witherable -import qualified Data.Foldable as F - --------------------------------------------------------------------------------- - --- | Lists of length at most three -data AtMostThree a = None - | One !a - | Two !a !a - | Three !a !a !a - deriving (Show,Eq,Ord,Functor,Foldable,Traversable) - -instance Filterable AtMostThree where - mapMaybe f = \case - None -> None - One x -> case f x of - Nothing -> None - Just x' -> One x' - Two x y -> case f x of - Nothing -> None - Just x' -> case f y of - Nothing -> One x' - Just y' -> Two x' y' - Three x y z -> case f x of - Nothing -> None - Just x' -> case f y of - Nothing -> One x' - Just y' -> case f z of - Nothing -> Two x' y' - Just z' -> Three x' y' z' - -instance Witherable AtMostThree - --- | Compute the common elements of two vectors of length 3 -commonElems :: Ord a - => Vector 3 a - -> Vector 3 a - -> AtMostThree a -commonElems u v = case Set.toAscList $ Set.intersection (Set.fromList $ F.toList u) - (Set.fromList $ F.toList v) of - [] -> None - [x] -> One x - [x,y] -> Two x y - [x,y,z] -> Three x y z - _ -> error "commonElems: absurd" diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs index 2214cb9ff..066469f8d 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs @@ -1,37 +1,24 @@ -{-# LANGUAGE UndecidableInstances #-} -------------------------------------------------------------------------------- -- | --- Module : HGeometry.LowerEnvelope.Connected +-- Module : HGeometry.Plane.LowerEnvelope.Connected -- Copyright : (C) Frank Staals -- License : see the LICENSE file -- Maintainer : Frank Staals -- --- A Representation of the Lower envelope of planes in Adjacency-list --- form. +-- A Representation of the Lower envelope of planes as a bunch of convex regions -- -------------------------------------------------------------------------------- module HGeometry.Plane.LowerEnvelope.Connected - ( LowerEnvelope'(LowerEnvelope) - , theUnboundedVertex, boundedVertices - - , singleton - , fromVertexForm' - - , BoundedVertexF(Vertex) - , location, definers, location2 - - , UnboundedVertex(UnboundedVertex) - , unboundedVertexId - , HasUnboundedEdges(..) - - , EdgeGeometry - , projectedEdgeGeometries, projectedEdgeGeometry + ( MinimizationDiagram + , asMap + , Region(..) + , CircularList + , module HGeometry.Plane.LowerEnvelope.Connected.Primitives + , module HGeometry.Plane.LowerEnvelope.Connected.Regions + , module HGeometry.Plane.LowerEnvelope.Connected.BruteForce ) where -import HGeometry.Plane.LowerEnvelope.Type -import HGeometry.Plane.LowerEnvelope.Connected.Type -import HGeometry.Plane.LowerEnvelope.Connected.FromVertexForm -import HGeometry.Vector.NonEmpty.Util () - - --------------------------------------------------------------------------------- +import HGeometry.Plane.LowerEnvelope.Connected.BruteForce +import HGeometry.Plane.LowerEnvelope.Connected.Primitives +import HGeometry.Plane.LowerEnvelope.Connected.Regions +import HGeometry.Plane.LowerEnvelope.Connected.Type diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs index 61fe67213..439c818ea 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/BruteForce.hs @@ -20,6 +20,7 @@ import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap import HGeometry.Plane.LowerEnvelope.Connected.Regions +import HGeometry.Plane.LowerEnvelope.Connected.Type import HGeometry.Point -------------------------------------------------------------------------------- diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/FromVertexForm.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/FromVertexForm.hs deleted file mode 100644 index 875446872..000000000 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/FromVertexForm.hs +++ /dev/null @@ -1,337 +0,0 @@ -{-# LANGUAGE UndecidableInstances #-} --------------------------------------------------------------------------------- --- | --- Module : HGeometry.Plane.LowerEnvelope.FromVertexForm --- Copyright : (C) Frank Staals --- License : see the LICENSE file --- Maintainer : Frank Staals --- --- A Representation of the Lower envelope of planes in Adjacency-list --- form. --- --------------------------------------------------------------------------------- -module HGeometry.Plane.LowerEnvelope.Connected.FromVertexForm - ( fromVertexForm' - ) where - -import Control.Lens -import qualified Data.Foldable as F -import Data.Foldable1 -import Data.Function (on) -import qualified Data.List as List -import Data.List.NonEmpty (NonEmpty(..)) -import qualified Data.List.NonEmpty as NonEmpty -import Data.Ord (comparing) -import qualified Data.Sequence as Seq -import qualified Data.Set as Set -import qualified Data.Vector as V -import Data.Vector.NonEmpty (NonEmptyVector) -import qualified Data.Vector.NonEmpty as NonEmptyV -import HGeometry.Combinatorial.Util -import HGeometry.Ext -import HGeometry.Foldable.Sort -import HGeometry.Foldable.Util -import HGeometry.HyperPlane.NonVertical -import HGeometry.Line -import HGeometry.Plane.LowerEnvelope.Connected.Type -import HGeometry.Plane.LowerEnvelope.Type -import qualified HGeometry.Plane.LowerEnvelope.VertexForm as VertexForm -import HGeometry.Point -import HGeometry.Properties -import HGeometry.Vector -import qualified HGeometry.Vector as Vec -import HGeometry.Vector.NonEmpty.Util () -import Witherable - -import Debug.Trace - --------------------------------------------------------------------------------- - --- | Given a Lower envelope in vertex form, construct the AdjacencyList representation out --- of it. --- --- pre: The set of vertices is non-empty --- --- \(O(n\log n)\) -fromVertexForm' :: forall plane r. (Plane_ plane r, Ord plane, Ord r, Fractional r - , Show plane, Show r - ) - => VertexForm.VertexForm plane -> LowerEnvelope' plane -fromVertexForm' lEnv = LowerEnvelope v0 boundedVs - where - v0 = UnboundedVertex unboundedEdges' - - -- TODO: we may want to sort the edges cyclically around every vertex. - - -- the bounded vertices. - boundedVs = Seq.zipWith ( \(_,v,_) outs -> v&incidentEdgesB .~ toSeq outs - ) boundedVs' boundedEsByOrigin - - toSeq = foldMap (Seq.singleton . snd) - - -- all edges leaving from bounded vertices, ordered by their VertexIds - boundedEsByOrigin = Seq.fromList . groupOnCheap fst . foldMap F.toList - $ allEdges - - -- computes all outgoing edges of all bounded vertices, they are grouped by face - allEdges :: [FaceEdges plane] - allEdges = traceShowWith ("allEdges",) . - fmap (faceToEdges . sortAlongBoundary) - . groupOnCheap definingPlane - $ foldMap (^._3) boundedVs' - - -- | construct the unbounded edges, by simply filtering and flipping the - -- edges towards infinity. - unboundedEdges' :: Seq.Seq (LEEdge plane) - unboundedEdges' = mapMaybe (\(u, e) -> if e^.destination == unboundedVertexId - then Just (flipEdge u e) else Nothing - -- flip every unbounded edge - ) - $ foldMap fromFoldable allEdges - - - -- | A sequence of bounded vertices, ordered by vertexID, together with a bunch of - -- copies; one for each face it will appear on - -- - -- TODO: i guess that in degenerate sitations it could be that there are copies that - -- wont actually proudce an edge (i.e. if there is some definer/plane containing the - -- vertex, that does not appear on the lower envelope). - -- - boundedVs' :: Seq.Seq ( VertexID, BoundedVertex plane, [IntermediateVertex plane]) - boundedVs' = ifoldMapOf (indexing (VertexForm.vertices'.withIndex)) mkVtx lEnv - - mkVtx i (v,defs) = let j = i+1 - in Seq.singleton ( j - , Vertex v defs mempty - , [ Vtx h j v defs | h <- F.toList defs ] - -- TODO: conceivably, we could delete the h from the - -- defs here already? instead of in the various delete h - -- cases - ) - -{- --- | Given a bunch of halflines that all share their starting point v, --- sort them cyclically around the starting point v. -sortCyclic :: (Foldable f, Plane_ plane r, Ord r, Fractional r) - => f (LEEdge plane) -> Seq.Seq (LEEdge plane) -sortCyclic = fromFoldable . sortBy @V.Vector cmpAround - where - cmpAround e1 e2 = ccwCmpAround origin (dir e1) (dir e2) - dir e = Point . negated $ edgeIntersectionLine e ^. direction - -- these direction vectors all point "inward", we negate them to make sure they are - -- pointing outward (towards +infty again). --} - - --- | For each bounded vertex (represented by their ID) the outgoing halfedge. -type FaceEdges plane = NonEmpty (VertexID, LEEdge plane) - - --- | Helper type for edgedefs -data EdgeDefs plane = EdgeDefs { common :: plane - , uNeigh :: plane - , vNeigh :: plane - } deriving (Show,Eq) - --- | Given a plane h, and vertices u (with its definers), and v (with its definers) that --- define an edge of h, computes: --- --- - plane h' that is on the other side of the edge from u to v, --- - the plane hu incident only to u that is adjacent to h, and --- - the plane hv incident only to v that is adjacent to h. -extractEdgeDefs :: (Ord plane - , Show plane, Show r - ) - => plane - -> Point 3 r -> VertexForm.Definers plane - -> Point 3 r -> VertexForm.Definers plane - -> Maybe (EdgeDefs plane) -extractEdgeDefs h u uDefs v vDefs - -- | traceShow ("extractEdgeDefs",h,u,uDefs,v,vDefs) False = undefined - -- | otherwise - = case traceShowWith ("commons",) commons of - [] -> Nothing - [h'] -> Just $ EdgeDefs h' hu hv - _ -> error "extractEdgeDefs: unhandled degeneracy. u and v have >2 planes in common." - where - commons = F.toList $ Set.delete h --- $ traceShowWith ("udefs intersect vDefs",) - (uDefs `Set.intersection` vDefs) - uOnlies = F.toList $ uDefs Set.\\ vDefs - vOnlies = F.toList $ vDefs Set.\\ uDefs - - hu = from' u uOnlies - hv = from' v vOnlies - - from' x hss - |traceShow ("from'",x,hss) False = undefined - | otherwise - = case hss of - [] -> error "extractEdgeDefs: absurd, too few definers" - [h'] -> h' - _hs -> error $ "extractEdgeDefs: unhandled degeneracy. More than 3 planes at a vertex. " - <> show ("extractEdgeDefs",h,u,v,uOnlies,vOnlies) - -- TODO we should either the neighbor of h in the order around the given - -- vertex here. - --------------------------------------------------------------------------------- - --- | Intermediate representation of a vertex. -data IntermediateVertex plane = Vtx { definingPlane :: !plane - , ivId :: {-# UNPACK #-} !VertexID - , ivLoc :: !(Point 3 (NumType plane)) - , _ivDefs :: VertexForm.Definers plane - } - --------------------------------------------------------------------------------- - --- | Sort the vertices of a (convex) face in counter clockwise order around its --- boundary. --- --- running time: \(O(n \log n)\) -sortAlongBoundary :: forall plane r. (Plane_ plane r, Ord r, Num r - , Show plane, Show r - ) - => NonEmptyVector (IntermediateVertex plane) - -> NonEmptyVector (IntermediateVertex plane) -sortAlongBoundary face = case mv0 of - Nothing -> face -- already sorted, since there is only one vertex - Just (Vtx _ _ v0 _) -> NonEmptyV.unsafeFromVector $ sortBy (cmpAround $ v1 .-. v0) face - -- TODO: I guess that in degenerate situations, there now may be consecutive - -- points (vertices) at the same location. We should group those and get rid of them? - where - -- we find the leftmost vertex v1 of the face (pick the lexicographically smallest - -- when there are multiple), and then compute its predecessor v0 in the order along - -- the boundary of the face (i.e. we find the point v for which the line through v1 - -- and v has maximum slope). - - -- To get the list of vertices in CCW order along the face, we then sort the vertices - -- around v1 - Vtx _ iv1 v1 _ = minimumBy (comparing ivLoc) face - v1' = projectPoint v1 - -- its predecessor v0, if it exists - mv0 :: Maybe (IntermediateVertex plane) - mv0 = foldr findPred Nothing face - findPred v acc - | ivId v == iv1 = acc -- skip v1 itself - | otherwise = Just $ case acc of - Nothing -> v -- anything is better than nothing - Just u -> maxBy cmpSlope' u v - - cmpSlope' :: IntermediateVertex plane -> IntermediateVertex plane -> Ordering - cmpSlope' u v = case ccw v1' (projectPoint $ ivLoc u) (projectPoint $ ivLoc v) of - CCW -> GT - CW -> LT - CoLinear -> EQ - - -- | sort the vertices around the direction wrt. (the downward projection of) w - cmpAround :: Vector 3 r -> IntermediateVertex plane -> IntermediateVertex plane -> Ordering - cmpAround w u v = - ccwCmpAroundWith (Vec.prefix w) v1' u' v' <> cmpByDistanceTo v1' u' v' - where - u' = projectPoint $ ivLoc u - v' = projectPoint $ ivLoc v - --- | given the vertices of the face in CCW order, compute all its edges of the face. -faceToEdges :: forall plane r. (Plane_ plane r, Ord r, Fractional r, Ord plane - , Show plane, Show r - ) - => NonEmptyVector (IntermediateVertex plane) -> FaceEdges plane -faceToEdges faceV = case toNonEmpty faceV of - u :| [] -> oneVertex u - u :| [v] -> twoVertices u v - _ -> toNonEmpty $ manyVertices - where - manyVertices = NonEmptyV.zipWith3 mkEdge (shiftR (-1) faceV) faceV (shiftR 1 faceV) - - mkEdge (Vtx _ _ up uDefs) (Vtx h vIdx vp vDefs) (Vtx _ wIdx wp wDefs) = - case extractEdgeDefs h vp vDefs wp wDefs of - Nothing -> -- vw are separated by the unbounded vertex - case extractEdgeDefs h up uDefs vp vDefs of - Nothing -> - error "mkEdge: absurd, u and v don't share a plane!?" - Just (EdgeDefs _hUV _hu hv) -> (vIdx, Edge unboundedVertexId h hv) - Just (EdgeDefs hVW _ _) -> (vIdx, Edge wIdx h hVW) - -- the main idea is that every pair of subsequent bounded vertices necessarily has - -- exactly one defining plane, except for when two bounded vertices are separted by - -- the vertex at infinity. - - --- | compute the face edges of the face of h, when v is the only vertex incident to h. -oneVertex :: (Plane_ plane r, Ord r, Fractional r, Ord plane - , Show plane, Show r - ) - => IntermediateVertex plane -> FaceEdges plane -oneVertex (Vtx h i v defs) = case List.sort outgoingEdges of - [ Left _hPred, Right hSucc ] -> NonEmpty.singleton (i, Edge unboundedVertexId h hSucc) - _ -> error "oneVertex. absurd. Other than a single Left, and Right" - where - otherPlanes = Set.delete h defs - -- | Tries to compute all outgoing edges involving h. Since every plane h appears - -- exactly once in the order around v, this should produce exactly two edges; one - -- left and one right. - outgoingEdges = foldMap (\h' -> let rest = toNonEmpty' $ Set.delete h' otherPlanes - in case outgoingUnboundedEdge v (Two h h') rest of - Just e -> [ tagOtherPlane $ e^.extra ] - Nothing -> [] - ) otherPlanes - tagOtherPlane (EdgeDefiners hl hr) = if h == hl then Right hr else Left hl - toNonEmpty' s = case NonEmpty.nonEmpty $ F.toList s of - Just xs -> xs - _ -> error "oneVertex. Absurd, there should be at least 3 definers" - --- | Compute the edges of the face incident to h, when there are only two vertices --- incident to that face. --- --- more or less a special case of the manyVertices scenario, in which we don't know --- the if edge should be oriented from u to v or from v to u. -twoVertices :: (Plane_ plane r, Ord r, Fractional r, Ord plane - , Show plane, Show r - ) - => IntermediateVertex plane -> IntermediateVertex plane -> FaceEdges plane -twoVertices (Vtx h ui up uDefs) (Vtx _ vi vp vDefs) = - case extractEdgeDefs h up uDefs vp vDefs >>= withUVOrder of - Nothing -> error "twoVertices. absurd, h and h' don't define an edge!?" - Just (EdgeDefs hr hu hv, uBeforeV) -- -> traceShowWith ("twoVertices",u,v,uBeforeV,) $ - | uBeforeV -> NonEmpty.fromList [ (vi, Edge ui h hr) - , (ui, Edge unboundedVertexId h hu) - ] - -- the edge must be oriented from v to u so that h is on the left - | otherwise -> NonEmpty.fromList [ (ui, Edge vi h hr) - , (vi, Edge unboundedVertexId h hv) - ] - -- the edge must be oriented from u to v - where - -- determine if u lies before v in the order of u and v along the intersection line of - -- h and hr. - withUVOrder e@(EdgeDefs hr _ _) = - ( \l -> let m = perpendicularTo l &anchorPoint .~ projectPoint up - in (e, projectPoint vp `onSide` m /= LeftSide) - -- if v lies on the left it lies further along commonLine. So u lies before - -- v if v does not lie in the left haflpalne - ) <$> intersectionLine' h hr - --------------------------------------------------------------------------------- --- * Convenience functions - --- | Group consecutive elements into non-empty groups. --- --- --- running time: \(O(n\log n)\) --- --- >>> groupOnCheap fst [ (1,"foo"), (2,"bar"), (2,"blaa"), (1,"boeez"), (1,"blap"), (4,"blax"), (4,"bleh"), (100,"floep")] --- [[(1,"foo"),(1,"boeez"),(1,"blap")],[(2,"bar"),(2,"blaa")],[(4,"blax"),(4,"bleh")],[(100,"floep")]] -groupOnCheap :: (Foldable f, Ord b) - => (a -> b) -> f a -> [NonEmptyV.NonEmptyVector a] -groupOnCheap f = fmap NonEmptyV.unsafeFromVector . V.groupBy ((==) `on` f) . sortOnCheap f - --- | returns the maximum using the given comparison function -maxBy :: (t -> t -> Ordering) -> t -> t -> t -maxBy cmp a b = case cmp a b of - LT -> b - _ -> a - - --- | shift the vector by d items to the right -shiftR :: Int -> NonEmptyVector v -> NonEmptyVector v -shiftR d v = let n = length v - in NonEmptyV.generate1 n $ \i -> v NonEmptyV.! ((i+n+d) `mod` n) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs index b0a76d930..52dc05062 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs @@ -21,7 +21,7 @@ import Data.Semigroup (First(..)) import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap -import HGeometry.Plane.LowerEnvelope.Connected.Regions +import HGeometry.Plane.LowerEnvelope.Connected.Type import HGeometry.Point import HGeometry.Vector @@ -45,8 +45,7 @@ instance (Ord r, Num r) => Ord (E r) where -- represented by its point, it stores a list of its outgoing edges, and some data. toPlaneGraph :: (Plane_ plane r, Num r, Ord r) => MinimizationDiagram r plane -> PlaneGraph (Point 2 r) (First r) (E r) -toPlaneGraph = mapWithKeyMerge toTriangulatedGr - +toPlaneGraph = mapWithKeyMerge toTriangulatedGr . asMap toTriangulatedGr :: (Plane_ plane r, Num r, Ord r) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Regions.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Regions.hs index cae262f5b..3a4585222 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Regions.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Regions.hs @@ -10,10 +10,7 @@ -- -------------------------------------------------------------------------------- module HGeometry.Plane.LowerEnvelope.Connected.Regions - ( MinimizationDiagram - , Region(..) - , CircularList - , fromVertexForm + ( fromVertexForm , Definers , fromCCWList @@ -37,106 +34,13 @@ import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical import HGeometry.Intersection import HGeometry.Line +import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap +import HGeometry.Plane.LowerEnvelope.Connected.Primitives +import HGeometry.Plane.LowerEnvelope.Connected.Type import HGeometry.Point import HGeometry.Properties import HGeometry.Vector -import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap -import HGeometry.Plane.LowerEnvelope.Connected.Primitives - --------------------------------------------------------------------------------- --- * The Minimization Diagram, i.e. the type that we use to represent our --- lower envelope - --- | A minimization daigram just maps every plane on the lower envelope to the region --- above which it is minimal. Every plane has at most one such a region. -type MinimizationDiagram r plane = Map plane (Region r (Point 2 r)) - -type instance NumType (MinimizationDiagram r plane) = r -type instance Dimension (MinimizationDiagram r plane) = 2 --- TODO: this is a bit of hack; maybe use a newtype instead. We only need these instances --- to actually print the thing in ipe. - --- | A region in the minimization diagram. The boundary is given in CCW order; i.e. the --- region is to the left of the boundary. -data Region r point = Bounded (CircularList point) - | Unbounded (Vector 2 r) - -- ^ vector indicating the direction of the unbounded edge - -- incident to the first vertex. Note that this vector - -- thus points INTO vertex v. - (NonEmpty point) - -- ^ the vertices in CCW order, - (Vector 2 r) - -- ^ the vector indicating the direction of the unbounded - -- edge incident to the last vertex. The vector points - -- away from the vertex (i.e. towards +infty). - deriving stock (Show,Eq,Functor,Foldable,Traversable) - -type instance NumType (Region r point) = r -type instance Dimension (Region r point) = Dimension point - --- | bounded regions are really circular lists, but we just represent them as lists for --- now. -type CircularList a = [a] - - --------------------------------------------------------------------------------- --- * The planes defining a vertex - --- | The vertices of a lower envelope is just a Map with every vertex its definers, --- i.e. the planes that define the vertex in CCW order around it. -type VertexForm r plane = Map (Point 3 r) (Definers plane) - ----------------------------------------- --- * The planes defining a vertex - --- | in CCW order, starting with the plane that is minimal at the vertical up direction --- from their common vertex. -newtype Definers plane = Definers [plane] - deriving stock (Show,Eq,Ord) - deriving newtype (Functor,Foldable) - --- | Given the planes in order, starting with the one that is closest in the up direction, --- construct the Definers. -fromCCWList :: [plane] -> Definers plane -fromCCWList = Definers - --- | Smart constructor for creating the definers of three planes -definers :: forall plane r.(Plane_ plane r, Ord r, Fractional r) - => Three plane -> Maybe (Point 3 r, Definers plane) -definers (Three h1@(Plane_ a1 b1 c1) h2 h3) = - do l12 <- intersectionLine h1 h2 - l13 <- intersectionLine h1 h3 - intersect l12 l13 >>= \case - Line_x_Line_Line _ -> Nothing - Line_x_Line_Point (Point2 x y) -> Just ( Point3 x y (a1 * x + b1* y + c1) - , Definers [hMin, hTwo, hThree] - ) - where - (hMin,h,h') = extractMinOn (evalAt $ Point2 x (y+1)) h1 h2 h3 - -- we compute the plane hMin that is cheapest directly above the vertex h and - -- h' are the other two planes. That Means hMin is the first definer (in the - -- CCW order). What remains is to determine the order in which the remaining - -- planes h and h' appear in the order. - (hTwo,hThree) = case ccwCmpAroundWith (Vector2 0 1) (origin :: Point 2 r) - (Point vMinH) (Point vMinH') of - LT -> (h,h') - EQ -> error "definers: weird degeneracy?" - GT -> (h',h) - - vMinH = fromMaybe err $ intersectionVector h hMin - vMinH' = fromMaybe err $ intersectionVector h' hMin - err = error "definers: absurd" - --- | given three elements, returns the minimum element in the first argument and the --- remaining two elements in the second and third argument (in arbitrary order). -extractMinOn :: Ord a => (c -> a) -> c -> c -> c -> (c, c, c) -extractMinOn f a b c = let (m,ab) = min' a b - (mi,c') = min' m c - in (mi, ab, c') - where - min' x y - | f x <= f y = (x,y) - | otherwise = (y,x) +import HGeometry.Plane.LowerEnvelope.Connected.VertexForm ---------------------------------------- @@ -231,7 +135,8 @@ mergeDefiners (Point3 x y _) = foldr (insertPlane $ Point2 x y) -- \(O(h\log h)\) assuming that the input is degenerate. fromVertexForm :: (Plane_ plane r, Ord plane, Ord r, Fractional r) => VertexForm r plane -> MinimizationDiagram r plane -fromVertexForm = Map.mapWithKey sortAroundBoundary . mapWithKeyMerge (\v defs -> +fromVertexForm = MinimizationDiagram + . Map.mapWithKey sortAroundBoundary . mapWithKeyMerge (\v defs -> Map.fromList [ (h, Set.singleton (v,defs)) | h <- F.toList defs ]) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs index f748a8891..611a91994 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs @@ -11,439 +11,93 @@ -- -------------------------------------------------------------------------------- module HGeometry.Plane.LowerEnvelope.Connected.Type - ( LowerEnvelope'(LowerEnvelope) - , theUnboundedVertex, boundedVertices - , traverseLowerEnvelope + ( MinimizationDiagram(..) + , asMap + , Region(..) + , CircularList - , singleton - , BoundedVertex - , BoundedVertexF(Vertex) - , location, definers, location2 + -- LowerEnvelope'(LowerEnvelope) + -- , theUnboundedVertex, boundedVertices + -- , traverseLowerEnvelope - , UnboundedVertex(UnboundedVertex) - , unboundedVertexId - , HasUnboundedEdges(..) + -- , singleton - , EdgeGeometry - , projectedEdgeGeometries, projectedEdgeGeometry + -- , BoundedVertex + -- , BoundedVertexF(Vertex) + -- , location, definers, location2 + -- , UnboundedVertex(UnboundedVertex) + -- , unboundedVertexId + -- , HasUnboundedEdges(..) - , outgoingUnboundedEdge - , edgeIntersectionLine - , intersectionLine' - , intersectionLineWithDefiners - , EdgeDefiners(..) + -- , EdgeGeometry + -- , projectedEdgeGeometries, projectedEdgeGeometry + + + -- , outgoingUnboundedEdge + -- , edgeIntersectionLine + -- , intersectionLine' + -- , intersectionLineWithDefiners + -- , EdgeDefiners(..) ) where -import Control.Applicative -import Control.Lens -import qualified Data.Foldable as F -import Data.Foldable1 -import Data.Functor.Apply (Apply, (<.*>)) -import qualified Data.Functor.Apply as Apply -import qualified Data.List.NonEmpty as NonEmpty -import Data.Semigroup.Traversable(traverse1Maybe) -import qualified Data.Sequence as Seq -import qualified Data.Set as Set -import qualified Data.Vector as V -import HGeometry.Box -import HGeometry.Combinatorial.Util -import HGeometry.Ext -import HGeometry.Foldable.Sort -import HGeometry.Foldable.Util -import HGeometry.HalfLine -import HGeometry.HyperPlane.Class -import HGeometry.HyperPlane.NonVertical -import HGeometry.Line -import HGeometry.Line.General -import HGeometry.LineSegment -import HGeometry.Plane.LowerEnvelope.Type -import HGeometry.Plane.LowerEnvelope.VertexForm (intersectionLine) -import qualified HGeometry.Plane.LowerEnvelope.VertexForm as VertexForm +import Control.Subcategory.Functor +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Map (Map) +import qualified Data.Map as Map import HGeometry.Point import HGeometry.Properties import HGeometry.Vector import HGeometry.Vector.NonEmpty.Util () -import Hiraffe.Graph -import Witherable - - --------------------------------------------------------------------------------- --- * Data type defining a connected lower envelope - --- | A connected lower envelope in adjacencylist form. --- --- invariant: there is always at least one bounded vertex. -data LowerEnvelope' plane = - LowerEnvelope !(UnboundedVertex plane) (Seq.Seq (BoundedVertex plane)) - -deriving instance (Show plane, Show (NumType plane)) => Show (LowerEnvelope' plane) -deriving instance (Eq plane, Eq (NumType plane)) => Eq (LowerEnvelope' plane) - -type instance NumType (LowerEnvelope' plane) = NumType plane -type instance Dimension (LowerEnvelope' plane) = 3 - -instance (Ord (NumType plane), Num (NumType plane)) => IsBoxable (LowerEnvelope' plane) where - -- ^ the bounding box contains all bounded vertices - boundingBox env = boundingBox . NonEmpty.fromList $ env^..boundedVertices.traverse.location - -- the fromList is safe since there is alwasy at least one vertex - --- | Traversal of the planes in the lower envelope. Since we need an Ord constraint we --- can't make LowerEnvelope' an instance of Traversable. --- --- Be aware that this may destroy some of the invariants. So use this function with care. -traverseLowerEnvelope :: ( Applicative f, NumType plane ~ NumType plane' - , Ord plane' - ) - => (plane -> f plane') - -> LowerEnvelope' plane - -> f (LowerEnvelope' plane') -traverseLowerEnvelope f (LowerEnvelope v0 vs) = - LowerEnvelope <$> traverse f v0 <*> traverse (traverseBoundedV f) vs - --- | lens to access the unbounded vertex -theUnboundedVertex :: Lens' (LowerEnvelope' plane) (UnboundedVertex plane) -theUnboundedVertex = lens (\(LowerEnvelope v _) -> v) - (\(LowerEnvelope _ vs) v -> LowerEnvelope v vs) - --- | Lens to access the sequence of bounded vertices. -boundedVertices :: Lens' (LowerEnvelope' plane) (Seq.Seq (BoundedVertex plane)) -boundedVertices = lens (\(LowerEnvelope _ vs) -> vs) - (\(LowerEnvelope u _ ) vs -> LowerEnvelope u vs) - --------------------------------------------------------------------------------- - --- | Given a single Vertex, construct a LowerEnvelope' out of it. -singleton :: (Plane_ plane r, Ord r, Fractional r, Ord plane) - => VertexForm.LEVertex plane -> LowerEnvelope' plane -singleton v = LowerEnvelope v0 (Seq.singleton v') - where - i = 1 -- the vertexID we are using for this vertex. - v' = fromLEVertex v - v0 = UnboundedVertex $ flipEdge i <$> view incidentEdgesB v' - -- do we need to reverse the sequenceo f edges? - --------------------------------------------------------------------------------- - --- | The unbounded vertex. -newtype UnboundedVertex plane = UnboundedVertex { _incidentEdgesU :: Seq.Seq (LEEdge plane) } - deriving (Show,Eq,Functor,Foldable,Traversable) - --- | access the incidentEdges of an unbounded vertex -incidentEdgesU :: Iso (UnboundedVertex plane) (UnboundedVertex plane') - (Seq.Seq (LEEdge plane)) (Seq.Seq (LEEdge plane')) -incidentEdgesU = coerced - --- | The vertexID of the unbounded vertex -unboundedVertexId :: VertexID -unboundedVertexId = 0 - --------------------------------------------------------------------------------- - --- | Vertices in of the lower envelope in adjacencylist form. -type BoundedVertex = BoundedVertexF Seq.Seq - --- | Given a vertex in vertex form, construct A BoundedVertex from --- it. Note that this essentially assumes the vertex is connected to --- the unbounded vertex with all its outgoing edgse. --- --- \(O(k \log k)\), where \(k\) is the number of definers. -fromLEVertex :: (Plane_ plane r, Ord r, Fractional r, Ord plane) - => VertexForm.LEVertex plane - -> BoundedVertex plane -fromLEVertex (VertexForm.LEVertex v defs) = Vertex v defs es - where - es = (\(_ :+ EdgeDefiners hl hr) -> Edge unboundedVertexId hl hr) <$> sortAroundStart es' - es' = mapMaybe (\t@(Two h1 h2) -> let defs' = toNonEmpty' $ Set.delete h1 $ Set.delete h2 defs - in outgoingUnboundedEdge v t defs' - ) $ uniquePairs defs - toNonEmpty' s = case NonEmpty.nonEmpty $ F.toList s of - Just xs -> xs - _ -> error "fromLEVertex: absurd, there should be at least 3 definers" - - --- | Given a bunch of halflines that all share their starting point v, --- sort them cyclically around the starting point v. -sortAroundStart :: (Foldable f, Ord r, Num r) - => f (HalfLine (Point 2 r) :+ e) -> Seq.Seq (HalfLine (Point 2 r) :+ e) -sortAroundStart = fromFoldable . sortBy @V.Vector compareAroundStart - where - compareAroundStart (HalfLine v d :+ _) - (HalfLine _ d' :+ _) = ccwCmpAround v (v .+^ d) (v .+^ d') - - - --- | Given a location of a vertex v, a pair of planes h1,h2 and the --- remaining defining planes of v, computes the outgoing half-line --- from v on which h1,h2 are the lowest (if such an halfline exists). -outgoingUnboundedEdge :: ( Plane_ plane r, Ord r, Fractional r - , Foldable1 f - ) - => Point 3 r -- ^ the location of the vertex v - -> Two plane -- ^ the pair of planes for which to compute - -- the halfine - -> f plane -- ^ the other planes intersecting at v - -> Maybe (HalfLine (Point 2 r) :+ EdgeDefiners plane) -outgoingUnboundedEdge v (Two h1 h2) h3s = - intersectionLineWithDefiners h1 h2 >>= toHalfLineFrom (projectPoint v) h3s - -- todo, if there are more planes, I guess we should check if the hl is not dominated by the other - -- planes either. - --- | Given : --- --- v : the projected location of the vertex --- hs : the remaining planes defining v (typically just one plane h3) --- l : the (projection of the) line l in which planes h1 and h2 intersect (containing v) --- --- we compute the half-line eminating from v in which h1 and h2 define --- an edge incident to v. -toHalfLineFrom :: (Plane_ plane r, Foldable1 f, Fractional r, Ord r) - => Point 2 r -- ^ vertex v - -> f plane -- ^ the remaining plane(s) hs - -> LinePV 2 r :+ EdgeDefiners plane -- ^ the line l - -> Maybe (HalfLine (Point 2 r) :+ EdgeDefiners plane) -toHalfLineFrom v hs ((LinePV _ w) :+ defs@(EdgeDefiners h1 h2)) = - validate w defs <|> validate ((-1) *^ w) (EdgeDefiners h2 h1) - -- We try both directions. Note that if we reverse the direction - -- of the line, what the left/right plane is changes. - -- - -- If neither direction works, then h1,h2 do not define a good - -- direction. This should happen only when there are more than 3 - -- planes intersecting in v, i.e. in degernate sitatuions - where - -- | test if direction d is a good direction, i.e. if towards direction d - -- h1 (and thus also h2) is actually lower than all remaining defining planes. - validate d defs' = let zVal = evalAt (v .+^ d) - in if all (\h3 -> zVal h1 < zVal h3) hs - then Just (HalfLine v d :+ defs') else Nothing - - - - --- | The planes left (above) and right (below) their intersection --- line. -data EdgeDefiners plane = EdgeDefiners { _leftPlane :: plane -- ^ above plane - , _rightPlane :: plane -- ^ below plane - } - --- | Computes the line in which the two planes intersect, and labels --- the halfplanes with the lowest plane. --- --- The returned line will have h to its left and h' to its right. -intersectionLineWithDefiners :: ( Plane_ plane r, Ord r, Fractional r) - => plane -> plane -> Maybe (LinePV 2 r :+ EdgeDefiners plane) -intersectionLineWithDefiners h h' = (:+ EdgeDefiners h h') <$> intersectionLine' h h' - - --- | Computes the line in which the two planes intersect. The returned line will have h to --- its left and h' to its right. --- --- -intersectionLine' :: ( Plane_ plane r, Ord r, Fractional r) - => plane -> plane -> Maybe (LinePV 2 r) -intersectionLine' h h' = intersectionLine h h' <&> \case - VerticalLineThrough x -> reorient (LinePV (Point2 x 0) (Vector2 0 1)) (Point2 (x-1) 0) - NonVertical l -> let l'@(LinePV p _) = fromLineEQ l - in reorient l' (p&yCoord %~ (+1)) - where - -- make sure h is to the left of the line - reorient l q = let f = evalAt q - in if f h <= f h' then l else l&direction %~ negated - fromLineEQ (LineEQ a b) = fromLinearFunction a b - --- | Computes the supportingline of the (downward projection of) the edge. -edgeIntersectionLine :: ( Plane_ plane r, Ord r, Fractional r) - => LEEdge plane -> LinePV 2 r -edgeIntersectionLine e = case intersectionLine' (e^.leftPlane) (e^.rightPlane) of - Just l -> l - Nothing -> error "edgeIntersectionLine: absurd. no intersection line !?" - - - --------------------------------------------------------------------------------- - --- | Types that have an IncidentEdges' field. -class HasIncidentEdges t where - -- | Lens to access the incident edges field. - incidentEdges' :: Lens' (t plane) (Seq.Seq (LEEdge plane)) - -instance HasIncidentEdges UnboundedVertex where - incidentEdges' = incidentEdgesU - -instance HasIncidentEdges BoundedVertex where - incidentEdges' = lens (view incidentEdgesB) (\(Vertex p d _) es -> Vertex p d es) - --- instance HasIncidentEdges (Vertex (LowerEnvelope' plane)) plane where --- incidentEdges' = undefined -- pick either incidentEdges on the left or on the right thing. - -------------------------------------------------------------------------------- --- * The Instances for HasVertices, HasDarts, HasFaces etc - -instance HasVertices' (LowerEnvelope' plane) where - type Vertex (LowerEnvelope' plane) = Either (UnboundedVertex plane) (BoundedVertex plane) - type VertexIx (LowerEnvelope' plane) = VertexID - - -- | note, trying to assign the unbounded vertex to something with index /= - -- unboundedVertexId is an error. - vertexAt :: VertexID - -> IndexedTraversal' VertexID (LowerEnvelope' plane) - (Vertex (LowerEnvelope' plane)) - vertexAt i - | i == unboundedVertexId = conjoined trav0 (itrav0 .indexed) - | otherwise = conjoined travVs (itravVs.indexed) - where - trav0 f (LowerEnvelope v0 vs) = flip LowerEnvelope vs . unLeft <$> f (Left v0) - - itrav0 :: Applicative f - => (VertexID -> Vertex (LowerEnvelope' plane) - -> f (Vertex (LowerEnvelope' plane))) - -> LowerEnvelope' plane -> f (LowerEnvelope' plane) - itrav0 f (LowerEnvelope v0 vs) = flip LowerEnvelope vs . unLeft <$> f 0 (Left v0) - - travVs :: Applicative f - => (Vertex (LowerEnvelope' plane) - -> f (Vertex (LowerEnvelope' plane))) - -> LowerEnvelope' plane -> f (LowerEnvelope' plane) - travVs f (LowerEnvelope v0 vs) = - LowerEnvelope v0 <$> (vs&ix (i-1) %%~ fmap unRight . f . Right) - itravVs f (LowerEnvelope v0 vs) = - LowerEnvelope v0 <$> (vs&ix (i-1) %%~ fmap unRight . f i . Right) - - unLeft = either id (error "LowerEnvelope'.vertexAt: trying to convert v_infty into a normal vertex is not allowed") - unRight = either (error "LowerEnvelope'.vertexAt: trying to convert a regular bounded vertex into v_infty is not allowed") id - {-# INLINE vertexAt #-} - - -instance HasVertices (LowerEnvelope' plane) (LowerEnvelope' plane') where - - vertices = conjoined traverse' (itraverse' . indexed) - where - traverse' :: (Apply f) - => (Vertex (LowerEnvelope' plane) -> f (Vertex (LowerEnvelope' plane'))) - -> LowerEnvelope' plane -> f (LowerEnvelope' plane') - traverse' f (LowerEnvelope v0 vs) = - LowerEnvelope <$> (fmap unLeft . f .Left) v0 - <.*> traverse1Maybe (fmap unRight . f . Right) vs +-- * The Minimization Diagram, i.e. the type that we use to represent our +-- lower envelope - unLeft = either id (error "LowerEnvelope'.vertices: trying to convert v_infty into a normal vertex is not allowed") - unRight = either (error "LowerEnvelope'.vertices: trying to convert a regular bounded vertex into v_infty is not allowed") id +-- | A minimization daigram just maps every plane on the lower envelope to the region +-- above which it is minimal. Every plane has at most one such a region. +newtype MinimizationDiagram r plane = MinimizationDiagram (Map plane (Region r (Point 2 r))) + deriving stock (Show,Eq) - itraverse' :: (Apply f) - => (VertexID -> Vertex (LowerEnvelope' plane) -> f (Vertex (LowerEnvelope' plane'))) - -> LowerEnvelope' plane -> f (LowerEnvelope' plane') - itraverse' f (LowerEnvelope v0 vs) = - LowerEnvelope <$> (fmap unLeft . f unboundedVertexId .Left) v0 - <.*> itraverse1Maybe (\i -> fmap unRight . f (i+1) . Right) vs - {-# INLINE vertices #-} - --- | Indexed version of 'traverse1Maybe'. -itraverse1Maybe :: (TraversableWithIndex i t, Apply f) - => (i -> a -> f b) -> t a -> Apply.MaybeApply f (t b) -itraverse1Maybe f = itraverse (\i -> Apply.MaybeApply . Left . f i) - ----------------------------------------- - -instance HasDarts' (LowerEnvelope' plane) where - type Dart (LowerEnvelope' plane) = LEEdge plane - type DartIx (LowerEnvelope' plane) = ( VertexIx (LowerEnvelope' plane) - , VertexIx (LowerEnvelope' plane) - ) - dartAt (u,v) = vertexAt u <.> beside l l - where - l :: HasIncidentEdges t => IndexedTraversal' VertexID (t plane) (LEEdge plane) - l = incidentEdges' .> first' v (view destination) - {-# INLINE dartAt #-} - --- | Helper function to access the an edge -first' :: Eq i => i -> (a -> i) - -> IndexedTraversal' i (Seq.Seq a) a -first' i getIdx f s = case Seq.findIndexL ((== i) . getIdx) s of - Nothing -> pure s - Just j -> s&ix j %%~ indexed f i - - -instance HasDarts (LowerEnvelope' plane) (LowerEnvelope' plane) where - darts = itraverse' . indexed - -- the traverse' variant that does not use the indices does not really look all that - -- useful. So I didn't bother writing the faster version. - where - itraverse' :: (Applicative f) - => (DartIx (LowerEnvelope' plane) -> LEEdge plane -> f (LEEdge plane)) - -> LowerEnvelope' plane -> f (LowerEnvelope' plane) - itraverse' f (LowerEnvelope v0 vs) = - LowerEnvelope <$> (v0&incidentEdgesU.traverse %%~ liftF f unboundedVertexId) - <*> itraverse (\i vtx -> vtx&incidentEdgesB.traverse %%~ liftF f (i+1)) vs - - liftF :: (DartIx (LowerEnvelope' plane) -> LEEdge plane -> f (LEEdge plane)) - -> VertexIx (LowerEnvelope' plane) - -> LEEdge plane -> f (LEEdge plane) - liftF f u e = f (u,e^.destination) e - {-# INLINE darts #-} - --------------------------------------------------------------------------------- +type instance NumType (MinimizationDiagram r plane) = r +type instance Dimension (MinimizationDiagram r plane) = 2 +instance Constrained (MinimizationDiagram r) where + type Dom (MinimizationDiagram r) plane = (Ord plane, NumType plane ~ r) -instance HasEdges' (LowerEnvelope' plane) where - type Edge (LowerEnvelope' plane) = LEEdge plane - type EdgeIx (LowerEnvelope' plane) = ( VertexIx (LowerEnvelope' plane) - , VertexIx (LowerEnvelope' plane) - ) - edgeAt = dartAt - {-# INLINE edgeAt #-} +instance CFunctor (MinimizationDiagram r) where + cmap f (MinimizationDiagram m) = MinimizationDiagram $ Map.mapKeys f m -instance HasEdges (LowerEnvelope' plane) (LowerEnvelope' plane) where - -- | Traversal of all edges in the graph; i.e. we traverse the edges only in the - -- direction from lower vertex id to higher vertexId. - edges = darts . ifiltered (\(u,v) _ -> u < v) - {-# INLINE edges #-} --- instance Graph_ (LowerEnvelope' plane) where -class HasUnboundedEdges t e | t -> e where - -- | Lens to access the unbounded edges. - unboundedEdges :: Lens' t (Seq.Seq e) +-- | Get the underlying Map that relates every plane in the envelope to its projected region +asMap :: MinimizationDiagram r plane -> Map plane (Region r (Point 2 r)) +asMap (MinimizationDiagram m) = m -instance HasUnboundedEdges (LowerEnvelope' plane) (LEEdge plane) where - unboundedEdges = theUnboundedVertex.incidentEdgesU --- | Edges are either halflines or Linesegments -type EdgeGeometry point = Either (HalfLine point) (ClosedLineSegment point) +-- | A region in the minimization diagram. The boundary is given in CCW order; i.e. the +-- region is to the left of the boundary. +data Region r point = Bounded (CircularList point) + | Unbounded (Vector 2 r) + -- ^ vector indicating the direction of the unbounded edge + -- incident to the first vertex. Note that this vector + -- thus points INTO vertex v. + (NonEmpty point) + -- ^ the vertices in CCW order, + (Vector 2 r) + -- ^ the vector indicating the direction of the unbounded + -- edge incident to the last vertex. The vector points + -- away from the vertex (i.e. towards +infty). + deriving stock (Show,Eq,Functor,Foldable,Traversable) --- | Get the projected geometries (halflines and edges line segments ) representing the --- edges. -projectedEdgeGeometries :: (Plane_ plane r, Ord r, Fractional r - , Show plane, Show r +type instance NumType (Region r point) = r +type instance Dimension (Region r point) = Dimension point - ) - => IndexedFold (EdgeIx (LowerEnvelope' plane)) - (LowerEnvelope' plane) (EdgeGeometry (Point 2 r)) -projectedEdgeGeometries = ifolding $ \env -> - (\(i,e) -> (i, projectedEdgeGeometry env i e)) <$> env^..edges.withIndex - -- TODO; this is kind of horrible, there must be a nicer way of writing this. +-- | bounded regions are really circular lists, but we just represent them as lists for +-- now. +type CircularList a = [a] --- | Computes the (projected) LineSegment or HalfLine representing a given edge. -projectedEdgeGeometry :: (Plane_ plane r, Ord r, Fractional r - , Show plane, Show r - ) - => LowerEnvelope' plane - -> EdgeIx (LowerEnvelope' plane) - -> Edge (LowerEnvelope' plane) - -> EdgeGeometry (Point 2 r) -projectedEdgeGeometry env (ui,vi) e = case env^?!vertexAt ui of - Left _unboundedVtx -> case env^?!vertexAt vi of - Left _ -> error "projectedEdgeGeometry: absurd edge from vInfty to vInfty" - Right v -> unboundedEdge v - Right u -> case env^?!vertexAt vi of - Left _unboundedVtx -> unboundedEdge u - Right v -> Right $ ClosedLineSegment (u^.location2) (v^.location2) - where - unboundedEdge v = let dir = edgeIntersectionLine e ^.direction.to negated - in Left $ HalfLine (v^.location2) dir - -- edge e is oriented from the lowerId towards the higherId, so in particular *from* - -- the unbounded vertex into the bounded vertex that means that to construct the - -- halfline we actually wish to flip the orientation of the line -------------------------------------------------------------------------------- diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/VertexForm.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/VertexForm.hs new file mode 100644 index 000000000..f1b438a31 --- /dev/null +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/VertexForm.hs @@ -0,0 +1,96 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : HGeometry.Plane.LowerEnvelope.VertexForm +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- The vertices of the lower envelope; i.e. the lower envelope in VertexForm +-- +-------------------------------------------------------------------------------- +module HGeometry.Plane.LowerEnvelope.Connected.VertexForm + ( VertexForm + , Definers(..) + , fromCCWList + , definers + ) where + +import qualified Data.Foldable as F +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Maybe (fromMaybe, listToMaybe) +import Data.Set (Set) +import qualified Data.Set as Set +import HGeometry.Combinatorial.Util +import HGeometry.HyperPlane.Class +import HGeometry.HyperPlane.NonVertical +import HGeometry.Intersection +import HGeometry.Line +import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap +import HGeometry.Plane.LowerEnvelope.Connected.Primitives +import HGeometry.Plane.LowerEnvelope.Connected.Type +import HGeometry.Point +import HGeometry.Properties +import HGeometry.Vector + +-------------------------------------------------------------------------------- +-- * The planes defining a vertex + +-- | The vertices of a lower envelope is just a Map with every vertex its definers, +-- i.e. the planes that define the vertex in CCW order around it. +type VertexForm r plane = Map (Point 3 r) (Definers plane) + +---------------------------------------- +-- * The planes defining a vertex + +-- | in CCW order, starting with the plane that is minimal at the vertical up direction +-- from their common vertex. +newtype Definers plane = Definers [plane] + deriving stock (Show,Eq,Ord) + deriving newtype (Functor,Foldable) + +-- | Given the planes in order, starting with the one that is closest in the up direction, +-- construct the Definers. +fromCCWList :: [plane] -> Definers plane +fromCCWList = Definers + +-- | Smart constructor for creating the definers of three planes +definers :: forall plane r.(Plane_ plane r, Ord r, Fractional r) + => Three plane -> Maybe (Point 3 r, Definers plane) +definers (Three h1@(Plane_ a1 b1 c1) h2 h3) = + do l12 <- intersectionLine h1 h2 + l13 <- intersectionLine h1 h3 + intersect l12 l13 >>= \case + Line_x_Line_Line _ -> Nothing + Line_x_Line_Point (Point2 x y) -> Just ( Point3 x y (a1 * x + b1* y + c1) + , Definers [hMin, hTwo, hThree] + ) + where + (hMin,h,h') = extractMinOn (evalAt $ Point2 x (y+1)) h1 h2 h3 + -- we compute the plane hMin that is cheapest directly above the vertex h and + -- h' are the other two planes. That Means hMin is the first definer (in the + -- CCW order). What remains is to determine the order in which the remaining + -- planes h and h' appear in the order. + (hTwo,hThree) = case ccwCmpAroundWith (Vector2 0 1) (origin :: Point 2 r) + (Point vMinH) (Point vMinH') of + LT -> (h,h') + EQ -> error "definers: weird degeneracy?" + GT -> (h',h) + + vMinH = fromMaybe err $ intersectionVector h hMin + vMinH' = fromMaybe err $ intersectionVector h' hMin + err = error "definers: absurd" + +-- | given three elements, returns the minimum element in the first argument and the +-- remaining two elements in the second and third argument (in arbitrary order). +extractMinOn :: Ord a => (c -> a) -> c -> c -> c -> (c, c, c) +extractMinOn f a b c = let (m,ab) = min' a b + (mi,c') = min' m c + in (mi, ab, c') + where + min' x y + | f x <= f y = (x,y) + | otherwise = (y,x) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/ConnectedNew.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/ConnectedNew.hs deleted file mode 100644 index c835fd909..000000000 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/ConnectedNew.hs +++ /dev/null @@ -1,19 +0,0 @@ --------------------------------------------------------------------------------- --- | --- Module : HGeometry.Plane.LowerEnvelope.ConnectedNew --- Copyright : (C) Frank Staals --- License : see the LICENSE file --- Maintainer : Frank Staals --- --- A Representation of the Lower envelope of planes as a bunch of convex regions --- --------------------------------------------------------------------------------- -module HGeometry.Plane.LowerEnvelope.ConnectedNew - ( module HGeometry.Plane.LowerEnvelope.Connected.Primitives - , module HGeometry.Plane.LowerEnvelope.Connected.Regions - , module HGeometry.Plane.LowerEnvelope.Connected.BruteForce - ) where - -import HGeometry.Plane.LowerEnvelope.Connected.Primitives -import HGeometry.Plane.LowerEnvelope.Connected.BruteForce -import HGeometry.Plane.LowerEnvelope.Connected.Regions diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs index dffc6b1d8..e67c8e821 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Naive.hs @@ -1,25 +1,28 @@ {-# LANGUAGE UndecidableInstances #-} module HGeometry.Plane.LowerEnvelope.Naive ( lowerEnvelope - , lowerEnvelopeVertexForm - -- , triangulatedLowerEnvelope - - , asVertex - , belowAll + , lowerEnvelopeWith ) where -------------------------------------------------------------------------------- -import Control.Monad (guard) +import Control.Lens +import Data.Either (partitionEithers) import Data.Foldable1 -import qualified Data.Set as Set -import HGeometry.Combinatorial.Util +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Maybe +import HGeometry.Ext import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical -import HGeometry.Plane.LowerEnvelope.AdjListForm (LowerEnvelope, fromVertexForm) -import HGeometry.Plane.LowerEnvelope.VertexForm +import HGeometry.Line +import HGeometry.Line.General +import qualified HGeometry.Line.LowerEnvelope as Line +import HGeometry.Plane.LowerEnvelope.Connected +import HGeometry.Plane.LowerEnvelope.Type import HGeometry.Point - +import HGeometry.Sequence.Alternating (firstWithNeighbors) +import HGeometry.Vector -------------------------------------------------------------------------------- -- | Brute force implementation that computes the lower envelope, by explicitly @@ -29,79 +32,104 @@ import HGeometry.Point -- -- -- running time: \(O(n^4 )\) -lowerEnvelope :: ( Plane_ plane r - , Ord r, Fractional r, Foldable1 f, Functor f, Ord plane - , Show plane, Show r - ) => f plane -> LowerEnvelope plane -lowerEnvelope hs = fromVertexForm hs $ lowerEnvelopeVertexForm hs - --------------------------------------------------------------------------------- - +lowerEnvelope :: ( Plane_ plane r + , Ord r, Fractional r, Foldable1 f, Functor f, Ord plane + , Show plane, Show r + ) => f plane -> LowerEnvelope plane +lowerEnvelope = lowerEnvelopeWith bruteForceLowerEnvelope -------------------------------------------------------------------------------- --- * Computing a lower envelope in vertex form --- | Brute force implementation that computes the vertex form of the --- lower envelope, by explicitly considering every triple of planes. --- --- pre: the input forms a *set* of planes, i.e. no duplicates --- +-- | Computes the lower envelope of the planes; possibly using the given function to +-- construct the minimization diagram. -- --- general position assumptions: None. Though keep in mind that --- e.g. if all planes are parallel there are no vertices. +-- In particular, we distinguish whether the planes define a connected lower envelope +-- (in that case, we actually call the algorithm on it). Or that the planes are degenerate +-- and are all parallel and thus define parallel strips. In that case, we use a two dimensional +-- lower envelope algorithm to compute the lower envelope of the (parallel) planes. -- +-- \(O(T(n) + n \log n)\). +lowerEnvelopeWith :: forall nonEmpty plane r. + ( Plane_ plane r + , Ord r, Fractional r, Foldable1 nonEmpty + ) + => (NonEmpty plane -> MinimizationDiagram r plane) + -> nonEmpty plane -> LowerEnvelope plane +lowerEnvelopeWith minimizationDiagram hs = case distinguish (toNonEmpty hs) of + Left lines' -> ParallelStrips . fromLines $ lines' + Right (Vector3 h1 h2 h3, hs') -> ConnectedEnvelope $ minimizationDiagram (h1 :| h2 : h3 : hs') + where + fromLines = firstWithNeighbors (\h _ h' -> fromMaybe err $ intersectionLine h h') + . fmap (view extra) . view Line._Alternating + . lowerEnvelope' + err = error "lowerEnvelopeWith. absurd: neighbouring planes must intersect" + + lowerEnvelope' :: NonEmpty (LineEQ r :+ plane) + -> Line.LowerEnvelope (Point 2 r) (LineEQ r :+ plane) + lowerEnvelope' = Line.lowerEnvelope + -- for whatever reason GHC 9.8 does not want to compile this without the type signature + -- things are fine with 9.10 though; so maybe at some point we an just inline this. + +-- | We distinguish whether the planes are all parallel (left), and thus actually just +-- define a 2D lower envelope of lines, or that the planes actually have a connected lower +-- envelope. In that case, we return three planes that certify that the lower envelope +-- will be connected, as well as the rest of the planes. -- +-- (Note that those three planes that we return do not necesarily define a vertex of the +-- envelope). -- --- running time: \(O(n^4 )\) -lowerEnvelopeVertexForm :: ( Plane_ plane r - , Ord r, Fractional r, Foldable f, Ord plane - ) => f plane -> VertexForm plane -lowerEnvelopeVertexForm hs = foldMap (maybe mempty singleton . asVertex hs) $ uniqueTriplets hs - --- | Given all planes, and a triple, computes if the triple defines a --- vertex of the lower envelope, and if so returns it. -asVertex :: (Plane_ plane r, Foldable f, Ord plane, Ord r, Fractional r) - => f plane -> Three plane -> Maybe (LEVertex plane) -asVertex hs t@(Three h1 h2 h3) = do v <- intersectionPoint t - guard (v `belowAll` hs) - pure $ LEVertex v (Set.fromList [h1,h2,h3]) - --- | 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 -> verticalSideTest v h /= GT) -{-# INLINE belowAll #-} +-- running time: \(O(n)\) +distinguish :: ( Plane_ plane r, Eq r, Fractional r + ) + => NonEmpty plane -> Either (NonEmpty (LineEQ r :+ plane)) + (Vector 3 plane, [plane]) +distinguish hs@(h :| rest) = case findDifferent lines' of + Left baseLine -> Left $ fmap (flip asLine baseLine) hs + Right (Vector2 h2 h3, rest') -> Right (Vector3 h h2 h3, rest') - - - -{- - -triangulatedLowerEnvelope :: ( Plane_ plane r - , Ord r, Fractional r, Foldable f - ) => f plane -> LowerEnvelope plane -triangulatedLowerEnvelope hs = undefined - --} --------------------------------------------------------------------------------- - -{- - --- TODO: attach the two defining halfplanes to the result - --- | Given two halfplanes h and h' computes the halfplane where h lies --- vertically below h'. -liesBelowIn :: (Plane_ plane r, Ord r, Fractional r) - => plane -> plane -> Maybe (HalfPlane r) -liesBelowIn (Plane_ a b c) (Plane_ a' b' c') = case b `compare` b' of - LT -> Just $ Above (LineEQ d e) - GT -> Just $ Below (LineEQ d e) - EQ -> case a `compare` a' of - LT -> Just $ RightOf f - GT -> Just $ LeftOf f - EQ -> Nothing where - d = (a-a') / (b - b') - e = (c-c') / (b - b') - f = (c-c') / (a - a') - --} + lines' = map (\h' -> (perpendicularTo' <$> intersectionLine h h') :+ h') rest + -- for each plane, compute (the downward projection of) its line of intersection with + -- with h. We actually all rotate these lines by 90 degrees. + -- + -- if these lines are all parallel, then our input was degenerate. We pick one of + -- those lines baseLine, and compute the intersection of the plane h with the vertical + -- wall through this + +perpendicularTo' :: (Fractional r, Eq r) => VerticalOrLineEQ r -> VerticalOrLineEQ r +perpendicularTo' = \case + VerticalLineThrough x -> NonVertical $ LineEQ 0 x + NonVertical (LineEQ 0 b) -> VerticalLineThrough b + NonVertical (LineEQ a b) -> NonVertical $ LineEQ (1/a) b + +-- | Given the planes tagged with their intersection line with some plane h0, returns +-- either a Left baseLine if all these lines are parallel, or two planes that have +-- non-paralel intersection lines (as well as the rest of the planes). +-- +-- the baseLine is one of the parallel lines (or vertical if there are no lines to begin with). +findDifferent :: (Eq r, Num r + ) + => [Maybe (VerticalOrLineEQ r) :+ plane] -> Either (VerticalOrLineEQ r) + (Vector 2 plane, [plane]) +findDifferent xs = case partitionEithers $ map distinguishParallel xs of + (_,[]) -> Left $ VerticalLineThrough 0 -- all planes happen to be paralell + (pars0, (l :+ h) : rest) -> case List.span (isParallelTo l) rest of + (_,[]) -> Left l -- degnerate; planes are all parralel + (pars1, (_ :+ h') : rest') -> Right ( Vector2 h h' + , pars0 <> map (view extra) (pars1 <> rest') + ) + where + distinguishParallel (m :+ h) = case m of + Nothing -> Left h + Just l -> Right (l :+ h) + +-- | given a plane h, and a line l, erect a vertical plane through l, and return the line +-- in which h intersects this vertical plane. +asLine :: (Plane_ plane r, Num r) + => plane -> VerticalOrLineEQ r -> LineEQ r :+ plane +asLine h@(Plane_ a b c) = \case + VerticalLineThrough x -> LineEQ b (a*x + c) :+ h + NonVertical (LineEQ a' b') -> LineEQ (a + a'*b) (b*b' + c) :+ h + -- we simply fill in the coordinates of the line into the equation for h. this gives us + -- line (either parameterized in terms of the y-coordinate, or in case of the + -- x-coordinate of R^3 in the second case). diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Type.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Type.hs index f3a7c8194..57397f412 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Type.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Type.hs @@ -1,127 +1,48 @@ {-# LANGUAGE UndecidableInstances #-} module HGeometry.Plane.LowerEnvelope.Type - ( VertexID - , BoundedVertexF(Vertex) - , location, definers, location2, incidentEdgesB - , traverseBoundedV - - , LEEdge(Edge) - -- , createEdge - , destination, leftPlane, rightPlane - , flipEdge + ( LowerEnvelope(..) + , pointLocateParallel ) where import Control.Lens import qualified Data.Set as Set +import qualified Data.Vector as Vector +import HGeometry.Algorithms.BinarySearch +import HGeometry.HyperPlane.Class +import HGeometry.HyperPlane.NonVertical +import HGeometry.Line.General +import HGeometry.Plane.LowerEnvelope.Connected import HGeometry.Point import HGeometry.Properties +import HGeometry.Sequence.Alternating (Alternating(..)) +import qualified HGeometry.Sequence.Alternating as Alternating -------------------------------------------------------------------------------- +-- * Data type defining a lower envelope --- | VertexIds are just ints -type VertexID = Int - --- | Orders the plane around their intersection point. -newtype AroundVertex plane = AroundVertex plane - deriving stock (Show,Eq,Functor,Foldable,Traversable) - --- instance Ord (AroundVertex plane) where --- (AroundVertex h) `compare` (AroundVertex h') = undefined --- where --- -- the vertex location; aquire this using reflection somehow --- v :: Point 2 (NumType plane) --- v = undefined - - -- I guess we somehow compute the bisector, and figure out if h lies above the bisector or below. - - --- | Data type representing the bounded vertices of the lower --- envelope. --- --- In the vertexform, the f will be Const (), i.e. we don't have --- adjacency information. But in case of Adjacencylist form we can --- implement it using a Sequence. -data BoundedVertexF f plane = Vertex { _location :: !(Point 3 (NumType plane)) - , _definers :: Set.Set plane - , _incidentEdgesB :: f (LEEdge plane) - -- ^ incident edges, in CCW order. - } - -deriving instance ( Show plane, Show (NumType plane) - , Show (f (LEEdge plane)) - ) => Show (BoundedVertexF f plane) -deriving instance ( Eq plane - , Eq (NumType plane) - , Eq (f (LEEdge plane)) - ) => Eq (BoundedVertexF f plane) - --- instance Functor (BoundedVertex ) - - - --- | Traverse the bounded vertex. --- --- \(O(k\log k)\), where \(k\) is the number of definers of v. -traverseBoundedV :: (Traversable f, Applicative g, NumType plane ~ NumType plane', Ord plane') - => (plane -> g plane') -> BoundedVertexF f plane -> g (BoundedVertexF f plane') -traverseBoundedV f = \case - Vertex l defs es -> Vertex l <$> traverseSet f defs <*> traverse (traverse f) es - --- | Traverse on a Set. -traverseSet :: (Ord b, Applicative f) => (a -> f b) -> Set.Set a -> f (Set.Set b) -traverseSet f = fmap Set.fromList . traverse f . Set.toAscList +-- | The lower enevelope of planes in R^3. (Or rather, its minimization diagram) +data LowerEnvelope plane = + ParallelStrips !(Alternating Vector.Vector (VerticalOrLineEQ (NumType plane)) plane) + | ConnectedEnvelope !(MinimizationDiagram (NumType plane) plane) +deriving instance (Show plane, Show (NumType plane)) => Show (LowerEnvelope plane) +deriving instance (Eq plane, Eq (NumType plane)) => Eq (LowerEnvelope plane) - --- | The location of the vertex -location :: (NumType plane ~ r) => Lens' (BoundedVertexF f plane) (Point 3 r) -location = lens _location (\v l -> v { _location = l }) - --- | Projected 2d location of the point -location2 :: (NumType plane ~ r) => Getter (BoundedVertexF f plane) (Point 2 r) -location2 = location . to projectPoint - --- | The three planes defining the vertex -definers :: Lens' (BoundedVertexF f plane) (Set.Set plane) -definers = lens _definers (\v ds -> v { _definers = ds }) - --- | Lens to access the incident edges of a vertex -incidentEdgesB :: Lens (BoundedVertexF f plane) - (BoundedVertexF g plane) (f (LEEdge plane)) (g (LEEdge plane)) -incidentEdgesB = lens _incidentEdgesB (\v es -> v { _incidentEdgesB = es }) - -------------------------------------------------------------------------------- --- | An (half)Edge in the Lower envelope -data LEEdge plane = Edge { _destination :: {-# UNPACK #-}!VertexID - , _leftPlane :: plane - , _rightPlane :: plane - } deriving stock (Show,Eq,Ord,Functor,Foldable,Traversable) - --- -- | Create an edge --- createEdge :: VertexID -- ^ the destination --- -> plane -- ^ the left plane --- -> plane -- ^ the right plane --- -> LEEdge plane --- createEdge u hl hr = Edge u hl hr - --- | Given some vertex u and an edge e from u towards some other --- vertex v, flip the edge e, s othat it is the edge from v to u. -flipEdge :: VertexID -> LEEdge plane -> LEEdge plane -flipEdge u (Edge _ hl hr) = Edge u hr hl - --- | Getter to access the destination field of an edge. -destination :: Getter (LEEdge plane) VertexID -destination = to _destination - - --- | Lens to access the plane left of/above the edge -leftPlane :: Lens' (LEEdge plane) plane -leftPlane = lens _leftPlane (\ed h -> ed { _leftPlane = h }) - --- | Lens to access the plane right of/below the edge -rightPlane :: Lens' (LEEdge plane) plane -rightPlane = lens _rightPlane (\ed h -> ed { _rightPlane = h }) - --------------------------------------------------------------------------------- +-- | Point locates the given point among the parralel strips. +-- +-- \(O(\log h)\), where \(h\) is the number of planes on the lower envelope +pointLocateParallel :: (Plane_ plane r, Point_ point 2 r, Ord r, Num r) + => point + -> Alternating Vector.Vector (VerticalOrLineEQ r) plane + -> plane +pointLocateParallel q (Alternating h0 hs) = case binarySearchIn (q `liesRightOf`) hs of + AllTrue (_,h) -> h + FlipsAt _ (_,h) -> h + AllFalse _ -> h0 + where + q' `liesRightOf` (sep, _) = case sep of + VerticalLineThrough x -> q^.xCoord > x + NonVertical l -> verticalSideTest q l == GT diff --git a/hgeometry/src/HGeometry/VoronoiDiagram.hs b/hgeometry/src/HGeometry/VoronoiDiagram.hs index dc20d0d86..18fbaf669 100644 --- a/hgeometry/src/HGeometry/VoronoiDiagram.hs +++ b/hgeometry/src/HGeometry/VoronoiDiagram.hs @@ -13,7 +13,7 @@ module HGeometry.VoronoiDiagram , VoronoiDiagram'(..) , voronoiDiagram , voronoiVertices - , edgeGeometries + , asMap ) where import HGeometry.VoronoiDiagram.ViaLowerEnvelope diff --git a/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs b/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs index 2744c1003..81bec6937 100644 --- a/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs +++ b/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs @@ -13,31 +13,39 @@ module HGeometry.VoronoiDiagram.ViaLowerEnvelope ( VoronoiDiagram(..) , VoronoiDiagram'(..) - , ColinearPoint + , asMap , voronoiDiagram , voronoiVertices - , edgeGeometries + -- , edgeGeometries ) where import Control.Lens +import Control.Subcategory.Functor import Data.Foldable1 +import qualified Data.List.NonEmpty as NonEmpty +import qualified Data.Map as Map +import Data.Set (Set) import qualified Data.Set as Set -import HGeometry.Box +import qualified Data.Vector as Vector import HGeometry.Duality import HGeometry.Ext import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical -import HGeometry.Plane.LowerEnvelope.AdjListForm -import HGeometry.Plane.LowerEnvelope.Naive (lowerEnvelopeVertexForm) -import HGeometry.Plane.LowerEnvelope.VertexForm (VertexForm, vertices') +import HGeometry.Line.General +import HGeometry.Plane.LowerEnvelope ( MinimizationDiagram, Region(..) + , lowerEnvelope, LowerEnvelope(..) + ) +import qualified HGeometry.Plane.LowerEnvelope as LowerEnvelope import HGeometry.Point import HGeometry.Properties +import HGeometry.Sequence.Alternating (Alternating(..)) -------------------------------------------------------------------------------- -- | A Voronoi diagram -data VoronoiDiagram point = AllColinear !(Set.Set (ColinearPoint point)) - | ConnectedVD !(VoronoiDiagram' point) +data VoronoiDiagram point = + AllColinear !(Alternating Vector.Vector (VerticalOrLineEQ (NumType point)) point) + | ConnectedVD !(VoronoiDiagram' point) deriving instance (Show point, Show (NumType point)) => Show (VoronoiDiagram point) deriving instance (Eq point, Eq (NumType point)) => Eq (VoronoiDiagram point) @@ -45,22 +53,11 @@ deriving instance (Eq point, Eq (NumType point)) => Eq (VoronoiDiagram poi type instance NumType (VoronoiDiagram point) = NumType point type instance Dimension (VoronoiDiagram point) = 2 -- Dimension point --- | Just a newtype around point; used to model a set of points that are all colinear in --- the Vornoi diagram. -newtype ColinearPoint point = ColinearPoint point - deriving (Show,Eq) - -instance Wrapped (ColinearPoint point) where - type Unwrapped (ColinearPoint point) = point - _Wrapped' = coerced - --- instance Rewrapped (ColinearPoint point) point -------------------------------------------------------------------------------- -- | A connected VoronoiDiagram -newtype VoronoiDiagram' point = - VoronoiDiagram (LowerEnvelope' (Plane (NumType point) :+ point)) +newtype VoronoiDiagram' point = VoronoiDiagram (MinimizationDiagram (NumType point) point) deriving instance (Show point, Show (NumType point)) => Show (VoronoiDiagram' point) deriving instance (Eq point, Eq (NumType point)) => Eq (VoronoiDiagram' point) @@ -70,58 +67,70 @@ type instance Dimension (VoronoiDiagram' point) = 2 -- Dimension point -- | Iso to Access the underlying LowerEnvelope _VoronoiDiagramLowerEnvelope :: Iso (VoronoiDiagram' point) (VoronoiDiagram' point') - (LowerEnvelope' (Plane (NumType point) :+ point)) - (LowerEnvelope' (Plane (NumType point') :+ point')) + (MinimizationDiagram (NumType point) point) + (MinimizationDiagram (NumType point') point') _VoronoiDiagramLowerEnvelope = coerced --------------------------------------------------------------------------------- +-- | Get, for each point, its Voronoi region +asMap :: (Point_ point 2 r, Ord point) + => VoronoiDiagram' point -> Map.Map point (Region r (Point 2 r)) +asMap = LowerEnvelope.asMap . view _VoronoiDiagramLowerEnvelope -instance (Ord (NumType point), Num (NumType point)) => IsBoxable (VoronoiDiagram' point) where - boundingBox vd = projectPoint <$> boundingBox (vd^._VoronoiDiagramLowerEnvelope) +-------------------------------------------------------------------------------- -------------------------------------------------------------------------------- -- | Computes the Voronoi Diagram, by lifting the points to planes, and computing -- the lower envelope of these planes. -- +-- O(n^4) ( we currently use the brute force implentation...) TODO: switch to something faster + -- \(O(n\log n)\) -voronoiDiagram :: ( Point_ point 2 r, Functor f, Ord point - , Ord r, Fractional r, Foldable1 f - , Show point, Show r - ) => f point -> VoronoiDiagram point -voronoiDiagram pts = case lowerEnvelope' . fmap (\p -> liftPointToPlane p :+ p) $ pts of - ParallelStrips strips -> AllColinear $ Set.mapMonotonic getPoint strips - ConnectedEnvelope env -> ConnectedVD $ VoronoiDiagram env - where - lowerEnvelope' hs = fromVertexForm hs $ upperEnvelopeVertexForm hs - getPoint = view (_Wrapped'.extra.to ColinearPoint) - --- | Computes all voronoi vertices -voronoiVertices :: ( Point_ point 2 r, Functor f, Ord point - , Ord r, Fractional r, Foldable f - ) => f point -> [Point 2 r] -voronoiVertices = map (projectPoint . fst) - . itoListOf vertices' - . upperEnvelopeVertexForm - . fmap (\p -> liftPointToPlane p :+ p) --- FIXME: get rid of the ord point constraint - --- | Computes the vertex form of the upper envelope. The z-coordinates are still flipped. -upperEnvelopeVertexForm :: ( Plane_ plane r - , Ord r, Fractional r, Foldable f, Functor f, Ord plane - ) => f plane -> VertexForm plane -upperEnvelopeVertexForm = lowerEnvelopeVertexForm . fmap flipZ +voronoiDiagram :: ( Point_ point 2 r, Functor f, Ord point + , Ord r, Fractional r, Foldable1 f + , Show point, Show r + ) => f point -> VoronoiDiagram point +voronoiDiagram = voronoiDiagramWith lowerEnvelope + +-- | Given a function to compute a lower envelope; construct use it to construct the +-- Voronoi diagram. +voronoiDiagramWith :: ( Point_ point 2 r, Functor nonEmpty, Ord point + , Ord r, Fractional r, Foldable1 nonEmpty + ) + => (nonEmpty (Plane r :+ point) -> LowerEnvelope (Plane r :+ point)) + + -> nonEmpty point + -> VoronoiDiagram point +voronoiDiagramWith lowerEnv pts = case lowerEnv . fmap (\p -> pointToPlane p :+ p) $ pts of + ParallelStrips strips -> AllColinear $ fmap (^.extra) strips + ConnectedEnvelope env -> ConnectedVD . VoronoiDiagram . cmap (^.extra) $ env where + -- lifts the point to a plane; so that the lower envelope corresponds to the VD + pointToPlane = flipZ . liftPointToPlane flipZ = over (hyperPlaneCoefficients.traverse) negate + +-- | Compute the vertices of the Voronoi diagram +voronoiVertices :: ( Point_ point 2 r, Functor f, Ord point + , Ord r, Fractional r, Foldable1 f + , Show point, Show r + , Ord point + ) => f point -> Set (Point 2 r) +voronoiVertices pts = case voronoiDiagram pts of + AllColinear _ -> mempty + ConnectedVD vd -> foldMap (\case + Bounded vs -> Set.fromList vs + Unbounded _ vs _ -> Set.fromList (NonEmpty.toList vs) + ) (asMap vd) + -------------------------------------------------------------------------------- --- | Get the halflines and line segments representing the VoronoiDiagram -edgeGeometries :: (Point_ point 2 r, Ord r, Fractional r +-- -- | Get the halflines and line segments representing the VoronoiDiagram +-- edgeGeometries :: (Point_ point 2 r, Ord r, Fractional r - , Show point, Show r - ) - => Fold (VoronoiDiagram' point) (EdgeGeometry (Point 2 r)) -edgeGeometries = _VoronoiDiagramLowerEnvelope.projectedEdgeGeometries --- TODO: figure out if this can be an indexed fold +-- , Show point, Show r +-- ) +-- => Fold (VoronoiDiagram' point) (EdgeGeometry (Point 2 r)) +-- edgeGeometries = _VoronoiDiagramLowerEnvelope.projectedEdgeGeometries +-- -- TODO: figure out if this can be an indexed fold diff --git a/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs b/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs index dcd6b4575..5e6bebfa3 100644 --- a/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs +++ b/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs @@ -7,6 +7,7 @@ module Plane.LowerEnvelopeSpec ) where import Control.Lens +import Data.Coerce import Data.Foldable import Data.Foldable1 import Data.List.NonEmpty (NonEmpty(..)) @@ -20,21 +21,24 @@ import HGeometry.Duality import HGeometry.Ext import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical +import HGeometry.Line import HGeometry.LineSegment import HGeometry.Number.Real.Rational +import HGeometry.Plane.LowerEnvelope.Connected import HGeometry.Plane.LowerEnvelope.Connected.Graph -import HGeometry.Plane.LowerEnvelope.ConnectedNew import HGeometry.Point import HGeometry.Polygon.Convex import HGeometry.Polygon.Simple +import HGeometry.Sequence.Alternating (separators) import HGeometry.Vector +import HGeometry.VoronoiDiagram +import qualified HGeometry.VoronoiDiagram as VD import Ipe import Ipe.Color import System.OsPath import Test.Hspec import Test.Hspec.WithTempFile import Test.QuickCheck.Instances () - -------------------------------------------------------------------------------- type R = RealNumber 5 @@ -75,13 +79,27 @@ instance (HasDefaultIpeOut point, Point_ point 2 r, Fractional r, Ord r -- 83.83267 217.24478 l -- h -instance ( Point_ point 2 r, Fractional r, Ord r +instance ( Point_ point 2 r, Fractional r, Ord r, Ord point + , Show point, Show r + ) + => HasDefaultIpeOut (VoronoiDiagram point) where + type DefaultIpeOut (VoronoiDiagram point) = Group + defIO = \case + AllColinear colinearPts -> let sites = toList colinearPts + bisectors = toList $ separators colinearPts + in ipeGroup . concat $ + [ [ iO $ defIO b | b <- bisectors ] + , [ iO $ defIO (p^.asPoint) | p <- sites ] + ] + + ConnectedVD vd -> defIO vd +instance ( Point_ point 2 r, Fractional r, Ord r, Ord point , Show point, Show r ) - => HasDefaultIpeOut (MinimizationDiagram r point) where - type DefaultIpeOut (MinimizationDiagram r point) = Group - defIO = ipeGroup . zipWith render (cycle $ drop 3 basicNamedColors) . Map.assocs + => HasDefaultIpeOut (VoronoiDiagram' point) where + type DefaultIpeOut (VoronoiDiagram' point) = Group + defIO = ipeGroup . zipWith render (cycle $ drop 3 basicNamedColors) . Map.assocs . VD.asMap where render color (site, voronoiRegion) = iO' $ ipeGroup [ iO $ defIO (site^.asPoint) ! attr SStroke color @@ -110,45 +128,18 @@ spec = describe "lower envelope tests" $ do testIpe [osp|degenerate2.ipe|] [osp|degenerate2_out|] - testIpeGraph [osp|foo.ipe|] - [osp|foo_graph_out|] - - - --- | Computes the vertex form of the upper envelope. The z-coordinates are still flipped. - -type VoronoiDiagram' r point = MinimizationDiagram r point + -- testIpeGraph [osp|foo.ipe|] + -- [osp|foo_graph_out|] -voronoiDiagram' :: ( Point_ point 2 r, Functor f, Ord point - , Ord r, Fractional r, Foldable1 f - , Show point, Show r - ) => f point -> VoronoiDiagram' r point -voronoiDiagram' = Map.mapKeysMonotonic (view extra) - . bruteForceLowerEnvelope - . fmap (\p -> flipZ (liftPointToPlane p) :+ p) - where - flipZ :: Num r => Plane r -> Plane r - flipZ = over (hyperPlaneCoefficients.traverse) negate - -voronoiVertices :: ( Point_ point 2 r, Functor f, Ord point - , Ord r, Fractional r, Foldable1 f - , Show point, Show r - , Ord point - ) => f point -> Set (Point 2 r) -voronoiVertices = foldMap (\case - Bounded pts -> Set.fromList pts - Unbounded _ pts _ -> Set.fromList (NonEmpty.toList pts) - ) . voronoiDiagram' - -- | Build voronoi diagrams on the input points testIpe :: OsPath -> OsPath -> Spec testIpe inFp outFp = do (points :: NonEmpty (Point 2 R :+ _)) <- runIO $ do inFp' <- getDataFileName ([osp|test-with-ipe/VoronoiDiagram/|] <> inFp) NonEmpty.fromList <$> readAllFrom inFp' - let vd = voronoiDiagram' $ view core <$> points + let vd = voronoiDiagram $ view core <$> points vv = voronoiVertices $ view core <$> points out = [ iO' points , iO' vd @@ -164,22 +155,24 @@ theEdges = ipeGroup . Map.foldMapWithKey (\v (adjs, _) -> ]) adjs) --- build a triangulated graph from the points in the input file -testIpeGraph :: OsPath -> OsPath -> Spec -testIpeGraph inFp outFp = do - (points :: NonEmpty (Point 2 R :+ _)) <- runIO $ do - inFp' <- getDataFileName ([osp|test-with-ipe/VoronoiDiagram/|] <> inFp) - NonEmpty.fromList <$> readAllFrom inFp' - let vd = voronoiDiagram' $ view core <$> points - vv = voronoiVertices $ view core <$> points - gr = toPlaneGraph $ Map.mapKeysMonotonic liftPointToPlane vd - out = [ iO' points - , iO' vd - , iO' $ theEdges gr - ] <> [ iO'' v $ attr SStroke red | v <- Set.toAscList vv ] - goldenWith [osp|data/test-with-ipe/Plane/LowerEnvelope/|] - (ipeFileGolden { name = outFp }) - (addStyleSheet opacitiesStyle $ singlePageFromContent out) +-- -- build a triangulated graph from the points in the input file +-- testIpeGraph :: OsPath -> OsPath -> Spec +-- testIpeGraph inFp outFp = do +-- (points :: NonEmpty (Point 2 R :+ _)) <- runIO $ do +-- inFp' <- getDataFileName ([osp|test-with-ipe/VoronoiDiagram/|] <> inFp) +-- NonEmpty.fromList <$> readAllFrom inFp' +-- case voronoiDiagram $ view core <$> points of +-- AllColinear _ -> fail "input points are colinear !?" +-- ConnectedVD vd -> do +-- let vv = voronoiVertices $ view core <$> points +-- gr = toPlaneGraph . Map.mapKeys liftPointToPlane $ VD.asMap vd -- TODO: we should flip the z no? +-- out = [ iO' points +-- , iO' vd +-- , iO' $ theEdges gr +-- ] <> [ iO'' v $ attr SStroke red | v <- Set.toAscList vv ] +-- goldenWith [osp|data/test-with-ipe/Plane/LowerEnvelope/|] +-- (ipeFileGolden { name = outFp }) +-- (addStyleSheet opacitiesStyle $ singlePageFromContent out) diff --git a/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs b/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs index fc4ae4ed2..c8a1b9bf1 100644 --- a/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs +++ b/hgeometry/test-with-ipe/test/VoronoiDiagram/VoronoiSpec.hs @@ -12,8 +12,8 @@ import Golden import HGeometry.Duality import HGeometry.Ext import HGeometry.Plane.LowerEnvelope -import HGeometry.Plane.LowerEnvelope.AdjListForm import qualified Data.Set as Set +import qualified Data.Vector as Vector import qualified Data.List.NonEmpty as NonEmpty import qualified Data.Sequence as Seq import HGeometry.Number.Real.Rational @@ -21,6 +21,7 @@ import Data.List.NonEmpty (NonEmpty(..)) import HGeometry.Point import HGeometry.Vector import HGeometry.Box +import HGeometry.Line.General import HGeometry.HalfLine import HGeometry.VoronoiDiagram import Ipe @@ -29,6 +30,8 @@ import Test.Hspec.WithTempFile import System.OsPath import Test.Hspec import Test.QuickCheck.Instances () +import Plane.LowerEnvelopeSpec () -- imports the ipe instances for Voronoi Diagram +import HGeometry.Sequence.Alternating (Alternating(..)) -- import Test.Util -------------------------------------------------------------------------------- @@ -40,22 +43,26 @@ spec = describe "Voronoi diagram tests" $ do -- prop "voronoi vertex is center disk" $ \c -> -- voronoiVertices inputs it "vertices of a trivial voronoi diagram" $ - voronoiVertices inputs `shouldBe` [Point2 5 5] + voronoiVertices inputs `shouldBe` (Set.fromList [Point2 5 5]) -- it "a trivial voronoi diagram" $ -- voronoiDiagram inputs `shouldBe` trivialVD - it "geometries of the trivial VD correct" $ - trivialVD^..edgeGeometries - `shouldBe` [Left (HalfLine (Point2 5 5) (Vector2 1 0)) - ,Left (HalfLine (Point2 5 5) (Vector2 0 (-1))) - ,Left (HalfLine (Point2 5 5) (Vector2 (-1) 1)) - ] - goldenWith [osp|data/test-with-ipe/golden/|] - (ipeContentGolden { name = [osp|trivalVoronoi|] }) + -- it "geometries of the trivial VD correct" $ + -- trivialVD^..edgeGeometries + -- `shouldBe` [Left (HalfLine (Point2 5 5) (Vector2 1 0)) + -- ,Left (HalfLine (Point2 5 5) (Vector2 0 (-1))) + -- ,Left (HalfLine (Point2 5 5) (Vector2 (-1) 1)) + -- ] + + + goldenWith [osp|data/test-with-ipe/VoronoiDiagram/|] + (ipeContentGolden { name = [osp|trivialVoronoi|] }) [ iO' inputs , iO' trivialVD ] + degenerateTests + testIpe [osp|trivial.ipe|] [osp|trivial_out|] testIpe [osp|simplest.ipe|] @@ -69,25 +76,31 @@ spec = describe "Voronoi diagram tests" $ do testIpe [osp|foo.ipe|] [osp|foo_out|] - - - degenerateTests :: Spec degenerateTests = describe "degnereate inputs" $ do it "single point diagram" $ voronoiDiagram (NonEmpty.singleton $ Point2 1 (2 :: R)) `shouldBe` - AllColinear undefined -- TODO + AllColinear (Alternating (Point2 1 2) mempty) it "two point diagram" $ voronoiDiagram (NonEmpty.fromList [Point2 1 (2 :: R), Point2 3 2]) `shouldBe` - AllColinear undefined -- TODO + AllColinear (Alternating (Point2 1 2) (Vector.fromList [(VerticalLineThrough 2, Point2 3 2)])) it "multiple parallel point diagram" $ voronoiDiagram (NonEmpty.fromList [ Point2 x (2 :: R) | x <- fromInteger <$> [1..10] ]) `shouldBe` - AllColinear undefined -- TODO + AllColinear (Alternating (Point2 1 2) . Vector.fromList $ + [(VerticalLineThrough 1.5,Point2 2 2) + ,(VerticalLineThrough 2.5,Point2 3 2) + ,(VerticalLineThrough 3.5,Point2 4 2) + ,(VerticalLineThrough 4.5,Point2 5 2) + ,(VerticalLineThrough 5.5,Point2 6 2) + ,(VerticalLineThrough 6.5,Point2 7 2) + ,(VerticalLineThrough 7.5,Point2 8 2) + ,(VerticalLineThrough 8.5,Point2 9 2) + ,(VerticalLineThrough 9.5,Point2 10 2)]) -- goldenWith [osp|data/test-with-ipe/golden/|] -- (ipeContentGolden { name = [osp|voronoi|] }) @@ -102,52 +115,55 @@ grow d (Box p q) = Box (p&coordinates %~ subtract d) (q&coordinates %~ (+d)) -instance (HasDefaultIpeOut point, Point_ point 2 r, Fractional r, Ord r - , Show r, Show point - ) - => HasDefaultIpeOut (VoronoiDiagram point) where - type DefaultIpeOut (VoronoiDiagram point) = Group - defIO = \case - AllColinear _pts -> ipeGroup [] - ConnectedVD vd -> defIO vd - - -instance (HasDefaultIpeOut point, Point_ point 2 r, Fractional r, Ord r - , Show r, Show point - ) - => HasDefaultIpeOut (VoronoiDiagram' point) where - type DefaultIpeOut (VoronoiDiagram' point) = Group - defIO vd = ipeGroup $ vd^..edgeGeometries.to render - where - bRect = boundingBox $ defaultBox :| [grow 1 $ boundingBox vd] - render = \case - Left hl -> iO $ ipeHalfLineIn bRect hl - Right seg -> iO' seg - -inputs :: [Point 2 R] -inputs = [origin, Point2 10 10, Point2 10 0] - - -trivialVD :: VoronoiDiagram' (Point 2 R) -trivialVD = VoronoiDiagram $ LowerEnvelope vInfty (Seq.fromList [bv]) - where - vInfty = UnboundedVertex $ Seq.fromList [Edge 1 h2 h3 - ,Edge 1 h3 h1 - ,Edge 1 h1 h2 - ] - bv = Vertex (Point3 5 5 0) - (Set.fromList planes) - (Seq.fromList $ - [ Edge 0 h2 h1 - , Edge 0 h3 h2 - , Edge 0 h1 h3 - ] - ) - planes = map (\p -> liftPointToPlane p :+ p) inputs - (h1,h2,h3) = case planes of - [h1',h2',h3'] -> (h1',h2',h3') - _ -> error "absurd" - -- order of the planes is incorrect, as is the z-coord. +-- instance (HasDefaultIpeOut point, Point_ point 2 r, Fractional r, Ord r +-- , Show r, Show point +-- ) +-- => HasDefaultIpeOut (VoronoiDiagram point) where +-- type DefaultIpeOut (VoronoiDiagram point) = Group +-- defIO = \case +-- AllColinear _pts -> ipeGroup [] +-- ConnectedVD vd -> defIO vd + + +-- instance (HasDefaultIpeOut point, Point_ point 2 r, Fractional r, Ord r +-- , Show r, Show point +-- ) +-- => HasDefaultIpeOut (VoronoiDiagram' point) where +-- type DefaultIpeOut (VoronoiDiagram' point) = Group +-- defIO vd = ipeGroup $ vd^..edgeGeometries.to render +-- where +-- bRect = boundingBox $ defaultBox :| [grow 1 $ boundingBox vd] +-- render = \case +-- Left hl -> iO $ ipeHalfLineIn bRect hl +-- Right seg -> iO' seg + + + +inputs :: NonEmpty (Point 2 R) +inputs = NonEmpty.fromList [origin, Point2 10 10, Point2 10 0] + +trivialVD :: VoronoiDiagram (Point 2 R) +trivialVD = voronoiDiagram inputs + + -- VoronoiDiagram $ LowerEnvelope vInfty (Seq.fromList [bv]) + -- where + -- vInfty = UnboundedVertex $ Seq.fromList [Edge 1 h2 h3 + -- ,Edge 1 h3 h1 + -- ,Edge 1 h1 h2 + -- ] + -- bv = Vertex (Point3 5 5 0) + -- (Set.fromList planes) + -- (Seq.fromList $ + -- [ Edge 0 h2 h1 + -- , Edge 0 h3 h2 + -- , Edge 0 h1 h3 + -- ] + -- ) + -- planes = map (\p -> liftPointToPlane p :+ p) inputs + -- (h1,h2,h3) = case planes of + -- [h1',h2',h3'] -> (h1',h2',h3') + -- _ -> error "absurd" + -- -- order of the planes is incorrect, as is the z-coord. testIpe :: OsPath -> OsPath -> Spec @@ -159,7 +175,7 @@ testIpe inFp outFp = do vv = voronoiVertices $ (view core) <$> points out = [ iO' points , iO' vd - ] <> [ iO'' v $ attr SStroke red | v <- vv ] + ] <> [ iO'' v $ attr SStroke red | v <- Set.toAscList vv ] goldenWith [osp|data/test-with-ipe/VoronoiDiagram/|] - (ipeContentGolden { name = outFp }) - out + (ipeFileGolden { name = outFp }) + (addStyleSheet opacitiesStyle $ singlePageFromContent out) diff --git a/hgeometry/test/LowerEnvelope/NaiveSpec.hs b/hgeometry/test/LowerEnvelope/NaiveSpec.hs deleted file mode 100644 index 31a87f7ab..000000000 --- a/hgeometry/test/LowerEnvelope/NaiveSpec.hs +++ /dev/null @@ -1,87 +0,0 @@ -module LowerEnvelope.NaiveSpec where - -import Control.Lens --- import qualified Data.Map as Map -import qualified Data.Set as Set --- import qualified Data.Vector as Boxed -import HGeometry.Combinatorial.Util -import HGeometry.HyperPlane.Class -import HGeometry.HyperPlane.NonVertical -import HGeometry.Instances () -import HGeometry.Plane.LowerEnvelope -import HGeometry.Plane.LowerEnvelope.VertexForm -import HGeometry.Number.Real.Rational -import HGeometry.Point -import Test.Hspec -import Test.Hspec.QuickCheck --- import Test.QuickCheck -import Test.QuickCheck.Instances () - --------------------------------------------------------------------------------- - -type R = RealNumber 5 - --- see also the test-with-ipe/test/VoronoiDiagram/VornoiSpec tests --- - --- TODO: I think fromVertexForm should also work for singleton inputs; i.e. just three planes intersecting in a vertex . That would be a good first testcase I think. --- for the rest we may watn to use the VoronoiDiagram tests to generate the edges of the VD - -spec :: Spec -spec = describe "lowerEnvelope tests" $ do - prop "intersection on plane" $ \h1 h2 (h3 :: Plane R) -> - verifyOnPlane h1 h2 h3 - it "asBoundedVertex" $ - case inputs of - [h1,h2,h3] -> asVertex inputs (Three h1 h2 h3) `shouldBe` - Just (LEVertex (Point3 10 10 10) (Set.fromList inputs)) - _ -> expectationFailure "input needs 3 planes" - prop "belowall" $ \h1 h2 h3 -> - case asVertex [h1,h2,h3] (Three h1 h2 (h3 :: Plane R)) of - Nothing -> True - Just v -> (v^.location) `belowAll` [h1,h2,h3] -{- - - it "vertices inputs" $ - let [h1,h2,h3] = inputs - p = Point3 10 10 10 - v = Vertex $ BoundedVertex p (Set.fromList inputs) - in - _boundedVertices (lowerEnvelope inputs) - `shouldBe` - (Map.fromList [(p, Set.fromList inputs)]) - it "halfEdges inputs" $ - let [h1,h2,h3] = inputs - v = Vertex $ BoundedVertex (Point3 10 10 10) (Set.fromList inputs) - in - _halfEdges (lowerEnvelope inputs) - `shouldBe` - [ HalfEdge v UnboundedVertex h1 - , HalfEdge UnboundedVertex v h2 - , HalfEdge v UnboundedVertex h3 - , HalfEdge UnboundedVertex v h1 - , HalfEdge v UnboundedVertex h2 - , HalfEdge UnboundedVertex v h3 - ] - -- there still seems to be s.t. wrong with the order of the leftplanes --} - --- | verify that the intersection point indeed lis on all three planes -verifyOnPlane :: (Fractional r, Ord r) - => Plane r -> Plane r -> Plane r -> Bool -verifyOnPlane h1 h2 h3 = case intersectionPoint (Three h1 h2 h3) of - Nothing -> True - Just (Point3 x y _) -> allEqual $ - (evalAt $ Point2 x y) <$> Three h1 h2 h3 - where - allEqual (Three a b c) = a == b && b == c - - -myEnv = lowerEnvelopeVertexForm inputs --- myTriEnv = triangulatedLowerEnvelope inputs - -inputs :: [Plane R] -inputs = [ Plane 1 0 0 - , Plane 0 1 0 - , Plane 0 0 10 - ] diff --git a/hgeometry/test/LowerEnvelope/RegionsSpec.hs b/hgeometry/test/LowerEnvelope/RegionsSpec.hs index 221ed0f15..f0d26ae37 100644 --- a/hgeometry/test/LowerEnvelope/RegionsSpec.hs +++ b/hgeometry/test/LowerEnvelope/RegionsSpec.hs @@ -7,7 +7,7 @@ import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical import HGeometry.Instances () import HGeometry.Number.Real.Rational -import HGeometry.Plane.LowerEnvelope.ConnectedNew +import HGeometry.Plane.LowerEnvelope.Connected import HGeometry.Point import HGeometry.Vector import Test.Hspec @@ -40,7 +40,7 @@ spec = describe "lowerEnvelope tests" $ do it "singleton diagram" $ do let v = Point2 10 10 :: Point 2 R [h1,h2,h3] <- pure inputs - bruteForceLowerEnvelope inputs `shouldBe` + (asMap $ bruteForceLowerEnvelope inputs) `shouldBe` Map.fromList [ (h1, Unbounded (Vector2 1 1) (NonEmpty.singleton v) (Vector2 0 1)) , (h2, Unbounded (Vector2 (-1) 0) (NonEmpty.singleton v) (Vector2 (-1) (-1)))