Python scipy.interpolate.RegularGridInterpolator() Examples

def __init__(self, fuel_map, full_load_curve):
        from scipy.interpolate import UnivariateSpline as Spl
        from scipy.interpolate import RegularGridInterpolator as Interpolator
        self.full_load_curve = full_load_curve
        self.fc = Interpolator(
            (fuel_map['speed'], fuel_map['power']), fuel_map['fuel'],
            bounds_error=False, fill_value=0
        )
        (s, p), fc = self.fc.grid, self.fc.values

        with np.errstate(divide='ignore', invalid='ignore'):
            e = np.maximum(0, p / fc)
        e[(p > full_load_curve(s)[:, None]) | (p < 0)] = np.nan
        b = ~np.isnan(e).all(1)
        (s, i), e = np.unique(s[b], return_index=True), e[b]
        b = ~np.isnan(e).all(0)
        p = p[b][np.nanargmax(e[:, b], 1)][i]

        func = Spl(s, p, w=1 / np.clip(p * .01, dfl.EPS, 1))
        s = np.unique(np.append(s, np.linspace(s.min(), s.max(), 1000)))
        p = func(s)
        self.max_power = p.max()
        self.speed_power = Spl(s, p, s=0)
        self.power_speed = Spl(*_invert(np.unique(p), s, p)[::-1], s=0, ext=3)
        self.idle_fc = self.fc((self.power_speed(0), 0)) 

def polychromatic(psfs, spectral_weights=None, interp_method='linear'):
        """Create a new PSF instance from an ensemble of monochromatic PSFs given spectral weights.

        The new PSF is the polychromatic PSF, assuming the wavelengths are
        sufficiently different that they do not interfere and the mode of
        imaging is incoherent.

        """
        if spectral_weights is None:
            spectral_weights = [1] * len(psfs)

        # find the most densely sampled PSF
        min_spacing = 1e99
        ref_idx = None
        ref_x = None
        ref_y = None
        ref_samples_x = None
        ref_samples_y = None
        for idx, psf in enumerate(psfs):
            if psf.sample_spacing < min_spacing:
                min_spacing = psf.sample_spacing
                ref_idx = idx
                ref_x = psf.x
                ref_y = psf.y
                ref_samples_x = psf.samples_x
                ref_samples_y = psf.samples_y

        merge_data = e.zeros((ref_samples_x, ref_samples_y, len(psfs)))
        for idx, psf in enumerate(psfs):
            # don't do anything to the reference PSF besides spectral scaling
            if idx is ref_idx:
                merge_data[:, :, idx] = psf.data * spectral_weights[idx]
            else:
                xv, yv = e.meshgrid(ref_x, ref_y)
                interpf = interpolate.RegularGridInterpolator((psf.y, psf.x), psf.data)
                merge_data[:, :, idx] = interpf((yv, xv), method=interp_method) * spectral_weights[idx]

        psf = PSF(data=merge_data.sum(axis=2), x=ref_x, y=ref_y)
        psf.spectral_weights = spectral_weights
        psf._renorm()
        return psf 

def interpolator(self, map_, station=None):
        maps = self.fetch_map(map_, station)
        nc = self.cell_count
        cc = (np.arange(nc[0]), np.arange(nc[1]), np.arange(nc[2]))
        return RegularGridInterpolator(cc, maps, bounds_error=False) 

