Source

Target

Commits (3)
Showing with 1576 additions and 141 deletions
+1576 -141
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
#!/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
This diff is collapsed.
...@@ -47,6 +47,7 @@ class camera: ...@@ -47,6 +47,7 @@ class camera:
self.imagesize = () # size of image (tuple) self.imagesize = () # size of image (tuple)
self.a = [] # a matrix from DLT, used to convert world to image space 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.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 = ""): def read_a(self,filename = ""):
""" """
...@@ -152,7 +153,7 @@ class camera: ...@@ -152,7 +153,7 @@ class camera:
return 0 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. 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) 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: ...@@ -201,13 +202,15 @@ class camera:
return 3 return 3
# If no xy0 are passed to the function they're calculated from the camera data # 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] xy0 = (self.position - (0.5*self.scale))[:2]
if len(imgIndx)==0: if len(imgIndx)==0:
print("WARNING : imgIndx is empty...") print("WARNING : imgIndx is empty...")
# Swap column to get (j,i) as j is x and i is y # Swap column to get (j,i) as j is x and i is y
JIimgIndx = np.array([imgIndx[:,1],imgIndx[:,0]]).T if ni == True:
ni = np.max(imgIndx[:,0])
JIimgIndx = np.array([imgIndx[:,1],ni-imgIndx[:,0]]).T
out = (JIimgIndx/res) + xy0 out = (JIimgIndx/res) + xy0
return out return out
...@@ -461,4 +464,85 @@ def read_grid(filename = ""): ...@@ -461,4 +464,85 @@ def read_grid(filename = ""):
return i,j 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
This diff is collapsed.