Source code for raidnr.util.ext_methods


# ---------------------------------------------------
# ¦¦¦¦¦¦¦    ¦¦¦¦¦   ¦¦¦¦  ¦¦¦¦¦¦¦¦  
#  ¦¦   ¦¦  ¦¦   ¦¦   ¦¦    ¦¦   ¦¦¦   
#  ¦¦¦¦¦¦   ¦¦¦¦¦¦¦   ¦¦    ¦¦   ¦¦¦   
#  ¦¦   ¦¦  ¦¦   ¦¦   ¦¦    ¦¦   ¦¦¦  
# ¦¦¦   ¦¦¦ ¦¦   ¦¦¦ ¦¦¦¦  ¦¦¦¦¦¦¦¦    
# ---------------------------------------------------
# Bryden Wood
# License: MIT, see full license in LICENSE.txt
# Web: repository-link
#-----------------------------------------------------
# Date: 30.09.2019
# Authors: Claudio Campanile, Konstantina Spanou
#-----------------------------------------------------
# Description
#-----------------------------------------------------
# Collection of extension functions for mixed use.
# Mainly geometrical and mathematical functions
#-----------------------------------------------------
# Docs Reference : link to module docs
#-----------------------------------------------------

# modules
import json
import math
import glob
from math import sqrt
# scientific modules
import numpy as np
import pandas as pd, numpy as np
from numpy.linalg import norm
from numpy import inf
from scipy import interpolate
import matplotlib.pyplot as plt
# geospatial modules
import geopandas as gpd, fiona, os, shapely
from shapely.geometry import LineString,Point,Polygon,shape, mapping, MultiLineString,MultiPoint
from shapely.ops import nearest_points,unary_union,cascaded_union
from shapely.ops import split
from pyproj import Proj, transform
import networkx as nx

# raidnr
from raidnr.config import CRS, DEBUG
from raidnr.bwgrid import grid_transformations as gt


speed=0