def get_mask_surface_data(self, radius=None):
        points_interf = np.vstack(
            (self.model._surface_points.df['X'].values, self.model._surface_points.df['Y'].values)).T
        points_orient = np.vstack((self.model._orientations.df['X'].values, self.model._orientations.df['Y'].values)).T

        mask_interf = self.get_data_within_extent(points_interf)
        mask_orient = self.get_data_within_extent(points_orient)

        xj = self.model._grid.topography.values_3D[:, :, 0][0, :]
        yj = self.model._grid.topography.values_3D[:, :, 1][:, 0]
        zj = self.model._grid.topography.values_3D[:, :, 2].T

        interpolate = RegularGridInterpolator((xj, yj), zj)

        Z_interf_interp = interpolate(points_interf[mask_interf])
        Z_orient_interp = interpolate(points_orient[mask_orient])

        if radius is None:
            radius = np.diff(zj).max()
        print(radius)

        dist_interf = np.abs(Z_interf_interp - self.model._surface_points.df['Z'].values[mask_interf])
        dist_orient = np.abs(Z_orient_interp - self.model._orientations.df['Z'].values[mask_orient])

        surfmask_interf = dist_interf < radius
        surfmask_orient = dist_orient < radius
        surf_indexes = np.flatnonzero(mask_interf)[surfmask_interf]
        orient_indexes = np.flatnonzero(mask_orient)[surfmask_orient]

        mask_surfpoints = np.zeros(points_interf.shape[0], dtype=bool)
        mask_orient = np.zeros(points_orient.shape[0], dtype=bool)

        mask_surfpoints[surf_indexes] = True
        mask_orient[orient_indexes] = True

        return mask_surfpoints, mask_orient 

def init_interpolator(self, x_grid, y_grid, z_grid, vel):

        pass
        # for ivel in range(3):
        #     self.interpolator.append(interpolate.RegularGridInterpolator((x_grid, y_grid, z_grid),
        #                                                         vel[ivel,:,:,:],
        #                                                         bounds_error=False,
        #                                                         fill_value=self.settings['u_out'][ivel])) 

def create_interpolator(self, data, x_grid, y_grid, z_grid, i_dim):
        interpolator = interpolate.RegularGridInterpolator((x_grid, y_grid, z_grid),
                                                            data,
                                                            bounds_error=False,
                                                            fill_value=0.0)
        return interpolator 

def run_sample(self, filepath):
        """
        Run sampling.
        """
        timer = common.Timer()
        Rs = self.get_views()

        # As rendering might be slower, we wait for rendering to finish.
        # This allows to run rendering and fusing in parallel (more or less).

        depths = common.read_hdf5(filepath)

        timer.reset()
        tsdf = self.fusion(depths, Rs)

        xs = np.linspace(-0.5, 0.5, tsdf.shape[0])
        ys = np.linspace(-0.5, 0.5, tsdf.shape[1])
        zs = np.linspace(-0.5, 0.5, tsdf.shape[2])
        tsdf_func = rgi((xs, ys, zs), tsdf)

        modelname = os.path.splitext(os.path.splitext(os.path.basename(filepath))[0])[0]
        points = self.get_random_points(tsdf)
        values = tsdf_func(points)
        t_loc, t_scale = self.get_transform(modelname)

        occupancy = (values <= 0.)
        out_file = self.get_outpath(filepath)
        np.savez(out_file, points=points, occupancy=occupancy, loc=t_loc, scale=t_scale)

        print('[Data] wrote %s (%f seconds)' % (out_file, timer.elapsed())) 

def latlon(self, indices):
        """
        Compute the latitude and longitude from the given (i,j) indices
        of the grid

        Parameters
        ----------
        indices : list of tuples
            i, j points to compute latitude and longitude

        Returns
        -------
        out : tuple of ndarray
            list of lat,lon points from the given i,j indices

        Examples
        --------
        >>> a = [(23.4, 16.5), (3.66, 22.43)]
        >>> idx = g.latlon(a)
        """
        from scipy.interpolate import RegularGridInterpolator

        lati = RegularGridInterpolator((self.I[0, :], self.J[:, 0]),
                                       self.lat_rho.T)
        loni = RegularGridInterpolator((self.I[0, :], self.J[:, 0]),
                                       self.lon_rho.T)

        return (lati(indices), loni(indices)) 

def _make_grad(self):
        # Get the gradient of the Time-crossing map
        dx, dy, dz = self._dist_gradient()
        standard_grid = (np.arange(self._t.shape[0]), np.arange(self._t.shape[1]),
                         np.arange(self._t.shape[2]))
        self._grad = (RegularGridInterpolator(standard_grid, dx),
                      RegularGridInterpolator(standard_grid, dy),
                      RegularGridInterpolator(standard_grid, dz)) 

