diff --git a/src/shapepy/analytic/bezier.py b/src/shapepy/analytic/bezier.py
index 2adb9fa6..13206f49 100644
--- a/src/shapepy/analytic/bezier.py
+++ b/src/shapepy/analytic/bezier.py
@@ -7,11 +7,11 @@
 from functools import lru_cache
 from typing import Iterable, Tuple, Union
 
-from ..rbool import SubSetR1, WholeR1, from_any
+from ..rbool import SubSetR1
 from ..scalar.quadrature import inner
 from ..scalar.reals import Math, Rational, Real
 from ..tools import Is, To
-from .polynomial import Polynomial
+from .polynomial import Polynomial, scale_coefs, shift_coefs
 
 
 @lru_cache(maxsize=None)
@@ -76,8 +76,13 @@ class Bezier(Polynomial):
     """
 
     def __init__(
-        self, coefs: Iterable[Real], domain: Union[None, SubSetR1] = None
+        self,
+        coefs: Iterable[Real],
+        reparam: Tuple[Real, Real] = (0, 1),
+        *,
+        domain: Union[None, SubSetR1] = None,
     ):
-        domain = WholeR1() if domain is None else from_any(domain)
-        poly_coefs = tuple(bezier2polynomial(coefs))
-        super().__init__(poly_coefs, domain)
+        coefs = tuple(bezier2polynomial(coefs))
+        knota, knotb = reparam
+        coefs = shift_coefs(scale_coefs(coefs, knotb - knota), knota)
+        super().__init__(coefs, domain=domain)
diff --git a/src/shapepy/analytic/polynomial.py b/src/shapepy/analytic/polynomial.py
index c462592d..47f9a56b 100644
--- a/src/shapepy/analytic/polynomial.py
+++ b/src/shapepy/analytic/polynomial.py
@@ -5,14 +5,36 @@
 from __future__ import annotations
 
 from numbers import Real
-from typing import Iterable, List, Union
+from typing import Iterable, Iterator, List, Union
 
 from ..rbool import IntervalR1, SubSetR1, WholeR1, from_any, move, scale
+from ..rbool.tools import is_continuous
 from ..scalar.reals import Math
 from ..tools import Is, To
 from .base import IAnalytic
 
 
+def scale_coefs(coefs: Iterable[Real], amount: Real) -> Iterator[Real]:
+    """Computes the polynomial p(t/A) from the coefficients of p(t)"""
+    if amount != 1:
+        inv = 1 / amount
+        coefs = (coef * inv**i for i, coef in enumerate(coefs))
+    yield from coefs
+
+
+def shift_coefs(coefs: Iterable[Real], amount: Real) -> Iterator[Real]:
+    """Computes the polynomial p(t-a) from the coefficients of p(t)"""
+    if amount != 0:
+        coefs = list(coefs)
+        for i, coef in enumerate(tuple(coefs)):
+            for j in range(i):
+                value = Math.binom(i, j) * (amount ** (i - j))
+                if (i + j) % 2:
+                    value *= -1
+                coefs[j] += coef * value
+    yield from coefs
+
+
 class Polynomial(IAnalytic):
     """
     Defines a polynomial with coefficients
@@ -36,7 +58,12 @@ class Polynomial(IAnalytic):
     5
     """
 
-    def __init__(self, coefs: Iterable[Real], domain: SubSetR1 = WholeR1()):
+    def __init__(
+        self, coefs: Iterable[Real], *, domain: Union[None, SubSetR1] = None
+    ):
+        domain = WholeR1() if domain is None else from_any(domain)
+        if not is_continuous(domain):
+            raise ValueError(f"Domain {domain} is not continuous")
         if not Is.iterable(coefs):
             raise TypeError("Expected an iterable of coefficients")
         coefs = tuple(coefs)
@@ -44,7 +71,7 @@ def __init__(self, coefs: Iterable[Real], domain: SubSetR1 = WholeR1()):
             raise ValueError("Cannot receive an empty tuple")
         degree = max((i for i, v in enumerate(coefs) if v * v > 0), default=0)
         self.__coefs = coefs[: degree + 1]
-        self.__domain = from_any(domain)
+        self.__domain = domain
 
     @property
     def domain(self) -> SubSetR1:
@@ -81,26 +108,28 @@ def __add__(self, other: Union[Real, Polynomial]) -> Polynomial:
         if not Is.instance(other, IAnalytic):
             coefs = list(self)
             coefs[0] += other
-            return Polynomial(coefs, self.domain)
+            return Polynomial(coefs, domain=self.domain)
         if not Is.instance(other, Polynomial):
-            return NotImplemented
+            raise NotImplementedError
         coefs = [0] * (1 + max(self.degree, other.degree))
         for i, coef in enumerate(self):
             coefs[i] += coef
         for i, coef in enumerate(other):
             coefs[i] += coef
-        return Polynomial(coefs, self.domain)
+        return Polynomial(coefs, domain=self.domain)
 
     def __mul__(self, other: Union[Real, Polynomial]) -> Polynomial:
         if not Is.instance(other, IAnalytic):
-            return Polynomial((other * coef for coef in self), self.domain)
+            return Polynomial(
+                (other * coef for coef in self), domain=self.domain
+            )
         if not Is.instance(other, Polynomial):
-            return NotImplemented
+            raise NotImplementedError
         coefs = [0 * self[0]] * (self.degree + other.degree + 1)
         for i, coefi in enumerate(self):
             for j, coefj in enumerate(other):
                 coefs[i + j] += coefi * coefj
-        return Polynomial(coefs, self.domain & other.domain)
+        return Polynomial(coefs, domain=self.domain & other.domain)
 
     def eval(self, node: Real, derivate: int = 0) -> Real:
         if node not in self.domain:
@@ -130,12 +159,12 @@ def eval(self, node: Real, derivate: int = 0) -> Real:
 
     def derivate(self, times=1):
         if self.degree < times:
-            return Polynomial([0 * self[0]], self.domain)
+            return Polynomial([0 * self[0]], domain=self.domain)
         coefs = (
             Math.factorial(n + times) // Math.factorial(n) * coef
             for n, coef in enumerate(self[times:])
         )
-        return Polynomial(coefs, self.domain)
+        return Polynomial(coefs, domain=self.domain)
 
     def integrate(self, domain):
         domain = from_any(domain)
@@ -161,21 +190,9 @@ def compose(self, function: IAnalytic) -> Polynomial:
         if function.degree != 1:
             raise ValueError("Only polynomials of degree = 1 are allowed")
         shift_amount, scale_amount = tuple(function)
-        coefs = list(self)
-        domain = self.domain
-        if scale_amount != 1:
-            inv = 1 / scale_amount
-            coefs = list(coef * inv**i for i, coef in enumerate(self))
-            domain = scale(domain, scale_amount)
-        if shift_amount != 0:
-            for i, coef in enumerate(tuple(coefs)):
-                for j in range(i):
-                    value = Math.binom(i, j) * (shift_amount ** (i - j))
-                    if (i + j) % 2:
-                        value *= -1
-                    coefs[j] += coef * value
-            domain = move(domain, shift_amount)
-        return Polynomial(coefs, domain)
+        coefs = shift_coefs(scale_coefs(self, scale_amount), shift_amount)
+        domain = move(scale(self.domain, scale_amount), shift_amount)
+        return Polynomial(coefs, domain=domain)
 
     def __repr__(self):
         return str(self.domain) + ": " + self.__str__()
diff --git a/src/shapepy/analytic/tools.py b/src/shapepy/analytic/tools.py
index 4a658c61..b723ea91 100644
--- a/src/shapepy/analytic/tools.py
+++ b/src/shapepy/analytic/tools.py
@@ -62,6 +62,11 @@ def find_minimum(
     raise NotExpectedError(f"Invalid analytic: {type(analytic)}")
 
 
+def is_constant(analytic: IAnalytic) -> bool:
+    """Tells if the given analytic function is constant"""
+    return Is.instance(analytic, Polynomial) and analytic.degree == 0
+
+
 class PolynomialFunctions:
     """Static class that stores static functions used for the generics
     functions above. This class specifics for Polynomial"""
diff --git a/src/shapepy/bool2d/boolean.py b/src/shapepy/bool2d/boolean.py
index 5d034e6e..adbc853a 100644
--- a/src/shapepy/bool2d/boolean.py
+++ b/src/shapepy/bool2d/boolean.py
@@ -6,7 +6,6 @@
 from __future__ import annotations
 
 from copy import copy
-from fractions import Fraction
 from typing import Dict, Iterable, Tuple, Union
 
 from shapepy.geometry.jordancurve import JordanCurve
@@ -370,12 +369,14 @@ def split_on_intersection(
                 jordansj = all_group_jordans[j]
                 for jordana in jordansi:
                     for jordanb in jordansj:
-                        intersection |= jordana.piecewise & jordanb.piecewise
+                        intersection |= (
+                            jordana.parametrize() & jordanb.parametrize()
+                        )
         intersection.evaluate()
         for jordans in all_group_jordans:
             for jordan in jordans:
-                split_knots = intersection.all_knots[id(jordan.piecewise)]
-                jordan.piecewise.split(split_knots)
+                split_knots = intersection.all_knots[id(jordan.parametrize())]
+                jordan.parametrize().split(split_knots)
 
     @staticmethod
     def pursue_path(
@@ -396,14 +397,14 @@ def pursue_path(
         We suppose there's no triple intersection
         """
         matrix = []
-        all_segments = [tuple(jordan.piecewise) for jordan in jordans]
+        all_segments = [tuple(jordan.parametrize()) for jordan in jordans]
         while True:
             index_segment %= len(all_segments[index_jordan])
             segment = all_segments[index_jordan][index_segment]
             if (index_jordan, index_segment) in matrix:
                 break
             matrix.append((index_jordan, index_segment))
-            last_point = segment(1)
+            last_point = segment(segment.knots[-1])
             possibles = []
             for i, jordan in enumerate(jordans):
                 if i == index_jordan:
