Interpoluj objętość 3D za pomocą numpy i lub scipy

The exact odpowiednik MATLAB 's interp3 będzie używał scipy' s interpn dla interpolacji jednorazowej:

import numpy as np
from scipy.interpolate import interpn

Vi = interpn((x,y,z), V, np.array([xi,yi,zi]).T)

Domyślną metodą dla MATLAB i scipy jest interpolacja liniowa, którą można zmienić za pomocą argumentu method. Zauważ, że tylko interpolacja liniowa i najbliższa jest obsługiwana przez interpn dla 3 wymiarów i powyżej, w przeciwieństwie do MATLAB, który obsługuje również interpolację sześcienną i splajnową.

Podczas wykonywania wielu wywołań interpolacyjnych na taką samą siatkę najlepiej zastosować obiekt interpolacji RegularGridInterpolator, jak w zaakceptowanej odpowiedzi powyżej . interpn używa RegularGridInterpolator wewnętrznie.

Zasadniczo, ndimage.map_coordinates działa we współrzędnych "indeksowych" (vel "voxel" lub "Pixel" współrzędne). Interfejs do niego wydaje się na początku trochę niezgrabny, ale daje dużo elastyczności.

Jeśli chcesz określić interpolowane współrzędne podobne do interp3 matlab, musisz przekonwertować współrzędne intput na współrzędne" indeksowe".

Istnieje również dodatkowa zmarszczka, która map_coordinates zawsze zachowuje dtype tablicy wejściowej na wyjściu. Jeśli Interpoluj tablicę liczb całkowitych, a otrzymasz wynik liczb całkowitych, który może być, ale nie musi być tym, czego chcesz. Dla poniższego fragmentu kodu, zakładam, że zawsze chcesz dane zmiennoprzecinkowe. (Jeśli nie, to w rzeczywistości jest to prostsze.)

Postaram się dodać więcej wyjaśnień później wieczorem (to raczej gęsty kod).

Podsumowując, Funkcja interp3, którą posiadam, jest bardziej złożona, niż może być potrzebna do Twoich dokładnych celów. Howver, powinno mniej więcej odwzorowywać zachowanie interp3 jak I zapamiętaj to (ignorując funkcję "zooming" interp3(data, zoom_factor), którą scipy.ndimage.zoom obsługuje.)

import numpy as np
from scipy.ndimage import map_coordinates

def main():
    data = np.arange(5*4*3).reshape(5,4,3)

    x = np.linspace(5, 10, data.shape[0])
    y = np.linspace(10, 20, data.shape[1])
    z = np.linspace(-100, 0, data.shape[2])

    # Interpolate at a single point
    print interp3(x, y, z, data, 7.5, 13.2, -27)

    # Interpolate a region of the x-y plane at z=-25
    xi, yi = np.mgrid[6:8:10j, 13:18:10j]
    print interp3(x, y, z, data, xi, yi, -25 * np.ones_like(xi))

def interp3(x, y, z, v, xi, yi, zi, **kwargs):
    """Sample a 3D array "v" with pixel corner locations at "x","y","z" at the
    points in "xi", "yi", "zi" using linear interpolation. Additional kwargs
    are passed on to ``scipy.ndimage.map_coordinates``."""
    def index_coords(corner_locs, interp_locs):
        index = np.arange(len(corner_locs))
        if np.all(np.diff(corner_locs) < 0):
            corner_locs, index = corner_locs[::-1], index[::-1]
        return np.interp(interp_locs, corner_locs, index)

    orig_shape = np.asarray(xi).shape
    xi, yi, zi = np.atleast_1d(xi, yi, zi)
    for arr in [xi, yi, zi]:
        arr.shape = -1

    output = np.empty(xi.shape, dtype=float)
    coords = [index_coords(*item) for item in zip([x, y, z], [xi, yi, zi])]

    map_coordinates(v, coords, order=1, output=output, **kwargs)

    return output.reshape(orig_shape)

main()

Pytanie jest stare, ale myślę, że wymaga wyjaśnienia, ponieważ nikt nie zauważył, że żądana operacja (interpolacja trójliniowa ) może być łatwo zaimplementowana od podstaw przy konsekwentnej oszczędności czasu obliczeniowego (około 10 razy szybciej) w.r.t. scipy.interpolate'S RegularGridInterpolator.

Kod

import numpy as np
from itertools import product

