import numpy as np
import matplotlib.pyplot as plt
import math
import operator

from classes.TW_Utility import TW_Utility
from classes.TW_Segment import TW_Point
from classes.TW_Segment import TW_Segment
from classes.TW_Linear_regression import TW_Linear_regression

import warnings

class TW_Segment_Regressor:
    def __init__(self, X, Y, thresold):
        self.segments = []
        self.currentSegment = None
        self.X = X
        self.Y = Y
        self.thresold = thresold

    # todo traiter les outlier
    # todo faire varier la taille de la sliding window
    # todo réfléchir à faire varier la di
    def process(self, X = None, Y = None):
        if not X is None:
            self.X = X

        if not Y is None:
            self.Y = Y
            
        X = self.X
        Y = self.Y

        if len(X) == 0:
            raise Exception("Empty data")
        
        if(len(X.shape)) > 2:
            raise Exception("3 dimensionals data not yet supported")
        
        segmentFirstPoint = None
        segmentNextPoint = None
        nextPointRegression = None

        for i in range (0, len(X)):
            if TW_Utility.empty_or_none(segmentFirstPoint):
                segmentFirstPoint = TW_Point(X[i], Y[i])
                currentSegment = TW_Segment(i)
                currentSegment.start = segmentFirstPoint
                currentSegment.addPoint(segmentFirstPoint)
            else:
                segmentNextPoint = TW_Point(X[i], Y[i])
                currentSegment.addPoint(segmentNextPoint)
                    
            nextPointRegression = None
            if not TW_Utility.empty_or_none(segmentNextPoint):
                #nextPointRegression = TW_2D_Linear_regression(np.array([segmentFirstPoint[0], segmentNextPoint[0]]),np.array([segmentFirstPoint[1], segmentNextPoint[1]]))
                nextPointRegression = TW_Linear_regression(np.array(currentSegment.getX()),np.array(currentSegment.getY()))
                nextPointRegression.process()
                currentSegment.regression = nextPointRegression

            if not TW_Utility.empty_or_none(segmentNextPoint) and len(X) > i +1:
                segmentLastPoint = TW_Point(X[i + 1], Y[i + 1])
                predictedLastPoint = nextPointRegression.predict([segmentLastPoint.X])

                diff = (abs(predictedLastPoint - segmentLastPoint.Y)  * 100) /  segmentLastPoint.Y

                if diff > self.thresold:
                    self.segments.append(currentSegment)
                    segmentFirstPoint = TW_Point(X[i], Y[i])
                    segmentNextPoint = None

                    currentSegment = TW_Segment(i)
                    currentSegment.start = segmentFirstPoint
                    currentSegment.addPoint(segmentFirstPoint)
            elif len(X) == i + 1:
                self.segments.append(currentSegment)

        return self.segments
    
    '''
    def _evaluateDistanceBetweenCordinates(self):
        points = None

        for i in range(0,len(self.segments) -1):
            currentSegment = self.segments[i]
            currentSegmentFirstPoint = currentSegment.points[0].getPoint()
            currentSegmentLastPoint = currentSegment.points[len(currentSegment.points) -1].getPoint()

            nextSegment = self.segments[i + 1]
            nextSegmentFirstPoint = nextSegment.points[0].getPoint()
            nextSegmentLastPoint = nextSegment.points[len(nextSegment.points) -1].getPoint()

            nbrColumns = TW_Utility.np_number_columns(currentSegmentFirstPoint)
            try:
                if nbrColumns == 2:
                    points = self.intersect_segments_2d(
                        currentSegmentFirstPoint, 
                        currentSegmentLastPoint,
                        nextSegmentFirstPoint,
                        nextSegmentLastPoint
                    )
                elif nbrColumns == 3:
                    points = self.intersect_segments_3d(
                        currentSegmentFirstPoint, 
                        currentSegmentLastPoint,
                        nextSegmentFirstPoint,
                        nextSegmentLastPoint
                    )
                else:
                    warnings.warn("Prediction on higher degrees than 3 might give wrong results")
                    points = self.intersect_segments_nd(
                        currentSegmentFirstPoint, 
                        currentSegmentLastPoint,
                        nextSegmentFirstPoint,
                        nextSegmentLastPoint
                    )
            except:
                points = None

        return currentSegment, nextSegment, points
    '''
    
    @staticmethod
    def getSegmentFittingCoordinates(segments, X):
        matchedSegment = None

        if TW_Utility.is_array(X) and len(X) > 1:
            raise Exception("X must be a unique feature with 1 Dimension or an array of one feature with N Dimensions")
            
        if not TW_Utility.is_array(X):
            X = [X]

        for i in range(0,len(segments)):
            segmentX = segments[i].getX()
            segmentXMin = segmentX[0]
            segmentXMax = segmentX[len(segmentX) -1]

            if len(X) < 2 and not TW_Utility.is_array(X[0]):
                if X[0] >= segmentXMin and X[0] <= segmentXMax:
                    matchedSegment = segments[i]
                    break
            else:
                valueToCompare = X[0]
                nbrDimensionsToConsider = len(valueToCompare)
                nbrDimensionsMatched = 0

                for j in range (0,len(valueToCompare)):    
                    if valueToCompare[j] >= segmentXMin[j] and valueToCompare[j] <= segmentXMax[j]:
                        nbrDimensionsMatched = nbrDimensionsMatched + 1

                if nbrDimensionsMatched == nbrDimensionsToConsider:
                    matchedSegment = segments[i]
                    break
        
        return matchedSegment
        
    def _getPointForDimension(self, nbrColumns, X):
        return X[nbrColumns -1]
    
        '''
        firstPoint = None
        for i in range(0,len(X)):
            firstPoint = X[0] if firstPoint is not None else firstPoint[0]

        return firstPoint
        '''

    def _isCoordinateBeforeFirstSegment(self, firstSegment, X):
        coordinateBeforeFirstSegment = False

        if TW_Utility.is_array(X) and len(X) > 1:
            raise Exception("X must be a unique feature with 1 Dimension or an array of one feature with N Dimensions")

        if not TW_Utility.is_array(X):
            X = [X]

        if len(X) < 2 and not TW_Utility.is_array(X[0]):
            if X[0] < firstSegment.getX()[0]:
                coordinateBeforeFirstSegment = True
        else:
            nbrColumns = TW_Utility.np_number_columns(X)
            valueToCompare = X[0]
            point = firstSegment.getX()[0]

            euclidianDistanceBetweenOriginAndFirstSegmentPoint = TW_Utility._euclidianDistance(np.zeros(nbrColumns), point)
            euclidianDistanceBetweenOriginAndCoordinate = TW_Utility._euclidianDistance(np.zeros(nbrColumns), valueToCompare)
            coordinateBeforeFirstSegment = euclidianDistanceBetweenOriginAndCoordinate < euclidianDistanceBetweenOriginAndFirstSegmentPoint

        return coordinateBeforeFirstSegment
    

    
    def _isCoordinateAfterLastSegment(self, lastSegment, X):
        coordinateAfterLastSegment = False

        lastSegmentPoints = lastSegment.getX()
        lastSegmentFirstPoint = lastSegmentPoints[0]
        lastSegmentLastPoint = lastSegmentPoints[len(lastSegmentPoints) -1]

        if TW_Utility.is_array(X) and len(X) > 1:
            raise Exception("X must be a unique feature with 1 Dimension or an array of one feature with N Dimensions")

        if not TW_Utility.is_array(X):
            X = [X]

        if len(X) < 2 and not TW_Utility.is_array(X[0]):
            if X[0] > lastSegmentLastPoint:
                coordinateAfterLastSegment = True
        else:
            nbrColumns = TW_Utility.np_number_columns(X)
            valueToCompare = X[0]

            origin = np.zeros(nbrColumns)
            lastSegmentFirstPoint = lastSegment.points[0].getX()
            lastSegmentLastPoint = lastSegment.points[len(lastSegment.points) -1].getX()


            euclidianDistanceBetweenOriginAndLastSegmentLastPoint = TW_Utility._euclidianDistance(origin, lastSegmentLastPoint)
            euclidianDistanceBetweenOriginAndPoint = TW_Utility._euclidianDistance(origin, valueToCompare)
            coordinateAfterLastSegment = euclidianDistanceBetweenOriginAndPoint > euclidianDistanceBetweenOriginAndLastSegmentLastPoint

        return coordinateAfterLastSegment

    def predict(self, X):
        if TW_Utility.is_array(X) and len(X) > 1:
            raise Exception("X must be a unique feature with 1 Dimension or an array of one feature with N Dimensions")

        if not TW_Utility.is_array(X):
            X = [X]

        predictedValue = None

        matchedSegment = TW_Segment_Regressor.getSegmentFittingCoordinates(self.segments, X)
        if matchedSegment is None:
            firstSegment = self.segments[0]
            lastSegment = self.segments[len(self.segments) -1]

            if self._isCoordinateBeforeFirstSegment(firstSegment, X):
                matchedSegment = firstSegment
            elif self._isCoordinateAfterLastSegment(lastSegment, X):
                matchedSegment = lastSegment
            else:
                euclidianDistances = {}
                for i in range(0,len(self.segments) -1):
                    currentSegment = self.segments[i]
                    currentSegmentPoints = currentSegment.getX()
                    currentSegmentFirstPoint = currentSegmentPoints[0]
                    currentSegmentLastPoint = currentSegmentPoints[len(currentSegmentPoints) -1]

                    nextSegment = self.segments[i + 1]
                    nextSegmentPoints = nextSegment.getX()
                    nextSegmentFirstPoint = nextSegmentPoints[0]
                    nextSegmentLastPoint = nextSegmentPoints[len(nextSegmentPoints) -1]

                    if len(X) < 2 and not TW_Utility.is_array(X[0]):
                        euclidianDistances[str(i) + '_currentSegmentStart'] = [
                            TW_Utility._euclidianDistance(X, currentSegmentFirstPoint),
                            currentSegment
                        ]
                        euclidianDistances[str(i) + '_currentSegmentEnd'] = [
                            TW_Utility._euclidianDistance(X, currentSegmentLastPoint),
                            currentSegment
                        ]
                        euclidianDistances[str(i) + '_nextSegmentStart'] = [
                            TW_Utility._euclidianDistance(X, nextSegmentFirstPoint),
                            nextSegment
                        ]
                        euclidianDistances[str(i) + '_nextSegmentEnd'] = [
                            TW_Utility._euclidianDistance(X, nextSegmentLastPoint),
                            nextSegment
                        ]
                    else:
                        valueToCompare = X[0]
                        euclidianDistances[str(i) + '_currentSegmentStart'] = [
                            TW_Utility._euclidianDistance(valueToCompare, currentSegmentFirstPoint),
                            currentSegment
                        ]
                        euclidianDistances[str(i) + '_currentSegmentEnd'] = [
                            TW_Utility._euclidianDistance(valueToCompare, currentSegmentLastPoint),
                            currentSegment
                        ]
                        euclidianDistances[str(i) + '_nextSegmentStart'] = [
                            TW_Utility._euclidianDistance(valueToCompare, nextSegmentFirstPoint),
                            nextSegment
                        ]
                        euclidianDistances[str(i) + '_nextSegmentEnd'] = [
                            TW_Utility._euclidianDistance(valueToCompare, nextSegmentLastPoint),
                            nextSegment
                        ]
                    
                euclidianDistanceMinRangeErrors = dict(sorted(euclidianDistances.items(), key=lambda item:item[1][0]))
                matchedSegment = euclidianDistances[next(iter(euclidianDistanceMinRangeErrors))][1]

        if matchedSegment is not None:
            predictedValue = matchedSegment.regression.predict(X)

        return matchedSegment, predictedValue
    
    '''
    def intersect_segments_2d_mine(p1, p2, p3, p4):
        Calcule le point d'intersection de deux segments en N dimensions, s'il existe.

        Args:
            p1 (np.array): Premier point du premier segment.
            p2 (np.array): Deuxième point du premier segment.
            p3 (np.array): Premier point du deuxième segment.
            p4 (np.array): Deuxième point du deuxième segment.

        Returns:
            np.array or None: Le point d'intersection si les segments se coupent,
                            sinon None.
        
        # Convertir les points en tableaux numpy pour faciliter les opérations
        p1 = np.array(p1, dtype=float)
        p2 = np.array(p2, dtype=float)
        p3 = np.array(p3, dtype=float)
        p4 = np.array(p4, dtype=float)

        # Vérifier que tous les points ont la même dimension
        if not (p1.shape == p2.shape == p3.shape == p4.shape):
            raise ValueError("Tous les points doivent avoir la même dimension.")
        
        # Vecteurs directeurs des droites
        v1 = p2 - p1  # Vecteur pour le premier segment
        v2 = p4 - p3  # Vecteur pour le deuxième segment

        #si v1 et v2 sont parallèle ou collinéaire  il existe un réel K tel que kV1 = V2
        k = v2[0] / v1[0]
        if v1[1] * k == v2[1]:
            return None

        #Y = AX + b
        Line1EquationA = p1[1] - p2[1] / (p1[0] - p2[0])
        Line1EquationB = p1[1] - (p1[0] * Line1EquationB[0])
        Line2EquationA = p2[1] - p2[1] / (p2[0] - p2[0])
        Line2EquationB = p2[1] - (p2[0] * Line2EquationB[0])

        #point d'intersection x se trouve à -b1 + b2 / (a1 -a2)
        intersectionX = (-Line1EquationB + Line2EquationB) / (Line1EquationA - Line2EquationA)
        intersectionY = (Line1EquationA * intersectionX) + Line1EquationB 

        return [intersectionX, intersectionY]
    '''
    

    def intersect_segments_2d(p1, p2, p3, p4, epsilon=1e-9):
        """
        Calcule le point d'intersection de deux segments 2D.
        Renvoie le point d'intersection [x, y] si une intersection unique existe.
        Renvoie None si les segments sont parallèles, colinéaires (non chevauchants) ou ne se coupent pas.
        Renvoie une liste de deux points si les segments sont colinéaires et se chevauchent (le segment d'intersection).
        """
        p1 = np.array(p1, dtype=float)
        p2 = np.array(p2, dtype=float)
        p3 = np.array(p3, dtype=float)
        p4 = np.array(p4, dtype=float)

        # Vecteurs directeurs des droites
        v1 = p2 - p1 # Direction du segment 1
        v2 = p4 - p3 # Direction du segment 2
        
        # Différence entre les points de départ
        dp = p3 - p1

        # Calcul du déterminant (produit croisé 2D) des vecteurs directeurs
        # Cela correspond à (v1_x * v2_y - v1_y * v2_x)
        det = v1[0] * v2[1] - v1[1] * v2[0]

        # --- Cas où les droites sont parallèles ou colinéaires ---
        if abs(det) < epsilon: # Les droites sont parallèles ou colinéaires
            # Vérifier si elles sont colinéaires (et potentiellement se chevauchent)
            # Si dp est parallèle à v1, alors elles sont colinéaires
            det_collinear = v1[0] * dp[1] - v1[1] * dp[0]
            
            if abs(det_collinear) < epsilon: # Les droites sont colinéaires
                # Elles sont colinéaires. Il faut vérifier s'il y a chevauchement des segments.
                # Convertir les points sur une dimension (ex: axe X si v1_x est grand, sinon axe Y)
                # ou projeter les points sur le vecteur directeur
                
                # Utiliser les projections sur les vecteurs directeurs pour vérifier le chevauchement
                # Projection des points du seg2 sur le seg1 (en prenant p1 comme origine)
                # t0 = (p3 - p1) . v1 / (v1 . v1)
                # t1 = (p4 - p1) . v1 / (v1 . v1)
                
                len_sq_v1 = np.dot(v1, v1)
                if len_sq_v1 < epsilon: # Segment 1 est un point
                    return p1 if np.linalg.norm(p3 - p1) < epsilon and np.linalg.norm(p4 - p1) < epsilon else None
                
                t_min_s1 = 0
                t_max_s1 = 1
                
                s_proj_p3 = np.dot(p3 - p1, v1) / len_sq_v1
                s_proj_p4 = np.dot(p4 - p1, v1) / len_sq_v1
                
                # Les intervalles de t pour le chevauchement
                # Si p3 < p4, alors min_s2 = s_proj_p3, max_s2 = s_proj_p4
                # Sinon, min_s2 = s_proj_p4, max_s2 = s_proj_p3
                min_s2 = min(s_proj_p3, s_proj_p4)
                max_s2 = max(s_proj_p3, s_proj_p4)

                # Calculer l'intervalle de chevauchement sur la ligne paramétrique de S1
                overlap_min = max(t_min_s1, min_s2)
                # vaut 1 si ||P3 - P1|| > ||P1 p2||


                overlap_max = min(t_max_s1, max_s2)
                # >= 1

                #si overlap_max > 1 alors il y'a chevauchement
                if overlap_min <= overlap_max + epsilon: # Il y a chevauchement
                    # Retourner le segment d'intersection (ses deux points d'extrémité)
                    # Il faut être prudent pour les cas de chevauchement ponctuel
                    if abs(overlap_min - overlap_max) < epsilon:
                        return [p1 + overlap_min * v1] # Un seul point d'intersection
                    else:
                        return [p1 + overlap_min * v1, p1 + overlap_max * v1]
                else:
                    return None # Colinéaires mais pas de chevauchement
            else:
                return None # Parallèles et distincts
        

        # --- Cas où les droites ne sont PAS parallèles (elles se coupent à un point unique) ---
        # Calcul des paramètres t et u
        t = (dp[0] * v2[1] - dp[1] * v2[0]) / det
        u = (dp[0] * v1[1] - dp[1] * v1[0]) / det

        # Vérifier si le point d'intersection se trouve sur les deux segments
        # t doit être entre 0 et 1 (inclus) pour le segment 1
        # u doit être entre 0 et 1 (inclus) pour le segment 2
        if 0 - epsilon <= t <= 1 + epsilon and 0 - epsilon <= u <= 1 + epsilon:
            # Le point d'intersection est sur les deux segments
            intersection_point = p1 + t * v1
            return intersection_point
        else:
            # Le point d'intersection est en dehors d'au moins un segment
            return None

    '''
    def intersect_segments_3dmine(gridX1, gridY1, gridZ1, gridX2, gridY2, gridZ2):

        #etape 1 trouver l equation du plan1
        #l'equation du plan 1 satisfait la matrice (X,Y,Z,1) * (A,B,C,D) = (0,0,0,0)

        #etape 2 trouver un vecteur normal à ce plan (ie le produit vectoriel de deux vecteurs non colinéaire)

        x1_values = gridX1[:n]
        y1_values = gridY1[:n]
        z1_values = gridZ1[:n]
        ones_column = np.ones(n)
        zero_column = np.zeros(n)
        firstPlanMatrix = np.column_stack((x1_values, y1_values, z1_values, ones_column))
        firstPlanParameters = np.linalg.solve(firstPlanMatrix, zero_column)
        firstPlanDirectorsVector = TW_Utility.get_vector_directors(firstPlanParameters)
        firstPlanNormalVector = np.cross(firstPlanDirectorsVector[0], firstPlanDirectorsVector[1])

        x2_values = gridX2[:n]
        y2_values = gridY2[:n]
        z2_values = gridZ2[:n]
        ones_column = np.ones(n)
        zero_column = np.zeros(n)
        secondPlanMatrix = np.column_stack((x2_values, y2_values, z2_values, ones_column))
        secondPlanParameters = np.linalg.solve(secondPlanMatrix, zero_column)
        secondPlanDirectorsVector = TW_Utility.get_vector_directors(secondPlanParameters)
        secondPlanNormalVector = np.cross(secondPlanDirectorsVector[0], secondPlanDirectorsVector[1])

        #etape 3 trouver un vecteur directeur sur la ligne d'intersection de ce plan 
        normalVector = np.cross(firstPlanNormalVector, secondPlanNormalVector)
        return normalVector
    '''
            

    def intersect_segments_3d(p1, p2, p3, p4):
        """
        Calcule le point d'intersection de deux segments 3D, s'il existe.

        Args:
            p1 (np.array): Premier point du premier segment [x1, y1, z1].
            p2 (np.array): Deuxième point du premier segment [x2, y2, z2].
            p3 (np.array): Premier point du deuxième segment [x3, y3, z3].
            p4 (np.array): Deuxième point du deuxième segment [x4, y4, z4].

        Returns:
            np.array or None: Le point d'intersection si les segments se coupent,
                            sinon None.
        """
        # Convertir les points en tableaux numpy pour faciliter les opérations
        p1 = np.array(p1)
        p2 = np.array(p2)
        p3 = np.array(p3)
        p4 = np.array(p4)

        # Vecteurs directeurs des droites
        v1 = p2 - p1  # Vecteur pour le premier segment
        v2 = p4 - p3  # Vecteur pour le deuxième segment

        # 1. Vérifier le parallélisme des droites
        # Si le produit vectoriel est nul ou très proche de zéro (en raison de la précision flottante)
        cross_product_v = np.cross(v1, v2)
        if np.linalg.norm(cross_product_v) < 1e-9:  # Presque parallèle
            # Les droites sont parallèles (ou colinéaires/identiques)
            # Vérifier si elles sont colinéaires (un point de l'une est sur l'autre)
            # Si P1-P3 est colinéaire à V1, alors elles sont colinéaires
            v_p1_p3 = p1 - p3
            collinear_check = np.cross(v_p1_p3, v1)
            if np.linalg.norm(collinear_check) < 1e-9:
                # Les droites sont colinéaires (identiques ou juste sur la même ligne)
                # Nous devons vérifier le chevauchement des segments sur cette ligne
                # C'est un problème d'intersection d'intervalles 1D.
                # Projection des points sur l'axe du vecteur v1 pour simplifier
                # (ou v2, peu importe car ils sont colinéaires)
                t0 = 0.0
                t1 = np.dot(p2 - p1, v1) / np.dot(v1, v1) if np.dot(v1, v1) != 0 else 0.0
                t_p3_proj = np.dot(p3 - p1, v1) / np.dot(v1, v1) if np.dot(v1, v1) != 0 else 0.0
                t_p4_proj = np.dot(p4 - p1, v1) / np.dot(v1, v1) if np.dot(v1, v1) != 0 else 0.0

                # Assurer que les intervalles sont triés
                t_seg1_min, t_seg1_max = sorted([t0, t1])
                t_seg2_min, t_seg2_max = sorted([t_p3_proj, t_p4_proj])

                # Calcul de l'intersection des intervalles
                overlap_min = max(t_seg1_min, t_seg2_min)
                overlap_max = min(t_seg1_max, t_seg2_max)

                if overlap_min <= overlap_max:
                    # Il y a un chevauchement. S'ils sont exactement au même point, on peut le retourner.
                    # Pour cet exemple, nous retournons le point de début de chevauchement s'il y en a un seul
                    # ou None si le chevauchement est un segment (pour simplifier).
                    # Dans un cas réel, vous pourriez vouloir retourner un segment pour l'intersection.
                    if overlap_min == overlap_max: # Les segments se touchent à un point
                        return p1 + overlap_min * v1
                    else: # Les segments se chevauchent sur une longueur
                        # Pour la simplicité de cette fonction, on considère que l'intersection est un point unique
                        # Si c'est un segment, la fonction retourne None (pas de point unique d'intersection)
                        # ou vous pouvez retourner (p1 + overlap_min * v1, p1 + overlap_max * v1)
                        # Je retourne None ici pour éviter de changer le type de retour.
                        return [p1 + overlap_min * v1, p1 + overlap_max * v1]
                else:
                    return None # Pas de chevauchement pour les segments colinéaires
            else:
                return None # Droites parallèles et distinctes, pas d'intersection

        # 2. Vérifier si les droites sont gauches (ne se coupent pas et ne sont pas parallèles)
        # Si le produit mixte n'est pas nul, les droites sont gauches
        # (p3 - p1) . (v1 x v2)
        mixed_product = np.dot((p3 - p1), cross_product_v)
        if abs(mixed_product) > 1e-9: # Utiliser une tolérance pour les nombres flottants
            return None # Droites gauches, pas d'intersection

        # 3. Calculer les paramètres t et u du point d'intersection
        # Système d'équations: P1 + t*V1 = P3 + u*V2
        # Soit: t*V1 - u*V2 = P3 - P1

        # On utilise les deux premières équations (x et y) pour trouver t et u,
        # puis on vérifie avec la troisième (z).
        # Matrice des coefficients pour [x, y]
        A = np.array([
            [v1[0], -v2[0]],
            [v1[1], -v2[1]]
        ])
        # Vecteur des constantes pour [x, y]
        B = np.array([
            p3[0] - p1[0],
            p3[1] - p1[1]
        ])

        # Résoudre le système A * [t, u] = B
        # Gérer le cas où le déterminant est proche de zéro (e.g., v1[0]*(-v2[1]) - (-v2[0])*v1[1] est nul)
        # ce qui signifie que les projections 2D sont parallèles, mais elles doivent être coplanaires en 3D
        # si elles ne sont pas parallèles en 3D.
        # On peut utiliser np.linalg.solve qui gère les cas singuliers ou pseudo-inverse si nécessaire.
        try:
            t_u = np.linalg.solve(A, B)
            t = t_u[0]
            u = t_u[1]
        except np.linalg.LinAlgError:
            # Cela peut arriver si la projection 2D est singulière, mais les droites sont coplanaires en 3D.
            # Dans ce cas, nous devons choisir d'autres axes.
            # Une méthode plus robuste utilise la pseudo-inverse ou un algorithme pour des systèmes surdéterminés.
            # Ou on peut choisir un autre sous-système 2x2.
            # Par exemple, si v1[0] et v2[0] sont petits, utilisez y et z.
            A_yz = np.array([
                [v1[1], -v2[1]],
                [v1[2], -v2[2]]
            ])
            B_yz = np.array([
                p3[1] - p1[1],
                p3[2] - p1[2]
            ])
            try:
                t_u = np.linalg.solve(A_yz, B_yz)
                t = t_u[0]
                u = t_u[1]
            except np.linalg.LinAlgError:
                # Si même yz est singulier, c'est un cas très dégénéré ou un bug.
                return None


        # Vérifier que les valeurs de t et u satisfont la troisième équation (pour la robustesse)
        # p1[2] + t * v1[2] doit être proche de p3[2] + u * v2[2]
        if not np.isclose(p1[2] + t * v1[2], p3[2] + u * v2[2]):
            return None # Incohérence, les droites ne se coupent pas ou problème de précision

        # 4. Vérifier si le point d'intersection se trouve sur les deux segments
        if 0 <= t <= 1 and 0 <= u <= 1:
            # Le point d'intersection est sur les deux segments
            intersection_point = p1 + t * v1
            return intersection_point
        else:
            # Le point d'intersection des droites n'est pas sur les segments
            return None
    
    def intersect_segments_nd(p1, p2, p3, p4):
        """
        Calcule le point d'intersection de deux segments en N dimensions, s'il existe.

        Args:
            p1 (np.array): Premier point du premier segment.
            p2 (np.array): Deuxième point du premier segment.
            p3 (np.array): Premier point du deuxième segment.
            p4 (np.array): Deuxième point du deuxième segment.

        Returns:
            np.array or None: Le point d'intersection si les segments se coupent,
                            sinon None.
        """
        # Convertir les points en tableaux numpy pour faciliter les opérations
        p1 = np.array(p1, dtype=float)
        p2 = np.array(p2, dtype=float)
        p3 = np.array(p3, dtype=float)
        p4 = np.array(p4, dtype=float)

        # Vérifier que tous les points ont la même dimension
        if not (p1.shape == p2.shape == p3.shape == p4.shape):
            raise ValueError("Tous les points doivent avoir la même dimension.")
        
        n_dim = p1.shape[0]

        # Vecteurs directeurs des droites
        v1 = p2 - p1  # Vecteur pour le premier segment
        v2 = p4 - p3  # Vecteur pour le deuxième segment

        # Construire la matrice du système linéaire A * [t, u] = B
        # L'équation est: p1 + t*v1 = p3 + u*v2
        # Réécrite: t*v1 - u*v2 = p3 - p1
        # Matrice A: [v1 | -v2] (où v1 et -v2 sont des colonnes)
        # Vecteur B: p3 - p1

        A = np.column_stack((v1, -v2))
        B = p3 - p1

        # Résoudre le système linéaire A @ [t, u] = B
        # Utiliser np.linalg.lstsq pour une solution plus robuste (moindres carrés)
        # qui peut gérer les systèmes surdéterminés (plus d'équations que d'inconnues, N_dim > 2)
        try:
            # lstsq retourne (solution, résidus, rang, valeurs singulières)
            # Nous sommes intéressés par la solution [t, u]
            t_u, residuals, rank, s = np.linalg.lstsq(A, B, rcond=None)
            t = t_u[0]
            u = t_u[1]
        except np.linalg.LinAlgError:
            # Cela indique un problème grave avec la matrice A, comme des vecteurs colinéaires
            # qui n'ont pas été gérés au préalable.
            return None

        # Vérifier le rang pour déterminer si les droites s'intersectent
        # Les droites s'intersectent si le système a une solution exacte
        # C'est-à-dire si le vecteur B est dans l'espace colonne de A.
        # Ceci est équivalent à dire que le résidu est proche de zéro.
        
        # La norme des résidus doit être très petite pour qu'il y ait une intersection exacte
        # (ce qui signifie que les droites sont coplanaires et se coupent)
        if residuals.size > 0 and np.linalg.norm(residuals) > 1e-9:
            return None # Les droites ne s'intersectent pas (gauches)

        # Vérification du parallélisme si le système a une solution exacte mais le rang est bas
        # C'est un cas spécial où v1 et v2 sont parallèles (colinéaires)
        # Si le rang est 1, cela signifie que v1 et v2 sont colinéaires.
        if rank < 2: # Cela signifie que v1 et v2 sont colinéaires
            # Les droites sont parallèles (ou colinéaires/identiques)
            # Vérifier si elles sont colinéaires (un point de l'une est sur l'autre)
            # Si P1-P3 est colinéaire à V1 (ou V2)
            v_p1_p3 = p1 - p3
            # Utiliser la norme du produit vectoriel pour 3D, ou simplement le produit scalaire
            # ou le test de colinéarité généralisé pour N-D (vérifier si les vecteurs sont multiples l'un de l'autre)
            # np.linalg.norm(np.cross(v_p1_p3, v1)) si n_dim == 3
            
            # Pour N-D: vérifier si (p1-p3) est colinéaire à v1.
            # Cela peut se faire en vérifiant si rank([v1, p1-p3]) == 1
            test_collinear_matrix = np.column_stack((v1, v_p1_p3))
            if np.linalg.matrix_rank(test_collinear_matrix, tol=1e-9) <= 1:
                # Les droites sont colinéaires (identiques ou juste sur la même ligne)
                # Nous devons vérifier le chevauchement des segments sur cette ligne.
                # Projection des points sur l'axe du vecteur v1
                
                # Pour éviter la division par zéro si v1 est nul (P1 == P2)
                v1_norm_sq = np.dot(v1, v1)
                if v1_norm_sq < 1e-9: # Segment 1 est un point
                    if np.linalg.norm(p3 - p1) < 1e-9 and np.linalg.norm(p4 - p1) < 1e-9: # Si P3 et P4 sont aussi P1
                        return p1 # Les deux segments sont le même point
                    elif np.linalg.norm(p3 - p1) < 1e-9: # P1 est P3
                        if 0 <= np.dot(p4 - p3, v2) / np.dot(v2, v2) <= 1 and np.linalg.norm(v2) > 1e-9:
                            return p1 # P1 est sur le segment P3P4
                    elif np.linalg.norm(p4 - p1) < 1e-9: # P1 est P4
                        if 0 <= np.dot(p3 - p4, -v2) / np.dot(-v2, -v2) <= 1 and np.linalg.norm(v2) > 1e-9:
                            return p1 # P1 est sur le segment P3P4
                    return None # Segment 1 est un point mais pas d'intersection
                
                t0 = 0.0
                t1 = np.dot(p2 - p1, v1) / v1_norm_sq
                t_p3_proj = np.dot(p3 - p1, v1) / v1_norm_sq
                t_p4_proj = np.dot(p4 - p1, v1) / v1_norm_sq

                t_seg1_min, t_seg1_max = sorted([t0, t1])
                t_seg2_min, t_seg2_max = sorted([t_p3_proj, t_p4_proj])

                overlap_min = max(t_seg1_min, t_seg2_min)
                overlap_max = min(t_seg1_max, t_seg2_max)

                if overlap_min <= overlap_max + 1e-9: # Tolérance pour les flottants
                    if abs(overlap_min - overlap_max) < 1e-9: # Les segments se touchent à un point
                        return p1 + overlap_min * v1
                    else: # Les segments se chevauchent sur une longueur
                        return None # Retourne None si l'intersection est un segment, pas un point unique
                else:
                    return None # Pas de chevauchement pour les segments colinéaires
            else:
                return None # Droites parallèles et distinctes, pas d'intersection
                
        # Si le système a été résolu avec succès et que les résidus sont proches de zéro,
        # les droites s'intersectent.
        # Vérifier si le point d'intersection se trouve sur les deux segments
        if 0 - 1e-9 <= t <= 1 + 1e-9 and 0 - 1e-9 <= u <= 1 + 1e-9: # Ajouter une tolérance pour les bornes
            intersection_point = p1 + t * v1
            return intersection_point
        else:
            # Le point d'intersection des droites n'est pas sur les segments
            return None
    
    def totalError(self):
        totalError = 0
        for i, segment in enumerate(self.segments):
            segmentError = self.segments[i].regression.lossFunction(self.segments[i].getX(), self.segments[i].getY())
            totalError = totalError + segmentError

        return totalError / len(self.segments)
    
    '''
    def intersect_segments_2d_mine(p1, p2, p3, p4):
        Calcule le point d'intersection de deux segments en N dimensions, s'il existe.

        Args:
            p1 (np.array): Premier point du premier segment.
            p2 (np.array): Deuxième point du premier segment.
            p3 (np.array): Premier point du deuxième segment.
            p4 (np.array): Deuxième point du deuxième segment.

        Returns:
            np.array or None: Le point d'intersection si les segments se coupent,
                            sinon None.
        
        # Convertir les points en tableaux numpy pour faciliter les opérations
        p1 = np.array(p1, dtype=float)
        p2 = np.array(p2, dtype=float)
        p3 = np.array(p3, dtype=float)
        p4 = np.array(p4, dtype=float)

        # Vérifier que tous les points ont la même dimension
        if not (p1.shape == p2.shape == p3.shape == p4.shape):
            raise ValueError("Tous les points doivent avoir la même dimension.")
        
        # Vecteurs directeurs des droites
        v1 = p2 - p1  # Vecteur pour le premier segment
        v2 = p4 - p3  # Vecteur pour le deuxième segment

        #si v1 et v2 sont parallèle ou collinéaire  il existe un réel K tel que kV1 = V2
        k = v2[0] / v1[0]
        if v1[1] * k == v2[1]:
            return None

        #Y = AX + b
        Line1EquationA = p1[1] - p2[1] / (p1[0] - p2[0])
        Line1EquationB = p1[1] - (p1[0] * Line1EquationB[0])
        Line2EquationA = p2[1] - p2[1] / (p2[0] - p2[0])
        Line2EquationB = p2[1] - (p2[0] * Line2EquationB[0])

        #point d'intersection x se trouve à -b1 + b2 / (a1 -a2)
        intersectionX = (-Line1EquationB + Line2EquationB) / (Line1EquationA - Line2EquationA)
        intersectionY = (Line1EquationA * intersectionX) + Line1EquationB 

        return [intersectionX, intersectionY]
    '''
    

    def intersect_segments_2d(p1, p2, p3, p4, epsilon=1e-9):
        """
        Calcule le point d'intersection de deux segments 2D.
        Renvoie le point d'intersection [x, y] si une intersection unique existe.
        Renvoie None si les segments sont parallèles, colinéaires (non chevauchants) ou ne se coupent pas.
        Renvoie une liste de deux points si les segments sont colinéaires et se chevauchent (le segment d'intersection).
        """
        p1 = np.array(p1, dtype=float)
        p2 = np.array(p2, dtype=float)
        p3 = np.array(p3, dtype=float)
        p4 = np.array(p4, dtype=float)

        # Vecteurs directeurs des droites
        v1 = p2 - p1 # Direction du segment 1
        v2 = p4 - p3 # Direction du segment 2
        
        # Différence entre les points de départ
        dp = p3 - p1

        # Calcul du déterminant (produit croisé 2D) des vecteurs directeurs
        # Cela correspond à (v1_x * v2_y - v1_y * v2_x)
        det = v1[0] * v2[1] - v1[1] * v2[0]

        # --- Cas où les droites sont parallèles ou colinéaires ---
        if abs(det) < epsilon: # Les droites sont parallèles ou colinéaires
            # Vérifier si elles sont colinéaires (et potentiellement se chevauchent)
            # Si dp est parallèle à v1, alors elles sont colinéaires
            det_collinear = v1[0] * dp[1] - v1[1] * dp[0]
            
            if abs(det_collinear) < epsilon: # Les droites sont colinéaires
                # Elles sont colinéaires. Il faut vérifier s'il y a chevauchement des segments.
                # Convertir les points sur une dimension (ex: axe X si v1_x est grand, sinon axe Y)
                # ou projeter les points sur le vecteur directeur
                
                # Utiliser les projections sur les vecteurs directeurs pour vérifier le chevauchement
                # Projection des points du seg2 sur le seg1 (en prenant p1 comme origine)
                # t0 = (p3 - p1) . v1 / (v1 . v1)
                # t1 = (p4 - p1) . v1 / (v1 . v1)
                
                len_sq_v1 = np.dot(v1, v1)
                if len_sq_v1 < epsilon: # Segment 1 est un point
                    return p1 if np.linalg.norm(p3 - p1) < epsilon and np.linalg.norm(p4 - p1) < epsilon else None
                
                t_min_s1 = 0
                t_max_s1 = 1
                
                s_proj_p3 = np.dot(p3 - p1, v1) / len_sq_v1
                s_proj_p4 = np.dot(p4 - p1, v1) / len_sq_v1
                
                # Les intervalles de t pour le chevauchement
                # Si p3 < p4, alors min_s2 = s_proj_p3, max_s2 = s_proj_p4
                # Sinon, min_s2 = s_proj_p4, max_s2 = s_proj_p3
                min_s2 = min(s_proj_p3, s_proj_p4)
                max_s2 = max(s_proj_p3, s_proj_p4)

                # Calculer l'intervalle de chevauchement sur la ligne paramétrique de S1
                overlap_min = max(t_min_s1, min_s2)
                # vaut 1 si ||P3 - P1|| > ||P1 p2||


                overlap_max = min(t_max_s1, max_s2)
                # >= 1

                #si overlap_max > 1 alors il y'a chevauchement
                if overlap_min <= overlap_max + epsilon: # Il y a chevauchement
                    # Retourner le segment d'intersection (ses deux points d'extrémité)
                    # Il faut être prudent pour les cas de chevauchement ponctuel
                    if abs(overlap_min - overlap_max) < epsilon:
                        return [p1 + overlap_min * v1] # Un seul point d'intersection
                    else:
                        return [p1 + overlap_min * v1, p1 + overlap_max * v1]
                else:
                    return None # Colinéaires mais pas de chevauchement
            else:
                return None # Parallèles et distincts
        

        # --- Cas où les droites ne sont PAS parallèles (elles se coupent à un point unique) ---
        # Calcul des paramètres t et u
        t = (dp[0] * v2[1] - dp[1] * v2[0]) / det
        u = (dp[0] * v1[1] - dp[1] * v1[0]) / det

        # Vérifier si le point d'intersection se trouve sur les deux segments
        # t doit être entre 0 et 1 (inclus) pour le segment 1
        # u doit être entre 0 et 1 (inclus) pour le segment 2
        if 0 - epsilon <= t <= 1 + epsilon and 0 - epsilon <= u <= 1 + epsilon:
            # Le point d'intersection est sur les deux segments
            intersection_point = p1 + t * v1
            return intersection_point
        else:
            # Le point d'intersection est en dehors d'au moins un segment
            return None

    '''
    def intersect_segments_3dmine(gridX1, gridY1, gridZ1, gridX2, gridY2, gridZ2):

        #etape 1 trouver l equation du plan1
        #l'equation du plan 1 satisfait la matrice (X,Y,Z,1) * (A,B,C,D) = (0,0,0,0)

        #etape 2 trouver un vecteur normal à ce plan (ie le produit vectoriel de deux vecteurs non colinéaire)

        x1_values = gridX1[:n]
        y1_values = gridY1[:n]
        z1_values = gridZ1[:n]
        ones_column = np.ones(n)
        zero_column = np.zeros(n)
        firstPlanMatrix = np.column_stack((x1_values, y1_values, z1_values, ones_column))
        firstPlanParameters = np.linalg.solve(firstPlanMatrix, zero_column)
        firstPlanDirectorsVector = TW_Utility.get_vector_directors(firstPlanParameters)
        firstPlanNormalVector = np.cross(firstPlanDirectorsVector[0], firstPlanDirectorsVector[1])

        x2_values = gridX2[:n]
        y2_values = gridY2[:n]
        z2_values = gridZ2[:n]
        ones_column = np.ones(n)
        zero_column = np.zeros(n)
        secondPlanMatrix = np.column_stack((x2_values, y2_values, z2_values, ones_column))
        secondPlanParameters = np.linalg.solve(secondPlanMatrix, zero_column)
        secondPlanDirectorsVector = TW_Utility.get_vector_directors(secondPlanParameters)
        secondPlanNormalVector = np.cross(secondPlanDirectorsVector[0], secondPlanDirectorsVector[1])

        #etape 3 trouver un vecteur directeur sur la ligne d'intersection de ce plan 
        normalVector = np.cross(firstPlanNormalVector, secondPlanNormalVector)
        return normalVector
    '''
            

    def intersect_segments_3d(p1, p2, p3, p4):
        """
        Calcule le point d'intersection de deux segments 3D, s'il existe.

        Args:
            p1 (np.array): Premier point du premier segment [x1, y1, z1].
            p2 (np.array): Deuxième point du premier segment [x2, y2, z2].
            p3 (np.array): Premier point du deuxième segment [x3, y3, z3].
            p4 (np.array): Deuxième point du deuxième segment [x4, y4, z4].

        Returns:
            np.array or None: Le point d'intersection si les segments se coupent,
                            sinon None.
        """
        # Convertir les points en tableaux numpy pour faciliter les opérations
        p1 = np.array(p1)
        p2 = np.array(p2)
        p3 = np.array(p3)
        p4 = np.array(p4)

        # Vecteurs directeurs des droites
        v1 = p2 - p1  # Vecteur pour le premier segment
        v2 = p4 - p3  # Vecteur pour le deuxième segment

        # 1. Vérifier le parallélisme des droites
        # Si le produit vectoriel est nul ou très proche de zéro (en raison de la précision flottante)
        cross_product_v = np.cross(v1, v2)
        if np.linalg.norm(cross_product_v) < 1e-9:  # Presque parallèle
            # Les droites sont parallèles (ou colinéaires/identiques)
            # Vérifier si elles sont colinéaires (un point de l'une est sur l'autre)
            # Si P1-P3 est colinéaire à V1, alors elles sont colinéaires
            v_p1_p3 = p1 - p3
            collinear_check = np.cross(v_p1_p3, v1)
            if np.linalg.norm(collinear_check) < 1e-9:
                # Les droites sont colinéaires (identiques ou juste sur la même ligne)
                # Nous devons vérifier le chevauchement des segments sur cette ligne
                # C'est un problème d'intersection d'intervalles 1D.
                # Projection des points sur l'axe du vecteur v1 pour simplifier
                # (ou v2, peu importe car ils sont colinéaires)
                t0 = 0.0
                t1 = np.dot(p2 - p1, v1) / np.dot(v1, v1) if np.dot(v1, v1) != 0 else 0.0
                t_p3_proj = np.dot(p3 - p1, v1) / np.dot(v1, v1) if np.dot(v1, v1) != 0 else 0.0
                t_p4_proj = np.dot(p4 - p1, v1) / np.dot(v1, v1) if np.dot(v1, v1) != 0 else 0.0

                # Assurer que les intervalles sont triés
                t_seg1_min, t_seg1_max = sorted([t0, t1])
                t_seg2_min, t_seg2_max = sorted([t_p3_proj, t_p4_proj])

                # Calcul de l'intersection des intervalles
                overlap_min = max(t_seg1_min, t_seg2_min)
                overlap_max = min(t_seg1_max, t_seg2_max)

                if overlap_min <= overlap_max:
                    # Il y a un chevauchement. S'ils sont exactement au même point, on peut le retourner.
                    # Pour cet exemple, nous retournons le point de début de chevauchement s'il y en a un seul
                    # ou None si le chevauchement est un segment (pour simplifier).
                    # Dans un cas réel, vous pourriez vouloir retourner un segment pour l'intersection.
                    if overlap_min == overlap_max: # Les segments se touchent à un point
                        return p1 + overlap_min * v1
                    else: # Les segments se chevauchent sur une longueur
                        # Pour la simplicité de cette fonction, on considère que l'intersection est un point unique
                        # Si c'est un segment, la fonction retourne None (pas de point unique d'intersection)
                        # ou vous pouvez retourner (p1 + overlap_min * v1, p1 + overlap_max * v1)
                        # Je retourne None ici pour éviter de changer le type de retour.
                        return [p1 + overlap_min * v1, p1 + overlap_max * v1]
                else:
                    return None # Pas de chevauchement pour les segments colinéaires
            else:
                return None # Droites parallèles et distinctes, pas d'intersection

        # 2. Vérifier si les droites sont gauches (ne se coupent pas et ne sont pas parallèles)
        # Si le produit mixte n'est pas nul, les droites sont gauches
        # (p3 - p1) . (v1 x v2)
        mixed_product = np.dot((p3 - p1), cross_product_v)
        if abs(mixed_product) > 1e-9: # Utiliser une tolérance pour les nombres flottants
            return None # Droites gauches, pas d'intersection

        # 3. Calculer les paramètres t et u du point d'intersection
        # Système d'équations: P1 + t*V1 = P3 + u*V2
        # Soit: t*V1 - u*V2 = P3 - P1

        # On utilise les deux premières équations (x et y) pour trouver t et u,
        # puis on vérifie avec la troisième (z).
        # Matrice des coefficients pour [x, y]
        A = np.array([
            [v1[0], -v2[0]],
            [v1[1], -v2[1]]
        ])
        # Vecteur des constantes pour [x, y]
        B = np.array([
            p3[0] - p1[0],
            p3[1] - p1[1]
        ])

        # Résoudre le système A * [t, u] = B
        # Gérer le cas où le déterminant est proche de zéro (e.g., v1[0]*(-v2[1]) - (-v2[0])*v1[1] est nul)
        # ce qui signifie que les projections 2D sont parallèles, mais elles doivent être coplanaires en 3D
        # si elles ne sont pas parallèles en 3D.
        # On peut utiliser np.linalg.solve qui gère les cas singuliers ou pseudo-inverse si nécessaire.
        try:
            t_u = np.linalg.solve(A, B)
            t = t_u[0]
            u = t_u[1]
        except np.linalg.LinAlgError:
            # Cela peut arriver si la projection 2D est singulière, mais les droites sont coplanaires en 3D.
            # Dans ce cas, nous devons choisir d'autres axes.
            # Une méthode plus robuste utilise la pseudo-inverse ou un algorithme pour des systèmes surdéterminés.
            # Ou on peut choisir un autre sous-système 2x2.
            # Par exemple, si v1[0] et v2[0] sont petits, utilisez y et z.
            A_yz = np.array([
                [v1[1], -v2[1]],
                [v1[2], -v2[2]]
            ])
            B_yz = np.array([
                p3[1] - p1[1],
                p3[2] - p1[2]
            ])
            try:
                t_u = np.linalg.solve(A_yz, B_yz)
                t = t_u[0]
                u = t_u[1]
            except np.linalg.LinAlgError:
                # Si même yz est singulier, c'est un cas très dégénéré ou un bug.
                return None


        # Vérifier que les valeurs de t et u satisfont la troisième équation (pour la robustesse)
        # p1[2] + t * v1[2] doit être proche de p3[2] + u * v2[2]
        if not np.isclose(p1[2] + t * v1[2], p3[2] + u * v2[2]):
            return None # Incohérence, les droites ne se coupent pas ou problème de précision

        # 4. Vérifier si le point d'intersection se trouve sur les deux segments
        if 0 <= t <= 1 and 0 <= u <= 1:
            # Le point d'intersection est sur les deux segments
            intersection_point = p1 + t * v1
            return intersection_point
        else:
            # Le point d'intersection des droites n'est pas sur les segments
            return None
    
    def intersect_segments_nd(p1, p2, p3, p4):
        """
        Calcule le point d'intersection de deux segments en N dimensions, s'il existe.

        Args:
            p1 (np.array): Premier point du premier segment.
            p2 (np.array): Deuxième point du premier segment.
            p3 (np.array): Premier point du deuxième segment.
            p4 (np.array): Deuxième point du deuxième segment.

        Returns:
            np.array or None: Le point d'intersection si les segments se coupent,
                            sinon None.
        """
        # Convertir les points en tableaux numpy pour faciliter les opérations
        p1 = np.array(p1, dtype=float)
        p2 = np.array(p2, dtype=float)
        p3 = np.array(p3, dtype=float)
        p4 = np.array(p4, dtype=float)

        # Vérifier que tous les points ont la même dimension
        if not (p1.shape == p2.shape == p3.shape == p4.shape):
            raise ValueError("Tous les points doivent avoir la même dimension.")
        
        n_dim = p1.shape[0]

        # Vecteurs directeurs des droites
        v1 = p2 - p1  # Vecteur pour le premier segment
        v2 = p4 - p3  # Vecteur pour le deuxième segment

        # Construire la matrice du système linéaire A * [t, u] = B
        # L'équation est: p1 + t*v1 = p3 + u*v2
        # Réécrite: t*v1 - u*v2 = p3 - p1
        # Matrice A: [v1 | -v2] (où v1 et -v2 sont des colonnes)
        # Vecteur B: p3 - p1

        A = np.column_stack((v1, -v2))
        B = p3 - p1

        # Résoudre le système linéaire A @ [t, u] = B
        # Utiliser np.linalg.lstsq pour une solution plus robuste (moindres carrés)
        # qui peut gérer les systèmes surdéterminés (plus d'équations que d'inconnues, N_dim > 2)
        try:
            # lstsq retourne (solution, résidus, rang, valeurs singulières)
            # Nous sommes intéressés par la solution [t, u]
            t_u, residuals, rank, s = np.linalg.lstsq(A, B, rcond=None)
            t = t_u[0]
            u = t_u[1]
        except np.linalg.LinAlgError:
            # Cela indique un problème grave avec la matrice A, comme des vecteurs colinéaires
            # qui n'ont pas été gérés au préalable.
            return None

        # Vérifier le rang pour déterminer si les droites s'intersectent
        # Les droites s'intersectent si le système a une solution exacte
        # C'est-à-dire si le vecteur B est dans l'espace colonne de A.
        # Ceci est équivalent à dire que le résidu est proche de zéro.
        
        # La norme des résidus doit être très petite pour qu'il y ait une intersection exacte
        # (ce qui signifie que les droites sont coplanaires et se coupent)
        if residuals.size > 0 and np.linalg.norm(residuals) > 1e-9:
            return None # Les droites ne s'intersectent pas (gauches)

        # Vérification du parallélisme si le système a une solution exacte mais le rang est bas
        # C'est un cas spécial où v1 et v2 sont parallèles (colinéaires)
        # Si le rang est 1, cela signifie que v1 et v2 sont colinéaires.
        if rank < 2: # Cela signifie que v1 et v2 sont colinéaires
            # Les droites sont parallèles (ou colinéaires/identiques)
            # Vérifier si elles sont colinéaires (un point de l'une est sur l'autre)
            # Si P1-P3 est colinéaire à V1 (ou V2)
            v_p1_p3 = p1 - p3
            # Utiliser la norme du produit vectoriel pour 3D, ou simplement le produit scalaire
            # ou le test de colinéarité généralisé pour N-D (vérifier si les vecteurs sont multiples l'un de l'autre)
            # np.linalg.norm(np.cross(v_p1_p3, v1)) si n_dim == 3
            
            # Pour N-D: vérifier si (p1-p3) est colinéaire à v1.
            # Cela peut se faire en vérifiant si rank([v1, p1-p3]) == 1
            test_collinear_matrix = np.column_stack((v1, v_p1_p3))
            if np.linalg.matrix_rank(test_collinear_matrix, tol=1e-9) <= 1:
                # Les droites sont colinéaires (identiques ou juste sur la même ligne)
                # Nous devons vérifier le chevauchement des segments sur cette ligne.
                # Projection des points sur l'axe du vecteur v1
                
                # Pour éviter la division par zéro si v1 est nul (P1 == P2)
                v1_norm_sq = np.dot(v1, v1)
                if v1_norm_sq < 1e-9: # Segment 1 est un point
                    if np.linalg.norm(p3 - p1) < 1e-9 and np.linalg.norm(p4 - p1) < 1e-9: # Si P3 et P4 sont aussi P1
                        return p1 # Les deux segments sont le même point
                    elif np.linalg.norm(p3 - p1) < 1e-9: # P1 est P3
                        if 0 <= np.dot(p4 - p3, v2) / np.dot(v2, v2) <= 1 and np.linalg.norm(v2) > 1e-9:
                            return p1 # P1 est sur le segment P3P4
                    elif np.linalg.norm(p4 - p1) < 1e-9: # P1 est P4
                        if 0 <= np.dot(p3 - p4, -v2) / np.dot(-v2, -v2) <= 1 and np.linalg.norm(v2) > 1e-9:
                            return p1 # P1 est sur le segment P3P4
                    return None # Segment 1 est un point mais pas d'intersection
                
                t0 = 0.0
                t1 = np.dot(p2 - p1, v1) / v1_norm_sq
                t_p3_proj = np.dot(p3 - p1, v1) / v1_norm_sq
                t_p4_proj = np.dot(p4 - p1, v1) / v1_norm_sq

                t_seg1_min, t_seg1_max = sorted([t0, t1])
                t_seg2_min, t_seg2_max = sorted([t_p3_proj, t_p4_proj])

                overlap_min = max(t_seg1_min, t_seg2_min)
                overlap_max = min(t_seg1_max, t_seg2_max)

                if overlap_min <= overlap_max + 1e-9: # Tolérance pour les flottants
                    if abs(overlap_min - overlap_max) < 1e-9: # Les segments se touchent à un point
                        return p1 + overlap_min * v1
                    else: # Les segments se chevauchent sur une longueur
                        return None # Retourne None si l'intersection est un segment, pas un point unique
                else:
                    return None # Pas de chevauchement pour les segments colinéaires
            else:
                return None # Droites parallèles et distinctes, pas d'intersection
                
        # Si le système a été résolu avec succès et que les résidus sont proches de zéro,
        # les droites s'intersectent.
        # Vérifier si le point d'intersection se trouve sur les deux segments
        if 0 - 1e-9 <= t <= 1 + 1e-9 and 0 - 1e-9 <= u <= 1 + 1e-9: # Ajouter une tolérance pour les bornes
            intersection_point = p1 + t * v1
            return intersection_point
        else:
            # Le point d'intersection des droites n'est pas sur les segments
            return None
        
    def _extendSegmentIfNecessary(self):
        return None