@@ -415,7 +416,7 @@ def pursue_path(
                 continue
             index_jordan = possibles[0]
             for j, segj in enumerate(all_segments[index_jordan]):
-                if segj(0) == last_point:
+                if segj(segj.knots[0]) == last_point:
                     index_segment = j
                     break
         return CyclicContainer(matrix)
@@ -434,7 +435,7 @@ def indexs_to_jordan(
         """
         beziers = []
         for index_jordan, index_segment in matrix_indexs:
-            new_bezier = jordans[index_jordan].piecewise[index_segment]
+            new_bezier = jordans[index_jordan].parametrize()[index_segment]
             new_bezier = copy(new_bezier)
             beziers.append(USegment(new_bezier))
         new_jordan = JordanCurve(beziers)
@@ -448,7 +449,7 @@ def follow_path(
         Returns a list of jordan curves which is the result
         of the intersection between 'jordansa' and 'jordansb'
         """
-        assert all(map(Is.jordan, jordans))
+        assert all(Is.instance(j, JordanCurve) for j in jordans)
         bez_indexs = []
         for ind_jord, ind_seg in start_indexs:
             indices_matrix = FollowPath.pursue_path(ind_jord, ind_seg, jordans)
@@ -480,7 +481,7 @@ def midpoints_one_shape(
         """
         for i, jordan in enumerate(shapea.jordans):
             for j, segment in enumerate(jordan.parametrize()):
-                mid_point = segment(Fraction(1, 2))
+                mid_point = segment((segment.knots[0] + segment.knots[-1]) / 2)
                 density = shapeb.density(mid_point)
                 mid_point_in = (float(density) > 0 and closed) or density == 1
                 if not inside ^ mid_point_in:
diff --git a/src/shapepy/bool2d/density.py b/src/shapepy/bool2d/density.py
index 54b19dc9..e9651588 100644
--- a/src/shapepy/bool2d/density.py
+++ b/src/shapepy/bool2d/density.py
@@ -32,13 +32,13 @@ def half_density_jordan(
         deltax = segment.xfunc - point.xcoord
         deltay = segment.yfunc - point.ycoord
         radius_square = deltax * deltax + deltay * deltay
-        minimal = find_minimum(radius_square, [0, 1])
+        minimal = find_minimum(radius_square, segment.domain)
         if minimal < 1e-6:
-            place = where_minimum(radius_square, [0, 1])
+            place = where_minimum(radius_square, segment.domain)
             if not Is.instance(place, SingleR1):
                 raise NotExpectedError(f"Not single value: {place}")
             parameter = To.finite(place.internal)
-            angle = segment(parameter, 1).angle
+            angle = segment.eval(parameter, 1).angle
             return line(angle)
     raise NotExpectedError("Not found minimum < 1e-6")
 
@@ -61,10 +61,10 @@ def lebesgue_density_jordan(
 
     segments = tuple(jordan.parametrize())
     for i, segmenti in enumerate(segments):
-        if point == segmenti(0):
+        if point == segmenti(segmenti.knots[0]):
             segmentj = segments[(i - 1) % len(segments)]
-            anglei = segmenti(0, 1).angle
-            anglej = segmentj(1, 1).angle
+            anglei = segmenti.eval(segmenti.knots[0], 1).angle
+            anglej = segmentj.eval(segmentj.knots[-1], 1).angle
             return sector(anglei, ~anglej)
 
     turns = IntegrateJordan.turns(jordan, point)
diff --git a/src/shapepy/geometry/base.py b/src/shapepy/geometry/base.py
index fc1dc92a..9420d80b 100644
--- a/src/shapepy/geometry/base.py
+++ b/src/shapepy/geometry/base.py
@@ -5,9 +5,11 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
-from typing import Iterable, Tuple
+from typing import Iterable, Tuple, Union
 
+from ..rbool import IntervalR1, WholeR1
 from ..scalar.reals import Real
+from ..tools import vectorize
 from .box import Box
 from .point import Point2D
 
@@ -72,19 +74,31 @@ class IParametrizedCurve(IGeometricCurve):
     Class interface for parametrized curves
     """
 
+    @property
+    @abstractmethod
+    def domain(self) -> Union[IntervalR1, WholeR1]:
+        """
+        The domain where the curve is defined.
+        """
+        raise NotImplementedError
+
     @property
     @abstractmethod
     def knots(self) -> Tuple[Real, ...]:
         """
-        The length of the curve
-        If the curve is not bounded, returns infinity
+        The subdivisions on the domain
         """
         raise NotImplementedError
 
     @abstractmethod
-    def __call__(self, node: Real, derivate: int = 0) -> Point2D:
+    def eval(self, node: Real, derivate: int = 0) -> Point2D:
+        """Evaluates the curve at given node"""
         raise NotImplementedError
 
+    @vectorize(1, 0)
+    def __call__(self, node: Real) -> Point2D:
+        return self.eval(node, 0)
+
     def __and__(self, other: IParametrizedCurve):
         return Future.intersect(self, other)
 
diff --git a/src/shapepy/geometry/concatenate.py b/src/shapepy/geometry/concatenate.py
index e4310bde..c3ead8d1 100644
--- a/src/shapepy/geometry/concatenate.py
+++ b/src/shapepy/geometry/concatenate.py
@@ -4,11 +4,9 @@
 
 from typing import Iterable, Union
 
-from ..analytic import Bezier
 from ..tools import Is, NotExpectedError
 from .base import IGeometricCurve
 from .piecewise import PiecewiseCurve
-from .point import cross
 from .segment import Segment
 from .unparam import UPiecewiseCurve, USegment
 
@@ -23,7 +21,7 @@ def concatenate(curves: Iterable[IGeometricCurve]) -> IGeometricCurve:
     if not all(Is.instance(curve, IGeometricCurve) for curve in curves):
         raise ValueError
     if all(Is.instance(curve, Segment) for curve in curves):
-        return concatenate_segments(curves)
+        return simplify_piecewise(PiecewiseCurve(curves))
     if all(Is.instance(curve, USegment) for curve in curves):
         return concatenate_usegments(curves)
     raise NotExpectedError(str(tuple(str(type(c)) for c in curves)))
@@ -35,32 +33,28 @@ def concatenate_usegments(
     """
     Concatenates all the unparametrized segments
     """
-    usegments = tuple(usegments)
-    assert all(Is.instance(useg, USegment) for useg in usegments)
-    union = concatenate_segments(useg.parametrize() for useg in usegments)
+    union = simplify_piecewise(UPiecewiseCurve(usegments).parametrize())
     return (
         USegment(union)
         if Is.instance(union, Segment)
-        else UPiecewiseCurve(map(USegment, union))
+        else UPiecewiseCurve(union)
     )
 
 
-def concatenate_segments(
-    segments: Iterable[Segment],
+def simplify_piecewise(
+    piecewise: PiecewiseCurve,
 ) -> Union[Segment, PiecewiseCurve]:
     """
     Concatenates all the segments
     """
-    segments = tuple(segments)
-    if len(segments) == 0:
-        raise ValueError(f"Number sizes: {len(segments)}")
     filtsegments = []
-    segments = iter(segments)
+    segments = iter(piecewise)
     segmenti = next(segments)
     for segmentj in segments:
-        try:
-            segmenti = bezier_and_bezier(segmenti, segmentj)
-        except ValueError:
+        if can_concatenate(segmenti, segmentj):
+            domain = segmenti.domain | segmentj.domain
+            segmenti = Segment(segmenti.xfunc, segmenti.yfunc, domain=domain)
+        else:
             filtsegments.append(segmenti)
             segmenti = segmentj
     filtsegments.append(segmenti)
@@ -71,25 +65,6 @@ def concatenate_segments(
     )
 
 
-def bezier_and_bezier(curvea: Segment, curveb: Segment) -> Segment:
-    """Computes the union of two bezier curves"""
-    if not Is.instance(curvea, Segment):
-        raise TypeError(f"Invalid type: {type(curvea)}")
-    if not Is.instance(curveb, Segment):
-        raise TypeError(f"Invalid type: {type(curveb)}")
-    if abs(cross(curvea(1, 1), curveb(0, 1))) > 1e-6:
-        raise ValueError
-    if curvea.xfunc.degree != curveb.xfunc.degree:
-        raise ValueError
-    if curvea.yfunc.degree != curveb.yfunc.degree:
-        raise ValueError
-    if curvea.xfunc.degree > 1 or curvea.yfunc.degree > 1:
-        raise ValueError
-
-    if curvea(1) != curveb(0):
-        raise ValueError
-    startpoint = curvea(0)
-    endpoint = curveb(1)
-    nxfunc = Bezier((startpoint[0], endpoint[0]))
-    nyfunc = Bezier((startpoint[1], endpoint[1]))
-    return Segment(nxfunc, nyfunc)
+def can_concatenate(curvea: Segment, curveb: Segment) -> bool:
+    """Tells if it's possible to concatenate two segments into a single one"""
+    return curvea.xfunc == curveb.xfunc and curvea.yfunc == curveb.yfunc
diff --git a/src/shapepy/geometry/factory.py b/src/shapepy/geometry/factory.py
index d3e40149..4d30c57f 100644
--- a/src/shapepy/geometry/factory.py
+++ b/src/shapepy/geometry/factory.py
@@ -11,7 +11,9 @@
 
 from ..analytic import Bezier
 from ..loggers import debug
-from ..tools import To
+from ..rbool import IntervalR1
+from ..scalar.reals import Real
+from ..tools import To, pairs
 from .jordancurve import JordanCurve
 from .point import Point2D, cartesian, rotate
 from .segment import Segment
@@ -26,7 +28,9 @@ class FactorySegment:
 
     @staticmethod
     @debug("shapepy.geometry.factory")
-    def bezier(ctrlpoints: Iterable[Point2D]) -> Segment:
+    def bezier(
+        ctrlpoints: Iterable[Point2D], limits: Tuple[Real, Real] = (0, 1)
+    ) -> Segment:
         """Initialize a bezier segment from a list of control points
 
         :param ctrlpoints: The list of control points
@@ -34,10 +38,11 @@ def bezier(ctrlpoints: Iterable[Point2D]) -> Segment:
         :return: The created segment
         :rtype: Segment
         """
+        domain = IntervalR1(limits[0], limits[1], True, True)
         ctrlpoints = tuple(map(To.point, ctrlpoints))
-        xfunc = Bezier((pt[0] for pt in ctrlpoints))
-        yfunc = Bezier((pt[1] for pt in ctrlpoints))
-        return Segment(xfunc, yfunc)
+        xfunc = Bezier((pt[0] for pt in ctrlpoints), limits)
+        yfunc = Bezier((pt[1] for pt in ctrlpoints), limits)
+        return Segment(xfunc, yfunc, domain=domain)
 
 
 class FactoryJordan:
@@ -65,14 +70,11 @@ def polygon(vertices: Tuple[Point2D, ...]) -> JordanCurve:
         ((0, 0), (4, 0), (0, 3))
 
         """
-        vertices = list(map(To.point, vertices))
-        nverts = len(vertices)
-        vertices.append(vertices[0])
-        beziers = [None] * nverts
-        for i in range(nverts):
-            ctrlpoints = vertices[i : i + 2]
-            new_bezier = FactorySegment.bezier(ctrlpoints)
-            beziers[i] = USegment(new_bezier)
+        vertices = tuple(map(To.point, vertices))
+        beziers = []
+        for i, points in enumerate(pairs(vertices, cyclic=True)):
+            new_bezier = FactorySegment.bezier(points, (i, i + 1))
+            beziers.append(USegment(new_bezier))
         return JordanCurve(beziers)
 
     @staticmethod
@@ -125,6 +127,8 @@ def circle(ndivangle: int):
             middle_point = rotate(middle_point, angle)
         end_point = all_ctrlpoints[0][0]
         all_ctrlpoints.append([start_point, middle_point, end_point])
-        return JordanCurve(
-            map(USegment, map(FactorySegment.bezier, all_ctrlpoints))
+        segments = (
+            FactorySegment.bezier(pts, [i, i + 1])
+            for i, pts in enumerate(all_ctrlpoints)
         )
+        return JordanCurve(segments)
diff --git a/src/shapepy/geometry/integral.py b/src/shapepy/geometry/integral.py
index 34817b58..f323cc5a 100644
--- a/src/shapepy/geometry/integral.py
+++ b/src/shapepy/geometry/integral.py
@@ -8,6 +8,7 @@
 
 from ..analytic.base import IAnalytic
 from ..analytic.tools import find_minimum
+from ..loggers import debug, get_logger
 from ..scalar.quadrature import AdaptativeIntegrator, IntegratorFactory
 from ..scalar.reals import Math
 from ..tools import Is, To
@@ -26,6 +27,7 @@ class IntegrateSegment:
     adaptative = AdaptativeIntegrator(direct, 1e-6)
 
     @staticmethod
+    @debug("shapepy.geometry.integral")
     def polynomial(curve: Segment, expx: int, expy: int):
         """
         Computes the integral
@@ -39,8 +41,10 @@ def polynomial(curve: Segment, expx: int, expy: int):
         pcrossdp = xfunc * yfunc.derivate()
         pcrossdp -= yfunc * xfunc.derivate()
         function = (xfunc**expx) * (yfunc**expy) * pcrossdp
+        logger = get_logger("shapepy.geometry.integral")
+        logger.debug(f"Integrating {function} over {curve.domain}")
         assert Is.instance(function, IAnalytic)
-        return function.integrate([0, 1]) / (expx + expy + 2)
+        return function.integrate(curve.domain) / (expx + expy + 2)
 
     @staticmethod
     def turns(curve: Segment, point: Point2D) -> float:
@@ -58,14 +62,14 @@ def turns(curve: Segment, point: Point2D) -> float:
         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:
+        if find_minimum(radius_square, curve.domain) < 1e-6:
             return To.rational(1, 2)
         crossf = deltax * deltay.derivate()
         crossf -= deltay * deltax.derivate()
         function = partial(
             lambda t, cf, rs: cf(t) / rs(t), cf=crossf, rs=radius_square
         )
-        radians = IntegrateSegment.adaptative.integrate(function, [0, 1])
+        radians = IntegrateSegment.adaptative.integrate(function, curve.domain)
         return radians / Math.tau
 
 
@@ -76,13 +80,14 @@ class IntegrateJordan:
     """
 
     @staticmethod
+    @debug("shapepy.geometry.integral")
     def polynomial(jordan: JordanCurve, expx: int, expy: int):
         """
         Computes the integral
 
         I = int x^expx * y^expy * ds
         """
-        assert Is.jordan(jordan)
+        assert Is.instance(jordan, JordanCurve)
         return sum(
             IntegrateSegment.polynomial(usegment.parametrize(), expx, expy)
             for usegment in jordan
diff --git a/src/shapepy/geometry/intersection.py b/src/shapepy/geometry/intersection.py
index 309d2fe2..df1665ab 100644
--- a/src/shapepy/geometry/intersection.py
+++ b/src/shapepy/geometry/intersection.py
@@ -14,6 +14,7 @@
 from fractions import Fraction
 from typing import Dict, Iterable, Set, Tuple, Union
 
+from ..loggers import debug, get_logger
 from ..rbool import (
     EmptyR1,
     SubSetR1,
@@ -24,7 +25,7 @@
 )
 from ..scalar.nodes_sample import NodeSampleFactory
 from ..scalar.reals import Real
-from ..tools import Is, NotExpectedError
+from ..tools import Is
 from .base import IGeometricCurve, IParametrizedCurve
 from .piecewise import PiecewiseCurve
 from .point import cross, inner
@@ -206,11 +207,19 @@ def param_and_param(
     """Computes the intersection between two parametrized curves"""
     assert Is.instance(curvea, IParametrizedCurve)
     assert Is.instance(curveb, IParametrizedCurve)
-    if Is.segment(curvea) and Is.segment(curveb):
-        return segment_and_segment(curvea, curveb)
-    if Is.piecewise(curvea) and Is.piecewise(curveb):
-        return intersect_piecewises(curvea, curveb)
-    raise NotExpectedError
+    if curvea.box() & curveb.box() is None:
+        return EmptyR1(), EmptyR1()
+    if Is.instance(curvea, PiecewiseCurve):
+        subseta, subsetb = EmptyR1(), EmptyR1()
+        for segmenta in curvea:
+            subb, suba = param_and_param(segmenta, curveb)
+            subseta |= suba
+            subsetb |= subb
+        return subseta, subsetb
+    if Is.instance(curveb, PiecewiseCurve):
+        # pylint: disable=arguments-out-of-order
+        return param_and_param(curveb, curvea)[::-1]
+    return segment_and_segment(curvea, curveb)
 
 
 def segment_is_linear(segment: Segment) -> bool:
@@ -224,25 +233,29 @@ def segment_and_segment(
     """Computes the intersection between two segment curves"""
     assert Is.instance(curvea, Segment)
     assert Is.instance(curveb, Segment)
-    if curvea.box() & curveb.box() is None:
-        return EmptyR1(), EmptyR1()
     if curvea == curveb:
-        return create_interval(0, 1), create_interval(0, 1)
+        return curvea.domain, curveb.domain
     if segment_is_linear(curvea) and segment_is_linear(curveb):
         return IntersectionSegments.lines(curvea, curveb)
     nptsa = max(curvea.xfunc.degree, curvea.yfunc.degree) + 4
     nptsb = max(curveb.xfunc.degree, curveb.yfunc.degree) + 4
-    usample = list(NodeSampleFactory.closed_linspace(nptsa))
-    vsample = list(NodeSampleFactory.closed_linspace(nptsb))
+    usample = [
+        curvea.knots[0] + u * (curvea.knots[-1] - curvea.knots[0])
+        for u in NodeSampleFactory.closed_linspace(nptsa)
+    ]
+    vsample = [
+        curveb.knots[0] + u * (curveb.knots[-1] - curveb.knots[0])
+        for u in NodeSampleFactory.closed_linspace(nptsb)
+    ]
     pairs = []
     for ui in usample:
         pairs += [(ui, vj) for vj in vsample]
     for _ in range(3):
         pairs = IntersectionSegments.bezier_and_bezier(curvea, curveb, pairs)
-        pairs.insert(0, (0, 0))
-        pairs.insert(0, (0, 1))
-        pairs.insert(0, (1, 0))
-        pairs.insert(0, (1, 1))
+        pairs.insert(0, (curvea.knots[0], curveb.knots[0]))
+        pairs.insert(0, (curvea.knots[0], curveb.knots[-1]))
+        pairs.insert(0, (curvea.knots[-1], curveb.knots[0]))
+        pairs.insert(0, (curvea.knots[-1], curveb.knots[-1]))
         # Filter values by distance of points
         tol_norm = 1e-6
         pairs = IntersectionSegments.filter_distance(
@@ -267,11 +280,15 @@ class IntersectionSegments:
 
     # pylint: disable=invalid-name, too-many-return-statements, too-many-locals
     @staticmethod
+    @debug("shapepy.geometry.intersection")
     def lines(curvea: Segment, curveb: Segment) -> Tuple[SubSetR1, SubSetR1]:
         """Finds the intersection of two line segments"""
         empty = EmptyR1()
-        A0, A1 = curvea(0), curvea(1)
-        B0, B1 = curveb(0), curveb(1)
+        logger = get_logger("shapepy.geometry.intersection")
+        A0 = curvea(curvea.knots[0])
+        A1 = curvea(curvea.knots[-1])
+        B0 = curveb(curveb.knots[0])
+        B1 = curveb(curveb.knots[-1])
         dA = A1 - A0
         dB = B1 - B0
         B0mA0 = B0 - A0
@@ -279,13 +296,19 @@ def lines(curvea: Segment, curveb: Segment) -> Tuple[SubSetR1, SubSetR1]:
         if dAxdB != 0:  # Lines are not parallel
             t0 = cross(B0mA0, dB) / dAxdB
             if t0 < 0 or 1 < t0:
+                logger.debug("1) Empty, Empty")
                 return empty, empty
             u0 = cross(B0mA0, dA) / dAxdB
             if u0 < 0 or 1 < u0:
+                logger.debug("2) Empty, Empty")
                 return empty, empty
+            t0 = curvea.knots[0] + t0 * (curvea.knots[-1] - curvea.knots[0])
+            u0 = curveb.knots[0] + u0 * (curveb.knots[-1] - curveb.knots[0])
+            logger.debug("3) Single, Single")
             return create_single(t0), create_single(u0)
         # Lines are parallel
         if cross(dA, B0mA0) != 0:
+            logger.debug("4) Empty, Empty")
             return empty, empty  # Parallel, but not colinear
         # Compute the projections
         dAodA = inner(dA, dA)
@@ -305,8 +328,14 @@ def lines(curvea: Segment, curveb: Segment) -> Tuple[SubSetR1, SubSetR1]:
         t1 = min(max(0, t1), 1)
         u0 = min(max(0, u0), 1)
         u1 = min(max(0, u1), 1)
+        t0 = curvea.knots[0] + t0 * (curvea.knots[-1] - curvea.knots[0])
+        u0 = curveb.knots[0] + u0 * (curveb.knots[-1] - curveb.knots[0])
+        t1 = curvea.knots[0] + t1 * (curvea.knots[-1] - curvea.knots[0])
+        u1 = curveb.knots[0] + u1 * (curveb.knots[-1] - curveb.knots[0])
         if t0 == t1 or u0 == u1:
+            logger.debug("6) Single, Single")
             return create_single(t0), create_single(u1)
+        logger.debug("7) Interval, Interval")
         return create_interval(t0, t1), create_interval(u0, u1)
 
     # pylint: disable=too-many-locals
@@ -345,8 +374,9 @@ def bezier_and_bezier(
                     continue
                 newu = u - (mat11 * vect0 - mat01 * vect1) / deter
                 newv = v - (mat00 * vect1 - mat01 * vect0) / deter
-                pair = (min(1, max(0, newu)), min(1, max(0, newv)))
-                new_pairs.add(pair)
+                newu = min(curvea.knots[-1], max(curvea.knots[0], newu))
+                newv = min(curveb.knots[-1], max(curveb.knots[0], newv))
+                new_pairs.add((newu, newv))
             pairs = list(new_pairs)
             for i, (ui, vi) in enumerate(pairs):
                 if Is.instance(ui, Fraction):
@@ -399,81 +429,3 @@ def filter_parameters(
                     j += 1
             index += 1
         return tuple(pairs)
-
-
-def intersect_piecewises(
-    curvea: PiecewiseCurve,
-    curveb: PiecewiseCurve,
-) -> Tuple[SubSetR1, SubSetR1]:
-    r"""Computes the intersection between two jordan curves
-
-    Finds the values of (:math:`a^{\star}`, :math:`b^{\star}`,
-    :math:`u^{\star}`, :math:`v^{\star}`) such
-
-    .. math::
-        S_{a^{\star}}(u^{\star}) == O_{b^{\star}}(v^{\star})
-
-    It computes the intersection between each pair of segments
-    from ``self`` and ``other`` and returns the matrix of coefficients
-
-    .. math::
-
-        \begin{bmatrix}
-        a_0 & b_0 & u_0 & v_0 \\
-        a_1 & b_1 & u_1 & v_1 \\
-        \vdots & \vdots & \vdots & \vdots \\
-        a_{n} & b_{n} & u_{n} & v_{n}
-        \end{bmatrix}
-
-    If two bezier curves are equal, then ``u_i = v_i = None``
-
-    * ``0 <= a_i < len(self.segments)``
-    * ``0 <= b_i < len(other.segments)``
-    * ``0 <= u_i <= 1`` or ``None``
-    * ``0 <= v_i <= 1`` or ``None``
-
-    Parameters
-    ----------
-    other : JordanCurve
-        The jordan curve which intersects ``self``
-    equal_beziers : bool, default = True
-        Flag to return (or not) when two segments are equal
-
-        If the flag ``equal_beziers`` are inactive,
-        then will remove when ``(ui, vi) == (None, None)``.
-
-    end_points : bool, default = True
-        Flag to return (or not) when jordans intersect at end points
-
-        If the flag ``end_points`` are inactive,
-        then will remove when ``(ui, vi)`` are
-        ``(0, 0)``, ``(0, 1)``, ``(1, 0)`` or ``(1, 1)``
-
-    :return: The matrix of coefficients ``[(ai, bi, ui, vi)]``
-                or an empty tuple in case of non-intersection
-    :rtype: tuple[(int, int, float, float)]
-
-
-    Example use
-    -----------
-    >>> from shapepy import JordanCurve
-    >>> vertices_a = [(0, 0), (2, 0), (2, 2), (0, 2)]
-    >>> jordan_a = FactoryJordan.polygon(vertices_a)
-    >>> vertices_b = [(1, 1), (3, 1), (3, 3), (1, 3)]
-    >>> jordan_b = FactoryJordan.polygon(vertices_b)
-    >>> jordan_a.intersection(jordan_b)
-    ((1, 0, 1/2, 1/2), (2, 3, 1/2, 1/2))
-
-    """
-    assert Is.piecewise(curvea)
-    assert Is.piecewise(curveb)
-
-    subseta, subsetb = EmptyR1(), EmptyR1()
-    for ai, sbezier in enumerate(curvea):
-        for bj, obezier in enumerate(curveb):
-            suba, subb = segment_and_segment(sbezier, obezier)
-            suba = suba.scale(curvea.knots[ai + 1] - curvea.knots[ai])
-            subseta |= suba.move(curvea.knots[ai])
-            subb = subb.scale(curveb.knots[bj + 1] - curveb.knots[bj])
-            subsetb |= subb.move(curveb.knots[bj])
-    return subseta, subsetb
diff --git a/src/shapepy/geometry/jordancurve.py b/src/shapepy/geometry/jordancurve.py
index abcf1e02..602d2453 100644
--- a/src/shapepy/geometry/jordancurve.py
+++ b/src/shapepy/geometry/jordancurve.py
@@ -6,76 +6,31 @@
 from __future__ import annotations
 
 from collections import deque
-from copy import copy
-from typing import Iterable, Iterator
+from typing import Iterable, Iterator, List, Union
 
 from ..analytic import IAnalytic
-from ..loggers import debug
+from ..loggers import debug, get_logger
 from ..scalar.reals import Real
 from ..tools import CyclicContainer, Is, pairs, reverse
-from .base import IGeometricCurve
-from .box import Box
-from .piecewise import PiecewiseCurve
-from .point import Point2D
+from .base import Future
+from .point import Point2D, cross
 from .segment import Segment
-from .unparam import UPiecewiseCurve, USegment, clean_usegment, self_intersect
+from .unparam import UPiecewiseCurve, USegment, self_intersect
 
 
-class JordanCurve(IGeometricCurve):
+class JordanCurve(UPiecewiseCurve):
     """
     Jordan Curve is an arbitrary closed curve which doesn't intersect itself.
     It stores a list of 'segments', each segment is a bezier curve
     """
 
-    def __init__(self, usegments: Iterable[USegment]):
-        usegments = tuple(usegments)
-        if not all(Is.instance(u, (USegment, Segment)) for u in usegments):
-            raise ValueError(f"Invalid tipos: {tuple(map(type, usegments))}")
-        usegments = tuple(
-            USegment(s) if Is.instance(s, Segment) else s for s in usegments
-        )
-        if any(map(self_intersect, usegments)):
-            raise ValueError("Segment must not self intersect")
-        for usegi, usegj in pairs(usegments, cyclic=True):
-            if usegi.end_point != usegj.start_point:
-                raise ValueError("The segments are not continuous")
-        self.__usegments = CyclicContainer(usegments)
-        self.__piecewise = UPiecewiseCurve(self.__usegments).parametrize()
+    def __init__(self, usegments: Iterable[Union[Segment, USegment]]):
+        super().__init__(clean_jordan(usegments))
+        for usegi in self:
+            if self_intersect(usegi):
+                raise ValueError(f"Segment self-intersect! {usegi}")
         self.__area = None
 
-    def __copy__(self) -> JordanCurve:
-        return self.__deepcopy__(None)
-
-    def __deepcopy__(self, memo) -> JordanCurve:
-        """Returns a deep copy of the jordan curve"""
-        return self.__class__(map(copy, self))
-
-    def box(self) -> Box:
-        """The box which encloses the jordan curve
-
-        :return: The box which encloses the jordan curve
-        :rtype: Box
-
-        Example use
-        -----------
-
-        >>> from shapepy import JordanCurve
-        >>> vertices = [(0, 0), (4, 0), (0, 3)]
-        >>> jordan = FactoryJordan.polygon(vertices)
-        >>> jordan.box()
-        Box with vertices (0, 0) and (4, 3)
-
-        """
-        box = None
-        for usegment in self:
-            box |= usegment.box()
-        return box
-
-    @property
-    def length(self) -> Real:
-        """The length of the curve"""
-        return sum(useg.length for useg in self)
-
     @property
     def area(self) -> Real:
         """The internal area"""
@@ -83,43 +38,6 @@ def area(self) -> Real:
             self.__area = compute_area(self)
         return self.__area
 
-    def parametrize(self) -> PiecewiseCurve:
-        """
-        Gives the piecewise curve
-        """
-        return self.__piecewise
-
-    @property
-    def piecewise(self) -> PiecewiseCurve:
-        """
-        Gives the internal piecewise curve
-        """
-        return self.__piecewise
-
-    def __iter__(self) -> Iterator[USegment]:
-        """Unparametrized Segments
-
-        When setting, it checks if the points are the same between
-        the junction of two segments to ensure a closed curve
-
-        :getter: Returns the tuple of connected planar beziers, not copy
-        :setter: Sets the segments of the jordan curve
-        :type: tuple[Segment]
-
-        Example use
-        -----------
-
-        >>> from shapepy import JordanCurve
-        >>> vertices = [(0, 0), (4, 0), (0, 3)]
-        >>> jordan = FactoryJordan.polygon(vertices)
-        >>> print(tuple(jordan))
-        (Segment (deg 1), Segment (deg 1), Segment (deg 1))
-        >>> print(tuple(jordan)[0])
-        Planar curve of degree 1 and control points ((0, 0), (4, 0))
-
-        """
-        yield from self.__usegments
-
     def vertices(self) -> Iterator[Point2D]:
         """Vertices
 
@@ -139,56 +57,47 @@ def vertices(self) -> Iterator[Point2D]:
         ((0, 0), (4, 0), (0, 3))
 
         """
-        yield from (useg.start_point for useg in self)
+        for useg in self:
+            seg = useg.parametrize()
+            yield seg.eval(seg.knots[0])
 
     def __str__(self) -> str:
-        nsegs = len(self.__usegments)
-        msg = f"Jordan Curve with {nsegs} segments and vertices\n"
+        msg = f"Jordan Curve with {len(self)} segments and vertices\n"
         msg += str(self.vertices())
         return msg
 
     def __repr__(self) -> str:
-        nsegs = len(self.__usegments)
         box = self.box()
-        return f"JC[{nsegs}:{box.lowpt},{box.toppt}]"
+        return f"JC[{len(self)}:{box.lowpt},{box.toppt}]"
 
     def __eq__(self, other: JordanCurve) -> bool:
+        logger = get_logger("shapepy.geometry.jordancurve")
+        logger.debug(f"     type: {type(other)}")
+        logger.debug(f"     box: {self.box() == other.box()}")
+        logger.debug(f"     len: {self.length == other.length}")
+        logger.debug(f"    area: {self.area == other.area}")
+        logger.debug(
+            f"    all1: {all(point in self for point in other.vertices())}"
+        )
+        logger.debug(
+            f"    all2: {all(point in other for point in self.vertices())}"
+        )
+        logger.debug(
+            f"    cycl: {CyclicContainer(self) == CyclicContainer(other)}"
+        )
         return (
             Is.instance(other, JordanCurve)
             and self.box() == other.box()
             and self.length == other.length
+            and self.area == other.area
             and all(point in self for point in other.vertices())
             and all(point in other for point in self.vertices())
-            and CyclicContainer(self.clean()) == CyclicContainer(other.clean())
+            and CyclicContainer(self) == CyclicContainer(other)
         )
 
-    @debug("shapepy.geometry.jordancurve")
-    def clean(self) -> JordanCurve:
-        """Cleans the jordan curve"""
-        usegments = list(map(clean_usegment, self))
-        index = 0
-        while index + 1 < len(usegments):
-            union = usegments[index] | usegments[index + 1]
-            if not Is.instance(union, USegment):
-                index += 1
-            else:
-                usegments.pop(index)
-                usegments[index] = union
-        while len(usegments) > 1:
-            union = usegments[-1] | usegments[0]
-            if not Is.instance(union, USegment):
-                break
-            usegments.pop(0)
-            usegments[-1] = union
-        return JordanCurve(usegments)
-
     def __invert__(self) -> JordanCurve:
         return JordanCurve(reverse(~useg for useg in self))
 
-    def __contains__(self, point: Point2D) -> bool:
-        """Tells if the point is on the boundary"""
-        return any(point in useg for useg in self)
-
 
 @debug("shapepy.geometry.jordancurve")
 def compute_area(jordan: JordanCurve) -> Real:
@@ -205,18 +114,20 @@ def compute_area(jordan: JordanCurve) -> Real:
         poly = xfunc * yfunc.derivate()
         poly -= yfunc * xfunc.derivate()
         assert Is.instance(poly, IAnalytic)
-        total += poly.integrate([0, 1])
+        total += poly.integrate(segment.domain)
     return total / 2
 
 
-def clean_jordan(jordan: JordanCurve) -> JordanCurve:
+@debug("shapepy.geometry.jordan")
+def clean_jordan(
+    usegments: Iterable[Union[Segment, USegment]],
+) -> Iterator[Union[Segment, USegment]]:
     """Cleans the jordan curve
 
-    Removes the uncessary nodes from jordan curve,
-    for example, after calling ``split`` function
+    Removes the uncessary nodes from jordan curve
 
-    :return: The same curve
-    :rtype: JordanCurve
+    :return: The set of segments
+    :rtype: A set of segments
 
     Example use
     -----------
@@ -229,31 +140,28 @@ def clean_jordan(jordan: JordanCurve) -> JordanCurve:
     ((0, 0), (4, 0), (0, 3))
 
     """
-    usegments = deque(map(clean_usegment, jordan))
-    for _ in range(len(usegments) + 1):
-        union = usegments[0] | usegments[1]
-        if Is.instance(union, USegment):
-            usegments.popleft()
-            usegments.popleft()
-            usegments.appendleft(union)
+    logger = get_logger("shapepy.geometry.jordan")
+    logger.debug("Segments = ")
+    usegments = deque(usegments)
+    endvectors = []
+    for i, usegment in enumerate(usegments):
+        if not Is.instance(usegment, (Segment, USegment)):
+            raise ValueError
+        segment = usegment.parametrize()
+        logger.debug(f"    {i}: {segment}")
+        vectora = segment.eval(segment.knots[0], 1)
+        vectorb = segment.eval(segment.knots[-1], 1)
+        endvectors.append((vectora, vectorb))
+    crosses: List[Real] = []
+    for (_, vectb), (vecta, _) in pairs(endvectors, cyclic=True):
+        crosses.append(cross(vectb, vecta))
+    logger.debug(f"End vectors = {endvectors}")
+    logger.debug(f"Crosses = {crosses}")
+    if not any(c == 0 for c in crosses):
+        return usegments
+    index = len(crosses) - 1
+    while index > 0 and crosses[index] == 0:
         usegments.rotate()
-    return JordanCurve(usegments)
-
-
-def is_jordan(obj: object) -> bool:
-    """
-    Checks if the parameter is a Jordan Curve
-
-    Parameters
-    ----------
-    obj : The object to be tested
-
-    Returns
-    -------
-    bool
-        True if the obj is a Jordan Curve, False otherwise
-    """
-    return Is.instance(obj, JordanCurve)
-
-
-Is.jordan = is_jordan
+        index -= 1
+    usegments = Future.concatenate(usegments)
+    return usegments
diff --git a/src/shapepy/geometry/piecewise.py b/src/shapepy/geometry/piecewise.py
index 094a30ba..7c23682a 100644
--- a/src/shapepy/geometry/piecewise.py
+++ b/src/shapepy/geometry/piecewise.py
@@ -5,11 +5,12 @@
 from __future__ import annotations
 
 from collections import defaultdict
-from typing import Iterable, Tuple, Union
+from typing import Iterable, Iterator, Tuple, Union
 
 from ..loggers import debug
+from ..rbool import IntervalR1
 from ..scalar.reals import Real
-from ..tools import Is, To, vectorize
+from ..tools import Is, To, pairs
 from .base import IParametrizedCurve
 from .box import Box
 from .point import Point2D
@@ -21,37 +22,39 @@ class PiecewiseCurve(IParametrizedCurve):
     Defines a piecewise curve that is the concatenation of several segments.
     """
 
-    def __init__(
-        self,
-        segments: Iterable[Segment],
-        knots: Union[None, Iterable[Real]] = None,
-    ):
+    def __init__(self, segments: Iterable[Segment]):
         segments = tuple(segments)
         if not all(Is.instance(seg, Segment) for seg in segments):
             raise ValueError("All segments must be instances of Segment")
-        if knots is None:
-            knots = tuple(map(To.rational, range(len(segments) + 1)))
-        else:
-            knots = tuple(sorted(map(To.finite, knots)))
-        for segi, segj in zip(segments, segments[1:]):
-            if segi(1) != segj(0):
-                raise ValueError("Not Continuous curve")
+        knots = list(segments[0].knots)
+        for segi, segj in pairs(segments):
+            knots.append(segj.knots[-1])
+            if segi.knots[-1] != segj.knots[0]:
+                raise ValueError(
+                    f"{segi.domain} and {segj.domain} not consecutive"
+                )
+            pointa = segi(segi.knots[-1])
+            pointb = segj(segj.knots[0])
+            if pointa != pointb:
+                raise ValueError(
+                    f"{segi} not continuous with {segj}: {pointa} != {pointb}"
+                )
+        self.__domain = IntervalR1(knots[0], knots[-1])
         self.__segments = segments
-        self.__knots = knots
+        self.__knots = tuple(knots)
 
     def __str__(self):
-        msgs = []
-        for i, segmenti in enumerate(self.__segments):
-            interval = [self.knots[i], self.knots[i + 1]]
-            msg = f"{interval}: {segmenti}"
-            msgs.append(msg)
-        return r"{" + ", ".join(msgs) + r"}"
+        return r"{" + ", ".join(map(str, self)) + r"}"
+
+    def __repr__(self):
+        return self.__str__()
 
     @property
-    def knots(self) -> Tuple[Real]:
-        """
-        Get the knots of the piecewise curve.
-        """
+    def domain(self):
+        return self.__domain
+
+    @property
+    def knots(self) -> Tuple[Real, ...]:
         return self.__knots
 
     @property
@@ -61,7 +64,7 @@ def length(self) -> Real:
         """
         return sum(seg.length for seg in self)
 
-    def __iter__(self):
+    def __iter__(self) -> Iterator[Segment]:
         yield from self.__segments
 
     def __getitem__(self, index: int) -> Segment:
@@ -94,8 +97,8 @@ def span(self, node: Real) -> Union[int, None]:
         """
         if not Is.real(node):
             raise ValueError
-        if node < self.knots[0] or self.knots[-1] < node:
-            return None
+        if node not in self.domain:
+            raise ValueError(f"Node {node} is not in {self.domain}")
         for i, knot in enumerate(self.knots[1:]):
             if node < knot:
                 return i
@@ -112,6 +115,7 @@ def box(self) -> Box:
             box |= bezier.box()
         return box
 
+    @debug("shapepy.geometry.piecewise")
     def split(self, nodes: Iterable[Real]) -> None:
         """
         Creates an opening in the piecewise curve
@@ -136,43 +140,15 @@ def split(self, nodes: Iterable[Real]) -> None:
             if i not in spansnodes:
                 newsegments.append(segmenti)
                 continue
-            knota, knotb = self.knots[i], self.knots[i + 1]
-            unit_nodes = (
-                (knot - knota) / (knotb - knota) for knot in spansnodes[i]
-            )
-            newsegments += list(segmenti.split(unit_nodes))
+            divisions = sorted(spansnodes[i] | set(segmenti.knots))
+            for ka, kb in pairs(divisions):
+                newsegments.append(segmenti.section([ka, kb]))
         self.__knots = tuple(sorted(list(self.knots) + list(nodes)))
         self.__segments = tuple(newsegments)
 
-    @vectorize(1, 0)
-    def __call__(self, node: float, derivate: int = 0) -> Point2D:
-        index = self.span(node)
-        if index is None:
-            raise ValueError(f"Node {node} is out of bounds")
-        knota, knotb = self.knots[index], self.knots[index + 1]
-        unitparam = (node - knota) / (knotb - knota)
-        segment = self[index]
-        return segment(unitparam, derivate)
+    def eval(self, node: float, derivate: int = 0) -> Point2D:
+        return self[self.span(node)].eval(node, derivate)
 
     def __contains__(self, point: Point2D) -> bool:
         """Tells if the point is on the boundary"""
         return any(point in bezier for bezier in self)
-
-
-def is_piecewise(obj: object) -> bool:
-    """
-    Checks if the parameter is a Piecewise curve
-
-    Parameters
-    ----------
-    obj : The object to be tested
-
-    Returns
-    -------
-    bool
-        True if the obj is a PiecewiseCurve, False otherwise
-    """
-    return Is.instance(obj, PiecewiseCurve)
-
-
-Is.piecewise = is_piecewise
diff --git a/src/shapepy/geometry/segment.py b/src/shapepy/geometry/segment.py
index a5c56fda..6daf41f4 100644
--- a/src/shapepy/geometry/segment.py
+++ b/src/shapepy/geometry/segment.py
@@ -13,16 +13,17 @@
 from __future__ import annotations
 
 from copy import copy
-from typing import Iterable, Optional, Tuple
+from typing import Optional, Tuple, Union
 
 from ..analytic.base import IAnalytic
 from ..analytic.bezier import Bezier
-from ..analytic.tools import find_minimum
+from ..analytic.tools import find_minimum, is_constant
 from ..loggers import debug
-from ..rbool import IntervalR1, from_any
+from ..rbool import IntervalR1, WholeR1, from_any, infimum, supremum
+from ..rbool.tools import is_continuous
 from ..scalar.quadrature import AdaptativeIntegrator, IntegratorFactory
 from ..scalar.reals import Math, Real
-from ..tools import Is, To, pairs, vectorize
+from ..tools import Is, To
 from .base import IParametrizedCurve
 from .box import Box
 from .point import Point2D, cartesian
@@ -34,18 +35,60 @@ class Segment(IParametrizedCurve):
     that contains a bezier curve inside it
     """
 
-    def __init__(self, xfunc: IAnalytic, yfunc: IAnalytic):
+    def __init__(
+        self,
+        xfunc: IAnalytic,
+        yfunc: IAnalytic,
+        *,
+        domain: Union[IntervalR1, WholeR1],
+    ):
         if not Is.instance(xfunc, IAnalytic):
             raise TypeError
         if not Is.instance(yfunc, IAnalytic):
             raise TypeError
+        domain = from_any(domain)
+        if not is_continuous(domain):
+            raise TypeError(f"Domain is not continuous: {domain}")
+        if domain not in (xfunc.domain & yfunc.domain):
+            raise ValueError(
+                f"Given domain must be in {xfunc.domain & yfunc.domain}"
+            )
         self.__length = None
-        self.__knots = (To.rational(0, 1), To.rational(1, 1))
+        self.__domain = domain
+        self.__knots = (infimum(self.domain), supremum(self.domain))
         self.__xfunc = xfunc
         self.__yfunc = yfunc
 
+    @property
+    def domain(self) -> Union[IntervalR1, WholeR1]:
+        return self.__domain
+
+    @property
+    def knots(self) -> Tuple[Real, Real]:
+        return self.__knots
+
+    @property
+    def xfunc(self) -> IAnalytic:
+        """
+        Gives the analytic function x(t) from p(t) = (x(t), y(t))
+        """
+        return self.__xfunc
+
+    @property
+    def yfunc(self) -> IAnalytic:
+        """
+        Gives the analytic function y(t) from p(t) = (x(t), y(t))
+        """
+        return self.__yfunc
+
+    @property
+    def length(self) -> Real:
+        if self.__length is None:
+            self.__length = compute_length(self)
+        return self.__length
+
     def __str__(self) -> str:
-        return f"BS{list(self.knots)}:({self.xfunc}, {self.yfunc})"
+        return f"BS{self.domain}:({self.xfunc}, {self.yfunc})"
 
     def __repr__(self) -> str:
         return str(self)
@@ -53,6 +96,7 @@ def __repr__(self) -> str:
     def __eq__(self, other: Segment) -> bool:
         return (
             Is.instance(other, Segment)
+            and self.domain == other.domain
             and self.xfunc == other.xfunc
             and self.yfunc == other.yfunc
         )
@@ -65,41 +109,13 @@ def __contains__(self, point: Point2D) -> bool:
         deltax = self.xfunc - point.xcoord
         deltay = self.yfunc - point.ycoord
         dist_square = deltax * deltax + deltay * deltay
-        return find_minimum(dist_square, [0, 1]) < 1e-12
+        return find_minimum(dist_square, self.domain) < 1e-12
 
-    @vectorize(1, 0)
-    def __call__(self, node: Real, derivate: int = 0) -> Point2D:
+    def eval(self, node: Real, derivate: int = 0) -> Point2D:
         xcoord = self.xfunc.eval(node, derivate)
         ycoord = self.yfunc.eval(node, derivate)
         return cartesian(xcoord, ycoord)
 
-    @property
-    def xfunc(self) -> IAnalytic:
-        """
-        Gives the analytic function x(t) from p(t) = (x(t), y(t))
-        """
-        return self.__xfunc
-
-    @property
-    def yfunc(self) -> IAnalytic:
-        """
-        Gives the analytic function y(t) from p(t) = (x(t), y(t))
-        """
-        return self.__yfunc
-
-    @property
-    def length(self) -> Real:
-        """
-        The length of the segment
-        """
-        if self.__length is None:
-            self.__length = compute_length(self)
-        return self.__length
-
-    @property
-    def knots(self) -> Tuple[Real, Real]:
-        return self.__knots
-
     def derivate(self, times: Optional[int] = 1) -> Segment:
         """
         Gives the first derivative of the curve
@@ -108,54 +124,43 @@ def derivate(self, times: Optional[int] = 1) -> Segment:
             raise ValueError(f"Times must be integer >= 1, not {times}")
         dxfunc = self.xfunc.derivate(times)
         dyfunc = self.yfunc.derivate(times)
-        return Segment(dxfunc, dyfunc)
+        return Segment(dxfunc, dyfunc, domain=self.domain)
 
     def box(self) -> Box:
         """Returns two points which defines the minimal exterior rectangle
 
         Returns the pair (A, B) with A[0] <= B[0] and A[1] <= B[1]
         """
-        xmin = find_minimum(self.xfunc, [0, 1])
-        xmax = -find_minimum(-self.xfunc, [0, 1])
-        ymin = find_minimum(self.yfunc, [0, 1])
-        ymax = -find_minimum(-self.yfunc, [0, 1])
+        xmin = find_minimum(self.xfunc, self.domain)
+        xmax = -find_minimum(-self.xfunc, self.domain)
+        ymin = find_minimum(self.yfunc, self.domain)
+        ymax = -find_minimum(-self.yfunc, self.domain)
         return Box(cartesian(xmin, ymin), cartesian(xmax, ymax))
 
     def __copy__(self) -> Segment:
         return self.__deepcopy__(None)
 
     def __deepcopy__(self, memo) -> Segment:
-        return Segment(copy(self.xfunc), copy(self.yfunc))
+        return Segment(copy(self.xfunc), copy(self.yfunc), domain=self.domain)
 
     def __invert__(self) -> Segment:
         """
         Inverts the direction of the curve.
         If the curve is clockwise, it becomes counterclockwise
         """
-        composition = Bezier([1, 0])
+        composition = Bezier(
+            [self.knots[-1], self.knots[0]], [self.knots[0], self.knots[-1]]
+        )
         xfunc = self.__xfunc.compose(composition)
         yfunc = self.__yfunc.compose(composition)
-        return Segment(xfunc, yfunc)
+        return Segment(xfunc, yfunc, domain=self.domain)
 
-    def split(self, nodes: Iterable[Real]) -> Tuple[Segment, ...]:
-        """
-        Splits the curve into more segments
-        """
-        nodes = (n for n in nodes if self.knots[0] <= n <= self.knots[-1])
-        nodes = sorted(set(nodes) | set(self.knots))
-        return tuple(self.extract([ka, kb]) for ka, kb in pairs(nodes))
-
-    def extract(self, interval: IntervalR1) -> Segment:
+    def section(self, domain: Union[IntervalR1, WholeR1]) -> Segment:
         """Extracts a subsegment from the given segment"""
-        interval = from_any(interval)
-        if not Is.instance(interval, IntervalR1):
-            raise TypeError
-        knota, knotb = interval[0], interval[1]
-        denom = 1 / (knotb - knota)
-        composition = Bezier([(-knota) * denom, (1 - knota) * denom])
-        nxfunc = self.xfunc.compose(composition)
-        nyfunc = self.yfunc.compose(composition)
-        return Segment(nxfunc, nyfunc)
+        domain = from_any(domain)
+        if domain not in self.domain:
+            raise ValueError(f"Given {domain} not in {self.domain}")
+        return Segment(self.xfunc, self.yfunc, domain=domain)
 
 
 @debug("shapepy.geometry.segment")
@@ -163,34 +168,15 @@ def compute_length(segment: Segment) -> Real:
     """
     Computes the length of the jordan curve
     """
-    domain = (0, 1)
     dpsquare = segment.xfunc.derivate() ** 2 + segment.yfunc.derivate() ** 2
     assert Is.instance(dpsquare, IAnalytic)
-    if dpsquare == dpsquare(0):  # Check if it's constant
-        return (domain[1] - domain[0]) * Math.sqrt(dpsquare(0))
+    if is_constant(dpsquare):  # Check if it's constant
+        knota, knotb = segment.knots
+        return (knotb - knota) * Math.sqrt(dpsquare((knota + knotb) / 2))
     integrator = IntegratorFactory.clenshaw_curtis(3)
     adaptative = AdaptativeIntegrator(integrator, 1e-9, 12)
 
     def function(node):
         return Math.sqrt(dpsquare(node))
 
-    return adaptative.integrate(function, domain)
-
-
-def is_segment(obj: object) -> bool:
-    """
-    Checks if the parameter is a Segment
-
-    Parameters
-    ----------
-    obj : The object to be tested
-
-    Returns
-    -------
-    bool
-        True if the obj is a Segment, False otherwise
-    """
-    return Is.instance(obj, Segment)
-
-
-Is.segment = is_segment
+    return adaptative.integrate(function, segment.domain)
diff --git a/src/shapepy/geometry/transform.py b/src/shapepy/geometry/transform.py
index 21e09742..c6025813 100644
--- a/src/shapepy/geometry/transform.py
+++ b/src/shapepy/geometry/transform.py
@@ -36,10 +36,9 @@ def move(curve: IGeometricCurve, vector: Point2D) -> IGeometricCurve:
     if Is.instance(curve, Segment):
         newxfunc = curve.xfunc + vector.xcoord
         newyfunc = curve.yfunc + vector.ycoord
-        return Segment(newxfunc, newyfunc)
+        return Segment(newxfunc, newyfunc, domain=curve.domain)
     if Is.instance(curve, PiecewiseCurve):
-        newsegs = (move(seg, vector) for seg in curve)
-        return PiecewiseCurve(newsegs, curve.knots)
+        return PiecewiseCurve(move(seg, vector) for seg in curve)
     if not Is.instance(curve, IParametrizedCurve):
         return curve.__class__(move(curve.parametrize(), vector))
     raise NotExpectedError(f"Invalid typo: {type(curve)}")
@@ -70,10 +69,9 @@ def scale(
     if Is.instance(curve, Segment):
         newxfunc = curve.xfunc * (amount if Is.real(amount) else amount[0])
         newyfunc = curve.yfunc * (amount if Is.real(amount) else amount[1])
-        return Segment(newxfunc, newyfunc)
+        return Segment(newxfunc, newyfunc, domain=curve.domain)
     if Is.instance(curve, PiecewiseCurve):
-        newsegs = (scale(seg, amount) for seg in curve)
-        return PiecewiseCurve(newsegs, curve.knots)
+        return PiecewiseCurve(scale(seg, amount) for seg in curve)
     if not Is.instance(curve, IParametrizedCurve):
         return curve.__class__(scale(curve.parametrize(), amount))
     raise NotExpectedError(f"Invalid typo: {type(curve)}")
@@ -104,10 +102,9 @@ def rotate(curve: IGeometricCurve, angle: Angle) -> IGeometricCurve:
         cos, sin = angle.cos(), angle.sin()
         newxfunc = cos * curve.xfunc - sin * curve.yfunc
         newyfunc = sin * curve.xfunc + cos * curve.yfunc
-        return Segment(newxfunc, newyfunc)
+        return Segment(newxfunc, newyfunc, domain=curve.domain)
     if Is.instance(curve, PiecewiseCurve):
-        newsegs = (rotate(seg, angle) for seg in curve)
-        return PiecewiseCurve(newsegs, curve.knots)
+        return PiecewiseCurve(rotate(seg, angle) for seg in curve)
     if not Is.instance(curve, IParametrizedCurve):
         return curve.__class__(rotate(curve.parametrize(), angle))
     raise NotExpectedError(f"Invalid typo: {type(curve)}")
diff --git a/src/shapepy/geometry/unparam.py b/src/shapepy/geometry/unparam.py
index 734ce143..4a4f0f91 100644
--- a/src/shapepy/geometry/unparam.py
+++ b/src/shapepy/geometry/unparam.py
@@ -5,8 +5,9 @@
 from __future__ import annotations
 
 from copy import copy
-from typing import Iterable
+from typing import Iterable, Iterator, Tuple, Union
 
+from ..analytic.polynomial import Polynomial
 from ..scalar.reals import Real
 from ..tools import Is
 from .base import IGeometricCurve
@@ -48,16 +49,6 @@ def length(self) -> Real:
         """
         return self.__segment.length
 
-    @property
-    def start_point(self) -> Point2D:
-        """Gives the start point of the USegment"""
-        return self.__segment(0)
-
-    @property
-    def end_point(self) -> Point2D:
-        """Gives the end point of the USegment"""
-        return self.__segment(1)
-
     def box(self) -> Box:
         """
         Gives the box that encloses the curve
@@ -73,25 +64,122 @@ def __eq__(self, other: IGeometricCurve) -> bool:
             raise TypeError
         segi = self.parametrize()
         segj = other.parametrize()
-        return segi(0) == segj(0) and segi(1) == segj(1)
+        return segi(segi.knots[0]) == segj(segj.knots[0]) and segi(
+            segi.knots[-1]
+        ) == segj(segj.knots[-1])
 
 
 class UPiecewiseCurve(IGeometricCurve):
     """Equivalent to PiecewiseCurve, but ignores the parametrization"""
 
-    def __init__(self, usegments: Iterable[USegment]):
+    def __init__(self, usegments: Iterable[Union[Segment, USegment]]):
         self.__usegments = tuple(usegments)
+        self.__piecewise = None
 
     @property
     def length(self) -> Real:
-        raise NotImplementedError
+        """The length of the curve"""
+        return sum(useg.length for useg in self)
+
+    def __iter__(self) -> Iterator[Union[Segment, USegment]]:
+        """Unparametrized Segments
+
+        When setting, it checks if the points are the same between
+        the junction of two segments to ensure a closed curve
+
+        :getter: Returns the tuple of connected planar beziers, not copy
+        :setter: Sets the segments of the jordan curve
+        :type: tuple[Segment]
+
+        Example use
+        -----------
+
+        >>> from shapepy import JordanCurve
+        >>> vertices = [(0, 0), (4, 0), (0, 3)]
+        >>> jordan = FactoryJordan.polygon(vertices)
+        >>> print(tuple(jordan))
+        (Segment (deg 1), Segment (deg 1), Segment (deg 1))
+        >>> print(tuple(jordan)[0])
+        Planar curve of degree 1 and control points ((0, 0), (4, 0))
+
+        """
+        yield from (
+            useg if Is.instance(useg, USegment) else USegment(useg)
+            for useg in self.__usegments
+        )
+
+    def __len__(self) -> int:
+        return len(self.__usegments)
+
+    def __contains__(self, point: Point2D) -> bool:
+        """Tells if the point is on the boundary"""
+        return any(point in useg for useg in self)
 
     def box(self) -> Box:
-        raise NotImplementedError
+        """The box which encloses the jordan curve
+
+        :return: The box which encloses the jordan curve
+        :rtype: Box
+
+        Example use
+        -----------
+
+        >>> from shapepy import JordanCurve
+        >>> vertices = [(0, 0), (4, 0), (0, 3)]
+        >>> jordan = FactoryJordan.polygon(vertices)
+        >>> jordan.box()
+        Box with vertices (0, 0) and (4, 3)
+
+        """
+        box = None
+        for usegment in self:
+            box |= usegment.box()
+        return box
 
     def parametrize(self) -> PiecewiseCurve:
         """Gives a parametrized curve"""
-        return PiecewiseCurve(useg.parametrize() for useg in self.__usegments)
+        if self.__piecewise is None:
+            newsegs = []
+            for i, usegment in enumerate(self.__usegments):
+                segment = usegment.parametrize()
+                pointa = segment(segment.knots[0])
+                pointb = segment(segment.knots[-1])
+                composition = find_linear_composition(
+                    segment.knots, [i, i + 1]
+                )
+                xfunc = segment.xfunc.compose(composition)
+                yfunc = segment.yfunc.compose(composition)
+                segment = Segment(xfunc, yfunc, domain=[i, i + 1])
+                assert segment(segment.knots[0]) == pointa
+                assert segment(segment.knots[-1]) == pointb
+                newsegs.append(segment)
+            self.__piecewise = PiecewiseCurve(newsegs)
+        return self.__piecewise
+
+
+def find_linear_composition(
+    domain: Tuple[Real, Real], image: Tuple[Real, Real]
+) -> Polynomial:
+    """
+    From a function f(x) defined in [a, b] (image),
+    we want to find a linear function g(t) in the domain [c, d] such
+
+    * h(t) = f(g(t))
+
+    Therefore
+        g(c) = a
+        g(d) = b
+    From hypothesis, g(t) = A + B * t
+        A + B * c = a
+        A + B * d = b
+    Therefore
+        B = (b-a)/(d-c)
+        A = (a*d-c*b)/(d-c)
+    """
+    denom = 1 / (domain[1] - domain[0])
+    const = denom * (image[0] * domain[1] - image[1] * domain[0])
+    slope = denom * (image[1] - image[0])
+    return Polynomial([const, slope])
 
 
 def clean_usegment(usegment: USegment) -> USegment:
diff --git a/src/shapepy/plot/plot.py b/src/shapepy/plot/plot.py
index 81b1535f..4052273d 100644
--- a/src/shapepy/plot/plot.py
+++ b/src/shapepy/plot/plot.py
@@ -31,7 +31,7 @@ def patch_segment(segment: Segment):
     commands = []
     xfunc, yfunc = segment.xfunc, segment.yfunc
     if xfunc.degree <= 1 and yfunc.degree <= 1:
-        vertices.append(segment(1))
+        vertices.append(segment(segment.knots[-1]))
         commands.append(Path.LINETO)
     elif xfunc.degree == 2 and yfunc.degree == 2:
         xfunc: Bezier = segment.xfunc
@@ -50,7 +50,7 @@ def path_shape(connected: ConnectedShape) -> Path:
     commands = []
     for jordan in connected.jordans:
         segments = tuple(useg.parametrize() for useg in jordan)
-        vertices.append(segments[0](0))
+        vertices.append(segments[0](segments[0].knots[0]))
         commands.append(Path.MOVETO)
         for segment in segments:
             verts, comms = patch_segment(segment)
@@ -67,7 +67,7 @@ def path_jordan(jordan: JordanCurve) -> Path:
     Creates the commands for matplotlib to plot the jordan curve
     """
     segments = tuple(useg.parametrize() for useg in jordan)
-    vertices = [segments[0](0)]
+    vertices = [segments[0](segments[0].knots[0])]
     commands = [Path.MOVETO]
     for segment in segments:
         verts, comms = patch_segment(segment)
diff --git a/src/shapepy/rbool/tools.py b/src/shapepy/rbool/tools.py
index df46dcc4..fd64769a 100644
--- a/src/shapepy/rbool/tools.py
+++ b/src/shapepy/rbool/tools.py
@@ -235,3 +235,27 @@ def subset_length(subset: SubSetR1) -> Real:
     if Is.instance(subset, DisjointR1):
         return sum(map(subset_length, subset.intervals))
     return 0
+
+
+def is_bounded(subset: SubSetR1) -> Real:
+    """
+    Tells if the given subset is limited, meaning it does not contain INF
+    """
+    subset = Future.convert(subset)
+    return not Is.instance(subset, WholeR1) and (
+        not (
+            Is.instance(subset, IntervalR1)
+            and not (Math.NEGINF < subset[0] and subset[1] < Math.POSINF)
+        )
+        and not (
+            Is.instance(subset, DisjointR1)
+            and not all(map(is_bounded, subset))
+        )
+    )
+
+
+def is_continuous(subset: SubSetR1) -> Real:
+    """
+    Tells if the given subset is continuous, there's no gaps
+    """
+    return not Is.instance(Future.convert(subset), DisjointR1)
diff --git a/src/shapepy/tools.py b/src/shapepy/tools.py
index 621dfc38..8efb53ab 100644
--- a/src/shapepy/tools.py
+++ b/src/shapepy/tools.py
@@ -113,9 +113,12 @@ def reverse(objs: Iterable[Any]) -> Iterable[Any]:
     return tuple(objs)[::-1]
 
 
+T = TypeVar("T")
+
+
 def pairs(
-    objs: Iterable[Any], /, *, cyclic: bool = False
-) -> Iterable[Tuple[Any, Any]]:
+    objs: Iterable[T], /, *, cyclic: bool = False
+) -> Iterable[Tuple[T, T]]:
     """Gives pairs of the objects in sequence
 
     Example
@@ -135,9 +138,6 @@ class NotExpectedError(Exception):
     """Raised when arrives in a section that were not expected"""
 
 
-T = TypeVar("T")
-
-
 class CyclicContainer(Generic[T]):
     """
     Class that allows checking if there's a circular similarity
diff --git a/tests/analytic/test_bezier.py b/tests/analytic/test_bezier.py
index 2a14f6ba..6f9b09bc 100644
--- a/tests/analytic/test_bezier.py
+++ b/tests/analytic/test_bezier.py
@@ -55,24 +55,24 @@ def test_matrices():
 @pytest.mark.dependency(depends=["test_build", "test_degree", "test_matrices"])
 def test_compare():
     domain = [0, 1]
-    bezier = Bezier([1], domain)
-    assert bezier == Polynomial([1], domain)
+    bezier = Bezier([1], domain=domain)
+    assert bezier == Polynomial([1], domain=domain)
     assert bezier == 1
 
-    bezier = Bezier([1, 1, 1], domain)
-    assert bezier == Polynomial([1], domain)
+    bezier = Bezier([1, 1, 1], domain=domain)
+    assert bezier == Polynomial([1], domain=domain)
     assert bezier == 1
 
-    bezier = Bezier([1, 2], domain)
-    assert bezier == Polynomial([1, 1], domain)
-    bezier = Bezier([1, 2, 3], domain)
-    assert bezier == Polynomial([1, 2], domain)
+    bezier = Bezier([1, 2], domain=domain)
+    assert bezier == Polynomial([1, 1], domain=domain)
+    bezier = Bezier([1, 2, 3], domain=domain)
+    assert bezier == Polynomial([1, 2], domain=domain)
 
     assert bezier != 1
     assert bezier != "asd"
 
     assert Bezier([1, 1, 1]) == Bezier([1])
-    assert Bezier([1, 1], [0, 1]) != Bezier([1], [-1, 2])
+    assert Bezier([1, 1], domain=[0, 1]) != Bezier([1], domain=[-1, 2])
 
 
 @pytest.mark.order(4)
@@ -115,6 +115,18 @@ def test_evaluate():
         assert bezier(0.5) == (a + 2 * b + c) / 4
         assert bezier(1) == c
 
+    bezier = Bezier([10, 20], [0, 1])
+    assert bezier(0) == 10
+    assert bezier(1) == 20
+
+    bezier = Bezier([10, 20], [-1, 1])
+    assert bezier(-1) == 10
+    assert bezier(1) == 20
+
+    bezier = Bezier([10, 20], [-1, 2])
+    assert bezier(-1) == 10
+    assert bezier(2) == 20
+
 
 @pytest.mark.order(4)
 @pytest.mark.dependency(depends=["test_build", "test_degree", "test_evaluate"])
diff --git a/tests/analytic/test_polynomial.py b/tests/analytic/test_polynomial.py
index 6afd97ca..950613a8 100644
--- a/tests/analytic/test_polynomial.py
+++ b/tests/analytic/test_polynomial.py
@@ -148,7 +148,7 @@ def test_print():
     assert str(poly) == "1 + 2 * t + 3 * t^2"
     repr(poly)
 
-    poly = Polynomial([1, 2, 3], [0, 1])
+    poly = Polynomial([1, 2, 3], domain=[0, 1])
     assert repr(poly) == "[0, 1]: 1 + 2 * t + 3 * t^2"
 
 
diff --git a/tests/bool2d/test_bool_overlap.py b/tests/bool2d/test_bool_overlap.py
index 8519f05c..d32dcb3f 100644
--- a/tests/bool2d/test_bool_overlap.py
+++ b/tests/bool2d/test_bool_overlap.py
@@ -7,6 +7,7 @@
 from shapepy.bool2d.base import EmptyShape, WholeShape
 from shapepy.bool2d.config import set_auto_clean
 from shapepy.bool2d.primitive import Primitive
+from shapepy.loggers import enable_logger
 
 
 @pytest.mark.order(43)
@@ -347,11 +348,12 @@ def test_sub(self):
         left = Primitive.circle(radius=3, center=(-10, 0))
         right = Primitive.circle(radius=3, center=(10, 0))
         with set_auto_clean(False):
-            shape = big - small | left ^ right
-            assert shape - shape is EmptyShape()
-            assert shape - (~shape) == shape
-            assert (~shape) - shape == ~shape
-            assert (~shape) - (~shape) is EmptyShape()
+            with enable_logger("shapepy.bool2d"):
+                shape = big - small  # | left ^ right
+                assert shape - shape is EmptyShape()
+                assert shape - (~shape) == shape
+                assert (~shape) - shape == ~shape
+                assert (~shape) - (~shape) is EmptyShape()
 
     @pytest.mark.order(43)
     @pytest.mark.timeout(40)
diff --git a/tests/bool2d/test_shape.py b/tests/bool2d/test_shape.py
index e56cc90a..ce4db3ef 100644
--- a/tests/bool2d/test_shape.py
+++ b/tests/bool2d/test_shape.py
@@ -45,13 +45,13 @@ def test_centered_rectangular(self):
                     rectangular.jordan, expx, expy
                 )
                 if expx % 2 or expy % 2:
-                    good = 0
+                    assert abs(test) < 1e-8
                 else:
                     good = width ** (expx + 1)
                     good *= height ** (expy + 1)
                     good /= (1 + expx) * (1 + expy)
                     good /= 2 ** (expx + expy)
-                assert abs(test - good) < 1e-9
+                    assert abs(test - good) < 1e-9 * abs(good)
 
     @pytest.mark.order(25)
     @pytest.mark.timeout(10)
@@ -78,7 +78,7 @@ def test_noncenter_rectangular(self):
                     center[1] - height / 2
                 ) ** (expy + 1)
                 good /= (1 + expx) * (1 + expy)
-                assert abs(test - good) < 1e-9 * abs(good)
+                assert abs(test - good) < 1e-9 * max(1, abs(good))
 
     @pytest.mark.order(25)
     @pytest.mark.timeout(10)
@@ -90,20 +90,20 @@ def test_noncenter_rectangular(self):
         ]
     )
     def test_centered_rombo(self):
