import math
import numbers
import traceback
from math import fabs
from typing import Union, Tuple, Optional, List

from geogst.core.mathematics.rotations import RotationAxis
from geogst.core.geometries.shapes.points import *
from geogst.core.utils.inspections.errors import Error
from geogst.core.utils.inspections.functions import caller_name
from geogst.core.mathematics.defaults import PRACTICAL_MIN_DIST, MAX_SCALAR_VALUE, MIN_SCALAR_VALUE
from geogst.core.mathematics.utils import are_close
from geogst.core.mathematics.vectors2d import Vect2D
from geogst.core.mathematics.vectors3d import Vect3D
from geogst.core.utils.types import check_type


class Segment(Shape):

    proper_space = 1

    @property
    def embedding_space(self) -> int:
        return self._coords.shape[1]

    def __init__(self,
                 start_pt: Point,
                 end_pt: Point
                 ):
        """
        Creates a segment instance.

        :param start_pt: the start point.
        :param end_pt: the end point.
        :return: the new segment instance.
        """

        super(Segment, self).__init__()

        self._coords = np.vstack((start_pt.coords, end_pt.coords))

    @property
    def coords(self) -> np.ndarray:

        return self._coords

    @property
    def start_pt(self):
        """Return start point"""

        return Point(*tuple(self._coords[0]))

    @property
    def end_pt(self):
        """Return end point"""

        return Point(*tuple(self._coords[1]))

    def __repr__(self) -> str:
        """
        Represents a Segment instance.

        :return: the Segment representation.
        """

        return f"Segment(start_pt={self.start_pt}, end_pt={self.end_pt})"

    def __iter__(self):
        """
        Return the elements of a Segment, i.e., start and end point.
        """

        return (Point(*tuple(self.coords[i])) for i in range(2))

    def clone(self) -> 'Segment':

        return Segment(self.start_pt, self.end_pt)

    def as_2d(self) -> 'Segment':
        """
        Creates a 2D copy of the original segment.

        Examples:
         >>> s3 = Segment(Point(1,2,3), Point(4, 5, 6))
         >>> s3.as_2d()
         Segment(start_pt=Point([1. 2.]), end_pt=Point([4. 5.]))
        """

        return Segment(self.start_pt.as_point2d(), self.end_pt.as_point2d())

    @property
    def length(self) -> numbers.Real:
        """Returns the length of the segment.

        Examples:
         >>> Segment(Point(0, 0, 0), Point(3, 4, 5)).length
         7.0710678118654755
        """

        return float( np.linalg.norm((self.end_pt - self.start_pt).coords) )

    def bounding_box(self) -> Union[None, Tuple[numbers.Real, numbers.Real, numbers.Real, numbers.Real]]:
        """Returns the bounding box of the segment."""

        try:

            xmin, xmax = self.x_range()
            ymin, ymax = self.y_range()

            return (xmin, ymin, xmax, ymax)

        except Exception as e:

            print(f"\n{e!r}")
            return None

    def length_2d(self) -> numbers.Real:
        """
        Returns the horizontal (2D) length of the segment.

        :return: the horizontal length of the segment.
        """

        return np.linalg.norm((self.end_pt - self.start_pt).coords[:2])

    def length_3d(self) -> numbers.Real:
        """
        Returns the 3D length of the segment.

        :return: the 3D length of the segment.
        """

        return np.linalg.norm((self.end_pt - self.start_pt).coords[:3])

    def area(self):
        """Calculate shape area"""

        return 0.0

    def delta_x(self) -> numbers.Real:
        """
        X delta between segment end point and start point.

        :return: the horizontal, x-parallel distance between segment end point and start point.
        """

        return self.end_pt.x - self.start_pt.x

    def delta_y(self) -> numbers.Real:
        """
        Y delta between segment end point and start point.

        :return: the horizontal, y-parallel distance between segment end point and start point.
        """

        return self.end_pt.y - self.start_pt.y

    def delta_z(self) -> Union[None, numbers.Real]:
        """
        Z delta between segment end point and start point.

        :return: the vertical distance between segment end point and start point.
        """

        try:
            return self.end_pt.z - self.start_pt.z
        except:
            return None

    def x_range(self) -> Tuple[numbers.Real, numbers.Real]:

        if self.start_pt.x < self.end_pt.x:
            return self.start_pt.x, self.end_pt.x
        else:
            return self.end_pt.x, self.start_pt.x

    def y_range(self) -> Tuple[numbers.Real, numbers.Real]:

        if self.start_pt.y < self.end_pt.y:
            return self.start_pt.y, self.end_pt.y
        else:
            return self.end_pt.y, self.start_pt.y

    def z_range(self) -> Tuple[numbers.Real, numbers.Real]:

        if self.start_pt.z < self.end_pt.z:
            return self.start_pt.z, self.end_pt.z
        else:
            return self.end_pt.z, self.start_pt.z

    def left_norm_vers2d(self) -> Vect2D:
        """
        Returns the left horizontal normal versor.

        :return: the left horizontal normal versor.
        """

        return Vect2D(
            x=-self.as_versor2d().y,
            y=self.as_versor2d().x
        )

    def right_norm_vers2d(self) -> Vect2D:
        """
        Returns the right horizontal normal versor.

        :return: the right horizontal normal versor.
        """

        return Vect2D(
            x=self.as_versor2d().y,
            y=-self.as_versor2d().x
        )

    def as_vector2d(self) -> Vect2D:
        """
        Convert a segment to a vector.
        """

        return Vect2D(
            x=self.delta_x(),
            y=self.delta_y()
        )

    def as_versor2d(self) -> Optional[Vect2D]:
        """
        Convert a segment to an optional versor.
        """

        return self.as_vector2d().versor()

    def as_vector3d(self) -> Vect3D:
        """
        Convert a segment to a vector.
        """

        return Vect3D(
            x=self.delta_x(),
            y=self.delta_y(),
            z=self.delta_z()
        )

    def as_versor3d(self) -> Vect3D:
        """
        Convert a segment to a versor.
        """

        vers, _ = self.as_vector3d().to_versor()

        return vers

    def conn_to_other(self,
                      another: 'Segment',
                      tol: numbers.Real = 1e-12
                      ) -> bool:
        """
        Check whether the first segment is sequentially connected to the second one.

        :param another: a segment to check for.
        :param tol: tolerance for distance between points.
        :return: whether the first segment is sequentially connected to the second one.

        Examples:
          >>> s1 = Segment(Point(0,0), Point(1,0))
          >>> s2 = Segment(Point(1,0), Point(2,0))
          >>> s1.conn_to_other(s2)
          True
          >>> s1 = Segment(Point(0,0,0), Point(1,0,0))
          >>> s2 = Segment(Point(1,0,0), Point(2,0,0))
          >>> s1.conn_to_other(s2)
          True
        """

        return self.end_pt.is_coincident(
            other=another.start_pt,
            tolerance=tol)

    def point_at_factor(self,
                        scale_factor: numbers.Real
                        ) -> Point:
        """
        Returns a point aligned with the segment
        and lying at given scale factor, where 1 is segment length
        and 0 is segment start.

        :param scale_factor: the scale factor, where 1 is the segment length.
        :return: Point at scale factor

        Examples:
          >>> s = Segment(Point(0,0), Point(1,0))
          >>> s.point_at_factor(0)
          Point([0. 0.])
          >>> s.point_at_factor(0.5)
          Point([0.5 0. ])
          >>> s.point_at_factor(1)
          Point([1. 0.])
          >>> s.point_at_factor(-1)
          Point([-1. 0.])
          >>> s.point_at_factor(-2)
          Point([-2.  0.])
          >>> s.point_at_factor(2)
          Point([2. 0.])
          >>> s = Segment(Point(0,0), Point(1,1))
          >>> s.point_at_factor(0.5)
          Point([0.5 0.5])
          >>> s = Segment(Point(0,0), Point(4,0))
          >>> s.point_at_factor(7.5)
          Point([30.  0.])
          >>> s = Segment(Point(0,0,0), Point(1,0,0))
          >>> s.point_at_factor(0)
          Point([0. 0. 0.])
          >>> s.point_at_factor(0.5)
          Point([0.5 0.  0. ])
          >>> s.point_at_factor(1)
          Point([1. 0. 0.])
          >>> s.point_at_factor(-1)
          Point([-1.  0.  0.])
          >>> s.point_at_factor(-2)
          Point([-2.  0.  0.])
          >>> s.point_at_factor(2)
          Point([2. 0. 0.])
          >>> s = Segment(Point(0,0,0), Point(0,0,1))
          >>> s.point_at_factor(0)
          Point([0. 0. 0.])
          >>> s.point_at_factor(0.5)
          Point([0.  0.  0.5])
          >>> s.point_at_factor(1)
          Point([0. 0. 1.])
          >>> s.point_at_factor(-1)
          Point([ 0.  0. -1.])
          >>> s.point_at_factor(-2)
          Point([ 0.  0. -2.])
          >>> s.point_at_factor(2)
          Point([0. 0. 2.])
          >>> s = Segment(Point(0,0,0), Point(1,1,1))
          >>> s.point_at_factor(0.5)
          Point([0.5 0.5 0.5])
          >>> s = Segment(Point(0,0,0), Point(4,0,0))
          >>> s.point_at_factor(7.5)
          Point([30.  0.  0.])
        """

        dx = self.delta_x() * scale_factor
        dy = self.delta_y() * scale_factor
        dz = None if self.delta_z() is None else self.delta_z() * scale_factor

        if dz is None:
            return Point(
                self.start_pt.x + dx,
                self.start_pt.y + dy
            )
        else:
            return Point(
                self.start_pt.x + dx,
                self.start_pt.y + dy,
                self.start_pt.z + dz
            )

    def midpoint(self) -> Point:
        """
        Returns the mid-point of the segment.

        Examples:
         >>> Segment(Point(1, 1), Point(3, 3)).midpoint()
         Point([2. 2.])
         >>> Segment(Point(0, 0), Point(0, 9)).midpoint()
         Point([0. 4.5])
         >>> Segment(Point(0, 0), Point(0, 0)).midpoint()
         Point([0. 0.])
         >>> Segment(Point(0, 0), Point(0, 20)).midpoint()
         Point([ 0. 10.])
        """

        return self.point_at_factor(0.5)

    def contains_pt(self,
                    pt: Point
                    ) -> bool:
        """
        Checks whether a point is contained in a segment.

        :param pt: the point for which to check containement.
        :return: whether a point is contained in a segment.

        Examples:
          >>> segment = Segment(Point(0, 0), Point(1, 0))
          >>> segment.contains_pt(Point(0, 0))
          True
          >>> segment.contains_pt(Point(1, 0))
          True
          >>> segment.contains_pt(Point(0.5, 0))
          True
          >>> segment.contains_pt(Point(0.5, 0.00001))
          False
          >>> segment.contains_pt(Point(1.00001, 0))
          False
          >>> segment.contains_pt(Point(0.000001, 0))
          True
          >>> segment.contains_pt(Point(-0.000001, 0))
          False
          >>> segment.contains_pt(Point(0.5, 1000))
          False
          >>> segment = Segment(Point(0, 0), Point(0, 1))
          >>> segment.contains_pt(Point(0, 0))
          True
          >>> segment.contains_pt(Point(0, 0.5))
          True
          >>> segment.contains_pt(Point(0, 1))
          True
          >>> segment.contains_pt(Point(0, 1.5))
          False
          >>> segment = Segment(Point(0, 0), Point(1, 1))
          >>> segment.contains_pt(Point(0.5, 0.5))
          True
          >>> segment.contains_pt(Point(1, 1))
          True
          >>> segment = Segment(Point(1, 2), Point(9, 8))
          >>> segment.contains_pt(segment.point_at_factor(0.745))
          True
          >>> segment.contains_pt(segment.point_at_factor(1.745))
          False
          >>> segment.contains_pt(segment.point_at_factor(-0.745))
          False
          >>> segment.contains_pt(segment.point_at_factor(0))
          True
          >>> segment = Segment(Point(0, 0, 0), Point(1, 0, 0))
          >>> segment.contains_pt(Point(0, 0, 0))
          True
          >>> segment.contains_pt(Point(1, 0, 0))
          True
          >>> segment.contains_pt(Point(0.5, 0, 0))
          True
          >>> segment.contains_pt(Point(0.5, 0.00001, 0))
          False
          >>> segment.contains_pt(Point(0.5, 0, 0.00001))
          False
          >>> segment.contains_pt(Point(1.00001, 0, 0))
          False
          >>> segment.contains_pt(Point(0.000001, 0, 0))
          True
          >>> segment.contains_pt(Point(-0.000001, 0, 0))
          False
          >>> segment.contains_pt(Point(0.5, 1000, 1000))
          False
          >>> segment = Segment(Point(0, 0, 0), Point(0, 1, 0))
          >>> segment.contains_pt(Point(0, 0, 0))
          True
          >>> segment.contains_pt(Point(0, 0.5, 0))
          True
          >>> segment.contains_pt(Point(0, 1, 0))
          True
          >>> segment.contains_pt(Point(0, 1.5, 0))
          False
          >>> segment = Segment(Point(0, 0, 0), Point(1, 1, 1))
          >>> segment.contains_pt(Point(0.5, 0.5, 0.5))
          True
          >>> segment.contains_pt(Point(1, 1, 1))
          True
          >>> segment = Segment(Point(1,2,3), Point(9,8,2))
          >>> segment.contains_pt(segment.point_at_factor(0.745))
          True
          >>> segment.contains_pt(segment.point_at_factor(1.745))
          False
          >>> segment.contains_pt(segment.point_at_factor(-0.745))
          False
          >>> segment.contains_pt(segment.point_at_factor(0))
          True
        """

        segment_length = self.length
        length_startpt_pt = self.start_pt.distance(pt)
        length_endpt_pt = self.end_pt.distance(pt)

        return are_close(
            a=segment_length,
            b=length_startpt_pt + length_endpt_pt
        )

    def point2d_at_hor_dist(self,
                            distance: numbers.Real
                            ) -> Point:
        """
        Returns a point aligned with the segment
        and lying at given distance from the segment start.

        :param distance: the distance from segment start
        :return: point at provided distance from segment start

        Examples:
          >>> s = Segment(Point(0,0), Point(1,0))
          >>> s.point2d_at_hor_dist(0)
          Point([0. 0.])
          >>> s.point2d_at_hor_dist(0.5)
          Point([0.5 0. ])
          >>> s.point2d_at_hor_dist(1)
          Point([1. 0.])
          >>> s.point2d_at_hor_dist(2)
          Point([2. 0.])
          >>> s = Segment(Point(0,0), Point(3, 4))
          >>> s.point2d_at_hor_dist(0)
          Point([0. 0.])
          >>> s.point2d_at_hor_dist(2.5)
          Point([1.5 2. ])
          >>> s.point2d_at_hor_dist(5)
          Point([3. 4.])
          >>> s.point2d_at_hor_dist(10)
          Point([6. 8.])
        """

        scale_factor_2d = distance / self.length_2d()

        dx = self.delta_x() * scale_factor_2d
        dy = self.delta_y() * scale_factor_2d

        return Point(
            self.start_pt.x + dx,
            self.start_pt.y + dy
        )

    def densify_as_line2d(self,
                          densify_distance: numbers.Real
                          ) -> 'Ln':
        """
        Densify a segment by adding additional points
        separated a distance equal to densify_distance.
        The result is no longer a Segment instance, instead it is a Ln instance.

        :param densify_distance: the densify distance
        :return: a line
        """

        return Ln.from_points(*self.densify_as_points2d(densify_distance=densify_distance))

    def densify_as_pts2d(self,
        densify_distance: numbers.Real,
        start_offset: numbers.Real = 0.0
    ) -> Tuple[Union[None, List[Point]], Error]:
        """
        Densify a segment as a list of points, by using the provided densify distance.

        :param densify_distance: the densify distance
        :param start_offset: the offset from the start point from which to densify.
        :return: an optional list of points and the error status.
        """

        try:

            length2d = self.length_2d()

            vers_2d = self.as_versor2d()

            step_vector = vers_2d.scale(densify_distance)

            # use the 2D start point as the fixed reference for the distance
            # check below: self.start_pt may be 3D (for a 3D segment), while
            # new_pt is always 2D, and comparing the two would otherwise raise
            # a shape-mismatch error
            start_pt_2d = self.start_pt.as_point2d()

            if start_offset == 0.0:
                start_point = start_pt_2d
            else:
                start_point = self.point2d_at_hor_dist(start_offset)

            pts = [start_point]

            n = 0
            while True:
                n += 1
                shift_vector = step_vector.scale(n)
                new_pt = start_point.shift(shift_vector.x, shift_vector.y)
                distance = start_pt_2d.distance(new_pt)
                if distance >= length2d:
                    break
                pts.append(new_pt)

            pts.append(self.end_pt.as_point2d())

            return pts, Error()

        except Exception as e:

            return None, Error(
                True,
                caller_name(),
                e,
                traceback.format_exc()
            )

    def densify_as_points2d(self,
                            densify_distance: numbers.Real,
                            start_offset: numbers.Real = 0.0
                            ) -> List[Point]:
        """
        DEPRECATED: use 'densify_as_pts2d'.

        Densify a segment as a list of points, by using the provided densify distance.

        :param densify_distance: the densify distance
        :param start_offset: the offset from the start point from which to densify.
        :return: a list of points
        """

        length2d = self.length_2d()

        vers_2d = self.as_versor2d()

        step_vector = vers_2d.scale(densify_distance)

        # use the 2D start point as the fixed reference for the distance
        # check below: self.start_pt may be 3D (for a 3D segment), while
        # new_pt is always 2D, and comparing the two would otherwise raise
        # a shape-mismatch error
        start_pt_2d = self.start_pt.as_point2d()

        if start_offset == 0.0:
            start_point = start_pt_2d
        else:
            start_point = self.point2d_at_hor_dist(start_offset)

        pts = [start_point]

        n = 0
        while True:
            n += 1
            shift_vector = step_vector.scale(n)
            new_pt = start_point.shift(shift_vector.x, shift_vector.y)
            distance = start_pt_2d.distance(new_pt)
            if distance >= length2d:
                break
            pts.append(new_pt)

        pts.append(self.end_pt.as_point2d())

        return pts

    def point_projection(self,
        point: Point
    ) -> Union[None, Point]:
        """
        Return the point projection on the segment.

        Examples:
          >>> s = Segment(start_pt=Point(0,0,0), end_pt=Point(1,0,0))
          >>> p = Point(0.5, 1, 4)
          >>> s.point_projection(p)
          Point([0.5 0.  0. ])
          >>> s = Segment(start_pt=Point(0,0,0), end_pt=Point(4,0,0))
          >>> p = Point(7.5, 19.2, -14.72)
          >>> s.point_projection(p)
          Point([7.5 0.  0. ])
        """

        other_segment = Segment(
            self.start_pt,
            point
        )

        scale_factor = self.as_vector3d().scalar_projection(other_segment.as_vector3d()) / self.length

        return self.point_at_factor(scale_factor)

    def distance_to_point(self,
                          point: Point
                          ) -> Union[None, numbers.Real]:
        """
        Returns the point distance to the segment.

        :param point: the point to calculate the distance with
        :return: the distance of the point to the segment

        Examples:
          >>> s = Segment(Point(0,0,0), Point(0,0,4))
          >>> s.distance_to_point(Point(-17.2, 0.0, -49))
          17.2
          >>> s.distance_to_point(Point(-17.2, 1.22, -49))
          17.24321315764553
        """

        pt_proj = self.point_projection(point)

        return point.distance(pt_proj)

    def segment_2d_m(self) -> Optional[numbers.Real]:

        denom = self.end_pt.x - self.start_pt.x

        if denom == 0.0:
            return None

        return (self.end_pt.y - self.start_pt.y) / denom

    def segment_2d_p(self) -> Optional[numbers.Real]:

        s2d_m = self.segment_2d_m()

        if s2d_m is None:
            return None

        return self.start_pt.y - s2d_m * self.start_pt.x

    def intersect_segments2d(
            self,
            other: 'Segment',
            tol: numbers.Real = PRACTICAL_MIN_DIST
    ) -> Union[None, Point, 'Segment']:
        """
        Determines the optional point or segment intersection between the segment pair.

        :param other: the second segment
        :param tol: the distance tolerance for collapsing a intersection segment into a point
        :return: the optional point or segment intersection between the segment pair.

        Examples:
          >>> s2 = Segment(Point(0,0), Point(1,0))
          >>> s1 = Segment(Point(0,0), Point(1,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0. 0.]), end_pt=Point([1. 0.]))
          >>> s1 = Segment(Point(-2,0), Point(-1,0))
          >>> s1.intersect_segments2d(s2) is None
          True
          >>> s1 = Segment(Point(-2,0), Point(0,0))
          >>> s1.intersect_segments2d(s2)
          Point([0. 0.])
          >>> s1 = Segment(Point(-2,0), Point(0.5,0.0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0. 0.]), end_pt=Point([0.5 0. ]))
          >>> s1 = Segment(Point(-2,0), Point(1,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0. 0.]), end_pt=Point([1. 0.]))
          >>> s1 = Segment(Point(-2,0), Point(2,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0. 0.]), end_pt=Point([1. 0.]))
          >>> s1 = Segment(Point(0,0), Point(0.5,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0. 0.]), end_pt=Point([0.5 0. ]))
          >>> s1 = Segment(Point(0.25,0), Point(0.75,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0.25 0.  ]), end_pt=Point([0.75 0.  ]))
          >>> s1 = Segment(Point(0.25,0), Point(1,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0.25 0.  ]), end_pt=Point([1. 0.]))
          >>> s1 = Segment(Point(0.25,0), Point(1.25,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0.25 0.  ]), end_pt=Point([1. 0.]))
          >>> s1 = Segment(Point(0,0), Point(1.25,0))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0. 0.]), end_pt=Point([1. 0.]))
          >>> s1 = Segment(Point(1,0), Point(1.25,0))
          >>> s1.intersect_segments2d(s2)
          Point([1. 0.])
          >>> s2 = Segment(Point(0,0), Point(1,1))
          >>> s1 = Segment(Point(0.25,0.25), Point(0.75,0.75))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0.25 0.25]), end_pt=Point([0.75 0.75]))
          >>> s1 = Segment(Point(0.25,0.25), Point(1.75,1.75))
          >>> s1.intersect_segments2d(s2)
          Segment(start_pt=Point([0.25 0.25]), end_pt=Point([1. 1.]))
          >>> s1 = Segment(Point(0.25,0.25), Point(1.75,0))
          >>> s1.intersect_segments2d(s2)
          Point([0.25 0.25])
          >>> s1 = Segment(Point(0.25,1), Point(0.75,0.75))
          >>> s1.intersect_segments2d(s2)
          Point([0.75 0.75])
          >>> s2 = Segment(Point(-1,-1), Point(1,1))
          >>> s1 = Segment(Point(-1,1), Point(1,-1))
          >>> s1.intersect_segments2d(s2)
          Point([0. 0.])
        """


        try:

            segment1 = self.as_2d()
            other = other.as_2d()

            s1_startpt_inside = segment1.segment_start_in(other)
            s2_startpt_inside = other.segment_start_in(segment1)

            s1_endpt_inside = segment1.segment_end_in(other)
            s2_endpt_inside = other.segment_end_in(segment1)

            elements = [s1_startpt_inside, s2_startpt_inside, s1_endpt_inside, s2_endpt_inside]

            if all(elements):
                return segment1.clone()

            if s1_startpt_inside and s1_endpt_inside:
                return segment1.clone()

            if s2_startpt_inside and s2_endpt_inside:
                return other.clone()

            if s1_startpt_inside and s2_startpt_inside:
                return point_or_segment2d(
                    segment1.start_pt,
                    other.start_pt,
                    tol=tol
                )

            if s1_startpt_inside and s2_endpt_inside:
                return point_or_segment2d(
                    segment1.start_pt,
                    other.end_pt,
                    tol=tol
                )

            if s1_endpt_inside and s2_startpt_inside:
                return point_or_segment2d(
                    other.start_pt,
                    segment1.end_pt,
                    tol=tol
                )

            if s1_endpt_inside and s2_endpt_inside:
                return point_or_segment2d(
                    segment1.end_pt,
                    other.end_pt,
                    tol=tol
                )

            if s1_startpt_inside:
                return segment1.start_pt

            if s1_endpt_inside:
                return segment1.end_pt

            if s2_startpt_inside:
                return other.start_pt

            if s2_endpt_inside:
                return other.end_pt

            shortest_segm_or_pt = shortest_segment_or_point2d(
                segment1,
                other,
                tol=tol
            )

            if not shortest_segm_or_pt:
                return None

            if not isinstance(shortest_segm_or_pt, Point):
                return None

            inters_pt = shortest_segm_or_pt

            if not segment1.contains_pt(inters_pt):
                return None

            if not other.contains_pt(inters_pt):
                return None

            return inters_pt

        except Exception as e:

            print(e)
            return None

    def intersection_2d_pt(self,
                           another: 'Segment'
                           ) -> Optional[Point]:
        """

        :param another:
        :return:
        """

        s_len2d = self.length_2d()
        a_len2d = another.length_2d()

        if s_len2d == 0.0 or a_len2d == 0.0:
            return None

        if self.start_pt.x == self.end_pt.x:  # self segment parallel to y axis
            x0 = self.start_pt.x
            m1, p1 = another.segment_2d_m(), another.segment_2d_p()
            if m1 is None:
                return None
            y0 = m1 * x0 + p1
        elif another.start_pt.x == another.end_pt.x:  # another segment parallel to y axis
            x0 = another.start_pt.x
            m1, p1 = self.segment_2d_m(), self.segment_2d_p()
            if m1 is None:
                return None
            y0 = m1 * x0 + p1
        else:  # no segment parallel to y axis
            m0, p0 = self.segment_2d_m(), self.segment_2d_p()
            m1, p1 = another.segment_2d_m(), another.segment_2d_p()
            if m0 is None or m1 is None:
                return None
            x0 = (p1 - p0) / (m0 - m1)
            y0 = m0 * x0 + p0

        return Point(x0, y0)

    def point_distance_from_start(self,
                                  point: Point
                                  ) -> Union[None, numbers.Real]:
        """
        Calculates the optional distance of the point along the segment.
        A zero value is for a point coinciding with the start point.
        Returns None if the point is not contained in the segment.

        :param point: the point to calculate the optional distance in the segment.
        :return: the optional distance of the point along the segment.
        """

        if not self.contains_pt(point):
            return None

        dist = self.start_pt.distance(point)

        return dist

    def subelement_distance_from_start(self,
                                       subelement: 'UnionPtSegment2D'
                                       ) -> Union[None, numbers.Real, Tuple[Union[None, numbers.Real], Union[None, numbers.Real]]]:
        """
        Calculates subelement distance (s) from segment start.
        Subelement is assumed to be fully contained in the segment (extremes comprised).
        """

        if isinstance(subelement, Point):
            return self.point_distance_from_start(subelement)
        elif isinstance(subelement, Segment):
            distance_start_segment = self.point_distance_from_start(subelement.start_pt)
            distance_end_segment = self.point_distance_from_start(subelement.end_pt)
            return distance_start_segment, distance_end_segment
        else:
            return None

    @classmethod
    def random(cls,
        lower_boundary: float = -MAX_SCALAR_VALUE,
        upper_boundary: float = MAX_SCALAR_VALUE):
        """
        Creates a random segment.

        :return: random segment
        """

        return cls(
            start_pt=Point.random(lower_boundary, upper_boundary),
            end_pt=Point.random(lower_boundary, upper_boundary)
        )

    def same_start(self,
                   another: 'Segment',
                   tol: numbers.Real = 1e-12
                   ) -> Union[None, bool]:
        """
        Check whether the two segments have the same start point.

        :param another: a segment to check for.
        :param tol: tolerance for distance between points.
        :return: whether the two segments have the same start point.

        Examples:
          >>> s1 = Segment(Point(0,0,0), Point(1,0,0))
          >>> s2 = Segment(Point(0,0,0), Point(0,1,0))
          >>> s1.same_start(s2)
          True
          >>> s1 = Segment(Point(0,0), Point(1,0))
          >>> s2 = Segment(Point(0,0), Point(0,1))
          >>> s1.same_start(s2)
          True
        """

        return self.start_pt.is_coincident(
            other=another.start_pt,
            tolerance=tol
        )

    def same_end(self,
                 another: 'Segment',
                 tol: numbers.Real = 1e-12
                 ) -> Union[None, bool]:
        """
        Check whether the two segments have the same end point.

        :param another: a segment to check for.
        :param tol: tolerance for distance between points.
        :return: whether the two segments have the same end point.

        Examples:
          >>> s1 = Segment(Point(0,0,0), Point(1,0,0))
          >>> s2 = Segment(Point(2,0,0), Point(1,0,0))
          >>> s1.same_end(s2)
          True
          >>> s1 = Segment(Point(0,0), Point(1,0))
          >>> s2 = Segment(Point(2,0), Point(1,0))
          >>> s1.same_end(s2)
          True
        """

        return self.end_pt.is_coincident(
            other=another.end_pt,
            tolerance=tol)

    def scale(self,
        scale_factor
    ) -> Union[None, 'Segment']:
        """
        Scale a segment by the given scale_factor.
        Start point does not change.

        :param scale_factor: the scale factor, where 1 is the segment length. May be negative.
        :return: Point at scale factor
        """

        end_pt = self.point_at_factor(scale_factor)

        return Segment(
            self.start_pt,
            end_pt)
    def rotate(self,
        rotation_axis: RotationAxis,
        center_point: 'Point' = None
        ) -> Union[type(None), 'Segment']:
        """
        Rotates a segment.

        :param rotation_axis:
        :param center_point:
        :return: the rotated segment

        Examples:
        >>> seg = Segment(Point(0,0,0), Point(0,0,1))
        >>> rot_ax = RotationAxis(0, 0, 90)
        >>> seg.rotate(rot_ax)
        Segment(start_pt=Point([0. 0. 0.]), end_pt=Point([1.00000000e+00 0.00000000e+00 2.22044605e-16]))
        >>> rot_ax = RotationAxis(0, 0, 180)
        >>> centr_pt = Point(0,0,0.5)
        >>> seg.rotate(rotation_axis=rot_ax, center_point=centr_pt)
        Segment(start_pt=Point([-6.123234e-17  0.000000e+00  1.000000e+00]), end_pt=Point([6.123234e-17 0.000000e+00 0.000000e+00]))
        >>> seg = Segment(Point(0,0,0), Point(1,1,0))
        >>> centr_pt = Point(1,0,0)
        >>> rot_ax = RotationAxis(0, 90, 90)
        >>> seg.rotate(rotation_axis=rot_ax, center_point=centr_pt)
        Segment(start_pt=Point([1.000000e+00 1.000000e+00 6.123234e-17]), end_pt=Point([ 2.00000000e+00  2.22044605e-16 -6.12323400e-17]))
        >>> seg = Segment(Point(1,1,1), Point(0,0,0))
        >>> rot_ax = RotationAxis(135, 0, 180)
        >>> centr_pt = Point(0.5,0.5,0.5)
        >>> seg.rotate(rotation_axis=rot_ax, center_point=centr_pt)
        Segment(start_pt=Point([1.66533454e-16 0.00000000e+00 1.11022302e-16]), end_pt=Point([1. 1. 1.]))
        """

        start_pt, end_pt = self

        rotated_start_pt = start_pt.rotate(
            rotation_axis=rotation_axis,
            center_point=center_point
        )

        if rotated_start_pt is None:
            return None

        rotated_end_pt = end_pt.rotate(
            rotation_axis=rotation_axis,
            center_point=center_point
        )

        if rotated_end_pt is None:
            return None

        return Segment(
            start_pt=rotated_start_pt,
            end_pt=rotated_end_pt
        )

    def shift(self,
              dx: numbers.Real,
              dy: numbers.Real
    ) -> 'Segment':
        """
        Shift a segment by dx and dy
        """

        return Segment(
            self.start_pt.shift(dx, dy),
            self.end_pt.shift(dx, dy)
        )

    def segment_start_in(self,
                         another: 'Segment'
                         ) -> bool:
        """
        Check whether the second segment contains the first segment start point.

        :param another: a segment to check for.
        :type another: Segment.
        :return: whether the second segment contains the first segment start point.
        :rtype: bool.

        Examples:
          >>> s1 = Segment(Point(0,0), Point(1,0))
          >>> s2 = Segment(Point(-0.5,0), Point(0.5,0))
          >>> s1.segment_start_in(s2)
          True
          >>> s1 = Segment(Point(0,0), Point(1,1))
          >>> s1.segment_start_in(s2)
          True
          >>> s1 = Segment(Point(0,1), Point(1,1))
          >>> s1.segment_start_in(s2)
          False
          >>> s1 = Segment(Point(-1,-1), Point(1,1))
          >>> s1.segment_start_in(s2)
          False
          >>> s1 = Segment(Point(0,0,0), Point(1,0,0))
          >>> s2 = Segment(Point(-0.5,0,0), Point(0.5,0,0))
          >>> s1.segment_start_in(s2)
          True
          >>> s1 = Segment(Point(0,0,0), Point(1,1,1))
          >>> s1.segment_start_in(s2)
          True
          >>> s1 = Segment(Point(0,1,0), Point(1,1,1))
          >>> s1.segment_start_in(s2)
          False
          >>> s1 = Segment(Point(-1,-1,-1), Point(1,1,1))
          >>> s1.segment_start_in(s2)
          False
        """

        return another.contains_pt(self.start_pt)

    def segment_end_in(self,
                       another: 'Segment'
                       ) -> bool:
        """
        Check whether the second segment contains the first segment end point.

        :param another: a segment to check for.
        :type another: Segment.
        :return: whether the second segment contains the first segment end point.
        :rtype: bool.

        Examples:
          >>> s1 = Segment(Point(0,0), Point(1,0))
          >>> s2 = Segment(Point(-0.5,0), Point(0.5,0))
          >>> s1.segment_end_in(s2)
          False
          >>> s1 = Segment(Point(0,0), Point(1,1))
          >>> s1.segment_end_in(s2)
          False
          >>> s1 = Segment(Point(0,1), Point(1,1))
          >>> s2 = Segment(Point(1,1), Point(0.5,0))
          >>> s1.segment_end_in(s2)
          True
          >>> s1 = Segment(Point(-1,-1), Point(1,1))
          >>> s2 = Segment(Point(0,2), Point(2,0))
          >>> s1.segment_end_in(s2)
          True
          >>> s1 = Segment(Point(0,0,0), Point(1,0,0))
          >>> s2 = Segment(Point(-0.5,0,0), Point(0.5,0,0))
          >>> s1.segment_end_in(s2)
          False
          >>> s1 = Segment(Point(0,0,0), Point(1,1,1))
          >>> s1.segment_end_in(s2)
          False
          >>> s1 = Segment(Point(0,1,0), Point(1,1,1))
          >>> s2 = Segment(Point(1,1,1), Point(0.5,0,0))
          >>> s1.segment_end_in(s2)
          True
          >>> s1 = Segment(Point(-1,-1,3), Point(1,1,3))
          >>> s2 = Segment(Point(0,2,3), Point(2,0,3))
          >>> s1.segment_end_in(s2)
          True
        """

        return another.contains_pt(self.end_pt)

    def tan_slope_angle(self) -> Union[type(None), numbers.Real]:
        """
        Calculates the delta z - delta s ratio of a segment.

        :return: the ratio.
        """

        len2d = self.length_2d()

        if len2d == 0.0:
            return None

        return self.delta_z() / len2d

    def slope_rad(self) -> Union[type(None), numbers.Real]:
        """
        Calculates the slope in radians of the segment.
        Positive is downward point, negative upward pointing.

        :return: the slope as radians.
        """

        delta_zs = self.tan_slope_angle()

        if delta_zs is None:
            return None
        else:
            return - math.atan(delta_zs)

