From e7cd5aa568701feac931b396cf03a25f84bafd74 Mon Sep 17 00:00:00 2001
From: Guillaume Bodart <guillaume.bodart@inrae.fr>
Date: Mon, 5 Jul 2021 14:10:18 +0200
Subject: [PATCH] MAJ latest dev
---
project.py | 269 ++++++++++++++++++++++++++++++++++++++++++++++-------
1 file changed, 236 insertions(+), 33 deletions(-)
diff --git a/project.py b/project.py
index 43ea412..6246b7a 100644
--- a/project.py
+++ b/project.py
@@ -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
--
GitLab