-        width, height = 3, 5
+        width, height = 1, 2
         rombo = Primitive.regular_polygon(4).scale((width, height))
         nx, ny = 5, 5
         for expx in range(nx):
             for expy in range(ny):
                 test = IntegrateJordan.polynomial(rombo.jordan, expx, expy)
                 if expx % 2 or expy % 2:
-                    good = 0
+                    assert abs(test) < 1e-6
                 else:
                     good = 4 * width ** (expx + 1) * height ** (expy + 1)
                     good *= math.factorial(expx) * math.factorial(expy)
                     good /= math.factorial(expx + expy)
                     good /= (1 + expx + expy) * (2 + expx + expy)
-                assert abs(test - good) < 1e-9
+                    assert abs(test - good) < 1e-9 * max(1, abs(good))
 
     @pytest.mark.order(25)
     @pytest.mark.timeout(10)
diff --git a/tests/geometry/test_jordan_polygon.py b/tests/geometry/test_jordan_polygon.py
index d7ad4383..58d29c26 100644
--- a/tests/geometry/test_jordan_polygon.py
+++ b/tests/geometry/test_jordan_polygon.py
@@ -8,7 +8,6 @@
 import pytest
 
 from shapepy.geometry.factory import FactoryJordan
-from shapepy.geometry.jordancurve import clean_jordan
 from shapepy.geometry.transform import move, rotate, scale
 from shapepy.scalar.angle import degrees, radians
 
