From 49350698640c7f4175aed11800da3eb5c42ec492 Mon Sep 17 00:00:00 2001 From: Theophile Terraz Date: Wed, 5 Jul 2023 15:49:30 +0200 Subject: [PATCH] add a basic st_to_m method --- src/Model/Geometry/Profile.py | 13 +++ src/Model/Geometry/ProfileXYZ.py | 25 +++++ src/Model/Geometry/Reach.py | 183 ++++++++++++++++++++++++++----- src/Scripts/plot_3DST.py | 2 +- 4 files changed, 196 insertions(+), 27 deletions(-) diff --git a/src/Model/Geometry/Profile.py b/src/Model/Geometry/Profile.py index f2678c68..407790bd 100644 --- a/src/Model/Geometry/Profile.py +++ b/src/Model/Geometry/Profile.py @@ -150,6 +150,19 @@ class Profile(object): 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(self, indexes: int): """Delete points at index diff --git a/src/Model/Geometry/ProfileXYZ.py b/src/Model/Geometry/ProfileXYZ.py index 08cf2cf4..611d9fdc 100644 --- a/src/Model/Geometry/ProfileXYZ.py +++ b/src/Model/Geometry/ProfileXYZ.py @@ -271,3 +271,28 @@ class ProfileXYZ(Profile): constant = station[index_profile_z_min[0]] 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 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 + diff --git a/src/Model/Geometry/Reach.py b/src/Model/Geometry/Reach.py index b2e35f7b..0830831f 100644 --- a/src/Model/Geometry/Reach.py +++ b/src/Model/Geometry/Reach.py @@ -459,7 +459,7 @@ class Reach: # Interpolate @timer - def interpolate_step(self, lim: list[int], guideline1: str, guideline2: str, step: float, longmi: int, linear: bool) + 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 @@ -476,9 +476,9 @@ class Reach: """ 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] >= len(self._profiles): + 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]: + if lim[1] < 0 or lim[0] > lim[1]: print('TODO: CRITICAL ERROR!') comp, inconmp = self.compute_guidelines() if not guideline1 in comp: @@ -487,29 +487,29 @@ class Reach: 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]) + 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 + 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) then + if (longmi == 1): d = max(d1, d2) - else if (longmi == 2) then + elif (longmi == 2): d = min(d1, d2) - else + else: d = 0.5*(d1 + d2) nbsti[i] = nint(d / pas) - 1 - if(nbsti[i].lt.0) nbsti[i] = 0 + if(nbsti[i] < 0): nbsti[i] = 0 - for i in range(lim[1]-1,len(self._profiles)-1) + for i in range(lim[1]-1,self.number_profiles-1): nbsti.append(0) return self.interpolate_nprof(nbsti, lineaire) @@ -547,10 +547,10 @@ class Reach: sections_out = [] kps = self.get_kp() for isect_in in range(len(nbsti)): - sections_out.append(self.section(isect_in)) + sections_out.append(self.profiles(isect_in)) dkp = kps[isect_in + 1] - kps[isect_in] - sect_in1 = self.section(isect_in) - sect_in2 = self.section(isect_in+1) + 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) @@ -565,11 +565,142 @@ class Reach: reach = self)) for i in range(len(x1)): # can do better... # add the points to the created section - sections_out[-1].add(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] - ) + 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.section(self.number_profiles()-1)) + 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] + 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]) + 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]) + diff --git a/src/Scripts/plot_3DST.py b/src/Scripts/plot_3DST.py index 2647a32a..d673c76f 100644 --- a/src/Scripts/plot_3DST.py +++ b/src/Scripts/plot_3DST.py @@ -32,7 +32,7 @@ def set_axes_equal(ax): ax.set_xlim3d([x_middle - plot_radius, x_middle + plot_radius]) ax.set_ylim3d([y_middle - plot_radius, y_middle + plot_radius]) - ax.set_zlim3d([z_middle - plot_radius, z_middle + plot_radius]) + ax.set_zlim3d([z_middle - 0.1*plot_radius, z_middle + 0.1*plot_radius]) st_file = sys.argv[1] my_reach = Reach(None)