Skip to content
GitLab

Source

Target

Commits (3)
Hide whitespace changes
Inline Side-by-side
Showing
with 1576 additions and 141 deletions
+1576 -141
__pycache__/project.cpython-38.pyc
View file @ 89c04e24
No preview for this file type
__pycache__/transform.cpython-37.pyc
View file @ 89c04e24
No preview for this file type
__pycache__/transform.cpython-38.pyc
View file @ 89c04e24
No preview for this file type
__pycache__/velocities.cpython-37.pyc
View file @ 89c04e24
No preview for this file type
__pycache__/velocities.cpython-38.pyc
View file @ 89c04e24
No preview for this file type
orthocyd.py 0 → 100755
View file @ 89c04e24
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 11 16:37:16 2022
################################
Curv2ortho (for LSPIV)
################################
Python library for ortho-rectification of curved surfaces
@author: guillaume.bodart
"""
import numpy as np
class polynom():
"""
Defines a polynom
Contains
---------------
- lstCoeff [npArray]
Coefficients of the polynom
- lim [tuple]
"""
def __init__(self,coeffs = np.zeros(0), lim = None):
if not isinstance(coeffs, np.ndarray):
print("\nERROR: coeffs must be a numpy array")
return 1
if not (isinstance(lim, tuple) or isinstance(lim, type(None))):
print("\nERROR: lim must be a tuple")
return 1
self.lstCoeff = coeffs
self.lim = lim
def set_coeff(self,coeffs,lim):
"""
Set the coefficient of the polynom
Parameters
----------
coeffs : npArray
Array containing the coefficients of the polynom sorted in descending order.
Returns
-------
errno
0 - OK
1 - error with the input.
"""
if not isinstance(coeffs, np.ndarray):
print("\nERROR: coeffs must be a numpy array")
return 1
if not isinstance(lim, tuple):
print("\nERROR: lim must be a tuple")
return 1
# Add coeff
self.lstCoeff = coeffs
self.lim = lim
return 0
def compute_Y(self,X):
"""
Compute the image of X by the polynomial function
Parameters
----------
X : float
Entry.
Returns
-------
out
Results.
"""
if len(self.lstCoeff)!=0:
if (isinstance(X,np.ndarray)):
# If X is an array or list, then adjust according to the limits
if isinstance(self.lim,tuple):
X = X[X > self.lim[0]]
X = X[X <= self.lim[1]]
else:
if isinstance(self.lim,tuple):
if X < self.lim[0] or X > self.lim[1]:
return np.nan
out = 0
for ind,el in enumerate(self.lstCoeff[::-1]):
out = out + el*X**(ind)
return out
else:
return 0
def compute_Yderive(self,X):
"""
Compute the image of X by the derivative of the polynomial function
Parameters
----------
X : float
Entry.
Returns
-------
out
Results.
"""
if len(self.lstCoeff)!=0:
# Compute the derivative
lstDeriv = np.zeros(0)
for ind,el in enumerate(self.lstCoeff[::-1]):
lstDeriv = np.append(lstDeriv,el*ind)
lstDeriv = lstDeriv[1:]
out = 0
for ind,el in enumerate(lstDeriv):
out = out + el*X**(ind)
return out
else:
return 0
def lenAB(poly,X):
"""
Compute the length of the polynomial function 'f' within the x sequence.
The function is defined as y=f(x). The step corresponds to the x sampling distance
Parameters
----------
poly : polynome()
polynomial function describing the curve of the flow surface.
X : npArray
x sequence used.
Returns
-------
float
length of the function whithin the given x sequence.
"""
s = 0
step = np.mean(X[1:]-X[:-1])
for el in X:
s = s + np.sqrt(1+poly.compute_Yderive(el)**2)*step
return np.abs(s)
def computeLongMesh(poly,X,nbV):
"""
Compute the longitudinal mesh according to the polynomial equation given and the x sequence.
nbV corresponds to the number of points along the curve (i.e. the number of pixels on the i axis of the ortho-rectified image)
Parameters
----------
poly : polynome
polynomial function describing the curve of the flow surface.
X : npArray
x sequence considered.
nbV : int
Number of points along the curve (i.e. the number of pixels on the i axis of the ortho-rectified image) .
Returns
-------
xy : TYPE
DESCRIPTION.
"""
# Create position vector
xy = np.empty((nbV,2))
# Compute the mesh step size
d = lenAB(poly,X)/nbV
# Compute first point
xy[0,:] = X[0],poly.compute_Y(X[0])
k = 1
for el in xy[:-1]:
# Compute vector director
xD = 1
yD = poly.compute_Yderive(el[0])
nD = np.sqrt(xD**2+yD**2)
xD = xD/nD
yD = yD/nD
# Compute position x2
x2 = el[0]+(d*xD)
y2 = el[1]+(d*yD)
xy[k,:] = x2,y2
k = k+1
return xy
class function():
"""
Defines a function made of polynoms
Contains
---------------
- poly [polynom()]
List of polynom() (Class containing the polynom)
- lim [tuple]
List containing the limits of each polynom
"""
def __init__(self):
self.poly = []
self.lim = []
def add_poly(self,poly):
if not isinstance(poly,polynom):
print("\nERROR: poly must be a polynom()\n")
return 1
self.poly.append(poly)
self.lim.append(poly.lim)
return 0
def compute_Y(self,X):
if isinstance(X,np.ndarray):
out = np.empty(0)
for el in self.poly:
out = np.concatenate((out,el.compute_Y(X)))
return out
else:
for el in self.poly:
if np.isfinite(el.compute_Y(X)):
out = el.compute_Y(X)
return out
def lenAB(self,X):
"""
Compute the length of the polynomial function 'f' within the x sequence.
The function is defined as y=f(x). The step corresponds to the x sampling distance
Parameters
----------
poly : polynome()
polynomial function describing the curve of the flow surface.
X : npArray
x sequence used.
Returns
-------
float
length of the function whithin the given x sequence.
"""
s = 0
step = np.mean(X[1:]-X[:-1])
for poly in self.poly:
for el in X:
if np.isfinite(poly.compute_Yderive(el)): #check out nan values
s = s + np.sqrt(1+poly.compute_Yderive(el)**2)*step
return np.abs(s)
def computeLongMesh(self,X,nbV):
"""
Compute the longitudinal mesh according to the polynomial equation given and the x sequence.
nbV corresponds to the number of points along the curve (i.e. the number of pixels on the i axis of the ortho-rectified image)
Parameters
----------
poly : polynome
polynomial function describing the curve of the flow surface.
X : npArray
x sequence considered.
nbV : int
Number of points along the curve (i.e. the number of pixels on the i axis of the ortho-rectified image) .
Returns
-------
xy : TYPE
DESCRIPTION.
"""
# Create position vector
xy = np.empty((nbV,2))
# Compute the mesh step size
d = self.lenAB(X)/nbV
# Compute first point
xy[0,:] = X[0],self.compute_Y(X[0])
k = 1
for el in xy[:-1]:
# Compute vector director
xD = 1
yD = self.compute_Yderive(el[0])
nD = np.sqrt(xD**2+yD**2)
xD = xD/nD
yD = yD/nD
# Compute position x2
x2 = el[0]+(d*xD)
y2 = el[1]+(d*yD)
xy[k,:] = x2,y2
k = k+1
return xy
project.py
View file @ 89c04e24
......@@ -16,13 +16,16 @@ import os
from tkinter import Tk
from tkinter.filedialog import askopenfilename
from zipfile import ZipFile
import glob
from lspiv import velocities,transform
from tqdm import tqdm
from tkinter.filedialog import askdirectory
import glob
from PIL import Image
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.backends.backend_pdf import PdfPages
import scipy.stats
class lspivProject():
"""
......@@ -137,12 +140,17 @@ class lspivProject():
self.PIV_param = velocities.pivParam() # PIV parameters
self.vel_raw = [] # list of class vel()
self.vel_filtre = [] # list of class vel()
self.vel_filtre_pp = []
self.vel_raw_pp = []
self.vel_moy = velocities.vel() # list of average velocities
self.contains = {} # Dictionnary resuming files in the project
self.img_src = [] # List of source images
self.img_transf = [] # List of transformed images
self.mask = [] # Mask of river. Needed for several routines
self.transect = [] # List of bathy()
self.Q = {}
self.vel_ens = [] # Ensemble velocity
self.rGrid = [] # correlation field (ensemble)
def readReport(self,filename = ""):
"""
......@@ -193,17 +201,17 @@ class lspivProject():
q_r = df.iat[5,15]
self.q_r = float(q_r[:-1])
# Array with results / transect
self.q_mes = np.zeros(int(df.iat[30,6]))
for k in range(int(df.iat[30,6])):
self.q_mes[k] = float(df.iat[62+k,15])
# U
self.u_m = np.zeros(int(df.iat[30,6]))
for k in range(int(df.iat[30,6])):
self.u_m[k] = float(df.iat[62+k,11])
# Qtot
self.q_tot = np.zeros(int(df.iat[30,6]))
for k in range(int(df.iat[30,6])):
self.q_tot[k] = float(df.iat[62+k,3])
# self.q_mes = np.zeros(int(df.iat[30,6]))
# for k in range(int(df.iat[30,6])):
# self.q_mes[k] = float(df.iat[62+k,15])
# # U
# self.u_m = np.zeros(int(df.iat[30,6]))
# for k in range(int(df.iat[30,6])):
# self.u_m[k] = float(df.iat[62+k,11])
# # Qtot
# self.q_tot = np.zeros(int(df.iat[30,6]))
# for k in range(int(df.iat[30,6])):
# self.q_tot[k] = float(df.iat[62+k,3])
# nb im
nb_img = df.iat[11,6]
self.nb_img = int(nb_img)
......@@ -271,8 +279,8 @@ class lspivProject():
self.PIV_param.alphaCoeff = float(df.iat[33,11])
self.PIV_param.rxyMax =(float(df.iat[32,11][:-1]),float(df.iat[32,15][:-1]))
self.PIV_param.dxp =(float(df.iat[31,11][:-1]),float(df.iat[31,15][:-1]))
self.PIV_param.rxyMax =(float(df.iat[31,11][:-1]),float(df.iat[31,15][:-1]))
self.PIV_param.dxp =(float(df.iat[32,11][:-1]),float(df.iat[32,15][:-1]))
self.notes = df.iat[96,1]
......@@ -480,19 +488,21 @@ class lspivProject():
dictFiles["vel_raw"] = len(lstraw) # Add entry to follow-up checking dictionnary
# Extract data, field by field
for fich in tqdm (lstraw,desc="Raw Velocities"):
buff = foldzip.read(fich)
tmp = buff.split()
data = np.reshape(np.asarray(tmp,dtype = float),(-1,5)) # Shaping data
buff = foldzip.read(fich).decode()
sLi = len((buff.split("\n")[0].split()))
tmp = buff.split()
data = np.reshape(np.asarray(tmp,dtype = float),(-1,sLi)) # Shaping data
# Create the vel() element containing the datas
elVelRaw = velocities.vel()
# Get data from extraction
elVelRaw.i = data[:,0].copy()# copy the values -> avoid references troubles
elVelRaw.i = self.PIV_param.ninj[0] - elVelRaw.i # correct i ax
#elVelRaw.i = self.PIV_param.ninj[0] - elVelRaw.i # correct i ax
elVelRaw.j = data[:,1].copy()
elVelRaw.di = data[:,2].copy()
elVelRaw.di = -elVelRaw.di # correct i ax
elVelRaw.dj = data[:,3].copy()
elVelRaw.r = data[:,4].copy()
elVelRaw.Ar = data[:,7].copy()
elVelRaw.vi = elVelRaw.di / self.PIV_param.dt
elVelRaw.vj = elVelRaw.dj / self.PIV_param.dt
# Compute other quantities
......@@ -545,7 +555,550 @@ class lspivProject():
self.contains = dictFiles
return 0
def display_raw(self,wd_p,Fname="/scatter_vxvy_raw.pdf",en=True):
"""
Scatter plot of raw velocities (all time step and all grid points)
Parameters
----------
wd_p : str
path to save file (without / at the end).
Fname : str , optional
File name. The default is "/scatter_vxvy_raw.pdf".
en : bool, optional
Figure's title in english (If false -> French). The default is True.
Returns
-------
int
errno.
"""
with PdfPages(wd_p + Fname) as pdf:
# Set cmap
cmap = matplotlib.cm.get_cmap("magma")
vx,vy,r = [],[],[]
for el in self.vel_raw:
vx = np.concatenate((vx,el.vx))
vy = np.concatenate((vy,el.vy))
r = np.concatenate((r,el.r))
cond = np.logical_and(np.isfinite(vx),np.isfinite(r))
cond = np.logical_and(cond,r!=-99.)
vx = vx[cond]
vy = vy[cond]
r = r[cond]
fig = plt.figure(figsize=(8.2, 8))
# Add a gridspec with two rows and two columns and a ratio of 2 to 7 between
# the size of the marginal axes and the main axes in both directions.
# Also adjust the subplot parameters for a square plot.
gs = fig.add_gridspec(2, 2, width_ratios=(7, 2), height_ratios=(2, 7),
left=0.1, right=0.9, bottom=0.1, top=0.9,
wspace=0.05, hspace=0.05)
ax = fig.add_subplot(gs[1, 0])
ax_histx = fig.add_subplot(gs[0, 0], sharex=ax)
ax_histy = fig.add_subplot(gs[1, 1], sharey=ax)
# no labels
ax_histx.tick_params(axis="x", labelbottom=False)
ax_histy.tick_params(axis="y", labelleft=False)
# the scatter plot:
ax.scatter(vx,
vy,
alpha = 0.35,
c = r,
cmap = "magma",
edgecolor = None,
linewidth = 0.1,
vmin = 0,
vmax = 1)
# Histograms
# Ensure one bin of the histogram represent one pixel
# -> useful to have an estimation of the displacement in pixels !
hx = ax_histx.hist(vx,
bins= self.PIV_param.SA["Sjp"]+self.PIV_param.SA["Sjm"],
linewidth = 0.5,density=True)
hy = ax_histy.hist(vy,
bins=self.PIV_param.SA["Sip"]+self.PIV_param.SA["Sim"],
orientation='horizontal',
linewidth = 0.5,density=True)
# Add colors on histograms
for ind,xl in enumerate(hx[1][:-1]):
xlm = hx[1][ind+1]
# Boucle sur les vitesses
r_ = []
for k,vv in enumerate(vx):
if vv >= xl and vv <= xlm:
r_.append(r[k])
rm = np.mean(r_)
# Apply color on hist
hx[2][ind]._facecolor = cmap(rm)
for ind,yl in enumerate(hy[1][:-1]):
ylm = hy[1][ind+1]
# Boucle sur les vitesses
r_ = []
for k,vv in enumerate(vy):
if vv >= yl and vv <= ylm:
r_.append(r[k])
rm = np.mean(r_)
# Apply color on hist
hy[2][ind]._facecolor = cmap(rm)
# Arrange axes
ax.grid(which='minor', linestyle=':', linewidth=0.5)
ax.minorticks_on()
ax.grid(which='major', linestyle='-', linewidth=0.75)
ax.tick_params(which='minor', bottom=True, left=True)
clb = matplotlib.cm.ScalarMappable(matplotlib.colors.Normalize(vmin=0, vmax=1),cmap = "magma")
if en:
cc = plt.colorbar(clb,label = "Correlation")
else:
cc = plt.colorbar(clb,label = "Corrélation")
# Compute data for title
IAm = np.round(self.PIV_param.IA * self.img_ref.res,3)
if en:
ax_histx.set_title("Data: $v_{raw}$"+" - IA: {} pix $\longleftrightarrow$ {} m - framestep: {}".format(self.PIV_param.IA,IAm,self.ind2))
else:
ax_histx.set_title("Données: $v_{brutes}$"+" - IA: {} pix $\longleftrightarrow$ {} m - framestep: {}".format(self.PIV_param.IA,IAm,self.ind2))
ax.annotate("SA :" ,(420,480),xycoords = "figure points",fontsize = 10)
ax.annotate("$Sim =" + str(self.PIV_param.SA["Sim"])+"pix$" ,(420,465),xycoords = "figure points",fontsize = 10)
vym = (self.PIV_param.SA["Sim"] * self.img_ref.res)/self.PIV_param.dt
ax.annotate(" $Sym =-" + str(np.round(vym,3)) +"m/s$" ,(500,465),xycoords = "figure points",fontsize = 10)
ax.annotate("$Sip =" + str(self.PIV_param.SA["Sip"])+"pix$" ,(420,450),xycoords = "figure points",fontsize = 10)
vyp = (self.PIV_param.SA["Sip"] * self.img_ref.res)/self.PIV_param.dt
ax.annotate(" $Syp =" + str(np.round(vyp,3)) +"m/s$" ,(500,450),xycoords = "figure points",fontsize = 10)
ax.annotate("$Sjm =" + str(self.PIV_param.SA["Sjm"])+"pix$" ,(420,435),xycoords = "figure points",fontsize = 10)
vxm = (self.PIV_param.SA["Sjm"] * self.img_ref.res)/self.PIV_param.dt
ax.annotate(" $Sxm =-" + str(np.round(vxm,3)) +"m/s$" ,(500,435),xycoords = "figure points",fontsize = 10)
ax.annotate("$Sjp =" + str(self.PIV_param.SA["Sjp"])+"pix$",(420,420),xycoords = "figure points",fontsize = 10)
vxp = (self.PIV_param.SA["Sjp"] * self.img_ref.res)/self.PIV_param.dt
ax.annotate(" $Sxp =" + str(np.round(vxp,3)) +"m/s$" ,(500,420),xycoords = "figure points",fontsize = 10)
ax.set_xlabel("$v_x [m/s]$")
ax.set_ylabel("$v_y [m/s]$")
pdf.savefig()
plt.close()
return 0
def display_filter(self,wd_p,Fname="/scatter_vxvy_filter.pdf",en=True):
"""
Scatter plot of raw velocities (all time step and all grid points)
Parameters
----------
wd_p : str
path to save file (without / at the end).
Fname : str , optional
File name. The default is "/scatter_vxvy_raw.pdf".
en : bool, optional
Figure's title in english (If false -> French). The default is True.
Returns
-------
int
errno.
"""
with PdfPages(wd_p + Fname) as pdf:
# Set cmap
cmap = matplotlib.cm.get_cmap("magma")
vx,vy,r = [],[],[]
for el in self.vel_filtre:
vx = np.concatenate((vx,el.vx))
vy = np.concatenate((vy,el.vy))
r = np.concatenate((r,el.r))
cond = np.logical_and(np.isfinite(vx),np.isfinite(r))
cond = np.logical_and(cond,r!=-99.)
vx = vx[cond]
vy = vy[cond]
r = r[cond]
fig = plt.figure(figsize=(8.2, 8))
# Add a gridspec with two rows and two columns and a ratio of 2 to 7 between
# the size of the marginal axes and the main axes in both directions.
# Also adjust the subplot parameters for a square plot.
gs = fig.add_gridspec(2, 2, width_ratios=(7, 2), height_ratios=(2, 7),
left=0.1, right=0.9, bottom=0.1, top=0.9,
wspace=0.05, hspace=0.05)
ax = fig.add_subplot(gs[1, 0])
ax_histx = fig.add_subplot(gs[0, 0], sharex=ax)
ax_histy = fig.add_subplot(gs[1, 1], sharey=ax)
# no labels
ax_histx.tick_params(axis="x", labelbottom=False)
ax_histy.tick_params(axis="y", labelleft=False)
# the scatter plot:
ax.scatter(vx,
vy,
alpha = 0.35,
c = r,
cmap = "magma",
edgecolor = None,
linewidth = 0.1,
vmin = 0,
vmax = 1)
# Histograms
# Ensure one bin of the histogram represent one pixel
# -> useful to have an estimation of the displacement in pixels !
hx = ax_histx.hist(vx,
bins= self.PIV_param.SA["Sjp"]+self.PIV_param.SA["Sjm"],
linewidth = 0.5,density=True)
hy = ax_histy.hist(vy,
bins=self.PIV_param.SA["Sip"]+self.PIV_param.SA["Sim"],
orientation='horizontal',
linewidth = 0.5,density=True)
# Add colors on histograms
for ind,xl in enumerate(hx[1][:-1]):
xlm = hx[1][ind+1]
# Boucle sur les vitesses
r_ = []
for k,vv in enumerate(vx):
if vv >= xl and vv <= xlm:
r_.append(r[k])
rm = np.mean(r_)
# Apply color on hist
hx[2][ind]._facecolor = cmap(rm)
for ind,yl in enumerate(hy[1][:-1]):
ylm = hy[1][ind+1]
# Boucle sur les vitesses
r_ = []
for k,vv in enumerate(vy):
if vv >= yl and vv <= ylm:
r_.append(r[k])
rm = np.mean(r_)
# Apply color on hist
hy[2][ind]._facecolor = cmap(rm)
# Arrange axes
ax.grid(which='minor', linestyle=':', linewidth=0.5)
ax.minorticks_on()
ax.grid(which='major', linestyle='-', linewidth=0.75)
ax.tick_params(which='minor', bottom=True, left=True)
clb = matplotlib.cm.ScalarMappable(matplotlib.colors.Normalize(vmin=0, vmax=1),cmap = "magma")
if en:
cc = plt.colorbar(clb,label = "Correlation")
else:
cc = plt.colorbar(clb,label = "Corrélation")
# # Compute data for title
# IAm = np.round(self.PIV_param.IA * self.img_ref.res,3)
# if en:
# ax_histx.set_title("Data: $v_{raw}$"+" - IA: {} pix $\longleftrightarrow$ {} m - framestep: {}".format(self.PIV_param.IA,IAm,self.ind2))
# else:
# ax_histx.set_title("Données: $v_{brutes}$"+" - IA: {} pix $\longleftrightarrow$ {} m - framestep: {}".format(self.PIV_param.IA,IAm,self.ind2))
# ax.annotate("SA :" ,(420,480),xycoords = "figure points",fontsize = 10)
# ax.annotate("$Sim =" + str(self.PIV_param.SA["Sim"])+"pix$" ,(420,465),xycoords = "figure points",fontsize = 10)
# vym = (self.PIV_param.SA["Sim"] * self.img_ref.res)/self.PIV_param.dt
# ax.annotate(" $Sym =-" + str(np.round(vym,3)) +"m/s$" ,(500,465),xycoords = "figure points",fontsize = 10)
# ax.annotate("$Sip =" + str(self.PIV_param.SA["Sip"])+"pix$" ,(420,450),xycoords = "figure points",fontsize = 10)
# vyp = (self.PIV_param.SA["Sip"] * self.img_ref.res)/self.PIV_param.dt
# ax.annotate(" $Syp =" + str(np.round(vyp,3)) +"m/s$" ,(500,450),xycoords = "figure points",fontsize = 10)
# ax.annotate("$Sjm =" + str(self.PIV_param.SA["Sjm"])+"pix$" ,(420,435),xycoords = "figure points",fontsize = 10)
# vxm = (self.PIV_param.SA["Sjm"] * self.img_ref.res)/self.PIV_param.dt
# ax.annotate(" $Sxm =-" + str(np.round(vxm,3)) +"m/s$" ,(500,435),xycoords = "figure points",fontsize = 10)
# ax.annotate("$Sjp =" + str(self.PIV_param.SA["Sjp"])+"pix$",(420,420),xycoords = "figure points",fontsize = 10)
# vxp = (self.PIV_param.SA["Sjp"] * self.img_ref.res)/self.PIV_param.dt
# ax.annotate(" $Sxp =" + str(np.round(vxp,3)) +"m/s$" ,(500,420),xycoords = "figure points",fontsize = 10)
# ax.set_xlabel("$v_x [m/s]$")
# ax.set_ylabel("$v_y [m/s]$")
pdf.savefig()
plt.close()
return 0
def computeAverage(self,boolMean = True):
# Create element
mean = self.vel_raw[0].copy()
# Run over grid points
for ind,el in enumerate(mean.vx):
vx,vy,dx,dy,n,r,crms,di,dj = np.empty(0),np.empty(0),np.empty(0),np.empty(0),np.empty(0),np.empty(0),np.empty(0),np.empty(0),np.empty(0)
# Extract results at point
for filtre in self.vel_filtre:
vx = np.append(vx,filtre.vx[ind])
vy = np.append(vy,filtre.vy[ind])
dx = np.append(dx,filtre.dx[ind])
dy = np.append(dy,filtre.dy[ind])
di = np.append(di,filtre.di[ind])
dj = np.append(dj,filtre.dj[ind])
n = np.append(n,filtre.n[ind])
r = np.append(r,filtre.r[ind])
#crms = np.append(crms,filtre.crms[ind])
# Remove filtered points
cond = np.logical_and(vx==0,vy ==0)
vx[cond] = np.nan
# Remove nan
di = di[np.isfinite(vx)]
dj = dj[np.isfinite(vx)]
vy = vy[np.isfinite(vx)]
dx = dx[np.isfinite(vx)]
dy = dy[np.isfinite(vx)]
n = n[np.isfinite(vx)]
r = r[np.isfinite(vx)]
vx = vx[np.isfinite(vx)]
#crms = crms[np.isfinite(crms)]
# Compute mean
if boolMean:
mean.vx[ind] = np.mean(vx)
mean.vy[ind] = np.mean(vy)
mean.dx[ind] = np.mean(dx)
mean.dy[ind] = np.mean(dy)
mean.n[ind] = np.mean(n)
mean.r[ind] = np.mean(r)
mean.di[ind] = np.mean(di)
mean.dj[ind] = np.mean(dj)
else:
mean.vx[ind] = np.median(vx)
mean.vy[ind] = np.median(vy)
mean.dx[ind] = np.median(dx)
mean.dy[ind] = np.median(dy)
mean.n[ind] = np.median(n)
mean.di[ind] = np.median(di)
mean.dj[ind] = np.median(dj)
self.vel_moy = mean
def writeAverageOut(self,wd):
"""
Write file average_vel.out in specified folder (wd).
This function is mainly used to compute discharge out of F-LSPIV as average_vel.out is used to obtain surface velocities on the transect
Parameters
----------
wd : str
Path to the folder were average_vel.out will be written.
Returns
-------
int
errno.
"""
if wd[-1] == "/":
wd = wd[:-1]
if not isinstance(self.vel_moy,list):
f = open(wd + "/average_vel.out","w")
for ind,el in enumerate(self.vel_moy.x):
x = np.round(self.vel_moy.x[ind],4)
y = np.round(self.vel_moy.y[ind],4)
vx = np.round(self.vel_moy.vx[ind],4)
vy = np.round(self.vel_moy.vy[ind],4)
li = "{} {} {} {}\n".format(x,y,vx,vy)
li.replace("nan","0.0")
f.write(li.replace("nan","0.0"))
f.close()
return 0
def writeFilterPiv(self,wd):
if wd[-1] == "/":
wd = wd[:-1]
if len(self.vel_filtre) > 0:
for ind,el in enumerate(self.vel_filtre):
f = open(wd + "/filter_piv%04d.dat" % (ind+1),"w")
for x,y,vx,vy,r in zip(el.x,el.y,el.vx,el.vy,el.r):
x = np.round(x,4)
y = np.round(y,4)
vx = np.round(vx,4)
vy = np.round(vy,4)
li = "{} {} {} {}\n".format(x,y,vx,vy)
li.replace("nan","0.0")
f.write(li.replace("nan","0.0"))
f.close()
else:
print("ERROR: no filter velocities")
return 1
return 0
def writeRawPiv(self,wd):
if wd[-1] == "/":
wd = wd[:-1]
if len(self.vel_raw) > 0:
for ind,el in enumerate(self.vel_raw):
f = open(wd + "/piv%04d.dat" % (ind+1),"w")
for x,y,vx,vy,r in zip(el.x,el.y,el.vx,el.vy,el.r):
x = np.round(x,4)
y = np.round(y,4)
vx = np.round(vx,4)
vy = np.round(vy,4)
li = "{} {} {} {}\n".format(x,y,vx,vy)
li.replace("nan","0.0")
f.write(li.replace("nan","0.0"))
f.close()
else:
print("ERROR: no filter velocities")
return 1
return 0
def genRawPp(self):
"""
Arrange data "per point".
Parameters
----------
Returns
-------
- erno [int]:
Error number
"""
if(len(self.vel_raw) > 1):
# Create list and get info
self.vel_raw_pp = []
nbPt = len(self.vel_raw[0].di)
# Initialize
for j in range(nbPt):
self.vel_raw_pp.append(velocities.vel(0))
# Run over results
for k,el in enumerate(self.vel_raw):
for j in range(nbPt):
self.vel_raw_pp[j].di = np.append(self.vel_raw_pp[j].di, el.di[j])
self.vel_raw_pp[j].dj = np.append(self.vel_raw_pp[j].dj, el.dj[j])
self.vel_raw_pp[j].vx = np.append(self.vel_raw_pp[j].vx, el.vx[j])
self.vel_raw_pp[j].vy = np.append(self.vel_raw_pp[j].vy, el.vy[j])
self.vel_raw_pp[j].n = np.append(self.vel_raw_pp[j].n, el.di[j])
self.vel_raw_pp[j].i = np.append(self.vel_raw_pp[j].i, el.i[j])
self.vel_raw_pp[j].j = np.append(self.vel_raw_pp[j].j, el.j[j])
self.vel_raw_pp[j].r = np.append(self.vel_raw_pp[j].r, el.r[j])
self.vel_raw_pp[j].x = np.append(self.vel_raw_pp[j].x, el.x[j])
self.vel_raw_pp[j].y = np.append(self.vel_raw_pp[j].y, el.y[j])
else:
print("No raw velocities, can't generate per point data...")
return 1
return 0
def genFilterPp(self):
"""
Arrange data "per point".
Parameters
----------
Returns
-------
- erno [int]:
Error number
"""
if(len(self.vel_filtre) > 1):
# Create list and get info
self.vel_filtre_pp = []
nbPt = len(self.vel_filtre[0].di)
for j in range(nbPt):
self.vel_filtre_pp.append(velocities.vel(0))
# Run over results
for k,el in enumerate(self.vel_filtre):
for j in range(nbPt):
self.vel_filtre_pp[j].di = np.append(self.vel_filtre_pp[j].di, el.di[j])
self.vel_filtre_pp[j].dj = np.append(self.vel_filtre_pp[j].dj, el.dj[j])
self.vel_filtre_pp[j].vx = np.append(self.vel_filtre_pp[j].vx, el.vx[j])
self.vel_filtre_pp[j].vy = np.append(self.vel_filtre_pp[j].vy, el.vy[j])
self.vel_filtre_pp[j].n = np.append(self.vel_filtre_pp[j].n, el.di[j])
self.vel_filtre_pp[j].i = np.append(self.vel_filtre_pp[j].i, el.i[j])
self.vel_filtre_pp[j].j = np.append(self.vel_filtre_pp[j].j, el.j[j])
self.vel_filtre_pp[j].r = np.append(self.vel_filtre_pp[j].r, el.r[j])
self.vel_filtre_pp[j].x = np.append(self.vel_filtre_pp[j].x, el.x[j])
self.vel_filtre_pp[j].y = np.append(self.vel_filtre_pp[j].y, el.y[j])
else:
print("No filtered velocities, can't generate per point data...")
return 1
return 0
def getGrid(self):
# Return grid
if len(self.grid) > 1:
return self.grid
else:
# Get grid from vel
if len(self.vel_raw) > 1:
i,j = self.vel_raw[0].i,self.vel_raw[0].j
grid = np.array([i,j]).T
return grid
else:
return None
def load_Q(self,wd):
# Get results from file Discharge.dat
try :
f = open(wd + "/Discharge.dat","r")
buff = f.read().split("\n")
f.close()
except Exception as e:
print("\nERROR: Problem opening Discharge.dat \n Exception: " + str(e))
return 1
# Extract discharge data
Q = {}
l = buff[0].split()
Q["Total"] = float(l[1])
Q["Area"] = float(l[2])
Q["Measured"] = float(l[4])
self.Q = Q
# Extract transect velocities
lst = buff[1:-1]
transect = bathy()
for el in lst:
tmp = el.split()
transect.X.append(float(tmp[2]))
transect.Y.append(float(tmp[3]))
transect.Z.append(float(tmp[4]))
transect.vx.append(float(tmp[5]))
transect.vy.append(float(tmp[6]))
transect.vn.append(float(tmp[7]))
transect.alpha.append(float(tmp[9]))
self.transect.append(transect)
return 0
def import_img_src(self,path = ""):
"""
Import source images to the lspiv project. path identifies the folder containing the images.
......@@ -750,7 +1303,7 @@ class lspivProject():
print("ERROR: no velocities found in vel_raw !")
return 1
# Create list and get info
vel_filter_pp = []
vel_filtre_pp = []
nbIt = len(self.vel_raw) # Number of intervals (or frame/step)
nbPt = len(self.vel_raw[0].di) # Number of points on grid
# Initialization of elements
......@@ -760,28 +1313,205 @@ class lspivProject():
# Fill in elements of list
for k,el in enumerate(self.vel_raw):
for j in range(nbPt):
vel_filter_pp[j].di[k] = el.di[j]
vel_filter_pp[j].dj[k] = el.dj[j]
vel_filter_pp[j].dx[k] = el.dx[j]
vel_filter_pp[j].dy[k] = el.dy[j]
vel_filter_pp[j].i[k] = el.i[j]
vel_filter_pp[j].j[k] = el.j[j]
vel_filter_pp[j].n[k] = el.n[j]
vel_filter_pp[j].nij[k] = el.nij[j]
vel_filter_pp[j].r[k] = el.r[j]
vel_filter_pp[j].vi[k] = el.vi[j]
vel_filter_pp[j].vj[k] = el.vj[j]
vel_filter_pp[j].vx[k] = el.vx[j]
vel_filter_pp[j].vy[k] = el.vy[j]
vel_filter_pp[j].theta[k] = el.theta[j]
vel_filter_pp[j].x[k] = el.x[j]
vel_filter_pp[j].y[k] = el.y[j]
vel_filtre_pp[j].di[k] = el.di[j]
vel_filtre_pp[j].dj[k] = el.dj[j]
vel_filtre_pp[j].dx[k] = el.dx[j]
vel_filtre_pp[j].dy[k] = el.dy[j]
vel_filtre_pp[j].i[k] = el.i[j]
vel_filtre_pp[j].j[k] = el.j[j]
vel_filtre_pp[j].n[k] = el.n[j]
vel_filtre_pp[j].nij[k] = el.nij[j]
vel_filtre_pp[j].r[k] = el.r[j]
vel_filtre_pp[j].vi[k] = el.vi[j]
vel_filtre_pp[j].vj[k] = el.vj[j]
vel_filtre_pp[j].vx[k] = el.vx[j]
vel_filtre_pp[j].vy[k] = el.vy[j]
vel_filtre_pp[j].theta[k] = el.theta[j]
vel_filtre_pp[j].x[k] = el.x[j]
vel_filtre_pp[j].y[k] = el.y[j]
# Add list to object lspivProject
self.vel_filter_pp = vel_filter_pp.copy()
self.vel_filtre_pp = vel_filtre_pp.copy()
return 0
def temporalFiltering(self,dest,k=3,tCoV = 0.4, tR = 0.25,CoV = True, dev = True):
if self.transect.a != 0 or self.transect.b !=0:
a,b = self.transect.a, self.transect.b
transect = True
if dest == "raw":
# Create filter velocities
self.vel_filtre = []
for el in self.vel_raw:
self.vel_filtre.append(el.copy())
else:
self.genFilterPp()
varA_,varA_2,fCoV_ = [],[],[]
# Run over grid points
for ind,el in tqdm(enumerate(self.vel_filtre_pp),desc="Temporal filtering..."):
di = el.di
dj = el.dj
# Force non measured points (where di==0 and dj==0) to nan before mean computation
cond = np.logical_or(di==0,dj ==0)
di[cond] = np.nan
dj[cond] = np.nan
if dev:
# ## FILTER 2 - DEVIATION FROM MEAN VELOCITY
# ##################################
# Evaluate mean and std at each point for di/dj
meanI = np.nanmean(di)
stdI = np.nanstd(di)
meanJ = np.nanmean(dj)
stdJ = np.nanstd(dj)
# Temporal filter (k*sigma tolerance)
mnI,mxI= meanI-(k*stdI),meanI+(k*stdI)
mnJ,mxJ= meanJ-(k*stdJ),meanJ+(k*stdJ)
# Apply filter
fI = np.logical_or(di>mxI,di<mnI)
fJ = np.logical_or(dj>mxJ,dj<mnJ)
stpI = np.where(fI==True)
stpJ = np.where(fJ==True)
# Apply on pp velocities
el.di[stpI] = np.nan
el.dj[stpI] = np.nan
el.r[stpI] = np.nan
el.vx[stpI] = np.nan
el.vy[stpI] = np.nan
el.di[stpJ] = np.nan
el.dj[stpJ] = np.nan
el.r[stpJ] = np.nan
el.vx[stpJ] = np.nan
el.vy[stpJ] = np.nan
for stp in stpI[0]:
self.vel_filtre[stp].vx[ind] = np.nan
self.vel_filtre[stp].vy[ind] = np.nan
self.vel_filtre[stp].dx[ind] = np.nan
self.vel_filtre[stp].dy[ind] = np.nan
self.vel_filtre[stp].di[ind] = np.nan
self.vel_filtre[stp].dj[ind] = np.nan
self.vel_filtre[stp].r[ind] = np.nan
for stp in stpJ[0]:
self.vel_filtre[stp].vx[ind] = np.nan
self.vel_filtre[stp].vy[ind] = np.nan
self.vel_filtre[stp].dx[ind] = np.nan
self.vel_filtre[stp].dy[ind] = np.nan
self.vel_filtre[stp].di[ind] = np.nan
self.vel_filtre[stp].dj[ind] = np.nan
self.vel_filtre[stp].r[ind] = np.nan
# ## FILTER 2 - COEFF. OF VARIATION
# ##################################
vx = el.vx
vy = el.vy
# Force non measured points (where di==0 and dj==0) to nan before mean computation
cond = np.logical_or(vx==0,vy ==0)
vx[cond] = np.nan
vy[cond] = np.nan
if transect:
# Project velocity on vector normal to transect
vp = -b*vx+a*vy
# Compute coefficient of variation on projected velocity
varVP = np.nanstd(vp)/np.nanmean(vp)
# Compute coefficient of variation on angle
aV = np.arctan2(vy,vx)
# Store in list to compute mean and std
meanA = scipy.stats.circmean(aV,high = np.pi,low=-np.pi,nan_policy="omit")
R = scipy.stats.circvar(aV,high = np.pi,low=-np.pi,nan_policy="omit")
# Center angle around mean
# aVC = (aV-meanA)
# aVCC = []
# # Account for angle out of [-pi,pi]
# for i in aVC:
# if (i < np.pi and i> -np.pi):
# aVCC.append(i)
# elif i >= np.pi:
# aVCC.append((np.pi*2)-i)
# else:
# aVCC.append(i%(2*np.pi))
# stdA = np.nanstd(aVCC)
# Filter out all velocities at the point based on CoV
if CoV and (varVP>tCoV or R>0.25):
el.di[:] = np.nan
el.dj[:] = np.nan
el.r[:] = np.nan
el.vx[:] = np.nan
el.vy[:] = np.nan
el.di[:] = np.nan
el.dj[:] = np.nan
el.r[:] = np.nan
el.vx[:] = np.nan
el.vy[:] = np.nan
for stp in np.arange(len(di)):
self.vel_filtre[stp].vx[ind] = np.nan
self.vel_filtre[stp].vy[ind] = np.nan
self.vel_filtre[stp].dx[ind] = np.nan
self.vel_filtre[stp].dy[ind] = np.nan
self.vel_filtre[stp].di[ind] = np.nan
self.vel_filtre[stp].dj[ind] = np.nan
self.vel_filtre[stp].r[ind] = np.nan
return 0
def applyOffsetGrid(self,x0,y0):
if len(self.vel_raw) > 0:
for el in self.vel_raw:
el.x = el.x + x0
el.y = el.y + y0
if len(self.vel_filtre) > 0:
for el in self.vel_filtre:
el.x = el.x + x0
el.y = el.y + y0
if isinstance(self.vel_moy,velocities.vel) > 0:
self.vel_moy.x = self.vel_moy.x + x0
self.vel_moy.y = self.vel_moy.y + y0
class bathy():
"""
Bathymetry info
"""
def __init__(self):
self.X = []
self.Y = []
self.Z = []
self.vx = []
self.vy = []
self.vn = []
self.xp = []
self.alpha = []
self.bExtrap = []
self.a = 0
self.b = 0
self.x0 = 0
self.y0 = 0
transform.py
View file @ 89c04e24
......@@ -47,6 +47,7 @@ class camera:
self.imagesize = () # size of image (tuple)
self.a = [] # a matrix from DLT, used to convert world to image space
self.a_inv = [] # a_inv matrix, used to convert from image to world coordinates
self.xy0 = () # Position of lowest image corner in world unit
def read_a(self,filename = ""):
"""
......@@ -152,7 +153,7 @@ class camera:
return 0
def image2space(self,imgIndx,xy0 = np.nan):
def image2space(self,imgIndx,xy0 = True,ni = True):
"""
Converts image indexes to world units with information of the camera and xy0.
Assume that the input image is already ortho-rectified (i.e. image can be assumed to be a plan XY)
......@@ -201,13 +202,15 @@ class camera:
return 3
# If no xy0 are passed to the function they're calculated from the camera data
if np.isnan(xy0):
if isinstance(xy0,bool) and self.scale != -99:
xy0 = (self.position - (0.5*self.scale))[:2]
if len(imgIndx)==0:
print("WARNING : imgIndx is empty...")
# Swap column to get (j,i) as j is x and i is y
JIimgIndx = np.array([imgIndx[:,1],imgIndx[:,0]]).T
# Swap column to get (j,i) as j is x and i is y
if ni == True:
ni = np.max(imgIndx[:,0])
JIimgIndx = np.array([imgIndx[:,1],ni-imgIndx[:,0]]).T
out = (JIimgIndx/res) + xy0
return out
......@@ -461,4 +464,85 @@ def read_grid(filename = ""):
return i,j
\ No newline at end of file
def bicubic_at(ip,jp,im):
i = int(ip)
j = int(jp)
sum_fx = 0
for k in range(4):
for l in range(4):
sx = abs(j+k-2-jp)
sy = abs(i+l-2-ip)
if ((sx >= 0) and (sx < 1)):
#Cx = un - deux*(sx*sx) + (sx*sx*sx)
Cx = 1. + (sx-2.)*sx*sx
elif ((sx >= 1) and (sx < 2)):
#Cx = four - huit*sx + cinq*(sx*sx) - (sx*sx*sx)
Cx = 4. + ((5. - sx)*sx - 8.)*sx
elif (sx > 2):
Cx = 0.
if ((sy >= 0) and (sy < 1)):
#Cy = un - deux*(sy*sy) + (sy*sy*sy)
Cy = 1 + (sy-2)*sy*sy
elif ((sy >= 1) and (sy < 2)):
#!Cy = four - huit*sy + cinq*(sy*sy) - (sy*sy*sy)
Cy = 4 + ((5. -sy)*sy - 8.)*sy
elif (sy > 2):
Cy = 0.
fx = im[i+l-1,j+k-1]*Cx*Cy
sum_fx = sum_fx+fx
if (int(sum_fx) > 255):
sum_fx = 255
if (int(sum_fx) < 0):
sum_fx = 0
return int(sum_fx)
def readGRPtable(filename):
# Open GRP file
try:
f = open(filename,"r")
buff = f.read()
f.close()
except Exception as e:
print("ERROR reading file coeff.dat : " + str(e))
return 1
# Get number of GRP
i0 = buff.find("\n")
i1 = i0 + buff[i0+1:].find("\n")
nbGRP = int(buff[i0+1:i1+1])
# Read file
if0 = buff.find("j")
tmp = np.array(buff[if0+2:].split(),float)
# Re-arrange data
lstGRP = tmp.reshape((nbGRP,5))
return lstGRP
def writeGRPtable(GRP,filename,ni,nj):
buff = "GRP V2.0 {} {}\n".format(int(ni),int(nj))
buff += "{}\nX\tY\tZ\ti\tj\n".format(int(np.shape(GRP)[0]))
for el in GRP:
buff += "%.3f\t%.3f\t%.3f\t%d\t%d\n"%(el[0],el[1],el[2],int(el[3]),int(el[4]))
# Open GRP file
try:
f = open(filename,"w")
f.write(buff)
f.close()
except Exception as e:
print("ERROR reading file coeff.dat : " + str(e))
return 1
return 0
\ No newline at end of file
velocities.py
View file @ 89c04e24
......@@ -7,13 +7,13 @@
Python LSPIV tools
Contains all definitions/functions that deals with
velocities, from determination to representation
velocities, from determination to representation
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import syri.project
#import syri.project
import lspiv.transform
from tkinter.filedialog import askdirectory,askopenfilename
from tkinter import Tk
......@@ -22,11 +22,13 @@ from tqdm import tqdm
import ctypes
from numpy.ctypeslib import ndpointer
import cv2 as cv
import os
import scipy
class vel:
""" Velocity class """
def __init__(self,size = 1):
def __init__(self,size = 0):
# Image Coordinate
self.i = np.empty((size), dtype = int)
self.j = np.empty((size), dtype = int)
......@@ -46,6 +48,7 @@ class vel:
self.theta = np.empty((size), dtype = float) # Angle
# Correlation (from piv algorithm)
self.r = np.empty((size), dtype = float)
self.Ar = np.empty((size), dtype = float) # Area > 0.8r_peak
def copy(self):
......@@ -69,6 +72,10 @@ class vel:
out.theta = self.theta.copy() # Angle
# Correlation (from piv algorithm)
out.r = self.r.copy()
if hasattr(self,'Ar'):
out.Ar = self.Ar.copy()
# Crms value (contrast in IA)
#out.crms = self.crms.copy()
return out
......@@ -92,6 +99,7 @@ class vel:
self.n = np.delete(self.n,f)
self.theta = np.delete(self.theta,f)
self.r = np.delete(self.r,f)
self.crms = np.delete(self.crms,f)
return
......@@ -99,6 +107,8 @@ class vel:
""" Apply filter on velocities """
dictVel = {"di" : self.di,
"dj" : self.dj,
"x" : self.x,
"y" : self.y,
"vi" : self.vi,
"vj" : self.vj,
"nij" : self.nij,
......@@ -108,17 +118,134 @@ class vel:
"vy" : self.vy,
"n" : self.n,
"theta" : self.theta,
"r" : self.r}
"r" : self.r,
"Ar" : self.Ar}
dictFiltre = {"di" : self.di,
"dj" : self.dj,
"vi" : self.vi,
"vj" : self.vj,
"nij" : self.nij,
"dx" : self.dx,
"dy" : self.dy,
"vx" : self.vx,
"vy" : self.vy,
"n" : self.n,
"theta" : self.theta,
"r" : self.r,
"Ar" : self.Ar}
# Checkings
if not qty in ["di","dj","vi","vj","nij","dx","dy","vx","vy","n","theta","r"]:
print("ERROR: qty isn't correct. expected : 'di','dj','vi','vj','nij','dx','dy','vx','vy','n','theta','r'")
if not qty in ["Ar","x","y","di","dj","vi","vj","nij","dx","dy","vx","vy","n","theta","r"]:
print("ERROR: qty isn't correct. expected : 'di','dj','vi','vj','nij','dx','dy','vx','vy','n','theta','r','Ar','x','y'")
return 1
# Apply condition and get logical mask
f = np.logical_or(dictVel[qty] < thresMin,dictVel[qty] > thresMax)
# Apply filter to each qty
for i in dictVel:
if np.shape(dictVel[i])== np.shape(f):
dictVel[i][f] = np.nan
for i in dictFiltre:
if np.shape(dictFiltre[i])== np.shape(f):
dictFiltre[i][f] = np.nan
return 0
def medianTest(self, dist_neighbor=10,nb_neighbor = 9, epsilon=0.1, thres=2,dest = "world",treeQuery = False):
"""
Apply median test (Westerweel and Scarano 2005) to remove spurious vectors.
KD-Tree algorithm (from scipy) is used to identify the neighboors of each points.
See help of kdTree function from scipy for more info.
Parameters
----------
dist_neighbor : int, optional
The maximum distance (in [m]) to consider a point as neighbor. The default is 10.
nb_neighbor : int, optional
Defines the maximum number of points identified as neighbor (n) such as n = nb_neighbor-1.
A value of 9 ensures that 8 neighbors will be identified.The default is 9.
epsilon : float, optional
Parameter of the median test, see (Westerweel and Scarano 2005). The default is 0.1.
thres : int, optional
Parameter of the median test, see (Westerweel and Scarano 2005). The default is 2.
direction : TYPE, optional
TODO... The default is None.
Returns
-------
int
DESCRIPTION.
"""
# Arrange grid and velocities in 2D arrays
if dest == "world":
grid = np.array([self.x, self.y]).T
vel = np.array([self.vx.copy(),self.vy.copy()]).T
else:
grid = np.array([self.i, self.j]).T
vel = np.array([self.di.copy(),self.dj.copy()]).T
# Use KDTree to identify neighbors
# Optimization : the tree is passed as argument
if isinstance(treeQuery,bool):
kdTree = scipy.spatial.KDTree(grid)
res = kdTree.query(grid,nb_neighbor,distance_upper_bound = dist_neighbor)
else:
res = treeQuery
# Look for neighbors
for ptInd,distInd in zip(res[1],res[0]):
# v0, velocity at point
v0x,v0y = vel[ptInd[0],0],vel[ptInd[0],1]
# Store neighbors velocities
vx,vy = [],[]
for el in ptInd:
if el < len(grid):
vx.append(vel[el,0])
vy.append(vel[el,1])
vx = np.array(vx)
vy = np.array(vy)
# X
# ####
# Calculate ref velocity
vxr = np.nanmedian(vx)
# Calculate residual
rix = []
for el in ptInd:
if el < len(grid):
rix.append(vel[el,0]-vxr)
rix = np.array(rix)
# Normalize target residual
r0x = np.abs(v0x-vxr)/(np.nanmedian(rix)+epsilon)
# Y
# ####
# Calculate ref velocity
vyr = np.nanmedian(vy)
# Calculate residual
riy = []
for el in ptInd:
if el < len(grid):
riy.append(vel[el,1]-vyr)
riy = np.array(riy)
# Normalize target residual
r0y = np.abs(v0y-vyr)/(np.nanmedian(riy)+epsilon)
# Filter data
if r0x > thres or r0y > thres:
self.di[ptInd[0]] = np.nan
self.dj[ptInd[0]] = np.nan
if len(self.dx) > 0:
self.dx[ptInd[0]] = np.nan
self.dy[ptInd[0]] = np.nan
self.vx[ptInd[0]] = np.nan
self.vy[ptInd[0]] = np.nan
self.n[ptInd[0]] = np.nan
self.r[ptInd[0]] = np.nan
return 0
......@@ -212,6 +339,10 @@ class vel:
cmapData = "norm",
cmap = None,
scale = None,
width = 0.0045,
headwidth = 4,
lw = 0.15,
ec="k",
**param_dict):
"""
Plot the data. Option plot allow to choose the quantity to be ploted: velocity (v) or displacement (d) and the units: world (XY) or image (IJ)
......@@ -263,6 +394,24 @@ class vel:
"""
# ---- CHECK INPUTS ----
# label of colorbar
if "lblClbr" in param_dict.keys():
lblClbr = param_dict["lblClbr"]
del param_dict["lblClbr"]
else:
lblClbr = "Velocity [m/s]"
# label of colorbar
if "vmin" in param_dict.keys():
vmin = param_dict["vmin"]
del param_dict["vmin"]
else:
vmin = None
# label of colorbar
if "vmax" in param_dict.keys():
vmax = param_dict["vmax"]
del param_dict["vmax"]
else:
vmax = None
# first check on size
if len(plot)!=3:
print("\nERROR: wrong format of variable 'plot'. See help file for more information.")
......@@ -270,12 +419,15 @@ class vel:
# then check contents
mtype = plot[:1].lower() # Case sensitivity
maxes = plot[1:].lower()
if np.logical_not(mtype in ["v","d"]) or np.logical_not(maxes in ["xy","ij"]):
print("\nERROR: wrong format of variable 'plot'. See help file for more information.")
return 1
if cmapData != "norm" and len(cmapData) != len(self.x):
if not isinstance(cmapData,str) and len(cmapData) != len(self.x):
print("\nERROR : Size of cmapData doesn't match data.")
return 1
# check axes
if ax == None: # Creates figure and ax
fig = plt.figure()
......@@ -284,6 +436,7 @@ class vel:
ax.set_title(title)
ax.set_xlabel(xlbl)
ax.set_ylabel(ylbl)
if type(ax) != matplotlib.axes._subplots.Subplot:
print("\nERROR: wrong format of variable ax. Expected : matplotlib.axes._subplots.Subplot.")
return 1
......@@ -300,7 +453,11 @@ class vel:
self.vy,
color = color,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
lw = lw,
ec = ec,
**param_dict)
except ValueError as e:
print("\nERROR : " + str(e))
......@@ -314,7 +471,13 @@ class vel:
else:
if cmapData == 'norm':
cmapData = norm
elif cmapData == "r":
cmapData = self.r
try:
if vmin == None:
vmin = cmapData[np.isfinite(cmapData)].min()
if vmax == None:
vmax = cmapData[np.isfinite(cmapData)].max()
qv = ax.quiver(self.x,
self.y,
self.vx,
......@@ -322,8 +485,15 @@ class vel:
cmapData,
cmap = cmap,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
clim = (vmin,vmax),
**param_dict)
#mapp = matplotlib.cm.ScalarMappable(matplotlib.colors.Normalize(vmin = vmin,vmax = vmax),cmap = cmap)
# Add colorbar
#plt.colorbar(mapp,ax=ax,label = lblClbr)
except ValueError as e:
print("\nERROR : " + str(e))
return 2
......@@ -335,17 +505,18 @@ class vel:
return 2
else:
# Compute norm
norm = np.sqrt(self.vi**2 + self.vj**2)
norm = np.sqrt(self.vx**2 + self.vy**2)
if cmap == None:
try:
qv = ax.quiver(self.j,
self.i,
self.vj,
self.vi,
norm,
self.vx,
self.vy,
color = color,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
**param_dict)
except ValueError as e:
print("\nERROR : " + str(e))
......@@ -359,16 +530,33 @@ class vel:
else:
if cmapData == 'norm':
cmapData = norm
elif cmapData == "r":
cmapData = self.r
cmapData[cmapData == -99] = np.nan
try:
if vmin == None:
vmin = cmapData[np.isfinite(cmapData)].min()
if vmax == None:
vmax = cmapData[np.isfinite(cmapData)].max()
qv = ax.quiver(self.j,
self.i,
self.vj,
self.vi,
self.vx,
self.vy,
cmapData,
cmap = cmap,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
clim = (vmin,vmax),
**param_dict)
if vmin == None:
vmin = cmapData[np.isfinite(cmapData)].min()
if vmax == None:
vmax = cmapData[np.isfinite(cmapData)].max()
#mapp = matplotlib.cm.ScalarMappable(matplotlib.colors.Normalize(vmin = vmin,vmax = vmax),cmap = cmap)
# Add colorbar
#plt.colorbar(mapp,ax=ax,label = lblClbr)
except ValueError as e:
print("\nERROR : " + str(e))
return 2
......@@ -377,7 +565,7 @@ class vel:
return 2
except:
print("\nERROR: Problems while plotting data.")
return 2
return 2
else:
if maxes == "xy":
# Compute norm
......@@ -391,7 +579,9 @@ class vel:
norm,
color = color,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
**param_dict)
except ValueError as e:
print("\nERROR : " + str(e))
......@@ -405,7 +595,14 @@ class vel:
else:
if cmapData == 'norm':
cmapData = norm
elif cmapData == "r":
cmapData = self.r
cmapData[cmapData == -99] = np.nan
try:
if vmin == None:
vmin = cmapData[np.isfinite(cmapData)].min()
if vmax == None:
vmax = cmapData[np.isfinite(cmapData)].max()
qv = ax.quiver(self.x,
self.y,
self.dx,
......@@ -413,8 +610,15 @@ class vel:
cmapData,
cmap = cmap,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
clim = (vmin,vmax),
**param_dict)
#mapp = matplotlib.cm.ScalarMappable(matplotlib.colors.Normalize(vmin = vmin,vmax = vmax),cmap = cmap)
# Add colorbar
#plt.colorbar(mapp,ax=ax,label = lblClbr)
except ValueError as e:
print("\nERROR : " + str(e))
return 2
......@@ -435,7 +639,9 @@ class vel:
self.di,
color = color,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
**param_dict)
except ValueError as e:
print("\nERROR : " + str(e))
......@@ -449,7 +655,14 @@ class vel:
else:
if cmapData == 'norm':
cmapData = norm
elif cmapData == "r":
cmapData = self.r
cmapData[cmapData == -99] = np.nan
try:
if vmin == None:
vmin = cmapData[np.isfinite(cmapData)].min()
if vmax == None:
vmax = cmapData[np.isfinite(cmapData)].max()
qv = ax.quiver(self.j,
self.i,
self.dj,
......@@ -457,8 +670,15 @@ class vel:
cmapData,
cmap = cmap,
scale = scale,
scale_units = "xy",
scale_units = "width",
width = width,
headwidth = headwidth,
clim = (vmin,vmax),
**param_dict)
mapp = matplotlib.cm.ScalarMappable(matplotlib.colors.Normalize(vmin = vmin,vmax = vmax),cmap = cmap)
# Add colorbar
plt.colorbar(mapp,ax=ax,label = lblClbr)
except ValueError as e:
print("\nERROR : " + str(e))
return 2
......@@ -470,74 +690,74 @@ class vel:
return 2
return qv
def toDataCl(self,qty = "vXY", camera = None):
"""
Transform velocity data from lspiv module to datacl of the syri module
# def toDataCl(self,qty = "vXY", camera = None):
# """
# Transform velocity data from lspiv module to datacl of the syri module
Parameters
----------
- qty [str (len:3)], optionnal
Data to be transformed (length = 3char).
First indicates the quantity: velocity (v) or displacement (d)
After indicates the units: meters (XY) or pixels (IJ)
The default is "vXY"
# Parameters
# ----------
# - qty [str (len:3)], optionnal
# Data to be transformed (length = 3char).
# First indicates the quantity: velocity (v) or displacement (d)
# After indicates the units: meters (XY) or pixels (IJ)
# The default is "vXY"
- camera [transform.camera()], optionnal
Use camera data to apply offset on x and y values of the lspiv values.
Most of the time it is used when simple scaling is applied under fudaa-lspiv.
If None then no offset is applied
# - camera [transform.camera()], optionnal
# Use camera data to apply offset on x and y values of the lspiv values.
# Most of the time it is used when simple scaling is applied under fudaa-lspiv.
# If None then no offset is applied
Returns
-------
- out [dataCl()]
dataCl class from the syri.project lib.
In case of wrong format for qty, out = 1
# Returns
# -------
# - out [dataCl()]
# dataCl class from the syri.project lib.
# In case of wrong format for qty, out = 1
"""
if len(qty)!=3:
print("\nERROR: wrong format of variable 'plot'. See help file for more information.")
return 1
# then check contents
mtype = qty[:1].lower() # Case sensitivity
maxes = qty[1:].lower()
if np.logical_not(mtype in ["v","p"]) or np.logical_not(maxes in ["xy","ij"]):
print("\nERROR: wrong format of variable 'plot'. See help file for more information.")
return 1
# Create dataCl()
out = syri.project.dataCl()
if type(camera) == lspiv.transform.camera:
xy0 = (camera.position - (0.5*camera.scale))[:2]
else:
xy0 = np.zeros((2,1))
# Fill with values
if mtype == "v":
if maxes == "xy":
tmp = np.array([self.x + xy0[0],self.y + xy0[1]]).T
out.pos = tmp
tmp = np.array([self.vx,self.vy]).T
out.vel = tmp
out.type = "vXY_lspiv"
else:
tmp = np.array([self.i,self.j]).T
out.pos = tmp
tmp = np.array([self.vi,self.vj]).T
out.vel = tmp
out.type = "vIJ_lspiv"
else:
if maxes == "xy":
tmp = np.array([self.x + xy0[0],self.y + xy0[1]]).T
out.pos = tmp
tmp = np.array([self.dx,self.dy]).T
out.vel = tmp
out.type = "dXY_lspiv"
else:
tmp = np.array([self.i,self.j]).T
out.pos = tmp
tmp = np.array([self.di,self.dj]).T
out.vel = tmp
out.type = "dIJ_lspiv"
# """
# if len(qty)!=3:
# print("\nERROR: wrong format of variable 'plot'. See help file for more information.")
# return 1
# # then check contents
# mtype = qty[:1].lower() # Case sensitivity
# maxes = qty[1:].lower()
# if np.logical_not(mtype in ["v","p"]) or np.logical_not(maxes in ["xy","ij"]):
# print("\nERROR: wrong format of variable 'plot'. See help file for more information.")
# return 1
# # Create dataCl()
# out = syri.project.dataCl()
# if type(camera) == lspiv.transform.camera:
# xy0 = (camera.position - (0.5*camera.scale))[:2]
# else:
# xy0 = np.zeros((2,1))
# # Fill with values
# if mtype == "v":
# if maxes == "xy":
# tmp = np.array([self.x + xy0[0],self.y + xy0[1]]).T
# out.pos = tmp
# tmp = np.array([self.vx,self.vy]).T
# out.vel = tmp
# out.type = "vXY_lspiv"
# else:
# tmp = np.array([self.i,self.j]).T
# out.pos = tmp
# tmp = np.array([self.vi,self.vj]).T
# out.vel = tmp
# out.type = "vIJ_lspiv"
# else:
# if maxes == "xy":
# tmp = np.array([self.x + xy0[0],self.y + xy0[1]]).T
# out.pos = tmp
# tmp = np.array([self.dx,self.dy]).T
# out.vel = tmp
# out.type = "dXY_lspiv"
# else:
# tmp = np.array([self.i,self.j]).T
# out.pos = tmp
# tmp = np.array([self.di,self.dj]).T
# out.vel = tmp
# out.type = "dIJ_lspiv"
return out
# return out
def warpImg(self,
im1,
......@@ -738,7 +958,38 @@ class vel:
else:
return np.array(np.clip(imout,0,255),dtype = np.uint8)
def rectifVel3D(self,ni,dt,
coeffInvFile = ""):
# Get coeff files
a_inv = lspiv.transform.read_a_inv(coeffInvFile)
if isinstance(a_inv,int):
print("ERROR: Problem while extracting coeff and coeff_inv")
return 1
# APPLY TRANSFORMATION
k = 0
x_,y_ = [],[]
for i,j,di,dj in zip(self.i,self.j,self.di,self.dj):
# Find displacement in meters
self.x[k],self.y[k] = lspiv.transform.CRT2space(a_inv,ni-i,j)
xt,yt = lspiv.transform.CRT2space(a_inv,ni-(i+di),j+dj)
x_.append(xt)
y_.append(yt)
self.dx[k],self.dy[k] = xt - self.x[k], yt - self.y[k]
self.vx[k],self.vy[k] = self.dx[k] / dt, self.dy[k] / dt
k = k+1
return 0
# def write(self,name,type = "filter"):
# # Open file
# f = open(name,"w")
# if type == "filter":
class pivParam:
""" PIV parameters class """
......@@ -760,19 +1011,44 @@ class pivParam:
self.alphaCoeff = -99.99
self.rxyMax =(-99.99,-99.99)
self.dxp = -99.99
def getPivParam(self,wd_p0):
# Checking
if not isinstance(wd_p0,str):
print("ERROR: wd_p0 isn't a string")
return 1
# Ensure path is ended by a '/'
if (wd_p0[-1:]) != "/":
wd_p0 = wd_p0 + "/"
# Open file
f = open(wd_p0 + "PIV_param.dat",'r')
buf = f.read().split("\n")
f.close()
self.IA = int(buf[1])
self.SA = {"Sim" : int(buf[3]),
"Sip" : int(buf[4]),
"Sjm" : int(buf[5]),
"Sjp" : int(buf[6])}
self.dt = float(buf[8])
self.rMaxMin =(float(buf[10]),float(buf[12]))
self.ninj =(int(buf[14]) ,int(buf[15]))
self.nMaxMin =(float(buf[17]),float(buf[18])) # norm filter
self.vxMaxMin =(float(buf[20]),float(buf[21]))
self.vyMaxMin =(float(buf[23]),float(buf[24]))
self.alphaCoeff = float(buf[26])
tmp = buf[28].split()
self.rxyMax =(float(tmp[0]),float(tmp[1]))
self.dxp = float(buf[30])
def getGrid(file,outFormat = "image"):
f = open(file,"r")
buf = f.read().split()
f.close()
f = open(file[:-8] + "PIV_param.dat","r")
for i in range(14):
buff = f.readline()
nj = int(f.readline())
ni = int(f.readline())
i,j = [],[]
for ind,l in enumerate(buf):
if ind%2 == 0:
......@@ -783,9 +1059,15 @@ def getGrid(file,outFormat = "image"):
i,j = np.array(i),np.array(j)
if outFormat == "image":
f = open(file[:-8] + "PIV_param.dat","r")
for i in range(14):
buff = f.readline()
nj = int(f.readline())
ni = int(f.readline())
i = ni - i
grid = np.array([i,j]).T
elif outFormat == "fudaa":
grid = np.array([i,j]).T
return grid
......@@ -880,6 +1162,8 @@ def dynExpans(im):
"""
mx, mn = im.max(),im.min()
return 255*((im-mn)/(mx-mn))
def NCC_at(im1,im2,ii,jj,IA,Sim,Sip,Sjm,Sjp):
"""
Perform normalized cross-correlation at point (ii, jj). IA is the side of the block to be matched, Sim-ip-jm-jp defines the searching area around point ii and jj
......@@ -918,7 +1202,7 @@ def NCC_at(im1,im2,ii,jj,IA,Sim,Sip,Sjm,Sjp):
if i1 >=0 and i2+1 < ni and j1 >=0 and j2+1 < nj and i1-Sim >=0 and j1-Sjm>=0 and i2+Sip+1 < ni and j2+Sjp+1 < nj :
imgIA = im1[i1:i2+1,j1:j2+1]
imgSA = im2[i1-Sim:i2+Sip+1,j1-Sjm:j2+Sjp+1]
res = cv.matchTemplate(imgSA,imgIA,cv.TM_CCOEFF_NORMED)
res = cv.matchTemplate(np.uint8(imgSA),np.uint8(imgIA),cv.TM_CCOEFF_NORMED)
else:
res = []
......@@ -1208,4 +1492,30 @@ def cub_conv1dx(dm,n,c,a = -1):
return 0
else:
return sum_fx
\ No newline at end of file
def setSAManualVel(mVx,mVy,resol,dt):
SA = {}
if np.abs(mVx) > np.abs(mVy):
# Flow on x
if mVx >0:
SA["Sjm"] = 1
SA["Sjp"] = int(np.ceil(((2*mVx)*dt)/resol))
else:
SA["Sjm"] = int(np.ceil(((2*np.abs(mVx))*dt)/resol))
SA["Sjp"] = 1
SA["Sim"] = int(np.ceil(((2*(np.abs(mVy)+0.1))*dt)/resol))
SA["Sip"] = int(np.ceil(((2*(np.abs(mVy)+0.1))*dt)/resol))
else:
# Flow on y
if mVy >0:
SA["Sim"] = 1
SA["Sip"] = int(np.ceil(((2*mVy)*dt)/resol))
else:
SA["Sim"] = int(np.ceil(((2*np.abs(mVx))*dt)/resol))
SA["Sip"] = 1
SA["Sjm"] = int(np.ceil(((2*(np.abs(mVx)+0.1))*dt)/resol))
SA["Sjp"] = int(np.ceil(((2*(np.abs(mVx)+0.1))*dt)/resol))
return SA
\ No newline at end of file