feat: add Density class

src/shapepy/__init__.py

@@ -21,6 +21,8 @@
 __version__ = importlib.metadata.version("shapepy")
 
 set_level("shapepy", level="INFO")
+# set_level("shapepy.bool2d", level="DEBUG")
+# set_level("shapepy.rbool", level="DEBUG")
 
 
 if __name__ == "__main__":

src/shapepy/bool2d/base.py

@@ -16,6 +16,7 @@
 from ..loggers import debug
 from ..scalar.angle import Angle
 from ..scalar.reals import Real
+from .density import Density
 
 
 class SubSetR2:
@@ -149,7 +150,7 @@ def rotate(self, angle: Angle) -> SubSetR2:
         raise NotImplementedError
 
     @abstractmethod
-    def density(self, center: Point2D) -> Real:
+    def density(self, center: Point2D) -> Density:
         """
         Computes the density of the subset around given point
 
@@ -237,8 +238,8 @@ def scale(self, _):
     def rotate(self, _):
         return self
 
-    def density(self, center: Point2D) -> Real:
-        return 0
+    def density(self, center: Point2D) -> Density:
+        return Density.zero
 
 
 class WholeShape(SubSetR2):
@@ -302,8 +303,8 @@ def scale(self, _):
     def rotate(self, _):
         return self
 
-    def density(self, center: Point2D) -> Real:
-        return 1
+    def density(self, center: Point2D) -> Density:
+        return Density.one
 
 
 # pylint: disable=duplicate-code

src/shapepy/bool2d/boolean.py

@@ -224,7 +224,7 @@ def midpoints_one_shape(
             for j, segment in enumerate(jordan.parametrize()):
                 mid_point = segment(Fraction(1, 2))
                 density = shapeb.density(mid_point)
-                mid_point_in = (density > 0 and closed) or density == 1
+                mid_point_in = (float(density) > 0 and closed) or density == 1
                 if not inside ^ mid_point_in:
                     yield (i, j)
 

src/shapepy/bool2d/density.py

@@ -0,0 +1,168 @@
+"""
+Defines the Density class that is used for the subsets to verify
+the density of a given point a
+
+It's the Lebesgue density.
+"""
+
+from __future__ import annotations
+
+from typing import Iterable, Optional, Union
+
+from ..analytic.tools import find_minimum, where_minimum
+from ..geometry.integral import IntegrateJordan
+from ..geometry.jordancurve import JordanCurve
+from ..geometry.point import Point2D
+from ..geometry.segment import Segment
+from ..loggers import debug
+from ..rbool import EmptyR1, SingleR1, SubSetR1, create_interval, subset_length
+from ..scalar.angle import Angle, degrees
+from ..scalar.reals import Real
+from ..tools import Is, NotExpectedError, To
+
+
+@debug("shapepy.bool2d.density")
+def half_density_jordan(
+    segments: Iterable[Segment], point: Point2D
+) -> Density:
+    """Computes the value of the density when the point is at
+    an smooth edge of any of the given segments
+    """
+    for segment in segments:
+        deltax = segment.xfunc - point.xcoord
+        deltay = segment.yfunc - point.ycoord
+        radius_square = deltax * deltax + deltay * deltay
+        minimal = find_minimum(radius_square, [0, 1])
+        if minimal < 1e-6:
+            place = where_minimum(radius_square, [0, 1])
+            if not Is.instance(place, SingleR1):
+                raise NotExpectedError(f"Not single value: {place}")
+            parameter = To.finite(place.internal)
+            angle = segment(parameter, 1).angle
+            return line(angle)
+    raise NotExpectedError("Not found minimum < 1e-6")
+
+
+@debug("shapepy.bool2d.density")
+def lebesgue_density_jordan(
+    jordan: JordanCurve, point: Optional[Point2D] = (0.0, 0.0)
+) -> Density:
+    """Computes the lebesgue density number from jordan curve
+
+    Returns a value in the interval [0, 1]:
+    * 0 -> means the point is outside the interior region
+    * 1 -> means the point is completly inside the interior
+    * between 0 and 1, it's on the boundary
+    """
+    point = To.point(point)
+    box = jordan.box()
+    if point not in box:
+        return Density.zero if jordan.area > 0 else Density.one
+
+    segments = tuple(jordan.parametrize())
+    for i, segmenti in enumerate(segments):
+        if point == segmenti(0):
+            segmentj = segments[(i - 1) % len(segments)]
+            anglei = segmenti(0, 1).angle
+            anglej = segmentj(1, 1).angle
+            return sector(anglei, ~anglej)
+
+    turns = IntegrateJordan.turns(jordan, point)
+    density = turns if jordan.area > 0 else 1 + turns
+    if density == 0.5:
+        return half_density_jordan(segments, point)
+    return Density.one if round(density) == 1 else Density.zero
+
+
+@debug("shapepy.bool2d.density")
+def line(angle: Angle) -> Density:
+    """Creates a Density of value 0.5 aligned with given angle"""
+    angle = To.angle(angle)
+    return sector(angle, angle + degrees(180))
+
+
+@debug("shapepy.bool2d.density")
+def sector(anglea: Angle, angleb: Angle) -> Density:
+    """Creates a Density instance within given two angles"""
+    uvala: Real = To.angle(anglea).turns % 1
+    uvalb: Real = To.angle(angleb).turns % 1
+    if uvalb == 0:
+        subset = create_interval(uvala, 1)
+    elif uvala < uvalb:
+        subset = create_interval(uvala, uvalb)
+    else:
+        subset = create_interval(0, uvalb) | create_interval(uvala, 1)
+    return Density(subset)
+
+
+class Density:
+    """
+    Density class that stores the sectors of circles that allows computing
+    the density of union and intersection of some subsets
+    """
+
+    zero: Density = None
+    one: Density = None
+
+    def __init__(self, subset: SubSetR1):
+        if not Is.instance(subset, SubSetR1):
+            raise TypeError
+        if subset not in create_interval(0, 1):
+            raise ValueError(f"{subset} not in [0, 1]")
+        self.subset = subset
+
+    def __float__(self) -> float:
+        value = subset_length(self.subset)
+        return float(value)
+
+    def __or__(self, value: Density):
+        if not Is.instance(value, Density):
+            raise TypeError
+        return Density(self.subset | value.subset)
+
+    def __and__(self, value: Density):
+        if not Is.instance(value, Density):
+            raise TypeError
+        return Density(self.subset & value.subset)
+
+    def __str__(self):  # pragma: no cover
+        return str(float(self))
+
+    def __repr__(self):  # pragma: no cover
+        return "D" + str(self)
+
+    @debug("shapepy.bool2d.density")
+    def __eq__(self, value: Union[Real, Density]) -> bool:
+        return abs(float(self) - float(value)) < 1e-6
+
+
+Density.zero = Density(EmptyR1())
+Density.one = Density(create_interval(0, 1))
+
+
+@debug("shapepy.bool2d.density")
+def unite_densities(densities: Iterable[Density]) -> Density:
+    """Computes the union of Density units"""
+    densities = iter(densities)
+    result = next(densities)
+    if not Is.instance(result, Density):
+        raise TypeError(f"Invalid {type(result)}")
+    for density in densities:
+        if not Is.instance(density, Density):
+            raise TypeError(f"Invalid {type(density)}")
+        result |= density
+    return result
+
+
+@debug("shapepy.bool2d.density")
+def intersect_densities(densities: Iterable[Density]) -> Density:
+    """Computes the intersection of Density units"""
+    densities = iter(densities)
+    result = next(densities)
+    if not Is.instance(result, Density):
+        raise TypeError(f"Invalid {type(result)}")
+    for density in densities:
+        if not Is.instance(density, Density):
+            raise TypeError(f"Invalid {type(density)}")
+        result &= density
+    return result

src/shapepy/bool2d/shape.py

@@ -13,13 +13,18 @@
 from typing import Iterable, Set, Tuple, Union
 
 from ..geometry.box import Box
-from ..geometry.integral import lebesgue_density_jordan
 from ..geometry.jordancurve import JordanCurve
 from ..geometry.point import Point2D
 from ..scalar.angle import Angle
 from ..scalar.reals import Real
-from ..tools import Is, To, prod
+from ..tools import Is, To
 from .base import EmptyShape, SubSetR2
+from .density import (
+    Density,
+    intersect_densities,
+    lebesgue_density_jordan,
+    unite_densities,
+)
 
 
 class SimpleShape(SubSetR2):
@@ -141,8 +146,7 @@ def __contains__(self, other: SubSetR2) -> bool:
         return self.__contains_point(other)
 
     def __contains_point(self, point: Point2D) -> bool:
-        density = self.density(point)
-
+        density = float(self.density(point))
         return density > 0 if self.boundary else density == 1
 
     def __contains_jordan(self, jordan: JordanCurve) -> bool:
@@ -210,7 +214,7 @@ def box(self) -> Box:
         """
         return self.jordan.box()
 
