orthocyd.py 8.14 KB
Newer Older
Guillaume  Bodart's avatar
Guillaume Bodart committed
#!/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