def interp_map_regular_grid(a, x, y, x_i, y_i, method="linear", debug=False, wrap=True):
    """Interpolates 2d fields from regular grid to another regular grid.

    wrap option: pads outer values/coordinates with other side of the array.
    Only works with lon/lat coordinates correctly.

    """
    # TODO these (interp_map*) should eventually end up in xgcm? Maybe not...
    # Pad borders to simulate wrap around coordinates
    # in global maps
    if wrap:

        x = x[[-1] + list(range(x.shape[0])) + [0]]
        y = y[[-1] + list(range(y.shape[0])) + [0]]

        x[0] = x[0] - 360
        x[-1] = x[-1] + 360

        y[0] = y[0] - 180
        y[-1] = y[-1] + 180

        a = a[:, [-1] + list(range(a.shape[1])) + [0]]
        a = a[[-1] + list(range(a.shape[0])) + [0], :]

    if debug:
        print("a shape", a.shape)
        print("x shape", x.shape)
        print("y shape", y.shape)
        print("x values", x[:])
        print("y values", y[:])
        print("x_i values", x_i[:])
        print("y_i values", y_i[:])

    xx_i, yy_i = np.meshgrid(x_i, y_i)
    f = spi.RegularGridInterpolator((x, y), a.T, method=method, bounds_error=False)
    int_points = np.vstack((xx_i.flatten(), yy_i.flatten())).T
    a_new = f(int_points)

    return a_new.reshape(xx_i.shape) 

def reload_calibration(self):
        '''
        reloads power calibration from file Windfreaktech_SynthHD.cal in instruments folder.
        '''
        try:
            data = np.loadtxt(qkit.cfg.get('user_insdir') + '/Windfreaktech_SynthHD.cal')
            f = data[1:, 0]
            p = data[0, 1:]
            values = data[1:, 1:]
            self._interp_amplitude = RegularGridInterpolator((f, p), values).__call__

        except IOError:
            raise IOError('Calibration file for WFT MW Source not found') 

def __init__(self, filename, verbose=1):
        """
        Initialisation of class to load templates from a file and create the interpolation
        objects

        Parameters
        ----------
        filename: string
            Location of Template file
        verbose: int
            Verbosity level,
            0 = no logging
            1 = File + interpolation point information
            2 = Detailed description of interpolation points
        """
        self.verbose = verbose
        if self.verbose:
            print("Loading lookup tables from", filename)

        grid, bins, template = self.parse_fits_table(filename)
        x_bins, y_bins = bins

        self.interpolator = interpolate.LinearNDInterpolator(
            grid, template, fill_value=0
        )
        self.nearest_interpolator = interpolate.NearestNDInterpolator(grid, template)

        self.grid_interp = interpolate.RegularGridInterpolator(
            (x_bins, y_bins),
            np.zeros([x_bins.shape[0], y_bins.shape[0]]),
            method="linear",
            bounds_error=False,
            fill_value=0,
        ) 

def __init__(self, x_norm, y_norm, z_map):
        self.xn = np.array(x_norm)
        self.yn = np.array(y_norm)
        self.zm = np.array(z_map)
        self.ip_x = Interpolate1D(pairs = zip(self.xn[:, 0], self.xn[:, 1]))
        self.ip_y = Interpolate1D(pairs = zip(self.yn[:, 0], self.yn[:, 1]))
        x_size, y_size = self.zm.shape
        self.ip_m = RegularGridInterpolator((np.arange(x_size), np.arange(y_size)), self.zm, method = "linear") 

def _make_interp_function_2d(self):
        """Generate a 2D interpolation function for this instance, used in sampling with exact_xy.

        Returns
        -------
        `scipy.interpolate.RegularGridInterpolator`
            interpolator instance.

        """
        if self.interpf_2d is None:
            self.interpf_2d = interpolate.RegularGridInterpolator((self.y, self.x), self.data)

        return self.interpf_2d 

