diff --git a/src/shapepy/__init__.py b/src/shapepy/__init__.py
index defee46b..53e63dbd 100644
--- a/src/shapepy/__init__.py
+++ b/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__":
diff --git a/src/shapepy/bool2d/base.py b/src/shapepy/bool2d/base.py
index 315f0d3d..53ca0feb 100644
--- a/src/shapepy/bool2d/base.py
+++ b/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
diff --git a/src/shapepy/bool2d/boolean.py b/src/shapepy/bool2d/boolean.py
index afd939c0..d4f0bc36 100644
--- a/src/shapepy/bool2d/boolean.py
+++ b/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)
 
diff --git a/src/shapepy/bool2d/density.py b/src/shapepy/bool2d/density.py
new file mode 100644
index 00000000..2a52a47e
--- /dev/null
+++ b/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
diff --git a/src/shapepy/bool2d/shape.py b/src/shapepy/bool2d/shape.py
index b7492489..d750c40e 100644
--- a/src/shapepy/bool2d/shape.py
+++ b/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(
diff --git a/src/shapepy/common.py b/src/shapepy/common.py
index 77e102d9..e91c3895 100644
--- a/src/shapepy/common.py
+++ b/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
     """
diff --git a/src/shapepy/geometry/concatenate.py b/src/shapepy/geometry/concatenate.py
index 72cb20b5..e4310bde 100644
--- a/src/shapepy/geometry/concatenate.py
+++ b/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
diff --git a/src/shapepy/geometry/integral.py b/src/shapepy/geometry/integral.py
index 7827d973..8ed81ba9 100644
--- a/src/shapepy/geometry/integral.py
+++ b/src/shapepy/geometry/integral.py
@@ -4,19 +4,16 @@
 
 from __future__ import annotations
 
+from copy import copy
 from functools import partial
-from typing import Optional, Union
 
 from ..analytic.base import IAnalytic
 from ..analytic.tools import find_minimum
-from ..common import derivate
-from ..loggers import debug
-from ..scalar.angle import arg
 from ..scalar.quadrature import AdaptativeIntegrator, IntegratorFactory
 from ..scalar.reals import Math
 from ..tools import Is, To
 from .jordancurve import JordanCurve
-from .point import Point2D, cross, inner
+from .point import Point2D
 from .segment import Segment
 
 
@@ -26,6 +23,9 @@ class IntegrateSegment:
     Defines methods to integrates over a segment
     """
 
+    direct = IntegratorFactory.closed_newton_cotes(3)
+    adaptative = AdaptativeIntegrator(direct, 1e-6)
+
     @staticmethod
     def polynomial(curve: Segment, expx: int, expy: int):
         """
@@ -43,6 +43,33 @@ def polynomial(curve: Segment, expx: int, expy: int):
         ipoly = function.integrate()
         return (ipoly(1) - ipoly(0)) / (expx + expy + 2)
 
+    @staticmethod
+    def turns(curve: Segment, point: Point2D) -> float:
+        """
+        Computes the integral
+
+        I = int_0^1 (u * v' - v * u') / (u * u + v * v) dt
+
+        Where
+
+        u = curve.x - point.x
+        v = curve.y - point.y
+        """
+        point = To.point(point)
+        deltax: IAnalytic = curve.xfunc - point.xcoord
+        deltay: IAnalytic = curve.yfunc - point.ycoord
+        radius_square = deltax * deltax + deltay * deltay
+        if find_minimum(radius_square, [0, 1]) < 1e-6:
+            return To.rational(1, 2)
+        crossf = (
+            deltax * copy(deltay).derivate() - deltay * copy(deltax).derivate()
+        )
+        function = partial(
+            lambda t, cf, rs: cf(t) / rs(t), cf=crossf, rs=radius_square
+        )
+        radians = IntegrateSegment.adaptative.integrate(function, [0, 1])
+        return radians / Math.tau
+
 
 class IntegrateJordan:
     """
@@ -63,49 +90,27 @@ def polynomial(jordan: JordanCurve, expx: int, expy: int):
             for usegment in jordan.usegments
         )
 
+    @staticmethod
+    def turns(jordan: JordanCurve, point: Point2D) -> float:
+        """
+        Computes the integral
 
-# pylint: disable=too-many-locals
-@debug("shapepy.geometry.integral")
-def lebesgue_density_jordan(
-    jordan: JordanCurve, point: Optional[Point2D] = (0.0, 0.0)
-) -> Union[int, float]:
-    """Computes the lebesgue density number from jordan curve
+        I = int_0^1 (u * v' - v * u') / (u * u + v * v) dt
 
-    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:
-        density = 0 if jordan.area > 0 else 1
-        return density
-
-    segments = tuple(jordan.parametrize())
-    for i, segmenti in enumerate(segments):
-        if point == segmenti(0):
-            segmentj = segments[(i - 1) % len(segments)]
-            deltapi = segmenti(0, 1)
-            deltapj = segmentj(1, 1)
-            innerval = inner(deltapi, deltapj)
-            crossval = cross(deltapi, deltapj)
-            angle = arg(-innerval, -crossval)
-            return angle.turns % 1
+        Where
 