@@ -152,7 +151,17 @@ def test_intersection(self):
         vertices1 = [(-1, 0), (1, 2), (3, 0), (1, -2)]
         square1 = FactoryJordan.polygon(vertices1)
         param0 = square0.parametrize()
+        assert param0(0) == (1, 0)
+        assert param0(1) == (-1, 2)
+        assert param0(2) == (-3, 0)
+        assert param0(3) == (-1, -2)
+        assert param0(4) == (1, 0)
         param1 = square1.parametrize()
+        assert param1(0) == (-1, 0)
+        assert param1(1) == (1, 2)
+        assert param1(2) == (3, 0)
+        assert param1(3) == (1, -2)
+        assert param1(4) == (-1, 0)
 
         inters = param0 & param0
         assert inters.all_subsets[id(param0)] == [0, 4]
@@ -162,9 +171,6 @@ def test_intersection(self):
         assert inters.all_knots[id(param1)] == {0, 1, 2, 3, 4}
 
         inters = param0 & param1
-        print(param0)
-        print(param1)
-        print(inters)
         assert inters.all_subsets[id(param0)] == {0.5, 3.5}
         assert inters.all_knots[id(param0)] == {0, 0.5, 1, 2, 3, 3.5, 4}
         assert inters.all_subsets[id(param1)] == {0.5, 3.5}
