Feat: add computation of Lebesgue density for all SubSetR2

docs/source/rst/theory.rst

@@ -36,7 +36,7 @@ and how they interact with each other.
     * :ref:`theory_bool_xor`
 * :ref:`theory_generalities`
     * :ref:`onedimen_integration`
-    * :ref:`winding_function`
+    * :ref:`lebesgue_density`
 
 -----------------------------------------------------------------------------
 
@@ -480,7 +480,7 @@ If this projected point is equal to the point itself, then the point is on the c
 Point in Shape
 --------------
 
-The :ref:`winding_function` is used to determine if the shape contains the point.
+The :ref:`lebesgue_density` is used to determine if the shape contains the point.
 
 Basically this function tells if a point is inside the shape, or outside, or at the boundary:
 
@@ -1124,12 +1124,12 @@ There are available schemas are bellow, with some nodes/weights depending on :ma
           - 1/2
           - 1
 
-.. _winding_function:
+.. _lebesgue_density:
 
-Winding function
-----------------
+Lebesgue Density function
+-------------------------
 
-The **Winding function** is a function on the plane, based on the a shape :math:`S`, that
+The **Density function** is a function on the plane, based on the a shape :math:`S`, that
 
 * Is equal to :math:`1` for interior points
 * Is equal to :math:`0` for exterior points

src/shapepy/__init__.py

@@ -10,7 +10,7 @@
 from .bool2d.base import EmptyShape, WholeShape
 from .bool2d.primitive import Primitive
 from .bool2d.shape import ConnectedShape, DisjointShape, SimpleShape
-from .common import move, rotate, scale
+from .common import lebesgue_density, move, rotate, scale
 from .geometry.integral import IntegrateJordan
 from .geometry.jordancurve import JordanCurve
 from .geometry.point import Point2D

src/shapepy/analytic/base.py