-    direct = IntegratorFactory.closed_newton_cotes(3)
-    integrator = AdaptativeIntegrator(direct, 1e-6)
-    radangle = 0
-    for segment in segments:
-        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
-        crossf = deltax * derivate(deltay) - deltay * derivate(deltax)
-        function = partial(
-            lambda t, cf, rs: cf(t) / rs(t), cf=crossf, rs=radius_square
-        )
-        radangle += integrator.integrate(function, [0, 1])
-    density = round(radangle / Math.tau)
-    return density if jordan.area > 0 else 1 + density
+        u = jordan.x - point.x
+        v = jordan.y - point.y
+
+        This function should return values other than
+        [-1, -0.5, 0, 0.5, 1]
+        It happens when the functions is discontinuous
+        """
+        result = 0
+        for usegment in jordan.usegments:
+            seg = usegment.parametrize()
+            delta_result = IntegrateSegment.turns(seg, point)
+            if delta_result == 0.5:
+                return 0.5 if jordan.area > 0 else -0.5
+            result += delta_result
+        return result
diff --git a/src/shapepy/geometry/intersection.py b/src/shapepy/geometry/intersection.py
index 21620bac..309d2fe2 100644
--- a/src/shapepy/geometry/intersection.py
+++ b/src/shapepy/geometry/intersection.py
@@ -16,9 +16,9 @@
 
 from ..rbool import (
     EmptyR1,
-    IntervalR1,
-    SingleR1,
     SubSetR1,
+    create_interval,
+    create_single,
     extract_knots,
     from_any,
 )
@@ -172,7 +172,7 @@ def __compute_two(
             return EmptyR1(), EmptyR1()
         if id(curvea) == id(curveb):  # Check if curves are equal
             curvea = curvea.parametrize()
-            subset = IntervalR1(curvea.knots[0], curvea.knots[-1])
+            subset = create_interval(curvea.knots[0], curvea.knots[-1])
             return subset, subset
         return curve_and_curve(curvea, curveb)
 
@@ -227,7 +227,7 @@ def segment_and_segment(
     if curvea.box() & curveb.box() is None:
         return EmptyR1(), EmptyR1()
     if curvea == curveb:
-        return IntervalR1(0, 1), IntervalR1(0, 1)
+        return create_interval(0, 1), create_interval(0, 1)
     if segment_is_linear(curvea) and segment_is_linear(curveb):
         return IntersectionSegments.lines(curvea, curveb)
     nptsa = max(curvea.xfunc.degree, curvea.yfunc.degree) + 4
@@ -283,7 +283,7 @@ def lines(curvea: Segment, curveb: Segment) -> Tuple[SubSetR1, SubSetR1]:
             u0 = cross(B0mA0, dA) / dAxdB
             if u0 < 0 or 1 < u0:
                 return empty, empty
-            return SingleR1(t0), SingleR1(u0)
+            return create_single(t0), create_single(u0)
         # Lines are parallel
         if cross(dA, B0mA0) != 0:
             return empty, empty  # Parallel, but not colinear
@@ -306,8 +306,8 @@ def lines(curvea: Segment, curveb: Segment) -> Tuple[SubSetR1, SubSetR1]:
         u0 = min(max(0, u0), 1)
         u1 = min(max(0, u1), 1)
         if t0 == t1 or u0 == u1:
-            return SingleR1(t0), SingleR1(u1)
-        return IntervalR1(t0, t1), IntervalR1(u0, u1)
+            return create_single(t0), create_single(u1)
+        return create_interval(t0, t1), create_interval(u0, u1)
 
     # pylint: disable=too-many-locals
     @staticmethod
diff --git a/src/shapepy/geometry/jordancurve.py b/src/shapepy/geometry/jordancurve.py
index e205d6d2..4e9f13d2 100644
--- a/src/shapepy/geometry/jordancurve.py
+++ b/src/shapepy/geometry/jordancurve.py
@@ -9,7 +9,6 @@
 from copy import copy
 from typing import Iterable, Iterator, Tuple, Union
 
-from ..common import clean
 from ..loggers import debug
 from ..scalar.angle import Angle
 from ..scalar.reals import Real
@@ -180,7 +179,7 @@ def __eq__(self, other: JordanCurve) -> bool:
             or not all(point in self for point in other.vertices())
         ):
             return False
-        return clean(self).usegments == clean(other).usegments
+        return copy(self).clean().usegments == copy(other).clean().usegments
 
     @debug("shapepy.geometry.jordancurve")
     def invert(self) -> JordanCurve:
diff --git a/src/shapepy/geometry/point.py b/src/shapepy/geometry/point.py
index b92f59f1..429cb8ab 100644
--- a/src/shapepy/geometry/point.py
+++ b/src/shapepy/geometry/point.py
@@ -117,7 +117,7 @@ def __neg__(self) -> Point2D:
             -self.xcoord,
             -self.ycoord,
             self.radius,
-            self.angle + degrees(180),
+            ~self.angle,
         )
 
     def __pos__(self) -> Point2D:
diff --git a/src/shapepy/loggers.py b/src/shapepy/loggers.py
index f1317c7e..a920ba72 100644
--- a/src/shapepy/loggers.py
+++ b/src/shapepy/loggers.py
@@ -121,6 +121,27 @@ def indent():
         IndentingLoggerAdapter.indent_level -= 1
 
 
+@contextmanager
+def enable_logger(base: str, /, *, level: logging._Level):
+    """Enables temporarily the given logger"""
+    current_enable = LogConfiguration.log_enabled
+    current_levels = {}
+    for name, logger in IndentingLoggerAdapter.instances.items():
+        if base in name:
+            current_levels[name] = logger.getEffectiveLevel()
+    try:
+        LogConfiguration.log_enabled = True
+        for name, logger in IndentingLoggerAdapter.instances.items():
+            if name in current_levels:
+                logger.setLevel(level)
+        yield
+    finally:
+        LogConfiguration.log_enabled = current_enable
+        for name, logger in IndentingLoggerAdapter.instances.items():
+            if name in current_levels:
+                logger.setLevel(current_levels[name])
+
+
 # Create decorator to use in functions
 def debug(name: Optional[str] = None, /, *, maxdepth: Optional[int] = None):
     """