@@ -443,28 +449,24 @@ def test_print(self):
     def test_clean(self):
         verticesa = [(-1, 0), (0, 0), (1, 0), (0, 1)]
         jordana = FactoryJordan.polygon(verticesa)
-        jordana = clean_jordan(jordana)
         verticesb = [(-1, 0), (1, 0), (0, 1)]
         jordanb = FactoryJordan.polygon(verticesb)
         assert jordana == jordanb
 
         verticesa = [(-1.0, 0.0), (0.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
         jordana = FactoryJordan.polygon(verticesa)
-        jordana = clean_jordan(jordana)
         verticesb = [(-1.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
         jordanb = FactoryJordan.polygon(verticesb)
         assert jordana == jordanb
 
         verticesa = [(0, 0), (1, 0), (0, 1), (-1, 0)]
         jordana = FactoryJordan.polygon(verticesa)
-        jordana = clean_jordan(jordana)
         verticesb = [(-1, 0), (1, 0), (0, 1)]
         jordanb = FactoryJordan.polygon(verticesb)
         assert jordana == jordanb
 
         verticesa = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (-1.0, 0.0)]
         jordana = FactoryJordan.polygon(verticesa)
-        jordana = clean_jordan(jordana)
         verticesb = [(-1.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
         jordanb = FactoryJordan.polygon(verticesb)
         assert jordana == jordanb
diff --git a/tests/geometry/test_piecewise.py b/tests/geometry/test_piecewise.py
index 529c87af..b05d1717 100644
--- a/tests/geometry/test_piecewise.py
+++ b/tests/geometry/test_piecewise.py
@@ -30,9 +30,10 @@ def test_build():
         ((1, 1), (0, 1)),
         ((0, 1), (0, 0)),
     ]
-    knots = range(len(points) + 1)
-    segments = tuple(map(FactorySegment.bezier, points))
-    PiecewiseCurve(segments, knots)
+    segments = (
+        FactorySegment.bezier(pts, [i, i + 1]) for i, pts in enumerate(points)
+    )
+    PiecewiseCurve(segments)
 
 
 @pytest.mark.order(14)
@@ -44,9 +45,10 @@ def test_box():
         ((1, 1), (0, 1)),
         ((0, 1), (0, 0)),
     ]
-    knots = range(len(points) + 1)
-    segments = tuple(map(FactorySegment.bezier, points))
-    piecewise = PiecewiseCurve(segments, knots)
+    segments = (
+        FactorySegment.bezier(pts, [i, i + 1]) for i, pts in enumerate(points)
+    )
+    piecewise = PiecewiseCurve(segments)
     box = piecewise.box()
     assert box.lowpt == (0, 0)
     assert box.toppt == (1, 1)
@@ -61,9 +63,10 @@ def test_evaluate():
         ((1, 1), (0, 1)),
         ((0, 1), (0, 0)),
     ]
-    knots = range(len(points) + 1)
-    segments = tuple(map(FactorySegment.bezier, points))
-    piecewise = PiecewiseCurve(segments, knots)
+    segments = (
+        FactorySegment.bezier(pts, [i, i + 1]) for i, pts in enumerate(points)
+    )
+    piecewise = PiecewiseCurve(segments)
     assert piecewise(0) == (0, 0)
     assert piecewise(1) == (1, 0)
     assert piecewise(2) == (1, 1)
@@ -85,9 +88,10 @@ def test_print():
         ((1, 1), (0, 1)),
         ((0, 1), (0, 0)),
     ]
-    knots = range(len(points) + 1)
-    segments = tuple(map(FactorySegment.bezier, points))
-    piecewise = PiecewiseCurve(segments, knots)
+    segments = (
+        FactorySegment.bezier(pts, [i, i + 1]) for i, pts in enumerate(points)
+    )
+    piecewise = PiecewiseCurve(segments)
     str(piecewise)
     repr(piecewise)
 
diff --git a/tests/geometry/test_segment.py b/tests/geometry/test_segment.py
index 6d4bcbdc..093b332e 100644
--- a/tests/geometry/test_segment.py
+++ b/tests/geometry/test_segment.py
@@ -34,6 +34,18 @@ def test_build():
     FactorySegment.bezier([(0, 0), (1, 0), (0, 1)])
 
 
+@pytest.mark.order(13)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin", "test_build"])
+def test_invert():
+    bezier = FactorySegment.bezier([(0, 0), (1, 0), (0, 1)])
+    invbez = ~bezier
+    assert bezier.domain == invbez.domain
+    assert bezier.knots == invbez.knots
+    assert bezier(bezier.knots[0]) == invbez(invbez.knots[-1])
+    assert bezier(bezier.knots[-1]) == invbez(invbez.knots[0])
+
+
 class TestDerivate:
     @pytest.mark.order(13)
     @pytest.mark.dependency(
@@ -127,11 +139,10 @@ def test_middle(self):
         half = Fraction(1, 2)
         points = [(0, 0), (1, 0)]
         curve = FactorySegment.bezier(points)
-        curvea = FactorySegment.bezier([(0, 0), (half, 0)])
-        curveb = FactorySegment.bezier([(half, 0), (1, 0)])
-        assert curve.extract([0, half]) == curvea
-        assert curve.extract([half, 1]) == curveb
-        assert curve.split([half]) == (curvea, curveb)
+        curvea = FactorySegment.bezier([(0, 0), (half, 0)], [0, 0.5])
+        curveb = FactorySegment.bezier([(half, 0), (1, 0)], [0.5, 1])
+        assert curve.section([0, half]) == curvea
+        assert curve.section([half, 1]) == curveb
 
         test = curvea | curveb
         assert test == curve
@@ -166,6 +177,7 @@ def test_print():
     depends=[
         "test_begin",
         "test_build",
+        "test_invert",
         "TestDerivate::test_end",
         "TestContains::test_end",
         "TestSplitUnite::test_end",
diff --git a/tests/rbool/test_limits.py b/tests/rbool/test_limits.py
index e8c856e2..8b1fb2db 100644
--- a/tests/rbool/test_limits.py
+++ b/tests/rbool/test_limits.py
@@ -1,6 +1,7 @@
 import pytest
 
-from shapepy.rbool import infimum, maximum, minimum, supremum
+from shapepy.rbool import from_any, infimum, maximum, minimum, supremum
+from shapepy.rbool.tools import is_bounded, is_continuous
 from shapepy.scalar.reals import Math
 
 
@@ -61,3 +62,85 @@ def test_inf_min_max_sup():
     assert minimum(subset) is None
     assert maximum(subset) is None
     assert supremum(subset) == 10
+
+
+@pytest.mark.order(17)
+@pytest.mark.timeout(1)
+@pytest.mark.dependency(
+    depends=[
+        "tests/rbool/test_build.py::test_all",
+        "tests/rbool/test_convert.py::test_all",
+    ],
+    scope="session",
+)
+def test_is_bounded():
+    subset = from_any(r"{}")
+    assert is_bounded(subset)
+
+    subset = from_any(r"(-inf, +inf)")
+    assert not is_bounded(subset)
+
+    subset = from_any(r"{1}")
+    assert is_bounded(subset)
+
+    subset = from_any(r"{1, 3}")
+    assert is_bounded(subset)
+
+    subset = from_any(r"[-10, 10]")
+    assert is_bounded(subset)
+
+    subset = from_any(r"(-10, 10]")
+    assert is_bounded(subset)
+
+    subset = from_any(r"[-10, 10)")
+    assert is_bounded(subset)
+
+    subset = from_any(r"(-10, 10)")
+    assert is_bounded(subset)
+
+    subset = from_any(r"(-inf, 10)")
+    assert not is_bounded(subset)
+
+    subset = from_any(r"(-10, +inf)")
+    assert not is_bounded(subset)
+
+
+@pytest.mark.order(17)
+@pytest.mark.timeout(1)
+@pytest.mark.dependency(
+    depends=[
+        "tests/rbool/test_build.py::test_all",
+        "tests/rbool/test_convert.py::test_all",
+    ],
+    scope="session",
+)
+def test_is_continuous():
+    subset = from_any(r"{}")
+    assert is_continuous(subset)
+
+    subset = from_any(r"(-inf, +inf)")
+    assert is_continuous(subset)
+
+    subset = from_any(r"{1}")
+    assert is_continuous(subset)
+
+    subset = from_any(r"{1, 3}")
+    assert not is_continuous(subset)
+
+    subset = from_any(r"[-10, 10]")
+    assert is_continuous(subset)
+
+    subset = from_any(r"(-10, 10]")
+    assert is_continuous(subset)
+
+    subset = from_any(r"[-10, 10)")
+    assert is_continuous(subset)
+
+    subset = from_any(r"(-10, 10)")
+    assert is_continuous(subset)
+
+    subset = from_any(r"(-inf, 10)")
+    assert is_continuous(subset)
+
+    subset = from_any(r"(-10, +inf)")
+    assert is_continuous(subset)