@@ -40,7 +40,7 @@ def pow_keys(exp: int) -> Set[int]:
     return pow_keys(exp // 2) | pow_keys(exp - exp // 2) | {exp}
 
 
-# pylint: disable=too-few-public-methods
+# pylint: disable=duplicate-code
 class IAnalytic(ABC):
     """
     Interface Class for Analytic classes

src/shapepy/bool2d/base.py

@@ -148,6 +148,33 @@ def rotate(self, angle: Angle) -> SubSetR2:
         """
         raise NotImplementedError
 
+    @abstractmethod
+    def density(self, center: Point2D) -> Real:
+        """
+        Computes the density of the subset around given point
+
+        Parameters
+        ----------
+
+        center : Point2D
+            The position to measure the density
+
+        :return: The density in the interval [0, 1]
+        :rtype: Real
+
+        Example use
+        -----------
+        >>> from shapepy import Primitive
+        >>> circle = Primitive.circle(radius=1)
+        >>> circle.density((0, 0))
+        1
+        >>> circle.density((5, 0))
+        0
+        >>> circle.density((1, 0))
+        0.5
+        """
+        raise NotImplementedError
+
 
 class EmptyShape(SubSetR2):
     """EmptyShape is a singleton class to represent an empty shape
@@ -210,6 +237,9 @@ def scale(self, _):
     def rotate(self, _):
         return self
 
+    def density(self, center: Point2D) -> Real:
+        return 0
+
 
 class WholeShape(SubSetR2):
     """WholeShape is a singleton class to represent all plane
@@ -272,7 +302,11 @@ def scale(self, _):
     def rotate(self, _):
         return self
 
+    def density(self, center: Point2D) -> Real:
+        return 1
+
 
+# pylint: disable=duplicate-code
 class Future:
     """
     Class that stores methods that are further defined.

src/shapepy/bool2d/boolean.py

@@ -223,8 +223,8 @@ def midpoints_one_shape(
         for i, jordan in enumerate(shapea.jordans):
             for j, segment in enumerate(jordan.piecewise):
                 mid_point = segment(Fraction(1, 2))
-                wind = shapeb.winding(mid_point)
-                mid_point_in = (wind > 0 and closed) or wind == 1
+                density = shapeb.density(mid_point)
+                mid_point_in = (density > 0 and closed) or density == 1
                 if not inside ^ mid_point_in:
                     yield (i, j)
 

src/shapepy/bool2d/shape.py

@@ -13,12 +13,12 @@
 from typing import Iterable, Set, Tuple, Union
 
 from ..geometry.box import Box
-from ..geometry.integral import winding_number
+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
+from ..tools import Is, To, prod
 from .base import EmptyShape, SubSetR2
 
 
@@ -141,9 +141,9 @@ def __contains__(self, other: SubSetR2) -> bool:
         return self.__contains_point(other)
 
     def __contains_point(self, point: Point2D) -> bool:
-        wind = self.winding(point)
+        density = self.density(point)
 
-        return wind > 0 if self.boundary else wind == 1
+        return density > 0 if self.boundary else density == 1
 
     def __contains_jordan(self, jordan: JordanCurve) -> bool:
         piecewise = jordan.parametrize()
@@ -210,15 +210,8 @@ def box(self) -> Box:
         """
         return self.jordan.box()
 
-    def winding(self, point: Point2D) -> Real:
-        """Gives the winding number.
-
-        0 means the point is outside the domain
-        1 means the point is inside the domain
-        between 0 and 1 means its on the boundary"""
-        point = To.point(point)
-        wind = winding_number(self.jordan, center=point)
-        return wind
+    def density(self, center: Point2D) -> Real:
+        return lebesgue_density_jordan(self.jordan, center)
 
 
 class ConnectedShape(SubSetR2):
@@ -351,17 +344,9 @@ def box(self) -> Box:
             box |= sub.jordan.box()
         return box
 
-    def winding(self, point: Point2D) -> Real:
-        """Gives the winding number.
-
-        0 means the point is outside the domain
-        1 means the point is inside the domain
-        between 0 and 1 means its on the boundary"""
-        point = To.point(point)
-        wind = 1
-        for subset in self.subshapes:
-            wind *= subset.winding(point)
-        return wind
+    def density(self, center: Point2D) -> Real:
+        center = To.point(center)
+        return prod(sub.density(center) for sub in self.subshapes)
 
 
 class DisjointShape(SubSetR2):
@@ -506,18 +491,10 @@ def box(self) -> Box:
             box |= sub.box()
         return box
 
-    def winding(self, point: Point2D) -> Real:
-        """Gives the winding number.
-
-        0 means the point is outside the domain
-        1 means the point is inside the domain
-        between 0 and 1 means its on the boundary"""
-        point = To.point(point)
-        for subset in self.subshapes:
-            wind = subset.winding(point)
-            if wind > 0:
-                return wind
-        return 0
+    def density(self, center: Point2D) -> Real:
+        center = To.point(center)
+        result = sum(sub.density(center) for sub in self.subshapes)
+        return min(result, 1)
 
 
 def divide_connecteds(

src/shapepy/common.py

@@ -3,6 +3,7 @@
 from copy import deepcopy
 from typing import Any, Tuple
 
+from .bool2d.base import SubSetR2
 from .scalar.angle import Angle
 from .scalar.reals import Real
 
@@ -88,3 +89,10 @@ def derivate(obj: Any) -> Any:
     Derivates the analytic function or the curve
     """
     return deepcopy(obj).derivate()
+
+
+def lebesgue_density(subset: SubSetR2, center: Tuple[Real, Real]) -> Real:
+    """
+    Calcules the density of given subset around given point
+    """
+    return subset.density(center)

src/shapepy/geometry/integral.py

@@ -66,22 +66,25 @@ def polynomial(jordan: JordanCurve, expx: int, expy: int):
 
 # pylint: disable=too-many-locals
 @debug("shapepy.geometry.integral")
-def winding_number(
-    jordan: JordanCurve, center: Optional[Point2D] = (0.0, 0.0)
+def lebesgue_density_jordan(
+    jordan: JordanCurve, point: Optional[Point2D] = (0.0, 0.0)
 ) -> Union[int, float]:
-    """Computes the winding number from jordan curve
+    """Computes the lebesgue density number from jordan curve
 
-    Returns [-1, -0.5, 0, 0.5 or 1]
+    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
     """
-    center = To.point(center)
+    point = To.point(point)
     box = jordan.box()
-    if center not in box:
-        wind = 0 if jordan.area > 0 else 1
-        return wind
+    if point not in box:
+        density = 0 if jordan.area > 0 else 1
+        return density
 
     segments = tuple(jordan.parametrize())
     for i, segmenti in enumerate(segments):
-        if center == segmenti(0):
+        if point == segmenti(0):
             segmentj = segments[(i - 1) % len(segments)]
             deltapi = segmenti(0, 1)
             deltapj = segmentj(1, 1)
@@ -94,8 +97,8 @@ def winding_number(
     integrator = AdaptativeIntegrator(direct, 1e-6)
     radangle = 0
     for segment in segments:
-        deltax: IAnalytic = segment.xfunc - center.xcoord
-        deltay: IAnalytic = segment.yfunc - center.ycoord
+        deltax: IAnalytic = segment.xfunc - point.xcoord
+        deltay: IAnalytic = segment.yfunc - point.ycoord
         radius_square = deltax * deltax + deltay * deltay
         if find_minimum(radius_square, [0, 1]) < 1e-6:
             return 0.5
@@ -104,5 +107,5 @@ def winding_number(
             lambda t, cf, rs: cf(t) / rs(t), cf=crossf, rs=radius_square
         )
         radangle += integrator.integrate(function, [0, 1])
-    wind = round(radangle / Math.tau)
-    return wind if jordan.area > 0 else 1 + wind
+    density = round(radangle / Math.tau)
+    return density if jordan.area > 0 else 1 + density

src/shapepy/tools.py

@@ -108,6 +108,21 @@ def wrapper(*args, **kwargs):
     return decorator
 
 
+def prod(values: Iterable[Any]) -> Any:
+    """Computes the product of given objects
+
+    Example
+    -------
+    >>> prod([3, 2, 5])
+    30
+    """
+    values = iter(values)
+    result = next(values)
+    for value in values:
+        result *= value
+    return result
+
+
 def reverse(objs: Iterable[Any]) -> Iterable[Any]:
     """Reverts the list/tuple"""
     return tuple(objs)[::-1]