diff --git a/src/shapepy/rbool.py b/src/shapepy/rbool.py
index bde18596..4cbc89a8 100644
--- a/src/shapepy/rbool.py
+++ b/src/shapepy/rbool.py
@@ -4,13 +4,57 @@
 
 import rbool
 
+from .loggers import debug
+from .scalar.reals import Real
+from .tools import Is
+
 EmptyR1: Type = rbool.Empty
 IntervalR1: Type = rbool.Interval
 SingleR1: Type = rbool.SingleValue
 SubSetR1: Type = rbool.SubSetR1
+DisjointR1: Type = rbool.Disjoint
 WholeR1: Type = rbool.Whole
 extract_knots: Callable[[Any], Iterator[Any]] = rbool.extract_knots
 from_any: Callable[[Any], object] = rbool.from_any
 shift = rbool.move
 scale = rbool.scale
 unite = rbool.unite
+infimum = rbool.infimum
+supremum = rbool.supremum
+
+
+def create_single(knot: Real) -> SingleR1:
+    """Creates a Subset on real-line that contains only one value"""
+    if not Is.finite(knot):
+        raise ValueError(f"Invalid {knot}")
+    return SingleR1(knot)
+
+
+def create_interval(knota: Real, knotb: Real) -> IntervalR1:
+    """Creates a closed interval [a, b] in the real line"""
+    if not Is.real(knota) or not Is.real(knotb):
+        raise TypeError(f"Invalid typos: {type(knota)}, {type(knotb)}")
+    if knotb <= knota:
+        raise ValueError(f"{knotb} <= {knota}")
+    return IntervalR1(knota, knotb)
+
+
+@debug("shapepy.rbool")
+def subset_length(subset: SubSetR1) -> Real:
+    """Computes the length of the subset
+
+    Example
+    -------
+    >>> subset_length([0, 1])
+    1
+    >>> subset_length([-3, 2])
+    5
+    >>> subset_length({})
+    0
+    """
+    subset = from_any(subset)
+    if Is.instance(subset, IntervalR1):
+        return subset[1] - subset[0]
+    if Is.instance(subset, DisjointR1):
+        return sum(map(subset_length, subset.intervals))
+    return 0
diff --git a/src/shapepy/scalar/angle.py b/src/shapepy/scalar/angle.py
index e2131b63..d5498f22 100644
--- a/src/shapepy/scalar/angle.py
+++ b/src/shapepy/scalar/angle.py
@@ -183,6 +183,9 @@ def __eq__(self, other: object) -> bool:
     def __float__(self):
         return float(self.radians)
 
+    def __invert__(self):
+        return Angle(self.direction + 2, self.part)
+
     def __add__(self, other: Angle) -> Angle:
         other: Angle = To.angle(other)
         return turns(self.turns + other.turns)
