Curve Tools 2: move to release T65825

17 files changed, 5979 insertions, 0 deletions

diff --git a/curve_tools/CurveIntersections.py b/curve_tools/CurveIntersections.py
new file mode 100644
index 00000000..06254701
--- /dev/null
+++ b/curve_tools/CurveIntersections.py
@@ -0,0 +1,830 @@
+import bpy
+from . import Math
+from . import Curves
+from . import Util
+
+from mathutils import Vector
+
+algoPOV = None
+algoDIR = None
+
+
+class BezierSegmentIntersectionPoint:
+    def __init__(self, segment, parameter, intersectionPoint):
+        self.segment = segment
+        self.parameter = parameter
+        self.intersectionPoint = intersectionPoint
+
+
+class BezierSegmentsIntersector:
+    def __init__(self, segment1, segment2, worldMatrix1, worldMatrix2):
+        self.segment1 = segment1
+        self.segment2 = segment2
+        self.worldMatrix1 = worldMatrix1
+        self.worldMatrix2 = worldMatrix2
+
+    def CalcFirstIntersection(self, nrSamples1, nrSamples2):
+        algorithm = bpy.context.scene.curvetools.IntersectCurvesAlgorithm
+
+        if algorithm == '3D':
+            return self.CalcFirstRealIntersection3D(nrSamples1, nrSamples2)
+
+        if algorithm == 'From View':
+            global algoDIR
+            if algoDIR is not None:
+                return self.CalcFirstRealIntersectionFromViewDIR(nrSamples1, nrSamples2)
+
+            global algoPOV
+            if algoPOV is not None:
+                return self.CalcFirstRealIntersectionFromViewPOV(nrSamples1, nrSamples2)
+
+        return None
+
+    def CalcFirstIntersection3D(self, nrSamples1, nrSamples2):
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        limitDistance = bpy.context.scene.curvetools.LimitDistance
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointLineSegments(P0, P1, Q0, Q1, limitDistance)
+                if intersectionPointData is None:
+                    continue
+
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                    intersectionSegment1Parameter,
+                                                                    intersectionPointData[2])
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                    intersectionSegment2Parameter,
+                                                                    intersectionPointData[3])
+
+                return [intersectionPoint1, intersectionPoint2]
+
+        return None
+
+    def CalcFirstRealIntersection3D(self, nrSamples1, nrSamples2):
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        limitDistance = bpy.context.scene.curvetools.LimitDistance
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointLineSegments(P0, P1, Q0, Q1, limitDistance)
+                if intersectionPointData is None:
+                    continue
+
+                # intersection point can't be an existing point
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / (fltNrSamples1))
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                if (Math.IsSamePoint(P0, worldPoint1, limitDistance)) or \
+                   (Math.IsSamePoint(P1, worldPoint1, limitDistance)):
+
+                    intersectionPoint1 = None
+                else:
+                    intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                        intersectionSegment1Parameter,
+                                                                        worldPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / (fltNrSamples2))
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                if (Math.IsSamePoint(Q0, worldPoint2, limitDistance)) or \
+                   (Math.IsSamePoint(Q1, worldPoint2, limitDistance)):
+
+                    intersectionPoint2 = None
+                else:
+                    intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                        intersectionSegment2Parameter,
+                                                                        worldPoint2)
+
+                return [intersectionPoint1, intersectionPoint2]
+
+        return None
+
+    def CalcFirstIntersectionFromViewDIR(self, nrSamples1, nrSamples2):
+        global algoDIR
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointsLineSegmentsDIR(P0, P1, Q0, Q1, algoDIR)
+                if intersectionPointData is None:
+                    continue
+
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                    intersectionSegment1Parameter,
+                                                                    worldPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                    intersectionSegment2Parameter,
+                                                                    worldPoint2)
+
+                return [intersectionPoint1, intersectionPoint2]
+
+        return None
+
+    def CalcFirstRealIntersectionFromViewDIR(self, nrSamples1, nrSamples2):
+        global algoDIR
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        limitDistance = bpy.context.scene.curvetools.LimitDistance
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointsLineSegmentsDIR(P0, P1, Q0, Q1, algoDIR)
+                if intersectionPointData is None:
+                    continue
+
+                # intersection point can't be an existing point
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / (fltNrSamples1))
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                if (Math.IsSamePoint(P0, worldPoint1, limitDistance)) or \
+                   (Math.IsSamePoint(P1, worldPoint1, limitDistance)):
+
+                    intersectionPoint1 = None
+                else:
+                    intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                        intersectionSegment1Parameter,
+                                                                        worldPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / (fltNrSamples2))
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                if (Math.IsSamePoint(Q0, worldPoint2, limitDistance)) or \
+                   (Math.IsSamePoint(Q1, worldPoint2, limitDistance)):
+
+                    intersectionPoint2 = None
+                else:
+                    intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                        intersectionSegment2Parameter,
+                                                                        worldPoint2)
+
+                return [intersectionPoint1, intersectionPoint2]
+
+        return None
+
+    def CalcFirstIntersectionFromViewPOV(self, nrSamples1, nrSamples2):
+        global algoPOV
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointsLineSegmentsPOV(P0, P1, Q0, Q1, algoPOV)
+                if intersectionPointData is None:
+                    continue
+
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                    intersectionSegment1Parameter,
+                                                                    worldPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                    intersectionSegment2Parameter,
+                                                                    worldPoint2)
+
+                return [intersectionPoint1, intersectionPoint2]
+
+        return None
+
+    def CalcFirstRealIntersectionFromViewPOV(self, nrSamples1, nrSamples2):
+        global algoPOV
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        limitDistance = bpy.context.scene.curvetools.LimitDistance
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointsLineSegmentsPOV(P0, P1, Q0, Q1, algoPOV)
+                if intersectionPointData is None:
+                    continue
+
+                # intersection point can't be an existing point
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                if (Math.IsSamePoint(P0, worldPoint1, limitDistance)) or \
+                   (Math.IsSamePoint(P1, worldPoint1, limitDistance)):
+
+                    intersectionPoint1 = None
+                else:
+                    intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                        intersectionSegment1Parameter,
+                                                                        worldPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                if (Math.IsSamePoint(Q0, worldPoint2, limitDistance)) or \
+                   (Math.IsSamePoint(Q1, worldPoint2, limitDistance)):
+
+                    intersectionPoint2 = None
+                else:
+                    intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                        intersectionSegment2Parameter,
+                                                                        worldPoint2)
+
+                return [intersectionPoint1, intersectionPoint2]
+
+        return None
+
+    def CalcIntersections(self, nrSamples1, nrSamples2):
+        algorithm = bpy.context.scene.curvetools.IntersectCurvesAlgorithm
+
+        if algorithm == '3D':
+            return self.CalcIntersections3D(nrSamples1, nrSamples2)
+
+        if algorithm == 'From View':
+            global algoDIR
+            if algoDIR is not None:
+                return self.CalcIntersectionsFromViewDIR(nrSamples1, nrSamples2)
+
+            global algoPOV
+            if algoPOV is not None:
+                return self.CalcIntersectionsFromViewPOV(nrSamples1, nrSamples2)
+
+        return [[], []]
+
+    def CalcIntersections3D(self, nrSamples1, nrSamples2):
+        rvIntersections1 = []
+        rvIntersections2 = []
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        limitDistance = bpy.context.scene.curvetools.LimitDistance
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointLineSegments(P0, P1, Q0, Q1, limitDistance)
+                if intersectionPointData is None:
+                    continue
+
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                    intersectionSegment1Parameter,
+                                                                    worldPoint1)
+                rvIntersections1.append(intersectionPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                    intersectionSegment2Parameter,
+                                                                    worldPoint2)
+                rvIntersections2.append(intersectionPoint2)
+
+        return [rvIntersections1, rvIntersections2]
+
+    def CalcIntersectionsFromViewDIR(self, nrSamples1, nrSamples2):
+        global algoDIR
+
+        rvIntersections1 = []
+        rvIntersections2 = []
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointsLineSegmentsDIR(P0, P1, Q0, Q1, algoDIR)
+                if intersectionPointData is None:
+                    continue
+
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                    intersectionSegment1Parameter,
+                                                                    worldPoint1)
+                rvIntersections1.append(intersectionPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                    intersectionSegment2Parameter,
+                                                                    worldPoint2)
+                rvIntersections2.append(intersectionPoint2)
+
+        return [rvIntersections1, rvIntersections2]
+
+    def CalcIntersectionsFromViewPOV(self, nrSamples1, nrSamples2):
+        global algoPOV
+
+        rvIntersections1 = []
+        rvIntersections2 = []
+
+        fltNrSamples1 = float(nrSamples1)
+        fltNrSamples2 = float(nrSamples2)
+
+        for iSample1 in range(nrSamples1):
+            segPar10 = float(iSample1) / fltNrSamples1
+            segPar11 = float(iSample1 + 1) / fltNrSamples1
+            P0 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar10)
+            P1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=segPar11)
+
+            for iSample2 in range(nrSamples2):
+                segPar20 = float(iSample2) / fltNrSamples2
+                segPar21 = float(iSample2 + 1) / fltNrSamples2
+                Q0 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar20)
+                Q1 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=segPar21)
+
+                intersectionPointData = Math.CalcIntersectionPointsLineSegmentsPOV(P0, P1, Q0, Q1, algoPOV)
+                if intersectionPointData is None:
+                    continue
+
+                intersectionSegment1Parameter = segPar10 + (intersectionPointData[0] / fltNrSamples1)
+                worldPoint1 = self.worldMatrix1 @ self.segment1.CalcPoint(parameter=intersectionSegment1Parameter)
+                intersectionPoint1 = BezierSegmentIntersectionPoint(self.segment1,
+                                                                    intersectionSegment1Parameter,
+                                                                    worldPoint1)
+                rvIntersections1.append(intersectionPoint1)
+
+                intersectionSegment2Parameter = segPar20 + (intersectionPointData[1] / fltNrSamples2)
+                worldPoint2 = self.worldMatrix2 @ self.segment2.CalcPoint(parameter=intersectionSegment2Parameter)
+                intersectionPoint2 = BezierSegmentIntersectionPoint(self.segment2,
+                                                                    intersectionSegment2Parameter,
+                                                                    worldPoint2)
+                rvIntersections2.append(intersectionPoint2)
+
+        return [rvIntersections1, rvIntersections2]
+
+
+class BezierSplineIntersectionPoint:
+    def __init__(self, spline, bezierSegmentIntersectionPoint):
+        self.spline = spline
+        self.bezierSegmentIntersectionPoint = bezierSegmentIntersectionPoint
+
+
+class BezierSplinesIntersector:
+    def __init__(self, spline1, spline2, worldMatrix1, worldMatrix2):
+        self.spline1 = spline1
+        self.spline2 = spline2
+        self.worldMatrix1 = worldMatrix1
+        self.worldMatrix2 = worldMatrix2
+
+    def CalcIntersections(self):
+        rvIntersections1 = []
+        rvIntersections2 = []
+
+        try:
+            nrSamplesPerSegment1 = int(self.spline1.resolution / self.spline1.nrSegments)
+        except:
+            nrSamplesPerSegment1 = 2
+        if nrSamplesPerSegment1 < 2:
+            nrSamplesPerSegment1 = 2
+
+        try:
+            nrSamplesPerSegment2 = int(self.spline2.resolution / self.spline2.nrSegments)
+        except:
+            nrSamplesPerSegment2 = 2
+        if nrSamplesPerSegment2 < 2:
+            nrSamplesPerSegment2 = 2
+
+        for segment1 in self.spline1.segments:
+            for segment2 in self.spline2.segments:
+                segmentsIntersector = BezierSegmentsIntersector(segment1, segment2,
+                                                                self.worldMatrix1, self.worldMatrix2)
+                segmentIntersections = segmentsIntersector.CalcIntersections(nrSamplesPerSegment1, nrSamplesPerSegment2)
+                if segmentIntersections is None:
+                    continue
+
+                segment1Intersections = segmentIntersections[0]
+                for segmentIntersection in segment1Intersections:
+                    splineIntersection = BezierSplineIntersectionPoint(self.spline1, segmentIntersection)
+                    rvIntersections1.append(splineIntersection)
+
+                segment2Intersections = segmentIntersections[1]
+                for segmentIntersection in segment2Intersections:
+                    splineIntersection = BezierSplineIntersectionPoint(self.spline2, segmentIntersection)
+                    rvIntersections2.append(splineIntersection)
+
+        return [rvIntersections1, rvIntersections2]
+
+
+class CurvesIntersector:
+    @staticmethod
+    def FromSelection():
+        selObjects = bpy.context.selected_objects
+        if len(selObjects) != 2:
+            raise Exception("len(selObjects) != 2")  # shouldn't be possible
+
+        blenderActiveCurve = bpy.context.active_object
+        blenderOtherCurve = selObjects[0]
+        if blenderActiveCurve == blenderOtherCurve:
+            blenderOtherCurve = selObjects[1]
+
+        aCurve = Curves.Curve(blenderActiveCurve)
+        oCurve = Curves.Curve(blenderOtherCurve)
+
+        return CurvesIntersector(aCurve, oCurve)
+
+    @staticmethod
+    def ResetGlobals():
+        global algoPOV
+        algoPOV = None
+        global algoDIR
+        algoDIR = None
+
+    @staticmethod
+    def InitGlobals():
+        CurvesIntersector.ResetGlobals()
+        global algoPOV
+        global algoDIR
+
+        algo = bpy.context.scene.curvetools.IntersectCurvesAlgorithm
+        if algo == 'From View':
+            regionView3D = Util.GetFirstRegionView3D()
+            if regionView3D is None:
+                print("### ERROR: regionView3D is None. Stopping.")
+                return
+
+            viewPerspective = regionView3D.view_perspective
+            print("--", "viewPerspective:", viewPerspective)
+
+            if viewPerspective == 'ORTHO':
+                viewMatrix = regionView3D.view_matrix
+                print("--", "viewMatrix:")
+                print(viewMatrix)
+
+                algoDIR = Vector((viewMatrix[2][0], viewMatrix[2][1], viewMatrix[2][2]))
+                print("--", "algoDIR:", algoDIR)
+
+            # ## TODO: doesn't work properly
+            if viewPerspective == 'PERSP':
+                viewMatrix = regionView3D.view_matrix
+                print("--", "viewMatrix:")
+                print(viewMatrix)
+
+                algoPOV = regionView3D.view_location.copy()
+                print("--", "algoPOV:", algoPOV)
+
+                otherPOV = Vector((viewMatrix[0][3], viewMatrix[1][3], viewMatrix[2][3]))
+                print("--", "otherPOV:", otherPOV)
+
+                localPOV = Vector((0, 0, 0))
+                globalPOV = viewMatrix * localPOV
+                print("--", "globalPOV:", globalPOV)
+
+                perspMatrix = regionView3D.perspective_matrix
+                print("--", "perspMatrix:")
+                print(perspMatrix)
+
+                globalPOVPersp = perspMatrix * localPOV
+                print("--", "globalPOVPersp:", globalPOVPersp)
+
+            if viewPerspective == 'CAMERA':
+                camera = bpy.context.scene.camera
+                if camera is None:
+                    print("### ERROR: camera is None. Stopping.")
+                    return
+
+                print("--", "camera:", camera)
+                cameraData = camera.data
+                print("--", "cameraData.type:", cameraData.type)
+
+                cameraMatrix = camera.matrix_world
+                print("--", "cameraMatrix:")
+                print(cameraMatrix)
+
+                if cameraData.type == 'ORTHO':
+                    cameraMatrix = camera.matrix_world
+                    # algoDIR = Vector((cameraMatrix[2][0], cameraMatrix[2][1], cameraMatrix[2][2]))
+                    algoDIR = Vector((- cameraMatrix[0][2], - cameraMatrix[1][2], - cameraMatrix[2][2]))
+                    print("--", "algoDIR:", algoDIR)
+
+                if cameraData.type == 'PERSP':
+                    algoPOV = camera.location.copy()
+                    print("--", "algoPOV:", algoPOV)
+
+    def __init__(self, activeCurve, otherCurve):
+        self.activeCurve = activeCurve
+        self.otherCurve = otherCurve
+
+        CurvesIntersector.InitGlobals()
+
+    def CalcIntersections(self):
+        rvIntersections1 = []
+        rvIntersections2 = []
+
+        worldMatrix1 = self.activeCurve.curve.matrix_world
+        worldMatrix2 = self.otherCurve.curve.matrix_world
+
+        for spline1 in self.activeCurve.splines:
+            for spline2 in self.otherCurve.splines:
+                splineIntersector = BezierSplinesIntersector(spline1, spline2, worldMatrix1, worldMatrix2)
+                splineIntersections = splineIntersector.CalcIntersections()
+                if splineIntersections is None:
+                    continue
+
+                spline1Intersections = splineIntersections[0]
+                for splineIntersection in spline1Intersections:
+                    rvIntersections1.append(splineIntersection)
+
+                spline2Intersections = splineIntersections[1]
+                for splineIntersection in spline2Intersections:
+                    rvIntersections2.append(splineIntersection)
+
+        return [rvIntersections1, rvIntersections2]
+
+    def CalcAndApplyIntersections(self):
+        mode = bpy.context.scene.curvetools.IntersectCurvesMode
+
+        if mode == 'Empty':
+            return self.CalcAndApplyEmptyAtIntersections()
+        if mode == 'Insert':
+            return self.CalcAndApplyInsertAtIntersections()
+        if mode == 'Split':
+            return self.CalcAndApplySplitAtIntersections()
+
+        return [0, 0]
+
+    def CalcAndApplyEmptyAtIntersections(self):
+        intersections = self.CalcIntersections()
+        intersectionsActive = intersections[0]
+        intersectionsOther = intersections[1]
+
+        nrActive = 0
+        nrOther = 0
+
+        affect = bpy.context.scene.curvetools.IntersectCurvesAffect
+
+        if (affect == 'Both') or (affect == 'Active'):
+            for splineIntersection in intersectionsActive:
+                iPoint = splineIntersection.bezierSegmentIntersectionPoint.intersectionPoint
+                bpy.ops.object.empty_add(type='PLAIN_AXES',
+                                         align='WORLD',
+                                         location=(iPoint.x, iPoint.y, iPoint.z), rotation=(0, 0, 0))
+                nrActive += 1
+
+        if (affect == 'Both') or (affect == 'Other'):
+            for splineIntersection in intersectionsOther:
+                iPoint = splineIntersection.bezierSegmentIntersectionPoint.intersectionPoint
+                bpy.ops.object.empty_add(type='PLAIN_AXES',
+                                         align='WORLD',
+                                         location=(iPoint.x, iPoint.y, iPoint.z), rotation=(0, 0, 0))
+                nrOther += 1
+
+        return [nrActive, nrOther]
+
+    def CalcAndApplyInsertAtIntersections(self):
+        nrActive = 0
+        nrOther = 0
+
+        affect = bpy.context.scene.curvetools.IntersectCurvesAffect
+        affectA = (affect == 'Both') or (affect == 'Active')
+        affectO = (affect == 'Both') or (affect == 'Other')
+
+        for iSplineA in range(len(self.activeCurve.splines)):
+            splineA = self.activeCurve.splines[iSplineA]
+            nrSegmentsA = len(splineA.segments)
+            resPerSegA = splineA.resolutionPerSegment
+
+            for iSplineO in range(len(self.otherCurve.splines)):
+                splineO = self.otherCurve.splines[iSplineO]
+                nrSegmentsO = len(splineO.segments)
+                resPerSegO = splineO.resolutionPerSegment
+
+                iSegA = 0
+                while True:
+                    segA = splineA.segments[iSegA]
+
+                    iSegO = 0
+                    while True:
+                        segO = splineO.segments[iSegO]
+
+                        segIntersector = BezierSegmentsIntersector(segA, segO,
+                                                                   self.activeCurve.worldMatrix,
+                                                                   self.otherCurve.worldMatrix)
+                        segFirstIntersection = segIntersector.CalcFirstIntersection(resPerSegA, resPerSegO)
+
+                        if segFirstIntersection is not None:
+                            intPointA = segFirstIntersection[0]
+                            intPointO = segFirstIntersection[1]
+                            # else does something weird if 1 of them is None..
+                            if (intPointA is not None) and (intPointO is not None):
+                                if affectA:
+                                    if intPointA is not None:
+                                        splineA.InsertPoint(segA, intPointA.parameter)
+
+                                        nrActive += 1
+                                        nrSegmentsA += 1
+
+                                if affectO:
+                                    if intPointO is not None:
+                                        splineO.InsertPoint(segO, intPointO.parameter)
+
+                                        nrOther += 1
+                                        nrSegmentsO += 1
+
+                        iSegO += 1
+                        if not (iSegO < nrSegmentsO):
+                            break
+
+                    iSegA += 1
+                    if not (iSegA < nrSegmentsA):
+                        break
+
+                if affectO:
+                    splineO.RefreshInScene()
+
+            if affectA:
+                splineA.RefreshInScene()
+
+        return [nrActive, nrOther]
+
+    def CalcAndApplySplitAtIntersections(self):
+        nrActive = 0
+        nrOther = 0
+
+        affect = bpy.context.scene.curvetools.IntersectCurvesAffect
+        affectA = (affect == 'Both') or (affect == 'Active')
+        affectO = (affect == 'Both') or (affect == 'Other')
+
+        nrSplinesA = len(self.activeCurve.splines)
+        nrSplinesO = len(self.otherCurve.splines)
+
+        iSplineA = 0
+        while True:
+            splineA = self.activeCurve.splines[iSplineA]
+            nrSegmentsA = len(splineA.segments)
+            resPerSegA = splineA.resolutionPerSegment
+
+            iSplineO = 0
+            while True:
+                splineO = self.otherCurve.splines[iSplineO]
+                nrSegmentsO = len(splineO.segments)
+                resPerSegO = splineO.resolutionPerSegment
+
+                iSegA = 0
+                while True:
+                    segA = splineA.segments[iSegA]
+
+                    iSegO = 0
+                    while True:
+                        segO = splineO.segments[iSegO]
+
+                        segIntersector = BezierSegmentsIntersector(segA, segO,
+                                                                   self.activeCurve.worldMatrix,
+                                                                   self.otherCurve.worldMatrix)
+                        segFirstIntersection = segIntersector.CalcFirstIntersection(resPerSegA, resPerSegO)
+
+                        if segFirstIntersection is not None:
+                            intPointA = segFirstIntersection[0]
+                            intPointO = segFirstIntersection[1]
+                            # else does something weird if 1 of them is None..
+                            if (intPointA is not None) and (intPointO is not None):
+                                if affectA:
+                                    if intPointA is not None:
+                                        print("--", "splineA.Split():")
+                                        newSplinesA = splineA.Split(segA, intPointA.parameter)
+                                        if newSplinesA is not None:
+                                            newResolutions = splineA.CalcDivideResolution(segA, intPointA.parameter)
+                                            newSplinesA[0].resolution = newResolutions[0]
+                                            newSplinesA[1].resolution = newResolutions[1]
+
+                                            splineA = newSplinesA[0]
+                                            self.activeCurve.splines[iSplineA] = splineA
+                                            self.activeCurve.splines.insert(iSplineA + 1, newSplinesA[1])
+
+                                            nrActive += 1
+
+                                if affectO:
+                                    if intPointO is not None:
+                                        print("--", "splineO.Split():")
+                                        newSplinesO = splineO.Split(segO, intPointO.parameter)
+                                        if newSplinesO is not None:
+                                            newResolutions = splineO.CalcDivideResolution(segO, intPointO.parameter)
+                                            newSplinesO[0].resolution = newResolutions[0]
+                                            newSplinesO[1].resolution = newResolutions[1]
+
+                                            splineO = newSplinesO[0]
+                                            self.otherCurve.splines[iSplineO] = splineO
+                                            self.otherCurve.splines.insert(iSplineO + 1, newSplinesO[1])
+
+                                            nrOther += 1
+
+                        nrSegmentsO = len(splineO.segments)
+                        iSegO += 1
+                        if not (iSegO < nrSegmentsO):
+                            break
+
+                    nrSegmentsA = len(splineA.segments)
+                    iSegA += 1
+                    if not (iSegA < nrSegmentsA):
+                        break
+
+                nrSplinesO = len(self.otherCurve.splines)
+                iSplineO += 1
+                if not (iSplineO < nrSplinesO):
+                    break
+
+            nrSplinesA = len(self.activeCurve.splines)
+            iSplineA += 1
+            if not (iSplineA < nrSplinesA):
+                break
+
+        if affectA:
+            print("")
+            print("--", "self.activeCurve.RebuildInScene():")
+            self.activeCurve.RebuildInScene()
+        if affectO:
+            print("")
+            print("--", "self.otherCurve.RebuildInScene():")
+            self.otherCurve.RebuildInScene()
+
+        return [nrActive, nrOther]
diff --git a/curve_tools/Curves.py b/curve_tools/Curves.py
new file mode 100644
index 00000000..e2608eeb
--- /dev/null
+++ b/curve_tools/Curves.py
@@ -0,0 +1,594 @@
+from . import Math
+
+import bpy
+
+
+class BezierPoint:
+    @staticmethod
+    def FromBlenderBezierPoint(blenderBezierPoint):
+        return BezierPoint(blenderBezierPoint.handle_left, blenderBezierPoint.co, blenderBezierPoint.handle_right)
+
+
+    def __init__(self, handle_left, co, handle_right):
+        self.handle_left = handle_left
+        self.co = co
+        self.handle_right = handle_right
+
+
+    def Copy(self):
+        return BezierPoint(self.handle_left.copy(), self.co.copy(), self.handle_right.copy())
+
+    def Reversed(self):
+        return BezierPoint(self.handle_right, self.co, self.handle_left)
+
+    def Reverse(self):
+        tmp = self.handle_left
+        self.handle_left = self.handle_right
+        self.handle_right = tmp
+
+
+class BezierSegment:
+    @staticmethod
+    def FromBlenderBezierPoints(blenderBezierPoint1, blenderBezierPoint2):
+        bp1 = BezierPoint.FromBlenderBezierPoint(blenderBezierPoint1)
+        bp2 = BezierPoint.FromBlenderBezierPoint(blenderBezierPoint2)
+
+        return BezierSegment(bp1, bp2)
+
+
+    def Copy(self):
+        return BezierSegment(self.bezierPoint1.Copy(), self.bezierPoint2.Copy())
+
+    def Reversed(self):
+        return BezierSegment(self.bezierPoint2.Reversed(), self.bezierPoint1.Reversed())
+
+    def Reverse(self):
+        # make a copy, otherwise neighboring segment may be affected
+        tmp = self.bezierPoint1.Copy()
+        self.bezierPoint1 = self.bezierPoint2.Copy()
+        self.bezierPoint2 = tmp
+        self.bezierPoint1.Reverse()
+        self.bezierPoint2.Reverse()
+
+
+    def __init__(self, bezierPoint1, bezierPoint2):
+        # bpy.types.BezierSplinePoint
+        # ## NOTE/TIP: copy() helps with repeated (intersection) action -- ??
+        self.bezierPoint1 = bezierPoint1.Copy()
+        self.bezierPoint2 = bezierPoint2.Copy()
+
+        self.ctrlPnt0 = self.bezierPoint1.co
+        self.ctrlPnt1 = self.bezierPoint1.handle_right
+        self.ctrlPnt2 = self.bezierPoint2.handle_left
+        self.ctrlPnt3 = self.bezierPoint2.co
+
+        self.coeff0 = self.ctrlPnt0
+        self.coeff1 = self.ctrlPnt0 * (-3.0) + self.ctrlPnt1 * (+3.0)
+        self.coeff2 = self.ctrlPnt0 * (+3.0) + self.ctrlPnt1 * (-6.0) + self.ctrlPnt2 * (+3.0)
+        self.coeff3 = self.ctrlPnt0 * (-1.0) + self.ctrlPnt1 * (+3.0) + self.ctrlPnt2 * (-3.0) + self.ctrlPnt3
+
+
+    def CalcPoint(self, parameter = 0.5):
+        parameter2 = parameter * parameter
+        parameter3 = parameter * parameter2
+
+        rvPoint = self.coeff0 + self.coeff1 * parameter + self.coeff2 * parameter2 + self.coeff3 * parameter3
+
+        return rvPoint
+
+
+    def CalcDerivative(self, parameter = 0.5):
+        parameter2 = parameter * parameter
+
+        rvPoint = self.coeff1 + self.coeff2 * parameter * 2.0 + self.coeff3 * parameter2 * 3.0
+
+        return rvPoint
+
+
+    def CalcLength(self, nrSamples = 2):
+        nrSamplesFloat = float(nrSamples)
+        rvLength = 0.0
+        for iSample in range(nrSamples):
+            par1 = float(iSample) / nrSamplesFloat
+            par2 = float(iSample + 1) / nrSamplesFloat
+
+            point1 = self.CalcPoint(parameter = par1)
+            point2 = self.CalcPoint(parameter = par2)
+            diff12 = point1 - point2
+
+            rvLength += diff12.magnitude
+
+        return rvLength
+
+
+    #http://en.wikipedia.org/wiki/De_Casteljau's_algorithm
+    def CalcSplitPoint(self, parameter = 0.5):
+        par1min = 1.0 - parameter
+
+        bez00 = self.ctrlPnt0
+        bez01 = self.ctrlPnt1
+        bez02 = self.ctrlPnt2
+        bez03 = self.ctrlPnt3
+
+        bez10 = bez00 * par1min + bez01 * parameter
+        bez11 = bez01 * par1min + bez02 * parameter
+        bez12 = bez02 * par1min + bez03 * parameter
+
+        bez20 = bez10 * par1min + bez11 * parameter
+        bez21 = bez11 * par1min + bez12 * parameter
+
+        bez30 = bez20 * par1min + bez21 * parameter
+
+        bezPoint1 = BezierPoint(self.bezierPoint1.handle_left, bez00, bez10)
+        bezPointNew = BezierPoint(bez20, bez30, bez21)
+        bezPoint2 = BezierPoint(bez12, bez03, self.bezierPoint2.handle_right)
+
+        return [bezPoint1, bezPointNew, bezPoint2]
+
+
+class BezierSpline:
+    @staticmethod
+    def FromSegments(listSegments):
+        rvSpline = BezierSpline(None)
+
+        rvSpline.segments = listSegments
+
+        return rvSpline
+
+
+    def __init__(self, blenderBezierSpline):
+        if not blenderBezierSpline is None:
+            if blenderBezierSpline.type != 'BEZIER':
+                print("## ERROR:", "blenderBezierSpline.type != 'BEZIER'")
+                raise Exception("blenderBezierSpline.type != 'BEZIER'")
+            if len(blenderBezierSpline.bezier_points) < 1:
+                if not blenderBezierSpline.use_cyclic_u:
+                    print("## ERROR:", "len(blenderBezierSpline.bezier_points) < 1")
+                    raise Exception("len(blenderBezierSpline.bezier_points) < 1")
+
+        self.bezierSpline = blenderBezierSpline
+
+        self.resolution = 12
+        self.isCyclic = False
+        if not self.bezierSpline is None:
+            self.resolution = self.bezierSpline.resolution_u
+            self.isCyclic = self.bezierSpline.use_cyclic_u
+
+        self.segments = self.SetupSegments()
+
+
+    def __getattr__(self, attrName):
+        if attrName == "nrSegments":
+            return len(self.segments)
+
+        if attrName == "bezierPoints":
+            rvList = []
+
+            for seg in self.segments: rvList.append(seg.bezierPoint1)
+            if not self.isCyclic: rvList.append(self.segments[-1].bezierPoint2)
+
+            return rvList
+
+        if attrName == "resolutionPerSegment":
+            try: rvResPS = int(self.resolution / self.nrSegments)
+            except: rvResPS = 2
+            if rvResPS < 2: rvResPS = 2
+
+            return rvResPS
+
+        if attrName == "length":
+            return self.CalcLength()
+
+        return None
+
+
+    def SetupSegments(self):
+        rvSegments = []
+        if self.bezierSpline is None: return rvSegments
+
+        nrBezierPoints = len(self.bezierSpline.bezier_points)
+        for iBezierPoint in range(nrBezierPoints - 1):
+            bezierPoint1 = self.bezierSpline.bezier_points[iBezierPoint]
+            bezierPoint2 = self.bezierSpline.bezier_points[iBezierPoint + 1]
+            rvSegments.append(BezierSegment.FromBlenderBezierPoints(bezierPoint1, bezierPoint2))
+        if self.isCyclic:
+            bezierPoint1 = self.bezierSpline.bezier_points[-1]
+            bezierPoint2 = self.bezierSpline.bezier_points[0]
+            rvSegments.append(BezierSegment.FromBlenderBezierPoints(bezierPoint1, bezierPoint2))
+
+        return rvSegments
+
+
+    def UpdateSegments(self, newSegments):
+        prevNrSegments = len(self.segments)
+        diffNrSegments = len(newSegments) - prevNrSegments
+        if diffNrSegments > 0:
+            newBezierPoints = []
+            for segment in newSegments: newBezierPoints.append(segment.bezierPoint1)
+            if not self.isCyclic: newBezierPoints.append(newSegments[-1].bezierPoint2)
+
+            self.bezierSpline.bezier_points.add(diffNrSegments)
+
+            for i, bezPoint in enumerate(newBezierPoints):
+                blBezPoint = self.bezierSpline.bezier_points[i]
+
+                blBezPoint.tilt = 0
+                blBezPoint.radius = 1.0
+
+                blBezPoint.handle_left_type = 'FREE'
+                blBezPoint.handle_left = bezPoint.handle_left
+                blBezPoint.co = bezPoint.co
+                blBezPoint.handle_right_type = 'FREE'
+                blBezPoint.handle_right = bezPoint.handle_right
+
+            self.segments = newSegments
+        else:
+            print("### WARNING: UpdateSegments(): not diffNrSegments > 0")
+
+
+    def Reversed(self):
+        revSegments = []
+
+        for iSeg in reversed(range(self.nrSegments)): revSegments.append(self.segments[iSeg].Reversed())
+
+        rvSpline = BezierSpline.FromSegments(revSegments)
+        rvSpline.resolution = self.resolution
+        rvSpline.isCyclic = self.isCyclic
+
+        return rvSpline
+
+
+    def Reverse(self):
+        revSegments = []
+
+        for iSeg in reversed(range(self.nrSegments)):
+            self.segments[iSeg].Reverse()
+            revSegments.append(self.segments[iSeg])
+
+        self.segments = revSegments
+
+
+    def CalcDivideResolution(self, segment, parameter):
+        if not segment in self.segments:
+            print("### WARNING: InsertPoint(): not segment in self.segments")
+            return None
+
+        iSeg = self.segments.index(segment)
+        dPar = 1.0 / self.nrSegments
+        splinePar = dPar * (parameter + float(iSeg))
+
+        res1 = int(splinePar * self.resolution)
+        if res1 < 2:
+            print("### WARNING: CalcDivideResolution(): res1 < 2 -- res1: %d" % res1, "-- setting it to 2")
+            res1 = 2
+
+        res2 = int((1.0 - splinePar) * self.resolution)
+        if res2 < 2:
+            print("### WARNING: CalcDivideResolution(): res2 < 2 -- res2: %d" % res2, "-- setting it to 2")
+            res2 = 2
+
+        return [res1, res2]
+        # return [self.resolution, self.resolution]
+
+
+    def CalcPoint(self, parameter):
+        nrSegs = self.nrSegments
+
+        segmentIndex = int(nrSegs * parameter)
+        if segmentIndex < 0: segmentIndex = 0
+        if segmentIndex > (nrSegs - 1): segmentIndex = nrSegs - 1
+
+        segmentParameter = nrSegs * parameter - segmentIndex
+        if segmentParameter < 0.0: segmentParameter = 0.0
+        if segmentParameter > 1.0: segmentParameter = 1.0
+
+        return self.segments[segmentIndex].CalcPoint(parameter = segmentParameter)
+
+
+    def CalcDerivative(self, parameter):
+        nrSegs = self.nrSegments
+
+        segmentIndex = int(nrSegs * parameter)
+        if segmentIndex < 0: segmentIndex = 0
+        if segmentIndex > (nrSegs - 1): segmentIndex = nrSegs - 1
+
+        segmentParameter = nrSegs * parameter - segmentIndex
+        if segmentParameter < 0.0: segmentParameter = 0.0
+        if segmentParameter > 1.0: segmentParameter = 1.0
+
+        return self.segments[segmentIndex].CalcDerivative(parameter = segmentParameter)
+
+
+    def InsertPoint(self, segment, parameter):
+        if not segment in self.segments:
+            print("### WARNING: InsertPoint(): not segment in self.segments")
+            return
+        iSeg = self.segments.index(segment)
+        nrSegments = len(self.segments)
+
+        splitPoints = segment.CalcSplitPoint(parameter = parameter)
+        bezPoint1 = splitPoints[0]
+        bezPointNew = splitPoints[1]
+        bezPoint2 = splitPoints[2]
+
+        segment.bezierPoint1.handle_right = bezPoint1.handle_right
+        segment.bezierPoint2 = bezPointNew
+
+        if iSeg < (nrSegments - 1):
+            nextSeg = self.segments[iSeg + 1]
+            nextSeg.bezierPoint1.handle_left = bezPoint2.handle_left
+        else:
+            if self.isCyclic:
+                nextSeg = self.segments[0]
+                nextSeg.bezierPoint1.handle_left = bezPoint2.handle_left
+
+
+        newSeg = BezierSegment(bezPointNew, bezPoint2)
+        self.segments.insert(iSeg + 1, newSeg)
+
+
+    def Split(self, segment, parameter):
+        if not segment in self.segments:
+            print("### WARNING: InsertPoint(): not segment in self.segments")
+            return None
+        iSeg = self.segments.index(segment)
+        nrSegments = len(self.segments)
+
+        splitPoints = segment.CalcSplitPoint(parameter = parameter)
+        bezPoint1 = splitPoints[0]
+        bezPointNew = splitPoints[1]
+        bezPoint2 = splitPoints[2]
+
+
+        newSpline1Segments = []
+        for iSeg1 in range(iSeg): newSpline1Segments.append(self.segments[iSeg1])
+        if len(newSpline1Segments) > 0: newSpline1Segments[-1].bezierPoint2.handle_right = bezPoint1.handle_right
+        newSpline1Segments.append(BezierSegment(bezPoint1, bezPointNew))
+
+        newSpline2Segments = []
+        newSpline2Segments.append(BezierSegment(bezPointNew, bezPoint2))
+        for iSeg2 in range(iSeg + 1, nrSegments): newSpline2Segments.append(self.segments[iSeg2])
+        if len(newSpline2Segments) > 1: newSpline2Segments[1].bezierPoint1.handle_left = newSpline2Segments[0].bezierPoint2.handle_left
+
+
+        newSpline1 = BezierSpline.FromSegments(newSpline1Segments)
+        newSpline2 = BezierSpline.FromSegments(newSpline2Segments)
+
+        return [newSpline1, newSpline2]
+
+
+    def Join(self, spline2, mode = 'At midpoint'):
+        if mode == 'At midpoint':
+            self.JoinAtMidpoint(spline2)
+            return
+
+        if mode == 'Insert segment':
+            self.JoinInsertSegment(spline2)
+            return
+
+        print("### ERROR: Join(): unknown mode:", mode)
+
+
+    def JoinAtMidpoint(self, spline2):
+        bezPoint1 = self.segments[-1].bezierPoint2
+        bezPoint2 = spline2.segments[0].bezierPoint1
+
+        mpHandleLeft = bezPoint1.handle_left.copy()
+        mpCo = (bezPoint1.co + bezPoint2.co) * 0.5
+        mpHandleRight = bezPoint2.handle_right.copy()
+        mpBezPoint = BezierPoint(mpHandleLeft, mpCo, mpHandleRight)
+
+        self.segments[-1].bezierPoint2 = mpBezPoint
+        spline2.segments[0].bezierPoint1 = mpBezPoint
+        for seg2 in spline2.segments: self.segments.append(seg2)
+
+        self.resolution += spline2.resolution
+        self.isCyclic = False    # is this ok?
+
+
+    def JoinInsertSegment(self, spline2):
+        self.segments.append(BezierSegment(self.segments[-1].bezierPoint2, spline2.segments[0].bezierPoint1))
+        for seg2 in spline2.segments: self.segments.append(seg2)
+
+        self.resolution += spline2.resolution    # extra segment will usually be short -- impact on resolution negligable
+
+        self.isCyclic = False    # is this ok?
+
+
+    def RefreshInScene(self):
+        bezierPoints = self.bezierPoints
+
+        currNrBezierPoints = len(self.bezierSpline.bezier_points)
+        diffNrBezierPoints = len(bezierPoints) - currNrBezierPoints
+        if diffNrBezierPoints > 0: self.bezierSpline.bezier_points.add(diffNrBezierPoints)
+
+        for i, bezPoint in enumerate(bezierPoints):
+            blBezPoint = self.bezierSpline.bezier_points[i]
+
+            blBezPoint.tilt = 0
+            blBezPoint.radius = 1.0
+
+            blBezPoint.handle_left_type = 'FREE'
+            blBezPoint.handle_left = bezPoint.handle_left
+            blBezPoint.co = bezPoint.co
+            blBezPoint.handle_right_type = 'FREE'
+            blBezPoint.handle_right = bezPoint.handle_right
+
+        self.bezierSpline.use_cyclic_u = self.isCyclic
+        self.bezierSpline.resolution_u = self.resolution
+
+
+    def CalcLength(self):
+        try: nrSamplesPerSegment = int(self.resolution / self.nrSegments)
+        except: nrSamplesPerSegment = 2
+        if nrSamplesPerSegment < 2: nrSamplesPerSegment = 2
+
+        rvLength = 0.0
+        for segment in self.segments:
+            rvLength += segment.CalcLength(nrSamples = nrSamplesPerSegment)
+
+        return rvLength
+
+
+    def GetLengthIsSmallerThan(self, threshold):
+        try: nrSamplesPerSegment = int(self.resolution / self.nrSegments)
+        except: nrSamplesPerSegment = 2
+        if nrSamplesPerSegment < 2: nrSamplesPerSegment = 2
+
+        length = 0.0
+        for segment in self.segments:
+            length += segment.CalcLength(nrSamples = nrSamplesPerSegment)
+            if not length < threshold: return False
+
+        return True
+
+
+class Curve:
+    def __init__(self, blenderCurve):
+        self.curve = blenderCurve
+        self.curveData = blenderCurve.data
+
+        self.splines = self.SetupSplines()
+
+
+    def __getattr__(self, attrName):
+        if attrName == "nrSplines":
+            return len(self.splines)
+
+        if attrName == "length":
+            return self.CalcLength()
+
+        if attrName == "worldMatrix":
+            return self.curve.matrix_world
+
+        if attrName == "location":
+            return self.curve.location
+
+        return None
+
+
+    def SetupSplines(self):
+        rvSplines = []
+        for spline in self.curveData.splines:
+            if spline.type != 'BEZIER':
+                print("## WARNING: only bezier splines are supported, atm; other types are ignored")
+                continue
+
+            try: newSpline = BezierSpline(spline)
+            except:
+                print("## EXCEPTION: newSpline = BezierSpline(spline)")
+                continue
+
+            rvSplines.append(newSpline)
+
+        return rvSplines
+
+
+    def RebuildInScene(self):
+        self.curveData.splines.clear()
+
+        for spline in self.splines:
+            blSpline = self.curveData.splines.new('BEZIER')
+            blSpline.use_cyclic_u = spline.isCyclic
+            blSpline.resolution_u = spline.resolution
+
+            bezierPoints = []
+            for segment in spline.segments: bezierPoints.append(segment.bezierPoint1)
+            if not spline.isCyclic: bezierPoints.append(spline.segments[-1].bezierPoint2)
+            #else: print("????", "spline.isCyclic")
+
+            nrBezierPoints = len(bezierPoints)
+            blSpline.bezier_points.add(nrBezierPoints - 1)
+
+            for i, blBezPoint in enumerate(blSpline.bezier_points):
+                bezPoint = bezierPoints[i]
+
+                blBezPoint.tilt = 0
+                blBezPoint.radius = 1.0
+
+                blBezPoint.handle_left_type = 'FREE'
+                blBezPoint.handle_left = bezPoint.handle_left
+                blBezPoint.co = bezPoint.co
+                blBezPoint.handle_right_type = 'FREE'
+                blBezPoint.handle_right = bezPoint.handle_right
+
+
+    def CalcLength(self):
+        rvLength = 0.0
+        for spline in self.splines:
+            rvLength += spline.length
+
+        return rvLength
+
+
+    def RemoveShortSplines(self, threshold):
+        splinesToRemove = []
+
+        for spline in self.splines:
+            if spline.GetLengthIsSmallerThan(threshold): splinesToRemove.append(spline)
+
+        for spline in splinesToRemove: self.splines.remove(spline)
+
+        return len(splinesToRemove)
+
+
+    def JoinNeighbouringSplines(self, startEnd, threshold, mode):
+        nrJoins = 0
+
+        while True:
+            firstPair = self.JoinGetFirstPair(startEnd, threshold)
+            if firstPair is None: break
+
+            firstPair[0].Join(firstPair[1], mode)
+            self.splines.remove(firstPair[1])
+
+            nrJoins += 1
+
+        return nrJoins
+
+
+    def JoinGetFirstPair(self, startEnd, threshold):
+        nrSplines = len(self.splines)
+
+        if startEnd:
+            for iCurrentSpline in range(nrSplines):
+                currentSpline = self.splines[iCurrentSpline]
+
+                for iNextSpline in range(iCurrentSpline + 1, nrSplines):
+                    nextSpline = self.splines[iNextSpline]
+
+                    currEndPoint = currentSpline.segments[-1].bezierPoint2.co
+                    nextStartPoint = nextSpline.segments[0].bezierPoint1.co
+                    if Math.IsSamePoint(currEndPoint, nextStartPoint, threshold): return [currentSpline, nextSpline]
+
+                    nextEndPoint = nextSpline.segments[-1].bezierPoint2.co
+                    currStartPoint = currentSpline.segments[0].bezierPoint1.co
+                    if Math.IsSamePoint(nextEndPoint, currStartPoint, threshold): return [nextSpline, currentSpline]
+
+            return None
+        else:
+            for iCurrentSpline in range(nrSplines):
+                currentSpline = self.splines[iCurrentSpline]
+
+                for iNextSpline in range(iCurrentSpline + 1, nrSplines):
+                    nextSpline = self.splines[iNextSpline]
+
+                    currEndPoint = currentSpline.segments[-1].bezierPoint2.co
+                    nextStartPoint = nextSpline.segments[0].bezierPoint1.co
+                    if Math.IsSamePoint(currEndPoint, nextStartPoint, threshold): return [currentSpline, nextSpline]
+
+                    nextEndPoint = nextSpline.segments[-1].bezierPoint2.co
+                    currStartPoint = currentSpline.segments[0].bezierPoint1.co
+                    if Math.IsSamePoint(nextEndPoint, currStartPoint, threshold): return [nextSpline, currentSpline]
+
+                    if Math.IsSamePoint(currEndPoint, nextEndPoint, threshold):
+                        nextSpline.Reverse()
+                        #print("## ", "nextSpline.Reverse()")
+                        return [currentSpline, nextSpline]
+
+                    if Math.IsSamePoint(currStartPoint, nextStartPoint, threshold):
+                        currentSpline.Reverse()
+                        #print("## ", "currentSpline.Reverse()")
+                        return [currentSpline, nextSpline]
+
+            return None
diff --git a/curve_tools/Math.py b/curve_tools/Math.py
new file mode 100644
index 00000000..4a61af4d
--- /dev/null
+++ b/curve_tools/Math.py
@@ -0,0 +1,145 @@
+from mathutils import *
+
+
+def IsSamePoint(v31, v32, limitDistance):
+    if (v31 - v32).magnitude < limitDistance: return True
+
+    return False
+
+
+class Plane:
+    @staticmethod
+    def XY():
+        p1 = Vector((0, 0, 0))
+        p2 = Vector((1, 0, 0))
+        p3 = Vector((0, 1, 0))
+
+        return Plane(p1, p2, p3)
+
+
+    # plane equation: (p - position).dot(normal) = 0
+    def __init__(self, P1, P2, P3):
+        self.normal = (P2 - P1).cross(P3 - P1)
+        self.normal.normalize()
+
+        self.position = P1
+
+
+    def CalcIntersectionPointLineSegment(self, PL1, PL2):
+        DL = PL2 - PL1
+
+        try: rvPar = ((self.position - PL1).dot(self.normal)) / (DL.dot(self.normal))
+        except: return None
+
+        return rvPar
+
+
+    def CalcNormalParameter(self, vector):
+        return (vector - self.position).dot(self.normal)
+
+
+    def CalcProjection(self, vector):
+        normalParameter = self.CalcNormalParameter(vector)
+
+        rvv3 = vector - (self.normal * normalParameter)
+
+        return [normalParameter, rvv3]
+
+
+
+# http://geomalgorithms.com/a07-_distance.html
+def CalcClosestPointLineSegments(v3P0, v3P1, v3Q0, v3Q1):
+    u = v3P1 - v3P0
+    v = v3Q1 - v3Q0
+
+    w0 = v3P0 - v3Q0
+    a = u.dot(u)
+    b = u.dot(v)
+    c = v.dot(v)
+    d = u.dot(w0)
+    e = v.dot(w0)
+
+
+    try: parP = (b * e - c * d) / (a * c - b * b)
+    except: return None
+
+    try: parQ = (a * e - b * d) / (a * c - b * b)
+    except: return None
+
+
+    return [parP, parQ]
+
+
+def CalcIntersectionPointLineSegments(v3P0, v3P1, v3Q0, v3Q1, limitDistance):
+    rvList = CalcClosestPointLineSegments(v3P0, v3P1, v3Q0, v3Q1)
+    if rvList is None: return None
+
+
+    parP = rvList[0]
+    if parP < 0.0: return None
+    if parP > 1.0: return None
+
+    parQ = rvList[1]
+    if parQ < 0.0: return None
+    if parQ > 1.0: return None
+
+
+    pointP = v3P0 + ((v3P1 - v3P0) * parP)
+    pointQ = v3Q0 + ((v3Q1 - v3Q0) * parQ)
+    if not IsSamePoint(pointP, pointQ, limitDistance): return None
+
+    return [parP, parQ, pointP, pointQ]
+
+
+def CalcIntersectionPointsLineSegmentsPOV(v3P0, v3P1, v3Q0, v3Q1, v3POV):
+    planeQ = Plane(v3POV, v3Q0, v3Q1)
+    parP = planeQ.CalcIntersectionPointLineSegment(v3P0, v3P1)
+    if parP is None: return None
+    if parP < 0.0: return None
+    if parP > 1.0: return None
+
+    planeP = Plane(v3POV, v3P0, v3P1)
+    parQ = planeP.CalcIntersectionPointLineSegment(v3Q0, v3Q1)
+    if parQ is None: return None
+    if parQ < 0.0: return None
+    if parQ > 1.0: return None
+
+    return [parP, parQ]
+
+
+def CalcIntersectionPointsLineSegmentsDIR(v3P0, v3P1, v3Q0, v3Q1, v3DIR):
+    v3POV = v3Q0 + v3DIR
+    planeQ = Plane(v3POV, v3Q0, v3Q1)
+    parP = planeQ.CalcIntersectionPointLineSegment(v3P0, v3P1)
+    if parP is None: return None
+    if parP < 0.0: return None
+    if parP > 1.0: return None
+
+    v3POV = v3P0 + v3DIR
+    planeP = Plane(v3POV, v3P0, v3P1)
+    parQ = planeP.CalcIntersectionPointLineSegment(v3Q0, v3Q1)
+    if parQ is None: return None
+    if parQ < 0.0: return None
+    if parQ > 1.0: return None
+
+    return [parP, parQ]
+
+
+
+def CalcRotationMatrix(v3From, v3To):
+    cross = v3From.cross(v3To)
+
+    try: angle = v3From.angle(v3To)
+    except: return Matrix.Identity(4)
+
+    return Matrix.Rotation(angle, 4, cross)        # normalize axis?
+
+
+def subdivide_cubic_bezier(p1, p2, p3, p4, t):
+    p12 = (p2 - p1) * t + p1
+    p23 = (p3 - p2) * t + p2
+    p34 = (p4 - p3) * t + p3
+    p123 = (p23 - p12) * t + p12
+    p234 = (p34 - p23) * t + p23
+    p1234 = (p234 - p123) * t + p123
+    return [p12, p123, p1234, p234, p34]
diff --git a/curve_tools/Operators.py b/curve_tools/Operators.py
new file mode 100644
index 00000000..6d7bd81d
--- /dev/null
+++ b/curve_tools/Operators.py
@@ -0,0 +1,634 @@
+import time
+import threading
+
+import bpy
+from bpy.props import *
+from bpy_extras import object_utils, view3d_utils
+from mathutils import  *
+from math import  *
+
+from . import Properties
+from . import Curves
+from . import CurveIntersections
+from . import Util
+from . import Surfaces
+from . import Math
+
+# 1 CURVE SELECTED
+# ################
+class OperatorCurveInfo(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorcurveinfo"
+    bl_label = "Info"
+    bl_description = "Displays general info about the active/selected curve"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1Curve()
+
+
+    def execute(self, context):
+        curve = Curves.Curve(context.active_object)
+
+        nrSplines = len(curve.splines)
+        nrSegments = 0
+        nrEmptySplines = 0
+        for spline in curve.splines:
+            nrSegments += spline.nrSegments
+            if spline.nrSegments < 1: nrEmptySplines += 1
+
+
+        self.report({'INFO'}, "nrSplines: %d; nrSegments: %d; nrEmptySplines: %d" % (nrSplines, nrSegments, nrEmptySplines))
+
+        return {'FINISHED'}
+
+
+
+class OperatorCurveLength(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorcurvelength"
+    bl_label = "Length"
+    bl_description = "Calculates the length of the active/selected curve"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1Curve()
+
+
+    def execute(self, context):
+        curve = Curves.Curve(context.active_object)
+
+        context.scene.curvetools.CurveLength = curve.length
+
+        return {'FINISHED'}
+
+
+
+class OperatorSplinesInfo(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsplinesinfo"
+    bl_label = "Info"
+    bl_description = "Displays general info about the splines of the active/selected curve"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1Curve()
+
+
+    def execute(self, context):
+        curve = Curves.Curve(context.active_object)
+        nrSplines = len(curve.splines)
+
+        print("")
+        print("OperatorSplinesInfo:", "nrSplines:", nrSplines)
+
+        nrEmptySplines = 0
+        for iSpline, spline in enumerate(curve.splines):
+            print("--", "spline %d of %d: nrSegments: %d" % (iSpline + 1, nrSplines, spline.nrSegments))
+
+            if spline.nrSegments < 1:
+                nrEmptySplines += 1
+                print("--", "--", "## WARNING: spline has no segments and will therefor be ignored in any further calculations")
+
+
+        self.report({'INFO'}, "nrSplines: %d; nrEmptySplines: %d" % (nrSplines, nrEmptySplines) + " -- more info: see console")
+
+        return {'FINISHED'}
+
+
+
+class OperatorSegmentsInfo(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsegmentsinfo"
+    bl_label = "Info"
+    bl_description = "Displays general info about the segments of the active/selected curve"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1Curve()
+
+
+    def execute(self, context):
+        curve = Curves.Curve(context.active_object)
+        nrSplines = len(curve.splines)
+        nrSegments = 0
+
+        print("")
+        print("OperatorSegmentsInfo:", "nrSplines:", nrSplines)
+
+        nrEmptySplines = 0
+        for iSpline, spline in enumerate(curve.splines):
+            nrSegmentsSpline = spline.nrSegments
+            print("--", "spline %d of %d: nrSegments: %d" % (iSpline + 1, nrSplines, nrSegmentsSpline))
+
+            if nrSegmentsSpline < 1:
+                nrEmptySplines += 1
+                print("--", "--", "## WARNING: spline has no segments and will therefor be ignored in any further calculations")
+                continue
+
+            for iSegment, segment in enumerate(spline.segments):
+                print("--", "--", "segment %d of %d coefficients:" % (iSegment + 1, nrSegmentsSpline))
+                print("--", "--", "--", "C0: %.6f, %.6f, %.6f" % (segment.coeff0.x, segment.coeff0.y, segment.coeff0.z))
+
+            nrSegments += nrSegmentsSpline
+
+        self.report({'INFO'}, "nrSplines: %d; nrSegments: %d; nrEmptySplines: %d" % (nrSplines, nrSegments, nrEmptySplines))
+
+        return {'FINISHED'}
+
+
+
+class OperatorOriginToSpline0Start(bpy.types.Operator):
+    bl_idname = "curvetools2.operatororigintospline0start"
+    bl_label = "OriginToSpline0Start"
+    bl_description = "Sets the origin of the active/selected curve to the starting point of the (first) spline. Nice for curve modifiers."
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1Curve()
+
+
+    def execute(self, context):
+        blCurve = context.active_object
+        blSpline = blCurve.data.splines[0]
+        newOrigin = blCurve.matrix_world @ blSpline.bezier_points[0].co
+
+        origOrigin = bpy.context.scene.cursor.location.copy()
+        print("--", "origOrigin: %.6f, %.6f, %.6f" % (origOrigin.x, origOrigin.y, origOrigin.z))
+        print("--", "newOrigin: %.6f, %.6f, %.6f" % (newOrigin.x, newOrigin.y, newOrigin.z))
+
+        bpy.context.scene.cursor.location = newOrigin
+        bpy.ops.object.origin_set(type='ORIGIN_CURSOR')
+        bpy.context.scene.cursor.location = origOrigin
+
+        self.report({'INFO'}, "TODO: OperatorOriginToSpline0Start")
+
+        return {'FINISHED'}
+
+
+
+# 2 CURVES SELECTED
+# #################
+class OperatorIntersectCurves(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorintersectcurves"
+    bl_label = "Intersect"
+    bl_description = "Intersects selected curves"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected2OrMoreCurves()
+
+
+    def execute(self, context):
+        print("### TODO: OperatorIntersectCurves.execute()")
+
+        algo = context.scene.curvetools.IntersectCurvesAlgorithm
+        print("-- algo:", algo)
+
+
+        mode = context.scene.curvetools.IntersectCurvesMode
+        print("-- mode:", mode)
+        # if mode == 'Split':
+            # self.report({'WARNING'}, "'Split' mode is not implemented yet -- <<STOPPING>>")
+            # return {'CANCELLED'}
+
+        affect = context.scene.curvetools.IntersectCurvesAffect
+        print("-- affect:", affect)
+
+        selected_objects = context.selected_objects
+        lenodjs = len(selected_objects)
+        print('lenodjs:', lenodjs)
+        for i in range(0, lenodjs):
+            for j in range(0, lenodjs):
+                if j != i:
+                    bpy.ops.object.select_all(action='DESELECT')
+                    selected_objects[i].select_set(True)
+                    selected_objects[j].select_set(True)
+        
+                    if selected_objects[i].type == 'CURVE' and selected_objects[j].type == 'CURVE':
+                        curveIntersector = CurveIntersections.CurvesIntersector.FromSelection()
+                        rvIntersectionNrs = curveIntersector.CalcAndApplyIntersections()
+
+                        self.report({'INFO'}, "Active curve points: %d; other curve points: %d" % (rvIntersectionNrs[0], rvIntersectionNrs[1]))
+        
+        for obj in selected_objects:
+            obj.select_set(True)
+        
+        return {'FINISHED'}
+
+
+class OperatorLoftCurves(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorloftcurves"
+    bl_label = "Loft"
+    bl_description = "Lofts selected curves"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected2Curves()
+
+
+    def execute(self, context):
+        #print("### TODO: OperatorLoftCurves.execute()")
+
+        loftedSurface = Surfaces.LoftedSurface.FromSelection()
+        loftedSurface.AddToScene()
+
+        self.report({'INFO'}, "OperatorLoftCurves.execute()")
+
+        return {'FINISHED'}
+
+
+
+class OperatorSweepCurves(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsweepcurves"
+    bl_label = "Sweep"
+    bl_description = "Sweeps the active curve along to other curve (rail)"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected2Curves()
+
+
+    def execute(self, context):
+        #print("### TODO: OperatorSweepCurves.execute()")
+
+        sweptSurface = Surfaces.SweptSurface.FromSelection()
+        sweptSurface.AddToScene()
+
+        self.report({'INFO'}, "OperatorSweepCurves.execute()")
+
+        return {'FINISHED'}
+
+
+
+# 3 CURVES SELECTED
+# #################
+class OperatorBirail(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorbirail"
+    bl_label = "Birail"
+    bl_description = "Generates a birailed surface from 3 selected curves -- in order: rail1, rail2 and profile"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected3Curves()
+
+
+    def execute(self, context):
+        birailedSurface = Surfaces.BirailedSurface.FromSelection()
+        birailedSurface.AddToScene()
+
+        self.report({'INFO'}, "OperatorBirail.execute()")
+
+        return {'FINISHED'}
+
+
+
+# 1 OR MORE CURVES SELECTED
+# #########################
+class OperatorSplinesSetResolution(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsplinessetresolution"
+    bl_label = "SplinesSetResolution"
+    bl_description = "Sets the resolution of all splines"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1OrMoreCurves()
+
+
+    def execute(self, context):
+        splRes = context.scene.curvetools.SplineResolution
+        selCurves = Util.GetSelectedCurves()
+
+        for blCurve in selCurves:
+            for spline in blCurve.data.splines:
+                spline.resolution_u = splRes
+
+        return {'FINISHED'}
+
+
+
+class OperatorSplinesRemoveZeroSegment(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsplinesremovezerosegment"
+    bl_label = "SplinesRemoveZeroSegment"
+    bl_description = "Removes splines with no segments -- they seem to creep up, sometimes.."
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1OrMoreCurves()
+
+
+    def execute(self, context):
+        selCurves = Util.GetSelectedCurves()
+
+        for blCurve in selCurves:
+            curve = Curves.Curve(blCurve)
+            nrSplines = curve.nrSplines
+
+            splinesToRemove = []
+            for spline in curve.splines:
+                if len(spline.segments) < 1: splinesToRemove.append(spline)
+            nrRemovedSplines = len(splinesToRemove)
+
+            for spline in splinesToRemove: curve.splines.remove(spline)
+
+            if nrRemovedSplines > 0: curve.RebuildInScene()
+
+            self.report({'INFO'}, "Removed %d of %d splines" % (nrRemovedSplines, nrSplines))
+
+        return {'FINISHED'}
+
+
+
+class OperatorSplinesRemoveShort(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsplinesremoveshort"
+    bl_label = "SplinesRemoveShort"
+    bl_description = "Removes splines with a length smaller than the threshold"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1OrMoreCurves()
+
+
+    def execute(self, context):
+        threshold = context.scene.curvetools.SplineRemoveLength
+        selCurves = Util.GetSelectedCurves()
+
+        for blCurve in selCurves:
+            curve = Curves.Curve(blCurve)
+            nrSplines = curve.nrSplines
+
+            nrRemovedSplines = curve.RemoveShortSplines(threshold)
+            if nrRemovedSplines > 0: curve.RebuildInScene()
+
+            self.report({'INFO'}, "Removed %d of %d splines" % (nrRemovedSplines, nrSplines))
+
+        return {'FINISHED'}
+
+
+
+class OperatorSplinesJoinNeighbouring(bpy.types.Operator):
+    bl_idname = "curvetools2.operatorsplinesjoinneighbouring"
+    bl_label = "SplinesJoinNeighbouring"
+    bl_description = "Joins neighbouring splines within a distance smaller than the threshold"
+
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1OrMoreCurves()
+
+
+    def execute(self, context):
+        selCurves = Util.GetSelectedCurves()
+
+        for blCurve in selCurves:
+            curve = Curves.Curve(blCurve)
+            nrSplines = curve.nrSplines
+
+            threshold = context.scene.curvetools.SplineJoinDistance
+            startEnd = context.scene.curvetools.SplineJoinStartEnd
+            mode = context.scene.curvetools.SplineJoinMode
+
+            nrJoins = curve.JoinNeighbouringSplines(startEnd, threshold, mode)
+            if nrJoins > 0: curve.RebuildInScene()
+
+            self.report({'INFO'}, "Applied %d joins on %d splines; resulting nrSplines: %d" % (nrJoins, nrSplines, curve.nrSplines))
+
+        return {'FINISHED'}
+        
+def SurfaceFromBezier(surfacedata, points, center):
+    
+    len_points = len(points) - 1
+    
+    if len_points % 2 == 0:
+        h = Math.subdivide_cubic_bezier(
+                        points[len_points].co, points[len_points].handle_right,
+                        points[0].handle_left, points[0].co, 0.5
+                        )
+        points.add(1)
+        len_points = len(points) - 1
+        points[len_points - 1].handle_right = h[0]
+        points[len_points].handle_left = h[1]
+        points[len_points].co =  h[2]
+        points[len_points].handle_right = h[3]
+        points[0].handle_left =  h[4]
+        
+    half = round((len_points + 1)/2) - 1
+    # 1
+    surfacespline1 = surfacedata.splines.new(type='NURBS')
+    surfacespline1.points.add(3)
+    surfacespline1.points[0].co = [points[0].co.x, points[0].co.y, points[0].co.z, 1]
+    surfacespline1.points[1].co = [points[0].handle_left.x, points[0].handle_left.y, points[0].handle_left.z, 1]
+    surfacespline1.points[2].co = [points[len_points].handle_right.x,points[len_points].handle_right.y, points[len_points].handle_right.z, 1]
+    surfacespline1.points[3].co = [points[len_points].co.x, points[len_points].co.y, points[len_points].co.z, 1]
+    for p in surfacespline1.points:
+        p.select = True
+    surfacespline1.use_endpoint_u = True
+    surfacespline1.use_endpoint_v = True
+
+    for i in range(0, half):
+     
+        if center:
+            # 2
+            surfacespline2 = surfacedata.splines.new(type='NURBS')
+            surfacespline2.points.add(3)
+            surfacespline2.points[0].co = [points[i].co.x, points[i].co.y, points[i].co.z, 1]
+            surfacespline2.points[1].co = [(points[i].co.x + points[len_points - i].co.x)/2,
+                                           (points[i].co.y + points[len_points - i].co.y)/2,
+                                           (points[i].co.z + points[len_points - i].co.z)/2, 1]
+            surfacespline2.points[2].co = [(points[len_points - i].co.x + points[i].co.x)/2,
+                                           (points[len_points - i].co.y + points[i].co.y)/2,
+                                           (points[len_points - i].co.z + points[i].co.z)/2, 1]
+            surfacespline2.points[3].co = [points[len_points - i].co.x, points[len_points - i].co.y, points[len_points - i].co.z, 1]
+            for p in surfacespline2.points:
+                p.select = True
+            surfacespline2.use_endpoint_u = True
+            surfacespline2.use_endpoint_v = True
+        
+        # 3
+        surfacespline3 = surfacedata.splines.new(type='NURBS')
+        surfacespline3.points.add(3)
+        surfacespline3.points[0].co = [points[i].handle_right.x, points[i].handle_right.y, points[i].handle_right.z, 1]
+        surfacespline3.points[1].co = [(points[i].handle_right.x + points[len_points - i].handle_left.x)/2,
+                                       (points[i].handle_right.y + points[len_points - i].handle_left.y)/2,
+                                       (points[i].handle_right.z + points[len_points - i].handle_left.z)/2, 1]
+        surfacespline3.points[2].co = [(points[len_points - i].handle_left.x + points[i].handle_right.x)/2,
+                                       (points[len_points - i].handle_left.y + points[i].handle_right.y)/2,
+                                       (points[len_points - i].handle_left.z + points[i].handle_right.z)/2, 1]
+        surfacespline3.points[3].co = [points[len_points - i].handle_left.x, points[len_points - i].handle_left.y, points[len_points - i].handle_left.z, 1]
+        for p in surfacespline3.points:
+            p.select = True
+        surfacespline3.use_endpoint_u = True
+        surfacespline3.use_endpoint_v = True
+    
+        # 4
+        surfacespline4 = surfacedata.splines.new(type='NURBS')
+        surfacespline4.points.add(3)
+        surfacespline4.points[0].co = [points[i + 1].handle_left.x, points[i + 1].handle_left.y, points[i + 1].handle_left.z, 1]
+        surfacespline4.points[1].co = [(points[i + 1].handle_left.x + points[len_points - i - 1].handle_right.x)/2,
+                                       (points[i + 1].handle_left.y + points[len_points - i - 1].handle_right.y)/2,
+                                       (points[i + 1].handle_left.z + points[len_points - i - 1].handle_right.z)/2, 1]
+        surfacespline4.points[2].co = [(points[len_points - i - 1].handle_right.x + points[i + 1].handle_left.x)/2,
+                                       (points[len_points - i - 1].handle_right.y + points[i + 1].handle_left.y)/2,
+                                       (points[len_points - i - 1].handle_right.z + points[i + 1].handle_left.z)/2, 1]
+        surfacespline4.points[3].co = [points[len_points - i - 1].handle_right.x, points[len_points - i - 1].handle_right.y, points[len_points - i - 1].handle_right.z, 1]
+        for p in surfacespline4.points:
+            p.select = True
+        surfacespline4.use_endpoint_u = True
+        surfacespline4.use_endpoint_v = True
+        
+        if center:
+            # 5
+            surfacespline5 = surfacedata.splines.new(type='NURBS')
+            surfacespline5.points.add(3)
+            surfacespline5.points[0].co = [points[i + 1].co.x, points[i + 1].co.y, points[i + 1].co.z, 1]
+            surfacespline5.points[1].co = [(points[i + 1].co.x + points[len_points - i - 1].co.x)/2,
+                                           (points[i + 1].co.y + points[len_points - i - 1].co.y)/2,
+                                           (points[i + 1].co.z + points[len_points - i - 1].co.z)/2, 1]
+            surfacespline5.points[2].co = [(points[len_points - i - 1].co.x + points[i + 1].co.x)/2,
+                                           (points[len_points - i - 1].co.y + points[i + 1].co.y)/2,
+                                           (points[len_points - i - 1].co.z + points[i + 1].co.z)/2, 1]
+            surfacespline5.points[3].co = [points[len_points - i - 1].co.x, points[len_points - i - 1].co.y, points[len_points - i - 1].co.z, 1]
+            for p in surfacespline5.points:
+                p.select = True
+            surfacespline5.use_endpoint_u = True
+            surfacespline5.use_endpoint_v = True
+        
+    # 6
+    surfacespline6 = surfacedata.splines.new(type='NURBS')
+    surfacespline6.points.add(3)
+    surfacespline6.points[0].co = [points[half].co.x, points[half].co.y, points[half].co.z, 1]
+    surfacespline6.points[1].co = [points[half].handle_right.x, points[half].handle_right.y, points[half].handle_right.z, 1]
+    surfacespline6.points[2].co = [points[half+1].handle_left.x, points[half+1].handle_left.y, points[half+1].handle_left.z, 1]
+    surfacespline6.points[3].co = [points[half+1].co.x, points[half+1].co.y, points[half+1].co.z, 1]
+    for p in surfacespline6.points:
+        p.select = True
+    surfacespline6.use_endpoint_u = True
+    surfacespline6.use_endpoint_v = True
+    
+    bpy.ops.object.mode_set(mode = 'EDIT') 
+    bpy.ops.curve.make_segment()
+        
+    for s in surfacedata.splines:
+        s.resolution_u = 4
+        s.resolution_v = 4
+        s.order_u = 4
+        s.order_v = 4
+        for p in s.points:
+            p.select = False
+
+class ConvertSelectedFacesToBezier(bpy.types.Operator):
+    bl_idname = "curvetools2.convert_selected_face_to_bezier"
+    bl_label = "Convert selected faces to Bezier"
+    bl_description = "Convert selected faces to Bezier"
+    bl_options = {'REGISTER', 'UNDO'}
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected1Mesh()
+
+    def execute(self, context):
+        # main function
+        bpy.ops.object.mode_set(mode = 'OBJECT')
+        active_object = context.active_object
+        meshdata = active_object.data
+        curvedata = bpy.data.curves.new('Curve' + active_object.name, type='CURVE')
+        curveobject = object_utils.object_data_add(context, curvedata)
+        curvedata.dimensions = '3D'
+        
+        for poly in meshdata.polygons:
+            if poly.select:
+                newSpline = curvedata.splines.new(type='BEZIER')
+                newSpline.use_cyclic_u = True
+                newSpline.bezier_points.add(poly.loop_total - 1)
+                npoint = 0
+                for loop_index in range(poly.loop_start, poly.loop_start + poly.loop_total):
+                    newSpline.bezier_points[npoint].co = meshdata.vertices[meshdata.loops[loop_index].vertex_index].co
+                    newSpline.bezier_points[npoint].handle_left_type = 'VECTOR'
+                    newSpline.bezier_points[npoint].handle_right_type = 'VECTOR'
+                    newSpline.bezier_points[npoint].select_control_point = True
+                    newSpline.bezier_points[npoint].select_left_handle = True
+                    newSpline.bezier_points[npoint].select_right_handle = True
+                    npoint += 1
+                                  
+        return {'FINISHED'}
+        
+class ConvertBezierToSurface(bpy.types.Operator):
+    bl_idname = "curvetools2.convert_bezier_to_surface"
+    bl_label = "Convert Bezier to Surface"
+    bl_description = "Convert Bezier to Surface"
+    bl_options = {'REGISTER', 'UNDO'}
+
+    Center : BoolProperty(
+            name="Center",
+            default=False,
+            description="Consider center points"
+            )
+            
+    Resolution_U: IntProperty(
+            name="Resolution_U",
+            default=4,
+            min=1, max=64,
+            soft_min=1,
+            description="Surface resolution U"
+            )
+            
+    Resolution_V: IntProperty(
+            name="Resolution_V",
+            default=4,
+            min=1, max=64,
+            soft_min=1,
+            description="Surface resolution V"
+            )
+            
+    def draw(self, context):
+        layout = self.layout
+
+         # general options
+        col = layout.column()
+        col.prop(self, 'Center')
+        col.prop(self, 'Resolution_U')
+        col.prop(self, 'Resolution_V')
+    
+    @classmethod
+    def poll(cls, context):
+        return (context.object is not None and
+                context.object.type == 'CURVE')
+
+    def execute(self, context):
+        # main function
+        bpy.ops.object.mode_set(mode = 'OBJECT') 
+        active_object = context.active_object
+        curvedata = active_object.data
+        
+        surfacedata = bpy.data.curves.new('Surface', type='SURFACE')
+        surfaceobject = object_utils.object_data_add(context, surfacedata)
+        surfaceobject.matrix_world = active_object.matrix_world
+        surfaceobject.rotation_euler = active_object.rotation_euler
+        surfacedata.dimensions = '3D'
+        surfaceobject.show_wire = True
+        surfaceobject.show_in_front = True
+        
+        for spline in curvedata.splines:
+            SurfaceFromBezier(surfacedata, spline.bezier_points, self.Center)
+            
+        for spline in surfacedata.splines:
+            len_p = len(spline.points)
+            len_devide_4 = round(len_p / 4) + 1
+            len_devide_2 = round(len_p / 2)
+            bpy.ops.object.mode_set(mode = 'EDIT')
+            for point_index in range(len_devide_4, len_p - len_devide_4):
+                if point_index != len_devide_2 and point_index != len_devide_2 - 1:
+                    spline.points[point_index].select = True
+                
+            surfacedata.resolution_u = self.Resolution_U
+            surfacedata.resolution_v = self.Resolution_V
+
+        return {'FINISHED'}
diff --git a/curve_tools/PathFinder.py b/curve_tools/PathFinder.py
new file mode 100644
index 00000000..801c5293
--- /dev/null
+++ b/curve_tools/PathFinder.py
@@ -0,0 +1,329 @@
+# ##### BEGIN GPL LICENSE BLOCK #####
+#
+#  This program is free software; you can redistribute it and / or
+#  modify it under the terms of the GNU General Public License
+#  as published by the Free Software Foundation; either version 2
+#  of the License, or (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program; if not, write to the Free Software Foundation,
+#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+# ##### END GPL LICENSE BLOCK #####
+
+bl_info = {
+    'name': 'PathFinder',
+    'author': 'Spivak Vladimir (http://cwolf3d.korostyshev.net)',
+    'version': (0, 5, 0),
+    'blender': (2, 80, 0),
+    'location': 'Curve Tools addon. (N) Panel',
+    'description': 'PathFinder - quick search, selection, removal of splines',
+    'warning': '', # used for warning icon and text in addons panel
+    'wiki_url': '',
+    'tracker_url': '',
+    'category': 'Curve'}
+    
+import time
+import threading
+
+import gpu
+import bgl
+from gpu_extras.batch import batch_for_shader
+
+import bpy
+from bpy.props import *
+from bpy_extras import object_utils, view3d_utils
+from mathutils import  *
+from math import  *
+
+from . import Properties
+from . import Curves
+from . import CurveIntersections
+from . import Util
+from . import Surfaces
+from . import Math
+
+def get_bezier_points(spline, matrix_world):
+    point_list = []
+    len_bezier_points = len(spline.bezier_points)
+    if len_bezier_points > 1:
+        for i in range(0, len_bezier_points - 1):
+            point_list.extend([matrix_world @ spline.bezier_points[i].co])
+            for t in range(0, 100, 2):
+                h = Math.subdivide_cubic_bezier(spline.bezier_points[i].co,
+                                           spline.bezier_points[i].handle_right,
+                                           spline.bezier_points[i + 1].handle_left,
+                                           spline.bezier_points[i + 1].co,
+                                           t/100)
+                point_list.extend([matrix_world @ h[2]])
+        if spline.use_cyclic_u and len_bezier_points > 2:
+            point_list.extend([matrix_world @ spline.bezier_points[len_bezier_points - 1].co])
+            for t in range(0, 100, 2):
+                h = Math.subdivide_cubic_bezier(spline.bezier_points[len_bezier_points - 1].co,
+                                           spline.bezier_points[len_bezier_points - 1].handle_right,
+                                           spline.bezier_points[0].handle_left,
+                                           spline.bezier_points[0].co,
+                                           t/100)
+                point_list.extend([matrix_world @ h[2]])
+            point_list.extend([matrix_world @ spline.bezier_points[0].co])
+
+    return point_list
+        
+def get_points(spline, matrix_world):
+    point_list = []
+    len_points = len(spline.points)
+    if len_points > 1:
+        for i in range(0, len_points - 1):
+            point_list.extend([matrix_world @ Vector((spline.points[i].co.x, spline.points[i].co.y, spline.points[i].co.z))])
+            for t in range(0, 100, 2):
+                x = (spline.points[i].co.x + t / 100 * spline.points[i + 1].co.x) / (1 + t / 100)
+                y = (spline.points[i].co.y + t / 100 * spline.points[i + 1].co.y) / (1 + t / 100)
+                z = (spline.points[i].co.z + t / 100 * spline.points[i + 1].co.z) / (1 + t / 100)
+                point_list.extend([matrix_world @ Vector((x, y, z))])
+        if spline.use_cyclic_u and len_points > 2:
+            point_list.extend([matrix_world @ Vector((spline.points[len_points - 1].co.x, spline.points[len_points - 1].co.y, spline.points[len_points - 1].co.z))])
+            for t in range(0, 100, 2):
+                x = (spline.points[len_points - 1].co.x + t / 100 * spline.points[0].co.x) / (1 + t / 100)
+                y = (spline.points[len_points - 1].co.y + t / 100 * spline.points[0].co.y) / (1 + t / 100)
+                z = (spline.points[len_points - 1].co.z + t / 100 * spline.points[0].co.z) / (1 + t / 100)
+                point_list.extend([matrix_world @ Vector((x, y, z))])
+            point_list.extend([matrix_world @ Vector((spline.points[0].co.x, spline.points[0].co.y, spline.points[0].co.z))])
+    return point_list
+
+def draw_bezier_points(self, context, spline, matrix_world, path_color, path_thickness):
+    
+    points = get_bezier_points(spline, matrix_world)
+    
+    shader = gpu.shader.from_builtin('3D_UNIFORM_COLOR')
+    batch = batch_for_shader(shader, 'LINES', {"pos": points})
+    
+    shader.bind()
+    shader.uniform_float("color", path_color)
+    bgl.glEnable(bgl.GL_BLEND)
+    bgl.glLineWidth(path_thickness)
+    batch.draw(shader)
+    
+def draw_points(self, context, spline, matrix_world, path_color, path_thickness):
+    
+    points = get_points(spline, matrix_world)
+    
+    shader = gpu.shader.from_builtin('3D_UNIFORM_COLOR')
+    batch = batch_for_shader(shader, 'LINES', {"pos": points})
+    
+    shader.bind()
+    shader.uniform_float("color", path_color)
+    bgl.glEnable(bgl.GL_BLEND)
+    bgl.glLineWidth(path_thickness)
+    batch.draw(shader)
+
+def near(location3D, point, radius):
+    factor = 0
+    if point.x > (location3D.x - radius):
+        factor += 1
+    if point.x < (location3D.x + radius):
+        factor += 1
+    if point.y > (location3D.y - radius):
+        factor += 1
+    if point.y < (location3D.y + radius):
+        factor += 1
+    if point.z > (location3D.z - radius):
+        factor += 1
+    if point.z < (location3D.z + radius):
+        factor += 1
+        
+    return factor
+
+def click(self, context, event):
+    bpy.ops.object.mode_set(mode = 'EDIT')
+    bpy.context.view_layer.update()
+    for object in context.selected_objects:
+        matrix_world = object.matrix_world
+        if object.type == 'CURVE':
+            curvedata = object.data
+            
+            radius = bpy.context.scene.curvetools.PathFinderRadius
+            
+            for spline in curvedata.splines:
+                len_bezier_points = len(spline.bezier_points)
+                factor_max = 0
+                for i in range(0, len_bezier_points):
+
+                    co = matrix_world @ spline.bezier_points[i].co
+                    factor = near(self.location3D, co, radius)
+                    if factor > factor_max:
+                                factor_max = factor
+                        
+                    if i < len_bezier_points - 1:
+                        for t in range(0, 100, 2):
+                            h = Math.subdivide_cubic_bezier(spline.bezier_points[i].co,
+                                           spline.bezier_points[i].handle_right,
+                                           spline.bezier_points[i + 1].handle_left,
+                                           spline.bezier_points[i + 1].co,
+                                           t/100)
+                            co = matrix_world @ h[2]
+                            factor = near(self.location3D, co, radius)
+                            if factor > factor_max:
+                                factor_max = factor
+
+                    if spline.use_cyclic_u and len_bezier_points > 2:
+                        for t in range(0, 100, 2):
+                            h = Math.subdivide_cubic_bezier(spline.bezier_points[len_bezier_points - 1].co,
+                                           spline.bezier_points[len_bezier_points - 1].handle_right,
+                                           spline.bezier_points[0].handle_left,
+                                           spline.bezier_points[0].co,
+                                           t/100)
+                            co = matrix_world @ h[2]
+                            factor = near(self.location3D, co, radius)
+                            if factor > factor_max:
+                                factor_max = factor
+
+                if factor_max == 6:
+                    args = (self, context, spline, matrix_world, self.path_color, self.path_thickness)
+                    self.handlers.append(bpy.types.SpaceView3D.draw_handler_add(draw_bezier_points, args, 'WINDOW', 'POST_VIEW'))
+
+                    for bezier_point in spline.bezier_points:
+                        bezier_point.select_control_point = True
+                        bezier_point.select_left_handle = True
+                        bezier_point.select_right_handle = True
+            
+            for spline in curvedata.splines:
+                len_points = len(spline.points)
+                factor_max = 0
+                for i in range(0, len_points):
+                    co = matrix_world @ Vector((spline.points[i].co.x, spline.points[i].co.y, spline.points[i].co.z))
+                    factor = near(self.location3D, co, radius)
+                    if factor > factor_max:
+                        factor_max = factor
+
+                    if i < len_bezier_points - 1:
+                        for t in range(0, 100, 2):
+                            x = (spline.points[i].co.x + t / 100 * spline.points[i + 1].co.x) / (1 + t / 100)
+                            y = (spline.points[i].co.y + t / 100 * spline.points[i + 1].co.y) / (1 + t / 100)
+                            z = (spline.points[i].co.z + t / 100 * spline.points[i + 1].co.z) / (1 + t / 100)
+                            co = matrix_world @ Vector((x, y, z))
+                            factor = near(self.location3D, co, radius)
+                            if factor > factor_max:
+                                factor_max = factor
+
+                    if spline.use_cyclic_u and len_points > 2:
+                        for t in range(0, 100, 2):
+                            x = (spline.points[len_points - 1].co.x + t / 100 * spline.points[0].co.x) / (1 + t / 100)
+                            y = (spline.points[len_points - 1].co.y + t / 100 * spline.points[0].co.y) / (1 + t / 100)
+                            z = (spline.points[len_points - 1].co.z + t / 100 * spline.points[0].co.z) / (1 + t / 100)
+                            co = matrix_world @ Vector((x, y, z))
+                            factor = near(self.location3D, co, radius)
+                            if factor > factor_max:
+                                factor_max = factor
+
+                if factor_max == 6:
+                    args = (self, context, spline, matrix_world, self.path_color, self.path_thickness)
+                    self.handlers.append(bpy.types.SpaceView3D.draw_handler_add(draw_points, args, 'WINDOW', 'POST_VIEW'))
+                        
+                    for point in spline.points:
+                        point.select = True    
+
+class PathFinder(bpy.types.Operator):
+    bl_idname = "curvetools2.pathfinder"
+    bl_label = "Path Finder"
+    bl_description = "Path Finder"
+    bl_options = {'REGISTER', 'UNDO'}
+    
+    x: IntProperty(name="x", description="x")
+    y: IntProperty(name="y", description="y")
+    location3D: FloatVectorProperty(name = "",
+                description = "Start location",
+                default = (0.0, 0.0, 0.0),
+                subtype = 'XYZ')
+                
+    handlers = []
+    
+    def __init__(self):
+        self.report({'INFO'}, "ESC or TAB - cancel")
+
+    def __del__(self):
+        self.report({'INFO'}, "PathFinder deactivated")
+        
+    def execute(self, context):
+        bpy.ops.object.mode_set(mode = 'EDIT')
+        
+        # color change in the panel
+        self.path_color = bpy.context.scene.curvetools.path_color
+        self.path_thickness = bpy.context.scene.curvetools.path_thickness
+
+    def modal(self, context, event):
+        context.area.tag_redraw()
+        
+        if event.type in {'ESC', 'TAB'}:  # Cancel
+            for handler in self.handlers:
+                try:
+                    bpy.types.SpaceView3D.draw_handler_remove(handler, 'WINDOW')
+                except:
+                    pass
+            for handler in self.handlers:
+                self.handlers.remove(handler)
+            return {'CANCELLED'}
+            
+        if event.type in {'X', 'DEL'}:  # Cancel
+            for handler in self.handlers:
+                try:
+                    bpy.types.SpaceView3D.draw_handler_remove(handler, 'WINDOW')
+                except:
+                    pass
+            for handler in self.handlers:
+                self.handlers.remove(handler)
+            bpy.ops.curve.delete(type='VERT') 
+            return {'RUNNING_MODAL'}         
+        
+        elif event.alt and event.type == 'LEFTMOUSE':
+            click(self, context, event)
+                                    
+        elif event.alt and event.type == 'RIGHTMOUSE':
+           click(self, context, event)
+            
+        elif event.type == 'A':
+            for handler in self.handlers:
+                try:
+                    bpy.types.SpaceView3D.draw_handler_remove(handler, 'WINDOW')
+                except:
+                    pass
+            for handler in self.handlers:
+                self.handlers.remove(handler)
+            bpy.ops.curve.select_all(action='DESELECT')
+            
+        elif event.type == 'MOUSEMOVE':  # 
+            self.x = event.mouse_x
+            self.y = event.mouse_y
+            region = bpy.context.region
+            rv3d = bpy.context.space_data.region_3d
+            self.location3D = view3d_utils.region_2d_to_location_3d(
+                region,
+                rv3d,
+                (event.mouse_region_x, event.mouse_region_y),
+                (0.0, 0.0, 0.0)
+                )       
+
+        return {'PASS_THROUGH'}
+
+    def invoke(self, context, event):
+        self.execute(context)
+        context.window_manager.modal_handler_add(self)
+        return {'RUNNING_MODAL'}
+        
+    @classmethod
+    def poll(cls, context):
+        return (context.object is not None and
+                context.object.type == 'CURVE')
+
+def register():
+    bpy.utils.register_class(PathFinder)
+
+def unregister():
+    bpy.utils.unregister_class(PathFinder)
+
+if __name__ == "__main__":
+    register()
\ No newline at end of file
diff --git a/curve_tools/Properties.py b/curve_tools/Properties.py
new file mode 100644
index 00000000..37cf72cc
--- /dev/null
+++ b/curve_tools/Properties.py
@@ -0,0 +1,89 @@
+import time
+
+import bpy
+from bpy.props import *
+
+
+
+class CurveTools2SelectedObjectHeader(bpy.types.Header):
+    bl_label = "Selection"
+    bl_space_type = "VIEW_3D"
+
+    def __init__(self):
+        self.update()
+
+
+    def update(self):
+        blenderSelectedObjects = bpy.context.selected_objects
+        selectedObjects = bpy.context.scene.curvetools.SelectedObjects
+
+        selectedObjectsToRemove = []
+        for selectedObject in selectedObjects:
+            if not selectedObject.IsElementOf(blenderSelectedObjects): selectedObjectsToRemove.append(selectedObject)
+        for selectedObject in selectedObjectsToRemove: selectedObjects.remove(selectedObject)
+
+        blenderObjectsToAdd = []
+        for blenderObject in blenderSelectedObjects:
+            if not CurveTools2SelectedObject.ListContains(selectedObjects, blenderObject): blenderObjectsToAdd.append(blenderObject)
+        for blenderObject in blenderObjectsToAdd:
+            newSelectedObject = CurveTools2SelectedObject(blenderObject)
+            selectedObjects.append(newSelectedObject)
+
+
+    def draw(self, context):
+        selectedObjects = bpy.context.scene.curvetools.SelectedObjects
+        nrSelectedObjects = len(selectedObjects)
+
+        layout = self.layout
+        row = layout.row()
+        row.label(text="Sel: " + str(nrSelectedObjects))
+
+
+class CurveTools2SelectedObject(bpy.types.PropertyGroup):
+    name: StringProperty(name = "name", default = "??")
+
+
+    @staticmethod
+    def UpdateThreadTarget(lock, sleepTime, selectedObjectNames, selectedBlenderObjectNames):
+        time.sleep(sleepTime)
+
+        newSelectedObjectNames = []
+
+        for name in selectedObjectNames:
+            if name in selectedBlenderObjectNames: newSelectedObjectNames.append(name)
+
+        for name in selectedBlenderObjectNames:
+            if not (name in selectedObjectNames): newSelectedObjectNames.append(name)
+
+        # sometimes it still complains about the context
+        try:
+            nrNewSelectedObjects = len(newSelectedObjectNames)
+            bpy.context.scene.curvetools.NrSelectedObjects = nrNewSelectedObjects
+
+            selectedObjects = bpy.context.scene.curvetools.SelectedObjects
+            selectedObjects.clear()
+            for i in range(nrNewSelectedObjects): selectedObjects.add()
+            for i, newSelectedObjectName in enumerate(newSelectedObjectNames):
+                selectedObjects[i].name = newSelectedObjectName
+        except: pass
+
+
+    @staticmethod
+    def GetSelectedObjectNames():
+        selectedObjects = bpy.context.scene.curvetools.SelectedObjects
+
+        rvNames = []
+        selectedObjectValues = selectedObjects.values()
+        for selectedObject in selectedObjectValues: rvNames.append(selectedObject.name)
+
+        return rvNames
+
+
+    @staticmethod
+    def GetSelectedBlenderObjectNames():
+        blenderSelectedObjects = bpy.context.selected_objects
+
+        rvNames = []
+        for blObject in blenderSelectedObjects: rvNames.append(blObject.name)
+
+        return rvNames
diff --git a/curve_tools/ShowCurveResolution.py b/curve_tools/ShowCurveResolution.py
new file mode 100644
index 00000000..d6fc9c22
--- /dev/null
+++ b/curve_tools/ShowCurveResolution.py
@@ -0,0 +1,128 @@
+# ##### BEGIN GPL LICENSE BLOCK #####
+#
+#  This program is free software; you can redistribute it and / or
+#  modify it under the terms of the GNU General Public License
+#  as published by the Free Software Foundation; either version 2
+#  of the License, or (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program; if not, write to the Free Software Foundation,
+#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110 - 1301, USA.
+#
+# ##### END GPL LICENSE BLOCK #####
+#
+
+
+# LOAD MODUL #
+import bpy
+from bpy import *
+from bpy.props import *
+from bpy.types import AddonPreferences, PropertyGroup
+
+import bgl
+import blf
+import gpu
+from gpu_extras.batch import batch_for_shader
+
+import math
+import mathutils
+from mathutils import Vector
+from mathutils.geometry import interpolate_bezier
+
+import bpy_extras
+from bpy_extras.view3d_utils import location_3d_to_region_2d as loc3d2d
+
+ 
+ 
+def get_points(spline):
+
+    bezier_points = spline.bezier_points
+
+    if len(bezier_points) < 2: 
+        return []
+
+    r = spline.resolution_u + 1
+    segments = len(bezier_points)
+    
+    if not spline.use_cyclic_u:
+        segments -= 1
+ 
+    point_list = []
+    for i in range(segments):
+        inext = (i + 1) % len(bezier_points)
+ 
+        bezier_points1 = bezier_points[i].co
+        handle1 = bezier_points[i].handle_right
+        handle2 = bezier_points[inext].handle_left
+        bezier_points2 = bezier_points[inext].co
+        
+        bezier = bezier_points1, handle1, handle2, bezier_points2, r
+        points = interpolate_bezier(*bezier)
+        point_list.extend(points)
+ 
+    return point_list
+ 
+def draw(self, context, spline, curve_vertcolor):
+    
+    points = get_points(spline)
+    
+    shader = gpu.shader.from_builtin('3D_UNIFORM_COLOR')
+    batch = batch_for_shader(shader, 'POINTS', {"pos": points})
+    
+    shader.bind()
+    shader.uniform_float("color", curve_vertcolor)
+    batch.draw(shader)
+
+
+class ShowCurveResolution(bpy.types.Operator):
+    bl_idname = "curve.show_resolution"
+    bl_label = "Show Curve Resolution"
+    bl_description = "Show curve Resolution / [ESC] - remove"
+    
+    handlers = []
+    
+    def modal(self, context, event):
+        context.area.tag_redraw()
+
+        if event.type in {'ESC'}: 
+            for handler in self.handlers:
+                try:
+                    bpy.types.SpaceView3D.draw_handler_remove(handler, 'WINDOW')
+                except:
+                    pass
+            for handler in self.handlers:
+                self.handlers.remove(handler)
+            return {'CANCELLED'}
+
+        return {'PASS_THROUGH'}
+
+
+    def invoke(self, context, event):
+
+        if context.area.type == 'VIEW_3D':
+        
+            # color change in the panel
+            curve_vertcolor = bpy.context.scene.curvetools.curve_vertcolor
+            
+            # the arguments we pass the the callback
+            
+            
+            splines = context.active_object.data.splines
+            for spline in splines:
+                args = (self, context, spline, curve_vertcolor)
+            
+                # Add the region OpenGL drawing callback
+                # draw in view space with 'POST_VIEW' and 'PRE_VIEW'
+                self.handlers.append(bpy.types.SpaceView3D.draw_handler_add(draw, args, 'WINDOW', 'POST_VIEW'))
+
+            context.window_manager.modal_handler_add(self)
+            return {'RUNNING_MODAL'}
+        else:
+            self.report({'WARNING'}, 
+            "View3D not found, cannot run operator")
+            return {'CANCELLED'}
diff --git a/curve_tools/Surfaces.py b/curve_tools/Surfaces.py
new file mode 100644
index 00000000..22b5119c
--- /dev/null
+++ b/curve_tools/Surfaces.py
@@ -0,0 +1,462 @@
+import bpy
+import bmesh
+
+from . import Math
+from . import Curves
+
+
+
+class LoftedSplineSurface:
+    def __init__(self, activeSpline, otherSpline, bMesh, vert0Index, resolution):
+        self.splineA = activeSpline
+        self.splineO = otherSpline
+
+        self.bMesh = bMesh
+        self.vert0Index = vert0Index
+        self.resolution = resolution
+
+
+    def Apply(self, worldMatrixA, worldMatrixO):
+        #deltaPar = 1.0 / float(self.resolution - 1)
+
+        par = 0.0
+        pointA = worldMatrixA @ self.splineA.CalcPoint(par)
+        pointO = worldMatrixO @ self.splineO.CalcPoint(par)
+        self.bMesh.verts[self.vert0Index].co = pointA
+        self.bMesh.verts[self.vert0Index + 1].co = pointO
+
+        fltResm1 = float(self.resolution - 1)
+        for i in range(1, self.resolution):
+            par = float(i) / fltResm1
+
+            pointA = worldMatrixA @ self.splineA.CalcPoint(par)
+            pointO = worldMatrixO @ self.splineO.CalcPoint(par)
+            self.bMesh.verts[self.vert0Index + 2 * i].co = pointA
+            self.bMesh.verts[self.vert0Index + 2 * i + 1].co = pointO
+
+
+    def AddFaces(self):
+        currIndexA = self.vert0Index
+        currIndexO = self.vert0Index + 1
+
+        bmVerts = self.bMesh.verts
+        bmVerts.ensure_lookup_table()
+
+        for i in range(1, self.resolution):
+            nextIndexA = self.vert0Index + 2 * i
+            nextIndexO = nextIndexA + 1
+
+            self.bMesh.faces.new([bmVerts[currIndexA], bmVerts[currIndexO], bmVerts[nextIndexO], bmVerts[nextIndexA]])
+
+            currIndexA = nextIndexA
+            currIndexO = nextIndexO
+
+
+class LoftedSurface:
+    @staticmethod
+    def FromSelection():
+        selObjects = bpy.context.selected_objects
+        if len(selObjects) != 2: raise Exception("len(selObjects) != 2") # shouldn't be possible
+
+        blenderActiveCurve = bpy.context.active_object
+        blenderOtherCurve = selObjects[0]
+        if blenderActiveCurve == blenderOtherCurve: blenderOtherCurve = selObjects[1]
+
+        aCurve = Curves.Curve(blenderActiveCurve)
+        oCurve = Curves.Curve(blenderOtherCurve)
+
+        name = "TODO: autoname"
+
+        return LoftedSurface(aCurve, oCurve, name)
+
+
+    def __init__(self, activeCurve, otherCurve, name = "LoftedSurface"):
+        self.curveA = activeCurve
+        self.curveO = otherCurve
+        self.name  = name
+
+        self.nrSplines = self.curveA.nrSplines
+        if self.curveO.nrSplines < self.nrSplines: self.nrSplines = self.curveO.nrSplines
+
+        self.bMesh = bmesh.new()
+
+        self.splineSurfaces = self.SetupSplineSurfaces()
+
+        self.Apply()
+
+
+    def SetupSplineSurfaces(self):
+        rvSplineSurfaces = []
+
+        currV0Index = 0
+        for i in range(self.nrSplines):
+            splineA = self.curveA.splines[i]
+            splineO = self.curveO.splines[i]
+
+            res = splineA.resolution
+            if splineO.resolution < res: res = splineO.resolution
+
+            for iv in range(2 * res): self.bMesh.verts.new()
+
+            splSurf = LoftedSplineSurface(splineA, splineO, self.bMesh, currV0Index, res)
+            splSurf.AddFaces()
+            rvSplineSurfaces.append(splSurf)
+
+            currV0Index += 2 * res
+
+        return rvSplineSurfaces
+
+
+    def Apply(self):
+        for splineSurface in self.splineSurfaces: splineSurface.Apply(self.curveA.worldMatrix, self.curveO.worldMatrix)
+
+
+    def AddToScene(self):
+        mesh = bpy.data.meshes.new("Mesh" + self.name)
+
+        self.bMesh.to_mesh(mesh)
+        mesh.update()
+
+        meshObject = bpy.data.objects.new(self.name, mesh)
+
+        bpy.context.collection.objects.link(meshObject)
+
+
+
+# active spline is swept over other spline (rail)
+class SweptSplineSurface:
+    def __init__(self, activeSpline, otherSpline, bMesh, vert0Index, resolutionA, resolutionO):
+        self.splineA = activeSpline
+        self.splineO = otherSpline
+
+        self.bMesh = bMesh
+        self.vert0Index = vert0Index
+        self.resolutionA = resolutionA
+        self.resolutionO = resolutionO
+
+
+    def Apply(self, worldMatrixA, worldMatrixO):
+        localPointsA = []
+        fltResAm1 = float(self.resolutionA - 1)
+        for i in range(self.resolutionA):
+            par = float(i) / fltResAm1
+            pointA = self.splineA.CalcPoint(par)
+            localPointsA.append(pointA)
+
+
+        worldPointsO = []
+        localDerivativesO = []
+        fltResOm1 = float(self.resolutionO - 1)
+        for i in range(self.resolutionO):
+            par = float(i) / fltResOm1
+
+            pointO = self.splineO.CalcPoint(par)
+            worldPointsO.append(worldMatrixO @ pointO)
+
+            derivativeO = self.splineO.CalcDerivative(par)
+            localDerivativesO.append(derivativeO)
+
+
+        currWorldMatrixA = worldMatrixA
+        worldMatrixOInv = worldMatrixO.inverted()
+        prevDerivativeO = localDerivativesO[0]
+        for iO in range(self.resolutionO):
+            currDerivativeO = localDerivativesO[iO]
+            localRotMatO = Math.CalcRotationMatrix(prevDerivativeO, currDerivativeO)
+
+            currLocalAToLocalO = worldMatrixOInv @ currWorldMatrixA
+            worldPointsA = []
+            for iA in range(self.resolutionA):
+                pointALocalToO = currLocalAToLocalO @ localPointsA[iA]
+                rotatedPointA = localRotMatO @ pointALocalToO
+                worldPointsA.append(worldMatrixO @ rotatedPointA)
+
+            worldOffsetsA = []
+            worldPoint0A = worldPointsA[0]
+            for i in range(self.resolutionA): worldOffsetsA.append(worldPointsA[i] - worldPoint0A)
+
+
+            for iA in range(self.resolutionA):
+                iVert = self.vert0Index + (self.resolutionA * iO) + iA
+                currVert = worldPointsO[iO] + worldOffsetsA[iA]
+                self.bMesh.verts[iVert].co = currVert
+
+            prevDerivativeO = currDerivativeO
+            currWorldMatrixA = worldMatrixO @ localRotMatO @ currLocalAToLocalO
+
+
+    def AddFaces(self):
+        bmVerts = self.bMesh.verts
+        bmVerts.ensure_lookup_table()
+
+        for iO in range(self.resolutionO - 1):
+            for iA in range(self.resolutionA - 1):
+                currIndexA1 = self.vert0Index + (self.resolutionA * iO) + iA
+                currIndexA2 = currIndexA1 + 1
+                nextIndexA1 = self.vert0Index + (self.resolutionA * (iO + 1)) + iA
+                nextIndexA2 = nextIndexA1 + 1
+
+                self.bMesh.faces.new([bmVerts[currIndexA1], bmVerts[currIndexA2], bmVerts[nextIndexA2], bmVerts[nextIndexA1]])
+
+
+
+class SweptSurface:
+    @staticmethod
+    def FromSelection():
+        selObjects = bpy.context.selected_objects
+        if len(selObjects) != 2: raise Exception("len(selObjects) != 2") # shouldn't be possible
+
+        blenderActiveCurve = bpy.context.active_object
+        blenderOtherCurve = selObjects[0]
+        if blenderActiveCurve == blenderOtherCurve: blenderOtherCurve = selObjects[1]
+
+        aCurve = Curves.Curve(blenderActiveCurve)
+        oCurve = Curves.Curve(blenderOtherCurve)
+
+        name = "TODO: autoname"
+
+        return SweptSurface(aCurve, oCurve, name)
+
+
+    def __init__(self, activeCurve, otherCurve, name = "SweptSurface"):
+        self.curveA = activeCurve
+        self.curveO = otherCurve
+        self.name  = name
+
+        self.nrSplines = self.curveA.nrSplines
+        if self.curveO.nrSplines < self.nrSplines: self.nrSplines = self.curveO.nrSplines
+
+        self.bMesh = bmesh.new()
+
+        self.splineSurfaces = self.SetupSplineSurfaces()
+
+        self.Apply()
+
+
+    def SetupSplineSurfaces(self):
+        rvSplineSurfaces = []
+
+        currV0Index = 0
+        for i in range(self.nrSplines):
+            splineA = self.curveA.splines[i]
+            splineO = self.curveO.splines[i]
+
+            resA = splineA.resolution
+            resO = splineO.resolution
+
+            for iv in range(resA * resO): self.bMesh.verts.new()
+
+            splSurf = SweptSplineSurface(splineA, splineO, self.bMesh, currV0Index, resA, resO)
+            splSurf.AddFaces()
+            rvSplineSurfaces.append(splSurf)
+
+            currV0Index += resA * resO
+
+        return rvSplineSurfaces
+
+
+    def Apply(self):
+        for splineSurface in self.splineSurfaces: splineSurface.Apply(self.curveA.worldMatrix, self.curveO.worldMatrix)
+
+
+    def AddToScene(self):
+        mesh = bpy.data.meshes.new("Mesh" + self.name)
+
+        self.bMesh.to_mesh(mesh)
+        mesh.update()
+
+        meshObject = bpy.data.objects.new(self.name, mesh)
+
+        bpy.context.collection.objects.link(meshObject)
+
+
+
+# profileSpline is swept over rail1Spline and scaled/rotated to have its endpoint on rail2Spline
+class BirailedSplineSurface:
+    def __init__(self, rail1Spline, rail2Spline, profileSpline, bMesh, vert0Index, resolutionRails, resolutionProfile):
+        self.rail1Spline = rail1Spline
+        self.rail2Spline = rail2Spline
+        self.profileSpline = profileSpline
+
+        self.bMesh = bMesh
+        self.vert0Index = vert0Index
+        self.resolutionRails = resolutionRails
+        self.resolutionProfile = resolutionProfile
+
+
+    def Apply(self, worldMatrixRail1, worldMatrixRail2, worldMatrixProfile):
+        localPointsProfile = []
+        fltResProfilem1 = float(self.resolutionProfile - 1)
+        for i in range(self.resolutionProfile):
+            par = float(i) / fltResProfilem1
+            pointProfile = self.profileSpline.CalcPoint(par)
+            localPointsProfile.append(pointProfile)
+
+
+        worldPointsRail1 = []
+        localDerivativesRail1 = []
+        worldPointsRail2 = []
+        fltResRailsm1 = float(self.resolutionRails - 1)
+        for i in range(self.resolutionRails):
+            par = float(i) / fltResRailsm1
+
+            pointRail1 = self.rail1Spline.CalcPoint(par)
+            worldPointsRail1.append(worldMatrixRail1 @ pointRail1)
+
+            derivativeRail1 = self.rail1Spline.CalcDerivative(par)
+            localDerivativesRail1.append(derivativeRail1)
+
+            pointRail2 = self.rail2Spline.CalcPoint(par)
+            worldPointsRail2.append(worldMatrixRail2 @ pointRail2)
+
+
+        currWorldMatrixProfile = worldMatrixProfile
+        worldMatrixRail1Inv = worldMatrixRail1.inverted()
+        prevDerivativeRail1 = localDerivativesRail1[0]
+        for iRail in range(self.resolutionRails):
+            currDerivativeRail1 = localDerivativesRail1[iRail]
+            localRotMatRail1 = Math.CalcRotationMatrix(prevDerivativeRail1, currDerivativeRail1)
+
+            currLocalProfileToLocalRail1 = worldMatrixRail1Inv @ currWorldMatrixProfile
+            worldPointsProfileRail1 = []
+            for iProfile in range(self.resolutionProfile):
+                pointProfileLocalToRail1 = currLocalProfileToLocalRail1 @ localPointsProfile[iProfile]
+                rotatedPointProfile = localRotMatRail1 @ pointProfileLocalToRail1
+                worldPointsProfileRail1.append(worldMatrixRail1 @ rotatedPointProfile)
+
+            worldOffsetsProfileRail1 = []
+            worldPoint0ProfileRail1 = worldPointsProfileRail1[0]
+            for iProfile in range(self.resolutionProfile): worldOffsetsProfileRail1.append(worldPointsProfileRail1[iProfile] - worldPoint0ProfileRail1)
+
+            worldStartPointProfileRail1 = worldPointsRail1[iRail]
+            worldEndPointProfileRail1 = worldStartPointProfileRail1 + worldOffsetsProfileRail1[-1]
+            v3From = worldEndPointProfileRail1 - worldStartPointProfileRail1
+            v3To = worldPointsRail2[iRail] - worldStartPointProfileRail1
+            if not v3From.magnitude == 0:
+                scaleFactorRail2 = v3To.magnitude / v3From.magnitude
+            else:
+                scaleFactorRail2 = 1
+            rotMatRail2 = Math.CalcRotationMatrix(v3From, v3To)
+
+            worldOffsetsProfileRail2 = []
+            for iProfile in range(self.resolutionProfile):
+                offsetProfileRail1 = worldOffsetsProfileRail1[iProfile]
+                worldOffsetsProfileRail2.append(rotMatRail2 @ (offsetProfileRail1 * scaleFactorRail2))
+
+
+            for iProfile in range(self.resolutionProfile):
+                iVert = self.vert0Index + (self.resolutionProfile * iRail) + iProfile
+                currVert = worldPointsRail1[iRail] + worldOffsetsProfileRail2[iProfile]
+                self.bMesh.verts[iVert].co = currVert
+
+            prevDerivativeRail1 = currDerivativeRail1
+            currWorldMatrixProfile = worldMatrixRail1 @ localRotMatRail1 @ currLocalProfileToLocalRail1
+
+
+    def AddFaces(self):
+        bmVerts = self.bMesh.verts
+        bmVerts.ensure_lookup_table()
+
+        for iRail in range(self.resolutionRails - 1):
+            for iProfile in range(self.resolutionProfile - 1):
+                currIndex1 = self.vert0Index + (self.resolutionProfile * iRail) + iProfile
+                currIndex2 = currIndex1 + 1
+                nextIndex1 = self.vert0Index + (self.resolutionProfile * (iRail + 1)) + iProfile
+                nextIndex2 = nextIndex1 + 1
+
+                self.bMesh.faces.new([bmVerts[currIndex1], bmVerts[currIndex2], bmVerts[nextIndex2], bmVerts[nextIndex1]])
+
+
+
+class BirailedSurface:
+    @staticmethod
+    def FromSelection():
+        nrSelectedObjects = bpy.context.scene.curvetools.NrSelectedObjects
+        if nrSelectedObjects != 3: raise Exception("nrSelectedObjects != 3") # shouldn't be possible
+
+
+        selectedObjects = bpy.context.scene.curvetools.SelectedObjects
+        selectedObjectValues = selectedObjects.values()
+
+        curveName = selectedObjectValues[0].name
+        rail1BlenderCurve = None
+        try: rail1BlenderCurve = bpy.data.objects[curveName]
+        except: rail1BlenderCurve = None
+        if rail1BlenderCurve is None: raise Exception("rail1BlenderCurve is None")
+
+        curveName = selectedObjectValues[1].name
+        rail2BlenderCurve = None
+        try: rail2BlenderCurve = bpy.data.objects[curveName]
+        except: rail2BlenderCurve = None
+        if rail2BlenderCurve is None: raise Exception("rail2BlenderCurve is None")
+
+        curveName = selectedObjectValues[2].name
+        profileBlenderCurve = None
+        try: profileBlenderCurve = bpy.data.objects[curveName]
+        except: profileBlenderCurve = None
+        if profileBlenderCurve is None: raise Exception("profileBlenderCurve is None")
+
+
+        rail1Curve = Curves.Curve(rail1BlenderCurve)
+        rail2Curve = Curves.Curve(rail2BlenderCurve)
+        profileCurve = Curves.Curve(profileBlenderCurve)
+
+        name = "TODO: autoname"
+
+        return BirailedSurface(rail1Curve, rail2Curve, profileCurve, name)
+
+
+    def __init__(self, rail1Curve, rail2Curve, profileCurve, name = "BirailedSurface"):
+        self.rail1Curve = rail1Curve
+        self.rail2Curve = rail2Curve
+        self.profileCurve = profileCurve
+        self.name  = name
+
+        self.nrSplines = self.rail1Curve.nrSplines
+        if self.rail2Curve.nrSplines < self.nrSplines: self.nrSplines = self.rail2Curve.nrSplines
+        if self.profileCurve.nrSplines < self.nrSplines: self.nrSplines = self.profileCurve.nrSplines
+
+        self.bMesh = bmesh.new()
+
+        self.splineSurfaces = self.SetupSplineSurfaces()
+
+        self.Apply()
+
+
+    def SetupSplineSurfaces(self):
+        rvSplineSurfaces = []
+
+        currV0Index = 0
+        for i in range(self.nrSplines):
+            splineRail1 = self.rail1Curve.splines[i]
+            splineRail2 = self.rail2Curve.splines[i]
+            splineProfile = self.profileCurve.splines[i]
+
+            resProfile = splineProfile.resolution
+            resRails = splineRail1.resolution
+            if splineRail2.resolution < resRails: resRails = splineRail2.resolution
+
+            for iv in range(resProfile * resRails): self.bMesh.verts.new()
+
+            splSurf = BirailedSplineSurface(splineRail1, splineRail2, splineProfile, self.bMesh, currV0Index, resRails, resProfile)
+            splSurf.AddFaces()
+            rvSplineSurfaces.append(splSurf)
+
+            currV0Index += resProfile * resRails
+
+        return rvSplineSurfaces
+
+
+    def Apply(self):
+        for splineSurface in self.splineSurfaces: splineSurface.Apply(self.rail1Curve.worldMatrix, self.rail2Curve.worldMatrix, self.profileCurve.worldMatrix)
+
+
+    def AddToScene(self):
+        mesh = bpy.data.meshes.new("Mesh" + self.name)
+
+        self.bMesh.to_mesh(mesh)
+        mesh.update()
+
+        meshObject = bpy.data.objects.new(self.name, mesh)
+
+        bpy.context.collection.objects.link(meshObject)
diff --git a/curve_tools/Util.py b/curve_tools/Util.py
new file mode 100644
index 00000000..30596697
--- /dev/null
+++ b/curve_tools/Util.py
@@ -0,0 +1,178 @@
+import bpy
+from mathutils import *
+
+
+def GetSelectedCurves():
+    rvList = []
+
+    for obj in bpy.context.selected_objects:
+        try:
+            if obj.type == "CURVE": rvList.append(obj)
+        except:
+            pass
+
+    return rvList
+    
+    
+def GetSelectedMeshes():
+    rvList = []
+
+    for obj in bpy.context.selected_objects:
+        try:
+            if obj.type == "MESH": rvList.append(obj)
+        except:
+            pass
+
+    return rvList
+
+
+def Selected1Curve():
+    try:
+        if len(GetSelectedCurves()) == 1:
+            return (bpy.context.active_object.type == "CURVE")
+    except:
+        pass
+
+    return False
+
+
+def Selected1Mesh():
+    try:
+        if len(GetSelectedMeshes()) == 1:
+            return (bpy.context.active_object.type == "MESH")
+    except:
+        pass
+
+    return False
+
+
+def Selected1SingleSplineCurve():
+    try:
+        if Selected1Curve():
+            return (len(bpy.context.active_object.data.splines) == 1)
+    except:
+        pass
+
+    return False
+
+
+def Selected2Curves():
+    try:
+        if len(GetSelectedCurves()) == 2:
+            return (bpy.context.active_object.type == "CURVE")
+    except:
+        pass
+
+    return False
+
+
+def Selected3Curves():
+    try:
+        if len(GetSelectedCurves()) == 3:
+            return (bpy.context.active_object.type == "CURVE")
+    except:
+        pass
+
+    return False
+
+
+def Selected1OrMoreCurves():
+    try:
+        if len(GetSelectedCurves()) > 0:
+            return (bpy.context.active_object.type == "CURVE")
+    except:
+        pass
+
+    return False
+    
+def Selected2OrMoreCurves():
+    try:
+        if len(GetSelectedCurves()) > 1:
+            return (bpy.context.active_object.type == "CURVE")
+    except:
+        pass
+
+    return False
+
+
+def Selected1OrMoreMesh():
+    try:
+        if len(GetSelectedMeshes()) > 0:
+            return (bpy.context.active_object.type == "MESH")
+    except:
+        pass
+
+    return False
+
+
+def GetToolsRegion():
+    for area in bpy.context.screen.areas:
+        if area.type == 'VIEW_3D':
+            for region in area.regions:
+                if region.type == 'TOOLS': return region
+
+    return None
+
+
+def GetFirstRegionView3D():
+    for area in bpy.context.screen.areas:
+        if area.type == 'VIEW_3D':
+            return area.spaces[0].region_3d
+
+    return None
+
+
+def LogFirstRegionView3D():
+    print("LogFirstRegionView3D()")
+    regionView3D = GetFirstRegionView3D()
+    if regionView3D is None:
+        print("--", "ERROR:", "regionView3D is None")
+        return
+
+    print("--", "view_matrix:")
+    print("--", "--", regionView3D.view_matrix)
+    print("--", "view_location:")
+    print("--", "--", regionView3D.view_location)
+
+
+class Intersection:
+    # listIP: list of BezierSplineIntersectionPoint
+    # return: list of splines
+    @staticmethod
+    def GetBezierSplines(listIP):
+        rvList = []
+
+        for ip in listIP:
+            if not (ip.spline in rvList): rvList.append(ip.spline)
+
+        return rvList
+
+
+    # listIP: list of BezierSplineIntersectionPoint
+    # return: list of segments
+    @staticmethod
+    def GetBezierSegments(listIP, spline):
+        rvList = []
+
+        for ip in listIP:
+            if not ip.spline is spline: continue
+
+            segIP = ip.bezierSegmentIntersectionPoint
+            if not (segIP.segment in rvList): rvList.append(segIP.segment)
+
+        return rvList
+
+
+    # listIP: list of BezierSplineIntersectionPoint
+    # return: list of floats (not necessarily ordered)
+    @staticmethod
+    def GetBezierSegmentParameters(listIP, segment):
+        rvList = []
+
+        for ip in listIP:
+            segIP = ip.bezierSegmentIntersectionPoint
+            if not segIP.segment is segment: continue
+
+            rvList.append(segIP.parameter)
+
+        return rvList
diff --git a/curve_tools/__init__.py b/curve_tools/__init__.py
new file mode 100644
index 00000000..04d63c09
--- /dev/null
+++ b/curve_tools/__init__.py
@@ -0,0 +1,563 @@
+# ##### BEGIN GPL LICENSE BLOCK #####
+#
+#  This program is free software; you can redistribute it and/or
+#  modify it under the terms of the GNU General Public License
+#  as published by the Free Software Foundation; either version 2
+#  of the License, or (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program; if not, write to the Free Software Foundation,
+#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+# ##### END GPL LICENSE BLOCK #####
+# Contributed to by guy lateur, Alexander Meißner (Lichtso),
+# Dealga McArdle (zeffii), Marvin.K.Breuer (MKB),
+# Spivak Vladimir (cwolf3d)
+# Origunally an addon by Mackraken
+
+
+bl_info = {
+    "name": "Curve Tools 2",
+    "description": "Adds some functionality for bezier/nurbs curve/surface modeling",
+    "author": "Mackraken",
+    "version": (0, 3, 3),
+    "blender": (2, 80, 0),
+    "location": "View3D > Tool Shelf > Addons Tab",
+    "warning": "WIP",
+    "wiki_url": "https://wiki.blender.org/index.php/Extensions:2.6/Py/"
+                "Scripts/Curve/Curve_Tools",
+    "tracker_url": "https://developer.blender.org/maniphest/task/edit/form/2/",
+    "category": "Add Curve"}
+
+
+import os, bpy, importlib, math
+from bpy.types import (
+        Operator,
+        Panel,
+        PropertyGroup,
+        )
+from bpy.props import (
+        BoolProperty,
+        IntProperty,
+        FloatProperty,
+        EnumProperty,
+        CollectionProperty,
+        StringProperty,
+        FloatVectorProperty,
+        )
+from . import Properties
+from . import Operators
+from . import auto_loft
+from . import curve_outline
+from . import curve_remove_doubles
+from . import PathFinder
+from . import ShowCurveResolution
+from . import internal, cad, toolpath, exports
+
+if 'internal' in locals():
+    importlib.reload(internal)
+    importlib.reload(cad)
+    importlib.reload(toolpath)
+    importlib.reload(exports)
+
+
+from bpy.types import (
+        AddonPreferences,
+        )
+
+
+def UpdateDummy(object, context):
+    scene = context.scene
+    SINGLEDROP = scene.UTSingleDrop
+    MOREDROP = scene.UTMOREDROP
+    LOFTDROP = scene.UTLoftDrop
+    ADVANCEDDROP = scene.UTAdvancedDrop
+    EXTENDEDDROP = scene.UTExtendedDrop
+    UTILSDROP = scene.UTUtilsDrop
+
+
+class SeparateOutline(Operator):
+    bl_idname = "object.sep_outline"
+    bl_label = "Separate Outline"
+    bl_options = {'REGISTER', 'UNDO'}
+    bl_description = "Makes 'Outline' separate mesh"
+
+    @classmethod
+    def poll(cls, context):
+        return (context.object is not None and
+                context.object.type == 'CURVE')
+
+    def execute(self, context):
+        bpy.ops.object.mode_set(mode = 'EDIT')
+        bpy.ops.curve.separate()
+
+        return {'FINISHED'}
+
+
+class CurveTools2Settings(PropertyGroup):
+    # selection
+    SelectedObjects: CollectionProperty(
+            type=Properties.CurveTools2SelectedObject
+            )
+    NrSelectedObjects: IntProperty(
+            name="NrSelectedObjects",
+            default=0,
+            description="Number of selected objects",
+            update=UpdateDummy
+            )
+    # curve
+    CurveLength: FloatProperty(
+            name="CurveLength",
+            default=0.0,
+            precision=6
+            )
+    # splines
+    SplineResolution: IntProperty(
+            name="SplineResolution",
+            default=64,
+            min=2, max=1024,
+            soft_min=2,
+            description="Spline resolution will be set to this value"
+            )
+    SplineRemoveLength: FloatProperty(
+            name="SplineRemoveLength",
+            default=0.001,
+            precision=6,
+            description="Splines shorter than this threshold length will be removed"
+            )
+    SplineJoinDistance: FloatProperty(
+            name="SplineJoinDistance",
+            default=0.001,
+            precision=6,
+            description="Splines with starting/ending points closer to each other "
+                        "than this threshold distance will be joined"
+            )
+    SplineJoinStartEnd: BoolProperty(
+            name="SplineJoinStartEnd",
+            default=False,
+            description="Only join splines at the starting point of one and the ending point of the other"
+            )
+    splineJoinModeItems = (
+            ('At midpoint', 'At midpoint', 'Join splines at midpoint of neighbouring points'),
+            ('Insert segment', 'Insert segment', 'Insert segment between neighbouring points')
+            )
+    SplineJoinMode: EnumProperty(
+            items=splineJoinModeItems,
+            name="SplineJoinMode",
+            default='At midpoint',
+            description="Determines how the splines will be joined"
+            )
+    # curve intersection
+    LimitDistance: FloatProperty(
+            name="LimitDistance",
+            default=0.0001,
+            precision=6,
+            description="Displays the result of the curve length calculation"
+            )
+
+    intAlgorithmItems = (
+            ('3D', '3D', 'Detect where curves intersect in 3D'),
+            ('From View', 'From View', 'Detect where curves intersect in the RegionView3D')
+            )
+    IntersectCurvesAlgorithm: EnumProperty(
+            items=intAlgorithmItems,
+            name="IntersectCurvesAlgorithm",
+            description="Determines how the intersection points will be detected",
+            default='3D'
+            )
+    intModeItems = (
+            ('Insert', 'Insert', 'Insert points into the existing spline(s)'),
+            ('Split', 'Split', 'Split the existing spline(s) into 2'),
+            ('Empty', 'Empty', 'Add empty at intersections')
+            )
+    IntersectCurvesMode: EnumProperty(
+            items=intModeItems,
+            name="IntersectCurvesMode",
+            description="Determines what happens at the intersection points",
+            default='Split'
+            )
+    intAffectItems = (
+            ('Both', 'Both', 'Insert points into both curves'),
+            ('Active', 'Active', 'Insert points into active curve only'),
+            ('Other', 'Other', 'Insert points into other curve only')
+            )
+    IntersectCurvesAffect: EnumProperty(
+            items=intAffectItems,
+            name="IntersectCurvesAffect",
+            description="Determines which of the selected curves will be affected by the operation",
+            default='Both'
+            )
+    PathFinderRadius: FloatProperty(
+            name="PathFinder detection radius",
+            default=0.2,
+            precision=6,
+            description="PathFinder detection radius"
+            )
+    curve_vertcolor: FloatVectorProperty(
+            name="OUT",
+            default=(0.2, 0.9, 0.9, 1),
+            size=4,
+            subtype="COLOR",
+            min=0,
+            max=1
+            )
+    path_color: FloatVectorProperty(
+            name="OUT",
+            default=(0.2, 0.9, 0.9, 0.1),
+            size=4,
+            subtype="COLOR",
+            min=0,
+            max=1
+            )
+    path_thickness: IntProperty(
+            name="Path thickness",
+            default=10,
+            min=1, max=1024,
+            soft_min=2,
+            description="Path thickness (px)"
+            )
+
+
+class VIEW3D_PT_CurvePanel(Panel):
+    bl_label = "Curve Tools 2"
+    bl_space_type = "VIEW_3D"
+    bl_region_type = "UI"
+    bl_options = {'DEFAULT_CLOSED'}
+    bl_category = "Tools"
+
+    @classmethod
+    def poll(cls, context):
+        return context.scene is not None
+
+    def draw(self, context):
+        scene = context.scene
+        SINGLEDROP = scene.UTSingleDrop
+        MOREDROP = scene.UTMOREDROP
+        LOFTDROP = scene.UTLoftDrop
+        ADVANCEDDROP = scene.UTAdvancedDrop
+        EXTENDEDDROP = scene.UTExtendedDrop
+        UTILSDROP = scene.UTUtilsDrop
+        layout = self.layout
+
+        # Single Curve options
+        box1 = self.layout.box()
+        col = box1.column(align=True)
+        row = col.row(align=True)
+        row.prop(scene, "UTSingleDrop", icon="TRIA_DOWN")
+        if SINGLEDROP:
+            # A. 1 curve
+            row = col.row(align=True)
+
+            # A.1 curve info/length
+            row.operator("curvetools2.operatorcurveinfo", text="Curve info")
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorcurvelength", text="Calc Length")
+            row.prop(context.scene.curvetools, "CurveLength", text="")
+
+            # A.2 splines info
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorsplinesinfo", text="Curve splines info")
+
+            # A.3 segments info
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorsegmentsinfo", text="Curve segments info")
+
+            # A.4 origin to spline0start
+            row = col.row(align=True)
+            row.operator("curvetools2.operatororigintospline0start", text="Set origin to spline start")
+
+        # Double Curve options
+        box2 = self.layout.box()
+        col = box2.column(align=True)
+        row = col.row(align=True)
+        row.prop(scene, "UTMOREDROP", icon="TRIA_DOWN")
+
+        if MOREDROP:
+            # B. 2 curves
+            row = col.row(align=True)
+
+            # B.1 curve intersections
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorintersectcurves", text="Intersect curves")
+
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "LimitDistance", text="LimitDistance")
+            # row.active = (context.scene.curvetools.IntersectCurvesAlgorithm == '3D')
+
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "IntersectCurvesAlgorithm", text="Algorithm")
+
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "IntersectCurvesMode", text="Mode")
+
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "IntersectCurvesAffect", text="Affect")
+
+        # Loft options
+        box1 = self.layout.box()
+        col = box1.column(align=True)
+        row = col.row(align=True)
+        row.prop(scene, "UTLoftDrop", icon="TRIA_DOWN")
+
+        if LOFTDROP:
+            # B.2 surface generation
+            wm = context.window_manager
+            scene = context.scene
+            layout = self.layout
+            layout.operator("curvetools2.create_auto_loft")
+            lofters = [o for o in scene.objects if "autoloft" in o.keys()]
+            for o in lofters:
+                layout.label(text=o.name)
+            # layout.prop(o, '["autoloft"]', toggle=True)
+            layout.prop(wm, "auto_loft", toggle=True)
+            layout.operator("curvetools2.update_auto_loft_curves")
+
+        # Advanced options
+        box1 = self.layout.box()
+        col = box1.column(align=True)
+        row = col.row(align=True)
+        row.prop(scene, "UTAdvancedDrop", icon="TRIA_DOWN")
+        if ADVANCEDDROP:
+            # C. 3 curves
+            row = col.row(align=True)
+            row.operator("object._curve_outline", text="Curve Outline")
+            row = col.row(align=True)
+            row.operator("object.sep_outline", text="Separate Outline")
+            row = col.row(align=True)
+            row.operator("curve.remove_doubles", text="Remove Doubles")
+            row = col.row(align=True)
+            row.operator("curve.bezier_points_fillet", text='Fillet')
+            row = col.row(align=True)
+            row.operator("curve.bezier_spline_divide", text='Divide')
+            row = col.row(align=True)
+            row.operator("curve.scale_reset", text='Scale Reset')
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorbirail", text="Birail")
+            row = col.row(align=True)        
+            row.operator("curvetools2.convert_selected_face_to_bezier", text="Convert selected faces to Bezier")
+            row = col.row(align=True)
+            row.operator("curvetools2.convert_bezier_to_surface", text="Convert Bezier to Surface")
+        
+        # Extended options
+        box1 = self.layout.box()
+        col = box1.column(align=True)
+        row = col.row(align=True)
+        row.prop(scene, "UTExtendedDrop", icon="TRIA_DOWN")
+        if EXTENDEDDROP:
+            for operator in cad.operators:
+                row = col.row(align=True)
+                row.operator(operator.bl_idname)
+            
+            for operator in toolpath.operators:
+                row = col.row(align=True)
+                row.operator(operator.bl_idname)
+
+        # Utils Curve options
+        box1 = self.layout.box()
+        col = box1.column(align=True)
+        row = col.row(align=True)
+        row.prop(scene, "UTUtilsDrop", icon="TRIA_DOWN")
+        if UTILSDROP:
+            # D.1 set spline resolution
+            row = col.row(align=True)
+            row.label(text="Show point Resolution:")
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorsplinessetresolution", text="Set resolution")
+            row.prop(context.scene.curvetools, "SplineResolution", text="")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "curve_vertcolor", text="")
+            row = col.row(align=True)
+            row.operator("curve.show_resolution", text="Run [ESC]")
+
+            # D.2 remove splines
+            row = col.row(align=True)
+            row.label(text="Remove splines:")
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorsplinesremovezerosegment", text="Remove 0-segments splines")
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorsplinesremoveshort", text="Remove short splines")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "SplineRemoveLength", text="Threshold remove")
+
+            # D.3 join splines
+            row = col.row(align=True)
+            row.label(text="Join splines:")
+            row = col.row(align=True)
+            row.operator("curvetools2.operatorsplinesjoinneighbouring", text="Join neighbouring splines")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "SplineJoinDistance", text="Threshold join")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "SplineJoinStartEnd", text="Only at start & end")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "SplineJoinMode", text="Join mode")
+
+            row = col.row(align=True)
+            row.label(text="PathFinder:")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "PathFinderRadius", text="PathFinder Radius")
+            row = col.row(align=True)
+            row.prop(context.scene.curvetools, "path_color", text="")
+            row.prop(context.scene.curvetools, "path_thickness", text="")
+            row = col.row(align=True)
+            row.operator("curvetools2.pathfinder", text="Run Path Finder [ESC]")
+            row = col.row(align=True)
+            row.label(text="ESC or TAB - exit from PathFinder")
+            row = col.row(align=True)
+            row.label(text="X or DEL - delete")
+            row = col.row(align=True)
+            row.label(text="Alt + mouse click - select spline")
+            row = col.row(align=True)
+            row.label(text="A - deselect all")
+            
+
+# Add-ons Preferences Update Panel
+
+# Define Panel classes for updating
+panels = (
+        VIEW3D_PT_CurvePanel,
+        )
+
+
+def update_panel(self, context):
+    message = "Curve Tools 2: Updating Panel locations has failed"
+    try:
+        for panel in panels:
+            if "bl_rna" in panel.__dict__:
+                bpy.utils.unregister_class(panel)
+
+        for panel in panels:
+            panel.bl_category = context.preferences.addons[__name__].preferences.category
+            bpy.utils.register_class(panel)
+
+    except Exception as e:
+        print("\n[{}]\n{}\n\nError:\n{}".format(__name__, message, e))
+        pass
+
+
+class CurveAddonPreferences(AddonPreferences):
+    # this must match the addon name, use '__package__'
+    # when defining this in a submodule of a python package.
+    bl_idname = __name__
+
+    category: StringProperty(
+            name="Tab Category",
+            description="Choose a name for the category of the panel",
+            default="Tools",
+            update=update_panel
+            )
+
+    def draw(self, context):
+        layout = self.layout
+
+        row = layout.row()
+        col = row.column()
+        col.label(text="Tab Category:")
+        col.prop(self, "category", text="")
+
+def menu_file_export(self, context):
+    for operator in exports.operators:
+        self.layout.operator(operator.bl_idname)
+
+def menu_file_import(self, context):
+    for operator in imports.operators:
+        self.layout.operator(operator.bl_idname)
+
+# REGISTER
+classes = cad.operators + toolpath.operators + exports.operators + [
+    Properties.CurveTools2SelectedObject,
+    CurveAddonPreferences,
+    CurveTools2Settings,
+    Operators.OperatorCurveInfo,
+    Operators.OperatorCurveLength,
+    Operators.OperatorSplinesInfo,
+    Operators.OperatorSegmentsInfo,
+    Operators.OperatorOriginToSpline0Start,
+    Operators.OperatorIntersectCurves,
+    Operators.OperatorLoftCurves,
+    Operators.OperatorSweepCurves,
+    Operators.OperatorBirail,
+    Operators.OperatorSplinesSetResolution,
+    Operators.OperatorSplinesRemoveZeroSegment,
+    Operators.OperatorSplinesRemoveShort,
+    Operators.OperatorSplinesJoinNeighbouring,
+    VIEW3D_PT_CurvePanel,
+    SeparateOutline,
+    Operators.ConvertSelectedFacesToBezier,
+    Operators.ConvertBezierToSurface,
+    PathFinder.PathFinder,
+    ShowCurveResolution.ShowCurveResolution,
+    ]
+
+def register():
+    bpy.types.Scene.UTSingleDrop = BoolProperty(
+            name="One Curve",
+            default=False,
+            description="One Curve"
+            )
+    bpy.types.Scene.UTMOREDROP = BoolProperty(
+            name="Curves",
+            default=False,
+            description="Curves"
+            )
+    bpy.types.Scene.UTLoftDrop = BoolProperty(
+            name="Two Curves Loft",
+            default=False,
+            description="Two Curves Loft"
+            )
+    bpy.types.Scene.UTAdvancedDrop = BoolProperty(
+            name="Advanced",
+            default=True,
+            description="Advanced"
+            )
+    bpy.types.Scene.UTExtendedDrop = BoolProperty(
+            name="Extended",
+            default=False,
+            description="Extended"
+            )
+    bpy.types.Scene.UTUtilsDrop = BoolProperty(
+            name="Curves Utils",
+            default=True,
+            description="Curves Utils"
+            )
+    
+    for cls in classes:
+        bpy.utils.register_class(cls)
+    
+    auto_loft.register()
+    
+    curve_outline.register()
+    
+    curve_remove_doubles.register()
+    
+    bpy.types.TOPBAR_MT_file_export.append(menu_file_export)
+    
+    bpy.types.Scene.curvetools = bpy.props.PointerProperty(type=CurveTools2Settings)
+
+
+def unregister():
+    del bpy.types.Scene.UTSingleDrop
+    del bpy.types.Scene.UTMOREDROP
+    del bpy.types.Scene.UTLoftDrop
+    del bpy.types.Scene.UTAdvancedDrop
+    del bpy.types.Scene.UTExtendedDrop
+    del bpy.types.Scene.UTUtilsDrop
+    
+    auto_loft.unregister()
+    
+    curve_outline.unregister()
+    
+    curve_remove_doubles.unregister()
+    
+    bpy.types.TOPBAR_MT_file_export.remove(menu_file_export)
+    
+    for cls in classes:
+        bpy.utils.unregister_class(cls)
+
+
+if __name__ == "__main__":
+    register()
diff --git a/curve_tools/auto_loft.py b/curve_tools/auto_loft.py
new file mode 100644
index 00000000..8cb89245
--- /dev/null
+++ b/curve_tools/auto_loft.py
@@ -0,0 +1,85 @@
+import bpy
+from bpy.props import BoolProperty
+from bpy.types import Operator, Panel
+from curve_tools.Surfaces import LoftedSurface
+from curve_tools.Curves import Curve
+from curve_tools import Util
+
+
+class OperatorAutoLoftCurves(Operator):
+    bl_idname = "curvetools2.create_auto_loft"
+    bl_label = "Loft"
+    bl_description = "Lofts selected curves"
+
+    @classmethod
+    def poll(cls, context):
+        return Util.Selected2Curves()
+
+    def execute(self, context):
+        #print("### TODO: OperatorLoftCurves.execute()")
+        mesh = bpy.data.meshes.new("LoftMesh")
+
+        curve0 = context.selected_objects[0]
+        curve1 = context.selected_objects[1]
+
+        ls = LoftedSurface(Curve(curve0), Curve(curve1), "AutoLoft")
+
+        ls.bMesh.to_mesh(mesh)
+
+        loftobj = bpy.data.objects.new(self.name, mesh)
+
+        context.collection.objects.link(loftobj)
+        loftobj["autoloft"] = True
+        if loftobj.get('_RNA_UI') is None:
+            loftobj['_RNA_UI'] = {}
+        loftobj['_RNA_UI']["autoloft"] = {
+                           "name": "Auto Loft",
+                           "description": "Auto loft from %s to %s" % (curve0.name, curve1.name),
+                           "curve0": curve0.name,
+                           "curve1": curve1.name}
+        #print(loftobj['_RNA_UI'].to_dict())
+        #self.report({'INFO'}, "OperatorAutoLoftCurves.execute()")
+
+        return {'FINISHED'}
+
+class AutoLoftModalOperator(Operator):
+    """Auto Loft"""
+    bl_idname = "curvetools2.update_auto_loft_curves"
+    bl_label = "Update Auto Loft"
+    bl_description = "Update Lofts selected curves"
+
+    _timer = None
+    @classmethod
+    def poll(cls, context):
+        # two curves selected.
+        return True
+
+    def execute(self, context):
+        scene = context.scene
+        lofters = [o for o in scene.objects if "autoloft" in o.keys()]
+        # quick hack
+        #print("TIMER", lofters)
+
+        for loftmesh in lofters:
+            #loftmesh.hide_select = True
+            rna = loftmesh['_RNA_UI']["autoloft"].to_dict()
+            curve0 = scene.objects.get(rna["curve0"])
+            curve1 = scene.objects.get(rna["curve1"])
+            if curve0 and curve1:
+                ls = LoftedSurface(Curve(curve0), Curve(curve1), loftmesh.name)
+                ls.bMesh.to_mesh(loftmesh.data)
+        return {'FINISHED'}
+
+def register():
+    bpy.utils.register_class(AutoLoftModalOperator)
+    bpy.utils.register_class(OperatorAutoLoftCurves)
+    bpy.types.WindowManager.auto_loft = BoolProperty(default=False,
+                                                     name="Auto Loft")
+    bpy.context.window_manager.auto_loft = False
+
+def unregister():
+    bpy.utils.unregister_class(AutoLoftModalOperator)
+    bpy.utils.unregister_class(OperatorAutoLoftCurves)
+
+if __name__ == "__main__":
+    register()
diff --git a/curve_tools/cad.py b/curve_tools/cad.py
new file mode 100644
index 00000000..e518875e
--- /dev/null
+++ b/curve_tools/cad.py
@@ -0,0 +1,248 @@
+#  ***** GPL LICENSE BLOCK *****
+#
+#  This program is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#  All rights reserved.
+#
+#  ***** GPL LICENSE BLOCK *****
+
+bl_info = {
+    'name': 'Curve CAD Tools',
+    'author': 'Alexander Meißner',
+    'version': (1, 0, 0),
+    'blender': (2, 80, 0),
+    'category': 'Curve',
+    'wiki_url': 'https://github.com/Lichtso/curve_cad',
+    'tracker_url': 'https://github.com/lichtso/curve_cad/issues'
+}
+
+import bpy
+from . import internal
+
+class Fillet(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_fillet'
+    bl_description = bl_label = 'Fillet'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    radius: bpy.props.FloatProperty(name='Radius', description='Radius of the rounded corners', unit='LENGTH', min=0.0, default=0.1)
+    chamfer_mode: bpy.props.BoolProperty(name='Chamfer', description='Cut off sharp without rounding', default=False)
+    limit_half_way: bpy.props.BoolProperty(name='Limit Half Way', description='Limits the segements to half their length in order to prevent collisions', default=False)
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        splines = internal.getSelectedSplines(True, True, True)
+        if len(splines) == 0:
+            self.report({'WARNING'}, 'Nothing selected')
+            return {'CANCELLED'}
+        for spline in splines:
+            internal.filletSpline(spline, self.radius, self.chamfer_mode, self.limit_half_way)
+            bpy.context.object.data.splines.remove(spline)
+        return {'FINISHED'}
+
+class Boolean(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_boolean'
+    bl_description = bl_label = 'Boolean'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    operation: bpy.props.EnumProperty(name='Type', items=[
+        ('UNION', 'Union', 'Boolean OR', 0),
+        ('INTERSECTION', 'Intersection', 'Boolean AND', 1),
+        ('DIFFERENCE', 'Difference', 'Active minus Selected', 2)
+    ])
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        if bpy.context.object.data.dimensions != '2D':
+            self.report({'WARNING'}, 'Can only be applied in 2D')
+            return {'CANCELLED'}
+        splines = internal.getSelectedSplines(True, True)
+        if len(splines) != 2:
+            self.report({'WARNING'}, 'Invalid selection')
+            return {'CANCELLED'}
+        bpy.ops.curve.spline_type_set(type='BEZIER')
+        splineA = bpy.context.object.data.splines.active
+        splineB = splines[0] if (splines[1] == splineA) else splines[1]
+        if not internal.bezierBooleanGeometry(splineA, splineB, self.operation):
+            self.report({'WARNING'}, 'Invalid selection')
+            return {'CANCELLED'}
+        return {'FINISHED'}
+
+class Intersection(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_intersection'
+    bl_description = bl_label = 'Intersection'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        segments = internal.bezierSegments(bpy.context.object.data.splines, True)
+        if len(segments) < 2:
+            self.report({'WARNING'}, 'Invalid selection')
+            return {'CANCELLED'}
+
+        internal.bezierMultiIntersection(segments)
+        return {'FINISHED'}
+
+class MergeEnds(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_merge_ends'
+    bl_description = bl_label = 'Merge Ends'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        points = []
+        selected_splines = []
+        is_last_point = []
+        for spline in bpy.context.object.data.splines:
+            if spline.type != 'BEZIER' or spline.use_cyclic_u:
+                continue
+            if spline.bezier_points[0].select_control_point:
+                points.append(spline.bezier_points[0])
+                selected_splines.append(spline)
+                is_last_point.append(False)
+            if spline.bezier_points[-1].select_control_point:
+                points.append(spline.bezier_points[-1])
+                selected_splines.append(spline)
+                is_last_point.append(True)
+
+        if len(points) != 2:
+            self.report({'WARNING'}, 'Invalid selection')
+            return {'CANCELLED'}
+
+        points[0].handle_left_type = 'FREE'
+        points[0].handle_right_type = 'FREE'
+        new_co = (points[0].co+points[1].co)*0.5
+
+        handle = (points[1].handle_left if is_last_point[1] else points[1].handle_right)+new_co-points[1].co
+        if is_last_point[0]:
+            points[0].handle_left += new_co-points[0].co
+            points[0].handle_right = handle
+        else:
+            points[0].handle_right += new_co-points[0].co
+            points[0].handle_left = handle
+        points[0].co = new_co
+
+        bpy.ops.curve.select_all(action='DESELECT')
+        points[1].select_control_point = True
+        bpy.ops.curve.delete()
+        selected_splines[0].bezier_points[-1 if is_last_point[0] else 0].select_control_point = True
+        selected_splines[1].bezier_points[-1 if is_last_point[1] else 0].select_control_point = True
+        bpy.ops.curve.make_segment()
+        bpy.ops.curve.select_all(action='DESELECT')
+        return {'FINISHED'}
+
+class Subdivide(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_subdivide'
+    bl_description = bl_label = 'Subdivide'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    params: bpy.props.StringProperty(name='Params', default='0.25 0.5 0.75')
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        segments = internal.bezierSegments(bpy.context.object.data.splines, True)
+        if len(segments) == 0:
+            self.report({'WARNING'}, 'Nothing selected')
+            return {'CANCELLED'}
+
+        cuts = []
+        for param in self.params.split(' '):
+            cuts.append({'param': max(0.0, min(float(param), 1.0))})
+        cuts.sort(key=(lambda cut: cut['param']))
+        for segment in segments:
+            segment['cuts'].extend(cuts)
+        internal.subdivideBezierSegments(segments)
+        return {'FINISHED'}
+
+class Array(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_array'
+    bl_description = bl_label = 'Array'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    offset: bpy.props.FloatVectorProperty(name='Offset', unit='LENGTH', description='Vector between to copies', subtype='DIRECTION', default=(0.0, 0.0, -1.0), size=3)
+    count: bpy.props.IntProperty(name='Count', description='Number of copies', min=1, default=2)
+    connect: bpy.props.BoolProperty(name='Connect', description='Concatenate individual copies', default=False)
+    serpentine: bpy.props.BoolProperty(name='Serpentine', description='Switch direction of every second copy', default=False)
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        splines = internal.getSelectedSplines(True, True)
+        if len(splines) == 0:
+            self.report({'WARNING'}, 'Nothing selected')
+            return {'CANCELLED'}
+        internal.arrayModifier(splines, self.offset, self.count, self.connect, self.serpentine)
+        return {'FINISHED'}
+
+class Circle(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_circle'
+    bl_description = bl_label = 'Circle'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        segments = internal.bezierSegments(bpy.context.object.data.splines, True)
+        if len(segments) != 1:
+            self.report({'WARNING'}, 'Invalid selection')
+            return {'CANCELLED'}
+
+        segment = internal.bezierSegmentPoints(segments[0]['beginPoint'], segments[0]['endPoint'])
+        circle = internal.circleOfBezier(segment)
+        if circle == None:
+            self.report({'WARNING'}, 'Not a circle')
+            return {'CANCELLED'}
+        bpy.context.scene.cursor.location = circle.center
+        bpy.context.scene.cursor.rotation_mode = 'QUATERNION'
+        bpy.context.scene.cursor.rotation_quaternion = circle.orientation.to_quaternion()
+        return {'FINISHED'}
+
+class Length(bpy.types.Operator):
+    bl_idname = 'curve.bezier_cad_length'
+    bl_description = bl_label = 'Length'
+
+    @classmethod
+    def poll(cls, context):
+        return internal.curveObject()
+
+    def execute(self, context):
+        segments = internal.bezierSegments(bpy.context.object.data.splines, True)
+        if len(segments) == 0:
+            self.report({'WARNING'}, 'Nothing selected')
+            return {'CANCELLED'}
+
+        length = 0
+        for segment in segments:
+            length += internal.bezierLength(internal.bezierSegmentPoints(segment['beginPoint'], segment['endPoint']))
+        self.report({'INFO'}, bpy.utils.units.to_string(bpy.context.scene.unit_settings.system, 'LENGTH', length))
+        return {'FINISHED'}
+
+operators = [Fillet, Boolean, Intersection, MergeEnds, Subdivide, Array, Circle, Length]
diff --git a/curve_tools/curve_outline.py b/curve_tools/curve_outline.py
new file mode 100644
index 00000000..6847e3c5
--- /dev/null
+++ b/curve_tools/curve_outline.py
@@ -0,0 +1,112 @@
+'''
+by Yann Bertrand, january 2014.
+
+BEGIN GPL LICENSE BLOCK
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software Foundation,
+Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+END GPL LICENCE BLOCK
+'''
+
+bl_info = {
+    "name": "Curve Outline",
+    "description": "creates an Outline",
+    "category": "Object",
+    "author": "Yann Bertrand (jimflim)",
+    "version": (0, 4),
+    "blender": (2, 69, 0),
+}
+
+import bpy
+from mathutils import Vector
+from mathutils.geometry import intersect_line_line
+
+
+def createOutline(curve, outline):
+
+    for spline in curve.data.splines[:]:
+        p = spline.bezier_points
+        out = []
+
+        n = ((p[0].handle_right-p[0].co).normalized()-(p[0].handle_left-p[0].co).normalized()).normalized()
+        n = Vector((-n[1], n[0], n[2]))
+        o = p[0].co+outline*n
+        out.append(o)
+
+        for i in range(1,len(p)):
+            n = ((p[i].handle_right-p[i].co).normalized()-(p[i].handle_left-p[i].co).normalized()).normalized()
+            n = Vector((-n[1], n[0], n[2]))
+            o = intersect_line_line(out[-1], (out[-1]+p[i].co-p[i-1].co), p[i].co, p[i].co+n)[0]
+            out.append(o)
+
+        curve.data.splines.new('BEZIER')
+        if spline.use_cyclic_u:
+            curve.data.splines[-1].use_cyclic_u = True
+        p_out = curve.data.splines[-1].bezier_points
+        p_out.add(len(out)-1)
+
+        for i in range(len(out)):
+            p_out[i].handle_left_type = 'FREE'
+            p_out[i].handle_right_type = 'FREE'
+
+            p_out[i].co = out[i]
+
+            if i<len(out)-1:
+                l = (p[i+1].co-p[i].co).length
+                l2 = (out[i]-out[i+1]).length
+
+            if i==0:
+                p_out[i].handle_left = out[i] + ((p[i].handle_left-p[i].co)*l2/l)
+            if i<len(out)-1:
+                p_out[i+1].handle_left = out[i+1] + ((p[i+1].handle_left-p[i+1].co)*l2/l)
+            p_out[i].handle_right = out[i] + ((p[i].handle_right-p[i].co)*l2/l)
+
+        for i in range(len(p)):
+            p_out[i].handle_left_type = p[i].handle_left_type
+            p_out[i].handle_right_type = p[i].handle_right_type
+
+    return
+
+
+class CurveOutline(bpy.types.Operator):
+    """Curve Outliner"""
+    bl_idname = "object._curve_outline"
+    bl_label = "Create Outline"
+    bl_options = {'REGISTER', 'UNDO'}
+    outline: bpy.props.FloatProperty(name="Amount", default=0.1)
+
+    @classmethod
+    def poll(cls, context):
+        return (context.object is not None and
+                context.object.type == 'CURVE')
+
+    def execute(self, context):
+        createOutline(context.object, self.outline)
+        return {'FINISHED'}
+
+    def invoke(self, context, event):
+        return context.window_manager.invoke_props_popup(self, event)
+
+def menu_func(self, context):
+    self.layout.operator(CurveOutline.bl_idname)
+
+def register():
+    bpy.utils.register_class(CurveOutline)
+
+def unregister():
+    bpy.utils.unregister_class(CurveOutline)
+
+if __name__ == "__main__":
+    register()
diff --git a/curve_tools/curve_remove_doubles.py b/curve_tools/curve_remove_doubles.py
new file mode 100644
index 00000000..ca26895b
--- /dev/null
+++ b/curve_tools/curve_remove_doubles.py
@@ -0,0 +1,110 @@
+import bpy, mathutils
+
+
+bl_info = {
+    'name': 'Curve Remove Doubles',
+    'author': 'Michael Soluyanov',
+    'version': (1, 1),
+    'blender': (2, 80, 0),
+    'location': 'View3D > Context menu (W/RMB) > Remove Doubles',
+    'description': 'Adds comand "Remove Doubles" for curves',
+    'category': 'Object'
+}
+
+def main(context, distance = 0.01):
+    
+    obj = context.active_object
+    dellist = []
+    
+    if bpy.ops.object.mode_set.poll():
+        bpy.ops.object.mode_set(mode='EDIT')
+    
+    for spline in obj.data.splines: 
+        if len(spline.bezier_points) > 1:
+            for i in range(0, len(spline.bezier_points)): 
+                
+                if i == 0:
+                    ii = len(spline.bezier_points) - 1
+                else:        
+                    ii = i - 1
+                    
+                dot = spline.bezier_points[i];
+                dot1 = spline.bezier_points[ii];   
+                    
+                while dot1 in dellist and i != ii:
+                    ii -= 1
+                    if ii < 0: 
+                        ii = len(spline.bezier_points)-1
+                    dot1 = spline.bezier_points[ii]
+                    
+                if dot.select_control_point and dot1.select_control_point and (i!=0 or spline.use_cyclic_u):   
+                    
+                    if (dot.co-dot1.co).length < distance:
+                        # remove points and recreate hangles
+                        dot1.handle_right_type = "FREE"
+                        dot1.handle_right = dot.handle_right
+                        dot1.co = (dot.co + dot1.co) / 2
+                        dellist.append(dot)
+                        
+                    else:
+                        # Handles that are on main point position converts to vector,
+                        # if next handle are also vector
+                        if dot.handle_left_type == 'VECTOR' and (dot1.handle_right - dot1.co).length < distance:
+                            dot1.handle_right_type = "VECTOR"
+                        if dot1.handle_right_type == 'VECTOR' and (dot.handle_left - dot.co).length < distance:
+                            dot.handle_left_type = "VECTOR"  
+                      
+                            
+    
+    bpy.ops.curve.select_all(action = 'DESELECT')
+
+    for dot in dellist:
+        dot.select_control_point = True
+        
+    count = len(dellist)
+    
+    bpy.ops.curve.delete(type = 'VERT')
+    
+    bpy.ops.curve.select_all(action = 'SELECT')
+    
+    return count
+    
+
+
+class CurveRemvDbs(bpy.types.Operator):
+    """Merge consecutive points that are near to each other"""
+    bl_idname = 'curve.remove_doubles'
+    bl_label = 'Remove Doubles'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    distance: bpy.props.FloatProperty(name = 'Distance', default = 0.01)
+
+    @classmethod
+    def poll(cls, context):
+        obj = context.active_object
+        return (obj and obj.type == 'CURVE')
+
+    def execute(self, context):
+        removed=main(context, self.distance)
+        self.report({'INFO'}, "Removed %d bezier points" % removed)
+        return {'FINISHED'}
+
+
+
+
+def menu_func(self, context):
+    self.layout.operator(CurveRemvDbs.bl_idname, text='Remove Doubles')
+
+def register():
+    bpy.utils.register_class(CurveRemvDbs)
+    bpy.types.VIEW3D_MT_edit_curve_context_menu.append(menu_func)
+
+def unregister():
+    bpy.utils.unregister_class(CurveRemvDbs)
+    bpy.types.VIEW3D_MT_edit_curve_context_menu.remove(menu_func)
+
+if __name__ == "__main__":
+    register()
+
+
+
diff --git a/curve_tools/exports.py b/curve_tools/exports.py
new file mode 100644
index 00000000..b3b7d900
--- /dev/null
+++ b/curve_tools/exports.py
@@ -0,0 +1,227 @@
+#  ***** GPL LICENSE BLOCK *****
+#
+#  This program is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#  All rights reserved.
+#
+#  ***** GPL LICENSE BLOCK *****
+
+import bpy, math
+from mathutils import Vector, Matrix
+from bpy_extras.io_utils import ExportHelper
+from . import internal
+
+class SvgExport(bpy.types.Operator, ExportHelper):
+    bl_idname = 'export_svg_format.svg'
+    bl_description = bl_label = 'Curves (.svg)'
+    filename_ext = '.svg'
+
+    selection_only: bpy.props.BoolProperty(name='Selection only', description='instead of exporting all visible curves')
+    absolute_coordinates: bpy.props.BoolProperty(name='Absolute coordinates', description='instead of relative coordinates')
+    viewport_projection: bpy.props.BoolProperty(name='Viewport projection', description='WYSIWYG instead of an local orthographic projection')
+    unit_name: bpy.props.EnumProperty(name='Unit', items=internal.units, default='mm')
+
+    def serialize_point(self, position, update_ref_position=True):
+        if self.transform:
+            position = self.transform@Vector((position[0], position[1], position[2], 1.0))
+            position *= 0.5/position.w
+        ref_position = self.origin if self.absolute_coordinates else self.ref_position
+        command = '{:.3f},{:.3f}'.format((position[0]-ref_position[0])*self.scale[0], (position[1]-ref_position[1])*self.scale[1])
+        if update_ref_position:
+            self.ref_position = position
+        return command
+
+    def serialize_point_command(self, point, drawing):
+        if self.absolute_coordinates:
+            return ('L' if drawing else 'M')+self.serialize_point(point.co)
+        else:
+            return ('l' if drawing else 'm')+self.serialize_point(point.co)
+
+    def serialize_curve_command(self, prev, next):
+        return ('C' if self.absolute_coordinates else 'c')+self.serialize_point(prev.handle_right, False)+' '+self.serialize_point(next.handle_left, False)+' '+self.serialize_point(next.co)
+
+    def serialize_spline(self, spline):
+        path = ''
+        points = spline.bezier_points if spline.type == 'BEZIER' else spline.points
+
+        for index, next in enumerate(points):
+            if index == 0:
+                path += self.serialize_point_command(next, False)
+            elif spline.type == 'BEZIER' and (points[index-1].handle_right_type != 'VECTOR' or next.handle_left_type != 'VECTOR'):
+                path += self.serialize_curve_command(points[index-1], next)
+            else:
+                path += self.serialize_point_command(next, True)
+
+        if spline.use_cyclic_u:
+            if spline.type == 'BEZIER' and (points[-1].handle_right_type != 'VECTOR' or points[0].handle_left_type != 'VECTOR'):
+                path += self.serialize_curve_command(points[-1], points[0])
+            else:
+                self.serialize_point(points[0].co)
+            path += 'Z' if self.absolute_coordinates else 'z'
+
+        return path
+
+    def serialize_object(self, obj):
+        if self.area:
+            self.transform = self.area.spaces.active.region_3d.perspective_matrix@obj.matrix_world
+            self.origin = Vector((-0.5, 0.5, 0, 0))
+        else:
+            self.transform = None
+            self.origin = Vector((obj.bound_box[0][0], obj.bound_box[7][1], obj.bound_box[0][2], 0))
+
+        xml = '\t<g id="'+obj.name+'">\n'
+        styles = {}
+        for spline in obj.data.splines:
+            style = 'none'
+            if obj.data.dimensions == '2D' and spline.use_cyclic_u:
+                if spline.material_index < len(obj.data.materials) and obj.data.materials[spline.material_index] != None:
+                    style = Vector(obj.data.materials[spline.material_index].diffuse_color)*255
+                else:
+                    style = Vector((0.8, 0.8, 0.8))*255
+                style = 'rgb({},{},{})'.format(round(style[0]), round(style[1]), round(style[2]))
+            if style in styles:
+                styles[style].append(spline)
+            else:
+                styles[style] = [spline]
+
+        for style, splines in styles.items():
+            style = 'fill:'+style+';'
+            if style == 'fill:none;':
+                style += 'stroke:black;'
+            xml += '\t\t<path style="'+style+'" d="'
+            self.ref_position = self.origin
+            for spline in splines:
+                xml += self.serialize_spline(spline)
+            xml += '"/>\n'
+
+        return xml+'\t</g>\n'
+
+    def execute(self, context):
+        objects = bpy.context.selected_objects if self.selection_only else bpy.context.visible_objects
+        curves = []
+        for obj in objects:
+            if obj.type == 'CURVE':
+                curves.append(obj)
+        if len(curves) == 0:
+            self.report({'WARNING'}, 'Nothing to export')
+            return {'CANCELLED'}
+
+        self.area = None
+        if self.viewport_projection:
+            for area in bpy.context.screen.areas:
+                if area.type == 'VIEW_3D':
+                    self.region = None
+                    for region in area.regions:
+                        if region.type == 'WINDOW':
+                            self.region = region
+                    if self.region == None:
+                        continue
+                    self.area = area
+                    self.bounds = Vector((self.region.width, self.region.height, 0))
+                    self.scale = Vector(self.bounds)
+                    if self.unit_name != 'px':
+                        self.unit_name = '-'
+
+        if self.area == None:
+            self.bounds = Vector((0, 0, 0))
+            for obj in curves:
+                self.bounds[0] = max(self.bounds[0], obj.bound_box[7][0]-obj.bound_box[0][0])
+                self.bounds[1] = max(self.bounds[1], obj.bound_box[7][1]-obj.bound_box[0][1])
+            self.scale = Vector((1, 1, 0))
+            for unit in internal.units:
+                if self.unit_name == unit[0]:
+                    self.scale *= 1.0/float(unit[2])
+                    break
+            self.scale *= context.scene.unit_settings.scale_length
+            self.bounds = Vector(a*b for a,b in zip(self.bounds, self.scale))
+
+        self.scale[1] *= -1
+        with open(self.filepath, 'w') as f:
+            svg_view = ('' if self.unit_name == '-' else 'width="{0:.3f}{2}" height="{1:.3f}{2}" ')+'viewBox="0 0 {0:.3f} {1:.3f}">\n'
+            f.write('''<?xml version="1.0" standalone="no"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<svg xmlns="http://www.w3.org/2000/svg" '''+svg_view.format(self.bounds[0], self.bounds[1], self.unit_name))
+            for obj in curves:
+                f.write(self.serialize_object(obj))
+            f.write('</svg>')
+
+        return {'FINISHED'}
+
+class GCodeExport(bpy.types.Operator, ExportHelper):
+    bl_idname = 'export_gcode_format.gcode'
+    bl_description = bl_label = 'Toolpath (.gcode)'
+    filename_ext = '.gcode'
+
+    speed: bpy.props.FloatProperty(name='Speed', description='Maximal speed in mm / minute', min=0, default=60)
+    step_angle: bpy.props.FloatProperty(name='Resolution', description='Smaller values make curves smoother by adding more vertices', unit='ROTATION', min=math.pi/128, default=math.pi/16)
+    local_coordinates: bpy.props.BoolProperty(name='Local coords', description='instead of global coordinates')
+    detect_circles: bpy.props.BoolProperty(name='Detect Circles', description='Export bezier circles and helixes as G02 and G03') # TODO: Detect polygon circles too, merge consecutive circle segments
+
+    @classmethod
+    def poll(cls, context):
+        obj = bpy.context.object
+        return obj != None and obj.type == 'CURVE' and len(obj.data.splines) == 1 and not obj.data.splines[0].use_cyclic_u
+
+    def execute(self, context):
+        self.scale = Vector((1, 1, 1))
+        self.scale *= context.scene.unit_settings.scale_length*1000.0
+        with open(self.filepath, 'w') as f:
+            f.write('G21\n') # Length is measured in millimeters
+            spline = bpy.context.object.data.splines[0]
+            if spline.use_cyclic_u:
+                return gcode
+            def transform(position):
+                result = Vector((position[0]*self.scale[0], position[1]*self.scale[1], position[2]*self.scale[2])) # , 1.0
+                return result if self.local_coordinates else bpy.context.object.matrix_world@result
+            points = spline.bezier_points if spline.type == 'BEZIER' else spline.points
+            prevSpeed = -1
+            for index, current in enumerate(points):
+                speed = self.speed*max(0.0, min(current.weight_softbody, 1.0))
+                if speed != prevSpeed and current.weight_softbody != 1.0:
+                    f.write('F{:.3f}\n'.format(speed))
+                    prevSpeed = speed
+                speed_code = 'G00' if current.weight_softbody == 1.0 else 'G01'
+                prev = points[index-1]
+                linear = spline.type != 'BEZIER' or index == 0 or (prev.handle_right_type == 'VECTOR' and current.handle_left_type == 'VECTOR')
+                position = transform(current.co)
+                if linear:
+                    f.write(speed_code+' X{:.3f} Y{:.3f} Z{:.3f}\n'.format(position[0], position[1], position[2]))
+                else:
+                    segment_points = internal.bezierSegmentPoints(prev, current)
+                    circle = None
+                    if self.detect_circles:
+                        for axis in range(0, 3):
+                            projected_points = []
+                            for point in segment_points:
+                                projected_point = Vector(point)
+                                projected_point[axis] = 0.0
+                                projected_points.append(projected_point)
+                            circle = internal.circleOfBezier(projected_points)
+                            if circle:
+                                normal = circle.orientation.col[2]
+                                center = transform(circle.center-prev.co)
+                                f.write('G{} G0{} I{:.3f} J{:.3f} K{:.3f} X{:.3f} Y{:.3f} Z{:.3f}\n'.format(19-axis, 3 if normal[axis] > 0.0 else 2, center[0], center[1], center[2], position[0], position[1], position[2]))
+                                break
+                    if circle == None:
+                        bezier_samples = 128
+                        prev_tangent = internal.bezierTangentAt(segment_points, 0).normalized()
+                        for t in range(1, bezier_samples+1):
+                            t /= bezier_samples
+                            tangent = internal.bezierTangentAt(segment_points, t).normalized()
+                            if t == 1 or math.acos(min(max(-1, prev_tangent@tangent), 1)) >= self.step_angle:
+                                position = transform(internal.bezierPointAt(segment_points, t))
+                                prev_tangent = tangent
+                                f.write(speed_code+' X{:.3f} Y{:.3f} Z{:.3f}\n'.format(position[0], position[1], position[2]))
+        return {'FINISHED'}
+
+operators = [SvgExport, GCodeExport]
diff --git a/curve_tools/internal.py b/curve_tools/internal.py
new file mode 100644
index 00000000..6741bb22
--- /dev/null
+++ b/curve_tools/internal.py
@@ -0,0 +1,958 @@
+#  ***** GPL LICENSE BLOCK *****
+#
+#  This program is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#  All rights reserved.
+#
+#  ***** GPL LICENSE BLOCK *****
+
+import bpy, math, cmath
+from mathutils import Vector, Matrix
+from collections import namedtuple
+
+units = [
+    ('-', 'None', '1.0', 0),
+    ('px', 'Pixel', '1.0', 1),
+    ('m', 'Meter', '1.0', 2),
+    ('dm', 'Decimeter', '0.1', 3),
+    ('cm', 'Centimeter', '0.01', 4),
+    ('mm', 'Millimeter', '0.001', 5),
+    ('yd', 'Yard', '0.9144', 6),
+    ('ft', 'Foot', '0.3048', 7),
+    ('in', 'Inch', '0.0254', 8)
+]
+
+param_tollerance = 0.0001
+AABB = namedtuple('AxisAlignedBoundingBox', 'center dimensions')
+Plane = namedtuple('Plane', 'normal distance')
+Circle = namedtuple('Circle', 'orientation center radius')
+
+def circleOfTriangle(a, b, c):
+    # https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates_from_cross-_and_dot-products
+    dirBA = a-b
+    dirCB = b-c
+    dirAC = c-a
+    normal = dirBA.cross(dirCB)
+    lengthBA = dirBA.length
+    lengthCB = dirCB.length
+    lengthAC = dirAC.length
+    lengthN = normal.length
+    if lengthN == 0:
+        return None
+    factor = -1/(2*lengthN*lengthN)
+    alpha = (dirBA@dirAC)*(lengthCB*lengthCB*factor)
+    beta = (dirBA@dirCB)*(lengthAC*lengthAC*factor)
+    gamma = (dirAC@dirCB)*(lengthBA*lengthBA*factor)
+    center = a*alpha+b*beta+c*gamma
+    radius = (lengthBA*lengthCB*lengthAC)/(2*lengthN)
+    tangent = (a-center).normalized()
+    orientation = Matrix.Identity(3)
+    orientation.col[2] = normal/lengthN
+    orientation.col[1] = (a-center).normalized()
+    orientation.col[0] = orientation.col[1].xyz.cross(orientation.col[2].xyz)
+    return Circle(orientation=orientation, center=center, radius=radius)
+
+def circleOfBezier(points, tollerance=0.000001, samples=16):
+    circle = circleOfTriangle(points[0], bezierPointAt(points, 0.5), points[3])
+    if circle == None:
+        return None
+    variance = 0
+    for t in range(0, samples):
+        variance += ((circle.center-bezierPointAt(points, (t+1)/(samples-1))).length/circle.radius-1) ** 2
+    variance /= samples
+    return None if variance > tollerance else circle
+
+def areaOfPolygon(vertices):
+    area = 0
+    for index, current in enumerate(vertices):
+        prev = vertices[index-1]
+        area += (current[0]+prev[0])*(current[1]-prev[1])
+    return area*0.5
+
+def linePointDistance(begin, dir, point):
+    return (point-begin).cross(dir.normalized()).length
+
+def linePlaneIntersection(origin, dir, plane):
+    det = dir@plane.normal
+    return float('nan') if det == 0 else (plane.distance-origin@plane.normal)/det
+
+def nearestPointOfLines(originA, dirA, originB, dirB, tollerance=0.0):
+    # https://en.wikipedia.org/wiki/Skew_lines#Nearest_Points
+    normal = dirA.cross(dirB)
+    normalA = dirA.cross(normal)
+    normalB = dirB.cross(normal)
+    divisorA = dirA@normalB
+    divisorB = dirB@normalA
+    if abs(divisorA) <= tollerance or abs(divisorB) <= tollerance:
+        return (float('nan'), float('nan'), None, None)
+    else:
+        paramA = (originB-originA)@normalB/divisorA
+        paramB = (originA-originB)@normalA/divisorB
+        return (paramA, paramB, originA+dirA*paramA, originB+dirB*paramB)
+
+def lineSegmentLineSegmentIntersection(beginA, endA, beginB, endB, tollerance=0.001):
+    dirA = endA-beginA
+    dirB = endB-beginB
+    intersection = nearestPointOfLines(beginA, dirA, beginB, dirB)
+    if math.isnan(intersection[0]) or (intersection[2]-intersection[3]).length > tollerance or \
+       intersection[0] < 0 or intersection[0] > 1 or intersection[1] < 0 or intersection[1] > 1:
+        return None
+    return intersection
+
+def aabbOfPoints(points):
+    min = Vector(points[0])
+    max = Vector(points[0])
+    for point in points:
+        for i in range(0, 3):
+            if min[i] > point[i]:
+                min[i] = point[i]
+            if max[i] < point[i]:
+                max[i] = point[i]
+    return AABB(center=(max+min)*0.5, dimensions=(max-min)*0.5)
+
+def aabbIntersectionTest(a, b, tollerance=0.0):
+    for i in range(0, 3):
+        if abs(a.center[i]-b.center[i]) > a.dimensions[i]+b.dimensions[i]+tollerance:
+            return False
+    return True
+
+def isPointInAABB(point, aabb, tollerance=0.0, ignore_axis=None):
+    for i in range(0, 3):
+        if i != ignore_axis and (point[i] < aabb.center[i]-aabb.dimensions[i]-tollerance or point[i] > aabb.center[i]+aabb.dimensions[i]+tollerance):
+            return False
+    return True
+
+def lineAABBIntersection(lineBegin, lineEnd, aabb):
+    intersections = []
+    for i in range(0, 3):
+        normal = [0, 0, 0]
+        normal = Vector(normal[0:i] + [1] + normal[i+1:])
+        for j in range(-1, 2, 2):
+            plane = Plane(normal=normal, distance=aabb.center[i]+j*aabb.dimensions[i])
+            param = linePlaneIntersection(lineBegin, lineEnd-lineBegin, plane)
+            if param < 0 or param > 1 or math.isnan(param):
+                continue
+            point = lineBegin+param*(lineEnd-lineBegin)
+            if isPointInAABB(point, aabb, 0.0, i):
+                intersections.append((param, point))
+    return intersections
+
+def bezierPointAt(points, t):
+    s = 1-t
+    return s*s*s*points[0] + 3*s*s*t*points[1] + 3*s*t*t*points[2] + t*t*t*points[3]
+
+def bezierTangentAt(points, t):
+    s = 1-t
+    return s*s*(points[1]-points[0])+2*s*t*(points[2]-points[1])+t*t*(points[3]-points[2])
+    # return s*s*points[0] + (s*s-2*s*t)*points[1] + (2*s*t-t*t)*points[2] + t*t*points[3]
+
+def bezierLength(points, beginT=0, endT=1, samples=1024):
+    # https://en.wikipedia.org/wiki/Arc_length#Finding_arc_lengths_by_integrating
+    vec = [points[1]-points[0], points[2]-points[1], points[3]-points[2]]
+    dot = [vec[0]@vec[0], vec[0]@vec[1], vec[0]@vec[2], vec[1]@vec[1], vec[1]@vec[2], vec[2]@vec[2]]
+    factors = [
+        dot[0],
+        4*(dot[1]-dot[0]),
+        6*dot[0]+4*dot[3]+2*dot[2]-12*dot[1],
+        12*dot[1]+4*(dot[4]-dot[0]-dot[2])-8*dot[3],
+        dot[0]+dot[5]+2*dot[2]+4*(dot[3]-dot[1]-dot[4])
+    ]
+    # https://en.wikipedia.org/wiki/Trapezoidal_rule
+    length = 0
+    prev_value = math.sqrt(factors[4]+factors[3]+factors[2]+factors[1]+factors[0])
+    for index in range(0, samples+1):
+        t = beginT+(endT-beginT)*index/samples
+        # value = math.sqrt(factors[4]*(t**4)+factors[3]*(t**3)+factors[2]*(t**2)+factors[1]*t+factors[0])
+        value = math.sqrt((((factors[4]*t+factors[3])*t+factors[2])*t+factors[1])*t+factors[0])
+        length += (prev_value+value)*0.5
+        prev_value = value
+    return length*3/samples
+
+# https://en.wikipedia.org/wiki/Root_of_unity
+# cubic_roots_of_unity = [cmath.rect(1, i/3*2*math.pi) for i in range(0, 3)]
+cubic_roots_of_unity = [complex(1, 0), complex(-1, math.sqrt(3))*0.5, complex(-1, -math.sqrt(3))*0.5]
+def bezierRoots(dists, tollerance=0.0001):
+    # https://en.wikipedia.org/wiki/Cubic_function
+    # y(t) = a*t^3 +b*t^2 +c*t^1 +d*t^0
+    a = 3*(dists[1]-dists[2])+dists[3]-dists[0]
+    b = 3*(dists[0]-2*dists[1]+dists[2])
+    c = 3*(dists[1]-dists[0])
+    d = dists[0]
+    if abs(a) > tollerance: # Cubic
+        E2 = a*c
+        E3 = a*a*d
+        A = (2*b*b-9*E2)*b+27*E3
+        B = b*b-3*E2
+        C = ((A+cmath.sqrt(A*A-4*B*B*B))*0.5) ** (1/3)
+        roots = []
+        for root in cubic_roots_of_unity:
+            root *= C
+            root = -1/(3*a)*(b+root+B/root)
+            if abs(root.imag) < tollerance and root.real > -param_tollerance and root.real < 1.0+param_tollerance:
+                roots.append(max(0.0, min(root.real, 1.0)))
+        # Remove doubles
+        roots.sort()
+        for index in range(len(roots)-1, 0, -1):
+            if abs(roots[index-1]-roots[index]) < param_tollerance:
+                roots.pop(index)
+        return roots
+    elif abs(b) > tollerance: # Quadratic
+        disc = c*c-4*b*d
+        if disc < 0:
+            return []
+        disc = math.sqrt(disc)
+        return [(-c-disc)/(2*b), (-c+disc)/(2*b)]
+    elif abs(c) > tollerance: # Linear
+        root = -d/c
+        return [root] if root >= 0.0 and root <= 1.0 else []
+    else: # Constant / Parallel
+        return [] if abs(d) > tollerance else float('inf')
+
+def xRaySplineIntersectionTest(spline, origin):
+    spline_points = spline.bezier_points if spline.type == 'BEZIER' else spline.points
+    cyclic_parallel_fix_flag = False
+    intersections = []
+
+    def areIntersectionsAdjacent(index):
+        if len(intersections) < 2:
+            return
+        prev = intersections[index-1]
+        current = intersections[index]
+        if prev[1] == current[0] and \
+           prev[2] > 1.0-param_tollerance and current[2] < param_tollerance and \
+           ((prev[3] < 0 and current[3] < 0) or (prev[3] > 0 and current[3] > 0)):
+            intersections.pop(index)
+
+    def appendIntersection(index, root, tangentY, intersectionX):
+        beginPoint = spline_points[index-1]
+        endPoint = spline_points[index]
+        if root == float('inf'): # Segment is parallel to ray
+            if index == 0 and spline.use_cyclic_u:
+                cyclic_parallel_fix_flag = True
+            if len(intersections) > 0 and intersections[-1][1] == beginPoint:
+                intersections[-1][1] = endPoint # Skip in adjacency test
+        elif intersectionX >= origin[0]:
+            intersections.append([beginPoint, endPoint, root, tangentY, intersectionX])
+            areIntersectionsAdjacent(len(intersections)-1)
+
+    if spline.type == 'BEZIER':
+        for index, endPoint in enumerate(spline.bezier_points):
+            if index == 0 and not spline.use_cyclic_u:
+                continue
+            beginPoint = spline_points[index-1]
+            points = (beginPoint.co, beginPoint.handle_right, endPoint.handle_left, endPoint.co)
+            roots = bezierRoots((points[0][1]-origin[1], points[1][1]-origin[1], points[2][1]-origin[1], points[3][1]-origin[1]))
+            if roots == float('inf'): # Intersection
+                appendIntersection(index, float('inf'), None, None)
+            else:
+                for root in roots:
+                    appendIntersection(index, root, bezierTangentAt(points, root)[1], bezierPointAt(points, root)[0])
+    elif spline.type == 'POLY':
+        for index, endPoint in enumerate(spline.points):
+            if index == 0 and not spline.use_cyclic_u:
+                continue
+            beginPoint = spline_points[index-1]
+            points = (beginPoint.co, endPoint.co)
+            if (points[0][0] < origin[0] and points[1][0] < origin[0]) or \
+               (points[0][1] < origin[1] and points[1][1] < origin[1]) or \
+               (points[0][1] > origin[1] and points[1][1] > origin[1]):
+                continue
+            diff = points[1]-points[0]
+            height = origin[1]-points[0][1]
+            if diff[1] == 0: # Parallel
+                if height == 0: # Intersection
+                    appendIntersection(index, float('inf'), None, None)
+            else: # Not parallel
+                root = height/diff[1]
+                appendIntersection(index, root, diff[1], points[0][0]+diff[0]*root)
+
+    if cyclic_parallel_fix_flag:
+        appendIntersection(0, float('inf'), None, None)
+    areIntersectionsAdjacent(0)
+    return intersections
+
+def isPointInSpline(point, spline):
+    return spline.use_cyclic_u and len(xRaySplineIntersectionTest(spline, point))%2 == 1
+
+def isSegmentLinear(points, tollerance=0.0001):
+    return 1.0-(points[1]-points[0]).normalized()@(points[3]-points[2]).normalized() < tollerance
+
+def bezierSegmentPoints(begin, end):
+    return [begin.co, begin.handle_right, end.handle_left, end.co]
+
+def deleteFromArray(item, array):
+    for index, current in enumerate(array):
+        if current is item:
+            array.pop(index)
+            break
+
+def copyAttributes(dst, src):
+    for attribute in dir(src):
+        try:
+            setattr(dst, attribute, getattr(src, attribute))
+        except:
+            pass
+
+def bezierSliceFromTo(points, minParam, maxParam):
+    fromP = bezierPointAt(points, minParam)
+    fromT = bezierTangentAt(points, minParam)
+    toP = bezierPointAt(points, maxParam)
+    toT = bezierTangentAt(points, maxParam)
+    paramDiff = maxParam-minParam
+    return [fromP, fromP+fromT*paramDiff, toP-toT*paramDiff, toP]
+
+def bezierIntersectionBroadPhase(solutions, pointsA, pointsB, aMin=0.0, aMax=1.0, bMin=0.0, bMax=1.0, depth=8, tollerance=0.001):
+    if aabbIntersectionTest(aabbOfPoints(bezierSliceFromTo(pointsA, aMin, aMax)), aabbOfPoints(bezierSliceFromTo(pointsB, bMin, bMax)), tollerance) == False:
+        return
+    if depth == 0:
+        solutions.append([aMin, aMax, bMin, bMax])
+        return
+    depth -= 1
+    aMid = (aMin+aMax)*0.5
+    bMid = (bMin+bMax)*0.5
+    bezierIntersectionBroadPhase(solutions, pointsA, pointsB, aMin, aMid, bMin, bMid, depth, tollerance)
+    bezierIntersectionBroadPhase(solutions, pointsA, pointsB, aMin, aMid, bMid, bMax, depth, tollerance)
+    bezierIntersectionBroadPhase(solutions, pointsA, pointsB, aMid, aMax, bMin, bMid, depth, tollerance)
+    bezierIntersectionBroadPhase(solutions, pointsA, pointsB, aMid, aMax, bMid, bMax, depth, tollerance)
+
+def bezierIntersectionNarrowPhase(broadPhase, pointsA, pointsB, tollerance=0.000001):
+    aMin = broadPhase[0]
+    aMax = broadPhase[1]
+    bMin = broadPhase[2]
+    bMax = broadPhase[3]
+    while (aMax-aMin > tollerance) or (bMax-bMin > tollerance):
+        aMid = (aMin+aMax)*0.5
+        bMid = (bMin+bMax)*0.5
+        a1 = bezierPointAt(pointsA, (aMin+aMid)*0.5)
+        a2 = bezierPointAt(pointsA, (aMid+aMax)*0.5)
+        b1 = bezierPointAt(pointsB, (bMin+bMid)*0.5)
+        b2 = bezierPointAt(pointsB, (bMid+bMax)*0.5)
+        a1b1Dist = (a1-b1).length
+        a2b1Dist = (a2-b1).length
+        a1b2Dist = (a1-b2).length
+        a2b2Dist = (a2-b2).length
+        minDist = min(a1b1Dist, a2b1Dist, a1b2Dist, a2b2Dist)
+        if a1b1Dist == minDist:
+            aMax = aMid
+            bMax = bMid
+        elif a2b1Dist == minDist:
+            aMin = aMid
+            bMax = bMid
+        elif a1b2Dist == minDist:
+            aMax = aMid
+            bMin = bMid
+        else:
+            aMin = aMid
+            bMin = bMid
+    return [aMin, bMin, minDist]
+
+def segmentIntersection(segmentA, segmentB, tollerance=0.001):
+    pointsA = bezierSegmentPoints(segmentA['beginPoint'], segmentA['endPoint'])
+    pointsB = bezierSegmentPoints(segmentB['beginPoint'], segmentB['endPoint'])
+    result = []
+    def addCut(paramA, paramB):
+        cutA = {'param': paramA, 'segment': segmentA}
+        cutB = {'param': paramB, 'segment': segmentB}
+        cutA['otherCut'] = cutB
+        cutB['otherCut'] = cutA
+        segmentA['cuts'].append(cutA)
+        segmentB['cuts'].append(cutB)
+        result.append([cutA, cutB])
+    if isSegmentLinear(pointsA) and isSegmentLinear(pointsB):
+        intersection = lineSegmentLineSegmentIntersection(pointsA[0], pointsA[3], pointsB[0], pointsB[3])
+        if intersection != None:
+            addCut(intersection[0], intersection[1])
+        return result
+    solutions = []
+    bezierIntersectionBroadPhase(solutions, pointsA, pointsB)
+    for index in range(0, len(solutions)):
+        solutions[index] = bezierIntersectionNarrowPhase(solutions[index], pointsA, pointsB)
+    for index in range(0, len(solutions)):
+        for otherIndex in range(0, len(solutions)):
+            if solutions[index][2] == float('inf'):
+                break
+            if index == otherIndex or solutions[otherIndex][2] == float('inf'):
+                continue
+            diffA = solutions[index][0]-solutions[otherIndex][0]
+            diffB = solutions[index][1]-solutions[otherIndex][1]
+            if diffA*diffA+diffB*diffB < 0.01:
+                if solutions[index][2] < solutions[otherIndex][2]:
+                    solutions[otherIndex][2] = float('inf')
+                else:
+                    solutions[index][2] = float('inf')
+    def areIntersectionsAdjacent(segmentA, segmentB, paramA, paramB):
+        return segmentA['endIndex'] == segmentB['beginIndex'] and paramA > 1-param_tollerance and paramB < param_tollerance
+    for solution in solutions:
+        if (solution[2] > tollerance) or \
+          (segmentA['spline'] == segmentB['spline'] and \
+          (areIntersectionsAdjacent(segmentA, segmentB, solution[0], solution[1]) or \
+           areIntersectionsAdjacent(segmentB, segmentA, solution[1], solution[0]))):
+            continue
+        addCut(solution[0], solution[1])
+    return result
+
+def bezierMultiIntersection(segments):
+    for index in range(0, len(segments)):
+        for otherIndex in range(index+1, len(segments)):
+            segmentIntersection(segments[index], segments[otherIndex])
+    prepareSegmentIntersections(segments)
+    subdivideBezierSegments(segments)
+
+def bezierSubivideAt(points, params):
+    if len(params) == 0:
+        return []
+    newPoints = []
+    newPoints.append(points[0]+(points[1]-points[0])*params[0])
+    for index, param in enumerate(params):
+        paramLeft = param
+        if index > 0:
+            paramLeft -= params[index-1]
+        paramRight = -param
+        if index == len(params)-1:
+            paramRight += 1.0
+        else:
+            paramRight += params[index+1]
+        point = bezierPointAt(points, param)
+        tangent = bezierTangentAt(points, param)
+        newPoints.append(point-tangent*paramLeft)
+        newPoints.append(point)
+        newPoints.append(point+tangent*paramRight)
+    newPoints.append(points[3]-(points[3]-points[2])*(1.0-params[-1]))
+    return newPoints
+
+def subdivideBezierSegment(segment):
+    # Blender only allows uniform subdivision. Use this method to subdivide at arbitrary params.
+    # NOTE: segment['cuts'] must be sorted by param
+    if len(segment['cuts']) == 0:
+        return
+
+    segment['beginPoint'] = segment['spline'].bezier_points[segment['beginIndex']]
+    segment['endPoint'] = segment['spline'].bezier_points[segment['endIndex']]
+    params = [cut['param'] for cut in segment['cuts']]
+    newPoints = bezierSubivideAt(bezierSegmentPoints(segment['beginPoint'], segment['endPoint']), params)
+    bpy.ops.curve.select_all(action='DESELECT')
+    segment['beginPoint'] = segment['spline'].bezier_points[segment['beginIndex']]
+    segment['beginPoint'].select_right_handle = True
+    segment['beginPoint'].handle_left_type = 'FREE'
+    segment['beginPoint'].handle_right_type = 'FREE'
+    segment['endPoint'] = segment['spline'].bezier_points[segment['endIndex']]
+    segment['endPoint'].select_left_handle = True
+    segment['endPoint'].handle_left_type = 'FREE'
+    segment['endPoint'].handle_right_type = 'FREE'
+
+    bpy.ops.curve.subdivide(number_cuts=len(params))
+    if segment['endIndex'] > 0:
+        segment['endIndex'] += len(params)
+    segment['beginPoint'] = segment['spline'].bezier_points[segment['beginIndex']]
+    segment['endPoint'] = segment['spline'].bezier_points[segment['endIndex']]
+    segment['beginPoint'].select_right_handle = False
+    segment['beginPoint'].handle_right = newPoints[0]
+    segment['endPoint'].select_left_handle = False
+    segment['endPoint'].handle_left = newPoints[-1]
+
+    for index, cut in enumerate(segment['cuts']):
+        cut['index'] = segment['beginIndex']+1+index
+        newPoint = segment['spline'].bezier_points[cut['index']]
+        newPoint.handle_left_type = 'FREE'
+        newPoint.handle_right_type = 'FREE'
+        newPoint.select_left_handle = False
+        newPoint.select_control_point = False
+        newPoint.select_right_handle = False
+        newPoint.handle_left = newPoints[index*3+1]
+        newPoint.co = newPoints[index*3+2]
+        newPoint.handle_right = newPoints[index*3+3]
+
+def prepareSegmentIntersections(segments):
+    def areCutsAdjacent(cutA, cutB):
+        return cutA['segment']['beginIndex'] == cutB['segment']['endIndex'] and \
+               cutA['param'] < param_tollerance and cutB['param'] > 1.0-param_tollerance
+    for segment in segments:
+        segment['cuts'].sort(key=(lambda cut: cut['param']))
+        for index in range(len(segment['cuts'])-1, 0, -1):
+            prev = segment['cuts'][index-1]
+            current = segment['cuts'][index]
+            if abs(prev['param']-current['param']) < param_tollerance and \
+               prev['otherCut']['segment']['spline'] == current['otherCut']['segment']['spline'] and \
+               (areCutsAdjacent(prev['otherCut'], current['otherCut']) or \
+                areCutsAdjacent(current['otherCut'], prev['otherCut'])):
+                deleteFromArray(prev['otherCut'], prev['otherCut']['segment']['cuts'])
+                deleteFromArray(current['otherCut'], current['otherCut']['segment']['cuts'])
+                segment['cuts'].pop(index-1 if current['otherCut']['param'] < param_tollerance else index)
+                current = segment['cuts'][index-1]['otherCut']
+                current['segment']['extraCut'] = current
+
+def subdivideBezierSegmentsOfSameSpline(segments):
+    # NOTE: segment['cuts'] must be sorted by param
+    indexOffset = 0
+    for segment in segments:
+        segment['beginIndex'] += indexOffset
+        if segment['endIndex'] > 0:
+            segment['endIndex'] += indexOffset
+        subdivideBezierSegment(segment)
+        indexOffset += len(segment['cuts'])
+    for segment in segments:
+        segment['beginPoint'] = segment['spline'].bezier_points[segment['beginIndex']]
+        segment['endPoint'] = segment['spline'].bezier_points[segment['endIndex']]
+
+def subdivideBezierSegments(segments):
+    # NOTE: segment['cuts'] must be sorted by param
+    groups = {}
+    for segment in segments:
+        spline = segment['spline']
+        if (spline in groups) == False:
+            groups[spline] = []
+        group = groups[spline]
+        group.append(segment)
+    for spline in groups:
+        subdivideBezierSegmentsOfSameSpline(groups[spline])
+
+def curveObject():
+    obj = bpy.context.object
+    return obj if obj != None and obj.type == 'CURVE' and obj.mode == 'EDIT' else None
+
+def bezierSegments(splines, selection_only):
+    segments = []
+    for spline in splines:
+        if spline.type != 'BEZIER':
+            continue
+        for index, current in enumerate(spline.bezier_points):
+            next = spline.bezier_points[(index+1) % len(spline.bezier_points)]
+            if next == spline.bezier_points[0] and not spline.use_cyclic_u:
+                continue
+            if not selection_only or (current.select_right_handle and next.select_left_handle):
+                segments.append({
+                    'spline': spline,
+                    'beginIndex': index,
+                    'endIndex': index+1 if index < len(spline.bezier_points)-1 else 0,
+                    'beginPoint': current,
+                    'endPoint': next,
+                    'cuts': []
+                })
+    return segments
+
+def getSelectedSplines(include_bezier, include_polygon, allow_partial_selection=False):
+    result = []
+    for spline in bpy.context.object.data.splines:
+        selected = not allow_partial_selection
+        if spline.type == 'BEZIER':
+            if not include_bezier:
+                continue
+            for index, point in enumerate(spline.bezier_points):
+                if point.select_left_handle == allow_partial_selection or \
+                   point.select_control_point == allow_partial_selection or \
+                   point.select_right_handle == allow_partial_selection:
+                    selected = allow_partial_selection
+                    break
+        elif spline.type == 'POLY':
+            if not include_polygon:
+                continue
+            for index, point in enumerate(spline.points):
+                if point.select == allow_partial_selection:
+                    selected = allow_partial_selection
+                    break
+        else:
+            continue
+        if selected:
+            result.append(spline)
+    return result
+
+def addObject(type, name):
+    bpy.ops.object.select_all(action='DESELECT')
+    if type == 'CURVE':
+        data = bpy.data.curves.new(name=name, type='CURVE')
+        data.dimensions = '3D'
+    elif type == 'MESH':
+        data = bpy.data.meshes.new(name=name, type='MESH')
+    obj = bpy.data.objects.new(name, data)
+    obj.location = bpy.context.scene.cursor.location
+    bpy.context.scene.collection.objects.link(obj)
+    obj.select_set(True)
+    bpy.context.view_layer.objects.active = obj
+    return obj
+
+def addPolygonSpline(obj, cyclic, vertices, weights=None, select=False):
+    spline = obj.data.splines.new(type='POLY')
+    spline.use_cyclic_u = cyclic
+    spline.points.add(len(vertices)-1)
+    for index, point in enumerate(spline.points):
+        point.co.xyz = vertices[index]
+        point.select = select
+        if weights:
+            point.weight_softbody = weights[index]
+    return spline
+
+def addBezierSpline(obj, cyclic, vertices, weights=None, select=False):
+    spline = obj.data.splines.new(type='BEZIER')
+    spline.use_cyclic_u = cyclic
+    spline.bezier_points.add(len(vertices)-1)
+    for index, point in enumerate(spline.bezier_points):
+        point.handle_left = vertices[index][0]
+        point.co = vertices[index][1]
+        point.handle_right = vertices[index][2]
+        if weights:
+            point.weight_softbody = weights[index]
+        point.select_left_handle = select
+        point.select_control_point = select
+        point.select_right_handle = select
+        if isSegmentLinear([vertices[index-1][1], vertices[index-1][2], vertices[index][0], vertices[index][1]]):
+            spline.bezier_points[index-1].handle_right_type = 'VECTOR'
+            point.handle_left_type = 'VECTOR'
+    return spline
+
+def polygonArcAt(center, radius, begin_angle, angle, step_angle, include_ends):
+    vertices = []
+    circle_samples = math.ceil(abs(angle)/step_angle)
+    for t in (range(0, circle_samples+1) if include_ends else range(1, circle_samples)):
+        t = begin_angle+angle*t/circle_samples
+        normal = Vector((math.cos(t), math.sin(t), 0))
+        vertices.append(center+normal*radius)
+    return vertices
+
+def bezierArcAt(tangent, normal, center, radius, angle, tollerance=0.99999):
+    transform = Matrix.Identity(4)
+    transform.col[0].xyz = tangent.cross(normal)*radius
+    transform.col[1].xyz = tangent*radius
+    transform.col[2].xyz = normal*radius
+    transform.col[3].xyz = center
+    segments = []
+    segment_count = math.ceil(abs(angle)/(math.pi*0.5)*tollerance)
+    angle /= segment_count
+    x0 = math.cos(angle*0.5)
+    y0 = math.sin(angle*0.5)
+    x1 = (4.0-x0)/3.0
+    y1 = (1.0-x0)*(3.0-x0)/(3.0*y0)
+    points = [
+        Vector((x0, -y0, 0)),
+        Vector((x1, -y1, 0)),
+        Vector((x1, y1, 0)),
+        Vector((x0, y0, 0))
+    ]
+    for i in range(0, segment_count):
+        rotation = Matrix.Rotation((i+0.5)*angle, 4, 'Z')
+        segments.append(list(map(lambda v: transform@(rotation@v), points)))
+    return segments
+
+def iterateSpline(spline, callback):
+    spline_points = spline.bezier_points if spline.type == 'BEZIER' else spline.points
+    for index, spline_point in enumerate(spline_points):
+        prev = spline_points[index-1]
+        current = spline_points[index]
+        next = spline_points[(index+1)%len(spline_points)]
+        if spline.type == 'BEZIER':
+            selected = current.select_control_point
+            prev_segment_points = bezierSegmentPoints(prev, current)
+            next_segment_points = bezierSegmentPoints(current, next)
+            prev_tangent = (prev_segment_points[3]-prev_segment_points[2]).normalized()
+            current_tangent = (next_segment_points[1]-next_segment_points[0]).normalized()
+            next_tangent = (next_segment_points[3]-next_segment_points[2]).normalized()
+        else:
+            selected = current.select
+            prev_segment_points = [prev.co.xyz, None, None, current.co.xyz]
+            next_segment_points = [current.co.xyz, None, None, next.co.xyz]
+            prev_tangent = (prev_segment_points[3]-prev_segment_points[0]).normalized()
+            current_tangent = next_tangent = (next_segment_points[3]-next_segment_points[0]).normalized()
+        normal = prev_tangent.cross(current_tangent).normalized()
+        angle = prev_tangent@current_tangent
+        angle = 0 if abs(angle-1.0) < 0.0001 else math.acos(angle)
+        is_first = (index == 0) and not spline.use_cyclic_u
+        is_last = (index == len(spline_points)-1) and not spline.use_cyclic_u
+        callback(prev_segment_points, next_segment_points, selected, prev_tangent, current_tangent, next_tangent, normal, angle, is_first, is_last)
+    return spline_points
+
+def offsetPolygonOfSpline(spline, offset, step_angle, round_line_join, bezier_samples=128, tollerance=0.000001):
+    def offsetVertex(position, tangent):
+        normal = Vector((-tangent[1], tangent[0], 0))
+        return position+normal*offset
+    vertices = []
+    def handlePoint(prev_segment_points, next_segment_points, selected, prev_tangent, current_tangent, next_tangent, normal, angle, is_first, is_last):
+        sign = math.copysign(1, normal[2])
+        angle *= sign
+        if is_last:
+            return
+        is_protruding = (abs(angle) > tollerance and abs(offset) > tollerance)
+        if is_protruding and not is_first and sign != math.copysign(1, offset): # Convex Corner
+            if round_line_join:
+                begin_angle = math.atan2(prev_tangent[1], prev_tangent[0])+math.pi*0.5
+                vertices.extend(polygonArcAt(next_segment_points[0], offset, begin_angle, angle, step_angle, False))
+            else:
+                distance = offset*math.tan(angle*0.5)
+                vertices.append(offsetVertex(next_segment_points[0], current_tangent)+current_tangent*distance)
+        if is_protruding or is_first:
+            vertices.append(offsetVertex(next_segment_points[0], current_tangent))
+        if spline.type == 'POLY' or isSegmentLinear(next_segment_points):
+            vertices.append(offsetVertex(next_segment_points[3], next_tangent))
+        else: # Trace Bezier Segment
+            prev_tangent = bezierTangentAt(next_segment_points, 0).normalized()
+            for t in range(1, bezier_samples+1):
+                t /= bezier_samples
+                tangent = bezierTangentAt(next_segment_points, t).normalized()
+                if t == 1 or math.acos(min(max(-1, prev_tangent@tangent), 1)) >= step_angle:
+                    vertices.append(offsetVertex(bezierPointAt(next_segment_points, t), tangent))
+                    prev_tangent = tangent
+    spline_points = iterateSpline(spline, handlePoint)
+
+    # Solve Self Intersections
+    original_area = areaOfPolygon([point.co for point in spline_points])
+    sign = -1 if offset < 0 else 1
+    i = (0 if spline.use_cyclic_u else 1)
+    while i < len(vertices):
+        j = i+2
+        while j < len(vertices) - (0 if i > 0 else 1):
+            intersection = lineSegmentLineSegmentIntersection(vertices[i-1], vertices[i], vertices[j-1], vertices[j])
+            if intersection == None:
+                j += 1
+                continue
+            intersection = (intersection[2]+intersection[3])*0.5
+            areaInner = sign*areaOfPolygon([intersection, vertices[i], vertices[j-1]])
+            areaOuter = sign*areaOfPolygon([intersection, vertices[j], vertices[i-1]])
+            if areaInner > areaOuter:
+                vertices = vertices[i:j]+[intersection]
+                i = (0 if spline.use_cyclic_u else 1)
+            else:
+                vertices = vertices[:i]+[intersection]+vertices[j:]
+            j = i+2
+        i += 1
+    new_area = areaOfPolygon(vertices)
+    return [vertices] if original_area*new_area >= 0 else []
+
+def filletSpline(spline, radius, chamfer_mode, limit_half_way, tollerance=0.0001):
+    vertices = []
+    distance_limit_factor = 0.5 if limit_half_way else 1.0
+    def handlePoint(prev_segment_points, next_segment_points, selected, prev_tangent, current_tangent, next_tangent, normal, angle, is_first, is_last):
+        distance = min((prev_segment_points[0]-prev_segment_points[3]).length*distance_limit_factor, (next_segment_points[0]-next_segment_points[3]).length*distance_limit_factor)
+        if not selected or is_first or is_last or angle == 0 or distance == 0 or \
+           (spline.type == 'BEZIER' and not (isSegmentLinear(prev_segment_points) and isSegmentLinear(next_segment_points))):
+            prev_handle = next_segment_points[0] if is_first else prev_segment_points[2] if spline.type == 'BEZIER' else prev_segment_points[0]
+            next_handle = next_segment_points[0] if is_last else next_segment_points[1] if spline.type == 'BEZIER' else next_segment_points[3]
+            vertices.append([prev_handle, next_segment_points[0], next_handle])
+            return
+        tan_factor = math.tan(angle*0.5)
+        offset = min(radius, distance/tan_factor)
+        distance = offset*tan_factor
+        circle_center = next_segment_points[0]+normal.cross(prev_tangent)*offset-prev_tangent*distance
+        segments = bezierArcAt(prev_tangent, normal, circle_center, offset, angle)
+        if chamfer_mode:
+            vertices.append([prev_segment_points[0], segments[0][0], segments[-1][3]])
+            vertices.append([segments[0][0], segments[-1][3], next_segment_points[3]])
+        else:
+            for i in range(0, len(segments)+1):
+                vertices.append([
+                    segments[i-1][2] if i > 0 else prev_segment_points[0],
+                    segments[i][0] if i < len(segments) else segments[i-1][3],
+                    segments[i][1] if i < len(segments) else next_segment_points[3]
+                ])
+    iterateSpline(spline, handlePoint)
+    i = 0 if spline.use_cyclic_u else 1
+    while(i < len(vertices)):
+        if (vertices[i-1][1]-vertices[i][1]).length < tollerance:
+            vertices[i-1][2] = vertices[i][2]
+            del vertices[i]
+        else:
+            i = i+1
+    return addBezierSpline(bpy.context.object, spline.use_cyclic_u, vertices)
+
+def discretizeCurve(spline, step_angle, samples):
+    vertices = []
+    def handlePoint(prev_segment_points, next_segment_points, selected, prev_tangent, current_tangent, next_tangent, normal, angle, is_first, is_last):
+        if is_last:
+            return
+        if isSegmentLinear(next_segment_points):
+            vertices.append(next_segment_points[3])
+        else:
+            prev_tangent = bezierTangentAt(next_segment_points, 0).normalized()
+            for t in range(1, samples+1):
+                t /= samples
+                tangent = bezierTangentAt(next_segment_points, t).normalized()
+                if t == 1 or math.acos(min(max(-1, prev_tangent@tangent), 1)) >= step_angle:
+                    vertices.append(bezierPointAt(next_segment_points, t))
+                    prev_tangent = tangent
+    iterateSpline(spline, handlePoint)
+    return vertices
+
+def bezierBooleanGeometry(splineA, splineB, operation):
+    if not splineA.use_cyclic_u or not splineB.use_cyclic_u:
+        return False
+    segmentsA = bezierSegments([splineA], False)
+    segmentsB = bezierSegments([splineB], False)
+
+    deletionFlagA = isPointInSpline(splineA.bezier_points[0].co, splineB)
+    deletionFlagB = isPointInSpline(splineB.bezier_points[0].co, splineA)
+    if operation == 'DIFFERENCE':
+        deletionFlagB = not deletionFlagB
+    elif operation == 'INTERSECTION':
+        deletionFlagA = not deletionFlagA
+        deletionFlagB = not deletionFlagB
+    elif operation != 'UNION':
+        return False
+
+    intersections = []
+    for segmentA in segmentsA:
+        for segmentB in segmentsB:
+            intersections.extend(segmentIntersection(segmentA, segmentB))
+    if len(intersections) == 0:
+        if deletionFlagA:
+            bpy.context.object.data.splines.remove(splineA)
+        if deletionFlagB:
+            bpy.context.object.data.splines.remove(splineB)
+        return True
+
+    prepareSegmentIntersections(segmentsA)
+    prepareSegmentIntersections(segmentsB)
+    subdivideBezierSegmentsOfSameSpline(segmentsA)
+    subdivideBezierSegmentsOfSameSpline(segmentsB)
+
+    def collectCuts(cuts, segments, deletionFlag):
+        for segmentIndex, segment in enumerate(segments):
+            if 'extraCut' in segment:
+                deletionFlag = not deletionFlag
+                segment['extraCut']['index'] = segment['beginIndex']
+                segment['extraCut']['deletionFlag'] = deletionFlag
+                cuts.append(segment['extraCut'])
+            else:
+                cuts.append(None)
+            cuts.extend(segments[segmentIndex]['cuts'])
+            segment['deletionFlag'] = deletionFlag
+            for cutIndex, cut in enumerate(segment['cuts']):
+                deletionFlag = not deletionFlag
+                cut['deletionFlag'] = deletionFlag
+    cutsA = []
+    cutsB = []
+    collectCuts(cutsA, segmentsA, deletionFlagA)
+    collectCuts(cutsB, segmentsB, deletionFlagB)
+
+    beginIndex = 0
+    for segment in segmentsA:
+        if segment['deletionFlag'] == False:
+            beginIndex = segment['beginIndex']
+            break
+        for cut in segment['cuts']:
+            if cut['deletionFlag'] == False:
+                beginIndex = cut['index']
+                break
+
+    cuts = cutsA
+    spline = splineA
+    index = beginIndex
+    backward = False
+    vertices = []
+    while True:
+        current = spline.bezier_points[index]
+        vertices.append([current.handle_left, current.co, current.handle_right])
+        if backward:
+            current.handle_left, current.handle_right = current.handle_right.copy(), current.handle_left.copy()
+        index += len(spline.bezier_points)-1 if backward else 1
+        index %= len(spline.bezier_points)
+        if spline == splineA and index == beginIndex:
+            break
+
+        cut = cuts[index]
+        if cut != None:
+            current = spline.bezier_points[index]
+            current_handle = current.handle_right if backward else current.handle_left
+            spline = splineA if spline == splineB else splineB
+            cuts = cutsA if spline == splineA else cutsB
+            index = cut['otherCut']['index']
+            backward = cut['otherCut']['deletionFlag']
+            next = spline.bezier_points[index]
+            if backward:
+                next.handle_right = current_handle
+            else:
+                next.handle_left = current_handle
+            if spline == splineA and index == beginIndex:
+                break
+
+    spline = addBezierSpline(bpy.context.object, True, vertices)
+    bpy.context.object.data.splines.remove(splineA)
+    bpy.context.object.data.splines.remove(splineB)
+    bpy.context.object.data.splines.active = spline
+    return True
+
+def truncateToFitBox(transform, spline, aabb):
+    spline_points = spline.points
+    aux = {
+        'traces': [],
+        'vertices': [],
+        'weights': []
+    }
+    def terminateTrace(aux):
+        if len(aux['vertices']) > 0:
+            aux['traces'].append((aux['vertices'], aux['weights']))
+        aux['vertices'] = []
+        aux['weights'] = []
+    for index, point in enumerate(spline_points):
+        begin = transform@point.co.xyz
+        end = spline_points[(index+1)%len(spline_points)]
+        inside = isPointInAABB(begin, aabb)
+        if inside:
+            aux['vertices'].append(begin)
+            aux['weights'].append(point.weight_softbody)
+        if index == len(spline_points)-1 and not spline.use_cyclic_u:
+            break
+        intersections = lineAABBIntersection(begin, transform@end.co.xyz, aabb)
+        if len(intersections) == 2:
+            terminateTrace(aux)
+            aux['traces'].append((
+                [intersections[0][1], intersections[1][1]],
+                [end.weight_softbody, end.weight_softbody]
+            ))
+        elif len(intersections) == 1:
+            aux['vertices'].append(intersections[0][1])
+            aux['weights'].append(end.weight_softbody)
+            if inside:
+                terminateTrace(aux)
+        elif inside and index == len(spline_points)-1 and spline.use_cyclic_u:
+            terminateTrace(aux)
+            aux['traces'][0] = (aux['traces'][-1][0]+aux['traces'][0][0], aux['traces'][-1][1]+aux['traces'][0][1])
+            aux['traces'].pop()
+    terminateTrace(aux)
+    return aux['traces']
+
+def arrayModifier(splines, offset, count, connect, serpentine):
+    if connect:
+        for spline in splines:
+            if spline.use_cyclic_u:
+                spline.use_cyclic_u = False
+                points = spline.points if spline.type == 'POLY' else spline.bezier_points
+                points.add(1)
+                copyAttributes(points[-1], points[0])
+    bpy.ops.curve.select_all(action='DESELECT')
+    for spline in splines:
+        if spline.type == 'BEZIER':
+            for point in spline.bezier_points:
+                point.select_left_handle = point.select_control_point = point.select_right_handle = True
+        elif spline.type == 'POLY':
+            for point in spline.points:
+                point.select = True
+    splines_at_layer = [splines]
+    for i in range(1, count):
+        bpy.ops.curve.duplicate()
+        bpy.ops.transform.translate(value=offset)
+        splines_at_layer.append(getSelectedSplines(True, True))
+        if serpentine:
+            bpy.ops.curve.switch_direction()
+    if connect:
+        for i in range(1, count):
+            prev_layer = splines_at_layer[i-1]
+            next_layer = splines_at_layer[i]
+            for j in range(0, len(next_layer)):
+                bpy.ops.curve.select_all(action='DESELECT')
+                if prev_layer[j].type == 'POLY':
+                    prev_layer[j].points[-1].select = True
+                else:
+                    prev_layer[j].bezier_points[-1].select_control_point = True
+                if next_layer[j].type == 'POLY':
+                    next_layer[j].points[0].select = True
+                else:
+                    next_layer[j].bezier_points[0].select_control_point = True
+                bpy.ops.curve.make_segment()
+    bpy.ops.curve.select_all(action='DESELECT')
diff --git a/curve_tools/toolpath.py b/curve_tools/toolpath.py
new file mode 100644
index 00000000..8353f7ac
--- /dev/null
+++ b/curve_tools/toolpath.py
@@ -0,0 +1,287 @@
+#  ***** GPL LICENSE BLOCK *****
+#
+#  This program is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#  All rights reserved.
+#
+#  ***** GPL LICENSE BLOCK *****
+
+import bpy, math, bmesh
+from mathutils import Vector, Matrix
+from . import internal
+
+class OffsetCurve(bpy.types.Operator):
+    bl_idname = 'curve.add_toolpath_offset_curve'
+    bl_description = bl_label = 'Offset Curve'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    offset: bpy.props.FloatProperty(name='Offset', description='Distace between the original and the first trace', unit='LENGTH', default=0.1)
+    pitch: bpy.props.FloatProperty(name='Pitch', description='Distace between two parallel traces', unit='LENGTH', default=0.1)
+    step_angle: bpy.props.FloatProperty(name='Resolution', description='Smaller values make curves smoother by adding more vertices', unit='ROTATION', min=math.pi/128, default=math.pi/16)
+    count: bpy.props.IntProperty(name='Count', description='Number of parallel traces', min=1, default=1)
+    round_line_join: bpy.props.BoolProperty(name='Round Line Join', description='Insert circle arcs at convex corners', default=True)
+
+    @classmethod
+    def poll(cls, context):
+        return bpy.context.object != None and bpy.context.object.type == 'CURVE'
+
+    def execute(self, context):
+        if bpy.context.object.mode == 'EDIT':
+            splines = internal.getSelectedSplines(True, True)
+        else:
+            splines = bpy.context.object.data.splines
+
+        if len(splines) == 0:
+            self.report({'WARNING'}, 'Nothing selected')
+            return {'CANCELLED'}
+
+        if bpy.context.object.mode != 'EDIT':
+            internal.addObject('CURVE', 'Offset Toolpath')
+            origin = bpy.context.scene.cursor.location
+        else:
+            origin = Vector((0.0, 0.0, 0.0))
+
+        for spline in splines:
+            spline_points = spline.bezier_points if spline.type == 'BEZIER' else spline.points
+            for spline_point in spline_points:
+                if spline_point.co.z != spline_points[0].co.z:
+                    self.report({'WARNING'}, 'Curves must be planar and in XY plane')
+                    return {'CANCELLED'}
+            for index in range(0, self.count):
+                traces = internal.offsetPolygonOfSpline(spline, self.offset+self.pitch*index, self.step_angle, self.round_line_join)
+                for trace in traces:
+                    internal.addPolygonSpline(bpy.context.object, spline.use_cyclic_u, [vertex-origin for vertex in trace])
+        return {'FINISHED'}
+
+class SliceMesh(bpy.types.Operator):
+    bl_idname = 'curve.add_toolpath_slice_mesh'
+    bl_description = bl_label = 'Slice Mesh'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    pitch_axis: bpy.props.FloatVectorProperty(name='Pitch & Axis', unit='LENGTH', description='Vector between to slices', subtype='DIRECTION', default=(0.0, 0.0, 0.1), size=3)
+    offset: bpy.props.FloatProperty(name='Offset', unit='LENGTH', description='Position of first slice along axis', default=-0.4)
+    slice_count: bpy.props.IntProperty(name='Count', description='Number of slices', min=1, default=9)
+
+    @classmethod
+    def poll(cls, context):
+        return bpy.context.object != None and bpy.context.object.mode == 'OBJECT'
+
+    def execute(self, context):
+        if bpy.context.object.type != 'MESH':
+            self.report({'WARNING'}, 'Active object must be a mesh')
+            return {'CANCELLED'}
+        depsgraph = context.evaluated_depsgraph_get()
+        mesh = bmesh.new()
+        mesh.from_object(bpy.context.object, depsgraph, deform=True, cage=False, face_normals=True)
+        mesh.transform(bpy.context.object.matrix_world)
+        toolpath = internal.addObject('CURVE', 'Slices Toolpath')
+        pitch_axis = Vector(self.pitch_axis)
+        axis = pitch_axis.normalized()
+        for i in range(0, self.slice_count):
+            aux_mesh = mesh.copy()
+            cut_geometry = bmesh.ops.bisect_plane(aux_mesh, geom=aux_mesh.edges[:]+aux_mesh.faces[:], dist=0, plane_co=pitch_axis*i+axis*self.offset, plane_no=axis, clear_outer=False, clear_inner=False)['geom_cut']
+            edge_pool = set([e for e in cut_geometry if isinstance(e, bmesh.types.BMEdge)])
+            while len(edge_pool) > 0:
+                current_edge = edge_pool.pop()
+                first_vertex = current_vertex = current_edge.verts[0]
+                vertices = [current_vertex.co]
+                follow_edge_loop = len(edge_pool) > 0
+                while follow_edge_loop:
+                    current_vertex = current_edge.other_vert(current_vertex)
+                    vertices.append(current_vertex.co)
+                    if current_vertex == first_vertex:
+                        break
+                    follow_edge_loop = False
+                    for edge in current_vertex.link_edges:
+                        if edge in edge_pool:
+                            current_edge = edge
+                            edge_pool.remove(current_edge)
+                            follow_edge_loop = True
+                            break
+                current_vertex = current_edge.other_vert(current_vertex)
+                vertices.append(current_vertex.co)
+                internal.addPolygonSpline(toolpath, False, vertices)
+            aux_mesh.free()
+        mesh.free()
+        return {'FINISHED'}
+
+class DiscretizeCurve(bpy.types.Operator):
+    bl_idname = 'curve.add_toolpath_discretize_curve'
+    bl_description = bl_label = 'Discretize Curve'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    step_angle: bpy.props.FloatProperty(name='Resolution', description='Smaller values make curves smoother by adding more vertices', unit='ROTATION', min=math.pi/512, default=math.pi/16)
+    samples: bpy.props.IntProperty(name='Sample Count', description='Number of samples to test per curve segment', min=1, default=128)
+
+    @classmethod
+    def poll(cls, context):
+        return bpy.context.object != None and bpy.context.object.type == 'CURVE'
+
+    def execute(self, context):
+        if bpy.context.object.mode == 'EDIT':
+            splines = internal.getSelectedSplines(True, False)
+        else:
+            splines = bpy.context.object.data.splines
+
+        if len(splines) == 0:
+            self.report({'WARNING'}, 'Nothing selected')
+            return {'CANCELLED'}
+
+        if bpy.context.object.mode != 'EDIT':
+            internal.addObject('CURVE', 'Discretized Curve')
+            origin = bpy.context.scene.cursor.location
+        else:
+            origin = Vector((0.0, 0.0, 0.0))
+
+        for spline in splines:
+            if spline.type != 'BEZIER':
+                continue
+            result = internal.discretizeCurve(spline, self.step_angle, self.samples)
+            internal.addPolygonSpline(bpy.context.object, spline.use_cyclic_u, [vertex-origin for vertex in result])
+        return {'FINISHED'}
+
+class Truncate(bpy.types.Operator):
+    bl_idname = 'curve.add_toolpath_truncate'
+    bl_description = bl_label = 'Truncate'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    min_dist: bpy.props.FloatProperty(name='Min Distance', unit='LENGTH', description='Remove vertices which are too close together', min=0.0, default=0.001)
+    z_hop: bpy.props.BoolProperty(name='Z Hop', description='Add movements to the ceiling at trace ends', default=True)
+
+    @classmethod
+    def poll(cls, context):
+        return bpy.context.object != None and bpy.context.object.mode == 'OBJECT'
+
+    def execute(self, context):
+        if bpy.context.object.type != 'EMPTY' or bpy.context.object.empty_display_type != 'CUBE':
+            self.report({'WARNING'}, 'Active object must be an empty of display type cube')
+            return {'CANCELLED'}
+        selection = bpy.context.selected_objects[:]
+        workspace = bpy.context.object
+        aabb = internal.AABB(center=Vector((0.0, 0.0, 0.0)), dimensions=Vector((1.0, 1.0, 1.0))*workspace.empty_display_size)
+        toolpath = internal.addObject('CURVE', 'Truncated Toolpath')
+        for curve in selection:
+            if curve.type == 'CURVE':
+                transform = workspace.matrix_world.inverted()@curve.matrix_world
+                inverse_transform = Matrix.Translation(-toolpath.location)@workspace.matrix_world
+                curve_traces = []
+                for spline in curve.data.splines:
+                    if spline.type == 'POLY':
+                        curve_traces += internal.truncateToFitBox(transform, spline, aabb)
+                for trace in curve_traces:
+                    i = len(trace[0])-1
+                    while i > 1:
+                        if (trace[0][i-1]-trace[0][i]).length < self.min_dist:
+                            trace[0].pop(i-1)
+                            trace[1].pop(i-1)
+                        i -= 1
+                    if self.z_hop:
+                        begin = Vector(trace[0][0])
+                        end = Vector(trace[0][-1])
+                        begin.z = end.z = workspace.empty_display_size
+                        trace[0].insert(0, begin)
+                        trace[1].insert(0, 1.0)
+                        trace[0].append(end)
+                        trace[1].append(1.0)
+                    internal.addPolygonSpline(toolpath, False, [inverse_transform@vertex for vertex in trace[0]], trace[1])
+        return {'FINISHED'}
+
+class RectMacro(bpy.types.Operator):
+    bl_idname = 'curve.add_toolpath_rect_macro'
+    bl_description = bl_label = 'Rect Macro'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    track_count: bpy.props.IntProperty(name='Number Tracks', description='How many tracks', min=1, default=10)
+    stride: bpy.props.FloatProperty(name='Stride', unit='LENGTH', description='Distance to previous track on the way back', min=0.0, default=0.5)
+    pitch: bpy.props.FloatProperty(name='Pitch', unit='LENGTH', description='Distance between two tracks', default=-1.0)
+    length: bpy.props.FloatProperty(name='Length', unit='LENGTH', description='Length of one track', default=10.0)
+    speed: bpy.props.FloatProperty(name='Speed', description='Stored in softbody goal weight', min=0.0, max=1.0, default=0.1)
+
+    def execute(self, context):
+        if not internal.curveObject():
+            internal.addObject('CURVE', 'Rect Toolpath')
+            origin = Vector((0.0, 0.0, 0.0))
+        else:
+            origin = bpy.context.scene.cursor.location
+        stride = math.copysign(self.stride, self.pitch)
+        length = self.length*0.5
+        vertices = []
+        weights = []
+        for i in range(0, self.track_count):
+            shift = i*self.pitch
+            flipped = -1 if (stride == 0 and i%2 == 1) else 1
+            vertices.append(origin+Vector((shift, -length*flipped, 0.0)))
+            weights.append(self.speed)
+            vertices.append(origin+Vector((shift, length*flipped, 0.0)))
+            weights.append(self.speed)
+            if stride != 0:
+                vertices.append(origin+Vector((shift-stride, length, 0.0)))
+                weights.append(self.speed)
+                vertices.append(origin+Vector((shift-stride, -length, 0.0)))
+                weights.append(1)
+        internal.addPolygonSpline(bpy.context.object, False, vertices, weights)
+        return {'FINISHED'}
+
+class DrillMacro(bpy.types.Operator):
+    bl_idname = 'curve.add_toolpath_drill_macro'
+    bl_description = bl_label = 'Drill Macro'
+    bl_options = {'REGISTER', 'UNDO'}
+
+    screw_count: bpy.props.FloatProperty(name='Screw Turns', description='How many screw truns', min=1.0, default=10.0)
+    spiral_count: bpy.props.FloatProperty(name='Spiral Turns', description='How many spiral turns', min=0.0, default=0.0)
+    vertex_count: bpy.props.IntProperty(name='Number Vertices', description = 'How many vertices per screw turn', min=3, default=32)
+    radius: bpy.props.FloatProperty(name='Radius', unit='LENGTH', description='Radius at tool center', min=0.0, default=5.0)
+    pitch: bpy.props.FloatProperty(name='Pitch', unit='LENGTH', description='Distance between two screw turns', min=0.0, default=1.0)
+    speed: bpy.props.FloatProperty(name='Speed', description='Stored in softbody goal weight', min=0.0, max=1.0, default=0.1)
+
+    def execute(self, context):
+        if not internal.curveObject():
+            internal.addObject('CURVE', 'Drill Toolpath')
+            origin = Vector((0.0, 0.0, 0.0))
+        else:
+            origin = bpy.context.scene.cursor.location
+        count = int(self.vertex_count*self.screw_count)
+        height = -count/self.vertex_count*self.pitch
+        vertices = []
+        weights = []
+        def addRadialVertex(param, radius, height):
+            angle = param*math.pi*2
+            vertices.append(origin+Vector((math.sin(angle)*radius, math.cos(angle)*radius, height)))
+            weights.append(self.speed)
+        if self.radius > 0:
+            if self.spiral_count > 0.0:
+                sCount = math.ceil(self.spiral_count*self.vertex_count)
+                for j in range(1, int(self.screw_count)+1):
+                    if j > 1:
+                        vertices.append(origin+Vector((0.0, 0.0, sHeight)))
+                        weights.append(self.speed)
+                    sHeight = max(-j*self.pitch, height)
+                    for i in range(0, sCount+1):
+                        sParam = i/self.vertex_count
+                        addRadialVertex(sParam, i/sCount*self.radius, sHeight)
+                    for i in range(0, self.vertex_count+1):
+                        addRadialVertex(sParam+(count+i)/self.vertex_count, self.radius, sHeight)
+            else:
+                for i in range(0, count):
+                    param = i/self.vertex_count
+                    addRadialVertex(param, self.radius, -param*self.pitch)
+                for i in range(0, self.vertex_count+1):
+                    addRadialVertex((count+i)/self.vertex_count, self.radius, height)
+            weights += [1, 1]
+        else:
+            weights += [self.speed, 1]
+        vertices += [origin+Vector((0.0, 0.0, height)), origin]
+        internal.addPolygonSpline(bpy.context.object, False, vertices, weights)
+        return {'FINISHED'}
+
+operators = [OffsetCurve, SliceMesh, DiscretizeCurve, Truncate, RectMacro, DrillMacro]