def trilinear_interpolation(x_volume, y_volume, z_volume, volume, x_needed, y_needed, z_needed):
    """
    Trilinear interpolation (from Wikipedia)

    :param x_volume: x points of the volume grid 
    :type crack_type: list or numpy.ndarray
    :param y_volume: y points of the volume grid 
    :type crack_type: list or numpy.ndarray
    :param x_volume: z points of the volume grid 
    :type crack_type: list or numpy.ndarray
    :param volume:   volume
    :type crack_type: list or numpy.ndarray
    :param x_needed: desired x coordinate of volume
    :type crack_type: float
    :param y_needed: desired y coordinate of volume
    :type crack_type: float
    :param z_needed: desired z coordinate of volume
    :type crack_type: float

    :return volume_needed: desired value of the volume, i.e. volume(x_needed, y_needed, z_needed)
    :type volume_needed: float
    """
    # dimensinoal check
    if np.shape(volume) != (len(x_volume), len(y_volume), len(z_volume)):
        raise ValueError(f'dimension mismatch, volume must be a ({len(x_volume)}, {len(y_volume)}, {len(z_volume)}) list or numpy.ndarray')
    # check of the indices needed for the correct control volume definition
    i = searchsorted(x_volume, x_needed)
    j = searchsorted(y_volume, y_needed)
    k = searchsorted(z_volume, z_needed)
    # control volume definition
    control_volume_coordinates = np.array(
        [[x_volume[i - 1], y_volume[j - 1], z_volume[k - 1]], [x_volume[i], y_volume[j], z_volume[k]]])
    xd = (np.array([x_needed, y_needed, z_needed]) - control_volume_coordinates[0]) / (control_volume_coordinates[1] - control_volume_coordinates[0])
    # interpolation along x
    c2 = [[0, 0], [0, 0]]
    for m, n in product([0, 1], [0, 1]):
        c2[m][n] = volume[i - 1][j - 1 + m][k - 1 + n] * (1 - xd[0]) + volume[i][j - 1 + m][k - 1 + n] * xd[0]
    # interpolation along y
    c1 = [0, 0]
    c1[0] = c2[0][0] * (1 - xd[1]) + c2[1][0] * xd[1]
    c1[1] = c2[0][1] * (1 - xd[1]) + c2[1][1] * xd[1]
    # interpolation along z
    volume_needed = c1[0] * (1 - xd[2]) + c1[1] * xd[2]
    return volume_needed

def searchsorted(l, x):
    for i in l:
        if i >= x: break
    return l.index(i)


from scipy.interpolate import RegularGridInterpolator
def trilin_interp_regular_grid(x_volume, y_volume, z_volume, volume, x_needed, y_needed, z_needed):
    # dimensinoal check
    if np.shape(volume) != (len(x_volume), len(y_volume), len(z_volume)):
        raise ValueError(f'dimension mismatch, volume must be a ({len(x_volume)}, {len(y_volume)}, {len(z_volume)}) list or numpy.ndarray')
    # trilinear interpolation on a regular grid
    fn = RegularGridInterpolator((x_volume,y_volume,z_volume), volume)
    volume_needed = fn(np.array([x_needed, y_needed, z_needed]))
    return volume_needed

Przykład

import numpy as np
import time

the_volume = np.array(
[[[0.902, 0.985, 1.12, 1.267, 1.366],
[0.822, 0.871, 0.959, 1.064, 1.141],
[0.744, 0.77, 0.824, 0.897, 0.954],
[0.669, 0.682, 0.715, 0.765, 0.806],
[0.597, 0.607, 0.631, 0.667, 0.695]],
[[1.059, 1.09, 1.384, 1.682, 1.881],
[0.948, 0.951, 1.079, 1.188, 1.251],
[0.792, 0.832, 0.888, 0.940, 0.971],
[0.726, 0.733, 0.754, 0.777, 0.792],
[0.642, 0.656, 0.675, 0.691, 0.700]]])

x_needed = np.linspace(100, 1000, 2)
y_needed = np.linspace(0.3, 1, 3)
z_needed = np.linspace(0, 1, 3)
start = time.time()
manual_trilint = []
for x in x_needed:
    for y in y_needed:
        for z in z_needed:
            manual_trilint.append(trilinear_interpolation([100, 1000], [0.2, 0.4, 0.6, 0.8, 1], [0, 0.2, 0.5, 0.8, 1], the_volume, x, y, z))
end = time.time()
print(end - start)

start = time.time()
auto_trilint = []
for x in x_needed:
    for y in y_needed:
        for z in z_needed:
            auto_trilint.append(trilin_interp_regular_grid([100, 1000], [0.2, 0.4, 0.6, 0.8, 1], [0, 0.2, 0.5, 0.8, 1], the_volume, x, y, z))
end = time.time()
print(end - start)