Как реализовать алгоритм трилинейной интерполяции для 3-х переменных?

у меня есть код, на с++, реализации билинейной интерполяции для двух переменных

class point{
public:
    double x,y,f;
    void init(double x1,double y1,double f1){
        x=x1;
        y=y1;
        f=f1;
    }
    point(){
    }
}*uk;
int n_x,n_y;
double interpolate(double x,double y){
    double h_x,h_y;
    h_x=uk[1].x-uk[0].x;
    h_y=uk[n_y].y-uk[0].y;
    int n,m;
    n=int(floor((x-uk[0].x)/h_x));
    m=int(floor((y-uk[0].y)/h_y));
    if((n<0)||(n>n_x)||(m<0)||(m>n_y)) return 0;
    if (n==n_x-1) n--;
    if (m==n_y-1) m--;
    double x_n,y_m;
    x_n=uk[n].x;
    y_m=uk[m*n_y].y;
    return  uk[n+m*n_y].f*(x_n+h_x-x)*(y_m+h_y-y)/h_x/h_y +
            uk[n+1+m*n_y].f*(x-x_n)*(y_m+h_y-y)/h_x/h_y +
            uk[n+(m+1)*n_y].f*(x_n+h_x-x)*(y-y_m)/h_x/h_y  +
            uk[n+1+(m+1)*n_y].f*(x-x_n)*(y-y_m)/h_x/h_y;
}

Может его можно будет как нибудь изменить? Нужно что бы в функцию передавались 3 значения X Y Z после отработки функции, выходные данные нужно будет сохранить в файл.
пример выходного файла:

5 -9 5.9
-1 6 1.1
-5 -7 5.5


Ответы (1 шт):

Автор решения: mr NAE

Что необходимо сделать так это заменить вектор с именованными компонентами на массив. Это может и не так удобно читать, но очень удобно программировать.

Пусть, inputChannels - размерность вектора, std::array<int, inputChannels> grids - размерность каждой компоненты вектора в массиве uk (далее data)

template <typename T, size_t inputChannels>
std::array<T, inputChannels> interpolate(const std::array<T, inputChannels> &value)
{
    //value is always in [0,1], so remap value before (after) call this
    std::array<T,  inputChannels> result;
    result.fill(-1);
    std::array<int, 1 << inputChannels> positions;
    positions.fill(0);
    std::array<T, positions.size()> weights;
    weights.fill(1);
    for (int i = 0; i < inputChannels; ++i) {
        // 1
        T reducedValue = (grids[i] - 1) * value[i];
        if (reducedValue < 0)
            return result;
        const bool hasNext = reducedValue + 1 < grids[i];
        const int pos = hasNext ? static_cast<int>(reducedValue) : (grids[i] - 1);
        const T gain = hasNext ? reducedValue - pos : 0;
        for (auto positionNumber = positions.size(); positionNumber-- > 0;) {
            const bool isPositionOdd = ((positionNumber >> i) & 1) > 0;
            positions[positionNumber] = grids[i] * positions[positionNumber] + pos + (hasNext && isPositionOdd ? 1 : 0);
            weights[positionNumber] *= isPositionOdd ? gain : 1 - gain;
        }
    }
    for (auto j =  inputChannels; j-- > 0;) {
        T sum = 0;
        for (auto i = positions.size(); i-- > 0;)
            sum += weights[i] * data[j +  inputChannels * positions[i]];
        result[j] = sum;
    }
    return result;
}
  1. Здесь предполагается, что каждая компонента нормирована к [0,1] включительно. Если это не так, то в отмеченном 1 месте следует задать своё отображение значений на индекс массива.

Остальное - это полный перебор комбинаций множителей для каждой компоненты соседних точек для заданного вектора. Например, для двумерного случая существуют четыре соседних точки: f00,f10,f01,f11. Интерполированное значение:

f00 * (1-a) * (1-b)+
f01 * (1-a) *    b +
f10 *    a  * (1-b)+
f11 *    a  *    b

Следовательно, один раз вычислив, скажем, b домножаем вес каждой точки на (1-b) или b в зависимости от номера точки. И так для каждой компоненты.

→ Ссылка