remove GL doublons

src/Model/Geometry/PointXYZ.py

@@ -16,7 +16,7 @@
 
 # -*- coding: utf-8 -*-
 
-from math import dist
+from math import dist, sqrt, acos
 import numpy as np
 
 from Model.Tools.PamhyrDB import SQLSubModel
@@ -204,3 +204,55 @@ class PointXYZ(Point, SQLSubModel):
             Euclidean distance between the two points
         """
         return dist((p1.x, p1.y, p1.z), (p2.x, p2.y, p2.z))
+
+    @staticmethod
+    def angle(p1, p2, p3):
+        """angle p1 p2 p3.
+
+        Args:
+            p1: A XYZ Point
+            p2: A XYZ Point
+            p3: A XYZ Point
+
+        Returns:
+            angle p1 p2 p3
+        """
+        a = PointXYZ.distance(p1, p2)
+        b = PointXYZ.distance(p2, p3)
+        c = PointXYZ.distance(p1, p3)
+        return acos((a**2 + b**2 - c**2) / (2 * a * b))
+
+
+    @staticmethod
+    def triangle_height(p1, p2, p3):
+        """triangle p1 p2 p3 height starting at p2.
+
+        Args:
+            p1: A XYZ Point
+            p2: A XYZ Point
+            p3: A XYZ Point
+
+        Returns:
+            triangle height starting at p2
+        """
+        a = PointXYZ.distance(p1, p2)
+        b = PointXYZ.distance(p2, p3)
+        c = PointXYZ.distance(p1, p3)
+        h2 = a**2 - ((c**2 - b**2 + a**2) / 2*c)**2
+        return sqrt(h2)
+
+    @staticmethod
+    def triangle_area(p1, p2, p3):
+        """triangle p1 p2 p3 area.
+
+        Args:
+            p1: A XYZ Point
+            p2: A XYZ Point
+            p3: A XYZ Point
+
+        Returns:
+            triangle height starting at p2
+        """
+        h = PointXYZ.triangle_height(p1, p2, p3)
+        d = PointXYZ.distance(p1, p3)
+        return h*d/2

src/Model/Geometry/Profile.py

@@ -214,6 +214,20 @@ class Profile(object):
         self.points.insert(index, point)
         self._status.modified()
 
+
+    def add_point(self, point:Point):
+        """append point at end of profile.
+
+        Args:
+            point: The point.
+
+        Returns:
+            Nothing.
+        """
+        self._points.append(point)
+        self._status.modified()
+
+
     def delete_i(self, indexes: list):
         self._points = list(
             map(

src/Model/Geometry/ProfileXYZ.py

@@ -440,6 +440,77 @@ class ProfileXYZ(Profile, SQLSubModel):
 
             return list(map(lambda s: s - constant, station))
 
+    @timer
+    def length(self) -> list:
+        l = []
+        x = self.x()
+        y = self.y()
+        z = self.z()
+        for i in range(len(self.points)-1):
+            dx = x[i+1] - x[i]
+            dy = y[i+1] - y[i]
+            dz = z[i+1] - z[i]
+            l.append( np.sqrt(dx*dx + dy*dy + dz*dz) )
+
+        return l
+
+
+    @timer
+    def cumul_length(self) -> list:
+        lcum = []
+        lc  = 0.0
+        l = self.length()
+        for i in l:
+            lc = lc + i
+            lcum.append(lc)
+
+        return lcum
+
+    @timer
+    def purge1(self, np) -> list:
+
+        points = self.points
+        if (self.number_points <= np): return points
+        while len(points) > np:
+            x1 = [point.x() for point in points]
+            y1 = [point.y() for point in points]
+            z1 = [point.z() for point in points]
+            l = [(x*x+y*y+z*z) for (x,y,z) in zip(x1,y1,z1)]
+            l2 = [(x + y) for (x, y) in zip(l[:-1], l[1:])]
+            imin = np.argmin(l2)+1
+            points.pop(imin)
+
+        return points
+
+    @timer
+    def purge2(self, np) -> list:
+
+        points = self.points
+        if (self.number_points <= np): return points
+        while len(points) > np:
+            l2 = []
+            for i in range(len(points)-2):
+                l2.append(PointXYZ.angle(points(i),points(i+1),points(i+2)))
+            imax = np.argmax(l2)+1
+            points.pop(imax)
+
+        return points
+
+    @timer
+    def purge3(self, np) -> list:
+
+        points = self.points
+        if (self.number_points <= np): return points
+        while len(points) > np:
+            l2 = []
+            for i in range(len(points)-2):
+                l2.append(PointXYZ.triangle_area(points(i),points(i+1),points(i+2)))
+            imin = np.argmin(l2)+1
+            points.pop(imin)
+
+        return points
+
+
     def width_approximation(self):
         rg = self.get_point_by_name("rg")
         rd = self.get_point_by_name("rd")

src/Model/Geometry/Reach.py

@@ -579,6 +579,266 @@ class Reach(SQLSubModel):
                 file_st.write("    999.9990     999.9990     999.9990")
                 file_st.write("\n")
 
+    # Interpolate
+
+    @timer
+    def interpolate_step(self, lim: list, guideline1: str, guideline2: str, step: float, longmi: int, linear: bool):
+        """interpolate with a given step between profile lim[0] and profile lim[1]
+        Args:
+            lim: index of profiles for interpolation
+            guideline1: first guideline for distance
+            guideline2: second guideline for distance
+            step: interpolation step
+            longmi: 1 if distance computed on min(dist(guideline1),dist(guideline2))
+                    2 if distance computed on max(dist(guideline1),dist(guideline2))
+                    else: distance computed on (dist(guideline1)+dist(guideline2))/2
+            linear: true if linear interpolation, false if spline interpolation
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+        if len(lim) != 2:
+            print('TODO: CRITICAL ERROR!') # error for now, only one interval
+        if lim[0] < 0 or lim[0] >= lim[1] or lim[1] >= self.number_profiles:
+            print('TODO: CRITICAL ERROR!')
+        if lim[1] < 0 or lim[0] > lim[1]:
+            print('TODO: CRITICAL ERROR!')
+        comp, inconmp = self.compute_guidelines()
+        if not guideline1 in comp:
+            print('TODO: CRITICAL ERROR!')
+        if not guideline2 in comp:
+            print('TODO: CRITICAL ERROR!')
+
+        nbsti = [] # number of profiles to add between self._profile[i] and self._profile[i+1]
+        for i in range(lim[0]):
+            nbsti.append(0)
+        for i in range(lim[0],lim[1]-1):
+            points = self.profile(i).points
+            point11 = next(filter(lambda p: p.name == guideline1, points))
+            point12 = next(filter(lambda p: p.name == guideline2, points))
+            points = self.profile(i+1).points
+            point21 = next(filter(lambda p: p.name == guideline1, points))
+            point22 = next(filter(lambda p: p.name == guideline2, points))
+            d1 = point11.dist(point21)
+            d2 = point12.dist(point22) # warning: 3D distance
+
+            if (longmi == 1):
+                d = max(d1, d2)
+            elif (longmi == 2):
+                d = min(d1, d2)
+            else:
+                d = 0.5*(d1 + d2)
+
+            nbsti[i] = nint(d / pas) - 1
+            if(nbsti[i] < 0): nbsti[i] = 0
+
+        for i in range(lim[1]-1,self.number_profiles-1):
+            nbsti.append(0)
+
+        return self.interpolate_nprof(nbsti, lineaire)
+
+    @timer
+    def interpolate_nprof(self, nbsti: List[int], linear: bool):
+        """add nbsti[i] profiles between self._profile[i] and self._profile[i+1]
+
+        Args:
+            nbsti: number of profiles to add between self._profile[i] and self._profile[i+1]
+            linear: true if linear interpolation, false if spline interpolation
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+
+        if(lineaire):
+            return self.interpolate_linear(nbsti)
+        else:
+            # we might add the spline interpolation later
+            print('TODO: CRITICAL ERROR!')
+
+    @timer
+    def interpolate_linear(self, nbsti: List[int]):
+        """add nbsti[i] profiles between self._profile[i] and self._profile[i+1]
+
+        Args:
+            nbsti: number of profiles to add between self._profile[i] and self._profile[i+1]
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+        if len(nbsti) != self.number_profiles()-1:
+            print('TODO: CRITICAL ERROR!') # error for now, only one interval
+        sections_out = []
+        kps = self.get_kp()
+        for isect_in in range(len(nbsti)):
+            sections_out.append(self.profiles(isect_in))
+            dkp = kps[isect_in + 1] - kps[isect_in]
+            sect_in1 = self.profiles(isect_in)
+            sect_in2 = self.profiles(isect_in+1)
+            for isect_loc in range(nbsti[isect_in]):
+                # RATIO entre les deux sections initiales
+                d=float(isect_loc)/float(nbsti(isect_in)+1)
+                x1 = sect_in1.x()
+                y1 = sect_in1.y()
+                z1 = sect_in1.z()
+                x2 = sect_in2.x()
+                y2 = sect_in2.y()
+                z2 = sect_in2.z()
+                tags = sect_in1.names()
+                section_out.append(ProfileXYZ(kp = kps[isect_in] + d*dkp,
+                                              reach = self))
+                for i in range(len(x1)): # can do better...
+                    # add the points to the created section
+                    sections_out[-1].add_point(PointXYZ(x = x1[i] + d*(x2[i] - x1[i]),
+                                                        y = y1[i] + d*(y2[i] - y1[i]),
+                                                        z = z1[i] + d*(z2[i] - z1[i]),
+                                                        name = tags[i]))
+        # don't forget the last one:
+        sections_out.append(self.profiles[self.number_profiles-1])
+        return sections_out
+
+    @timer
+    def st_to_m(self):
+        """convert a ST bief to a M bief (same number of points on every section)
+
+        Args:
+            nbsti: number of profiles to add between self._profile[i] and self._profile[i+1]
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+        isectmax = 0
+        npoint_max = 0
+        new_profiles = []
+        np = []
+        nprof = len(self.number_profiles)
+        # find the profile with more points
+        for isect in range(nprof):
+            new_profiles.append(self.profile(isect))
+            points = new_profiles[isect].points
+            np.append(len(points))
+            if(npoint_max < len(points)):
+                npoint_max = len(points)
+                isectmax = isect
+
+        # let's go ! up from isectmax
+        for isect in range(isectmax+1,nprof):
+            # TODO case new_profiles[isect].length() = 0.0
+            # TODO case new_profiles[isect-1].length() = 0.0
+            l = new_profiles[isect-1].length() # length of each segment
+            lcum = new_profiles[isect-1].cumul_length() # lcum[-1] is length
+            alpha = [i/lcum[-1] for i in lcum]
+            l = new_profiles[isect].length()
+            lcum = new_profiles[isect].cumul_length()
+            beta = [i/lcum[-1] for i in lcum]
+            for i in range(np[isect], npoints_max):
+                beta.append(1.0)
+
+
+
+            for k in range(np[isect]): # boucle sur les points à rajouter
+                trouve = False
+                for j0 in range(np[isect]-1):
+                    if(trouve):break
+                    undeplus=True
+                    for j1 in range(j0, np[isect]-1):
+                        if(undeplus):
+                            if(beta[j1].le.alpha[j1]):
+                                undeplus=False
+                                break
+
+                    if (undeplus):
+                        if (j0 == np[isect]-2):
+                            trouve=True
+                        else:
+                            j1=j0+1
+                            if (beta[j0]-alpha[j0] > alpha[j1]-beta[j0]): # beta(j0) > (alpha(j1)+alpha(j2))/2
+                                trouve = True
+
+
+                    if (trouve):
+                    # on rajoute un point
+                        ratio=(alpha[j0]-beta[j0-1])/(beta[j0]-beta[j0-1])
+                        if (ratio < 0.0):
+                            # on double le point a gauche
+                            new_profiles[isect].insert_point(j0,new_profiles[isect].points[j0])
+                            beta[j0]=beta[j0-1]
+                            new_profiles[isect].points[j0+1].name('')
+                        else:
+                            p = PointXYZ(x = (1.0-ratio)*new_profiles[isect].points[j0].x + ratio*new_profiles[isect].points[j0].x,
+                                         y = (1.0-ratio)*new_profiles[isect].points[j0].y + ratio*new_profiles[isect].points[j0].y,
+                                         z = (1.0-ratio)*new_profiles[isect].points[j0].z + ratio*new_profiles[isect].points[j0].z
+                                        )
+                            new_profiles[isect].insert_point(j0,p)
+
+            n = len(new_profiles[isect].points)
+            if n <= npoint_max: # we add the missing points at the end of the profile
+                for k in range(n,npoints_max): new_profiles[isect].add_point(new_profiles[isect].points[-1])
+
+
+        # let's go ! down from isectmax
+        for isect in range(isectmax-1,-1,-1):
+            # TODO case new_profiles[isect].length() = 0.0
+            # TODO case new_profiles[isect-1].length() = 0.0
+            l = new_profiles[isect+1].length() # length of each segment
+            lcum = new_profiles[isect+1].cumul_length() # lcum[-1] is length
+            alpha = [i/lcum[-1] for i in lcum]
+            l = new_profiles[isect].length()
+            lcum = new_profiles[isect].cumul_length()
+            beta = [i/lcum[-1] for i in lcum]
+            for i in range(np[isect], npoints_max):
+                beta.append(1.0)
+
+
+
+            for k in range(np[isect]): # boucle sur les points à rajouter
+                trouve = False
+                for j0 in range(np[isect]-1):
+                    if(trouve):break
+                    undeplus=True
+                    for j1 in range(j0, np[isect]-1):
+                        if(undeplus):
+                            if(beta[j1].le.alpha[j1]):
+                                undeplus=False
+                                break
+
+                    if (undeplus):
+                        if (j0 == np[isect]-2):
+                            trouve=True
+                        else:
+                            j1=j0+1
+                            if (beta[j0]-alpha[j0] > alpha[j1]-beta[j0]): # beta(j0) > (alpha(j1)+alpha(j2))/2
+                                trouve = True
+
+
+                    if (trouve):
+                    # on rajoute un point
+                        ratio=(alpha[j0]-beta[j0-1])/(beta[j0]-beta[j0-1])
+                        if (ratio < 0.0):
+                            # on double le point a gauche
+                            new_profiles[isect].insert_point(j0,new_profiles[isect].points[j0])
+                            new_profiles[isect].points[j0+1].name('')
+                            beta[j0]=beta[j0-1]
+                        else:
+                            p = PointXYZ(x = (1.0-ratio)*new_profiles[isect].points[j0].x + ratio*new_profiles[isect].points[j0].x,
+                                         y = (1.0-ratio)*new_profiles[isect].points[j0].y + ratio*new_profiles[isect].points[j0].y,
+                                         z = (1.0-ratio)*new_profiles[isect].points[j0].z + ratio*new_profiles[isect].points[j0].z,
+                                        )
+                            new_profiles[isect].insert_point(j0,p)
+
+            n = len(new_profiles[isect].points)
+            if n <= npoint_max: # we add the missing points at the end of the profile
+                for k in range(n,npoints_max): new_profiles[isect].add_point(new_profiles[isect].points[-1])
+
+    @timer
+    def purge(self, np: int) -> list:
+        profiles = []
+        for profile in self.profiles():
+            points = profile.purge2(np)
+            profiles.append(ProfileXYZ())
+            profiles[-1].import_points(points)
+
+        return profiles
+
     def get_incline(self):
         first = self.profile(0)
         last = self.profile(len(self) - 1)