def uniform_cart_to_polar(x, y, data):
    """Interpolate data uniformly sampled in cartesian coordinates to polar coordinates.

    Parameters
    ----------
    x : `numpy.ndarray`
        sorted 1D array of x sample pts
    y : `numpy.ndarray`
        sorted 1D array of y sample pts
    data : `numpy.ndarray`
        data sampled over the (x,y) coordinates

    Returns
    -------
    rho : `numpy.ndarray`
        samples for interpolated values
    phi : `numpy.ndarray`
        samples for interpolated values
    f(rho,phi) : `numpy.ndarray`
        data uniformly sampled in (rho,phi)

    """
    # create a set of polar coordinates to interpolate onto
    xmin, xmax = x.min(), x.max()
    ymin, ymax = y.min(), y.max()

    _max = max(abs(e.asarray([xmin, xmax, ymin, ymax])))

    rho = e.linspace(0, _max, len(x))
    phi = e.linspace(0, 2 * e.pi, len(y))
    rv, pv = e.meshgrid(rho, phi)

    # map points to x, y and make a grid for the original samples
    xv, yv = polar_to_cart(rv, pv)

    # interpolate the function onto the new points
    f = interpolate.RegularGridInterpolator((y, x), data, bounds_error=False, fill_value=0)
    return rho, phi, f((yv, xv), method='linear') 

def interpolate_var_to_points(self,
                                  points,
                                  variable,
                                  method='linear',
                                  indices=None,
                                  slices=None,
                                  mask=None,
                                  **kwargs):
        try:
            from scipy.interpolate import RegularGridInterpolator
        except ImportError:
            raise ImportError("The scipy package is required to use "
                              "Grid_R.interpolate_var_to_points\n"
                              " -- interpolating a regular grid")
        points = np.asarray(points, dtype=np.float64)
        just_one = (points.ndim == 1)
        points = points.reshape(-1, 2)
        if slices is not None:
            variable = variable[slices]
            if np.ma.isMA(variable):
                variable = variable.filled(0) #eventually should use Variable fill value
        x = self.node_lon if variable.shape[0] == len(self.node_lon) else self.node_lat
        y = self.node_lat if x is self.node_lon else self.node_lon
        interp_func = RegularGridInterpolator((x, y),
                                              variable,
                                              method=method,
                                              bounds_error=False,
                                              fill_value=0)
        if x is self.node_lon:
            vals = interp_func(points, method=method)
        else:
            vals = interp_func(points[:, ::-1], method=method)
        if just_one:
            return vals[0]
        else:
            return vals 

def nearest_2D_interpolator(lats_o, lons_o, values):
    '''
    Produces a 2D interpolator function using the nearest value interpolation
    method. If the grids are regular grids, uses the
    scipy.interpolate.RegularGridInterpolator,
    otherwise, scipy.intepolate.griddata

    Values can be interpolated from the returned function as follows:
       f = nearest_2D_interpolator(lat_origin, lon_origin, values_origin)
       interp_values = f(lat_interp, lon_interp)


    Parameters
    -----------
    lats_o: numpy.ndarray
        Latitude coordinates of the values usde by the interpolator
    lons_o: numpy.ndarray
        Longitude coordinates of the values usde by the interpolator
    values: numpy.ndarray
        Values usde by the interpolator


    Returns
    --------
    interpolator: function
        Nearest neighbour interpolator function
    '''
    # Determine if we are dealing with a regular grid
    if is_regular_grid(lats_o[2:-2, 2:-2], lons_o[2:-2, 2:-2]):
        return _nearest_2D_interpolator_reg(lats_o, lons_o, values)
    else:
        return _nearest_2D_interpolator_arb(lats_o, lons_o, values) 

def expected_LD(alpha, velocity, cl, lift_to_drag, airFrame):
    # data = data[data[:,0].argsort()]
    expected_value = 0
    # V = 12*45
    V = 1
    pdfs = []
    alpha = np.sort(np.unique(alpha.ravel()).T)
    velocity = np.sort(np.unique(velocity.ravel()).T)
    lift_to_drag = np.reshape(lift_to_drag, [len(alpha), len(velocity)])
    f_interpolation = RegularGridInterpolator((alpha, velocity), lift_to_drag)
    for i in range(len(airFrame.samples)):
        pdf = airFrame.pdf.score_samples(airFrame.samples[i,:])
        pdf = np.exp(pdf)
        # print(alpha.ravel().shape)
        # print(velocity.ravel().shape)
        # print(lift_to_drag.ravel().shape)
        # print(pdf)
        # print(airFrame.samples[i,:][0])
        
        data_i = C172.denormalize(np.array(airFrame.samples[i,:][0])).T
        try:
            LD_interpolated = f_interpolation(data_i)
        except(ValueError):
            print(data_i)
        # print(pdf, LD_interpolated)
        expected_value += LD_interpolated[0]
        pdfs.append(pdf)
    total_pdf = sum(pdfs)
    # print(total_pdf, expected_value, expected_value/len(airFrame.samples))
    expected_value = expected_value/len(airFrame.samples)
    return(expected_value) 

