Skip to content

GitLab

Commit e7cd5aa5 authored by Guillaume Bodart's avatar Guillaume Bodart
Browse files

MAJ latest dev

parent a29843f5
Hide whitespace changes
Inline Side-by-side
project.py
View file @ e7cd5aa5
......@@ -25,6 +25,8 @@ import matplotlib as mpl
import plotlib
import cv2 as cv
from tkinter.filedialog import askdirectory
import ctypes
from numpy.ctypeslib import ndpointer
if 'monStyle' in plt.style.available:
plt.style.use('monStyle')
......@@ -770,7 +772,19 @@ class dataCl():
lblClbr = param_dict["lblClbr"]
del param_dict["lblClbr"]
else:
lblClbr = "Velocity difference [m/s]"
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
# ---- PLOT ----
if mtype == "vel":
......@@ -809,7 +823,11 @@ class dataCl():
scale_units = "xy",
**param_dict)
mapp = mpl.cm.ScalarMappable(mpl.colors.Normalize(vmin = norm.min(),vmax = norm.max()),cmap = cmap)
if vmin == None:
vmin = norm.min()
if vmax == None:
vmax = norm.min()
mapp = mpl.cm.ScalarMappable(mpl.colors.Normalize(vmin = vmin,vmax = vmax),cmap = cmap)
# Add colorbar
plt.colorbar(mapp,ax=ax,label = lblClbr)
......@@ -1068,6 +1086,11 @@ class dataCl():
- showCdf [int], optional
Flag that indicates whether the cumulative density function (cdf) has to be plotted or not
Default is False
--- ADDITIONARY ARGUMENT FOR HIST
- legend [bool], optionnal
Indicates if legend has to be plotted, mainly for histograms and statistics.
--- ADDITIONARY ARGUMENT FOR EE_PIX :
......@@ -1104,7 +1127,9 @@ class dataCl():
camera = param_dict["camera"]
del param_dict["camera"]
else:
camera = None
if method in ["EE_pix","AE_pix","AE_2D_pix"]:
print("\nERROR: camera data is mandatory for this method ! Please add it to kwargs")
return 1
# ---- PLOT ----
dictPlotFct = {"dist" : plotlib.dist,
......@@ -1130,8 +1155,12 @@ class dataCl():
selfDsp = self.copy()
compDsp = compData.copy()
# Converts now to pixel displacement
selfDsp.vel = velSpace2dispImage(selfDsp.vel, camera)
compDsp.vel = velSpace2dispImage(compDsp.vel, camera)
try:
selfDsp.vel = camera.velSpace2dispImage(selfDsp.vel)
compDsp.vel = camera.velSpace2dispImage(compDsp.vel)
except Exception as e:
print("\nERROR: Problem while changing from velocity to pixels displacement... \n Exception: " + str(e))
return 1
# Plot
meth = method[:-4]
mplt = dictPlotFct[meth](ax,selfDsp.vel,compDsp.vel,**param_dict)
......@@ -1286,8 +1315,12 @@ class dataCl():
selfDsp = self.copy()
compDsp = compData.copy()
# Converts now to pixel displacement
selfDsp.vel = velSpace2dispImage(selfDsp.vel, camera)
compDsp.vel = velSpace2dispImage(compDsp.vel, camera)
try:
selfDsp.vel = camera.velSpace2dispImage(selfDsp.vel)
compDsp.vel = camera.velSpace2dispImage(compDsp.vel)
except Exception as e:
print("\nERROR: Problem while changing from velocity to pixels displacement... \n Exception: " + str(e))
return 1
# compute diff
ee = compDsp.vel - selfDsp.vel
# Compute norm
......@@ -1636,10 +1669,10 @@ class dataCl():
self.vel[Flspiv_data.vel == 0] = np.nan
Flspiv_data.vel[Flspiv_data.vel == 0] = np.nan
elif case ==2:
isFiltre = np.array([(Flspiv_data.vx == 0),(Flspiv_data.vy == 0)]).T
isFiltre = np.logical_and(Flspiv_data.vx == 0,Flspiv_data.vy == 0)
self.vel[isFiltre] = np.nan
Flspiv_data.vx[isFiltre[:,0]] = np.nan
Flspiv_data.vy[isFiltre[:,0]] = np.nan
Flspiv_data.vx[isFiltre] = np.nan
Flspiv_data.vy[isFiltre] = np.nan
return
......@@ -1780,7 +1813,109 @@ class dataCl():
mapp = mpl.cm.ScalarMappable(mpl.colors.Normalize(vmin = vmin,vmax = vmax),cmap = cmap)
# Add colorbar
fig.colorbar(mapp,ax=ax,label = lblClbr)
def filter(self):
""" Remove filtered velocities from list """
f = np.logical_not(np.isfinite(self.vel))
self.pos = np.array([np.delete(self.pos[:,0],f[:,0]),np.delete(self.pos[:,1],f[:,1])]).T
self.vel = np.array([np.delete(self.vel[:,0],f[:,0]),np.delete(self.vel[:,1],f[:,1])]).T
if self.normals != []:
self.normals = np.array([np.delete(self.normals[:,0],f[:,0]),np.delete(self.normals[:,1],f[:,1])]).T
return
def warpImg(self,
camera,
im1,
im2,
outRMS = True):
"""
Apply velocity field to warp image2 back into position of image1.
This tool is used for Backend Error Projection.
Output is the warped image and the RMS error computed over the velocity field
Parameters
----------
- im1 : [npArray]
First image (could be ndg or color image)
- im2 : [npArray]
Second image image (could be ndg or color image)
- outRMS : [bool]
If True then the RMS error between im1 and im1_warp is returned.
Error is computed only within the warped area.
Returns
-------
- warpedImage : [npArray]
Image 2 "warped" with the velocity field of image 1
- RMS : [float]
RMS error computed over all image
- errno : [int]
- 1 - Wrong format for inputs
"""
#### Checkings
if np.shape(im1) != np.shape(im2):
print("\nERROR: Shape of images mismatch")
return 1
# image format -> ndg or color
if len(np.shape(im1)) != 2:
print("\nERROR: Wrong image format")
return 1
#### Get data
# shape
ny,nx = np.shape(im1)[:2]
# Remove nan from arrays
self.filter()
# indexes
tmp = space2img(self.pos, camera)
ind = np.array(tmp,dtype = int)
# displacement [pixels]
disp = camera.velSpace2dispImage(self.vel)
# Prepare grid for interpolation
ij = np.meshgrid(np.arange(ny), np.arange(nx),indexing = "ij") # indexes of image (/!\ meshgrid gives [0] -> x (j) and [1] -> y (i))
gridij = np.array([ij[0],ij[1]]) # Create 2D grid containing i and j indexes (stored [i,j] -> [y,x])
gridij = np.moveaxis(gridij,0,2) # Reshape the data for compatibility -> [ni,nj,2] instead of [2,ni,nj]
# Interpolate displacement on all image
uv_i = griddata(ind, disp, gridij, method="cubic")
# Put velocities in grid and change nan to 0 (for interpolation function)
u = np.nan_to_num(uv_i[:,:,1], copy=True, nan=0.0, posinf=None, neginf=None)
v = np.nan_to_num(uv_i[:,:,0], copy=True, nan=0.0, posinf=None, neginf=None)
#### Warp image
# Data have to be double for interpolation function
imout = np.array(im1,dtype = np.float64)
imin = np.array(im2,dtype = np.float64)
# Use ctypes to wrap C function to Python (much much faster!!)
lib = ctypes.CDLL('/home/guillaume.bodart/Documents/04_Dev/03_Lib/Image/01_bicubic_c_func/bicubic_interpolation.so')
# Select the function in the dyn. lib
bicint = lib.bicubic_interpolation_warp
# Arrange input and outputs
bicint.restype = None
bicint.argtypes = [ndpointer(ctypes.c_double),
ndpointer(ctypes.c_double),
ndpointer(ctypes.c_double),
ndpointer(ctypes.c_double),
ctypes.c_int,
ctypes.c_int,
ctypes.c_bool]
try:
bicint(imin,u,v,imout,nx,ny,True)
except Exception as e:
print("\nERROR: Problem while warping image... \n Exception : " + str(e))
return 2
# Compute difference between im1 and backprojection
# First apply dynamic expansion based on the data within the warped area
ppout = dynExpans(imout[np.isfinite(uv_i[:,:,0])])
pp1 = dynExpans(im1[np.isfinite(uv_i[:,:,0])].astype(float))
# comput diff
diff = ppout - pp1
RMS = np.sqrt(np.sum([x**2 for x in diff])/len(diff))
if outRMS:
return np.array(np.clip(imout,0,255),dtype = np.uint8),RMS
else:
return np.array(np.clip(imout,0,255),dtype = np.uint8)
class filteredCl(dataCl):
......@@ -1990,15 +2125,17 @@ def image2space(imgIndx,camera,xy0 = np.nan):
- 1 - Information missing on the camera : no resolution and no scaling
- 2 - Image isn't squarred, problem with scaling...
- 3 - ResolutionInfo of the camera isn't recognize
- 4 - Shape of grid isn't correct'
"""
if camera.resolutionInfo == "":
# Case no resolution is entered, use simple scaling (typical blender case with othographic camera)
if camera.scale == -99:
print("ERROR : Information missing on the camera : no resolution and no scaling")
print("ERROR: Information missing on the camera : no resolution and no scaling")
return 1
# Image size is assumed to be squared
if camera.imagesize[0]!=camera.imagesize[1]:
print("ERROR : Image isn't squarred, problem with scaling...")
print("ERROR: Image isn't squarred, problem with scaling...")
return 1
res = camera.imagesize[0] / camera.scale
elif camera.resolutionInfo == "m/pix":
......@@ -2008,23 +2145,25 @@ def image2space(imgIndx,camera,xy0 = np.nan):
# Case resolution is given in pixels/meter
res = camera.resolution
else:
print("ERROR : resolution information of the camera isn't recognize. Expected : 'm/pix' or 'pix/m'")
print("ERROR: resolution information of the camera isn't recognize. Expected : 'm/pix' or 'pix/m'")
return 3
# Checkings
if len(np.shape(imgIndx)) > 2:
print("ERROR: Shape of grid isn't correct. Expected : [;,2]")
return 4
# If no xy0 are passed to the function they're calculated from the camera data
if np.isnan(xy0):
xy0 = (camera.position - (0.5*camera.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
print("WARNING: imgIndx is empty...")
# Swap column to get (j,i) as j is x and i is y and inverse i ax to be in mathematical coordinate space
JIimgIndx = np.array([imgIndx[:,1],camera.imagesize[0] - imgIndx[:,0]]).T
out = (JIimgIndx/res) + xy0
return out
def velSpace2dispImage(spaceVel,camera,xy0 = np.nan):
def space2img(spaceIndx,camera,xy0 = np.nan):
"""
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)
......@@ -2052,17 +2191,14 @@ def velSpace2dispImage(spaceVel,camera,xy0 = np.nan):
- 2 - Image isn't squarred, problem with scaling...
- 3 - ResolutionInfo of the camera isn't recognize
"""
# if isinstance(camera,transform.camera):
# print("\nERROR : wrong format of camera passed to the function. Expected : lspiv.transform.camera ")
# return 1
if camera.resolutionInfo == "":
# Case no resolution is entered, use simple scaling (typical blender case with othographic camera)
if camera.scale == -99:
print("\nERROR : Information missing on the camera : no resolution and no scaling")
print("ERROR : Information missing on the camera : no resolution and no scaling")
return 1
# Image size is assumed to be squared
if camera.imagesize[0]!=camera.imagesize[1]:
print("\nERROR : Image isn't squarred, problem with scaling...")
print("ERROR : Image isn't squarred, problem with scaling...")
return 1
res = camera.imagesize[0] / camera.scale
elif camera.resolutionInfo == "m/pix":
......@@ -2072,23 +2208,90 @@ def velSpace2dispImage(spaceVel,camera,xy0 = np.nan):
# Case resolution is given in pixels/meter
res = camera.resolution
else:
print("\nERROR : resolution information of the camera isn't recognize. Expected : 'm/pix' or 'pix/m'")
print("ERROR : resolution information of the camera isn't recognize. Expected : 'm/pix' or 'pix/m'")
return 3
# If no xy0 are passed to the function they're calculated from the camera data
if np.isnan(xy0):
xy0 = (camera.position - (0.5*camera.scale))[:2]
if len(spaceVel)==0:
print("\nWARNING : Values are empty...")
if len(spaceIndx)==0:
print("WARNING : imgIndx is empty...")
# Swap column to get (x,y) as j is x and i is y
JIimgIndx = (spaceIndx - xy0) * res # Because here the resolution is in pix/meters
out = np.array([camera.imagesize[0]-JIimgIndx[:,1],JIimgIndx[:,0]]).T
return out
def dynExpans(val):
mx, mn = val.max(),val.min()
return 255*((val-mn)/(mx-mn))
# def velSpace2dispImage(spaceVel,camera,xy0 = np.nan):
# """
# 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)
# Parameters
# ----------
# - imgIndx : [npArrayInt]
# Indexes of the images to be transformed. --> shape : [(ind),(i,j)]
# WARNING : image coordinate are expected.
# Values from Fudaa-LSPIV are usually given in a mathematical coordinate space with (0,0) located down left
# So the i index need to be "invert" --> correctiIndex = imageSize(i) - iIndexFudaa
# - camera : [lspiv.transform.camera()]
# Contains the information of the camera
# - xyz0 : [npArrayFloat]
# Offset (in meters) apply on the values X, Y for transformation
# Returns
# -------
# - SpaceCoord [npArrayFloat]
# X,Y values in meters
# - errno [int]
# In case of an error this function returns an integer
# - 1 - Information missing on the camera : no resolution and no scaling
# - 2 - Image isn't squarred, problem with scaling...
# - 3 - ResolutionInfo of the camera isn't recognize
# """
# # if isinstance(camera,transform.camera):
# # print("\nERROR : wrong format of camera passed to the function. Expected : lspiv.transform.camera ")
# # return 1
# if camera.resolutionInfo == "":
# # Case no resolution is entered, use simple scaling (typical blender case with othographic camera)
# if camera.scale == -99:
# print("\nERROR : Information missing on the camera : no resolution and no scaling")
# return 1
# # Image size is assumed to be squared
# if camera.imagesize[0]!=camera.imagesize[1]:
# print("\nERROR : Image isn't squarred, problem with scaling...")
# return 1
# res = camera.imagesize[0] / camera.scale
# elif camera.resolutionInfo == "m/pix":
# # Case resolution is given in meter/pixels (typical units from Fudaa-LSPIV project)
# res = 1/camera.resolution
# elif camera.resolutionInfo == "pix/m":
# # Case resolution is given in pixels/meter
# res = camera.resolution
# else:
# print("\nERROR : resolution information of the camera isn't recognize. Expected : 'm/pix' or 'pix/m'")
# return 3
# # If no xy0 are passed to the function they're calculated from the camera data
# if np.isnan(xy0):
# xy0 = (camera.position - (0.5*camera.scale))[:2]
# if len(spaceVel)==0:
# print("\nWARNING : Values are empty...")
# Compute displacement based on camera info
XYdisp = spaceVel / camera.fps
JIspaceV = (XYdisp)*res
# Swap column to get (i,j) as x is j and y is i
IJindex = np.array([JIspaceV[:,1],JIspaceV[:,0]]).T
# # Compute displacement based on camera info
# XYdisp = spaceVel / camera.fps
# JIspaceV = (XYdisp)*res
# # Swap column to get (i,j) as x is j and y is i
# IJindex = np.array([JIspaceV[:,1],JIspaceV[:,0]]).T
return IJindex
# return IJindex
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment