orthocyd.py

+#!/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

transform.py

@@ -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