def expected_MonteCarlos(data, airFrame):
    # data = data[data[:,0].argsort()]
    alpha = data[:,0]
    velocity = data[:,1]
    cl = data[:,2]
    lift_to_drag = data[:,3]
    expected_value = 0
    # V = 12*45
    V = 1
    pdfs = []
    f_interpolation = RegularGridInterpolator((np.unique(alpha.ravel()).T, np.unique(velocity.ravel()).T), np.reshape(lift_to_drag, [200,200]))
    for i in range(len(airFrame.samples)):
        pdf = airFrame.pdf.score_samples(airFrame.samples[i,:])
        pdf = np.exp(pdf)
        # print(alpha.ravel().shape)
        # print(velocity.ravel().shape)
        # print(lift_to_drag.ravel().shape)
        # print(pdf)
        # print(airFrame.samples[i,:][0])
        
        data_i = C172.denormalize(np.array(airFrame.samples[i,:][0])).T
        try:
            LD_interpolated = f_interpolation(data_i)
        except(ValueError):
            print(data_i)
        # print(pdf, LD_interpolated)
        expected_value += (V/1)*pdf[0]*LD_interpolated[0]
        pdfs.append(pdf)
    total_pdf = sum(pdfs)
    expected_value /= total_pdf
    return(expected_value) 

def _nearest_2D_interpolator_reg(lats_o, lons_o, values, lats_d, lons_d):
    f = RegularGridInterpolator((np.flipud(lats_o[:, 0]), lons_o[0, :]),
                                np.flipud(values), method='nearest',
                                bounds_error=False)
    return f 

def make_nearest_interpolator_separated(field, grid=None):
	'''Make a nearest interpolator for a field on a separated grid.

	Parameters
	----------
	field : Field or ndarray
		The field to interpolate.
	grid : Grid or None
		The grid of the field. If it is given, the grid of `field` is replaced by this grid.

	Returns
	-------
	Field generator
		The interpolator, as a Field generator. The grid on which this field generator will evaluated, does
		not have to be separated.
	'''
	if grid is None:
		grid = field.grid
	else:
		field = Field(field, grid)

	axes_reversed = np.array(grid.separated_coords)
	interp = RegularGridInterpolator(axes_reversed, field.shaped, 'nearest', False)

	def interpolator(evaluated_grid):
		evaluated_coords = np.flip(np.array(evaluated_grid.coords), 0)
		res = interp(evaluated_coords.T)
		return Field(res.ravel(), evaluated_grid)

	return interpolator 

def bilinear_2D_interpolator(lats_o, lons_o, values):
    '''
    Produces a 2D interpolator function using the bilinear interpolation
    method. If the grids are regular grids, uses the
    scipy.interpolate.RegularGridInterpolator,
    otherwise, scipy.intepolate.griddata

    Values can be interpolated from the returned function as follows:
       f = nearest_2D_interpolator(lat_origin, lon_origin, values_origin)
       interp_values = f(lat_interp, lon_interp)


    Parameters
    -----------
    lats_o: numpy.ndarray
        Latitude coordinates of the values usde by the interpolator
    lons_o: numpy.ndarray
        Longitude coordinates of the values usde by the interpolator
    values: numpy.ndarray
        Values usde by the interpolator


    Returns
    --------
    interpolator: function
        Bilinear interpolator function
    '''
    if is_regular_grid(lats_o[2:-2, 2:-2], lons_o[2:-2, 2:-2]):
        return _bilinear_2D_interpolator_reg(lats_o, lons_o, values)
    else:
        return _bilinear_2D_interpolator_arb(lats_o, lons_o, values) 