[docs]def resample_pts(pt1, pt2, n): """ Return a list of points in between two ND-points at equal distance. Parameters ---------- pt1 : tuple initial point pt2 : tuple final point Returns ------- pts : list<tuple> """ if n < 1: return [] # from no. of pts to no. of segments n += 1 pts = [] for i in range(1,n): # note: coord1 means the coord of pt1. its iterations through all the coord make this methods suitable for ND points (2D, 3D, etc) pt = tuple([coord1*float(1.-(i/float(n))) + coord2*float(i/float(n)) for coord1,coord2 in zip(pt1,pt2)]) pts.append(pt) return pts
[docs]def pt_dist_sq(pt1, pt2): """ Calculate the squared distance between two ND-points. Parameters ---------- pt1 : tuple initial point pt2 : tuple final point Returns ------- float squared distance """ return sum([(c1-c2)**2 for c1,c2 in zip(list(pt1), list(pt2))])
[docs]def resample_poly(linestring, limit_distance, rmv_duplicates=True, tolerance = 0.01): """ Resample a shapely.geometry.Linestring to avoid segments above a length threshold. Parameters ---------- linestring : shapely.geometry.Linestring the linestring to resample. limit_distance : float maximum distance a segment can be. rmv_duplicates : bool remove duplicate points. tolerance: float distance below which two points are considered duplicates. Returns ------- shapely.geometry.Linestring resampled Linestring object """ pts = [] for i,pt in enumerate(linestring.coords[:-1]): if i == 0: pts.append(pt) pt2= linestring.coords[i+1] d = pt_dist_sq(pt,pt2) # working with squared of distances if d < tolerance and rmv_duplicates: #remove duplicate points continue if d > limit_distance**2.: #insert mid-points to address limit_distance mids = resample_pts(pt,pt2, math.ceil(d / limit_distance**2)) else: mids = [] pts += mids + [pt2] return LineString(pts)
[docs]def reconstruct_poly(linestring, step_distance): """ Resample curve by interpolation. Parameters ---------- linestring: Shapely.geometry.LineString Polyline to interpolate step_distance: float granularity on resampled, interpolated curve Returns ------- np.array link: https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html#spline-interpolation """ tck, u = interpolate_poly(linestring) delta_u = step_distance / linestring.length unew = np.arange(0, 1.0 + delta_u, delta_u) return points_from_spline(tck, unew)
[docs]def interpolate_poly(linestring): """ Resample curve by interpolation. Parameters ---------- linestring: Shapely.geometry.LineString Polyline to interpolate step_distance: float granularity on resampled, interpolated curve Returns ------- np.array link: https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html#spline-interpolation """ pts = np.array(linestring.coords) tck, u = interpolate.splprep(pts.transpose(), s=0) return tck, u
[docs]def bezier(poly, s): """ Deprecated ref: https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html#spline-interpolation """ pts = poly tck, u = interpolate.splprep(pts.transpose(), s=10000.) return tck, u
[docs]def points_from_spline(spline, u): """ Extract points at u parameters. Parameters ---------- spline : spline3d numpy spline function u : list<float> parameters at which extract points Returns ------- nd-array points along a spline """ return np.array(interpolate.splev(u, spline)).transpose()
[docs]def tangent_vector(pt, tan, unitize=True, mult=1.0): """ create a tangent vector from a point and a direction. Parameters ---------- pt : numpy.array vector application point tan : numpy.array vector Returns ------- numpy.array (start, end) positions, where start and end are (x,y[,z]) arrays """ if unitize: norm = np.linalg.norm(tan) tan = tan/norm delta = pt + tan*mult return np.array([pt, delta], dtype=float).transpose()
[docs]def circle_3pt(A,B,C): """ Calculate a circple by three points. Parameters ---------- A, B, C : array-like points to which the circle passes through, as [x,y,z] link: https://stackoverflow.com/questions/20314306/find-arc-circle-equation-given-three-points-in-space-3d """ if not (isinstance(A,np.ndarray) and isinstance(B,np.ndarray) and isinstance(C,np.ndarray)): raise ValueError("A,B,C must be numpy arrays") if not ((len(A)==3 and len(B)==3 and len(C)==3) or (len(A)==2 and len(B)==2 and len(C)==2)): raise ValueError("A,B,C must be either 2D or 3D points") # Triangle Lengths a = np.linalg.norm(B - A) b = np.linalg.norm(C - B) c = np.linalg.norm(C - A) try: # Heron's formula for triangle's surface k = np.sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c))) / 4. den = a*b*c # Denumerator; make sure there is no division by zero. if den == 0.0: # Very unlikely, but just to be sure curvature = 0.01 else: curvature = 4.*k / den return 1./curvature except Exception as e: print("ERROR: {}".format(e)) return 0
[docs]def moving_average(a, n=3) : """ Moving average of an array Parameters ---------- a : array-like array to calculate the moving average of Returns ------- array modified array """ if not isinstance(a, np.ndarray): a = np.array(a) if n > 0: # head s = [] for ii,i in enumerate(a[:n-1]): r = np.sum(a[max(0,n-ii):ii+n]) / ((ii+n) - max(0,n-ii)) s.append(r) s = np.array(s) #rest ret = np.cumsum(a, dtype=float) ret[n:] = ret[n:] - ret[:-n] return np.append(s, ret[n - 1:] / n) else: return a
[docs]def speedOnDiffCant(R, cant_max=150, cant_small=100,cant_xsmall=50,cant_min=25,max_speed=200): """ """ global speed if (R >= 200) & (R <= 5000): speed = sqrt(R*cant_max/11.82)*(5/18) elif R < 100: speed = sqrt(R*cant_min/11.82)*(5/18) elif (R < 200) & (R >= 150): speed = sqrt(R*cant_small/11.82)*(5/18) elif (R < 150) & (R >= 100): speed = sqrt(R*cant_xsmall/11.82)*(5/18) elif R > 5000: speed = max_speed*(5/18) return speed - (speed%5)
[docs]def repeatForEach(elements, times): """ link: https://github.com/matplotlib/matplotlib/issues/8484 """ return [element for element in elements for _ in range(times)]
[docs]def renderColorsForQuiver3d(colors): """ Used to fix the colors of a quiver 3d. link: https://github.com/matplotlib/matplotlib/issues/8484 """ # filter out 0-length vector corresponding to rgb=(0,0,0) that won't be drawn colors = list(filter(lambda x: x!=(0.,0.,0.), colors)) # set up rgb color for each lines of arrows return colors + repeatForEach(colors, 2)
[docs]def smoothen_poly(poly, tolerance, as_ratio=False): """ Simplify a polyline to not exceed a maximum distance from the input. Parameters ---------- poly: numpy.array a set of 2D or 3D points tolerance: float maximum allowed distance Returns ------- numpy.array array (x,y[,z]) points """ s = poly.shape assert (len(s) == 2), "poly is not a 2-d numpy.array" assert (poly.shape[1] == 2 or poly.shape[1] == 3), "data points must be either 2d or 3d, not {}-d".format(poly.shape[1]) # calc furthest point from start-end line i,d = max_var(poly) d = d if not math.isnan(d) else 0.0 p1, p2 = poly[0], poly[-1] if as_ratio: ratio = d/np.linalg.norm(p2-p1) if (ratio <= tolerance) or (len(poly)<3) or (i==0) or (i==(len(poly)-1)): return np.array([p1, p2]) else: return np.concatenate((smoothen_poly(poly[:i], tolerance, True)[:-1], smoothen_poly(poly[i:], tolerance, True))) if (d<=tolerance) or (len(poly)<3) or (i==0) or (i==(len(poly)-1)): return np.array([p1, p2]) else: # split to debug r1 = smoothen_poly(poly[:i], tolerance)[:-1] r2 = smoothen_poly(poly[i:], tolerance) r = np.concatenate((r1, r2)) return r
[docs]def max_var(poly): """ Given a polyline, finds the furthest point from the line constructed as (start-end) Parameters ---------- poly: numpy.array array (x,y[,z]) points Returns ------- tuple (index, distance) of the firthest point """ p1 = poly[0] p2 = poly[-1] d = [np.linalg.norm(np.cross(p2-p1,p3-p1)/np.linalg.norm(p2-p1)) for p3 in poly] i = np.arange(0,len(poly),1) r = sorted(zip(i,d), key=lambda x:x[1], reverse=True)[0] return r
[docs]def angle_deviation_3pt(p0,p1=np.array([0,0]),p2=None): """ Compute angle (in degrees) for p0p1p2 corner Parameters ---------- p0,p1,p2,p3 : array-like points in the form of [x,y] Returns ------- float angle """ if p2 is None: p2 = p1 + np.array([1, 0]) v0 = np.array(p0) - np.array(p1) v1 = np.array(p2) - np.array(p1) angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1)) return np.degrees(angle)
[docs]def extract_elevation(pt, dtm_grid, layer='DTM'): """ Extract elevation of an alignment in relation to a DTM Parameters ---------- pt : numpy.array 3-d array in world coordinates dtm_grid : raidnr.bwgrid.grid.Grid Grid object layer : str [optional] the layer to extract from the grid. Note: this function can be used to extract other types of layers (e.g. Woodland) Returns ------- float average of the four closest values, weighted by distance ratios """ # assert isinstance(dtm_grid, Grid), "dtm_grid must be a Grid object, given {}".format(type(dtm_grid)) # calc raster coords pt=dtm_grid.to_raster(pt, to_int=False) # find 4 closest coordinates j0=math.floor(pt[1]) j1=math.ceil(pt[1]) i0=math.floor(pt[0]) i1=math.ceil(pt[0]) # calc 4 distances squared d1 = (pt[0]-i0)**2 + (pt[1]-j0)**2 d2 = (pt[0]-i0)**2 + (pt[1]-j1)**2 d3 = (pt[0]-i1)**2 + (pt[1]-j0)**2 d4 = (pt[0]-i1)**2 + (pt[1]-j1)**2 dt = d1+d2+d3+d4 # take altitudes e1 = dtm_grid[layer][i0,j0] e2 = dtm_grid[layer][i0,j1] e3 = dtm_grid[layer][i1,j0] e4 = dtm_grid[layer][i1,j1] #return weighted average elevation = (e1*d1+e2*d2+e3*d3+e4*d4)/dt if(elevation<0): return 0 return elevation
[docs]def distance_deviation_3pt(p1,p2,p): """ WIP """ d = norm(np.cross(p2-p1, p1-p))/norm(p2-p1) return d
# IO
[docs]def list_files(path): """ List all files in a directory Parameters ---------- path : str folder path """ fs=[] for f in os.listdir(path): if os.path.isfile(os.path.join(path,f)): fs.append(f) return fs
[docs]def group_files(fns): """ group files by filename. Useful for shapefiles Parameters ---------- fn : list<str> list of filenames """ groups = {} for fn in fns: name = str(os.path.splitext(os.path.split(fn)[1])[0]) if not name in groups.keys(): groups[name] = [] groups[name].append(fn) return groups
[docs]def filter_files(path, ext='shp', no_of_groups=None): """ Filter groups of files by extension. Useful to filter .shp Parameters ---------- path : str path to search files into ext : str [optional] file extension. default='.shp' no_of_groups: int or None [optional] check that only n files exists with that extension (i.e. n groups) Returns ------- list<str> file paths selected """ fns = list_files(path) groups = group_files(fns) if no_of_groups: assert no_of_groups == len(groups), "Found wrong number of file groups. Check if your process is saving the files in the wrong location" # extract the file with required extension files = glob.glob(path + "/*." + ext) return files
[docs]def fetch_shp(shpPath,name): """ WIP """ directory=os.path.join(shpPath,name) files=os.listdir(directory) for i in files: if i.split(".")[1]=="shp": return os.path.join(directory,i)
[docs]def load_speed_limits(fileName, output_full_df=False): """ WIP """ railway_lines=gpd.read_file(fileName) # distances = list(railway_lines['geometry'].apply(lambda x: int(x.length) if int(x.length)>0 else 1)) speeds=np.array(railway_lines['Speed_mph']) speeds[0]=30 speeds=speeds.astype('float64') speeds=speeds/2.237 railway_lines['speed_ms']=speeds # railway_lines['geometry'] = railway_lines['geometry'].to_crs(epsg=27000) segments=railway_lines['geometry'] if output_full_df: return railway_lines, speeds else: return segments, speeds
[docs]def purge_poly(fn): """ Purge a collection of segments and build up one polyline as LineString object Parameters ---------- fn : str shapefile filename Returns ------- LineString re-assembled polyine """ rail=gpd.read_file(fn) #Visulise the points of that form the lines if DEBUG: print(rail[:5]) points=[] for i in range(len(rail)): line=rail.geometry[i] for j in range(len(line.coords[:])): points.append(Point(line.coords[:][j])) fig, ax = plt.subplots(figsize=(12,9)) gpd.GeoDataFrame(points,columns={'geometry':[]}).plot(ax=ax,markersize=5) x,y=points[0].x,points[0].y # the first point in the geodataframe plt.scatter(x,y) plt.show() count=0 fromLine=[] starts, ends, ids = [], [], [] lines = rail['geometry'] for line in rail['geometry']: starts.append(Point(line.coords[:][0])) ends.append(Point(line.coords[:][-1])) ids.append(count) count+=1 #this dataframe will be used to order the line segments fromLine=pd.DataFrame(data={'geometry':rail['geometry'],'pt1':starts,'pt2':ends,'id':ids}) # #this dataframe will be used to order the line segments # fromLine=pd.DataFrame(fromLine,columns={'geometry':[],'pt1':[],'pt2':[],'id':[]}) G=nx.read_shp(fn) #the maximum number of edges per node - identify if there 'junctions' max_node_degree=max(list(zip(*list(G.degree)))[1]) # if max_node_degree==2: # print('There are no branches') # print(len(fromLine)) #choose a random line i=np.random.randint(0,len(starts)) origin=fromLine['geometry'][i] #create an empty list to store the ordered coords lineCoords1=[] #the ids idList=list(fromLine['id']) #take the first subset with point pt1 firstPoint=fromLine['pt1'][i] #the second will search for the point pt2 secondPoint=fromLine['pt2'][i] #create the list to store ordered coordinates - add the first point lineCoords1=firstPoint.coords[:] #remove the origin line from the list idList.remove(fromLine['id'][i]) #create a new list sample=fromLine.loc[fromLine['id'].apply(lambda x: x in idList)].reset_index(drop=True) samplePoint=firstPoint a=1 while a==1: if samplePoint in list(sample['pt1']): nextLine=sample.loc[sample['pt1'].apply(lambda x: x==samplePoint )].reset_index(drop=True) lineCoords1=lineCoords1 + nextLine['geometry'][0].coords[:] samplePoint=nextLine['pt2'][0] idList.remove(nextLine['id'][0]) sample=sample.loc[sample['id'].apply(lambda x: x in idList)].reset_index(drop=True) a=1 elif samplePoint in list(sample['pt2']): nextLine=sample.loc[sample['pt2'].apply(lambda x: x==samplePoint )].reset_index(drop=True) lineCoords1=lineCoords1 + list(reversed(nextLine['geometry'][0].coords[:])) samplePoint=nextLine['pt1'][0] idList.remove(nextLine['id'][0]) sample=sample.loc[sample['id'].apply(lambda x: x in idList)].reset_index(drop=True) a=1 else: samplePoint=Point(100,100) # print('No more') a=0 samplePoint=secondPoint lineCoords2=secondPoint.coords[:] #add the second point of the reference line a=1 while a==1: if samplePoint in list(sample['pt1']): nextLine=sample.loc[sample['pt1'].apply(lambda x: x==samplePoint )].reset_index(drop=True) lineCoords2=lineCoords2 + nextLine['geometry'][0].coords[:] samplePoint=nextLine['pt2'][0] idList.remove(nextLine['id'][0]) sample=sample.loc[sample['id'].apply(lambda x: x in idList)].reset_index(drop=True) a=1 elif samplePoint in list(sample['pt2']): nextLine=sample.loc[sample['pt2'].apply(lambda x: x==samplePoint )].reset_index(drop=True) lineCoords2=lineCoords2 + list(reversed(nextLine['geometry'][0].coords[:])) samplePoint=nextLine['pt1'][0] idList.remove(nextLine['id'][0]) sample=sample.loc[sample['id'].apply(lambda x: x in idList)].reset_index(drop=True) a=1 else: samplePoint=Point(100,100) # print('No more') a=0 lineCoords=[] if firstPoint.coords[:][0] in lineCoords1: lineCoords=list(reversed(lineCoords1)) + origin.coords[:] + lineCoords2 else: lineCoords=list(reversed(lineCoords2)) + origin.coords[:] + lineCoords1 oneLine=LineString(lineCoords) # test line created and initial line length (assuming no double lines exist) evaluation=rail['geometry'].length.sum() #the length of the original shapefile # if round(evaluation,3) == round(oneLine.length,3): # print('Line lengths match!') # else: # print('WARNING: line length mismatch') oneLine=gpd.GeoDataFrame(geometry=gpd.GeoSeries(oneLine)) oneLine.crs=CRS return oneLine
[docs]def snap_bbox(bbox, path=None): """ Snap a bounding box to the w3w grid Parameters ---------- bbox : array-like [min_x, min_y, max_x, max_y] fn : str folder where to save the bbox shapefile Returns ------- array-link [min_x, min_y, max_x, max_y] snapped to the w3w grid """ # user defined test_bbox=Polygon([[bbox[0],bbox[1]],[bbox[2],bbox[1]],[bbox[2],bbox[3]],[bbox[0],bbox[3]]]) # build a gpd df bbox_gpd = gpd.GeoDataFrame(geometry=[test_bbox]) # reproject bbox_gpd.crs = {'init' :'epsg:27700'} # align with the grid NW,SE = gt.get_UG_grid_bounds(bbox_gpd,corners=True) # adjust the new bounds minx,maxy=NW[0][0],NW[1][0] maxx,miny=SE[0][2],SE[1][2] # same procedure bbox = [minx, miny, maxx, maxy] if path: test_bbox=Polygon([[bbox[0], bbox[1]], [bbox[2], bbox[1]], [bbox[2], bbox[3]], [bbox[0], bbox[3]]]) bbox_gpd = gpd.GeoDataFrame(geometry=[test_bbox]) bbox_gpd.crs = {'init' :'epsg:27700'} bbox_gpd.to_file(os.path.join(path,'bbox.shp')) return bbox
[docs]def concatenate_arrays(myarray): """ """ lenOne=0 for i in range(len(myarray)): if len(myarray[i]) == 1: lenOne+=1 for i in range(len(myarray)-lenOne): if (len(myarray[i]) == 1): myarray[i+1] = np.concatenate((myarray[i],myarray[i+1])) del myarray[i] lastIndex = len(myarray)-1 if len(myarray[lastIndex]) == 1: myarray[lastIndex-1] = np.concatenate((myarray[lastIndex-1],myarray[lastIndex])) del myarray[lastIndex] return myarray
[docs]def slope(X,Y): """ Calculate the slope as ratio 1:X meters and as percentage X: ndarray of the end distance of the segment with x gradient Y: ndarray of the gradient values of each segment Returns ndarray of ratio and pct values of slope """ diff_ = abs(np.diff(np.array(X)))/np.diff(np.array(Y)) ratio = np.nan_to_num(diff_,0) pct = (1/ratio) * 100 return ratio,pct
[docs]def distance_mask(dists, lim_dist, lim_n): """ Extract subsets of a list of distances, given a limit threshold. Used for """ groups = {'attached':[],'detached':[],'mask':[]} tmp = [0] p = 0 for i,(d0,d1) in enumerate(zip(dists[:-1],dists[1:])): if d0>lim_dist and d1>lim_dist: # above lim p+=1 elif d0<=lim_dist and d1<=lim_dist: # below lim p+=1 elif d0<=lim_dist and d1>lim_dist: # detach from alignment groups['attached'].append(tmp+[i]) groups['mask'].append(0) p=0 tmp=[i+1] else: # d0>lim_dist and d1<=lim_dist: # attach to alignment if p > lim_n: #check that we are re-attaching a subset longer than lim_n groups['detached'].append(tmp+[i]) groups['mask'].append(1) p=0 tmp=[i+1] else: # skip this detached group and add it to the attached last = groups['attached'][-1] delta = last[1] - last[0] tmp = [last[0]] del groups['attached'][-1] del groups['mask'][-1] p = p + delta if i == len(dists)-2: if d0<lim_dist and d1<lim_dist: groups['attached'].append(tmp+[i+1]) groups['mask'].append(0) else: groups['detached'].append(tmp+[i+1]) groups['mask'].append(1) return groups
[docs]def pick_from_arrays(a1, a2, mask): """ Sukhbinder 5 April 2017 Based on: https://github.com/sukhbinder/intersection/blob/master/intersection.py """ result=[] count_1 = 0 count_2 = 0 for m in mask: if m: result.append(a1[count_1]) count_1+=1 else: result.append(a2[count_2]) count_2+=1 return result
def _rect_inter_inner(x1,x2): """ """ n1=x1.shape[0]-1 n2=x2.shape[0]-1 X1=np.c_[x1[:-1],x1[1:]] X2=np.c_[x2[:-1],x2[1:]] S1=np.tile(X1.min(axis=1),(n2,1)).T S2=np.tile(X2.max(axis=1),(n1,1)) S3=np.tile(X1.max(axis=1),(n2,1)).T S4=np.tile(X2.min(axis=1),(n1,1)) return S1,S2,S3,S4 def _rectangle_intersection_(x1,y1,x2,y2): """ """ S1,S2,S3,S4=_rect_inter_inner(x1,x2) S5,S6,S7,S8=_rect_inter_inner(y1,y2) C1=np.less_equal(S1,S2) C2=np.greater_equal(S3,S4) C3=np.less_equal(S5,S6) C4=np.greater_equal(S7,S8) ii,jj=np.nonzero(C1 & C2 & C3 & C4) return ii,jj
[docs]def intersection(x1,y1,x2,y2): """ Intersections of curves.Computes the (x,y) locations where two curves intersect. The curves can be broken with NaNs or have vertical segments. Usage ----- x,y=intersection(x1,y1,x2,y2) Example ------- a, b = 1, 2 phi = np.linspace(3, 10, 100) x1 = a*phi - b*np.sin(phi) y1 = a - b*np.cos(phi) x2=phi y2=np.sin(phi)+2 x,y=intersection(x1,y1,x2,y2) plt.plot(x1,y1,c='r') plt.plot(x2,y2,c='g') plt.plot(x,y,'*k') plt.show() """ ii,jj=_rectangle_intersection_(x1,y1,x2,y2) n=len(ii) dxy1=np.diff(np.c_[x1,y1],axis=0) dxy2=np.diff(np.c_[x2,y2],axis=0) T=np.zeros((4,n)) AA=np.zeros((4,4,n)) AA[0:2,2,:]=-1 AA[2:4,3,:]=-1 AA[0::2,0,:]=dxy1[ii,:].T AA[1::2,1,:]=dxy2[jj,:].T BB=np.zeros((4,n)) BB[0,:]=-x1[ii].ravel() BB[1,:]=-x2[jj].ravel() BB[2,:]=-y1[ii].ravel() BB[3,:]=-y2[jj].ravel() for i in range(n): try: T[:,i]=np.linalg.solve(AA[:,:,i],BB[:,i]) except: T[:,i]=np.NaN in_range= (T[0,:] >=0) & (T[1,:] >=0) & (T[0,:] <=1) & (T[1,:] <=1) xy0=T[2:,in_range] xy0=xy0.T return xy0[:,0],xy0[:,1]
[docs]def line_intersection(line1, line2): """ Intersection between two lines Parameters ---------- line1 : array-like give two points as ((A_x,A_y),(B_x,B_y)) line2 : array-like give two points as ((A_x,A_y),(B_x,B_y)) Returns ------- tuple intersection point as (x,y) ref: https://stackoverflow.com/questions/20677795/how-do-i-compute-the-intersection-point-of-two-lines """ assert(len(line1[0]) == len(line2[0])), "lines does not share the same coordinate system dimension" if len(line1[0]) is not 2: # in case of 3d point is_3d1 = sum([line1[i][2] for i in range(len(line1))]) != 0. is_3d2 = sum([line2[i][2] for i in range(len(line2))]) != 0. if is_3d1: print("WARNING: line1: {} is 3d, will be projected on XY".format(line1)) if is_3d2: print("WARNING: line2: {} is 3d, will be projected on XY".format(line2)) try: line1 = [pt[:2] for pt in line1] line1 = [pt[:2] for pt in line1] except Exception as e: print("ERROR: line projection to XY failed, exception: {}".format(e)) xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0]) ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1]) def det(a, b): return a[0] * b[1] - a[1] * b[0] div = det(xdiff, ydiff) if div == 0: raise Exception('lines do not intersect') d = (det(*line1), det(*line2)) x = det(d, xdiff) / div y = det(d, ydiff) / div return (x, y)
[docs]def ellipse_pts(): """ (x-1)^2 + (y-1)^2 + 2cxy = 1 c = 1/2 ref1: https://math.stackexchange.com/questions/252368/defining-constructing-an-ellipse ref2: https://math.stackexchange.com/questions/2645689/what-is-the-parametric-equation-of-a-rotated-ellipse-given-the-angle-of-rotatio ref3: go to wolfram alpha, insert the equation and click on properties for the centre and radii """ # parameter s=0.1 #step alpha = np.arange(0,2*np.pi+s,s) # rotation angle t = (3./4.)*np.pi # centre cx = 0.666666667 cy = 0.666666667 # radii r_x = 0.816497 # major r_y = 0.471405 # minor xe = [r_x*math.cos(a)*math.cos(t) - r_y*math.sin(a)*math.sin(t) + cx for a in alpha] ye = [r_x*math.cos(a)*math.sin(t) + r_y*math.sin(a)*math.cos(t) + cy for a in alpha] return xe,ye
[docs]def circle_pts(): """ (x-1)^2 + (y-1)^2 = 1 Returns ------- array list of points of the circle sector (0 to pi/2/) """ cx,cy = 1.,1. r=1. theta = np.arange(np.pi,np.pi+np.pi*0.5,0.1) xe = [r*math.cos(t)+cx for t in theta] ye = [r*math.sin(t)+cy for t in theta] return xe, ye
[docs]def recover_homogenous_affine_transformation(p, p_prime): ''' Find the unique homogeneous affine transformation that maps a set of 3 points to another set of 3 points in 3D space: p_prime == np.dot(p, R) + t where `R` is an unknown rotation matrix, `t` is an unknown translation vector, and `p` and `p_prime` are the original and transformed set of points stored as row vectors: p = np.array((p1, p2, p3)) p_prime = np.array((p1_prime, p2_prime, p3_prime)) The result of this function is an augmented 4-by-4 matrix `A` that represents this affine transformation: np.column_stack((p_prime, (1, 1, 1))) == \ np.dot(np.column_stack((p, (1, 1, 1))), A) Source: https://math.stackexchange.com/a/222170 (robjohn) ''' # construct intermediate matrix Q = p[1:] - p[0] Q_prime = p_prime[1:] - p_prime[0] # calculate rotation matrix R = np.dot(np.linalg.inv(np.row_stack((Q, np.cross(*Q)))), np.row_stack((Q_prime, np.cross(*Q_prime)))) # calculate translation vector t = p_prime[0] - np.dot(p[0], R) # calculate affine transformation matrix return np.column_stack((np.row_stack((R, t)), (0, 0, 0, 1)))
[docs]def transform_pt(point, trans_mat): """ Apply affine matrix to points Parameters ---------- points : array-like array of points as [x,y,z] trans_mat : array-like affine transformation matrix. Use `recover_homogenous_affine_transformation` from this module to retrieve a transformation matrix from two sets of three points """ a = np.array([point[0], point[1], point[2], 1]) ap = np.dot(a, trans_mat)[:3] return [ap[0], ap[1], ap[2]]
[docs]def apply_circle(p1, p2, p3): """ Apply affine transformation to a circle sector to pass by three points p1 : array-like point on alignment 1 p2 : array-like intersection point of tangents p1 : array-like point on alignment 2 Returns ------- array array of points """ p = np.array([[1.,0.,0.],[0.,0.,0.],[0.,1.,0.]]) p_prime = np.array([p1,p2,p3]) mat = recover_homogenous_affine_transformation(p, p_prime) # calc ellipse x,y = circle_pts() pts = np.array([x, y, [0.]*len(x) ]).transpose() pts_ = [transform_pt(pt,mat) for pt in pts] return np.array(pts_)