UnionPtSegment2D = Union[Point, Segment]


def shortest_segment_or_point3d(
    first_segment: Segment,
    second_segment: Segment,
    tol: numbers.Real = PRACTICAL_MIN_DIST
) -> Optional[Union[Segment, Point]]:

    """
    Calculates the optional shortest segment - or the intersection point - between two lines represented by two segments.

    Adapted from:
        http://paulbourke.net/geometry/pointlineplane/

    C code from:
        http://paulbourke.net/geometry/pointlineplane/lineline.c
[
    typedef struct {
    double x,y,z;
    } XYZ;

    /*
    Calculate the line segment PaPb that is the shortest route between
    two lines P1P2 and P3P4. Calculate also the values of mua and mub where
      Pa = P1 + mua (P2 - P1)
      Pb = P3 + mub (P4 - P3)
    Return FALSE if no solution exists.
    */
    int LineLineIntersect(
    XYZ p1,XYZ p2,XYZ p3,XYZ p4,XYZ *pa,XYZ *pb,
    double *mua, double *mub)
    {
    XYZ p13,p43,p21;
    double d1343,d4321,d1321,d4343,d2121;
    double numer,denom;

    p13.x = p1.x - p3.x;
    p13.y = p1.y - p3.y;
    p13.z = p1.z - p3.z;
    p43.x = p4.x - p3.x;
    p43.y = p4.y - p3.y;
    p43.z = p4.z - p3.z;
    if (ABS(p43.x) < EPS && ABS(p43.y) < EPS && ABS(p43.z) < EPS)
      return(FALSE);
    p21.x = p2.x - p1.x;
    p21.y = p2.y - p1.y;
    p21.z = p2.z - p1.z;
    if (ABS(p21.x) < EPS && ABS(p21.y) < EPS && ABS(p21.z) < EPS)
      return(FALSE);

    d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z;
    d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z;
    d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z;
    d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z;
    d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z;

    denom = d2121 * d4343 - d4321 * d4321;
    if (ABS(denom) < EPS)
      return(FALSE);
    numer = d1343 * d4321 - d1321 * d4343;

    *mua = numer / denom;
    *mub = (d1343 + d4321 * (*mua)) / d4343;

    pa->x = p1.x + *mua * p21.x;
    pa->y = p1.y + *mua * p21.y;
    pa->z = p1.z + *mua * p21.z;
    pb->x = p3.x + *mub * p43.x;
    pb->y = p3.y + *mub * p43.y;
    pb->z = p3.z + *mub * p43.z;

    return(TRUE);
    }

    :param first_segment: the first segment
    :param second_segment: the second segment
    :param tol: tolerance value for collapsing a segment into the midpoint.
    :return: the optional shortest segment or an intersection point.
    """

    check_type(second_segment, "Second Cartesian line", Segment)

    p1 = first_segment.start_pt
    p2 = first_segment.end_pt

    p3 = second_segment.start_pt
    p4 = second_segment.end_pt

    p13 = Point(
        p1.x - p3.x,
        p1.y - p3.y,
        p1.z - p3.z
    )

    p43 = Point(
        p4.x - p3.x,
        p4.y - p3.y,
        p4.z - p3.z
    )

    if p43.as_vector().is_close_to_zero():
        return None

    p21 = Point(
        p2.x - p1.x,
        p2.y - p1.y,
        p2.z - p1.z,
    )

    if p21.as_vector().is_close_to_zero():
        return None

    d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z
    d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z
    d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z
    d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z
    d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z

    denom = d2121 * d4343 - d4321 * d4321

    if fabs(denom) < MIN_SCALAR_VALUE:
        return None

    numer = d1343 * d4321 - d1321 * d4343

    mua = numer / denom
    mub = (d1343 + d4321 * mua) / d4343

    pa = Point(
        p1.x + mua * p21.x,
        p1.y + mua * p21.y,
        p1.z + mua * p21.z
    )

    pb = Point(
        p3.x + mub * p43.x,
        p3.y + mub * p43.y,
        p3.z + mub * p43.z
    )

    intersection = point_or_segment3d(
        point1=pa,
        point2=pb,
        tol=tol
    )

    return intersection