def _bilinear_2D_interpolator_reg(lats_o, lons_o, values):
    f = RegularGridInterpolator((np.flipud(lats_o[:, 0]), lons_o[0, :]),
                                np.flipud(values), method='linear',
                                bounds_error=False)
    return f 

def make_linear_interpolator_separated(field, grid=None, fill_value=np.nan):
	'''Make a linear interpolators for a separated grid.

	Parameters
	----------
	field : Field or ndarray
		The field to interpolate.
	grid : Grid or None
		The grid of the field. If it is given, the grid of `field` is replaced by this grid.
	fill_value : scalar
		The value to use for points outside of the domain of the input field. If this is None, the values
		outside the domain are extrapolated.

	Returns
	-------
	Field generator
		The interpolator, as a Field generator. The grid on which this field generator will evaluated, does
		not have to be separated.
	'''
	if grid is None:
		grid = field.grid
	else:
		field = Field(field, grid)

	axes_reversed = np.array(grid.separated_coords)
	interp = RegularGridInterpolator(axes_reversed, field.shaped, 'linear', False, fill_value)

	def interpolator(evaluated_grid):
		evaluated_coords = np.flip(np.array(evaluated_grid.coords), 0)
		res = interp(evaluated_coords.T)
		return Field(res.ravel(), evaluated_grid)

	return interpolator 

def __init__(self, points, values, delta_list, extrapolate=False):
        # type: (Sequence[np.multiarray.ndarray], np.multiarray.ndarray, List[float], bool) -> None
        input_range = [(pvec[0], pvec[-1]) for pvec in points]
        DiffFunction.__init__(self, input_range, delta_list=delta_list)
        self._points = points
        self._extrapolate = extrapolate
        self.fun = interp.RegularGridInterpolator(points, values, method='linear',
                                                  bounds_error=not extrapolate,
                                                  fill_value=None) 

def interpolate_from_raster(input_raster, xy_array_to_interpolate,
                            method='nearest'):
    '''
    Interpolate input_raster to xy array.

    Parameters
    ----------
    input_raster : str
        filepath to raster
    xy_array : ndarray, shape (num_x, num_y)
        Node values to interpolate z at
    
    Returns
    -------
    interp_z : ndarray, shape (len(xy_array))
        Interpolated z values as xy_array (x, y)

    See Also
    --------

    See scipy.interpolate.RegularGridInterpolator documentation for further
    details.
    '''
    x_raster, y_raster, raster_grid = raster_XYZ(input_raster)
    interp_f = RegularGridInterpolator((y_raster, x_raster),
                                       raster_grid[2],
                                       bounds_error=False,
                                       fill_value=np.nan)
    # Need to flip
    xy_flipped = np.fliplr(xy_array_to_interpolate)
    interp_z = interp_f(xy_flipped, method=method)

    return interp_z 

def __getitem__(self, idx):
        """Get the random cutout data cube corresponding to idx.

        Args:
          idx: int, index of the crop to return. must be smaller than len(self).

        Returns:
          space_time_crop_hres (*optional): array of shape [4, nt_hres, nz_hres, nx_hres],
          where 4 are the phys channels pbuw.
          space_time_crop_lres: array of shape [4, nt_lres, nz_lres, nx_lres], where 4 are the phys
          channels pbuw.
          point_coord: array of shape [n_samp_pts_per_crop, 3], where 3 are the t, x, z dims.
                       CAUTION - point_coord are normalized to (0, 1) for the relative window.
          point_value: array of shape [n_samp_pts_per_crop, 4], where 4 are the phys channels pbuw.
        """
        t_id, z_id, x_id = self.rand_start_id[idx]
        space_time_crop_hres = self.data[:,
                                         t_id:t_id+self.nt_hres,
                                         z_id:z_id+self.nz_hres,
                                         x_id:x_id+self.nx_hres]  # [c, t, z, x]

        # create low res grid from hi res space time crop
        # apply filter
        space_time_crop_hres_fil = self.filter(space_time_crop_hres)

        interp = RegularGridInterpolator(
            (np.arange(self.nt_hres), np.arange(self.nz_hres), np.arange(self.nx_hres)),
            values=space_time_crop_hres_fil.transpose(1, 2, 3, 0), method=self.lres_interp)

        lres_coord = np.stack(np.meshgrid(np.linspace(0, self.nt_hres-1, self.nt_lres),
                                          np.linspace(0, self.nz_hres-1, self.nz_lres),
                                          np.linspace(0, self.nx_hres-1, self.nx_lres),
                                          indexing='ij'), axis=-1)
        space_time_crop_lres = interp(lres_coord).transpose(3, 0, 1, 2)  # [c, t, z, x]

        # create random point samples within space time crop
        point_coord = np.random.rand(self.n_samp_pts_per_crop, 3) * (self.scale_hres - 1)
        point_value = interp(point_coord)
        point_coord = point_coord / (self.scale_hres - 1)

        if self.normalize_output:
            space_time_crop_lres = self.normalize_grid(space_time_crop_lres)
            point_value = self.normalize_points(point_value)
        if self.normalize_hres:
            space_time_crop_hres = self.normalize_grid(space_time_crop_hres)

        return_tensors = [space_time_crop_lres, point_coord, point_value]

        # cast everything to float32
        return_tensors = [t.astype(np.float32) for t in return_tensors]

        if self.return_hres:
            return_tensors = [space_time_crop_hres] + return_tensors
        return tuple(return_tensors) 

