Source code for ecpi.common.num.interpol_tools

'''
#TODO: to remove
'''

import logging
from  scipy.interpolate import RectBivariateSpline, splev, splrep
import numpy as np 
import matplotlib.pyplot as plt

logger = logging.getLogger(__name__) 


[docs]class InterpolInversMonotoneFunction(object): """Interpolation of a monotone function. """ def __init__(self, a_x, a_f_x, degree_spline=1): assert a_x.size == a_f_x.size # do spline, and inverse x , y self.x = a_x self.y = a_f_x self.spl = splrep(a_x, a_f_x, k=degree_spline) logger.debug(f'(knot vector, bspline coeff, bspline degree): {self.spl[0]}') logger.debug(f'weighted sum of residuals: {self.spl[1]}') if self.spl[2] <= 0: # cf. scipy.interpolate.splrep.html logger.info(f'bspline fitting status: {self.spl[2]}') else: logger.error(f'bspline fitting status: {self.spl[2]}') def __call__(self, a_y): return splev(a_y, self.spl)
[docs]class InterpolLin2DCenter(object): """ URL: https://scipython.com/book/chapter-8-scipy/examples/ two-dimensional-interpolation-with-scipyinterpolaterectbivariatespline/ input : image and vector translation in pixel unit output : interpol image only in domain center after translation """ def __init__(self, p_ima, p_ox, p_oy): # coherence test if (p_ima.shape[0] % 2 == 1): raise if (p_ox % 2 == 1) or (p_oy % 2 == 1): raise if (p_ox > p_ima.shape[0]) or (p_oy > p_ima.shape[1]): raise self.sox = p_ox self.soy = p_oy self.vmin = np.min(p_ima) self.vmax = np.max(p_ima) # 0 padding to obtain 0 out of domain with interpol method ima_0 = np.zeros((p_ima.shape[0] + 2, p_ima.shape[1] + 2), dtype=np.float32) ima_0[1:-1, 1:-1] = np.array(p_ima, dtype=np.float32) # plan interpol ix = np.arange(0, ima_0.shape[0], dtype=np.float32) + 0.5 iy = np.arange(0, ima_0.shape[1], dtype=np.float32) + 0.5 self.o_interpol = RectBivariateSpline(ix, iy, ima_0, kx=1, ky=1, s=0) # create center array self.aix = np.arange(0, p_ox, dtype=np.float32) + (ima_0.shape[0] - p_ox) / 2 + 0.5 self.aiy = np.arange(0, p_oy, dtype=np.float32) + (ima_0.shape[1] - p_oy) / 2 + 0.5
[docs] def interpol(self, p_dx, p_dy): # add translate to center array logger.debug(f'offset: ({p_dx}, {p_dy})') aixd = self.aix + p_dx aiyd = self.aiy + p_dy # interpol aimx, aimy = np.meshgrid(aixd, aiyd, indexing='ij') out_inter = self.o_interpol.ev(aimx, aimy) self.last = [aimx, aimy, out_inter] return out_inter
[docs] def plot_interpol(self, p_title=""): # pragma: no cover plt.figure() plt.title(p_title) plt.imshow(self.last[2], vmin=self.vmin, vmax=self.vmax) # plt.imshow(self.last[2] , origin='lower') plt.grid(linewidth=1) plt.colorbar()