def point_or_segment2d(
        point1: Point,
        point2: Point,
        tol: numbers.Real = PRACTICAL_MIN_DIST
) -> Union[Point, Segment]:
    """
    Creates a 2D point or segment based on the points distance.

    :param point1: first input point.
    :param point2: second input point.
    :param tol: distance tolerance between the two points.
    :return: point or segment based on their distance.
    :raise: Exception.
    """

    check_type(point1, "First point", Point)
    check_type(point2, "Second point", Point)

    point1 = point1.as_point2d()
    point2 = point2.as_point2d()

    if point1.distance(point2) <= tol:
        return Point(
            (point1.x + point2.x) / 2,
            (point1.y + point2.y) / 2
        )
    else:
        return Segment(
            start_pt=point1,
            end_pt=point2
        )


def shortest_segment_or_point2d(
    first_segment: Segment,
    second_segment: Segment,
    tol: numbers.Real = PRACTICAL_MIN_DIST
) -> Optional[Union[Segment, Point]]:
    """
    TODO: check correct implementation for 2D case, since it derives from 3D implementation.

    Calculates the optional shortest segment - or intersection point - between two segments.

    Adapted from:
        http://paulbourke.net/geometry/pointlineplane/

    C code from:
        http://paulbourke.net/geometry/pointlineplane/lineline.c
[
    typedef struct {
    double x,y,z;
    } XYZ;

    /*
    Calculate the line segment PaPb that is the shortest route between
    two lines P1P2 and P3P4. Calculate also the values of mua and mub where
      Pa = P1 + mua (P2 - P1)
      Pb = P3 + mub (P4 - P3)
    Return FALSE if no solution exists.
    */
    int LineLineIntersect(
    XYZ p1,XYZ p2,XYZ p3,XYZ p4,XYZ *pa,XYZ *pb,
    double *mua, double *mub)
    {
    XYZ p13,p43,p21;
    double d1343,d4321,d1321,d4343,d2121;
    double numer,denom;

    p13.x = p1.x - p3.x;
    p13.y = p1.y - p3.y;
    p13.z = p1.z - p3.z;
    p43.x = p4.x - p3.x;
    p43.y = p4.y - p3.y;
    p43.z = p4.z - p3.z;
    if (ABS(p43.x) < EPS && ABS(p43.y) < EPS && ABS(p43.z) < EPS)
      return(FALSE);
    p21.x = p2.x - p1.x;
    p21.y = p2.y - p1.y;
    p21.z = p2.z - p1.z;
    if (ABS(p21.x) < EPS && ABS(p21.y) < EPS && ABS(p21.z) < EPS)
      return(FALSE);

    d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z;
    d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z;
    d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z;
    d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z;
    d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z;

    denom = d2121 * d4343 - d4321 * d4321;
    if (ABS(denom) < EPS)
      return(FALSE);
    numer = d1343 * d4321 - d1321 * d4343;

    *mua = numer / denom;
    *mub = (d1343 + d4321 * (*mua)) / d4343;

    pa->x = p1.x + *mua * p21.x;
    pa->y = p1.y + *mua * p21.y;
    pa->z = p1.z + *mua * p21.z;
    pb->x = p3.x + *mub * p43.x;
    pb->y = p3.y + *mub * p43.y;
    pb->z = p3.z + *mub * p43.z;

    return(TRUE);
    }

    :param first_segment: the first segment
    :param second_segment: the second segment
    :param tol: tolerance value for collapsing a segment into the midpoint.
    :return: the optional shortest segment or an intersection point.
    """

    check_type(second_segment, "Second Cartesian line", Segment)

    p1 = first_segment.start_pt
    p2 = first_segment.end_pt

    p3 = second_segment.start_pt
    p4 = second_segment.end_pt

    p13 = Point(
        p1.x - p3.x,
        p1.y - p3.y
    )

    p43 = Point(
        p4.x - p3.x,
        p4.y - p3.y
    )

    if p43.is_coincident(Point(0, 0)):
        return None

    p21 = Point(
        p2.x - p1.x,
        p2.y - p1.y
    )

    if p21.is_coincident(Point(0, 0)):
        return None

    d1343 = p13.x * p43.x + p13.y * p43.y
    d4321 = p43.x * p21.x + p43.y * p21.y
    d1321 = p13.x * p21.x + p13.y * p21.y
    d4343 = p43.x * p43.x + p43.y * p43.y
    d2121 = p21.x * p21.x + p21.y * p21.y

    denom = d2121 * d4343 - d4321 * d4321

    if fabs(denom) < MIN_SCALAR_VALUE:
        return None

    numer = d1343 * d4321 - d1321 * d4343

    mua = numer / denom
    mub = (d1343 + d4321 * mua) / d4343

    pa = Point(
        p1.x + mua * p21.x,
        p1.y + mua * p21.y
    )

    pb = Point(
        p3.x + mub * p43.x,
        p3.y + mub * p43.y
    )

    intersection = point_or_segment2d(
        point1=pa,
        point2=pb,
        tol=tol
    )

    return intersection


def point_or_segment3d(
        point1: Point,
        point2: Point,
        tol: numbers.Real = PRACTICAL_MIN_DIST
) -> Union[Point, Segment]:
    """
    Creates a point or segment based on the points distance.

    :param point1: first input point.
    :type point1: Point.
    :param point2: second input point.
    :type point2: Point.
    :param tol: distance tolerance between the two points.
    :type tol: numbers.Real.
    :return: point or segment based on their distance.
    :rtype: PointOrSegment.
    :raise: Exception.
    """

    check_type(point1, "First point", Point)
    check_type(point2, "Second point", Point)

    if point1.distance(point2) <= tol:
        return Point(
            (point1.x + point2.x) / 2,
            (point1.y + point2.y) / 2,
            (point1.z + point2.z) / 2
        )
    else:
        return Segment(
            start_pt=point1,
            end_pt=point2
        )