def elastic_transform_3d(img_numpy, labels=None, alpha=1, sigma=20, c_val=0.0, method="linear"):
    """
    :param img_numpy: 3D medical image modality
    :param labels: 3D medical image labels
    :param alpha: scaling factor of gaussian filter
    :param sigma: standard deviation of random gaussian filter
    :param c_val: fill value
    :param method: interpolation method. supported methods : ("linear", "nearest")
    :return: deformed image and/or label
    """
    assert img_numpy.ndim == 3, 'Wrong img shape, provide 3D img'
    if labels is not None:
        assert img_numpy.shape == labels.shape, "Shapes of img and label do not much!"
    shape = img_numpy.shape

    # Define 3D coordinate system
    coords = np.arange(shape[0]), np.arange(shape[1]), np.arange(shape[2])

    # Interpolated img
    im_intrps = RegularGridInterpolator(coords, img_numpy,
                                        method=method,
                                        bounds_error=False,
                                        fill_value=c_val)

    # Get random elastic deformations
    dx = gaussian_filter((np.random.rand(*shape) * 2 - 1), sigma,
                         mode="constant", cval=0.) * alpha
    dy = gaussian_filter((np.random.rand(*shape) * 2 - 1), sigma,
                         mode="constant", cval=0.) * alpha
    dz = gaussian_filter((np.random.rand(*shape) * 2 - 1), sigma,
                         mode="constant", cval=0.) * alpha

    # Define sample points
    x, y, z = np.mgrid[0:shape[0], 0:shape[1], 0:shape[2]]
    indices = np.reshape(x + dx, (-1, 1)), \
              np.reshape(y + dy, (-1, 1)), \
              np.reshape(z + dz, (-1, 1))

    # Interpolate 3D image image
    img_numpy = im_intrps(indices).reshape(shape)

    # Interpolate labels
    if labels is not None:
        lab_intrp = RegularGridInterpolator(coords, labels,
                                            method="nearest",
                                            bounds_error=False,
                                            fill_value=0)

        labels = lab_intrp(indices).reshape(shape).astype(labels.dtype)
        return img_numpy, labels

    return img_numpy 