-    def density(self, center: Point2D) -> Real:
+    def density(self, center: Point2D) -> Density:
         return lebesgue_density_jordan(self.jordan, center)
 
 
@@ -344,9 +348,10 @@ def box(self) -> Box:
             box |= sub.jordan.box()
         return box
 
-    def density(self, center: Point2D) -> Real:
+    def density(self, center: Point2D) -> Density:
         center = To.point(center)
-        return prod(sub.density(center) for sub in self.subshapes)
+        densities = (sub.density(center) for sub in self.subshapes)
+        return intersect_densities(densities)
 
 
 class DisjointShape(SubSetR2):
@@ -493,8 +498,8 @@ def box(self) -> Box:
 
     def density(self, center: Point2D) -> Real:
         center = To.point(center)
-        result = sum(sub.density(center) for sub in self.subshapes)
-        return min(result, 1)
+        densities = (sub.density(center) for sub in self.subshapes)
+        return unite_densities(densities)
 
 
 def divide_connecteds(

src/shapepy/common.py

@@ -4,6 +4,7 @@
 from typing import Any, Tuple
 
 from .bool2d.base import SubSetR2
+from .bool2d.density import Density
 from .scalar.angle import Angle
 from .scalar.reals import Real
 
@@ -91,7 +92,7 @@ def derivate(obj: Any) -> Any:
     return deepcopy(obj).derivate()
 
 
-def lebesgue_density(subset: SubSetR2, center: Tuple[Real, Real]) -> Real:
+def lebesgue_density(subset: SubSetR2, center: Tuple[Real, Real]) -> Density:
     """
     Calcules the density of given subset around given point
     """

src/shapepy/geometry/concatenate.py

@@ -64,7 +64,6 @@ def concatenate_segments(
             filtsegments.append(segmenti)
             segmenti = segmentj
     filtsegments.append(segmenti)
-    print("filtsegments = ", filtsegments)
     return (
         PiecewiseCurve(filtsegments)
         if len(filtsegments) > 1