@@ -308,6 +311,8 @@ def turns(self) -> Real:
         >>> degrees(360).turns
         0
         """
+        if self.direction == 0 and self.part < 0:
+            return To.rational(1) + self.part
         return self.part + To.rational(self.direction, 4)
 
     def sin(self) -> Real:
diff --git a/src/shapepy/tools.py b/src/shapepy/tools.py
index f4441daf..621dfc38 100644
--- a/src/shapepy/tools.py
+++ b/src/shapepy/tools.py
@@ -8,7 +8,7 @@
 import types
 from collections import deque
 from functools import wraps
-from typing import Any, Generic, Iterable, Tuple, TypeVar
+from typing import Any, Generic, Iterable, Iterator, Tuple, TypeVar
 
 import numpy as np
 
@@ -108,21 +108,6 @@ 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]
@@ -167,7 +152,7 @@ class CyclicContainer(Generic[T]):
     def __init__(self, values: Iterable[T]):
         self.__values = tuple(values)
 
-    def __iter__(self) -> Iterable[T]:
+    def __iter__(self) -> Iterator[T]:
         yield from self.__values
 
     def __getitem__(self, index):
diff --git a/tests/bool2d/test_contains.py b/tests/bool2d/test_contains.py
index 04194926..f28a4c9b 100644
--- a/tests/bool2d/test_contains.py
+++ b/tests/bool2d/test_contains.py
@@ -13,7 +13,7 @@
 from shapepy.geometry.factory import FactoryJordan
 
 
-@pytest.mark.order(23)
+@pytest.mark.order(24)
 @pytest.mark.dependency(
     depends=[
         "tests/geometry/test_integral.py::test_all",
@@ -31,12 +31,12 @@ class TestObjectsInEmptyWhole:
     Test relative to special cases, a empty shape and whole domain
     """
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["test_begin"])
     def test_begin(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInEmptyWhole::test_begin"])
     def test_singleton(self):
         empty0 = EmptyShape()
@@ -49,7 +49,7 @@ def test_singleton(self):
         assert whole0 is whole1
         assert empty1 != whole1
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -62,7 +62,7 @@ def test_empty(self):
         assert empty in empty
         assert empty in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -75,7 +75,7 @@ def test_whole(self):
         assert whole not in empty
         assert whole in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -89,7 +89,7 @@ def test_point(self):
             assert point not in empty
             assert point in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -110,7 +110,7 @@ def test_jordan(self):
         assert jordan not in empty
         assert jordan in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -125,7 +125,7 @@ def test_simple_shape(self):
         assert shape not in empty
         assert shape in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -142,7 +142,7 @@ def test_connected_shape(self):
         assert shape not in empty
         assert shape in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -159,7 +159,7 @@ def test_disjoint_shape(self):
         assert shape not in empty
         assert shape in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInEmptyWhole::test_begin",
@@ -182,7 +182,7 @@ class TestObjectsInJordan:
     Tests the respective position
     """
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "test_begin",
@@ -192,7 +192,7 @@ class TestObjectsInJordan:
     def test_begin(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.timeout(10)
     @pytest.mark.dependency(
         depends=[
@@ -221,7 +221,7 @@ def test_boundary_point(self):
         assert (0.5, 1) in square
         assert (0, 0.5) in square
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.timeout(10)
     @pytest.mark.dependency(
         depends=[
@@ -238,7 +238,7 @@ def test_interior_point(self):
         square = FactoryJordan.polygon(vertices)
         assert (1, 1) not in square
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.timeout(10)
     @pytest.mark.dependency(
         depends=[
@@ -255,7 +255,7 @@ def test_exterior_point(self):
         assert (-1, -1) not in square
         assert (2, 2) not in square
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.timeout(10)
     @pytest.mark.dependency(
         depends=[
@@ -274,7 +274,7 @@ class TestObjectsInSimple:
     Tests the respective position
     """
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "test_begin",
@@ -285,7 +285,7 @@ class TestObjectsInSimple:
     def test_begin(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInSimple::test_begin"])
     def test_empty(self):
         empty = EmptyShape()
@@ -293,7 +293,7 @@ def test_empty(self):
         assert empty in square
         assert square not in empty
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInSimple::test_begin"])
     def test_whole(self):
         whole = WholeShape()
@@ -301,7 +301,7 @@ def test_whole(self):
         assert whole not in square
         assert square in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -325,7 +325,7 @@ def test_keep_type(self):
             assert len(test_types) == len(good_types)
             assert test_types == good_types
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -358,7 +358,7 @@ def test_point(self):
         assert (3, 3) not in square
         assert (-3, -3) not in square
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -392,7 +392,7 @@ def test_jordan(self):
         assert ~(big_square.jordan) in (~small_square)
         assert ~(big_square.jordan) in (~big_square)
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -428,7 +428,7 @@ def test_simple(self):
         assert (~left) not in (~right)
         assert (~right) not in (~left)
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -450,7 +450,7 @@ def test_connected(self):
         assert connected in (~small_square)
         assert connected not in (~big_square)
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -486,7 +486,7 @@ def test_disjoint(self):
         square = Primitive.square(side=7, center=(0, 0))
         assert disj_shape not in square
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInSimple::test_begin",
@@ -508,7 +508,7 @@ class TestObjectsInConnected:
     Tests the respective position
     """
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "test_begin",
@@ -520,7 +520,7 @@ class TestObjectsInConnected:
     def test_begin(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInConnected::test_begin"])
     def test_empty(self):
         empty = EmptyShape()
@@ -528,7 +528,7 @@ def test_empty(self):
         assert empty in square
         assert square not in empty
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInConnected::test_begin"])
     def test_whole(self):
         whole = WholeShape()
@@ -536,7 +536,7 @@ def test_whole(self):
         assert whole not in square
         assert square in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInConnected::test_begin",
@@ -570,7 +570,7 @@ def test_point(self):
         assert (1.5, 1.5) in connected
         assert (1.5, -1.5) in connected
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInConnected::test_begin",
@@ -596,7 +596,7 @@ def test_jordan(self):
         assert (~big_jordan) in connected
         assert (~big_jordan) in connected
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInConnected::test_begin",
@@ -616,7 +616,7 @@ def test_simple(self):
         assert (~small_square) not in connected
         assert (~big_square) not in connected
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInConnected::test_begin",
@@ -637,7 +637,7 @@ def test_connected(self):
         assert (~small_square) not in (~big_square)
         assert (~big_square) in (~big_square)
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInConnected::test_begin",
@@ -652,7 +652,7 @@ def test_connected(self):
     def test_disjoint(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInConnected::test_begin",
@@ -674,7 +674,7 @@ class TestObjectsInDisjoint:
     Tests the respective position
     """
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "test_begin",
@@ -687,7 +687,7 @@ class TestObjectsInDisjoint:
     def test_begin(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInDisjoint::test_begin"])
     def test_empty(self):
         empty = EmptyShape()
@@ -698,7 +698,7 @@ def test_empty(self):
         assert empty in disj_shape
         assert disj_shape not in empty
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(depends=["TestObjectsInDisjoint::test_begin"])
     def test_whole(self):
         whole = WholeShape()
@@ -709,7 +709,7 @@ def test_whole(self):
         assert whole not in disj_shape
         assert disj_shape in whole
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInDisjoint::test_begin",
@@ -727,7 +727,7 @@ def test_point(self):
         assert (-3, 0) in disj_shape
         assert (3, 0) in disj_shape
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInDisjoint::test_begin",
@@ -739,7 +739,7 @@ def test_point(self):
     def test_jordan(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInDisjoint::test_begin",
@@ -752,7 +752,7 @@ def test_jordan(self):
     def test_simple(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInDisjoint::test_begin",
@@ -766,7 +766,7 @@ def test_simple(self):
     def test_connected(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInDisjoint::test_begin",
@@ -781,7 +781,7 @@ def test_connected(self):
     def test_disjoint(self):
         pass
 
-    @pytest.mark.order(23)
+    @pytest.mark.order(24)
     @pytest.mark.dependency(
         depends=[
             "TestObjectsInDisjoint::test_begin",
@@ -798,7 +798,7 @@ def test_end(self):
         pass
 
 
-@pytest.mark.order(23)
+@pytest.mark.order(24)
 @pytest.mark.dependency(
     depends=[
         "TestObjectsInEmptyWhole::test_end",
diff --git a/tests/bool2d/test_density.py b/tests/bool2d/test_density.py
index 42cef717..2bc0b577 100644
--- a/tests/bool2d/test_density.py
+++ b/tests/bool2d/test_density.py
@@ -3,21 +3,27 @@
 of an object with respect to another
 """
 
+from fractions import Fraction as frac
+
 import pytest
 
 from shapepy import lebesgue_density
 from shapepy.bool2d.base import EmptyShape, WholeShape
+from shapepy.bool2d.density import lebesgue_density_jordan
 from shapepy.bool2d.primitive import Primitive
 from shapepy.bool2d.shape import ConnectedShape, DisjointShape
+from shapepy.geometry.factory import FactoryJordan
 from shapepy.geometry.point import polar
-from shapepy.scalar.angle import degrees
+from shapepy.scalar.angle import degrees, turns
 
 
-@pytest.mark.order(22)
+@pytest.mark.order(23)
 @pytest.mark.dependency(
     depends=[
         "tests/geometry/test_integral.py::test_all",
+        "tests/geometry/test_jordan_polygon.py::test_all",
         "tests/bool2d/test_empty_whole.py::test_end",
+        "tests/bool2d/test_primitive.py::test_end",
     ],
     scope="session",
 )
@@ -25,8 +31,230 @@ def test_begin():
     pass
 
 
-@pytest.mark.order(22)
-@pytest.mark.dependency(depends=["test_begin"])
+class TestJordan:
+    """
+    Tests the respective position
+    """
+
+    @pytest.mark.order(23)
+    @pytest.mark.dependency(depends=["test_begin"])
+    def test_begin(self):
+        pass
+
+    @pytest.mark.order(23)
+    @pytest.mark.dependency(
+        depends=[
+            "TestJordan::test_begin",
+        ]
+    )
+    def test_standard_square(self):
+        vertices = [(-1, -1), (1, -1), (1, 1), (-1, 1)]
+        jordan = FactoryJordan.polygon(vertices)
+        assert jordan.length == 8
+        assert jordan.area == 4
+
+        interiors = {
+            (0, 0),
+            (0.5, 0.5),
+            (-0.5, 0.5),
+            (-0.5, -0.5),
+            (0.5, -0.5),
+        }
+        exteriors = {
+            (2, 1),
+            (2, 2),
+            (3, 3),
+            (-1, -2),
+        }
+        mid_edges = {
+            (1, 0),
+            (0, 1),
+            (-1, 0),
+            (0, -1),
+        }
+        corners = {
+            (-1, -1): 0.25,
+            (1, -1): 0.25,
+            (1, 1): 0.25,
+            (-1, 1): 0.25,
+        }
+        for point in interiors:
+            assert lebesgue_density_jordan(jordan, point) == 1
+        for point in exteriors:
+            assert lebesgue_density_jordan(jordan, point) == 0
+        for point in mid_edges:
+            assert lebesgue_density_jordan(jordan, point) == 0.5
+        for point, density in corners.items():
+            assert lebesgue_density_jordan(jordan, point) == density
+
+    @pytest.mark.order(23)
+    @pytest.mark.dependency(
+        depends=[
+            "TestJordan::test_begin",
+            "TestJordan::test_standard_square",
+        ]
+    )
+    def test_inverted_square(self):
+        vertices = [(-1, -1), (-1, 1), (1, 1), (1, -1)]
+        jordan = FactoryJordan.polygon(vertices)
+        assert jordan.length == 8
+        assert jordan.area == -4
+
+        interiors = {
+            (2, 1),
+            (2, 2),
+            (3, 3),
+            (-1, -2),
+        }
+        exteriors = {(0, 0)}
+        mid_edges = {
+            (1, 0),
+            (0, 1),
+            (-1, 0),
+            (0, -1),
+        }
+        corners = {
+            (-1, -1): 0.75,
+            (1, -1): 0.75,
+            (1, 1): 0.75,
+            (-1, 1): 0.75,
+        }
+        for point in interiors:
+            assert lebesgue_density_jordan(jordan, point) == 1
+        for point in exteriors:
+            assert lebesgue_density_jordan(jordan, point) == 0
+        for point in mid_edges:
+            assert lebesgue_density_jordan(jordan, point) == 0.5
+        for point, density in corners.items():
+            assert lebesgue_density_jordan(jordan, point) == density
+
+    @pytest.mark.order(23)
+    @pytest.mark.dependency(
+        depends=[
+            "TestJordan::test_begin",
+            "TestJordan::test_standard_square",
+            "TestJordan::test_inverted_square",
+        ]
+    )
+    def test_standard_triangle(self):
+        vertices = [(0, 0), (3, 0), (0, 3)]
+        jordan = FactoryJordan.polygon(vertices)
+        assert jordan.area == 9 / 2
+
+        interiors = {(1, 1)}
+        exteriors = {
+            (2, 2),
+            (3, 3),
+            (-1, -1),
+        }
+        mid_edges = {
+            (1, 0),
+            (2, 0),
+            (2, 1),
+            (1, 2),
+            (0, 2),
+            (0, 1),
+        }
+        corners = {
+            (0, 0): 0.25,
+            (3, 0): 0.125,
+            (0, 3): 0.125,
+        }
+        for point in interiors:
+            assert lebesgue_density_jordan(jordan, point) == 1
+        for point in exteriors:
+            assert lebesgue_density_jordan(jordan, point) == 0
+        for point in mid_edges:
+            assert lebesgue_density_jordan(jordan, point) == 0.5
+        for point, density in corners.items():
+            assert lebesgue_density_jordan(jordan, point) == density
+
+    @pytest.mark.order(23)
+    @pytest.mark.dependency(
+        depends=[
+            "TestJordan::test_begin",
+            "TestJordan::test_inverted_square",
+            "TestJordan::test_standard_triangle",
+        ]
+    )
+    def test_inverted_triangle(self):
+        vertices = [(0, 0), (0, 3), (3, 0)]
+        jordan = FactoryJordan.polygon(vertices)
+        assert jordan.area == -9 / 2
+
+        interiors = {
+            (2, 2),
+            (3, 3),
+            (-1, -1),
+        }
+        exteriors = {(1, 1)}
+        mid_edges = {
+            (1, 0),
+            (2, 0),
+            (2, 1),
+            (1, 2),
+            (0, 2),
+            (0, 1),
+        }
+        corners = {
+            (0, 0): 0.75,
+            (3, 0): 0.875,
+            (0, 3): 0.875,
+        }
+        for point in interiors:
+            assert lebesgue_density_jordan(jordan, point) == 1
+        for point in exteriors:
+            assert lebesgue_density_jordan(jordan, point) == 0
+        for point in mid_edges:
+            assert lebesgue_density_jordan(jordan, point) == 0.5
+        for point, density in corners.items():
+            assert lebesgue_density_jordan(jordan, point) == density
+
+    @pytest.mark.order(23)
+    @pytest.mark.timeout(10)
+    @pytest.mark.dependency(
+        depends=[
+            "TestJordan::test_begin",
+            "TestJordan::test_standard_square",
+            "TestJordan::test_inverted_square",
+            "TestJordan::test_standard_triangle",
+            "TestJordan::test_inverted_triangle",
+        ]
+    )
+    def test_regular_polygon(self):
+        for nsides in range(3, 10):
+            angles = (turns(frac(i, nsides)) for i in range(nsides))
+            ctrlpoints = tuple((a.cos(), a.sin()) for a in angles)
+
+            # Counter clockwise
+            vertices = tuple(ctrlpoints)
+            jordancurve = FactoryJordan.polygon(vertices)
+            density = lebesgue_density_jordan(jordancurve)
+            assert density == 1
+
+            # Clockwise
+            vertices = tuple(ctrlpoints[::-1])
+            jordancurve = FactoryJordan.polygon(vertices)
+            density = lebesgue_density_jordan(jordancurve)
+            assert density == 0
+
+    @pytest.mark.order(23)
+    @pytest.mark.dependency(
+        depends=[
+            "TestJordan::test_begin",
+            "TestJordan::test_standard_square",
+            "TestJordan::test_inverted_square",
+            "TestJordan::test_standard_triangle",
+            "TestJordan::test_inverted_triangle",
+            "TestJordan::test_regular_polygon",
+        ]
+    )
+    def test_all(self):
+        pass
+
+
+@pytest.mark.order(23)
+@pytest.mark.dependency(depends=["test_begin", "TestJordan::test_all"])
 def test_empty_whole():
     empty = EmptyShape()
     whole = WholeShape()
@@ -45,8 +273,10 @@ def test_empty_whole():
         assert whole.density(point) == 1
 
 
-@pytest.mark.order(22)
-@pytest.mark.dependency(depends=["test_begin", "test_empty_whole"])
+@pytest.mark.order(23)
+@pytest.mark.dependency(
+    depends=["test_begin", "test_empty_whole", "TestJordan::test_all"]
+)
 def test_simple_shape():
     shape = Primitive.triangle(3)
     # Corners
@@ -82,12 +312,13 @@ def test_simple_shape():
         assert shape.density(point) == value
 
 
-@pytest.mark.order(22)
+@pytest.mark.order(23)
 @pytest.mark.dependency(
     depends=[
         "test_begin",
         "test_empty_whole",
         "test_simple_shape",
+        "TestJordan::test_all",
     ]
 )
 def test_connected_shape():
@@ -139,12 +370,13 @@ def test_connected_shape():
         assert shape.density(point) == value
 
 
-@pytest.mark.order(22)
+@pytest.mark.order(23)
 @pytest.mark.dependency(
     depends=[
         "test_begin",
         "test_simple_shape",
         "test_connected_shape",
+        "TestJordan::test_all",
     ]
 )
 def test_disjoint_shape():
@@ -186,10 +418,11 @@ def test_disjoint_shape():
         assert shape.density(point) == value
 
 
-@pytest.mark.order(22)
+@pytest.mark.order(23)
 @pytest.mark.dependency(
     depends=[
         "test_begin",
+        "TestJordan::test_all",
         "test_empty_whole",
         "test_simple_shape",
         "test_connected_shape",
diff --git a/tests/geometry/test_integral.py b/tests/geometry/test_integral.py
index f47600f5..f0a27957 100644
--- a/tests/geometry/test_integral.py
+++ b/tests/geometry/test_integral.py
@@ -3,7 +3,6 @@
 import pytest
 
 from shapepy.geometry.factory import FactoryJordan, FactorySegment
-from shapepy.geometry.integral import lebesgue_density_jordan
 
 
 @pytest.mark.order(15)
@@ -64,193 +63,6 @@ def test_area():
     assert jordan.area == -6
 
 
-class TestDensity:
-    """
-    Tests the respective position
-    """
-
-    @pytest.mark.order(15)
-    @pytest.mark.dependency(depends=["test_begin", "test_area"])
-    def test_begin(self):
-        pass
-
-    @pytest.mark.order(15)
-    @pytest.mark.dependency(
-        depends=[
-            "TestDensity::test_begin",
-        ]
-    )
-    def test_standard_square(self):
-        vertices = [(-1, -1), (1, -1), (1, 1), (-1, 1)]
-        jordan = FactoryJordan.polygon(vertices)
-        assert jordan.length == 8
-        assert jordan.area == 4
-
-        interiors = {(0, 0)}
-        exteriors = {
-            (2, 1),
-            (2, 2),
-            (3, 3),
-            (-1, -2),
-        }
-        mid_edges = {
-            (1, 0),
-            (0, 1),
-            (-1, 0),
-            (0, -1),
-        }
-        corners = {
-            (-1, -1): 0.25,
-            (1, -1): 0.25,
-            (1, 1): 0.25,
-            (-1, 1): 0.25,
-        }
-        for point in interiors:
-            assert lebesgue_density_jordan(jordan, point) == 1
-        for point in exteriors:
-            assert lebesgue_density_jordan(jordan, point) == 0
-        for point in mid_edges:
-            assert lebesgue_density_jordan(jordan, point) == 0.5
-        for point, density in corners.items():
-            assert lebesgue_density_jordan(jordan, point) == density
-
-    @pytest.mark.order(15)
-    @pytest.mark.dependency(
-        depends=[
-            "TestDensity::test_begin",
-            "TestDensity::test_standard_square",
-        ]
-    )
-    def test_inverted_square(self):
-        vertices = [(-1, -1), (-1, 1), (1, 1), (1, -1)]
-        jordan = FactoryJordan.polygon(vertices)
-        assert jordan.length == 8
-        assert jordan.area == -4
-
-        interiors = {
-            (2, 1),
-            (2, 2),
-            (3, 3),
-            (-1, -2),
-        }
-        exteriors = {(0, 0)}
-        mid_edges = {
-            (1, 0),
-            (0, 1),
-            (-1, 0),
-            (0, -1),
-        }
-        corners = {
-            (-1, -1): 0.75,
-            (1, -1): 0.75,
-            (1, 1): 0.75,
-            (-1, 1): 0.75,
-        }
-        for point in interiors:
-            assert lebesgue_density_jordan(jordan, point) == 1
-        for point in exteriors:
-            assert lebesgue_density_jordan(jordan, point) == 0
-        for point in mid_edges:
-            assert lebesgue_density_jordan(jordan, point) == 0.5
-        for point, density in corners.items():
-            assert lebesgue_density_jordan(jordan, point) == density
-
-    @pytest.mark.order(15)
-    @pytest.mark.dependency(
-        depends=[
-            "TestDensity::test_begin",
-            "TestDensity::test_standard_square",
-            "TestDensity::test_inverted_square",
-        ]
-    )
-    def test_standard_triangle(self):
-        vertices = [(0, 0), (3, 0), (0, 3)]
-        jordan = FactoryJordan.polygon(vertices)
-        assert jordan.area == 9 / 2
-
-        interiors = {(1, 1)}
-        exteriors = {
-            (2, 2),
-            (3, 3),
-            (-1, -1),
-        }
-        mid_edges = {
-            (1, 0),
-            (2, 0),
-            (2, 1),
-            (1, 2),
-            (0, 2),
-            (0, 1),
-        }
-        corners = {
-            (0, 0): 0.25,
-            (3, 0): 0.125,
-            (0, 3): 0.125,
-        }
-        for point in interiors:
-            assert lebesgue_density_jordan(jordan, point) == 1
-        for point in exteriors:
-            assert lebesgue_density_jordan(jordan, point) == 0
-        for point in mid_edges:
-            assert lebesgue_density_jordan(jordan, point) == 0.5
-        for point, density in corners.items():
-            assert lebesgue_density_jordan(jordan, point) == density
-
-    @pytest.mark.order(15)
-    @pytest.mark.dependency(
-        depends=[
-            "TestDensity::test_begin",
-            "TestDensity::test_inverted_square",
-            "TestDensity::test_standard_triangle",
-        ]
-    )
-    def test_inverted_triangle(self):
-        vertices = [(0, 0), (0, 3), (3, 0)]
-        jordan = FactoryJordan.polygon(vertices)
-        assert jordan.area == -9 / 2
-
-        interiors = {
-            (2, 2),
-            (3, 3),
-            (-1, -1),
-        }
-        exteriors = {(1, 1)}
-        mid_edges = {
-            (1, 0),
-            (2, 0),
-            (2, 1),
-            (1, 2),
-            (0, 2),
-            (0, 1),
-        }
-        corners = {
-            (0, 0): 0.75,
-            (3, 0): 0.875,
-            (0, 3): 0.875,
-        }
-        for point in interiors:
-            assert lebesgue_density_jordan(jordan, point) == 1
-        for point in exteriors:
-            assert lebesgue_density_jordan(jordan, point) == 0
-        for point in mid_edges:
-            assert lebesgue_density_jordan(jordan, point) == 0.5
-        for point, density in corners.items():
-            assert lebesgue_density_jordan(jordan, point) == density
-
-    @pytest.mark.order(15)
-    @pytest.mark.dependency(
-        depends=[
-            "TestDensity::test_begin",
-            "TestDensity::test_standard_square",
-            "TestDensity::test_inverted_square",
-            "TestDensity::test_standard_triangle",
-            "TestDensity::test_inverted_triangle",
-        ]
-    )
-    def test_all(self):
-        pass
-
-
 @pytest.mark.order(15)
 @pytest.mark.timeout(10)
 @pytest.mark.dependency(
@@ -259,7 +71,6 @@ def test_all(self):
         "test_segment_length",
         "test_jordan_length",
         "test_area",
-        "TestDensity::test_all",
     ]
 )
 def test_all():
diff --git a/tests/geometry/test_jordan_polygon.py b/tests/geometry/test_jordan_polygon.py
index 4aa56975..3aa929cf 100644
--- a/tests/geometry/test_jordan_polygon.py
+++ b/tests/geometry/test_jordan_polygon.py
@@ -7,8 +7,8 @@
 import numpy as np
 import pytest
 
+from shapepy.bool2d.density import lebesgue_density_jordan
 from shapepy.geometry.factory import FactoryJordan
-from shapepy.geometry.integral import lebesgue_density_jordan
 from shapepy.geometry.jordancurve import clean_jordan
 from shapepy.scalar.angle import degrees, radians
 
@@ -332,26 +332,6 @@ class TestIntegrateJordan:
     def test_begin(self):
         pass
 
-    @pytest.mark.order(15)
-    @pytest.mark.timeout(10)
-    @pytest.mark.dependency(depends=["TestIntegrateJordan::test_begin"])
-    def test_density_regular_polygon(self):
-        # Counter clockwise
-        for nsides in range(3, 10):
-            angles = np.linspace(0, math.tau, nsides + 1)
-            ctrlpoints = np.vstack([np.cos(angles), np.sin(angles)]).T
-            jordancurve = FactoryJordan.polygon(ctrlpoints)
-            density = lebesgue_density_jordan(jordancurve)
-            assert abs(density - 1) < 1e-9
-
-        # Clockwise
-        for nsides in range(3, 10):
-            angles = np.linspace(math.tau, 0, nsides + 1)
-            ctrlpoints = np.vstack([np.cos(angles), np.sin(angles)]).T
-            jordancurve = FactoryJordan.polygon(ctrlpoints)
-            density = lebesgue_density_jordan(jordancurve)
-            assert abs(density - 0) < 1e-9
-
     @pytest.mark.order(15)
     @pytest.mark.timeout(10)
     @pytest.mark.dependency(depends=["TestIntegrateJordan::test_begin"])
@@ -423,7 +403,6 @@ def test_lenght_regular_polygon(self):
     @pytest.mark.dependency(
         depends=[
             "TestIntegrateJordan::test_begin",
-            "TestIntegrateJordan::test_density_regular_polygon",
             "TestIntegrateJordan::test_lenght_triangle",
             "TestIntegrateJordan::test_lenght_square",
             "TestIntegrateJordan::test_lenght_regular_polygon",
diff --git a/tests/scalar/test_angle.py b/tests/scalar/test_angle.py
index 63f79086..e66c38e0 100644
--- a/tests/scalar/test_angle.py
+++ b/tests/scalar/test_angle.py
@@ -1,5 +1,6 @@
 import math
 import random
+from fractions import Fraction
 
 import pytest
 
@@ -107,6 +108,12 @@ def test_directions():
     for degval in range(226, 315):
         assert degrees(degval).direction % 4 == 3
 
+    assert degrees(0).part == 0
+    assert degrees(90).part == 0
+    assert degrees(180).part == 0
+    assert degrees(270).part == 0
+    assert degrees(45).part == 0.125
+
 
 @pytest.mark.order(2)
 @pytest.mark.timeout(1)
@@ -137,6 +144,19 @@ def test_compare():
 @pytest.mark.order(2)
 @pytest.mark.timeout(1)
 @pytest.mark.dependency(depends=["test_directions", "test_compare"])
+def test_compare_all():
+    for degval in range(0, 360):
+        angle = degrees(degval)
+        assert angle.degrees == degval
+        assert angle.turns == Fraction(degval, 360)
+        assert angle.radians == math.tau * Fraction(degval, 360)
+
+
+@pytest.mark.order(2)
+@pytest.mark.timeout(1)
+@pytest.mark.dependency(
+    depends=["test_directions", "test_compare", "test_compare_all"]
+)
 def test_evaluate_arg():
     assert arg(0, 0) == degrees(0)
     assert arg(1, 0) == degrees(0)
@@ -209,6 +229,12 @@ def test_evaluate_operations():
 
         assert 2 * anglea == anglea + anglea
 
+    for _ in range(1000):
+        a = random.randint(0, 720)
+        anglea = degrees(a)
+        angleb = degrees(a + 180)
+        assert ~angleb == anglea
+
 
 @pytest.mark.order(2)
 @pytest.mark.timeout(1)
@@ -249,6 +275,7 @@ def test_print():
         "test_build_arg",
         "test_directions",
         "test_compare",
+        "test_compare_all",
         "test_evaluate_arg",
         "test_evaluate_sincos",
         "test_evaluate_operations",