def read_grid(self, i_grid, i_cache=0):
        """
        This function returns an interpolator list of size 3 made of `scipy.interpolate.RegularGridInterpolator`
        objects.
        """
        velocities = ['ux', 'uy', 'uz']
        interpolator = list()
        for i_dim in range(3):
            file_name = self.grid_data['grid'][i_grid][velocities[i_dim]]['file']
            if i_cache == 0:
                if not self.settings['store_field']:
                    # load file, but dont copy it
                    self.vel_holder0[i_dim] = np.memmap(self.route + '/' + file_name,
                                                   # dtype='float64',
                                                   dtype=self.grid_data['grid'][i_grid][velocities[i_dim]]['Precision'],
                                                   shape=(self.grid_data['dimensions'][2],
                                                          self.grid_data['dimensions'][1],
                                                          self.grid_data['dimensions'][0]),
                                                   order='F')
                else:
                    # load and store file
                    self.vel_holder0[i_dim] = (np.fromfile(open(self.route + '/' + file_name, 'rb'),
                                                          dtype=self.grid_data['grid'][i_grid][velocities[i_dim]]['Precision']).\
                                                          reshape((self.grid_data['dimensions'][2],
                                                                   self.grid_data['dimensions'][1],
                                                                   self.grid_data['dimensions'][0]),
                                                                   order='F'))

                interpolator.append(self.create_interpolator(self.vel_holder0[i_dim],
                                                        self.grid_data['initial_x_grid'],
                                                        self.grid_data['initial_y_grid'],
                                                        self.grid_data['initial_z_grid'],
                                                        i_dim=i_dim))
            elif i_cache == 1:
                if not self.settings['store_field']:
                    # load file, but dont copy it
                    self.vel_holder1[i_dim] = np.memmap(self.route + '/' + file_name,
                                                   # dtype='float64',
                                                   dtype=self.grid_data['grid'][i_grid][velocities[i_dim]]['Precision'],
                                                   shape=(self.grid_data['dimensions'][2],
                                                          self.grid_data['dimensions'][1],
                                                          self.grid_data['dimensions'][0]),
                                                   order='F')
                else:
                    # load and store file
                    self.vel_holder1[i_dim] = (np.fromfile(open(self.route + '/' + file_name, 'rb'),
                                                          dtype=self.grid_data['grid'][i_grid][velocities[i_dim]]['Precision']).\
                                                          reshape((self.grid_data['dimensions'][2],
                                                                   self.grid_data['dimensions'][1],
                                                                   self.grid_data['dimensions'][0]),
                                                                   order='F'))

                interpolator.append(self.create_interpolator(self.vel_holder1[i_dim],
                                                        self.grid_data['initial_x_grid'],
                                                        self.grid_data['initial_y_grid'],
                                                        self.grid_data['initial_z_grid'],
                                                        i_dim=i_dim))
            else:
                raise Error('i_cache has to be 0 or 1')
        return interpolator 

def nonmax(img, sigma=2, threshold=0):
    '''
    Finds directional local maxima
    in a gradient array, as used in the Canny edge detector, but made
    separately accessible here for greater flexibility. The result is a
    logical array with the value true where the gradient magnitude is a
    local maximum along the gradient direction.
    '''

    # Get normalised gaussian gradients
    eps = 1e-12
    gx = gaussian_filter1d(img, sigma, axis=0, order=1)
    gy = gaussian_filter1d(img, sigma, axis=1, order=1)
    gz = gaussian_filter1d(img, sigma, axis=2, order=1)
    gmag = np.sqrt(gx**2 + gy**2 + gz**2)

    gx = gx / (gmag + eps)
    gy = gy / (gmag + eps)
    gz = gz / (gmag + eps)
    standard_grid = (np.arange(gmag.shape[0]), np.arange(gmag.shape[1]),
                     np.arange(gmag.shape[2]))
    ginterp = RegularGridInterpolator(standard_grid, gmag, bounds_error=False)

    # Interpolate the graident magnitudes
    idx = np.argwhere(
        img > threshold
    )  # Double-check if the original image should be used to check
    xidx = idx[:, 0]
    yidx = idx[:, 1]
    zidx = idx[:, 2]
    dx = gx[xidx, yidx, zidx]
    dy = gy[xidx, yidx, zidx]
    dz = gz[xidx, yidx, zidx]
    gmag_0 = gmag[xidx, yidx, zidx]
    gmag_1 = ginterp(np.stack((xidx + dx, yidx + dy, zidx + dz), axis=-1))
    gmag_2 = ginterp(np.stack((xidx - dx, yidx - dy, zidx - dz), axis=-1))

    # Suppress nonmax voxels
    keep = np.logical_and(gmag_0 > gmag_1, gmag_0 > gmag_2)
    gmag.fill(0)
    gmag[xidx, yidx, zidx] = keep.astype('float')

    return gmag